Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
vtkMitkOpenGLVolumeTextureMapper3D.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #ifdef _OPENMP
18 #include <omp.h>
19 #endif
20 
21 #include "vtkWindows.h"
22 #include "mitkCommon.h"
24 
25 #define GPU_INFO MITK_INFO("mapper.vr")
26 #define GPU_WARN MITK_WARN("mapper.vr")
27 
28 #include "vtkCamera.h"
29 #include "vtkDataArray.h"
30 #include "vtkImageData.h"
31 #include "vtkLight.h"
32 #include "vtkLightCollection.h"
33 #include "vtkMath.h"
34 #include "vtkMatrix4x4.h"
35 #include "vtkObjectFactory.h"
36 #include "vtkOpenGLExtensionManager.h"
37 #include "vtkPlane.h"
38 #include "vtkPlaneCollection.h"
39 #include "vtkPointData.h"
40 #include "vtkRenderWindow.h"
41 #include "vtkRenderer.h"
42 #include "vtkTimerLog.h"
43 #include "vtkTransform.h"
44 #include "vtkVolumeProperty.h"
45 #include "vtkgl.h"
46 
47 #include "vtkOpenGLRenderWindow.h"
48 
49 #define myGL_COMPRESSED_RGB_S3TC_DXT1_EXT 0x83F0
50 #define myGL_COMPRESSED_LUMINANCE_ALPHA_LATC2_EXT 0x8C72
51 #define myGL_COMPRESSED_RGBA_S3TC_DXT5_EXT 0x83F3
52 
54 
55  //# We need some temporary variables
56  "TEMP index, normal, finalColor;\n"
57  "TEMP temp,temp1, temp2, temp3,temp4; \n"
58  "TEMP sampleColor;\n"
59  "TEMP ndotl, ndoth, ndotv; \n"
60  "TEMP lightInfo, lightResult;\n"
61 
62  //# We are going to use the first
63  //# texture coordinate
64  "ATTRIB tex0 = fragment.texcoord[0];\n"
65 
66  //# This is the lighting information
67  "PARAM lightDirection = program.local[0];\n"
68  "PARAM halfwayVector = program.local[1];\n"
69  "PARAM coefficient = program.local[2];\n"
70  "PARAM lightDiffColor = program.local[3]; \n"
71  "PARAM lightSpecColor = program.local[4]; \n"
72  "PARAM viewVector = program.local[5];\n"
73  "PARAM constants = program.local[6];\n"
74 
75  //# This is our output color
76  "OUTPUT out = result.color;\n"
77 
78  //# Look up the gradient direction
79  //# in the third volume
80  "TEX temp2, tex0, texture[0], 3D;\n"
81 
82  //# This normal is stored 0 to 1, change to -1 to 1
83  //# by multiplying by 2.0 then adding -1.0.
84  "MAD normal, temp2, constants.x, constants.y;\n"
85 
86  "DP3 temp4, normal, normal;\n"
87  "RSQ temp, temp4.x;\n"
88  "MUL normal, normal, temp;\n"
89 
90  //"RCP temp4,temp.x;\n"
91 
92  //"MUL temp2.w,temp2.w,temp4.x;\n"
93 
94  //"MUL_SAT temp2.w,temp2.w,6.0;\n"
95 
96  "TEX sampleColor, tex0, texture[1], 3D;\n"
97 
98  //# Take the dot product of the light
99  //# direction and the normal
100  "DP3 ndotl, normal, lightDirection;\n"
101 
102  //# Take the dot product of the halfway
103  //# vector and the normal
104  "DP3 ndoth, normal, halfwayVector;\n"
105 
106  "DP3 ndotv, normal, viewVector;\n"
107 
108  //# flip if necessary for two sided lighting
109  "MUL temp3, ndotl, constants.y; \n"
110  "CMP ndotl, ndotv, ndotl, temp3;\n"
111  "MUL temp3, ndoth, constants.y; \n"
112  "CMP ndoth, ndotv, ndoth, temp3;\n"
113 
114  //# put the pieces together for a LIT operation
115  "MOV lightInfo.x, ndotl.x; \n"
116  "MOV lightInfo.y, ndoth.x; \n"
117  "MOV lightInfo.w, coefficient.w; \n"
118 
119  //# compute the lighting
120  "LIT lightResult, lightInfo;\n"
121 
122  //# COLOR FIX
123  "MUL lightResult, lightResult, 4.0;\n"
124 
125  //# This is the ambient contribution
126  "MUL finalColor, coefficient.x, sampleColor;\n"
127 
128  //# This is the diffuse contribution
129  "MUL temp3, lightDiffColor, sampleColor;\n"
130  "MUL temp3, temp3, lightResult.y;\n"
131  "ADD finalColor, finalColor, temp3;\n"
132 
133  //# This is th specular contribution
134  "MUL temp3, lightSpecColor, lightResult.z; \n"
135 
136  //# Add specular into result so far, and replace
137  //# with the original alpha.
138  "ADD out, finalColor, temp3;\n"
139  "MOV out.w, temp2.w;\n"
140 
141  "END\n";
142 
144 
145  //# This is the fragment program for one
146  //# component data with shading
147 
148  //# We need some temporary variables
149  "TEMP index, normal, finalColor;\n"
150  "TEMP temp,temp1, temp2, temp3,temp4; \n"
151  "TEMP sampleColor;\n"
152  "TEMP ndotl, ndoth, ndotv; \n"
153  "TEMP lightInfo, lightResult;\n"
154 
155  //# We are going to use the first
156  //# texture coordinate
157  "ATTRIB tex0 = fragment.texcoord[0];\n"
158 
159  //# This is the lighting information
160  "PARAM lightDirection = program.local[0];\n"
161  "PARAM halfwayVector = program.local[1];\n"
162  "PARAM coefficient = program.local[2];\n"
163  "PARAM lightDiffColor = program.local[3]; \n"
164  "PARAM lightSpecColor = program.local[4]; \n"
165  "PARAM viewVector = program.local[5];\n"
166  "PARAM constants = program.local[6];\n"
167 
168  //# This is our output color
169  "OUTPUT out = result.color;\n"
170 
171  //# Look up the gradient direction
172  //# in the third volume
173  "TEX temp2, tex0, texture[0], 3D;\n"
174 
175  // Gradient Compution
176 
177  //# Look up the scalar value / gradient
178  //# magnitude in the first volume
179  //"TEX temp1, tex0, texture[0], 3D;\n"
180 
181  /*
182 
183  "ADD temp3,tex0,{-0.005,0,0};\n"
184  "TEX temp2,temp3, texture[0], 3D;\n"
185 
186  //"ADD temp3,tex0,{ 0.005,0,0};\n"
187  //"TEX temp1,temp3, texture[0], 3D;\n"
188 
189  "SUB normal.x,temp2.y,temp1.y;\n"
190 
191  "ADD temp3,tex0,{0,-0.005,0};\n"
192  "TEX temp2,temp3, texture[0], 3D;\n"
193 
194  //"ADD temp3,tex0,{0, 0.005,0};\n"
195  //"TEX temp1,temp3, texture[0], 3D;\n"
196 
197  "SUB normal.y,temp2.y,temp1.y;\n"
198 
199  "ADD temp3,tex0,{0,0,-0.005};\n"
200  "TEX temp2,temp3, texture[0], 3D;\n"
201 
202  //"ADD temp3,tex0,{0,0, 0.005};\n"
203  //"TEX temp1,temp3, texture[0], 3D;\n"
204 
205  "SUB normal.z,temp2.y,temp1.y;\n"
206 
207  */
208 
209  //"MOV normal,{1,1,1};\n"
210 
211  "MOV index.x,temp2.a;\n"
212 
213  //# This normal is stored 0 to 1, change to -1 to 1
214  //# by multiplying by 2.0 then adding -1.0.
215  "MAD normal, temp2, constants.x, constants.y;\n"
216 
217  //# Swizzle this to use (a,r) as texture
218  //# coordinates
219  //"SWZ index, temp1, a, r, 1, 1;\n"
220 
221  //# Use this coordinate to look up a
222  //# final color in the third texture
223  //# (this is a 2D texture)
224 
225  "DP3 temp4, normal, normal;\n"
226 
227  "RSQ temp, temp4.x;\n"
228 
229  "RCP temp4,temp.x;\n"
230 
231  "MUL normal, normal, temp;\n"
232 
233  "MOV index.y, temp4.x;\n"
234 
235  "TEX sampleColor, index, texture[1], 2D;\n"
236 
237  //"MUL sampleColor.w,sampleColor.w,temp4.x;\n"
238 
239  //# Take the dot product of the light
240  //# direction and the normal
241  "DP3 ndotl, normal, lightDirection;\n"
242 
243  //# Take the dot product of the halfway
244  //# vector and the normal
245  "DP3 ndoth, normal, halfwayVector;\n"
246 
247  "DP3 ndotv, normal, viewVector;\n"
248 
249  //# flip if necessary for two sided lighting
250  "MUL temp3, ndotl, constants.y; \n"
251  "CMP ndotl, ndotv, ndotl, temp3;\n"
252  "MUL temp3, ndoth, constants.y; \n"
253  "CMP ndoth, ndotv, ndoth, temp3;\n"
254 
255  //# put the pieces together for a LIT operation
256  "MOV lightInfo.x, ndotl.x; \n"
257  "MOV lightInfo.y, ndoth.x; \n"
258  "MOV lightInfo.w, coefficient.w; \n"
259 
260  //# compute the lighting
261  "LIT lightResult, lightInfo;\n"
262 
263  //# COLOR FIX
264  "MUL lightResult, lightResult, 4.0;\n"
265 
266  //# This is the ambient contribution
267  "MUL finalColor, coefficient.x, sampleColor;\n"
268 
269  //# This is the diffuse contribution
270  "MUL temp3, lightDiffColor, sampleColor;\n"
271  "MUL temp3, temp3, lightResult.y;\n"
272  "ADD finalColor, finalColor, temp3;\n"
273 
274  //# This is th specular contribution
275  "MUL temp3, lightSpecColor, lightResult.z; \n"
276 
277  //# Add specular into result so far, and replace
278  //# with the original alpha.
279  "ADD out, finalColor, temp3;\n"
280  "MOV out.w, sampleColor.w;\n"
281 
282  "END\n";
283 
284 //#ifndef VTK_IMPLEMENT_MESA_CXX
286 //#endif
287 
289 {
290  // GPU_INFO << "vtkMitkOpenGLVolumeTextureMapper3D";
291 
292  this->Initialized = 0;
293  this->Volume1Index = 0;
294  this->Volume2Index = 0;
295  this->Volume3Index = 0;
296  this->ColorLookupIndex = 0;
297  this->AlphaLookupIndex = 0;
298  this->RenderWindow = NULL;
299  this->SupportsCompressedTexture = false;
300 
302  prgRGBAShade = 0;
303 }
304 
306 {
307  // GPU_INFO << "~vtkMitkOpenGLVolumeTextureMapper3D";
309  vtkgl::DeleteProgramsARB(1, &prgOneComponentShade);
310 
311  if (prgRGBAShade)
312  vtkgl::DeleteProgramsARB(1, &prgRGBAShade);
313 }
314 
315 // Release the graphics resources used by this texture.
317 {
318  // GPU_INFO << "ReleaseGraphicsResources";
319 
320  if ((this->Volume1Index || this->Volume2Index || this->Volume3Index || this->ColorLookupIndex) && renWin)
321  {
322  static_cast<vtkRenderWindow *>(renWin)->MakeCurrent();
323 #ifdef GL_VERSION_1_1
324  // free any textures
325  this->DeleteTextureIndex(&this->Volume1Index);
326  this->DeleteTextureIndex(&this->Volume2Index);
327  this->DeleteTextureIndex(&this->Volume3Index);
328  this->DeleteTextureIndex(&this->ColorLookupIndex);
329  this->DeleteTextureIndex(&this->AlphaLookupIndex);
330 #endif
331  }
332  this->Volume1Index = 0;
333  this->Volume2Index = 0;
334  this->Volume3Index = 0;
335  this->ColorLookupIndex = 0;
336  this->RenderWindow = NULL;
337  this->SupportsCompressedTexture = false;
338  this->SupportsNonPowerOfTwoTextures = false;
339 
340  this->Modified();
341 }
342 
343 // Release the graphics resources used by this texture.
345 {
346  // GPU_INFO << "ReleaseGraphicsResources";
347 
348  vtkWindow *renWin = renderer->GetVtkRenderer()->GetRenderWindow();
349 
350  if ((this->Volume1Index || this->Volume2Index || this->Volume3Index || this->ColorLookupIndex) && renWin)
351  {
352  static_cast<vtkRenderWindow *>(renWin)->MakeCurrent();
353 #ifdef GL_VERSION_1_1
354  // free any textures
355  this->DeleteTextureIndex(&this->Volume1Index);
356  this->DeleteTextureIndex(&this->Volume2Index);
357  this->DeleteTextureIndex(&this->Volume3Index);
358  this->DeleteTextureIndex(&this->ColorLookupIndex);
359  this->DeleteTextureIndex(&this->AlphaLookupIndex);
360 #endif
361  }
362  this->Volume1Index = 0;
363  this->Volume2Index = 0;
364  this->Volume3Index = 0;
365  this->ColorLookupIndex = 0;
366  this->RenderWindow = NULL;
367  this->SupportsCompressedTexture = false;
368  this->SupportsNonPowerOfTwoTextures = false;
369 
370  this->Modified();
371 }
372 
373 void vtkMitkOpenGLVolumeTextureMapper3D::Render(vtkRenderer *ren, vtkVolume *vol)
374 {
375  // GPU_INFO << "Render";
376 
377  ren->GetRenderWindow()->MakeCurrent();
378 
379  if (!this->Initialized)
380  {
381  // this->Initialize();
382  this->Initialize(ren);
383  }
384 
385  if (!this->RenderPossible)
386  {
387  vtkErrorMacro("required extensions not supported");
388  return;
389  }
390 
391  vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
392  vtkPlaneCollection *clipPlanes;
393  vtkPlane *plane;
394  int numClipPlanes = 0;
395  double planeEquation[4];
396 
397  // build transformation
398  vol->GetMatrix(matrix);
399  matrix->Transpose();
400 
401  glPushAttrib(GL_ENABLE_BIT | GL_COLOR_BUFFER_BIT | GL_STENCIL_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_POLYGON_BIT |
402  GL_TEXTURE_BIT);
403 
404  int i;
405 
406  // Use the OpenGL clip planes
407  clipPlanes = this->ClippingPlanes;
408  if (clipPlanes)
409  {
410  numClipPlanes = clipPlanes->GetNumberOfItems();
411  if (numClipPlanes > 6)
412  {
413  vtkErrorMacro(<< "OpenGL guarantees only 6 additional clipping planes");
414  }
415 
416  for (i = 0; i < numClipPlanes; i++)
417  {
418  glEnable(static_cast<GLenum>(GL_CLIP_PLANE0 + i));
419 
420  plane = static_cast<vtkPlane *>(clipPlanes->GetItemAsObject(i));
421 
422  planeEquation[0] = plane->GetNormal()[0];
423  planeEquation[1] = plane->GetNormal()[1];
424  planeEquation[2] = plane->GetNormal()[2];
425  planeEquation[3] = -(planeEquation[0] * plane->GetOrigin()[0] + planeEquation[1] * plane->GetOrigin()[1] +
426  planeEquation[2] * plane->GetOrigin()[2]);
427  glClipPlane(static_cast<GLenum>(GL_CLIP_PLANE0 + i), planeEquation);
428  }
429  }
430 
431  // insert model transformation
432  glMatrixMode(GL_MODELVIEW);
433  glPushMatrix();
434  glMultMatrixd(matrix->Element[0]);
435 
436  glColor4f(1.0, 1.0, 1.0, 1.0);
437 
438  // Turn lighting off - the polygon textures already have illumination
439  glDisable(GL_LIGHTING);
440 
441  // vtkGraphicErrorMacro(ren->GetRenderWindow(),"Before actual render method");
442 
443  this->RenderFP(ren, vol);
444 
445  // pop transformation matrix
446  glMatrixMode(GL_MODELVIEW);
447  glPopMatrix();
448 
449  matrix->Delete();
450  glPopAttrib();
451 }
452 
453 void vtkMitkOpenGLVolumeTextureMapper3D::RenderFP(vtkRenderer *ren, vtkVolume *vol)
454 {
455  // GPU_INFO << "RenderFP";
456 
457  /*
458  glAlphaFunc (GL_GREATER, static_cast<GLclampf>(1.0/255.0));
459  glEnable (GL_ALPHA_TEST);
460  */
461 
462  glEnable(GL_BLEND);
463  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
464 
465  int components = this->GetInput()->GetNumberOfScalarComponents();
466  switch (components)
467  {
468  case 1:
469  this->RenderOneIndependentShadeFP(ren, vol);
470  break;
471 
472  case 4:
473  this->RenderRGBAShadeFP(ren, vol);
474  break;
475  }
476 
477  vtkgl::ActiveTexture(vtkgl::TEXTURE2);
478  glDisable(GL_TEXTURE_2D);
479  glDisable(vtkgl::TEXTURE_3D);
480 
481  vtkgl::ActiveTexture(vtkgl::TEXTURE1);
482  glDisable(GL_TEXTURE_2D);
483  glDisable(vtkgl::TEXTURE_3D);
484 
485  vtkgl::ActiveTexture(vtkgl::TEXTURE0);
486  glDisable(GL_TEXTURE_2D);
487  glDisable(vtkgl::TEXTURE_3D);
488 
489  glDisable(GL_BLEND);
490 }
491 
493 {
494  // GPU_INFO << "DeleteTextureIndex";
495 
496  if (glIsTexture(*index))
497  {
498  GLuint tempIndex;
499  tempIndex = *index;
500  glDeleteTextures(1, &tempIndex);
501  *index = 0;
502  }
503 }
504 
506 {
507  // GPU_INFO << "CreateTextureIndex";
508 
509  GLuint tempIndex = 0;
510  glGenTextures(1, &tempIndex);
511  *index = static_cast<long>(tempIndex);
512 }
513 
514 void vtkMitkOpenGLVolumeTextureMapper3D::RenderPolygons(vtkRenderer *ren, vtkVolume *vol, int stages[4])
515 {
516  // GPU_INFO << "RenderPolygons";
517 
518  vtkRenderWindow *renWin = ren->GetRenderWindow();
519 
520  if (renWin->CheckAbortStatus())
521  {
522  return;
523  }
524 
525  double bounds[27][6];
526  float distance2[27];
527 
528  int numIterations;
529  int i, j, k;
530 
531  // No cropping case - render the whole thing
532  if (!this->Cropping)
533  {
534  // Use the input data bounds - we'll take care of the volume's
535  // matrix during rendering
536  this->GetInput()->GetBounds(bounds[0]);
537  numIterations = 1;
538  }
539  // Simple cropping case - render the subvolume
540  else if (this->CroppingRegionFlags == 0x2000)
541  {
542  this->GetCroppingRegionPlanes(bounds[0]);
543  numIterations = 1;
544  }
545  // Complex cropping case - render each region in back-to-front order
546  else
547  {
548  // Get the camera position
549  double camPos[4];
550  ren->GetActiveCamera()->GetPosition(camPos);
551 
552  double volBounds[6];
553  this->GetInput()->GetBounds(volBounds);
554 
555  // Pass camera through inverse volume matrix
556  // so that we are in the same coordinate system
557  vtkMatrix4x4 *volMatrix = vtkMatrix4x4::New();
558  vol->GetMatrix(volMatrix);
559  camPos[3] = 1.0;
560  volMatrix->Invert();
561  volMatrix->MultiplyPoint(camPos, camPos);
562  volMatrix->Delete();
563  if (camPos[3])
564  {
565  camPos[0] /= camPos[3];
566  camPos[1] /= camPos[3];
567  camPos[2] /= camPos[3];
568  }
569 
570  // These are the region limits for x (first four), y (next four) and
571  // z (last four). The first region limit is the lower bound for
572  // that axis, the next two are the region planes along that axis, and
573  // the final one in the upper bound for that axis.
574  float limit[12];
575  for (i = 0; i < 3; i++)
576  {
577  limit[i * 4] = volBounds[i * 2];
578  limit[i * 4 + 1] = this->CroppingRegionPlanes[i * 2];
579  limit[i * 4 + 2] = this->CroppingRegionPlanes[i * 2 + 1];
580  limit[i * 4 + 3] = volBounds[i * 2 + 1];
581  }
582 
583  // For each of the 27 possible regions, find out if it is enabled,
584  // and if so, compute the bounds and the distance from the camera
585  // to the center of the region.
586  int numRegions = 0;
587  int region;
588  for (region = 0; region < 27; region++)
589  {
590  int regionFlag = 1 << region;
591 
592  if (this->CroppingRegionFlags & regionFlag)
593  {
594  // what is the coordinate in the 3x3x3 grid
595  int loc[3];
596  loc[0] = region % 3;
597  loc[1] = (region / 3) % 3;
598  loc[2] = (region / 9) % 3;
599 
600  // compute the bounds and center
601  float center[3];
602  for (i = 0; i < 3; i++)
603  {
604  bounds[numRegions][i * 2] = limit[4 * i + loc[i]];
605  bounds[numRegions][i * 2 + 1] = limit[4 * i + loc[i] + 1];
606  center[i] = (bounds[numRegions][i * 2] + bounds[numRegions][i * 2 + 1]) / 2.0;
607  }
608 
609  // compute the distance squared to the center
610  distance2[numRegions] = (camPos[0] - center[0]) * (camPos[0] - center[0]) +
611  (camPos[1] - center[1]) * (camPos[1] - center[1]) +
612  (camPos[2] - center[2]) * (camPos[2] - center[2]);
613 
614  // we've added one region
615  numRegions++;
616  }
617  }
618 
619  // Do a quick bubble sort on distance
620  for (i = 1; i < numRegions; i++)
621  {
622  for (j = i; j > 0 && distance2[j] > distance2[j - 1]; j--)
623  {
624  float tmpBounds[6];
625  float tmpDistance2;
626 
627  for (k = 0; k < 6; k++)
628  {
629  tmpBounds[k] = bounds[j][k];
630  }
631  tmpDistance2 = distance2[j];
632 
633  for (k = 0; k < 6; k++)
634  {
635  bounds[j][k] = bounds[j - 1][k];
636  }
637  distance2[j] = distance2[j - 1];
638 
639  for (k = 0; k < 6; k++)
640  {
641  bounds[j - 1][k] = tmpBounds[k];
642  }
643  distance2[j - 1] = tmpDistance2;
644  }
645  }
646 
647  numIterations = numRegions;
648  }
649 
650  // loop over all regions we need to render
651  for (int loop = 0; loop < numIterations; loop++)
652  {
653  // Compute the set of polygons for this region
654  // according to the bounds
655  this->ComputePolygons(ren, vol, bounds[loop]);
656 
657  // Loop over the polygons
658  for (i = 0; i < this->NumberOfPolygons; i++)
659  {
660  if (renWin->CheckAbortStatus())
661  {
662  return;
663  }
664 
665  float *ptr = this->PolygonBuffer + 36 * i;
666 
667  glBegin(GL_TRIANGLE_FAN);
668 
669  for (j = 0; j < 6; j++)
670  {
671  if (ptr[0] < 0.0)
672  {
673  break;
674  }
675 
676  for (k = 0; k < 4; k++)
677  {
678  if (stages[k])
679  {
680  vtkgl::MultiTexCoord3fv(vtkgl::TEXTURE0 + k, ptr);
681  }
682  }
683  glVertex3fv(ptr + 3);
684 
685  ptr += 6;
686  }
687  glEnd();
688  }
689  }
690 }
691 
692 // This method moves the scalars from the input volume into volume1 (and
693 // possibly volume2) which are the 3D texture maps used for rendering.
694 //
695 // In the case where our volume is a power of two, the copy is done
696 // directly. If we need to resample, then trilinear interpolation is used.
697 //
698 // A shift/scale is applied to the input scalar value to produce an 8 bit
699 // value for the texture volume.
700 //
701 // When the input data is one component, the scalar value is placed in the
702 // second component of the two component volume1. The first component is
703 // filled in later with the gradient magnitude.
704 //
705 // When the input data is two component non-independent, the first component
706 // of the input data is placed in the first component of volume1, and the
707 // second component of the input data is placed in the third component of
708 // volume1. Volume1 has three components - the second is filled in later with
709 // the gradient magnitude.
710 //
711 // When the input data is four component non-independent, the first three
712 // components of the input data are placed in volume1 (which has three
713 // components), and the fourth component is placed in the second component
714 // of volume2. The first component of volume2 is later filled in with the
715 // gradient magnitude.
716 
717 template <class T>
718 class ScalarGradientCompute
719 {
720  T *dataPtr;
721  unsigned char *tmpPtr;
722  unsigned char *tmpPtr2;
723  int sizeX;
724  int sizeY;
725  int sizeZ;
726  int sizeXY;
727  int sizeXm1;
728  int sizeYm1;
729  int sizeZm1;
730  int fullX;
731  int fullY;
732  int fullZ;
733  int fullXY;
734  int currentChunkStart;
735  int currentChunkEnd;
736 
737  int offZ;
738 
739  float offset;
740  float scale;
741 
742 public:
743  ScalarGradientCompute(T *_dataPtr,
744  unsigned char *_tmpPtr,
745  unsigned char *_tmpPtr2,
746  int _sizeX,
747  int _sizeY,
748  int _sizeZ,
749  int _fullX,
750  int _fullY,
751  int _fullZ,
752  float _offset,
753  float _scale)
754  {
755  dataPtr = _dataPtr;
756  tmpPtr = _tmpPtr;
757  tmpPtr2 = _tmpPtr2;
758  sizeX = _sizeX;
759  sizeY = _sizeY;
760  sizeZ = _sizeZ;
761  fullX = _fullX;
762  fullY = _fullY;
763  fullZ = _fullZ;
764  offset = _offset;
765  scale = _scale;
766 
767  sizeXY = sizeX * sizeY;
768  sizeXm1 = sizeX - 1;
769  sizeYm1 = sizeY - 1;
770  sizeZm1 = sizeZ - 1;
771 
772  fullXY = fullX * fullY;
773  }
774 
775  inline float sample(int x, int y, int z) { return float(dataPtr[x + y * sizeX + z * sizeXY]); }
776  inline void fill(int x, int y, int z)
777  {
778  int doff = x + y * fullX + (z - offZ) * fullXY;
779 
780  tmpPtr[doff * 4 + 0] = 0;
781  tmpPtr[doff * 4 + 1] = 0;
782  tmpPtr[doff * 4 + 2] = 0;
783  tmpPtr[doff * 4 + 3] = 0;
784  /*
785 tmpPtr2[doff*3+0]= 0;
786 tmpPtr2[doff*3+1]= 0;
787 tmpPtr2[doff*3+2]= 0;
788 */
789  }
790 
791  inline int clamp(int x)
792  {
793  if (x < 0)
794  x = 0;
795  else if (x > 255)
796  x = 255;
797  return x;
798  }
799 
800  inline void write(int x, int y, int z, float grayValue, float gx, float gy, float gz)
801  {
802  /*
803  gx /= aspect[0];
804  gy /= aspect[1];
805  gz /= aspect[2];
806  */
807  // Compute the gradient magnitude
808 
809  int iGrayValue = static_cast<int>((grayValue + offset) * scale + 0.5f);
810 
811  gx *= scale;
812  gy *= scale;
813  gz *= scale;
814 
815  float t = sqrtf(gx * gx + gy * gy + gz * gz);
816 
817  if (t > 0.01f)
818  {
819  if (t < 2.0f)
820  {
821  float fac = 2.0f / t;
822  gx *= fac;
823  gy *= fac;
824  gz *= fac;
825  }
826  else if (t > 255.0f)
827  {
828  float fac = 255.0f / t;
829  gx *= fac;
830  gy *= fac;
831  gz *= fac;
832  }
833  }
834  else
835  {
836  gx = gy = gz = 0.0f;
837  }
838 
839  int nx = static_cast<int>(0.5f * gx + 127.5f);
840  int ny = static_cast<int>(0.5f * gy + 127.5f);
841  int nz = static_cast<int>(0.5f * gz + 127.5f);
842 
843  int doff = x + y * fullX + (z - offZ) * fullXY;
844 
845  // tmpPtr[doff*2+0]= 0;
846 
847  tmpPtr[doff * 4 + 0] = clamp(nx);
848  tmpPtr[doff * 4 + 1] = clamp(ny);
849  tmpPtr[doff * 4 + 2] = clamp(nz);
850  tmpPtr[doff * 4 + 3] = clamp(iGrayValue);
851 
852  /*
853  if( z == fullZ/2 )
854  if( y == fullY/2 )
855  MITK_INFO << x << " " << y << " " << z << " : " << iGrayValue << " : " << iGradient;
856  */
857  }
858 
859  inline void compute(int x, int y, int z)
860  {
861  float grayValue = sample(x, y, z);
862  float gx, gy, gz;
863 
864  gx = sample(x + 1, y, z) - sample(x - 1, y, z);
865  gy = sample(x, y + 1, z) - sample(x, y - 1, z);
866  gz = sample(x, y, z + 1) - sample(x, y, z - 1);
867 
868  write(x, y, z, grayValue, gx, gy, gz);
869  }
870 
871  inline void computeClamp(int x, int y, int z)
872  {
873  float grayValue = sample(x, y, z);
874  float gx, gy, gz;
875 
876  if (x == 0)
877  gx = 2.0f * (sample(x + 1, y, z) - grayValue);
878  else if (x == sizeXm1)
879  gx = 2.0f * (grayValue - sample(x - 1, y, z));
880  else
881  gx = sample(x + 1, y, z) - sample(x - 1, y, z);
882 
883  if (y == 0)
884  gy = 2.0f * (sample(x, y + 1, z) - grayValue);
885  else if (y == sizeYm1)
886  gy = 2.0f * (grayValue - sample(x, y - 1, z));
887  else
888  gy = sample(x, y + 1, z) - sample(x, y - 1, z);
889 
890  if (z == 0)
891  gz = 2.0f * (sample(x, y, z + 1) - grayValue);
892  else if (z == sizeZm1)
893  gz = 2.0f * (grayValue - sample(x, y, z - 1));
894  else
895  gz = sample(x, y, z + 1) - sample(x, y, z - 1);
896 
897  write(x, y, z, grayValue, gx, gy, gz);
898  }
899 
900  inline void compute1D(int y, int z)
901  {
902  int x;
903 
904  x = 0;
905  computeClamp(x, y, z);
906  x++;
907 
908  while (x < sizeX - 1)
909  {
910  compute(x, y, z);
911  x++;
912  }
913 
914  if (x < sizeX)
915  {
916  computeClamp(x, y, z);
917  x++;
918  }
919 
920  while (x < fullX)
921  {
922  fill(x, y, z);
923  x++;
924  }
925  }
926 
927  inline void fill1D(int y, int z)
928  {
929  int x;
930 
931  x = 0;
932  while (x < fullX)
933  {
934  fill(x, y, z);
935  x++;
936  }
937  }
938 
939  inline void computeClamp1D(int y, int z)
940  {
941  int x;
942 
943  x = 0;
944 
945  while (x < sizeX)
946  {
947  computeClamp(x, y, z);
948  x++;
949  }
950 
951  while (x < fullX)
952  {
953  fill(x, y, z);
954  x++;
955  }
956  }
957 
958  inline void computeClamp2D(int z)
959  {
960  int y;
961 
962  y = 0;
963 
964  while (y < sizeY)
965  {
966  computeClamp1D(y, z);
967  y++;
968  }
969 
970  while (y < fullY)
971  {
972  fill1D(y, z);
973  y++;
974  }
975  }
976 
977  inline void compute2D(int z)
978  {
979  int y;
980 
981  y = 0;
982  computeClamp1D(y, z);
983  y++;
984 
985  while (y < sizeY - 1)
986  {
987  compute1D(y, z);
988  y++;
989  }
990 
991  if (y < sizeY)
992  {
993  computeClamp1D(y, z);
994  y++;
995  }
996 
997  while (y < fullY)
998  {
999  fill1D(y, z);
1000  y++;
1001  }
1002  }
1003 
1004  inline void fill2D(int z)
1005  {
1006  int y;
1007 
1008  y = 0;
1009  while (y < fullY)
1010  {
1011  fill1D(y, z);
1012  y++;
1013  }
1014  }
1015 
1016  inline void fillSlices(int currentChunkStart, int currentChunkEnd)
1017  {
1018  offZ = currentChunkStart;
1019 
1020 /*
1021  int num = omp_get_num_procs();
1022  MITK_INFO << "omp uses " << num << " processors";
1023 */
1024 
1025 #pragma omp parallel for
1026  for (int z = currentChunkStart; z <= currentChunkEnd; z++)
1027  {
1028  if (z == 0 || z == sizeZ - 1)
1029  computeClamp2D(z);
1030  else if (z >= sizeZ)
1031  fill2D(z);
1032  else
1033  compute2D(z);
1034  }
1035  }
1036 };
1037 
1038 template <class T>
1040  T *dataPtr, vtkMitkVolumeTextureMapper3D *me, float offset, float scale, GLuint volume1, GLuint /*volume2*/)
1041 {
1042  // T *inPtr;
1043  // unsigned char *outPtr, *outPtr2;
1044  // int i, j, k;
1045  // int idx;
1046 
1047  int inputDimensions[3];
1048  double inputSpacing[3];
1049  vtkImageData *input = me->GetInput();
1050 
1051  input->GetDimensions(inputDimensions);
1052  input->GetSpacing(inputSpacing);
1053 
1054  int outputDimensions[3];
1055  float outputSpacing[3];
1056  me->GetVolumeDimensions(outputDimensions);
1057  me->GetVolumeSpacing(outputSpacing);
1058 
1059  // int components = input->GetNumberOfScalarComponents();
1060 
1061  // double wx, wy, wz;
1062  // double fx, fy, fz;
1063  // int x, y, z;
1064 
1065  // double sampleRate[3];
1066  // sampleRate[0] = outputSpacing[0] / static_cast<double>(inputSpacing[0]);
1067  // sampleRate[1] = outputSpacing[1] / static_cast<double>(inputSpacing[1]);
1068  // sampleRate[2] = outputSpacing[2] / static_cast<double>(inputSpacing[2]);
1069 
1070  int fullX = outputDimensions[0];
1071  int fullY = outputDimensions[1];
1072  int fullZ = outputDimensions[2];
1073 
1074  int sizeX = inputDimensions[0];
1075  int sizeY = inputDimensions[1];
1076  int sizeZ = inputDimensions[2];
1077 
1078  int chunkSize = 64;
1079 
1080  if (fullZ < chunkSize)
1081  chunkSize = fullZ;
1082 
1083  int numChunks = (fullZ + (chunkSize - 1)) / chunkSize;
1084 
1085  // inPtr = dataPtr;
1086 
1087  unsigned char *tmpPtr = new unsigned char[fullX * fullY * chunkSize * 4];
1088  unsigned char *tmpPtr2 = 0; // new unsigned char[fullX*fullY*chunkSize*3];
1089 
1090  // For each Chunk
1091  {
1092  ScalarGradientCompute<T> sgc(dataPtr, tmpPtr, tmpPtr2, sizeX, sizeY, sizeZ, fullX, fullY, fullZ, offset, scale);
1093 
1094  int currentChunk = 0;
1095 
1096  while (currentChunk < numChunks)
1097  {
1098  int currentChunkStart = currentChunk * chunkSize;
1099  int currentChunkEnd = currentChunkStart + chunkSize - 1;
1100 
1101  if (currentChunkEnd > (fullZ - 1))
1102  currentChunkEnd = (fullZ - 1);
1103 
1104  int currentChunkSize = currentChunkEnd - currentChunkStart + 1;
1105 
1106  sgc.fillSlices(currentChunkStart, currentChunkEnd);
1107 
1108  glBindTexture(vtkgl::TEXTURE_3D, volume1);
1109  vtkgl::TexSubImage3D(vtkgl::TEXTURE_3D,
1110  0,
1111  0,
1112  0,
1113  currentChunkStart,
1114  fullX,
1115  fullY,
1116  currentChunkSize,
1117  GL_RGBA,
1118  GL_UNSIGNED_BYTE,
1119  tmpPtr);
1120  currentChunk++;
1121  }
1122  }
1123 
1124  delete[] tmpPtr;
1125 }
1126 
1127 class RGBACompute
1128 {
1129  unsigned char *dataPtr;
1130  unsigned char *tmpPtr;
1131  unsigned char *tmpPtr2;
1132  int sizeX;
1133  int sizeY;
1134  int sizeZ;
1135  int sizeXY;
1136  int sizeXm1;
1137  int sizeYm1;
1138  int sizeZm1;
1139  int fullX;
1140  int fullY;
1141  int fullZ;
1142  int fullXY;
1143  // int currentChunkStart;
1144  // int currentChunkEnd;
1145 
1146  int offZ;
1147 
1148 public:
1149  RGBACompute(unsigned char *_dataPtr,
1150  unsigned char *_tmpPtr,
1151  unsigned char *_tmpPtr2,
1152  int _sizeX,
1153  int _sizeY,
1154  int _sizeZ,
1155  int _fullX,
1156  int _fullY,
1157  int _fullZ)
1158  {
1159  dataPtr = _dataPtr;
1160  tmpPtr = _tmpPtr;
1161  tmpPtr2 = _tmpPtr2;
1162  sizeX = _sizeX;
1163  sizeY = _sizeY;
1164  sizeZ = _sizeZ;
1165  fullX = _fullX;
1166  fullY = _fullY;
1167  fullZ = _fullZ;
1168 
1169  sizeXY = sizeX * sizeY;
1170  sizeXm1 = sizeX - 1;
1171  sizeYm1 = sizeY - 1;
1172  sizeZm1 = sizeZ - 1;
1173 
1174  fullXY = fullX * fullY;
1175  }
1176 
1177  inline int sample(int x, int y, int z) { return dataPtr[(x + y * sizeX + z * sizeXY) * 4 + 3]; }
1178  inline void fill(int x, int y, int z)
1179  {
1180  int doff = x + y * fullX + (z - offZ) * fullXY;
1181 
1182  tmpPtr[doff * 4 + 0] = 0;
1183  tmpPtr[doff * 4 + 1] = 0;
1184  tmpPtr[doff * 4 + 2] = 0;
1185  tmpPtr[doff * 4 + 3] = 0;
1186 
1187  tmpPtr2[doff * 3 + 0] = 0;
1188  tmpPtr2[doff * 3 + 1] = 0;
1189  tmpPtr2[doff * 3 + 2] = 0;
1190  }
1191 
1192  inline int clamp(int x)
1193  {
1194  if (x < 0)
1195  x = 0;
1196  else if (x > 255)
1197  x = 255;
1198  return x;
1199  }
1200 
1201  inline void write(int x, int y, int z, int iGrayValue, int gx, int gy, int gz)
1202  {
1203  /*
1204  gx /= aspect[0];
1205  gy /= aspect[1];
1206  gz /= aspect[2];
1207  */
1208  int nx = static_cast<int>(0.5f * gx + 127.5f);
1209  int ny = static_cast<int>(0.5f * gy + 127.5f);
1210  int nz = static_cast<int>(0.5f * gz + 127.5f);
1211 
1212  int doff = x + y * fullX + (z - offZ) * fullXY;
1213 
1214  // tmpPtr[doff*2+0]= 0;
1215 
1216  tmpPtr[doff * 4 + 0] = clamp(nx);
1217  tmpPtr[doff * 4 + 1] = clamp(ny);
1218  tmpPtr[doff * 4 + 2] = clamp(nz);
1219  tmpPtr[doff * 4 + 3] = clamp(iGrayValue);
1220 
1221  int soff = x + y * sizeX + z * sizeXY;
1222 
1223  tmpPtr2[doff * 3 + 0] = dataPtr[soff * 4 + 0];
1224  tmpPtr2[doff * 3 + 1] = dataPtr[soff * 4 + 1];
1225  tmpPtr2[doff * 3 + 2] = dataPtr[soff * 4 + 2];
1226 
1227  /*
1228  if( z == fullZ/2 )
1229  if( y == fullY/2 )
1230  MITK_INFO << x << " " << y << " " << z << " : " << iGrayValue << " : " << iGradient;
1231  */
1232  }
1233 
1234  inline void compute(int x, int y, int z)
1235  {
1236  int grayValue = sample(x, y, z);
1237  int gx, gy, gz;
1238 
1239  gx = sample(x + 1, y, z) - sample(x - 1, y, z);
1240  gy = sample(x, y + 1, z) - sample(x, y - 1, z);
1241  gz = sample(x, y, z + 1) - sample(x, y, z - 1);
1242 
1243  write(x, y, z, grayValue, gx, gy, gz);
1244  }
1245 
1246  inline void computeClamp(int x, int y, int z)
1247  {
1248  int grayValue = sample(x, y, z);
1249  int gx, gy, gz;
1250 
1251  if (x == 0)
1252  gx = 2 * (sample(x + 1, y, z) - grayValue);
1253  else if (x == sizeXm1)
1254  gx = 2 * (grayValue - sample(x - 1, y, z));
1255  else
1256  gx = sample(x + 1, y, z) - sample(x - 1, y, z);
1257 
1258  if (y == 0)
1259  gy = 2 * (sample(x, y + 1, z) - grayValue);
1260  else if (y == sizeYm1)
1261  gy = 2 * (grayValue - sample(x, y - 1, z));
1262  else
1263  gy = sample(x, y + 1, z) - sample(x, y - 1, z);
1264 
1265  if (z == 0)
1266  gz = 2 * (sample(x, y, z + 1) - grayValue);
1267  else if (z == sizeZm1)
1268  gz = 2 * (grayValue - sample(x, y, z - 1));
1269  else
1270  gz = sample(x, y, z + 1) - sample(x, y, z - 1);
1271 
1272  write(x, y, z, grayValue, gx, gy, gz);
1273  }
1274 
1275  inline void compute1D(int y, int z)
1276  {
1277  int x = 0;
1278 
1279  computeClamp(x, y, z);
1280  x++;
1281 
1282  while (x < sizeX - 1)
1283  {
1284  compute(x, y, z);
1285  x++;
1286  }
1287 
1288  if (x < sizeX)
1289  {
1290  computeClamp(x, y, z);
1291  x++;
1292  }
1293 
1294  while (x < fullX)
1295  {
1296  fill(x, y, z);
1297  x++;
1298  }
1299  }
1300 
1301  inline void fill1D(int y, int z)
1302  {
1303  int x = 0;
1304 
1305  while (x < fullX)
1306  {
1307  fill(x, y, z);
1308  x++;
1309  }
1310  }
1311 
1312  inline void computeClamp1D(int y, int z)
1313  {
1314  int x = 0;
1315 
1316  while (x < sizeX)
1317  {
1318  computeClamp(x, y, z);
1319  x++;
1320  }
1321 
1322  while (x < fullX)
1323  {
1324  fill(x, y, z);
1325  x++;
1326  }
1327  }
1328 
1329  inline void computeClamp2D(int z)
1330  {
1331  int y = 0;
1332 
1333  while (y < sizeY)
1334  {
1335  computeClamp1D(y, z);
1336  y++;
1337  }
1338 
1339  while (y < fullY)
1340  {
1341  fill1D(y, z);
1342  y++;
1343  }
1344  }
1345 
1346  inline void compute2D(int z)
1347  {
1348  int y = 0;
1349 
1350  computeClamp1D(y, z);
1351  y++;
1352 
1353  while (y < sizeY - 1)
1354  {
1355  compute1D(y, z);
1356  y++;
1357  }
1358 
1359  if (y < sizeY)
1360  {
1361  computeClamp1D(y, z);
1362  y++;
1363  }
1364 
1365  while (y < fullY)
1366  {
1367  fill1D(y, z);
1368  y++;
1369  }
1370  }
1371 
1372  inline void fill2D(int z)
1373  {
1374  int y = 0;
1375 
1376  while (y < fullY)
1377  {
1378  fill1D(y, z);
1379  y++;
1380  }
1381  }
1382 
1383  inline void fillSlices(int currentChunkStart, int currentChunkEnd)
1384  {
1385  offZ = currentChunkStart;
1386 
1387 #pragma omp parallel for
1388  for (int z = currentChunkStart; z <= currentChunkEnd; z++)
1389  {
1390  if (z == 0 || z == sizeZ - 1)
1391  computeClamp2D(z);
1392  else if (z >= sizeZ)
1393  fill2D(z);
1394  else
1395  compute2D(z);
1396  }
1397  }
1398 };
1399 
1400 void vtkVolumeTextureMapper3DComputeRGBA(unsigned char *dataPtr,
1402  GLuint volume1,
1403  GLuint volume2)
1404 {
1405  // unsigned char *inPtr;
1406  // unsigned char *outPtr, *outPtr2;
1407  // int i, j, k;
1408  // int idx;
1409 
1410  int inputDimensions[3];
1411  double inputSpacing[3];
1412  vtkImageData *input = me->GetInput();
1413 
1414  input->GetDimensions(inputDimensions);
1415  input->GetSpacing(inputSpacing);
1416 
1417  int outputDimensions[3];
1418  float outputSpacing[3];
1419  me->GetVolumeDimensions(outputDimensions);
1420  me->GetVolumeSpacing(outputSpacing);
1421 
1422  int components = input->GetNumberOfScalarComponents();
1423 
1424  MITK_INFO << "components are " << components;
1425 
1426  // double wx, wy, wz;
1427  // double fx, fy, fz;
1428  // int x, y, z;
1429 
1430  // double sampleRate[3];
1431  // sampleRate[0] = outputSpacing[0] / static_cast<double>(inputSpacing[0]);
1432  // sampleRate[1] = outputSpacing[1] / static_cast<double>(inputSpacing[1]);
1433  // sampleRate[2] = outputSpacing[2] / static_cast<double>(inputSpacing[2]);
1434 
1435  int fullX = outputDimensions[0];
1436  int fullY = outputDimensions[1];
1437  int fullZ = outputDimensions[2];
1438 
1439  int sizeX = inputDimensions[0];
1440  int sizeY = inputDimensions[1];
1441  int sizeZ = inputDimensions[2];
1442 
1443  int chunkSize = 64;
1444 
1445  if (fullZ < chunkSize)
1446  chunkSize = fullZ;
1447 
1448  int numChunks = (fullZ + (chunkSize - 1)) / chunkSize;
1449 
1450  // inPtr = dataPtr;
1451 
1452  unsigned char *tmpPtr = new unsigned char[fullX * fullY * chunkSize * 4];
1453  unsigned char *tmpPtr2 = new unsigned char[fullX * fullY * chunkSize * 3];
1454 
1455  // For each Chunk
1456  {
1457  RGBACompute sgc(dataPtr, tmpPtr, tmpPtr2, sizeX, sizeY, sizeZ, fullX, fullY, fullZ);
1458 
1459  int currentChunk = 0;
1460 
1461  while (currentChunk < numChunks)
1462  {
1463  // MITK_INFO << "processing chunk " << currentChunk;
1464 
1465  int currentChunkStart = currentChunk * chunkSize;
1466  int currentChunkEnd = currentChunkStart + chunkSize - 1;
1467 
1468  if (currentChunkEnd > (fullZ - 1))
1469  currentChunkEnd = (fullZ - 1);
1470 
1471  int currentChunkSize = currentChunkEnd - currentChunkStart + 1;
1472 
1473  sgc.fillSlices(currentChunkStart, currentChunkEnd);
1474 
1475  glBindTexture(vtkgl::TEXTURE_3D, volume1);
1476  vtkgl::TexSubImage3D(vtkgl::TEXTURE_3D,
1477  0,
1478  0,
1479  0,
1480  currentChunkStart,
1481  fullX,
1482  fullY,
1483  currentChunkSize,
1484  GL_RGBA,
1485  GL_UNSIGNED_BYTE,
1486  tmpPtr);
1487 
1488  glBindTexture(vtkgl::TEXTURE_3D, volume2);
1489  vtkgl::TexSubImage3D(vtkgl::TEXTURE_3D,
1490  0,
1491  0,
1492  0,
1493  currentChunkStart,
1494  fullX,
1495  fullY,
1496  currentChunkSize,
1497  GL_RGB,
1498  GL_UNSIGNED_BYTE,
1499  tmpPtr2);
1500 
1501  currentChunk++;
1502  }
1503  }
1504 
1505  delete[] tmpPtr;
1506  delete[] tmpPtr2;
1507 }
1508 
1509 //-----------------------------------------------------------------------------
1511 {
1512  // Get the image data
1513  vtkImageData *input = this->GetInput();
1514 
1515  // How big does the Volume need to be?
1516  int dim[3];
1517  input->GetDimensions(dim);
1518 
1519  int powerOfTwoDim[3];
1520 
1522  {
1523  for (int i = 0; i < 3; i++)
1524  powerOfTwoDim[i] = (dim[i] + 1) & ~1;
1525 
1526  // MITK_INFO << "using non-power-two even textures (" <<
1527  // (1.0-double(dim[0]*dim[1]*dim[2])/double(powerOfTwoDim[0]*powerOfTwoDim[1]*powerOfTwoDim[2])) * 100.0 << "%
1528  // memory wasted)";
1529  }
1530  else
1531  {
1532  for (int i = 0; i < 3; i++)
1533  {
1534  powerOfTwoDim[i] = 4;
1535  while (powerOfTwoDim[i] < dim[i])
1536  powerOfTwoDim[i] *= 2;
1537  }
1538 
1539  MITK_WARN << "using power-two textures ("
1540  << (1.0 -
1541  double(dim[0] * dim[1] * dim[2]) / double(powerOfTwoDim[0] * powerOfTwoDim[1] * powerOfTwoDim[2])) *
1542  100.0
1543  << "% memory wasted)";
1544  }
1545 
1546  // Save the volume size
1547  this->VolumeDimensions[0] = powerOfTwoDim[0];
1548  this->VolumeDimensions[1] = powerOfTwoDim[1];
1549  this->VolumeDimensions[2] = powerOfTwoDim[2];
1550 
1551  // What is the spacing?
1552  double spacing[3];
1553  input->GetSpacing(spacing);
1554 
1555  // Compute the new spacing
1556  this->VolumeSpacing[0] = (dim[0] - 1.01) * spacing[0] / static_cast<double>(this->VolumeDimensions[0] - 1);
1557  this->VolumeSpacing[1] = (dim[1] - 1.01) * spacing[1] / static_cast<double>(this->VolumeDimensions[1] - 1);
1558  this->VolumeSpacing[2] = ((dim[2]) - 1.01) * spacing[2] / static_cast<double>(this->VolumeDimensions[2] - 1);
1559 }
1560 
1561 //-----------------------------------------------------------------------------
1563 {
1564  // Get the image data
1565  vtkImageData *input = this->GetInput();
1566  // input->Update(); //VTK6_TODO
1567 
1568  bool needUpdate = false;
1569 
1570  // Has the volume changed in some way?
1571  if (this->SavedTextureInput != input || this->SavedTextureMTime.GetMTime() < input->GetMTime())
1572  needUpdate = true;
1573 
1574  // Do we have any volume on the gpu already?
1575  if (!this->Volume1Index)
1576  needUpdate = true;
1577 
1578  if (!needUpdate)
1579  return true;
1580 
1582 
1583  int components = input->GetNumberOfScalarComponents();
1584 
1585  // Find the scalar range
1586  double scalarRange[2];
1587  input->GetPointData()->GetScalars()->GetRange(scalarRange, components - 1);
1588 
1589  // Is the difference between max and min less than 4096? If so, and if
1590  // the data is not of float or double type, use a simple offset mapping.
1591  // If the difference between max and min is 4096 or greater, or the data
1592  // is of type float or double, we must use an offset / scaling mapping.
1593  // In this case, the array size will be 4096 - we need to figure out the
1594  // offset and scale factor.
1595  float offset;
1596  float scale;
1597 
1598  int arraySizeNeeded;
1599 
1600  int scalarType = input->GetScalarType();
1601 
1602  if (scalarType == VTK_FLOAT || scalarType == VTK_DOUBLE || scalarRange[1] - scalarRange[0] > 255)
1603  {
1604  arraySizeNeeded = 256;
1605  offset = -scalarRange[0];
1606  scale = 255.0 / (scalarRange[1] - scalarRange[0]);
1607  }
1608  else
1609  {
1610  arraySizeNeeded = static_cast<int>(scalarRange[1] - scalarRange[0] + 1);
1611  offset = -scalarRange[0];
1612  scale = 1.0;
1613  }
1614 
1615  this->ColorTableSize = arraySizeNeeded;
1616  this->ColorTableOffset = offset;
1617  this->ColorTableScale = scale;
1618 
1619  // Allocating volume on gpu
1620  {
1621  // Deleting old textures
1622  this->DeleteTextureIndex(&this->Volume1Index);
1623  this->DeleteTextureIndex(&this->Volume2Index);
1624  this->DeleteTextureIndex(&this->Volume3Index);
1625 
1626  this->CreateTextureIndex(&this->Volume1Index);
1627  // this->CreateTextureIndex(&this->Volume2Index);
1628 
1629  int dim[3];
1630  this->GetVolumeDimensions(dim);
1631 
1632  vtkgl::ActiveTexture(vtkgl::TEXTURE0);
1633 
1634  MITK_INFO << "allocating volume on gpu";
1635 
1636  GLint gradientScalarTextureFormat = GL_RGBA8;
1637 
1639  gradientScalarTextureFormat = myGL_COMPRESSED_RGBA_S3TC_DXT5_EXT;
1640 
1641  glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1642  vtkgl::TexImage3D(
1643  vtkgl::TEXTURE_3D, 0, gradientScalarTextureFormat, dim[0], dim[1], dim[2], 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);
1644  this->Setup3DTextureParameters(true);
1645  }
1646 
1647  // Transfer the input volume to the RGBA volume
1648  void *dataPtr = input->GetScalarPointer();
1649 
1650  switch (scalarType)
1651  {
1652  vtkTemplateMacro(vtkVolumeTextureMapper3DComputeScalars(
1653  static_cast<VTK_TT *>(dataPtr), this, offset, scale, this->Volume1Index, this->Volume2Index));
1654  }
1655 
1656  this->SavedTextureInput = input;
1657  this->SavedTextureMTime.Modified();
1658 
1659  return true;
1660 }
1661 
1662 //-----------------------------------------------------------------------------
1664 {
1665  // Get the image data
1666  vtkImageData *input = this->GetInput();
1667  // input->Update(); //VTK6_TODO
1668 
1669  bool needUpdate = false;
1670 
1671  // Has the volume changed in some way?
1672  if (this->SavedTextureInput != input || this->SavedTextureMTime.GetMTime() < input->GetMTime())
1673  needUpdate = true;
1674 
1675  // Do we have any volume on the gpu already?
1676  if (!this->Volume1Index)
1677  needUpdate = true;
1678 
1679  if (!needUpdate)
1680  return true;
1681 
1682  MITK_INFO << "updating rgba volume";
1683 
1685 
1686  // Allocating volume on gpu
1687  {
1688  // Deleting old textures
1689  this->DeleteTextureIndex(&this->Volume1Index);
1690  this->DeleteTextureIndex(&this->Volume2Index);
1691  this->DeleteTextureIndex(&this->Volume3Index);
1692 
1693  this->CreateTextureIndex(&this->Volume1Index);
1694  this->CreateTextureIndex(&this->Volume2Index);
1695 
1696  int dim[3];
1697  this->GetVolumeDimensions(dim);
1698 
1699  MITK_INFO << "allocating volume on gpu";
1700 
1701  GLint gradientScalarTextureFormat = GL_RGBA8;
1702  GLint colorTextureFormat = GL_RGB8;
1703 
1705  {
1706  gradientScalarTextureFormat = myGL_COMPRESSED_RGBA_S3TC_DXT5_EXT;
1707  colorTextureFormat = myGL_COMPRESSED_RGB_S3TC_DXT1_EXT;
1708  }
1709 
1710  vtkgl::ActiveTexture(vtkgl::TEXTURE0);
1711  glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1712  vtkgl::TexImage3D(
1713  vtkgl::TEXTURE_3D, 0, gradientScalarTextureFormat, dim[0], dim[1], dim[2], 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);
1714  this->Setup3DTextureParameters(true);
1715 
1716  glBindTexture(vtkgl::TEXTURE_3D, this->Volume2Index);
1717  vtkgl::TexImage3D(vtkgl::TEXTURE_3D, 0, colorTextureFormat, dim[0], dim[1], dim[2], 0, GL_RGB, GL_UNSIGNED_BYTE, 0);
1718  this->Setup3DTextureParameters(true);
1719  }
1720 
1721  // Transfer the input volume to the RGBA volume
1722  unsigned char *dataPtr = (unsigned char *)input->GetScalarPointer();
1723  vtkVolumeTextureMapper3DComputeRGBA(dataPtr, this, this->Volume1Index, this->Volume2Index);
1724 
1725  this->SavedTextureInput = input;
1726  this->SavedTextureMTime.Modified();
1727 
1728  return true;
1729 }
1730 
1732 {
1733  // GPU_INFO << "Setup3DTextureParameters";
1734 
1735  if (linear)
1736  {
1737  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1738  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
1739  }
1740  else
1741  {
1742  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
1743  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
1744  }
1745 
1746  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
1747  glTexParameterf(vtkgl::TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
1748 }
1749 
1750 void vtkMitkOpenGLVolumeTextureMapper3D::SetupOneIndependentTextures(vtkRenderer *vtkNotUsed(ren), vtkVolume *vol)
1751 {
1752  // Update the volume containing the 2 byte scalar / gradient magnitude
1753  this->UpdateVolumes(vol);
1754 
1755  // Update the dependent 2D color table mapping scalar value and
1756  // gradient magnitude to RGBA
1757  if (this->UpdateColorLookup(vol) || !this->ColorLookupIndex)
1758  {
1759  this->DeleteTextureIndex(&this->ColorLookupIndex);
1760  this->DeleteTextureIndex(&this->AlphaLookupIndex);
1761 
1762  this->CreateTextureIndex(&this->ColorLookupIndex);
1763 
1764  vtkgl::ActiveTexture(vtkgl::TEXTURE1);
1765  glBindTexture(GL_TEXTURE_2D, this->ColorLookupIndex);
1766 
1767  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
1768  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
1769  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
1770  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
1771 
1772  // MITK_INFO << "uploading transferfunction";
1773 
1774  GLint colorLookupTextureFormat = GL_RGBA8;
1775 
1777  colorLookupTextureFormat = myGL_COMPRESSED_RGBA_S3TC_DXT5_EXT;
1778 
1779  glTexImage2D(GL_TEXTURE_2D, 0, colorLookupTextureFormat, 256, 256, 0, GL_RGBA, GL_UNSIGNED_BYTE, this->ColorLookup);
1780  }
1781 }
1782 
1783 void vtkMitkOpenGLVolumeTextureMapper3D::SetupRGBATextures(vtkRenderer *vtkNotUsed(ren), vtkVolume *vol)
1784 {
1785  MITK_INFO << "SetupFourDependentTextures";
1786 
1787  this->UpdateVolumesRGBA(vol);
1788 
1789  /*
1790 vtkgl::ActiveTexture( vtkgl::TEXTURE0 );
1791 glDisable( GL_TEXTURE_2D );
1792 glEnable( vtkgl::TEXTURE_3D );
1793 
1794 vtkgl::ActiveTexture( vtkgl::TEXTURE1 );
1795 glDisable( GL_TEXTURE_2D );
1796 glEnable( vtkgl::TEXTURE_3D );
1797 
1798 vtkgl::ActiveTexture( vtkgl::TEXTURE2 );
1799 glDisable( GL_TEXTURE_2D );
1800 glEnable( vtkgl::TEXTURE_3D );
1801 
1802 // Update the volume containing the 3 byte scalars / gradient magnitude
1803 if ( this->UpdateVolumes( vol ) || !this->Volume1Index ||
1804  !this->Volume2Index || !this->Volume3Index )
1805  {
1806  int dim[3];
1807  this->GetVolumeDimensions(dim);
1808 
1809  vtkgl::ActiveTexture( vtkgl::TEXTURE0 );
1810  glBindTexture(vtkgl::TEXTURE_3D,0);
1811  this->DeleteTextureIndex(&this->Volume1Index);
1812  this->CreateTextureIndex(&this->Volume1Index);
1813  glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1814  vtkgl::TexImage3D(vtkgl::TEXTURE_3D,0,this->InternalRGB,dim[0],dim[1],
1815  dim[2],0,GL_RGB,GL_UNSIGNED_BYTE,this->Volume1);
1816 
1817  vtkgl::ActiveTexture( vtkgl::TEXTURE1 );
1818  glBindTexture(vtkgl::TEXTURE_3D,0);
1819  this->DeleteTextureIndex(&this->Volume2Index);
1820  this->CreateTextureIndex(&this->Volume2Index);
1821  glBindTexture(vtkgl::TEXTURE_3D, this->Volume2Index);
1822  vtkgl::TexImage3D(vtkgl::TEXTURE_3D,0,this->InternalLA,dim[0],dim[1],
1823  dim[2],0,GL_LUMINANCE_ALPHA,GL_UNSIGNED_BYTE,
1824  this->Volume2);
1825 
1826  vtkgl::ActiveTexture( vtkgl::TEXTURE2 );
1827  glBindTexture(vtkgl::TEXTURE_3D,0);
1828  this->DeleteTextureIndex(&this->Volume3Index);
1829  this->CreateTextureIndex(&this->Volume3Index);
1830  glBindTexture(vtkgl::TEXTURE_3D, this->Volume3Index);
1831  vtkgl::TexImage3D(vtkgl::TEXTURE_3D,0,this->InternalRGB,dim[0],dim[1],
1832  dim[2],0,GL_RGB,GL_UNSIGNED_BYTE,this->Volume3);
1833  }
1834 
1835 vtkgl::ActiveTexture( vtkgl::TEXTURE0 );
1836 glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1837 this->Setup3DTextureParameters( true );
1838 
1839 vtkgl::ActiveTexture( vtkgl::TEXTURE1 );
1840 glBindTexture(vtkgl::TEXTURE_3D, this->Volume2Index);
1841 this->Setup3DTextureParameters( true );
1842 
1843 vtkgl::ActiveTexture( vtkgl::TEXTURE2 );
1844 glBindTexture(vtkgl::TEXTURE_3D_EXT, this->Volume3Index);
1845 this->Setup3DTextureParameters( true );
1846 
1847 vtkgl::ActiveTexture( vtkgl::TEXTURE3 );
1848 glEnable( GL_TEXTURE_2D );
1849 glDisable( vtkgl::TEXTURE_3D );
1850 
1851 // Update the dependent 2D table mapping scalar value and
1852 // gradient magnitude to opacity
1853 if ( this->UpdateColorLookup( vol ) || !this->AlphaLookupIndex )
1854  {
1855  this->DeleteTextureIndex(&this->ColorLookupIndex);
1856 
1857  vtkgl::ActiveTexture( vtkgl::TEXTURE3 );
1858  glBindTexture(GL_TEXTURE_2D,0);
1859  this->DeleteTextureIndex(&this->AlphaLookupIndex);
1860  this->CreateTextureIndex(&this->AlphaLookupIndex);
1861  glBindTexture(GL_TEXTURE_2D, this->AlphaLookupIndex);
1862 
1863  glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
1864  glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST );
1865  glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP );
1866  glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP );
1867 
1868  //MITK_INFO << "uploading transferfunction";
1869 
1870  glTexImage2D(GL_TEXTURE_2D,0,this->InternalAlpha, 256, 256, 0,
1871  GL_ALPHA, GL_UNSIGNED_BYTE, this->AlphaLookup );
1872  }
1873 
1874 vtkgl::ActiveTexture( vtkgl::TEXTURE3 );
1875 glBindTexture(GL_TEXTURE_2D, this->AlphaLookupIndex);
1876 
1877 */
1878 }
1879 
1881 {
1882  // GPU_INFO << "RenderOneIndependentShadeFP";
1883 
1884  this->SetupOneIndependentTextures(ren, vol);
1885 
1886  glEnable(vtkgl::FRAGMENT_PROGRAM_ARB);
1887 
1888  vtkgl::BindProgramARB(vtkgl::FRAGMENT_PROGRAM_ARB, prgOneComponentShade);
1889 
1890  this->SetupProgramLocalsForShadingFP(ren, vol);
1891 
1892  // Bind Textures
1893  {
1894  vtkgl::ActiveTexture(vtkgl::TEXTURE0);
1895  glDisable(GL_TEXTURE_2D);
1896  glEnable(vtkgl::TEXTURE_3D);
1897  glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1898 
1899  vtkgl::ActiveTexture(vtkgl::TEXTURE1);
1900  glEnable(GL_TEXTURE_2D);
1901  glDisable(vtkgl::TEXTURE_3D);
1902  glBindTexture(GL_TEXTURE_2D, this->ColorLookupIndex);
1903 
1904  vtkgl::ActiveTexture(vtkgl::TEXTURE2);
1905  glDisable(GL_TEXTURE_2D);
1906  glEnable(vtkgl::TEXTURE_3D);
1907  glBindTexture(vtkgl::TEXTURE_3D, this->Volume2Index);
1908  }
1909 
1910  int stages[4] = {1, 1, 1, 0};
1911  this->RenderPolygons(ren, vol, stages);
1912 
1913  glDisable(vtkgl::FRAGMENT_PROGRAM_ARB);
1914 }
1915 
1916 void vtkMitkOpenGLVolumeTextureMapper3D::RenderRGBAShadeFP(vtkRenderer *ren, vtkVolume *vol)
1917 {
1918  this->SetupRGBATextures(ren, vol);
1919 
1920  glEnable(vtkgl::FRAGMENT_PROGRAM_ARB);
1921 
1922  vtkgl::BindProgramARB(vtkgl::FRAGMENT_PROGRAM_ARB, prgRGBAShade);
1923 
1924  this->SetupProgramLocalsForShadingFP(ren, vol);
1925 
1926  // Bind Textures
1927  {
1928  vtkgl::ActiveTexture(vtkgl::TEXTURE0);
1929  glDisable(GL_TEXTURE_2D);
1930  glEnable(vtkgl::TEXTURE_3D);
1931  glBindTexture(vtkgl::TEXTURE_3D, this->Volume1Index);
1932 
1933  vtkgl::ActiveTexture(vtkgl::TEXTURE1);
1934  glDisable(GL_TEXTURE_2D);
1935  glEnable(vtkgl::TEXTURE_3D);
1936  glBindTexture(vtkgl::TEXTURE_3D, this->Volume2Index);
1937  }
1938 
1939  int stages[4] = {1, 1, 1, 0};
1940  this->RenderPolygons(ren, vol, stages);
1941 
1942  glDisable(vtkgl::FRAGMENT_PROGRAM_ARB);
1943 }
1944 
1946  vtkVolume *vol,
1947  GLfloat lightDirection[2][4],
1948  GLfloat lightDiffuseColor[2][4],
1949  GLfloat lightSpecularColor[2][4],
1950  GLfloat halfwayVector[2][4],
1951  GLfloat ambientColor[4])
1952 {
1953  // GPU_INFO << "GetLightInformation";
1954 
1955  float ambient = vol->GetProperty()->GetAmbient();
1956  float diffuse = vol->GetProperty()->GetDiffuse();
1957  float specular = vol->GetProperty()->GetSpecular();
1958 
1959  vtkTransform *volumeTransform = vtkTransform::New();
1960 
1961  volumeTransform->SetMatrix(vol->GetMatrix());
1962  volumeTransform->Inverse();
1963 
1964  vtkLightCollection *lights = ren->GetLights();
1965  lights->InitTraversal();
1966 
1967  vtkLight *light[2];
1968  light[0] = lights->GetNextItem();
1969  light[1] = lights->GetNextItem();
1970 
1971  int lightIndex = 0;
1972 
1973  double cameraPosition[3];
1974  double cameraFocalPoint[3];
1975 
1976  ren->GetActiveCamera()->GetPosition(cameraPosition);
1977  ren->GetActiveCamera()->GetFocalPoint(cameraFocalPoint);
1978 
1979  double viewDirection[3];
1980 
1981  volumeTransform->TransformPoint(cameraPosition, cameraPosition);
1982  volumeTransform->TransformPoint(cameraFocalPoint, cameraFocalPoint);
1983 
1984  viewDirection[0] = cameraFocalPoint[0] - cameraPosition[0];
1985  viewDirection[1] = cameraFocalPoint[1] - cameraPosition[1];
1986  viewDirection[2] = cameraFocalPoint[2] - cameraPosition[2];
1987 
1988  vtkMath::Normalize(viewDirection);
1989 
1990  ambientColor[0] = 0.0;
1991  ambientColor[1] = 0.0;
1992  ambientColor[2] = 0.0;
1993  ambientColor[3] = 0.0;
1994 
1995  for (lightIndex = 0; lightIndex < 2; lightIndex++)
1996  {
1997  float dir[3] = {0, 0, 0};
1998  float half[3] = {0, 0, 0};
1999 
2000  if (light[lightIndex] == NULL || light[lightIndex]->GetSwitch() == 0)
2001  {
2002  lightDiffuseColor[lightIndex][0] = 0.0;
2003  lightDiffuseColor[lightIndex][1] = 0.0;
2004  lightDiffuseColor[lightIndex][2] = 0.0;
2005  lightDiffuseColor[lightIndex][3] = 0.0;
2006 
2007  lightSpecularColor[lightIndex][0] = 0.0;
2008  lightSpecularColor[lightIndex][1] = 0.0;
2009  lightSpecularColor[lightIndex][2] = 0.0;
2010  lightSpecularColor[lightIndex][3] = 0.0;
2011  }
2012  else
2013  {
2014  float lightIntensity = light[lightIndex]->GetIntensity();
2015  double lightColor[3];
2016 
2017  light[lightIndex]->GetDiffuseColor(lightColor);
2018 
2019  double lightPosition[3];
2020  double lightFocalPoint[3];
2021  light[lightIndex]->GetTransformedPosition(lightPosition);
2022  light[lightIndex]->GetTransformedFocalPoint(lightFocalPoint);
2023 
2024  volumeTransform->TransformPoint(lightPosition, lightPosition);
2025  volumeTransform->TransformPoint(lightFocalPoint, lightFocalPoint);
2026 
2027  dir[0] = lightPosition[0] - lightFocalPoint[0];
2028  dir[1] = lightPosition[1] - lightFocalPoint[1];
2029  dir[2] = lightPosition[2] - lightFocalPoint[2];
2030 
2031  vtkMath::Normalize(dir);
2032 
2033  lightDiffuseColor[lightIndex][0] = lightColor[0] * diffuse * lightIntensity;
2034  lightDiffuseColor[lightIndex][1] = lightColor[1] * diffuse * lightIntensity;
2035  lightDiffuseColor[lightIndex][2] = lightColor[2] * diffuse * lightIntensity;
2036  lightDiffuseColor[lightIndex][3] = 1.0;
2037 
2038  lightSpecularColor[lightIndex][0] = lightColor[0] * specular * lightIntensity;
2039  lightSpecularColor[lightIndex][1] = lightColor[1] * specular * lightIntensity;
2040  lightSpecularColor[lightIndex][2] = lightColor[2] * specular * lightIntensity;
2041  lightSpecularColor[lightIndex][3] = 0.0;
2042 
2043  half[0] = dir[0] - viewDirection[0];
2044  half[1] = dir[1] - viewDirection[1];
2045  half[2] = dir[2] - viewDirection[2];
2046 
2047  vtkMath::Normalize(half);
2048 
2049  ambientColor[0] += ambient * lightColor[0];
2050  ambientColor[1] += ambient * lightColor[1];
2051  ambientColor[2] += ambient * lightColor[2];
2052  }
2053 
2054  lightDirection[lightIndex][0] = (dir[0] + 1.0) / 2.0;
2055  lightDirection[lightIndex][1] = (dir[1] + 1.0) / 2.0;
2056  lightDirection[lightIndex][2] = (dir[2] + 1.0) / 2.0;
2057  lightDirection[lightIndex][3] = 0.0;
2058 
2059  halfwayVector[lightIndex][0] = (half[0] + 1.0) / 2.0;
2060  halfwayVector[lightIndex][1] = (half[1] + 1.0) / 2.0;
2061  halfwayVector[lightIndex][2] = (half[2] + 1.0) / 2.0;
2062  halfwayVector[lightIndex][3] = 0.0;
2063  }
2064 
2065  volumeTransform->Delete();
2066 }
2067 
2069 {
2070  // GPU_INFO << "SetupProgramLocalsForShadingFP";
2071 
2072  GLfloat lightDirection[2][4];
2073  GLfloat lightDiffuseColor[2][4];
2074  GLfloat lightSpecularColor[2][4];
2075  GLfloat halfwayVector[2][4];
2076  GLfloat ambientColor[4];
2077 
2078  float ambient = vol->GetProperty()->GetAmbient();
2079  float diffuse = vol->GetProperty()->GetDiffuse();
2080  float specular = vol->GetProperty()->GetSpecular();
2081  float specularPower = vol->GetProperty()->GetSpecularPower();
2082 
2083  vtkTransform *volumeTransform = vtkTransform::New();
2084 
2085  volumeTransform->SetMatrix(vol->GetMatrix());
2086  volumeTransform->Inverse();
2087 
2088  vtkLightCollection *lights = ren->GetLights();
2089  lights->InitTraversal();
2090 
2091  vtkLight *light[2];
2092  light[0] = lights->GetNextItem();
2093  light[1] = lights->GetNextItem();
2094 
2095  int lightIndex = 0;
2096 
2097  double cameraPosition[3];
2098  double cameraFocalPoint[3];
2099 
2100  ren->GetActiveCamera()->GetPosition(cameraPosition);
2101  ren->GetActiveCamera()->GetFocalPoint(cameraFocalPoint);
2102 
2103  volumeTransform->TransformPoint(cameraPosition, cameraPosition);
2104  volumeTransform->TransformPoint(cameraFocalPoint, cameraFocalPoint);
2105 
2106  double viewDirection[4];
2107 
2108  viewDirection[0] = cameraFocalPoint[0] - cameraPosition[0];
2109  viewDirection[1] = cameraFocalPoint[1] - cameraPosition[1];
2110  viewDirection[2] = cameraFocalPoint[2] - cameraPosition[2];
2111  viewDirection[3] = 0.0;
2112 
2113  vtkMath::Normalize(viewDirection);
2114 
2115  ambientColor[0] = 0.0;
2116  ambientColor[1] = 0.0;
2117  ambientColor[2] = 0.0;
2118  ambientColor[3] = 0.0;
2119 
2120  for (lightIndex = 0; lightIndex < 2; lightIndex++)
2121  {
2122  float dir[3] = {0, 0, 0};
2123  float half[3] = {0, 0, 0};
2124 
2125  if (light[lightIndex] == NULL || light[lightIndex]->GetSwitch() == 0)
2126  {
2127  lightDiffuseColor[lightIndex][0] = 0.0;
2128  lightDiffuseColor[lightIndex][1] = 0.0;
2129  lightDiffuseColor[lightIndex][2] = 0.0;
2130  lightDiffuseColor[lightIndex][3] = 0.0;
2131 
2132  lightSpecularColor[lightIndex][0] = 0.0;
2133  lightSpecularColor[lightIndex][1] = 0.0;
2134  lightSpecularColor[lightIndex][2] = 0.0;
2135  lightSpecularColor[lightIndex][3] = 0.0;
2136  }
2137  else
2138  {
2139  float lightIntensity = light[lightIndex]->GetIntensity();
2140  double lightColor[3];
2141 
2142  light[lightIndex]->GetDiffuseColor(lightColor);
2143 
2144  double lightPosition[3];
2145  double lightFocalPoint[3];
2146  light[lightIndex]->GetTransformedPosition(lightPosition);
2147  light[lightIndex]->GetTransformedFocalPoint(lightFocalPoint);
2148 
2149  volumeTransform->TransformPoint(lightPosition, lightPosition);
2150  volumeTransform->TransformPoint(lightFocalPoint, lightFocalPoint);
2151 
2152  dir[0] = lightPosition[0] - lightFocalPoint[0];
2153  dir[1] = lightPosition[1] - lightFocalPoint[1];
2154  dir[2] = lightPosition[2] - lightFocalPoint[2];
2155 
2156  vtkMath::Normalize(dir);
2157 
2158  lightDiffuseColor[lightIndex][0] = lightColor[0] * diffuse * lightIntensity;
2159  lightDiffuseColor[lightIndex][1] = lightColor[1] * diffuse * lightIntensity;
2160  lightDiffuseColor[lightIndex][2] = lightColor[2] * diffuse * lightIntensity;
2161  lightDiffuseColor[lightIndex][3] = 0.0;
2162 
2163  lightSpecularColor[lightIndex][0] = lightColor[0] * specular * lightIntensity;
2164  lightSpecularColor[lightIndex][1] = lightColor[1] * specular * lightIntensity;
2165  lightSpecularColor[lightIndex][2] = lightColor[2] * specular * lightIntensity;
2166  lightSpecularColor[lightIndex][3] = 0.0;
2167 
2168  half[0] = dir[0] - viewDirection[0];
2169  half[1] = dir[1] - viewDirection[1];
2170  half[2] = dir[2] - viewDirection[2];
2171 
2172  vtkMath::Normalize(half);
2173 
2174  ambientColor[0] += ambient * lightColor[0];
2175  ambientColor[1] += ambient * lightColor[1];
2176  ambientColor[2] += ambient * lightColor[2];
2177  }
2178 
2179  lightDirection[lightIndex][0] = dir[0];
2180  lightDirection[lightIndex][1] = dir[1];
2181  lightDirection[lightIndex][2] = dir[2];
2182  lightDirection[lightIndex][3] = 0.0;
2183 
2184  halfwayVector[lightIndex][0] = half[0];
2185  halfwayVector[lightIndex][1] = half[1];
2186  halfwayVector[lightIndex][2] = half[2];
2187  halfwayVector[lightIndex][3] = 0.0;
2188  }
2189 
2190  volumeTransform->Delete();
2191 
2192  vtkgl::ProgramLocalParameter4fARB(vtkgl::FRAGMENT_PROGRAM_ARB,
2193  0,
2194  lightDirection[0][0],
2195  lightDirection[0][1],
2196  lightDirection[0][2],
2197  lightDirection[0][3]);
2198 
2199  vtkgl::ProgramLocalParameter4fARB(
2200  vtkgl::FRAGMENT_PROGRAM_ARB, 1, halfwayVector[0][0], halfwayVector[0][1], halfwayVector[0][2], halfwayVector[0][3]);
2201 
2202  vtkgl::ProgramLocalParameter4fARB(vtkgl::FRAGMENT_PROGRAM_ARB, 2, ambient, diffuse, specular, specularPower);
2203 
2204  vtkgl::ProgramLocalParameter4fARB(vtkgl::FRAGMENT_PROGRAM_ARB,
2205  3,
2206  lightDiffuseColor[0][0],
2207  lightDiffuseColor[0][1],
2208  lightDiffuseColor[0][2],
2209  lightDiffuseColor[0][3]);
2210 
2211  vtkgl::ProgramLocalParameter4fARB(vtkgl::FRAGMENT_PROGRAM_ARB,
2212  4,
2213  lightSpecularColor[0][0],
2214  lightSpecularColor[0][1],
2215  lightSpecularColor[0][2],
2216  lightSpecularColor[0][3]);
2217 
2218  vtkgl::ProgramLocalParameter4fARB(
2219  vtkgl::FRAGMENT_PROGRAM_ARB, 5, viewDirection[0], viewDirection[1], viewDirection[2], viewDirection[3]);
2220 
2221  vtkgl::ProgramLocalParameter4fARB(vtkgl::FRAGMENT_PROGRAM_ARB, 6, 2.0, -1.0, 0.0, 0.0);
2222 }
2223 
2224 int vtkMitkOpenGLVolumeTextureMapper3D::IsRenderSupported(vtkRenderer *renderer, vtkVolumeProperty * /*property*/)
2225 {
2226  // GPU_INFO << "IsRenderSupported";
2227 
2228  if (!this->Initialized)
2229  {
2230  // this->Initialize();
2231  this->Initialize(renderer);
2232  }
2233 
2234  if (!this->RenderPossible)
2235  {
2236  MITK_WARN << "vtkMitkOpenGLVolumeTextureMapper3D::IsRenderSupported Rendering not possible";
2237  return 0;
2238  }
2239 
2240  if (!this->GetInput())
2241  {
2242  MITK_WARN << "vtkMitkOpenGLVolumeTextureMapper3D::IsRenderSupported No input available";
2243  return 0;
2244  }
2245 
2246  return 1;
2247 }
2248 
2250 {
2251  // GPU_INFO << "Initialize";
2252 
2253  this->Initialized = 1;
2254  // vtkOpenGLExtensionManager * extensions = vtkOpenGLExtensionManager::New();
2255  // extensions->SetRenderWindow(NULL); // set render window to the current one.
2256  vtkOpenGLExtensionManager *extensions =
2257  static_cast<vtkOpenGLRenderWindow *>(renderer->GetRenderWindow())->GetExtensionManager();
2258 
2259  int supports_texture3D = extensions->ExtensionSupported("GL_VERSION_1_2");
2260  if (supports_texture3D)
2261  {
2262  extensions->LoadExtension("GL_VERSION_1_2");
2263  }
2264  else
2265  {
2266  supports_texture3D = extensions->ExtensionSupported("GL_EXT_texture3D");
2267  if (supports_texture3D)
2268  {
2269  extensions->LoadCorePromotedExtension("GL_EXT_texture3D");
2270  }
2271  }
2272 
2273  int supports_multitexture = extensions->ExtensionSupported("GL_VERSION_1_3");
2274  if (supports_multitexture)
2275  {
2276  extensions->LoadExtension("GL_VERSION_1_3");
2277  }
2278  else
2279  {
2280  supports_multitexture = extensions->ExtensionSupported("GL_ARB_multitexture");
2281  if (supports_multitexture)
2282  {
2283  extensions->LoadCorePromotedExtension("GL_ARB_multitexture");
2284  }
2285  }
2286 
2287  this->SupportsCompressedTexture = extensions->ExtensionSupported("GL_VERSION_1_3") == 1;
2288  if (!this->SupportsCompressedTexture)
2289  {
2290  this->SupportsCompressedTexture = extensions->ExtensionSupported("GL_ARB_texture_compression") == 1;
2291  if (this->SupportsCompressedTexture)
2292  {
2293  extensions->LoadCorePromotedExtension("GL_ARB_texture_compression");
2294  }
2295  }
2296  // GPU_INFO(this->SupportsCompressedTexture) << "supporting compressed textures";
2297 
2298  this->SupportsNonPowerOfTwoTextures = extensions->ExtensionSupported("GL_VERSION_2_0") ||
2299  extensions->ExtensionSupported("GL_ARB_texture_non_power_of_two");
2300 
2301  // GPU_INFO << "np2: " << (this->SupportsNonPowerOfTwoTextures?1:0);
2302 
2303  int supports_GL_ARB_fragment_program = extensions->ExtensionSupported("GL_ARB_fragment_program");
2304  if (supports_GL_ARB_fragment_program)
2305  {
2306  extensions->LoadExtension("GL_ARB_fragment_program");
2307  }
2308 
2309  int supports_GL_ARB_vertex_program = extensions->ExtensionSupported("GL_ARB_vertex_program");
2310  if (supports_GL_ARB_vertex_program)
2311  {
2312  extensions->LoadExtension("GL_ARB_vertex_program");
2313  }
2314 
2315  RenderPossible = 0;
2316 
2317  if (supports_texture3D && supports_multitexture && supports_GL_ARB_fragment_program &&
2318  supports_GL_ARB_vertex_program && vtkgl::TexImage3D && vtkgl::ActiveTexture && vtkgl::MultiTexCoord3fv &&
2319  vtkgl::GenProgramsARB && vtkgl::DeleteProgramsARB && vtkgl::BindProgramARB && vtkgl::ProgramStringARB &&
2320  vtkgl::ProgramLocalParameter4fARB)
2321  {
2322  RenderPossible = 1;
2323  }
2324  else
2325  {
2326  std::string errString =
2327  "no gpu-acceleration possible cause following extensions/methods are missing or unsupported:";
2328  if (!supports_texture3D)
2329  errString += " EXT_TEXTURE3D";
2330  if (!supports_multitexture)
2331  errString += " EXT_MULTITEXTURE";
2332  if (!supports_GL_ARB_fragment_program)
2333  errString += " ARB_FRAGMENT_PROGRAM";
2334  if (!supports_GL_ARB_vertex_program)
2335  errString += " ARB_VERTEX_PROGRAM";
2336  if (!vtkgl::TexImage3D)
2337  errString += " glTexImage3D";
2338  if (!vtkgl::ActiveTexture)
2339  errString += " glActiveTexture";
2340  if (!vtkgl::MultiTexCoord3fv)
2341  errString += " glMultiTexCoord3fv";
2342  if (!vtkgl::GenProgramsARB)
2343  errString += " glGenProgramsARB";
2344  if (!vtkgl::DeleteProgramsARB)
2345  errString += " glDeleteProgramsARB";
2346  if (!vtkgl::BindProgramARB)
2347  errString += " glBindProgramARB";
2348  if (!vtkgl::ProgramStringARB)
2349  errString += " glProgramStringARB";
2350  if (!vtkgl::ProgramLocalParameter4fARB)
2351  errString += " glProgramLocalParameter4fARB";
2352  GPU_WARN << errString;
2353  };
2354 
2355  if (RenderPossible)
2356  {
2357  vtkgl::GenProgramsARB(1, &prgOneComponentShade);
2358  vtkgl::BindProgramARB(vtkgl::FRAGMENT_PROGRAM_ARB, prgOneComponentShade);
2359  vtkgl::ProgramStringARB(vtkgl::FRAGMENT_PROGRAM_ARB,
2360  vtkgl::PROGRAM_FORMAT_ASCII_ARB,
2361  static_cast<GLsizei>(strlen(vtkMitkVolumeTextureMapper3D_OneComponentShadeFP)),
2363 
2364  vtkgl::GenProgramsARB(1, &prgRGBAShade);
2365  vtkgl::BindProgramARB(vtkgl::FRAGMENT_PROGRAM_ARB, prgRGBAShade);
2366  vtkgl::ProgramStringARB(vtkgl::FRAGMENT_PROGRAM_ARB,
2367  vtkgl::PROGRAM_FORMAT_ASCII_ARB,
2368  static_cast<GLsizei>(strlen(vtkMitkVolumeTextureMapper3D_FourDependentShadeFP)),
2370  }
2371 }
2372 
2373 // ----------------------------------------------------------------------------
2374 // Print the vtkMitkOpenGLVolumeTextureMapper3D
2375 void vtkMitkOpenGLVolumeTextureMapper3D::PrintSelf(ostream &os, vtkIndent indent)
2376 {
2377  // vtkOpenGLExtensionManager * extensions = vtkOpenGLExtensionManager::New();
2378  // extensions->SetRenderWindow(NULL); // set render window to current render window
2379 
2380  os << indent << "Initialized " << this->Initialized << endl;
2381  /* if ( this->Initialized )
2382  {
2383  os << indent << "Supports GL_VERSION_1_2:"
2384  << extensions->ExtensionSupported( "GL_VERSION_1_2" ) << endl;
2385  os << indent << "Supports GL_EXT_texture3D:"
2386  << extensions->ExtensionSupported( "GL_EXT_texture3D" ) << endl;
2387  os << indent << "Supports GL_VERSION_1_3:"
2388  << extensions->ExtensionSupported( "GL_VERSION_1_3" ) << endl;
2389  os << indent << "Supports GL_ARB_multitexture: "
2390  << extensions->ExtensionSupported( "GL_ARB_multitexture" ) << endl;
2391  os << indent << "Supports GL_NV_texture_shader2: "
2392  << extensions->ExtensionSupported( "GL_NV_texture_shader2" ) << endl;
2393  os << indent << "Supports GL_NV_register_combiners2: "
2394  << extensions->ExtensionSupported( "GL_NV_register_combiners2" )
2395  << endl;
2396  os << indent << "Supports GL_ATI_fragment_shader: "
2397  << extensions->ExtensionSupported( "GL_ATI_fragment_shader" ) << endl;
2398  os << indent << "Supports GL_ARB_fragment_program: "
2399  << extensions->ExtensionSupported( "GL_ARB_fragment_program" ) << endl;
2400  os << indent << "Supports GL_ARB_texture_compression: "
2401  << extensions->ExtensionSupported( "GL_ARB_texture_compression" )
2402  << endl;
2403  os << indent << "Supports GL_VERSION_2_0:"
2404  << extensions->ExtensionSupported( "GL_VERSION_2_0" )
2405  << endl;
2406  os << indent << "Supports GL_ARB_texture_non_power_of_two:"
2407  << extensions->ExtensionSupported( "GL_ARB_texture_non_power_of_two" )
2408  << endl;
2409  }
2410  extensions->Delete();
2411  */
2412 
2413  if (this->RenderWindow != 0)
2414  {
2415  vtkOpenGLExtensionManager *extensions =
2416  static_cast<vtkOpenGLRenderWindow *>(this->RenderWindow)->GetExtensionManager();
2417 
2418  if (this->Initialized)
2419  {
2420  os << indent << "Supports GL_VERSION_1_2:" << extensions->ExtensionSupported("GL_VERSION_1_2") << endl;
2421  os << indent << "Supports GL_EXT_texture3D:" << extensions->ExtensionSupported("GL_EXT_texture3D") << endl;
2422  os << indent << "Supports GL_VERSION_1_3:" << extensions->ExtensionSupported("GL_VERSION_1_3") << endl;
2423  os << indent << "Supports GL_ARB_multitexture: " << extensions->ExtensionSupported("GL_ARB_multitexture") << endl;
2424  os << indent << "Supports GL_NV_texture_shader2: " << extensions->ExtensionSupported("GL_NV_texture_shader2")
2425  << endl;
2426  os << indent
2427  << "Supports GL_NV_register_combiners2: " << extensions->ExtensionSupported("GL_NV_register_combiners2")
2428  << endl;
2429  os << indent << "Supports GL_ATI_fragment_shader: " << extensions->ExtensionSupported("GL_ATI_fragment_shader")
2430  << endl;
2431  os << indent << "Supports GL_ARB_fragment_program: " << extensions->ExtensionSupported("GL_ARB_fragment_program")
2432  << endl;
2433  os << indent
2434  << "Supports GL_ARB_texture_compression: " << extensions->ExtensionSupported("GL_ARB_texture_compression")
2435  << endl;
2436  os << indent << "Supports GL_VERSION_2_0:" << extensions->ExtensionSupported("GL_VERSION_2_0") << endl;
2437  os << indent << "Supports GL_ARB_texture_non_power_of_two:"
2438  << extensions->ExtensionSupported("GL_ARB_texture_non_power_of_two") << endl;
2439  }
2440  }
2441 
2442  this->Superclass::PrintSelf(os, indent);
2443 }
void RenderPolygons(vtkRenderer *ren, vtkVolume *vol, int stages[4])
const char * vtkMitkVolumeTextureMapper3D_FourDependentShadeFP
#define MITK_INFO
Definition: mitkLogMacros.h:22
const char * vtkMitkVolumeTextureMapper3D_OneComponentShadeFP
Organizes the rendering process.
void vtkVolumeTextureMapper3DComputeScalars(T *dataPtr, vtkMitkVolumeTextureMapper3D *me, float offset, float scale, GLuint volume1, GLuint)
void SetupProgramLocalsForShadingFP(vtkRenderer *ren, vtkVolume *vol)
#define myGL_COMPRESSED_RGB_S3TC_DXT1_EXT
void ComputePolygons(vtkRenderer *ren, vtkVolume *vol, double bounds[6])
#define myGL_COMPRESSED_RGBA_S3TC_DXT5_EXT
void PrintSelf(ostream &os, vtkIndent indent) override
static Vector3D offset
void RenderRGBAShadeFP(vtkRenderer *ren, vtkVolume *vol)
#define MITK_WARN
Definition: mitkLogMacros.h:23
int IsRenderSupported(vtkRenderer *ren, vtkVolumeProperty *) override
void GetLightInformation(vtkRenderer *ren, vtkVolume *vol, GLfloat lightDirection[2][4], GLfloat lightDiffuseColor[2][4], GLfloat lightSpecularColor[2][4], GLfloat halfwayVector[2][4], GLfloat *ambient)
vtkStandardNewMacro(vtkMitkOpenGLVolumeTextureMapper3D)
void Normalize(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer im2, mitk::Image::Pointer mask1, std::string output)
Definition: CLBrainMask.cpp:40
virtual void RenderFP(vtkRenderer *ren, vtkVolume *vol)
virtual void Render(vtkRenderer *ren, vtkVolume *vol) override
vtkRenderer * GetVtkRenderer() const
void RenderOneIndependentShadeFP(vtkRenderer *ren, vtkVolume *vol)
void vtkVolumeTextureMapper3DComputeRGBA(unsigned char *dataPtr, vtkMitkVolumeTextureMapper3D *me, GLuint volume1, GLuint volume2)
void SetupRGBATextures(vtkRenderer *ren, vtkVolume *vol)
void SetupOneIndependentTextures(vtkRenderer *ren, vtkVolume *vol)
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.