RayCasting with PyOpenGL

  1import numpy, wx
  2from numpy import arange, array, float32, int8
  3
  4from VolRenderSkel import *
  5
  6fragment_shader_src = '''
  7        uniform sampler1D TransferFunction;
  8        uniform sampler3D VolumeData;
  9        uniform sampler2D RayEnd;
 10        uniform sampler2D RayStart;
 11        
 12        float stepSize = 0.001; 
 13        
 14        void main(void)
 15        {
 16        
 17            // Get the end point of the ray (from the front-culled faces rendering)
 18            vec3 rayStart = texture2D(RayStart, gl_TexCoord[0].st).xyz;
 19            vec3 rayEnd = texture2D(RayEnd, gl_TexCoord[0].st).xyz;
 20            
 21            // Get a vector from back to front
 22            vec3 traverseVector = rayEnd - rayStart;
 23
 24            // The maximum length of the ray
 25            float maxLength = length(traverseVector);
 26              
 27            // Construct a ray in the correct direction and of step length
 28            vec3 step = stepSize * normalize(traverseVector);
 29            vec3 ray = step;
 30       
 31            // The color accumulation buffer
 32            vec4 acc = vec4(0.0, 0.0, 0.0, 0.0);
 33
 34            // Holds current voxel color
 35            vec4 voxelColor;
 36
 37            // Advance ray
 38            for (int i = 0; i < int(1/stepSize); ++i)
 39            {
 40                if (length(ray) >= maxLength || acc.a >= 0.99) 
 41                {
 42                    acc.a = 1.0;
 43                    break;
 44                }
 45
 46                voxelColor = texture1D(TransferFunction, texture3D(VolumeData, ray + rayStart).w);
 47
 48                // Accumulate RGB : acc.rgb = voxelColor.rgb*voxelColor.a + (1.0 - voxelColor.a)*acc.rgb;
 49                acc.rgb = mix(acc.rgb, voxelColor.rgb, voxelColor.a);
 50
 51                // Accumulate Opacity: acc.a = acc.a + (1.0 - acc.a)*voxelColor.a;
 52                acc.a = mix(voxelColor.a, 1.0, acc.a);
 53
 54                ray += step;
 55
 56            }
 57        
 58            gl_FragColor = acc;
 59            return;
 60        }
 61        '''
 62
 63class RayCaster(VolumeRenderSkeleton):
 64       
 65    def __init__(self, parent):
 66        
 67        VolumeRenderSkeleton.__init__(self, parent)
 68
 69        self.fragment_src_file = 'raycast.f.c'
 70        self.vertex_src_file = 'raycast.v.c'
 71
 72        self.data_scale = array([1.0, 1.0, 1.0], dtype=float32)
 73        self.iso_value = 0.0
 74
 75        # Shader sources
 76        self.vertex_color_shader = '''
 77                void main(void)
 78                {
 79                    gl_FrontColor = vec4(gl_Vertex.xyz, 1.0);
 80                    gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
 81                }
 82                '''
 83        
 84        # Default
 85        self.fragment_shader_src = fragment_shader_src
 86        self.vertex_shader_src = '''
 87            void main(void)
 88            {
 89                gl_TexCoord[0] = gl_Vertex;
 90                gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
 91            }
 92        '''
 93
 94
 95    def InitGL(self):
 96        VolumeRenderSkeleton.InitGL(self)
 97        self.color_program = compile_program(self.vertex_color_shader, None)        
 98        self.BuildGeometry()
 99
