1 /*=========================================================================
2 Program:   Visualization Toolkit
3 Module:    TestGPURayCastClippingUserTransform.cxx
4 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
5 All rights reserved.
6 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
7 This software is distributed WITHOUT ANY WARRANTY; without even
8 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 PURPOSE.  See the above copyright notice for more information.
10 =========================================================================*/
11 
12 // Description
13 // This test creates a vtkImageData with two components.
14 // The data is volume rendered considering the two components as independent.
15 #include <fstream>
16 #include <iostream>
17 
18 #include "vtkActor.h"
19 #include "vtkCamera.h"
20 #include "vtkColorTransferFunction.h"
21 #include "vtkCommand.h"
22 #include "vtkGPUVolumeRayCastMapper.h"
23 #include "vtkImageData.h"
24 #include "vtkImageReader2.h"
25 #include "vtkInteractorStyleImage.h"
26 #include "vtkInteractorStyleTrackballCamera.h"
27 #include "vtkMatrix4x4.h"
28 #include "vtkNew.h"
29 #include "vtkOutlineFilter.h"
30 #include "vtkPiecewiseFunction.h"
31 #include "vtkPlane.h"
32 #include "vtkPlaneCollection.h"
33 #include "vtkPointData.h"
34 #include "vtkPolyDataMapper.h"
35 #include "vtkRenderWindow.h"
36 #include "vtkRenderWindowInteractor.h"
37 #include "vtkRenderer.h"
38 #include "vtkUnsignedShortArray.h"
39 #include "vtkVolume.h"
40 #include "vtkVolumeProperty.h"
41 #include "vtksys/FStream.hxx"
42 
43 #include "vtkRegressionTestImage.h"
44 #include "vtkTestUtilities.h"
45 
ComputeNormal(double * reference,bool flipSign)46 static double* ComputeNormal(double* reference, bool flipSign)
47 {
48   double* normal = new double[3];
49   if (flipSign)
50   {
51     normal[0] = -reference[0];
52     normal[1] = -reference[1];
53     normal[2] = -reference[2];
54   }
55   else
56   {
57     normal[0] = reference[0];
58     normal[1] = reference[1];
59     normal[2] = reference[2];
60   }
61 
62   return normal;
63 }
64 
ComputeOrigin(double * focalPoint,double * reference,double distance,bool flipSign)65 static double* ComputeOrigin(double* focalPoint, double* reference, double distance, bool flipSign)
66 {
67 
68   double* origin = new double[3];
69   if (flipSign)
70   {
71     origin[0] = focalPoint[0] - distance * reference[0];
72     origin[1] = focalPoint[1] - distance * reference[1];
73     origin[2] = focalPoint[2] - distance * reference[2];
74   }
75   else
76   {
77     origin[0] = focalPoint[0] + distance * reference[0];
78     origin[1] = focalPoint[1] + distance * reference[1];
79     origin[2] = focalPoint[2] + distance * reference[2];
80   }
81 
82   return origin;
83 }
84 
UpdateFrontClippingPlane(vtkPlane * frontClippingPlane,double * normal,double * focalPoint,double slabThickness)85 static void UpdateFrontClippingPlane(
86   vtkPlane* frontClippingPlane, double* normal, double* focalPoint, double slabThickness)
87 {
88   // The front plane is the start of ray cast
89   // The front normal should be in the same direction as the camera
90   // direction (opposite to the plane facing direction)
91   double* frontNormal = ComputeNormal(normal, true);
92 
93   // Get the origin of the front clipping plane
94   double halfSlabThickness = slabThickness / 2;
95   double* frontOrigin = ComputeOrigin(focalPoint, normal, halfSlabThickness, false);
96 
97   // Set the normal and origin of the front clipping plane
98   frontClippingPlane->SetNormal(frontNormal);
99   frontClippingPlane->SetOrigin(frontOrigin);
100 
101   delete[] frontNormal;
102   delete[] frontOrigin;
103 }
104 
UpdateRearClippingPlane(vtkPlane * rearClippingPlane,double * normal,double * focalPoint,double slabThickness)105 static void UpdateRearClippingPlane(
106   vtkPlane* rearClippingPlane, double* normal, double* focalPoint, double slabThickness)
107 {
108 
109   // The rear normal is the end of ray cast
110   // The rear normal should be in the opposite direction to the
111   // camera direction (same as the plane facing direction)
112   double* rearNormal = ComputeNormal(normal, false);
113 
114   // Get the origin of the rear clipping plane
115   double halfSlabThickness = slabThickness / 2;
116   double* rearOrigin = ComputeOrigin(focalPoint, normal, halfSlabThickness, true);
117 
118   // Set the normal and origin of the rear clipping plane
119   rearClippingPlane->SetNormal(rearNormal);
120   rearClippingPlane->SetOrigin(rearOrigin);
121 
122   delete[] rearNormal;
123   delete[] rearOrigin;
124 }
125 
126 class vtkInteractorStyleCallback : public vtkCommand
127 {
128 public:
New()129   static vtkInteractorStyleCallback* New() { return new vtkInteractorStyleCallback; }
130 
Execute(vtkObject * caller,unsigned long,void *)131   void Execute(vtkObject* caller, unsigned long, void*) override
132   {
133     vtkInteractorStyle* style = reinterpret_cast<vtkInteractorStyle*>(caller);
134 
135     vtkCamera* camera = style->GetCurrentRenderer()->GetActiveCamera();
136     // vtkCamera *camera = reinterpret_cast<vtkCamera*>(caller);
137 
138     // Get the normal and focal point of the camera
139     double* normal = camera->GetViewPlaneNormal();
140     double* focalPoint = camera->GetFocalPoint();
141 
142     // Fixed slab thickness
143     slabThickness = 3.0;
144     UpdateFrontClippingPlane(frontClippingPlane, normal, focalPoint, slabThickness);
145     UpdateRearClippingPlane(rearClippingPlane, normal, focalPoint, slabThickness);
146   }
147 
148   vtkInteractorStyleCallback() = default;
149 
SetFrontClippingPlane(vtkPlane * fcPlane)150   void SetFrontClippingPlane(vtkPlane* fcPlane) { this->frontClippingPlane = fcPlane; }
151 
SetRearClippingPlane(vtkPlane * rcPlane)152   void SetRearClippingPlane(vtkPlane* rcPlane) { this->rearClippingPlane = rcPlane; }
153 
154   double slabThickness;
155   vtkPlane* frontClippingPlane;
156   vtkPlane* rearClippingPlane;
157 };
158 
TestGPURayCastClippingUserTransform(int argc,char * argv[])159 int TestGPURayCastClippingUserTransform(int argc, char* argv[])
160 {
161   int width = 256;
162   int height = 256;
163   int depth = 148;
164   double spacing[3] = { 1.4844, 1.4844, 1.2 };
165 
166   // Read the image
167   std::streampos size;
168   char* memblock;
169 
170   char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/MagnitudeImage_256x256x148");
171 
172   vtksys::ifstream file(fname, ios::in | ios::binary | ios::ate);
173   if (file.is_open())
174   {
175     size = file.tellg();
176     memblock = new char[size];
177     file.seekg(0, ios::beg);
178     file.read(memblock, size);
179     file.close();
180   }
181   else
182   {
183     cout << "Unable to open file";
184     return 1;
185   }
186 
187   // Convert to short
188   unsigned short* shortData = new unsigned short[size / 2];
189   int idx = 0;
190   int idx2 = 0;
191   for (idx = 0; idx < size / 2; idx++)
192   {
193     idx2 = idx * 2;
194     shortData[idx] = (short)(((memblock[idx2] & 0xFF) << 8) | (memblock[idx2 + 1] & 0xFF));
195   }
196 
197   //
198   int volumeSizeInSlice = width * height * depth;
199   vtkNew<vtkUnsignedShortArray> dataArrayMag;
200   dataArrayMag->Allocate(volumeSizeInSlice, 0);
201   dataArrayMag->SetNumberOfComponents(1);
202   dataArrayMag->SetNumberOfTuples(volumeSizeInSlice);
203   dataArrayMag->SetArray(shortData, volumeSizeInSlice, 1);
204 
205   vtkNew<vtkImageData> imageData;
206   imageData->SetDimensions(width, height, depth);
207   imageData->SetSpacing(spacing);
208   imageData->GetPointData()->SetScalars(dataArrayMag);
209 
210   // Create a clipping plane
211   vtkNew<vtkPlane> frontClippingPlane;
212   vtkNew<vtkPlane> rearClippingPlane;
213 
214   // Create a clipping plane collection
215   vtkNew<vtkPlaneCollection> clippingPlaneCollection;
216   clippingPlaneCollection->AddItem(frontClippingPlane);
217   clippingPlaneCollection->AddItem(rearClippingPlane);
218 
219   // Create a mapper
220   vtkNew<vtkGPUVolumeRayCastMapper> volumeMapper;
221   // volumeMapper->SetInputConnection(reader->GetOutputPort());
222   volumeMapper->SetInputData(imageData);
223   volumeMapper->SetBlendModeToMaximumIntensity();
224   volumeMapper->AutoAdjustSampleDistancesOff();
225   volumeMapper->SetSampleDistance(1.0);
226   volumeMapper->SetImageSampleDistance(1.0);
227   volumeMapper->SetClippingPlanes(clippingPlaneCollection);
228 
229   // Create volume scale opacity
230   vtkNew<vtkPiecewiseFunction> volumeScalarOpacity;
231   volumeScalarOpacity->AddPoint(0, 0.0);
232   volumeScalarOpacity->AddPoint(32767, 1.0);
233   volumeScalarOpacity->ClampingOn();
234 
235   // Create a property
236   vtkNew<vtkVolumeProperty> volumeProperty;
237   volumeProperty->SetInterpolationTypeToLinear();
238   volumeProperty->ShadeOff();
239   volumeProperty->SetAmbient(1.0);
240   volumeProperty->SetDiffuse(0.0);
241   volumeProperty->SetSpecular(0.0);
242   volumeProperty->IndependentComponentsOn();
243   volumeProperty->SetScalarOpacity(volumeScalarOpacity);
244   volumeProperty->SetColor(volumeScalarOpacity);
245 
246   // Create a volume
247   vtkNew<vtkVolume> volume;
248   volume->SetMapper(volumeMapper);
249   volume->SetProperty(volumeProperty);
250   volume->PickableOff();
251 
252   // Rotate the blue props
253   double rowVector[3] = { 0.0, 0.0, -1.0 };
254   double columnVector[3] = { 1.0, 0.0, 0.0 };
255   double normalVector[3];
256   vtkMath::Cross(rowVector, columnVector, normalVector);
257   double position[3] = { 0.0, 0.0, 0.0 };
258 
259   vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
260   matrix->Identity();
261   matrix->SetElement(0, 0, rowVector[0]);
262   matrix->SetElement(0, 1, rowVector[1]);
263   matrix->SetElement(0, 2, rowVector[2]);
264   matrix->SetElement(0, 3, position[0]);
265   matrix->SetElement(1, 0, columnVector[0]);
266   matrix->SetElement(1, 1, columnVector[1]);
267   matrix->SetElement(1, 2, columnVector[2]);
268   matrix->SetElement(1, 3, position[1]);
269   matrix->SetElement(2, 0, normalVector[0]);
270   matrix->SetElement(2, 1, normalVector[1]);
271   matrix->SetElement(2, 2, normalVector[2]);
272   matrix->SetElement(2, 3, position[2]);
273 
274   volume->SetUserMatrix(matrix);
275 
276   // Create a outline filter
277   vtkNew<vtkOutlineFilter> outlineFilter;
278   outlineFilter->SetInputData(imageData);
279 
280   // Create an outline mapper
281   vtkNew<vtkPolyDataMapper> outlineMapper;
282   outlineMapper->SetInputConnection(outlineFilter->GetOutputPort());
283 
284   vtkNew<vtkActor> outline;
285   outline->SetMapper(outlineMapper);
286   outline->PickableOff();
287 
288   // Create a renderer
289   vtkNew<vtkRenderer> ren;
290   ren->AddViewProp(volume);
291   ren->AddViewProp(outline);
292 
293   // Get the center of volume
294   double* center = volume->GetCenter();
295   double cameraFocal[3];
296   cameraFocal[0] = center[0];
297   cameraFocal[1] = center[1];
298   cameraFocal[2] = center[2];
299 
300   double cameraViewUp[3] = { 0.00, -1.00, 0.00 };
301   double cameraNormal[3] = { 0.00, 0.00, -1.00 };
302   double cameraDistance = 1000.0;
303 
304   double cameraPosition[3];
305   cameraPosition[0] = cameraFocal[0] + cameraDistance * cameraNormal[0];
306   cameraPosition[1] = cameraFocal[1] + cameraDistance * cameraNormal[1];
307   cameraPosition[2] = cameraFocal[2] + cameraDistance * cameraNormal[2];
308 
309   // Update clipping planes
310   UpdateFrontClippingPlane(frontClippingPlane, cameraNormal, cameraFocal, 3.0);
311   UpdateRearClippingPlane(rearClippingPlane, cameraNormal, cameraFocal, 3.0);
312 
313   // Get the active camera
314   vtkCamera* camera = ren->GetActiveCamera();
315   camera->ParallelProjectionOn();
316   camera->SetParallelScale(250);
317   camera->SetPosition(cameraPosition);
318   camera->SetFocalPoint(cameraFocal);
319   camera->SetViewUp(cameraViewUp);
320 
321   // Create a render window
322   vtkNew<vtkRenderWindow> renWin;
323   renWin->SetSize(500, 500);
324   renWin->AddRenderer(ren);
325 
326   // Create a style
327   vtkNew<vtkInteractorStyleImage> style;
328   style->SetInteractionModeToImage3D();
329 
330   // Create a interactor style callback
331   vtkNew<vtkInteractorStyleCallback> interactorStyleCallback;
332   interactorStyleCallback->frontClippingPlane = frontClippingPlane;
333   interactorStyleCallback->rearClippingPlane = rearClippingPlane;
334   style->AddObserver(vtkCommand::InteractionEvent, interactorStyleCallback);
335 
336   // Create an interactor
337   vtkNew<vtkRenderWindowInteractor> iren;
338   iren->SetInteractorStyle(style);
339   iren->SetRenderWindow(renWin);
340 
341   // Start
342   iren->Initialize();
343   renWin->Render();
344 
345   int retVal = vtkRegressionTestImageThreshold(renWin, 70);
346 
347   if (retVal == vtkRegressionTester::DO_INTERACTOR)
348   {
349     iren->Start();
350   }
351 
352   delete[] memblock;
353   delete[] shortData;
354   delete[] fname;
355 
356   return !retVal;
357 }
358