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