100        return
101    
102    def SetupUniforms(self):
103        VolumeRenderSkeleton.SetupUniforms(self)
104
105        glActiveTexture(GL_TEXTURE1)
106        glBindTexture(GL_TEXTURE_3D, self.vol_data)
107        glUniform1i(glGetUniformLocation(self.program, "VolumeData"), 1)
108        glUniform1i(glGetUniformLocation(self.program, "RayEnd"), 2)
109        glUniform1i(glGetUniformLocation(self.program, "RayStart"), 3)
110
111    def OnDraw(self,event):
112        self.SetCurrent()
113        if not self.init:
114            self.InitGL()
115            self.init = True
116
117        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
118        glLoadIdentity()
119
120        glTranslate(0.0, 0.0, -2.6)
121        glPushMatrix()
122        glRotate(self.rotation_y, 0.0, 1.0, 0.0)
123        glRotate(self.rotation_x, 1.0, 0.0, 0.0)
124        glTranslate(-0.5, -0.5, -0.5)
125
126        # Render the back texture
127        glEnable(GL_CULL_FACE)
128        glCullFace(GL_FRONT)
129        glUseProgram(self.color_program)
130        glBindBuffer(GL_ARRAY_BUFFER, self.cube_vbo)
131        glInterleavedArrays(GL_V3F, 0, None)
132        glDrawArrays(GL_QUADS, 0, 4*8)
133        glBindBuffer(GL_ARRAY_BUFFER, 0)
134
135        # Get the texture data
136        glActiveTexture(GL_TEXTURE2)
137        glBindTexture(GL_TEXTURE_2D, self.back_texture)
138        glCopyTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, 0, 0, 512, 512)
139
140        # Render the front texture
141        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
142        glCullFace(GL_BACK)
143        glBindBuffer(GL_ARRAY_BUFFER, self.cube_vbo)
144        glInterleavedArrays(GL_V3F, 0, None)
145        glDrawArrays(GL_QUADS, 0, 4*8)
146        # Get the texture data
147        glActiveTexture(GL_TEXTURE3)
148        glBindTexture(GL_TEXTURE_2D, self.front_texture)
149        glCopyTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, 0, 0, 512, 512)
150
151        # Clear the screen
152        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
153        glDisable(GL_CULL_FACE)
154        glBindBuffer(GL_ARRAY_BUFFER, 0)
155        glUseProgram(None)
156
157        # Render the wireframe box
158        glDisableClientState(GL.GL_COLOR_ARRAY)
159        glColor(0.0, 1.0, 0.0)
160        glDisable(GL_LIGHTING)
161        glDisable(GL_TEXTURE_2D)
162        glVertexPointerf(self.box)
163        glDrawArrays(GL_LINES, 0, len(self.box))
164        glPopMatrix()
165        
166        # Render the cube
167        glUseProgram(self.program)
168        
169        # Setup the texture information
170        self.SetupUniforms()
171
172        # Render flat quad and fudge it to the center 
173        glTranslate(-0.58, -0.58, -0.5)
174        glColor(0.0, 0.0, 1.0)
175        glVertexPointerf(self.flat_plane)
176        glDrawArrays(GL_QUADS, 0, 4)
177        
178        self.SwapBuffers()
179        
180        return
181
182    def BuildGeometry(self):
183        
184        cube_panels = []
185
186        for t in ['xy', 'xz', 'yz']:
187            for p in [0.0, 1.0]:
188                plane_verts = gen_plane(t=t, p=p)
189
190                if not int(p):
191                    plane_verts.reverse()
192
193                for pv in plane_verts:
194                    # Colors
195                    cube_panels.extend(pv)
196         
197        self.cube_vbo = simple.GLuint(0)
198        self.cube_data = numpy.array(cube_panels, dtype=numpy.float32)
199        glGenBuffers(1, self.cube_vbo)
200        glBindBuffer(GL_ARRAY_BUFFER, self.cube_vbo)
201        glBufferData(GL_ARRAY_BUFFER, ADT.arrayByteCount(self.cube_data), 
202                     ADT.voidDataPointer(self.cube_data), GL_STATIC_DRAW_ARB)
203        
204        self.flat_plane = numpy.array(gen_plane(t='xy', p=2.0));
205        
206        self.box = [[0.0, 0.0, 0.0],[0.0, 0.0, 1.0],
207                    [1.0, 0.0, 0.0],[1.0, 0.0, 1.0],
208                    [1.0, 1.0, 0.0],[1.0, 1.0, 1.0],
209                    [0.0, 1.0, 0.0],[0.0, 1.0, 1.0]]
210        
211        self.box.extend(box_side())
212        self.box.extend(box_side(z=1.0))
213
214    def LoadVolumeData(self):
215
216        # Load the volume data
217        buffer = open('data/MRI_head.dat', 'rb').read()
218        data = numpy.frombuffer(buffer, numpy.int8)
219        self.vol_data = glGenTextures(1)
220        glPixelStorei(GL_UNPACK_ALIGNMENT,1)
221        glBindTexture(GL_TEXTURE_3D, self.vol_data)
222        glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP)
223        glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP)
224        glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP)
225        glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
226        glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
227    
228        glTexImage3D(GL_TEXTURE_3D, 0, GL_ALPHA, 128, 128, 128, 0, GL_ALPHA, GL_UNSIGNED_BYTE, data)
229    
230        # Create a texture for the ends of the rays
231        self.back_texture = glGenTextures(1)
232        glPixelStorei(GL_UNPACK_ALIGNMENT,1)
233        glBindTexture(GL_TEXTURE_2D, self.back_texture)
234        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP)
235        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP)
236        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
237        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
238    
239        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, 512, 512, 0, GL_RGB, GL_UNSIGNED_BYTE, None)
240        
241        # Create a texture for the start of the rays
242        self.front_texture = glGenTextures(1)
243        glPixelStorei(GL_UNPACK_ALIGNMENT,1)
244        glBindTexture(GL_TEXTURE_2D, self.front_texture)
245        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP)
246        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP)
247        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
248        glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
249    
250        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, 512, 512, 0, GL_RGB, GL_UNSIGNED_BYTE, None)
251        return
252
253if __name__ == '__main__':
254    app = wx.App()
255    glutInit()
256    frame = wx.Frame(None, -1, 'Volume RayCasting with PyOpenGL', wx.DefaultPosition, wx.Size(600, 600))
257    canvas = RayCaster(frame)
258
259    frame.Show()
260    app.MainLoop()
261
262    # Cleanup