Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkExtractSliceFilterTest.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 (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #include <itkImage.h>
14 #include <itkImageRegionIterator.h>
15 #include <mitkExtractSliceFilter.h>
16 #include <mitkIOUtil.h>
17 #include <mitkITKImageImport.h>
18 #include <mitkImageAccessByItk.h>
19 #include <mitkImageCast.h>
21 #include <mitkInteractionConst.h>
22 #include <mitkNumericTypes.h>
23 #include <mitkRotationOperation.h>
25 #include <mitkTestingMacros.h>
26 
27 #include <cstdlib>
28 #include <ctime>
29 #include <cmath>
30 
31 #include <mitkGeometry3D.h>
32 
33 #include <vtkActor.h>
34 #include <vtkImageActor.h>
35 #include <vtkImageData.h>
36 #include <vtkImageMapToColors.h>
37 #include <vtkImageMapper.h>
38 #include <vtkImageReslice.h>
39 #include <vtkInteractorStyleImage.h>
40 #include <vtkLookupTable.h>
41 #include <vtkPlaneSource.h>
42 #include <vtkPolyDataMapper.h>
43 #include <vtkRenderWindow.h>
44 #include <vtkRenderWindowInteractor.h>
45 #include <vtkRenderer.h>
46 #include <vtkSmartPointer.h>
47 #include <vtkTexture.h>
48 
49 // use this to create the test volume on the fly
50 #define CREATE_VOLUME
51 
52 // use this to save the created volume
53 //#define SAVE_VOLUME
54 
55 // use this to calculate the error from the sphere mathematical model to our pixel based one
56 //#define CALC_TESTFAILURE_DEVIATION
57 
58 // use this to render an oblique slice through a specified image
59 //#define SHOW_SLICE_IN_RENDER_WINDOW
60 
61 // use this to have infos printed in mbilog
62 //#define EXTRACTOR_DEBUG
63 
64 /*these are the deviations calculated by the function CalcTestFailureDeviation (see for details)*/
65 #define Testfailure_Deviation_Mean_128 0.853842
66 
67 #define Testfailure_Deviation_Volume_128 0.145184
68 #define Testfailure_Deviation_Diameter_128 1.5625
69 
70 #define Testfailure_Deviation_Mean_256 0.397693
71 
72 #define Testfailure_Deviation_Volume_256 0.0141357
73 #define Testfailure_Deviation_Diameter_256 0.78125
74 
75 #define Testfailure_Deviation_Mean_512 0.205277
76 
77 #define Testfailure_Deviation_Volume_512 0.01993
78 #define Testfailure_Deviation_Diameter_512 0.390625
79 
80 class mitkExtractSliceFilterTestClass
81 {
82 public:
83  static void TestSlice(mitk::PlaneGeometry *planeGeometry, std::string testname)
84  {
85  TestPlane = planeGeometry;
86  TestName = testname;
87 
88  mitk::ScalarType centerCoordValue = TestvolumeSize / 2.0;
89  mitk::ScalarType center[3] = {centerCoordValue, centerCoordValue, centerCoordValue};
90  mitk::Point3D centerIndex(center);
91 
92  double radius = TestvolumeSize / 4.0;
93  if (TestPlane->Distance(centerIndex) >= radius)
94  return; // outside sphere
95 
96  // feed ExtractSliceFilter
98 
99  slicer->SetInput(TestVolume);
100  slicer->SetWorldGeometry(TestPlane);
101  slicer->Update();
102 
103  MITK_TEST_CONDITION_REQUIRED(slicer->GetOutput() != nullptr, "Extractor returned a slice");
104 
105  mitk::Image::Pointer reslicedImage = slicer->GetOutput();
106 
107  AccessFixedDimensionByItk(reslicedImage, TestSphereRadiusByItk, 2);
108  AccessFixedDimensionByItk(reslicedImage, TestSphereAreaByItk, 2);
109 
110  /*
111  double devArea, devDiameter;
112  if(TestvolumeSize == 128.0){ devArea = Testfailure_Deviation_Volume_128; devDiameter =
113  Testfailure_Deviation_Diameter_128; }
114  else if(TestvolumeSize == 256.0){devArea = Testfailure_Deviation_Volume_256; devDiameter =
115  Testfailure_Deviation_Diameter_256;}
116  else if (TestvolumeSize == 512.0){devArea = Testfailure_Deviation_Volume_512; devDiameter =
117  Testfailure_Deviation_Diameter_512;}
118  else{devArea = Testfailure_Deviation_Volume_128; devDiameter = Testfailure_Deviation_Diameter_128;}
119  */
120 
121  std::string areatestName = TestName.append(" area");
122  std::string diametertestName = TestName.append(" testing diameter");
123 
124  // TODO think about the deviation, 1% makes no sense at all
125  MITK_TEST_CONDITION(std::abs(100 - testResults.percentageAreaCalcToPixel) < 1, areatestName);
126  MITK_TEST_CONDITION(std::abs(100 - testResults.percentageRadiusToPixel) < 1, diametertestName);
127 
128 #ifdef EXTRACTOR_DEBUG
129  MITK_INFO << TestName << " >>> "
130  << "planeDistanceToSphereCenter: " << testResults.planeDistanceToSphereCenter;
131  MITK_INFO << "area in pixels: " << testResults.areaInPixel << " <-> area in mm: " << testResults.areaCalculated
132  << " = " << testResults.percentageAreaCalcToPixel << "%";
133 
134  MITK_INFO << "calculated diameter: " << testResults.diameterCalculated
135  << " <-> diameter in mm: " << testResults.diameterInMM
136  << " <-> diameter in pixel: " << testResults.diameterInPixel << " = "
137  << testResults.percentageRadiusToPixel << "%";
138 #endif
139  }
140 
141  /*
142  * get the radius of the slice of a sphere based on pixel distance from edge to edge of the circle.
143  */
144  template <typename TPixel, unsigned int VImageDimension>
145  static void TestSphereRadiusByItk(itk::Image<TPixel, VImageDimension> *inputImage)
146  {
147  typedef itk::Image<TPixel, VImageDimension> InputImageType;
148 
149  // set the index to the middle of the image's edge at x and y axis
150  typename InputImageType::IndexType currentIndexX;
151  currentIndexX[0] = (int)(TestvolumeSize / 2.0);
152  currentIndexX[1] = 0;
153 
154  typename InputImageType::IndexType currentIndexY;
155  currentIndexY[0] = 0;
156  currentIndexY[1] = (int)(TestvolumeSize / 2.0);
157 
158  // remember the last pixel value
159  double lastValueX = inputImage->GetPixel(currentIndexX);
160  double lastValueY = inputImage->GetPixel(currentIndexY);
161 
162  // storage for the index marks
163  std::vector<typename InputImageType::IndexType> indicesX;
164  std::vector<typename InputImageType::IndexType> indicesY;
165 
166  /*Get four indices on the edge of the circle*/
167  while (currentIndexX[1] < TestvolumeSize && currentIndexX[0] < TestvolumeSize)
168  {
169  // move x direction
170  currentIndexX[1] += 1;
171 
172  // move y direction
173  currentIndexY[0] += 1;
174 
175  if (inputImage->GetPixel(currentIndexX) > lastValueX)
176  {
177  // mark the current index
178  typename InputImageType::IndexType markIndex;
179  markIndex[0] = currentIndexX[0];
180  markIndex[1] = currentIndexX[1];
181 
182  indicesX.push_back(markIndex);
183  }
184  else if (inputImage->GetPixel(currentIndexX) < lastValueX)
185  {
186  // mark the current index
187  typename InputImageType::IndexType markIndex;
188  markIndex[0] = currentIndexX[0];
189  markIndex[1] = currentIndexX[1] - 1; // value inside the sphere
190 
191  indicesX.push_back(markIndex);
192  }
193 
194  if (inputImage->GetPixel(currentIndexY) > lastValueY)
195  {
196  // mark the current index
197  typename InputImageType::IndexType markIndex;
198  markIndex[0] = currentIndexY[0];
199  markIndex[1] = currentIndexY[1];
200 
201  indicesY.push_back(markIndex);
202  }
203  else if (inputImage->GetPixel(currentIndexY) < lastValueY)
204  {
205  // mark the current index
206  typename InputImageType::IndexType markIndex;
207  markIndex[0] = currentIndexY[0];
208  markIndex[1] = currentIndexY[1] - 1; // value inside the sphere
209 
210  indicesY.push_back(markIndex);
211  }
212 
213  // found both marks?
214  if (indicesX.size() == 2 && indicesY.size() == 2)
215  break;
216 
217  // the new 'last' values
218  lastValueX = inputImage->GetPixel(currentIndexX);
219  lastValueY = inputImage->GetPixel(currentIndexY);
220  }
221 
222  /*
223  *If we are here we found the four marks on the edge of the circle.
224  *For the case our plane is rotated and shifted, we have to calculate the center of the circle,
225  *else the center is the intersection of both straight lines between the marks.
226  *When we have the center, the diameter of the circle will be checked to the reference value(math!).
227  */
228 
229  // each distance from the first mark of each direction to the center of the straight line between the marks
230  double distanceToCenterX = std::abs(indicesX[0][1] - indicesX[1][1]) / 2.0;
231  // double distanceToCenterY = std::abs(indicesY[0][0] - indicesY[1][0]) / 2.0;
232 
233  // the center of the straight lines
234  typename InputImageType::IndexType centerX;
235  // typename InputImageType::IndexType centerY;
236 
237  centerX[0] = indicesX[0][0];
238  centerX[1] = indicesX[0][1] + distanceToCenterX;
239  // TODO think about implicit cast to int. this is not the real center of the image, which could be between two
240  // pixels
241 
242  // centerY[0] = indicesY[0][0] + distanceToCenterY;
243  // centerY[1] = inidcesY[0][1];
244 
245  typename InputImageType::IndexType currentIndex(centerX);
246  lastValueX = inputImage->GetPixel(currentIndex);
247 
248  long sumpixels = 0;
249 
250  std::vector<typename InputImageType::IndexType> diameterIndices;
251 
252  // move up
253  while (currentIndex[1] < TestvolumeSize)
254  {
255  currentIndex[1] += 1;
256 
257  if (inputImage->GetPixel(currentIndex) != lastValueX)
258  {
259  typename InputImageType::IndexType markIndex;
260  markIndex[0] = currentIndex[0];
261  markIndex[1] = currentIndex[1] - 1;
262 
263  diameterIndices.push_back(markIndex);
264  break;
265  }
266 
267  sumpixels++;
268  lastValueX = inputImage->GetPixel(currentIndex);
269  }
270 
271  currentIndex[1] -= sumpixels; // move back to center to go in the other direction
272  lastValueX = inputImage->GetPixel(currentIndex);
273 
274  // move down
275  while (currentIndex[1] >= 0)
276  {
277  currentIndex[1] -= 1;
278 
279  if (inputImage->GetPixel(currentIndex) != lastValueX)
280  {
281  typename InputImageType::IndexType markIndex;
282  markIndex[0] = currentIndex[0];
283  markIndex[1] = currentIndex[1]; // outside sphere because we want to calculate the distance from edge to edge
284 
285  diameterIndices.push_back(markIndex);
286  break;
287  }
288 
289  sumpixels++;
290  lastValueX = inputImage->GetPixel(currentIndex);
291  }
292 
293  /*
294  *Now sumpixels should be the apromximate diameter of the circle. This is checked with the calculated diameter from
295  *the plane transformation(math).
296  */
297  mitk::Point3D volumeCenter;
298  volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0;
299 
300  double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter);
301 
302  double sphereRadius = TestvolumeSize / 4.0;
303 
304  // calculate the radius of the circle cut from the sphere by the plane
305  double diameter = 2.0 * std::sqrt(std::pow(sphereRadius, 2) - std::pow(planeDistanceToSphereCenter, 2));
306 
307  double percentageRadiusToPixel = 100 / diameter * sumpixels;
308 
309  /*
310  *calculate the radius in mm by the both marks of the center line by using the world coordinates
311  */
312  // get the points as 3D coordinates
313  mitk::Vector3D diameterPointRight, diameterPointLeft;
314 
315  diameterPointRight[2] = diameterPointLeft[2] = 0.0;
316 
317  diameterPointLeft[0] = diameterIndices[0][0];
318  diameterPointLeft[1] = diameterIndices[0][1];
319 
320  diameterPointRight[0] = diameterIndices[1][0];
321  diameterPointRight[1] = diameterIndices[1][1];
322 
323  // transform to worldcoordinates
324  TestVolume->GetGeometry()->IndexToWorld(diameterPointLeft, diameterPointLeft);
325  TestVolume->GetGeometry()->IndexToWorld(diameterPointRight, diameterPointRight);
326 
327  // euklidian distance
328  double diameterInMM = ((diameterPointLeft * -1.0) + diameterPointRight).GetNorm();
329 
330  testResults.diameterInMM = diameterInMM;
331  testResults.diameterCalculated = diameter;
332  testResults.diameterInPixel = sumpixels;
333  testResults.percentageRadiusToPixel = percentageRadiusToPixel;
334  testResults.planeDistanceToSphereCenter = planeDistanceToSphereCenter;
335  }
336 
337  /*brute force the area pixel by pixel*/
338  template <typename TPixel, unsigned int VImageDimension>
339  static void TestSphereAreaByItk(itk::Image<TPixel, VImageDimension> *inputImage)
340  {
341  typedef itk::Image<TPixel, VImageDimension> InputImageType;
342  typedef itk::ImageRegionConstIterator<InputImageType> ImageIterator;
343 
344  ImageIterator imageIterator(inputImage, inputImage->GetLargestPossibleRegion());
345  imageIterator.GoToBegin();
346 
347  int sumPixelsInArea = 0;
348 
349  while (!imageIterator.IsAtEnd())
350  {
351  if (inputImage->GetPixel(imageIterator.GetIndex()) == pixelValueSet)
352  sumPixelsInArea++;
353  ++imageIterator;
354  }
355 
356  mitk::Point3D volumeCenter;
357  volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0;
358 
359  double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter);
360 
361  double sphereRadius = TestvolumeSize / 4.0;
362 
363  // calculate the radius of the circle cut from the sphere by the plane
364  double radius = std::sqrt(std::pow(sphereRadius, 2) - std::pow(planeDistanceToSphereCenter, 2));
365 
366  double areaInMM = 3.14159265358979 * std::pow(radius, 2);
367 
368  testResults.areaCalculated = areaInMM;
369  testResults.areaInPixel = sumPixelsInArea;
370  testResults.percentageAreaCalcToPixel = 100 / areaInMM * sumPixelsInArea;
371  }
372 
373  /*
374  * random a voxel. define plane through this voxel. reslice at the plane. compare the pixel vaues of the voxel
375  * in the volume with the pixel value in the resliced image.
376  * there are some indice shifting problems which causes the test to fail for oblique planes. seems like the chosen
377  * worldcoordinate is not corrresponding to the index in the 2D image. and so the pixel values are not the same as
378  * expected.
379  */
380  static void PixelvalueBasedTest()
381  {
382  /* setup itk image */
383  typedef itk::Image<unsigned short, 3> ImageType;
384 
385  typedef itk::ImageRegionConstIterator<ImageType> ImageIterator;
386 
387  ImageType::Pointer image = ImageType::New();
388 
389  ImageType::IndexType start;
390  start[0] = start[1] = start[2] = 0;
391 
392  ImageType::SizeType size;
393  size[0] = size[1] = size[2] = 32;
394 
395  ImageType::RegionType imgRegion;
396  imgRegion.SetSize(size);
397  imgRegion.SetIndex(start);
398 
399  image->SetRegions(imgRegion);
400  image->SetSpacing(1.0);
401  image->Allocate();
402 
403  ImageIterator imageIterator(image, image->GetLargestPossibleRegion());
404  imageIterator.GoToBegin();
405 
406  unsigned short pixelValue = 0;
407 
408  // fill the image with distinct values
409  while (!imageIterator.IsAtEnd())
410  {
411  image->SetPixel(imageIterator.GetIndex(), pixelValue);
412  ++imageIterator;
413  ++pixelValue;
414  }
415  /* end setup itk image */
416 
417  mitk::Image::Pointer imageInMitk;
418  CastToMitkImage(image, imageInMitk);
419 
420  /*mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New();
421  writer->SetInput(imageInMitk);
422  std::string file = "C:\\Users\\schroedt\\Desktop\\cube.nrrd";
423  writer->SetFileName(file);
424  writer->Update();*/
425 
426  PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Frontal);
427  PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Sagittal);
428  PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Axial);
429  }
430 
431  static void PixelvalueBasedTestByPlane(mitk::Image *imageInMitk, mitk::PlaneGeometry::PlaneOrientation orientation)
432  {
433  typedef itk::Image<unsigned short, 3> ImageType;
434 
435  // set the seed of the rand function
436  srand((unsigned)time(nullptr));
437 
438  /* setup a random orthogonal plane */
439  int sliceindex = 17; // rand() % 32;
440  bool isFrontside = true;
441  bool isRotated = false;
442 
443  if (orientation == mitk::PlaneGeometry::Axial)
444  {
445  /*isFrontside = false;
446  isRotated = true;*/
447  }
448 
450 
451  plane->InitializeStandardPlane(imageInMitk->GetGeometry(), orientation, sliceindex, isFrontside, isRotated);
452 
453  mitk::Point3D origin = plane->GetOrigin();
454  mitk::Vector3D normal;
455  normal = plane->GetNormal();
456 
457  normal.Normalize();
458 
459  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
460 
461  plane->SetOrigin(origin);
462 
463  // we dont need this any more, because we are only testing orthogonal planes
464  /*mitk::Vector3D rotationVector;
465  rotationVector[0] = randFloat();
466  rotationVector[1] = randFloat();
467  rotationVector[2] = randFloat();
468 
469 
470  float degree = randFloat() * 180.0;
471 
472  mitk::RotationOperation* op = new mitk::RotationOperation(mitk::OpROTATE, plane->GetCenter(), rotationVector,
473  degree);
474  plane->ExecuteOperation(op);
475  delete op;*/
476 
477  /* end setup plane */
478 
479  /* define a point in the 3D volume.
480  * add the two axis vectors of the plane (each multiplied with a
481  * random number) to the origin. now the two random numbers
482  * become our index coordinates in the 2D image, because the
483  * length of the axis vectors is 1.
484  */
485  mitk::Point3D planeOrigin = plane->GetOrigin();
486  mitk::Vector3D axis0, axis1;
487  axis0 = plane->GetAxisVector(0);
488  axis1 = plane->GetAxisVector(1);
489  axis0.Normalize();
490  axis1.Normalize();
491 
492  unsigned char n1 = 7; // rand() % 32;
493  unsigned char n2 = 13; // rand() % 32;
494 
495  mitk::Point3D testPoint3DInWorld;
496  testPoint3DInWorld = planeOrigin + (axis0 * n1) + (axis1 * n2);
497 
498  // get the index of the point in the 3D volume
499  ImageType::IndexType testPoint3DInIndex;
500  imageInMitk->GetGeometry()->WorldToIndex(testPoint3DInWorld, testPoint3DInIndex);
501 
502  itk::Index<3> testPoint2DInIndex;
503 
504  /* end define a point in the 3D volume.*/
505 
506  // do reslicing at the plane
508  slicer->SetInput(imageInMitk);
509  slicer->SetWorldGeometry(plane);
510 
511  slicer->Update();
512 
513  mitk::Image::Pointer slice = slicer->GetOutput();
514 
515  // Get TestPoiont3D as Index in Slice
516  slice->GetGeometry()->WorldToIndex(testPoint3DInWorld, testPoint2DInIndex);
517 
518  mitk::Point3D p, sliceIndexToWorld, imageIndexToWorld;
519  p[0] = testPoint2DInIndex[0];
520  p[1] = testPoint2DInIndex[1];
521  p[2] = testPoint2DInIndex[2];
522  slice->GetGeometry()->IndexToWorld(p, sliceIndexToWorld);
523 
524  p[0] = testPoint3DInIndex[0];
525  p[1] = testPoint3DInIndex[1];
526  p[2] = testPoint3DInIndex[2];
527  imageInMitk->GetGeometry()->IndexToWorld(p, imageIndexToWorld);
528 
529  itk::Index<2> testPoint2DIn2DIndex;
530  testPoint2DIn2DIndex[0] = testPoint2DInIndex[0];
531  testPoint2DIn2DIndex[1] = testPoint2DInIndex[1];
532 
533  typedef mitk::ImagePixelReadAccessor<unsigned short, 3> VolumeReadAccessorType;
534  typedef mitk::ImagePixelReadAccessor<unsigned short, 2> SliceReadAccessorType;
535  VolumeReadAccessorType VolumeReadAccessor(imageInMitk);
536  SliceReadAccessorType SliceReadAccessor(slice);
537 
538  // compare the pixelvalues of the defined point in the 3D volume with the value of the resliced image
539  unsigned short valueAt3DVolume = VolumeReadAccessor.GetPixelByIndex(testPoint3DInIndex);
540  unsigned short valueAtSlice = SliceReadAccessor.GetPixelByIndex(testPoint2DIn2DIndex);
541 
542  // valueAt3DVolume == valueAtSlice is not always working. because of rounding errors
543  // indices are shifted
544  MITK_TEST_CONDITION(valueAt3DVolume == valueAtSlice, "comparing pixelvalues for orthogonal plane");
545 
546  vtkSmartPointer<vtkImageData> imageInVtk = imageInMitk->GetVtkImageData();
547  vtkSmartPointer<vtkImageData> sliceInVtk = slice->GetVtkImageData();
548 
549  double PixelvalueByMitkOutput = sliceInVtk->GetScalarComponentAsDouble(n1, n2, 0, 0);
550  // double valueVTKinImage = imageInVtk->GetScalarComponentAsDouble(testPoint3DInIndex[0], testPoint3DInIndex[1],
551  // testPoint3DInIndex[2], 0);
552 
553  /* Test that everything is working equally if vtkoutput is used instead of the default output
554  * from mitk ImageToImageFilter
555  */
557  slicerWithVtkOutput->SetInput(imageInMitk);
558  slicerWithVtkOutput->SetWorldGeometry(plane);
559  slicerWithVtkOutput->SetVtkOutputRequest(true);
560 
561  slicerWithVtkOutput->Update();
562  vtkSmartPointer<vtkImageData> vtkImageByVtkOutput = slicerWithVtkOutput->GetVtkOutput();
563  double PixelvalueByVtkOutput = vtkImageByVtkOutput->GetScalarComponentAsDouble(n1, n2, 0, 0);
564 
565  MITK_TEST_CONDITION(PixelvalueByMitkOutput == PixelvalueByVtkOutput,
566  "testing convertion of image output vtk->mitk by reslicer");
567 
568 /*================ mbilog outputs ===========================*/
569 #ifdef EXTRACTOR_DEBUG
570  MITK_INFO << "\n"
571  << "TESTINFO index: " << sliceindex << " orientation: " << orientation << " frontside: " << isFrontside
572  << " rotated: " << isRotated;
573  MITK_INFO << "\n"
574  << "slice index to world: " << sliceIndexToWorld;
575  MITK_INFO << "\n"
576  << "image index to world: " << imageIndexToWorld;
577 
578  MITK_INFO << "\n"
579  << "vtk: slice: " << PixelvalueByMitkOutput << ", image: " << valueVTKinImage;
580 
581  MITK_INFO << "\n"
582  << "testPoint3D InWorld" << testPoint3DInWorld << " is " << testPoint2DInIndex << " in 2D";
583  MITK_INFO << "\n"
584  << "randoms: " << ((int)n1) << ", " << ((int)n2);
585  MITK_INFO << "\n"
586  << "point is inside plane: " << plane->IsInside(testPoint3DInWorld)
587  << " and volume: " << imageInMitk->GetGeometry()->IsInside(testPoint3DInWorld);
588 
589  MITK_INFO << "\n"
590  << "volume idx: " << testPoint3DInIndex << " = " << valueAt3DVolume;
591  MITK_INFO << "\n"
592  << "volume world: " << testPoint3DInWorld << " = " << valueAt3DVolumeByWorld;
593  MITK_INFO << "\n"
594  << "slice idx: " << testPoint2DInIndex << " = " << valueAtSlice;
595 
596  itk::Index<3> curr;
597  curr[0] = curr[1] = curr[2] = 0;
598 
599  for (int i = 0; i < 32; ++i)
600  {
601  for (int j = 0; j < 32; ++j)
602  {
603  ++curr[1];
604  if (SliceReadAccessor.GetPixelByIndex(curr) == valueAt3DVolume)
605  {
606  MITK_INFO << "\n" << valueAt3DVolume << " MATCHED mitk " << curr;
607  }
608  }
609  curr[1] = 0;
610  ++curr[0];
611  }
612 
613  typedef itk::Image<unsigned short, 2> Image2DType;
614 
615  Image2DType::Pointer img = Image2DType::New();
616  CastToItkImage(slice, img);
617  typedef itk::ImageRegionConstIterator<Image2DType> Iterator2D;
618 
619  Iterator2D iter(img, img->GetLargestPossibleRegion());
620  iter.GoToBegin();
621  while (!iter.IsAtEnd())
622  {
623  if (img->GetPixel(iter.GetIndex()) == valueAt3DVolume)
624  MITK_INFO << "\n" << valueAt3DVolume << " MATCHED itk " << iter.GetIndex();
625 
626  ++iter;
627  }
628 #endif // EXTRACTOR_DEBUG
629  }
630 
631  /* random a float value */
632  static float randFloat()
633  {
634  return (((float)rand() + 1.0) / ((float)RAND_MAX + 1.0)) +
635  (((float)rand() + 1.0) / ((float)RAND_MAX + 1.0)) / ((float)RAND_MAX + 1.0);
636  }
637 
638  /* create a sphere with the size of the given testVolumeSize*/
639  static void InitializeTestVolume()
640  {
641 #ifdef CREATE_VOLUME
642 
643  // do sphere creation
644  ItkVolumeGeneration();
645 
646 #ifdef SAVE_VOLUME
647  // save in file
649  writer->SetInput(TestVolume);
650 
651  std::string file;
652 
653  std::ostringstream filename;
654  filename << "C:\\home\\schroedt\\MITK\\Modules\\ImageExtraction\\Testing\\Data\\sphere_";
655  filename << TestvolumeSize;
656  filename << ".nrrd";
657 
658  file = filename.str();
659 
660  writer->SetFileName(file);
661  writer->Update();
662 #endif // SAVE_VOLUME
663 
664 #endif
665 
666 #ifndef CREATE_VOLUME // read from file
667 
669 
670  std::string filename = locator->FindFile("sphere_512.nrrd.mhd", "Modules/ImageExtraction/Testing/Data");
671 
672  TestVolume = mitk::IOUtil::Load<mitk::Image>(filename);
673 
674 #endif
675 
676 #ifdef CALC_TESTFAILURE_DEVIATION
677  // get the TestFailureDeviation in %
678  AccessFixedDimensionByItk(TestVolume, CalcTestFailureDeviation, 3);
679 #endif
680  }
681 
682  // the test result of the sphere reslice
683  struct SliceProperties
684  {
685  double planeDistanceToSphereCenter;
686  double diameterInMM;
687  double diameterInPixel;
688  double diameterCalculated;
689  double percentageRadiusToPixel;
690  double areaCalculated;
691  double areaInPixel;
692  double percentageAreaCalcToPixel;
693  };
694 
695  static mitk::Image::Pointer TestVolume;
696  static double TestvolumeSize;
697  static mitk::PlaneGeometry::Pointer TestPlane;
698  static std::string TestName;
699  static unsigned char pixelValueSet;
700  static SliceProperties testResults;
701  static double TestFailureDeviation;
702 
703 private:
704  /*
705  * Generate a sphere with a radius of TestvolumeSize / 4.0
706  */
707  static void ItkVolumeGeneration()
708  {
709  typedef itk::Image<unsigned char, 3> TestVolumeType;
710 
711  typedef itk::ImageRegionConstIterator<TestVolumeType> ImageIterator;
712 
713  TestVolumeType::Pointer sphereImage = TestVolumeType::New();
714 
715  TestVolumeType::IndexType start;
716  start[0] = start[1] = start[2] = 0;
717 
718  TestVolumeType::SizeType size;
719  size[0] = size[1] = size[2] = TestvolumeSize;
720 
721  TestVolumeType::RegionType imgRegion;
722  imgRegion.SetSize(size);
723  imgRegion.SetIndex(start);
724 
725  sphereImage->SetRegions(imgRegion);
726  sphereImage->SetSpacing(1.0);
727  sphereImage->Allocate();
728 
729  sphereImage->FillBuffer(0);
730 
731  mitk::Vector3D center;
732  center[0] = center[1] = center[2] = TestvolumeSize / 2.0;
733 
734  double radius = TestvolumeSize / 4.0;
735 
736  double pixelValue = pixelValueSet;
737 
738  ImageIterator imageIterator(sphereImage, sphereImage->GetLargestPossibleRegion());
739  imageIterator.GoToBegin();
740 
741  mitk::Vector3D currentVoxelInIndex;
742 
743  while (!imageIterator.IsAtEnd())
744  {
745  currentVoxelInIndex[0] = imageIterator.GetIndex()[0];
746  currentVoxelInIndex[1] = imageIterator.GetIndex()[1];
747  currentVoxelInIndex[2] = imageIterator.GetIndex()[2];
748 
749  double distanceToCenter = (center + (currentVoxelInIndex * -1.0)).GetNorm();
750 
751  // if distance to center is smaller then the radius of the sphere
752  if (distanceToCenter < radius)
753  {
754  sphereImage->SetPixel(imageIterator.GetIndex(), pixelValue);
755  }
756 
757  ++imageIterator;
758  }
759 
760  CastToMitkImage(sphereImage, TestVolume);
761  }
762 
763  /* calculate the devation of the voxel object to the mathematical sphere object.
764  * this is use to make a statement about the accuracy of the resliced image, eg. the circle's diameter or area.
765  */
766  template <typename TPixel, unsigned int VImageDimension>
767  static void CalcTestFailureDeviation(itk::Image<TPixel, VImageDimension> *inputImage)
768  {
769  typedef itk::Image<TPixel, VImageDimension> InputImageType;
770  typedef itk::ImageRegionConstIterator<InputImageType> ImageIterator;
771 
772  ImageIterator iterator(inputImage, inputImage->GetLargestPossibleRegion());
773  iterator.GoToBegin();
774 
775  int volumeInPixel = 0;
776 
777  while (!iterator.IsAtEnd())
778  {
779  if (inputImage->GetPixel(iterator.GetIndex()) == pixelValueSet)
780  volumeInPixel++;
781  ++iterator;
782  }
783 
784  double diameter = TestvolumeSize / 2.0;
785  double volumeCalculated = (1.0 / 6.0) * 3.14159265358979 * std::pow(diameter, 3);
786 
787  double volumeDeviation = std::abs(100 - (100 / volumeCalculated * volumeInPixel));
788 
789  typename InputImageType::IndexType index;
790  index[0] = index[1] = TestvolumeSize / 2.0;
791  index[2] = 0;
792 
793  int sumpixels = 0;
794  while (index[2] < TestvolumeSize)
795  {
796  if (inputImage->GetPixel(index) == pixelValueSet)
797  sumpixels++;
798  index[2] += 1;
799  }
800 
801  double diameterDeviation = std::abs(100 - (100 / diameter * sumpixels));
802 #ifdef DEBUG
803  MITK_INFO << "volume deviation: " << volumeDeviation << " diameter deviation:" << diameterDeviation;
804 #endif
805  mitkExtractSliceFilterTestClass::TestFailureDeviation = (volumeDeviation + diameterDeviation) / 2.0;
806  }
807 };
808 /*================ #END class ================*/
809 
810 /*================#BEGIN Instanciation of members ================*/
811 mitk::Image::Pointer mitkExtractSliceFilterTestClass::TestVolume = mitk::Image::New();
812 double mitkExtractSliceFilterTestClass::TestvolumeSize = 256.0;
813 mitk::PlaneGeometry::Pointer mitkExtractSliceFilterTestClass::TestPlane = mitk::PlaneGeometry::New();
814 std::string mitkExtractSliceFilterTestClass::TestName = "";
815 unsigned char mitkExtractSliceFilterTestClass::pixelValueSet = 255;
816 mitkExtractSliceFilterTestClass::SliceProperties mitkExtractSliceFilterTestClass::testResults = {
817  -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
818 double mitkExtractSliceFilterTestClass::TestFailureDeviation = 0.0;
819 /*================ #END Instanciation of members ================*/
820 
821 /*================ #BEGIN test main ================*/
822 int mitkExtractSliceFilterTest(int /*argc*/, char * /*argv*/ [])
823 {
824  MITK_TEST_BEGIN("mitkExtractSliceFilterTest")
825 
826  // pixelvalue based testing
827  mitkExtractSliceFilterTestClass::PixelvalueBasedTest();
828 
829  // initialize sphere test volume
830  mitkExtractSliceFilterTestClass::InitializeTestVolume();
831 
832  mitk::Vector3D spacing = mitkExtractSliceFilterTestClass::TestVolume->GetGeometry()->GetSpacing();
833 
834  // the center of the sphere = center of image
835  double sphereCenter = mitkExtractSliceFilterTestClass::TestvolumeSize / 2.0;
836 
837  double planeSize = mitkExtractSliceFilterTestClass::TestvolumeSize;
838 
839  /* axial plane */
841  geometryAxial->InitializeStandardPlane(
842  planeSize, planeSize, spacing, mitk::PlaneGeometry::Axial, sphereCenter, false, true);
843  geometryAxial->ChangeImageGeometryConsideringOriginOffset(true);
844 
845  mitk::Point3D origin = geometryAxial->GetOrigin();
846  mitk::Vector3D normal;
847  normal = geometryAxial->GetNormal();
848 
849  normal.Normalize();
850 
851  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
852 
853  // geometryAxial->SetOrigin(origin);
854 
855  mitkExtractSliceFilterTestClass::TestSlice(geometryAxial, "Testing axial plane");
856  /* end axial plane */
857 
858  /* sagittal plane */
860  geometrySagital->InitializeStandardPlane(
861  planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, sphereCenter, true, false);
862  geometrySagital->ChangeImageGeometryConsideringOriginOffset(true);
863 
864  origin = geometrySagital->GetOrigin();
865  normal = geometrySagital->GetNormal();
866 
867  normal.Normalize();
868 
869  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
870 
871  // geometrySagital->SetOrigin(origin);
872 
873  mitkExtractSliceFilterTestClass::TestSlice(geometrySagital, "Testing sagittal plane");
874  /* sagittal plane */
875 
876  /* sagittal shifted plane */
877  mitk::PlaneGeometry::Pointer geometrySagitalShifted = mitk::PlaneGeometry::New();
878  geometrySagitalShifted->InitializeStandardPlane(
879  planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, (sphereCenter - 14), true, false);
880  geometrySagitalShifted->ChangeImageGeometryConsideringOriginOffset(true);
881 
882  origin = geometrySagitalShifted->GetOrigin();
883  normal = geometrySagitalShifted->GetNormal();
884 
885  normal.Normalize();
886 
887  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
888 
889  // geometrySagitalShifted->SetOrigin(origin);
890 
891  mitkExtractSliceFilterTestClass::TestSlice(geometrySagitalShifted, "Testing sagittal plane shifted");
892  /* end sagittal shifted plane */
893 
894  /* coronal plane */
896  geometryCoronal->InitializeStandardPlane(
897  planeSize, planeSize, spacing, mitk::PlaneGeometry::Frontal, sphereCenter, true, false);
898  geometryCoronal->ChangeImageGeometryConsideringOriginOffset(true);
899 
900  origin = geometryCoronal->GetOrigin();
901  normal = geometryCoronal->GetNormal();
902 
903  normal.Normalize();
904 
905  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
906 
907  // geometryCoronal->SetOrigin(origin);
908 
909  mitkExtractSliceFilterTestClass::TestSlice(geometryCoronal, "Testing coronal plane");
910  /* end coronal plane */
911 
912  /* oblique plane */
914  obliquePlane->InitializeStandardPlane(
915  planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, sphereCenter, true, false);
916  obliquePlane->ChangeImageGeometryConsideringOriginOffset(true);
917 
918  origin = obliquePlane->GetOrigin();
919  normal = obliquePlane->GetNormal();
920 
921  normal.Normalize();
922 
923  origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5
924 
925  // obliquePlane->SetOrigin(origin);
926 
927  mitk::Vector3D rotationVector;
928  rotationVector[0] = 0.2;
929  rotationVector[1] = 0.4;
930  rotationVector[2] = 0.62;
931 
932  float degree = 37.0;
933 
935  new mitk::RotationOperation(mitk::OpROTATE, obliquePlane->GetCenter(), rotationVector, degree);
936  obliquePlane->ExecuteOperation(op);
937  delete op;
938 
939  mitkExtractSliceFilterTestClass::TestSlice(obliquePlane, "Testing oblique plane");
940 /* end oblique plane */
941 
942 #ifdef SHOW_SLICE_IN_RENDER_WINDOW
943  /*================ #BEGIN vtk render code ================*/
944 
945  // set reslicer for renderwindow
946 
947  mitk::Image::Pointer pic = mitk::IOUtil::Load<mitk::Image>(filename);
948  vtkSmartPointer<vtkImageReslice> slicer = vtkSmartPointer<vtkImageReslice>::New();
949 
950  slicer->SetInput(pic->GetVtkImageData());
951 
953  obliquePl->InitializeStandardPlane(
954  pic->GetGeometry(), mitk::PlaneGeometry::Sagittal, pic->GetGeometry()->GetCenter()[0], true, false);
955  obliquePl->ChangeImageGeometryConsideringOriginOffset(true);
956 
957  mitk::Point3D origin2 = obliquePl->GetOrigin();
958  mitk::Vector3D n;
959  n = obliquePl->GetNormal();
960 
961  n.Normalize();
962 
963  origin2 += n * 0.5; // pixelspacing is 1, so half the spacing is 0.5
964 
965  obliquePl->SetOrigin(origin2);
966 
968  rotation[0] = 0.534307;
969  rotation[1] = 0.000439605;
970  rotation[2] = 0.423017;
971  MITK_INFO << rotation;
972 
973  mitk::RotationOperation *operation =
974  new mitk::RotationOperation(mitk::OpROTATE, obliquePl->GetCenter(), rotationVector, degree);
975  obliquePl->ExecuteOperation(operation);
976  delete operation;
977 
978  double origin[3];
979  origin[0] = obliquePl->GetOrigin()[0];
980  origin[1] = obliquePl->GetOrigin()[1];
981  origin[2] = obliquePl->GetOrigin()[2];
982  slicer->SetResliceAxesOrigin(origin);
983 
984  mitk::Vector3D right, bottom, normal;
985  right = obliquePl->GetAxisVector(0);
986  bottom = obliquePl->GetAxisVector(1);
987  normal = obliquePl->GetNormal();
988 
989  right.Normalize();
990  bottom.Normalize();
991  normal.Normalize();
992 
993  double cosines[9];
994 
995  mitk::vnl2vtk(right.GetVnlVector(), cosines); // x
996 
997  mitk::vnl2vtk(bottom.GetVnlVector(), cosines + 3); // y
998 
999  mitk::vnl2vtk(normal.GetVnlVector(), cosines + 6); // n
1000 
1001  slicer->SetResliceAxesDirectionCosines(cosines);
1002 
1003  slicer->SetOutputDimensionality(2);
1004  slicer->Update();
1005 
1006  // set vtk renderwindow
1007  vtkSmartPointer<vtkPlaneSource> vtkPlane = vtkSmartPointer<vtkPlaneSource>::New();
1008  vtkPlane->SetOrigin(0.0, 0.0, 0.0);
1009 
1010  // These two points define the axes of the plane in combination with the origin.
1011  // Point 1 is the x-axis and point 2 the y-axis.
1012  // Each plane is transformed according to the view (axial, coronal and saggital) afterwards.
1013  vtkPlane->SetPoint1(1.0, 0.0, 0.0); // P1: (xMax, yMin, depth)
1014  vtkPlane->SetPoint2(0.0, 1.0, 0.0); // P2: (xMin, yMax, depth)
1015  // these are not the correct values for all slices, only a square plane by now
1016 
1017  vtkSmartPointer<vtkPolyDataMapper> imageMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
1018  imageMapper->SetInputConnection(vtkPlane->GetOutputPort());
1019 
1020  vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
1021 
1022  // built a default lookuptable
1023  lookupTable->SetRampToLinear();
1024  lookupTable->SetSaturationRange(0.0, 0.0);
1025  lookupTable->SetHueRange(0.0, 0.0);
1026  lookupTable->SetValueRange(0.0, 1.0);
1027  lookupTable->Build();
1028  // map all black values to transparent
1029  lookupTable->SetTableValue(0, 0.0, 0.0, 0.0, 0.0);
1030  lookupTable->SetRange(-255.0, 255.0);
1031  // lookupTable->SetRange(-1022.0, 1184.0);//pic3D range
1032 
1033  vtkSmartPointer<vtkTexture> texture = vtkSmartPointer<vtkTexture>::New();
1034 
1035  texture->SetInput(slicer->GetOutput());
1036 
1037  texture->SetLookupTable(lookupTable);
1038 
1039  texture->SetMapColorScalarsThroughLookupTable(true);
1040 
1041  vtkSmartPointer<vtkActor> imageActor = vtkSmartPointer<vtkActor>::New();
1042  imageActor->SetMapper(imageMapper);
1043  imageActor->SetTexture(texture);
1044 
1045  // Setup renderers
1046  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
1047  renderer->AddActor(imageActor);
1048 
1049  // Setup render window
1050  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
1051  renderWindow->AddRenderer(renderer);
1052 
1053  // Setup render window interactor
1054  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
1055  vtkSmartPointer<vtkInteractorStyleImage> style = vtkSmartPointer<vtkInteractorStyleImage>::New();
1056  renderWindowInteractor->SetInteractorStyle(style);
1057 
1058  // Render and start interaction
1059  renderWindowInteractor->SetRenderWindow(renderWindow);
1060  // renderer->AddViewProp(imageActor);
1061 
1062  renderWindow->Render();
1063  renderWindowInteractor->Start();
1064 // always end with this!
1065 /*================ #END vtk render code ================*/
1066 #endif // SHOW_SLICE_IN_RENDER_WINDOW
1067 
1068  MITK_TEST_END()
1069 }
void IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const
Convert (continuous or discrete) index coordinates of a vector vec_units to world coordinates (in mm)...
Gives locked and index-based read access for a particular image part. The class provides several set-...
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define AccessFixedDimensionByItk(mitkImage, itkImageTypeFunction, dimension)
Access a mitk-image with known dimension by an itk-image.
itk::Image< unsigned char, 3 > ImageType
double ScalarType
bool IsInside(const mitk::Point3D &p) const
Test whether the point p (world coordinates in mm) is inside the bounding box.
#define MITK_TEST_CONDITION_REQUIRED(COND, MSG)
static Pointer New()
virtual vtkImageData * GetVtkImageData(int t=0, int n=0)
Get a volume at a specific time t of channel n as a vtkImageData.
Definition: mitkImage.cpp:217
section GeneralTestsDeprecatedOldTestingStyle Deprecated macros All tests with MITK_TEST_BEGIN()
Constants for most interaction classes, due to the generic StateMachines.
static Matrix3D rotation
#define MITK_TEST_CONDITION(COND, MSG)
void vnl2vtk(const vnl_vector< Tin > &in, Tout *out)
Image class for storing images.
Definition: mitkImage.h:72
mitk::Image::Pointer image
static Pointer New()
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
static Pointer New()
Describes a two-dimensional, rectangular plane.
static Pointer New()
and MITK_TEST_END()
Operation, that holds everything necessary for an rotation operation on mitk::BaseData.
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:138
int mitkExtractSliceFilterTest(int, char *[])
static StandardFileLocations * GetInstance()