License New BSD license
Lines 473
Keywords
OpenGL (15) PyOpenGL (9) PyQt (6) Qt (4) Volume Rendering (9)
Included in this Library
Permissions
Owner: Stou S.
Viewable by Everyone
Editable by All Siafoo Users
Hide
Writing an article is easy - try our reStructured Text demo Join Siafoo Now or Learn More

PyQt4 and PyOpenGL Volume Rendering Skeleton Atom Feed 0

In Brief Keys:... more
# 's
  1# Created on Nov 8, 2009
2#
3# Copyright (c) 2010, Stou Sandalski (stou@icapsid.net)
4# All rights reserved.
5#
6# Redistribution and use in source and binary forms, with or without
7# modification, are permitted provided that the following conditions are met:
8# * Redistributions of source code must retain the above copyright
9# notice, this list of conditions and the following disclaimer.
10# * Redistributions in binary form must reproduce the above copyright
11# notice, this list of conditions and the following disclaimer in the
12# documentation and/or other materials provided with the distribution.
13#
14# THIS SOFTWARE IS PROVIDED BY THE AUTHORS ''AS IS'' AND ANY
15# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
16# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
17# DISCLAIMED. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY
18# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
19# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
20# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
21# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
23# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25import logging, sys
26
27from numpy import arange, append, array, dot, isfinite, hstack, float32, matrix, vdot, vstack, zeros
28from PyQt4.QtCore import Qt, QSettings, QVariant, SIGNAL
29from PyQt4.QtOpenGL import QGLWidget
30
31try:
32 from OpenGL.GL import *
33 from OpenGL.GLU import gluPerspective
34except ImportError, e:
35 print e
36 raise Exception('You should install PyOpenGL')
37
38try:
39 from OpenGL.GL.shaders import *
40except ImportError:
41 pass
42
43log = logging.getLogger(__name__)
44
45# HACK:
46from numpy import *
47
48# TODO: Cleanup
49def box_side(w=1.0, z=0.0):
50 return [[0.0, 0.0, z], [w, 0.0, z], [w, 0.0, z], [w, w, z],
51 [w, w, z], [0.0, w, z], [0.0, w, z], [0.0, 0.0, z]]
52
53box = [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0],
54 [1.0, 0.0, 0.0], [1.0, 0.0, 1.0],
55 [1.0, 1.0, 0.0], [1.0, 1.0, 1.0],
56 [0.0, 1.0, 0.0], [0.0, 1.0, 1.0]]
57box.extend(box_side())
58box.extend(box_side(z=1.0))
59
60class VolumePanelGL(QGLWidget):
61 '''
62 Widget for drawing two spirals.
63 '''
64
65 # Triangulation indeces
66 tri_ix = [0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 5]
67
68 # We need a Nx4 array because the modelview matrix is 4x4
69 box_array = hstack((array(box, dtype=float32), zeros((len(box), 1), dtype=float32) + 1.0))
70
71 def __init__(self, parent, core):
72 super(VolumePanelGL, self).__init__(parent)
73 self.setMinimumSize(500, 500)
74
75 self.core = core
76 self.settings = QSettings()
77
78 self.fragment_shader_src = '''
79 void main(void)
80 {
81 gl_FragColor = vec4(1.0, 0.0, 0.0, 0.0);
82 }
83 '''
84
85 self.vertex_shader_src = '''
86 uniform mat4 OldModelView;
87 varying vec4 ecPosition;
88 varying vec3 ecPosition3;
89
90 void main(void)
91 {
92 gl_TexCoord[0] = OldModelView*gl_Vertex;
93 ecPosition = gl_ModelViewMatrix * gl_Vertex;
94 ecPosition3 = vec3(ecPosition) / ecPosition.w;
95 gl_Position = ftransform();
96 }
97 '''
98
99 self.fragment_src_file = 'src/shaders/default.glf'
100 self.vertex_src_file = 'src/shaders/default.glv'
101
102 # Flags and Settings
103 self.lighting = self.settings.value('volume_rendering/lighting', QVariant(False)).toBool()
104
105 self.plane_count = 512
106
107 self.rotation = [0.0, 0.0]
108
109 # Is OpenGL initialized?
110 self.init = False
111 # Have the slices been created
112 self.init_slices = False
113
114 self.program = None
115 self.transfer_functions = None
116 self.volume_textures = None
117
118 # Connect signals
119 self.connect(core, SIGNAL('imageLoaded'), self.LoadVolumeData)
120 self.connect(core, SIGNAL('imageUnloaded'), self.LoadVolumeData)
121 self.connect(core, SIGNAL('transferFunctionUpdate'), self.UpdateTransferFunctions)
122
123 self.setFocusPolicy(Qt.ClickFocus)
124
125 def paintGL(self):
126 '''
127 Drawing routine
128 '''
129
130 if not self.init:
131 self.initializeGL()
132 self.init = True
133
134 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
135 glLoadIdentity()
136 glEnable(GL_BLEND)
137 glEnable(GL_POLYGON_SMOOTH)
138
139 glPushMatrix()
140 glTranslate(0.0, 0.0, -2.0)
141 glRotate(self.rotation[0], 1.0, 0.0, 0.0)
142 glRotate(self.rotation[1], 0.0, 1.0, 0.0)
143 glTranslate(-0.5, -0.5, -0.5)
144
145 # Draw the box
146 glUseProgram(0)
147 glColor(0.0, 1.0, 0.0)
148 glDisable(GL_LIGHTING)
149 glVertexPointerf(self.box_array[:, :3])
150 glDrawArrays(GL_LINES, 0, self.box_array.shape[0])
151
152 model_view = matrix(glGetFloatv(GL_MODELVIEW_MATRIX))
153
154 glPopMatrix()
155
156 if not self.core.images:
157 return
158
159# # Adjust light position
160# light_pos = array(model_view.T*array([0.5, 0.5, -1.0, 0.0]).reshape(-1, 1)).flatten()
161# light_pos = light_pos.tolist()[3:]
162# glLight(GL_LIGHT0, GL_POSITION, light_pos)
163
164 # Draw the slice planes
165 glUseProgram(self.program)
166
167 # We need the inverse of the old model_view matrix to generate
168 # texture coordinates
169 glUniformMatrix4fv(glGetUniformLocation(self.program, "OldModelView"),
170 1, True, model_view.getI())
171 self.SetupUniforms()
172
173 if not self.init_slices:
174 self.planes = self.BuildViewAlignedPlanes(model_view)
175
176 # Draw plane intersections
177 glColor(1.0, 1.0, 0.0)
178 glVertexPointerf(self.planes)
179 glDrawArrays(GL_TRIANGLES, 0, self.planes.shape[0])
180
181 glBindBuffer(GL_ARRAY_BUFFER, 0)
182
183 def resizeGL(self, w, h):
184 '''
185 Resize the GL window
186 '''
187
188 glViewport(0, 0, w, h)
189 self.initializeGL()
190
191 def initializeGL(self):
192 '''
193 Initialize GL
194 '''
195
196 glEnable(GL_DEPTH_TEST)
197 glEnable(GL_BLEND)
198 glEnableClientState(GL_VERTEX_ARRAY)
199
200 glDepthFunc(GL_LESS)
201 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
202 glShadeModel(GL_SMOOTH)
203
204 glClearColor(0.0, 0.0, 0.0, 1.0)
205 glClearDepth(1.0)
206
207 glMatrixMode(GL_PROJECTION)
208 glLoadIdentity()
209 gluPerspective(45.0, self.width()/float(self.height()), 0.1, 1000.0)
210
211 glMatrixMode(GL_MODELVIEW)
212 glLoadIdentity()
213
214 self.SetupLighting()
215 self.LoadVolumeData()
216
217 self.init = True
218
219 def keyPressEvent(self, event):
220 ''' Key press event '''
221
222 if event.text() == 'r':
223 self.compile_program()
224 self.update()
225
226 def mousePressEvent(self, event):
227 self.prev_mouse = (event.y(), event.x())
228
229 def mouseMoveEvent(self, event):
230
231 # TODO: Actually calculate these
232 radius = 3.0
233 screen_width = 180
234 screen_height = 180
235
236 dimensions = (screen_height, screen_width)
237 pos = (event.y(), event.x())
238
239 for i in range(2):
240 self.rotation[i] += math.asin((pos[i] - self.prev_mouse[i])/float(dimensions[i]))*180/3.14
241 self.rotation[i] %= 360.0
242
243 self.prev_mouse = pos
244
245 self.init_slices = False
246
247 self.update()
248 event.accept()
249
250 def toggleLighting(self, enabled):
251 ''' Toggle the lighting '''
252
253 self.lighting = enabled
254
255 self.settings.setValue('volume_rendering/lighting', QVariant(self.lighting))
256 self.update()
257
258 def compile_program(self):
259 '''
260 Compile a Shader program given the vertex
261 and fragment sources
262 '''
263
264 params={'imageCount':len(self.core.images)}
265
266 # Load the Shader sources from the files
267 self.LoadShaderSources()
268
269 # Delete the previous program
270# if self.program is not None:
271# glDeleteProgram(self.program)
272
273 self.program = glCreateProgram()
274
275 shaders = []
276
277 for shader_type, src in ((GL_VERTEX_SHADER, self.vertex_shader_src % (params or {})),
278 (GL_FRAGMENT_SHADER, self.fragment_shader_src % (params or {}))):
279 shader = glCreateShader(shader_type)
280 glShaderSource(shader, src)
281 glCompileShader(shader)
282
283 shaders.append(shader)
284
285 status = glGetShaderiv(shader, GL_COMPILE_STATUS)
286
287 if not status:
288 if glGetShaderiv(shader, GL_INFO_LOG_LENGTH) > 0:
289 log = glGetShaderInfoLog(shader)
290 print >> sys.stderr, log
291 glDeleteShader(shader)
292 raise ValueError, 'Shader compilation failed'
293
294 glAttachShader(self.program, shader)
295
296 glLinkProgram(self.program)
297
298 for shader in shaders:
299 glDeleteShader(shader)
300
301 def SetupLighting(self):
302 '''
303 Initialize default lighting
304 '''
305
306 glLight(GL_LIGHT0, GL_AMBIENT, (0.4, 0.4, 0.4))
307 glLight(GL_LIGHT0, GL_DIFFUSE, (0.7, 0.7, 0.7))
308 glLight(GL_LIGHT0, GL_SPECULAR, (0.2, 0.2, 0.2))
309
310 glLight(GL_LIGHT0, GL_POSITION, (-1.0, -1.0, -1.0))
311 glEnable(GL_LIGHT0)
312
313 def SetupUniforms(self):
314 set_uniform_i = lambda k, v: glUniform1i(glGetUniformLocation(self.program, k), v)
315 set_uniform_f = lambda k, v: glUniform1f(glGetUniformLocation(self.program, k), v)
316
317 image_count = len(self.core.images)
318
319 data_scales = []
320
321 # Texture coordinate transformation matrix
322 tex_transform = []
323
324 # Init the textures
325 for i in range(image_count):
326 glActiveTexture(GL_TEXTURE0 + 2*i)
327 glBindTexture(GL_TEXTURE_1D, self.transfer_functions[i])
328
329 glActiveTexture(GL_TEXTURE0 + 2*i + 1)
330 glBindTexture(GL_TEXTURE_3D, self.volume_textures[i])
331 image = self.core.images[i]
332
333 scale = list(image.image_size/float(image.image_size.max()))
334
335 data_scales.append((scale[2], scale[1], scale[0]))
336
337# print 'Image stuff' , image.ijk_to_world
338# print 'Image scale' , image.ijk_to_world / [[2], [1], [1], [1]]
339
340 k = image.ijk_to_world.flatten()
341 tex_transform.append(k)
342
343
344 glUniform1iv(glGetUniformLocation(self.program, 'TransferFunctions[0]'),
345 image_count, range(0, len(self.transfer_functions) + 1, 2))
346 glUniform1iv(glGetUniformLocation(self.program, 'Volumes[0]'),
347 image_count, range(1, len(self.transfer_functions) + 2, 2))
348 glUniform3fv(glGetUniformLocation(self.program, 'DataScales[0]'),
349 image_count, data_scales)
350
351 glUniformMatrix4fv(glGetUniformLocation(self.program, 'TexTransform[0]'),
352 image_count, False, require(tex_transform, float32, 'C'))
353
354 # Flags
355 set_uniform_i('EnableLighting', self.lighting)
356 set_uniform_i('VolumeCount', image_count)
357 set_uniform_f('CellSize', 1/256.0)
358
359 def BuildViewAlignedPlanes(self, model_view):
360 '''Generate the geometry for the
361 view aligned slice planes
362 '''
363
364 box_array = self.box_array
365
366 # Get bounding box in eye coordinates
367 points = array(model_view.T*box_array.T)
368
369 # Slice plane normal
370 N_plane = array([0.0, 0.0, 1.0], dtype=float32).reshape(3, 1)
371
372 z_min = points[2,].min()
373 z_max = points[2,].max()
374
375 # Generate Z coordinates for slice planes
376 samples_z = arange(z_min, z_max, (z_max - z_min)/self.plane_count).reshape(-1, 1)
377 # Points on slice plane centers
378 plane_pts = hstack((zeros(self.plane_count*2).reshape(-1, 2), samples_z))
379
380 # Storage for the triangles
381 triangles = array([], dtype=float32)
382
383 # Test intersections
384 for plane_pt in plane_pts:
385 plane_pt = plane_pt.reshape(-1, 1)
386
387 ints = []
388
389 p1 = points[:3, range(0, box_array.shape[0], 2)]
390 p2 = points[:3, arange(0, box_array.shape[0], 2) + 1]
391
392 u = (dot(N_plane.T, plane_pt - p1)/dot(N_plane.T, p2 - p1))
393
394 # Make sure we didn't divide by zero and that the intersection points
395 # are within U
396 valid_us = vstack((isfinite(u), (u <= 1.0), (u >= 0.0))).all(0)
397
398 p1 = p1[ : , valid_us]
399 p2 = p2[ : , valid_us]
400 u = u[0, valid_us]
401
402 ints = (p1 + u*(p2 - p1)).T
403
404 # Center of polygon
405 center = ints.mean(axis=0)
406
407 ints_c = ints - center
408
409 x = ints_c[:,0]
410 y = ints_c[:,1]
411
412 # Get the pseudo-angles
413 angles = y/(abs(x) + abs(y))
414 x_lt_0 = x < 0.0
415 x_gt_0 = x >= 0.0
416
417 angles[x_lt_0] = 2.0 - angles[x_lt_0]
418 angles[x_gt_0] = 4.0 + angles[x_gt_0]
419
420 # Sort based on the angles
421 ints = ints[angles.argsort()]
422
423 # Triangulate
424 try:
425 # TODO: Is this appending really bad?
426 data = ints[self.tri_ix[ : 3*(len(ints) - 2)]]
427 triangles = append(triangles, data)
428 except IndexError:
429 print 'Unexpected number of intersection points', len(ints)
430
431 self.init_slices = True
432
433 return triangles.reshape((-1, 3))
434
435 def _init_transfer_functions(self, image):
436 return image.t_function.get_map().flatten()
437
438 def LoadTransferFunctions(self):
439
440 if not self.core.images:
441 return
442
443 # Create Texture
444 # Free the previous transfer functions
445 glDeleteTextures(self.transfer_functions)
446
447 #TODO: Free These textures
448 self.transfer_functions = glGenTextures(len(self.core.images))
449
450 if isinstance(self.transfer_functions, long):
451 self.transfer_functions = [self.transfer_functions]
452
453 for i in range(len(self.core.images)):
454 image = self.core.images[i]
455
456 glBindTexture(GL_TEXTURE_1D, self.transfer_functions[i])
457 data = self._init_transfer_functions(image)
458
459 glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
460 glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP)
461 glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
462 glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
463 glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
464 glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
465 glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA, 256, 0, GL_RGBA, GL_FLOAT, data)
466
467 return
468
469 def UpdateTransferFunctions(self):
470
471 for i in range(len(self.core.images)):
472 image = self.core.images[i]
473 data = self._init_transfer_functions(image)
474 glBindTexture(GL_TEXTURE_1D, self.transfer_functions[i])
475 glTexSubImage1D(GL_TEXTURE_1D, 0, 0, 256, GL_RGBA, GL_FLOAT, data)
476
477 self.update()
478
479 def LoadVolumeData(self):
480
481 if not (self.core.images and self.init):
482 self.update()
483 return
484
485 # These have to be loaded everytime a new image is loaded/unloaded
486 self.LoadTransferFunctions()
487
488 # Recompile the program because it requires the image count
489 self.compile_program()
490
491 # Free the previous volume textures
492 glDeleteTextures(self.volume_textures)
493
494 # Create some textures
495 self.volume_textures = glGenTextures(len(self.core.images))
496
497 if isinstance(self.volume_textures, long):
498 self.volume_textures = [self.volume_textures]
499
500 for image, texture in zip(self.core.images, self.volume_textures):
501 # Create Texture
502 glPixelStorei(GL_UNPACK_ALIGNMENT,1)
503 glBindTexture(GL_TEXTURE_3D, texture)
504
505 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP)
506 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP)
507 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP)
508
509 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
510 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
511
512 glTexImage3D(GL_TEXTURE_3D, 0, GL_ALPHA,
513 image.data.shape[2], image.data.shape[1], image.data.shape[0],
514 0, GL_ALPHA, GL_UNSIGNED_BYTE, image.data)
515
516 self.update()
517 return
518
519 def LoadShaderSources(self):
520 try:
521 self.fragment_shader_src = open(self.fragment_src_file).read()
522 except IOError, e:
523 print 'Fragment source not found, using default'
524
525 try:
526 self.vertex_shader_src = open(self.vertex_src_file).read()
527 except IOError, e:
528 print 'Vertex source not found, using default'
529

Keys:

r - Reload shader sources

Warning

This needs to be finished

This code uses view-aligned slices