Medical Imaging Interaction Toolkit  2018.4.99-ef453c4b
Medical Imaging Interaction Toolkit
mitkGeometry3DTest.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 "mitkGeometry3D.h"
14 
15 #include <vnl/vnl_quaternion.h>
16 #include <vnl/vnl_quaternion.hxx>
17 
18 #include "mitkInteractionConst.h"
19 #include "mitkRotationOperation.h"
20 #include <mitkImageCast.h>
21 #include <mitkMatrixConvert.h>
22 
23 #include "mitkTestingMacros.h"
24 #include <fstream>
25 #include <mitkNumericTypes.h>
26 
28 {
29  int direction;
30  for (direction = 0; direction < 3; ++direction)
31  {
32  mitk::Vector3D frontToBack(0.);
33  switch (direction)
34  {
35  case 0:
36  frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(true, false, false);
37  break; // 7-3
38  case 1:
39  frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, true, false);
40  break; // 7-5
41  case 2:
42  frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, false, true);
43  break; // 7-2
44  }
45  std::cout << "Testing GetAxisVector(int) vs GetAxisVector(bool, bool, bool): ";
46  if (mitk::Equal(geometry->GetAxisVector(direction), frontToBack) == false)
47  {
48  std::cout << "[FAILED]" << std::endl;
49  return false;
50  }
51  std::cout << "[PASSED]" << std::endl;
52  }
53  return true;
54 }
55 
57 {
58  int direction;
59  for (direction = 0; direction < 3; ++direction)
60  {
61  if (mitk::Equal(geometry->GetAxisVector(direction).GetNorm(), geometry->GetExtentInMM(direction)) == false)
62  {
63  std::cout << "[FAILED]" << std::endl;
64  return false;
65  }
66  std::cout << "[PASSED]" << std::endl;
67  }
68  return true;
69 }
70 
71 // a part of the test requires axis-parallel coordinates
73 {
74  MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems: ");
75  mitk::Point3D origin = geometry3d->GetOrigin();
76  mitk::Point3D dummy;
77 
78  MITK_TEST_OUTPUT(<< " Testing index->world->index conversion consistency");
79  geometry3d->WorldToIndex(origin, dummy);
80  geometry3d->IndexToWorld(dummy, dummy);
81  MITK_TEST_CONDITION_REQUIRED(dummy == origin, "");
82 
83  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin, mitk::Point3D)==(0,0,0)");
84  mitk::Point3D globalOrigin;
85  mitk::FillVector3D(globalOrigin, 0, 0, 0);
86 
87  mitk::Point3D originContinuousIndex;
88  geometry3d->WorldToIndex(origin, originContinuousIndex);
89  MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, "");
90 
91  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin, itk::Index)==(0,0,0)");
92  itk::Index<3> itkindex;
93  geometry3d->WorldToIndex(origin, itkindex);
94  itk::Index<3> globalOriginIndex;
95  mitk::vtk2itk(globalOrigin, globalOriginIndex);
96  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
97 
98  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin-0.5*spacing, itk::Index)==(0,0,0)");
99  mitk::Vector3D halfSpacingStep = geometry3d->GetSpacing() * 0.5;
101  mitk::Point3D originOffCenter = origin - halfSpacingStep;
102  geometry3d->WorldToIndex(originOffCenter, itkindex);
103  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
104 
105  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)");
106  originOffCenter = origin + halfSpacingStep;
107  originOffCenter -= 0.0001;
108  geometry3d->WorldToIndex(originOffCenter, itkindex);
109  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
110 
111  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)");
112  originOffCenter = origin + halfSpacingStep;
113  itk::Index<3> global111;
114  mitk::FillVector3D(global111, 1, 1, 1);
115  geometry3d->WorldToIndex(originOffCenter, itkindex);
116  MITK_TEST_CONDITION_REQUIRED(itkindex == global111, "");
117 
118  MITK_TEST_OUTPUT(<< " Testing WorldToIndex(GetCenter())==BoundingBox.GetCenter: ");
119  mitk::Point3D center = geometry3d->GetCenter();
120  mitk::Point3D centerContIndex;
121  geometry3d->WorldToIndex(center, centerContIndex);
122  mitk::BoundingBox::ConstPointer boundingBox = geometry3d->GetBoundingBox();
123  mitk::BoundingBox::PointType centerBounds = boundingBox->GetCenter();
124  // itk assumes corner based geometry. If our geometry is center based (imageGoe == true), everything needs to be
125  // shifted
126  if (geometry3d->GetImageGeometry())
127  {
128  centerBounds[0] -= 0.5;
129  centerBounds[1] -= 0.5;
130  centerBounds[2] -= 0.5;
131  }
132  MITK_TEST_CONDITION_REQUIRED(mitk::Equal(centerContIndex, centerBounds), "");
133 
134  MITK_TEST_OUTPUT(<< " Testing GetCenter()==IndexToWorld(BoundingBox.GetCenter): ");
135  center = geometry3d->GetCenter();
136  mitk::Point3D centerBoundsInWorldCoords;
137  geometry3d->IndexToWorld(centerBounds, centerBoundsInWorldCoords);
138  MITK_TEST_CONDITION_REQUIRED(mitk::Equal(center, centerBoundsInWorldCoords), "");
139 
140  return EXIT_SUCCESS;
141 }
142 
144 {
145  MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems for vectors: ");
146  mitk::Vector3D xAxisMM = geometry3d->GetAxisVector(0);
147  mitk::Vector3D xAxisContinuousIndex;
148  mitk::Vector3D xAxisContinuousIndexDeprecated;
149 
150  mitk::Point3D p, pIndex, origin;
151  origin = geometry3d->GetOrigin();
152  p[0] = xAxisMM[0];
153  p[1] = xAxisMM[1];
154  p[2] = xAxisMM[2];
155 
156  geometry3d->WorldToIndex(p, pIndex);
157 
158  geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndexDeprecated);
159  geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndex);
160  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == pIndex[0], "");
161  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == pIndex[1], "");
162  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == pIndex[2], "");
163 
164  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == xAxisContinuousIndexDeprecated[0], "");
165  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == xAxisContinuousIndexDeprecated[1], "");
166  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == xAxisContinuousIndexDeprecated[2], "");
167 
168  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == pIndex[0], "");
169  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == pIndex[1], "");
170  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == pIndex[2], "");
171 
172  geometry3d->IndexToWorld(xAxisContinuousIndex, xAxisContinuousIndex);
173  geometry3d->IndexToWorld(xAxisContinuousIndexDeprecated, xAxisContinuousIndexDeprecated);
174  geometry3d->IndexToWorld(pIndex, p);
175 
176  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex == xAxisMM, "");
177  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == p[0], "");
178  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == p[1], "");
179  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == p[2], "");
180 
181  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated == xAxisMM, "");
182  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == p[0], "");
183  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == p[1], "");
184  MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == p[2], "");
185 
186  return EXIT_SUCCESS;
187 }
188 
190 {
191  MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems: ");
192 
193  // creating testing data
194  itk::Index<4> itkIndex4, itkIndex4b;
195  itk::Index<3> itkIndex3, itkIndex3b;
196  itk::Index<2> itkIndex2, itkIndex2b;
197  itk::Index<3> mitkIndex, mitkIndexb;
198 
199  itkIndex4[0] = itkIndex4[1] = itkIndex4[2] = itkIndex4[3] = 4;
200  itkIndex3[0] = itkIndex3[1] = itkIndex3[2] = 6;
201  itkIndex2[0] = itkIndex2[1] = 2;
202  mitkIndex[0] = mitkIndex[1] = mitkIndex[2] = 13;
203 
204  // check for constistency
205  mitk::Point3D point;
206  geometry3d->IndexToWorld(itkIndex2, point);
207  geometry3d->WorldToIndex(point, itkIndex2b);
208 
209  MITK_TEST_CONDITION_REQUIRED(((itkIndex2b[0] == itkIndex2[0]) && (itkIndex2b[1] == itkIndex2[1])),
210  "Testing itk::index<2> for IndexToWorld/WorldToIndex consistency");
211 
212  geometry3d->IndexToWorld(itkIndex3, point);
213  geometry3d->WorldToIndex(point, itkIndex3b);
214 
216  ((itkIndex3b[0] == itkIndex3[0]) && (itkIndex3b[1] == itkIndex3[1]) && (itkIndex3b[2] == itkIndex3[2])),
217  "Testing itk::index<3> for IndexToWorld/WorldToIndex consistency");
218 
219  geometry3d->IndexToWorld(itkIndex4, point);
220  geometry3d->WorldToIndex(point, itkIndex4b);
221 
222  MITK_TEST_CONDITION_REQUIRED(((itkIndex4b[0] == itkIndex4[0]) && (itkIndex4b[1] == itkIndex4[1]) &&
223  (itkIndex4b[2] == itkIndex4[2]) && (itkIndex4b[3] == 0)),
224  "Testing itk::index<3> for IndexToWorld/WorldToIndex consistency");
225 
226  geometry3d->IndexToWorld(mitkIndex, point);
227  geometry3d->WorldToIndex(point, mitkIndexb);
228 
230  ((mitkIndexb[0] == mitkIndex[0]) && (mitkIndexb[1] == mitkIndex[1]) && (mitkIndexb[2] == mitkIndex[2])),
231  "Testing mitk::Index for IndexToWorld/WorldToIndex consistency");
232 
233  return EXIT_SUCCESS;
234 }
235 
236 #include <itkImage.h>
237 
239 {
240  MITK_TEST_OUTPUT(<< "Testing whether itk::Image coordinates are center-based.");
241  typedef itk::Image<int, 3> ItkIntImage3D;
242  ItkIntImage3D::Pointer itkintimage = ItkIntImage3D::New();
243  ItkIntImage3D::SizeType size;
244  size.Fill(10);
245  mitk::Point3D origin;
246  mitk::FillVector3D(origin, 2, 3, 7);
247  itkintimage->Initialize();
248  itkintimage->SetRegions(size);
249  itkintimage->SetOrigin(origin);
250  std::cout << "[PASSED]" << std::endl;
251 
252  MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToContinuousIndex(origin)==(0,0,0)");
253  mitk::Point3D globalOrigin;
254  mitk::FillVector3D(globalOrigin, 0, 0, 0);
255 
256  itk::ContinuousIndex<mitk::ScalarType, 3> originContinuousIndex;
257  itkintimage->TransformPhysicalPointToContinuousIndex(origin, originContinuousIndex);
258  MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, "");
259 
260  MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin)==(0,0,0)");
261  itk::Index<3> itkindex;
262  itkintimage->TransformPhysicalPointToIndex(origin, itkindex);
263  itk::Index<3> globalOriginIndex;
264  mitk::vtk2itk(globalOrigin, globalOriginIndex);
265  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
266 
267  MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin-0.5*spacing)==(0,0,0)");
268  mitk::Vector3D halfSpacingStep = itkintimage->GetSpacing() * 0.5;
270  mitk::Point3D originOffCenter = origin - halfSpacingStep;
271  itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex);
272  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
273 
275  << " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)");
276  originOffCenter = origin + halfSpacingStep;
277  originOffCenter -= 0.0001;
278  itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex);
279  MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, "");
280 
281  MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)");
282  originOffCenter = origin + halfSpacingStep;
283  itk::Index<3> global111;
284  mitk::FillVector3D(global111, 1, 1, 1);
285  itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex);
286  MITK_TEST_CONDITION_REQUIRED(itkindex == global111, "");
287 
288  MITK_TEST_OUTPUT(<< "=> Yes, itk::Image coordinates are center-based.");
289 
290  return EXIT_SUCCESS;
291 }
292 
293 int testGeometry3D(bool imageGeometry)
294 {
295  // Build up a new image Geometry
297  float bounds[] = {-10.0, 17.0, -12.0, 188.0, 13.0, 211.0};
298 
299  MITK_TEST_OUTPUT(<< "Initializing");
300  geometry3d->Initialize();
301 
302  MITK_TEST_OUTPUT(<< "Setting ImageGeometry to " << imageGeometry);
303  geometry3d->SetImageGeometry(imageGeometry);
304 
305  MITK_TEST_OUTPUT(<< "Setting bounds by SetFloatBounds(): " << bounds);
306  geometry3d->SetFloatBounds(bounds);
307 
308  MITK_TEST_OUTPUT(<< "Testing AxisVectors");
309  if (testGetAxisVectorVariants(geometry3d) == false)
310  return EXIT_FAILURE;
311 
312  if (testGetAxisVectorExtent(geometry3d) == false)
313  return EXIT_FAILURE;
314 
315  MITK_TEST_OUTPUT(<< "Creating an AffineTransform3D transform");
316  mitk::AffineTransform3D::MatrixType matrix;
317  matrix.SetIdentity();
318  matrix(1, 1) = 2;
319  mitk::AffineTransform3D::Pointer transform;
320  transform = mitk::AffineTransform3D::New();
321  transform->SetMatrix(matrix);
322 
323  MITK_TEST_OUTPUT(<< "Testing a SetIndexToWorldTransform");
324  geometry3d->SetIndexToWorldTransform(transform);
325 
326  MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing");
327  const mitk::Vector3D &spacing1 = geometry3d->GetSpacing();
328  mitk::Vector3D expectedSpacing;
329  expectedSpacing.Fill(1.0);
330  expectedSpacing[1] = 2;
331  if (mitk::Equal(spacing1, expectedSpacing) == false)
332  {
333  MITK_TEST_OUTPUT(<< " [FAILED]");
334  return EXIT_FAILURE;
335  }
336 
337  MITK_TEST_OUTPUT(<< "Testing a Compose(transform)");
338  geometry3d->Compose(transform);
339 
340  MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing");
341  const mitk::Vector3D &spacing2 = geometry3d->GetSpacing();
342  expectedSpacing[1] = 4;
343  if (mitk::Equal(spacing2, expectedSpacing) == false)
344  {
345  MITK_TEST_OUTPUT(<< " [FAILED]");
346  return EXIT_FAILURE;
347  }
348 
349  MITK_TEST_OUTPUT(<< "Testing correctness of SetSpacing");
350  mitk::Vector3D newspacing;
351  mitk::FillVector3D(newspacing, 1.5, 2.5, 3.5);
352  geometry3d->SetSpacing(newspacing);
353  const mitk::Vector3D &spacing3 = geometry3d->GetSpacing();
354  if (mitk::Equal(spacing3, newspacing) == false)
355  {
356  MITK_TEST_OUTPUT(<< " [FAILED]");
357  return EXIT_FAILURE;
358  }
359 
360  // Seperate Test function for Index and World consistency
361  testIndexAndWorldConsistency(geometry3d);
364 
365  MITK_TEST_OUTPUT(<< "Testing a rotation of the geometry");
366  double angle = 35.0;
367  mitk::Vector3D rotationVector;
368  mitk::FillVector3D(rotationVector, 1, 0, 0);
369  mitk::Point3D center = geometry3d->GetCenter();
370  auto op = new mitk::RotationOperation(mitk::OpROTATE, center, rotationVector, angle);
371  geometry3d->ExecuteOperation(op);
372 
373  MITK_TEST_OUTPUT(<< "Testing mitk::GetRotation() and success of rotation");
375  mitk::GetRotation(geometry3d, rotation);
376  mitk::Vector3D voxelStep = rotation * newspacing;
377  mitk::Vector3D voxelStepIndex;
378  geometry3d->WorldToIndex(voxelStep, voxelStepIndex);
379  mitk::Vector3D expectedVoxelStepIndex;
380  expectedVoxelStepIndex.Fill(1);
381  MITK_TEST_CONDITION_REQUIRED(mitk::Equal(voxelStepIndex, expectedVoxelStepIndex), "");
382  delete op;
383  std::cout << "[PASSED]" << std::endl;
384 
385  MITK_TEST_OUTPUT(<< "Testing that ImageGeometry is still " << imageGeometry);
386  MITK_TEST_CONDITION_REQUIRED(geometry3d->GetImageGeometry() == imageGeometry, "");
387 
388  // Test if the translate function moves the origin correctly.
389  mitk::Point3D oldOrigin = geometry3d->GetOrigin();
390 
391  // use some random values for translation
392  mitk::Vector3D translationVector;
393  translationVector.SetElement(0, 17.5f);
394  translationVector.SetElement(1, -32.3f);
395  translationVector.SetElement(2, 4.0f);
396  // compute ground truth
397  mitk::Point3D tmpResult = geometry3d->GetOrigin() + translationVector;
398  geometry3d->Translate(translationVector);
399  MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), tmpResult), "Testing if origin was translated.");
400 
401  translationVector *= -1; // vice versa
402  geometry3d->Translate(translationVector);
403 
404  MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), oldOrigin),
405  "Testing if the translation could be done vice versa.");
406 
407  return EXIT_SUCCESS;
408 }
409 
411 {
412  // Epsilon. Allowed difference for rotationvalue
413  float eps = 0.0001;
414 
415  // Cast ITK and MITK images and see if geometry stays
416  typedef itk::Image<double, 2> Image2DType;
417  typedef itk::Image<double, 3> Image3DType;
418 
419  // Create 3D ITK Image from Scratch, cast to 3D MITK image, compare Geometries
420  Image3DType::Pointer image3DItk = Image3DType::New();
421  Image3DType::RegionType myRegion;
422  Image3DType::SizeType mySize;
423  Image3DType::IndexType myIndex;
424  Image3DType::SpacingType mySpacing;
425  Image3DType::DirectionType myDirection, rotMatrixX, rotMatrixY, rotMatrixZ;
426  mySpacing[0] = 31;
427  mySpacing[1] = 0.1;
428  mySpacing[2] = 2.9;
429  myIndex[0] = -15;
430  myIndex[1] = 15;
431  myIndex[2] = 12;
432  mySize[0] = 10;
433  mySize[1] = 2;
434  mySize[2] = 555;
435  myRegion.SetSize(mySize);
436  myRegion.SetIndex(myIndex);
437  image3DItk->SetSpacing(mySpacing);
438  image3DItk->SetRegions(myRegion);
439  image3DItk->Allocate();
440  image3DItk->FillBuffer(0);
441 
442  myDirection.SetIdentity();
443  rotMatrixX.SetIdentity();
444  rotMatrixY.SetIdentity();
445  rotMatrixZ.SetIdentity();
446 
447  mitk::Image::Pointer mitkImage;
448 
449  // direction [row] [coloum]
450  MITK_TEST_OUTPUT(<< "Casting a rotated 3D ITK Image to a MITK Image and check if Geometry is still same");
451  for (double rotX = 0; rotX < (itk::Math::pi * 2); rotX += 0.4)
452  {
453  // Set Rotation X
454  rotMatrixX[1][1] = cos(rotX);
455  rotMatrixX[1][2] = -sin(rotX);
456  rotMatrixX[2][1] = sin(rotX);
457  rotMatrixX[2][2] = cos(rotX);
458 
459  for (double rotY = 0; rotY < (itk::Math::pi * 2); rotY += 0.33)
460  {
461  // Set Rotation Y
462  rotMatrixY[0][0] = cos(rotY);
463  rotMatrixY[0][2] = sin(rotY);
464  rotMatrixY[2][0] = -sin(rotY);
465  rotMatrixY[2][2] = cos(rotY);
466 
467  for (double rotZ = 0; rotZ < (itk::Math::pi * 2); rotZ += 0.5)
468  {
469  // Set Rotation Z
470  rotMatrixZ[0][0] = cos(rotZ);
471  rotMatrixZ[0][1] = -sin(rotZ);
472  rotMatrixZ[1][0] = sin(rotZ);
473  rotMatrixZ[1][1] = cos(rotZ);
474 
475  // Multiply matrizes
476  myDirection = myDirection * rotMatrixX * rotMatrixY * rotMatrixZ;
477  image3DItk->SetDirection(myDirection);
478  mitk::CastToMitkImage(image3DItk, mitkImage);
479  const mitk::AffineTransform3D::MatrixType &matrix =
480  mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix();
481 
482  for (int row = 0; row < 3; row++)
483  {
484  for (int col = 0; col < 3; col++)
485  {
486  double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col];
487  double itkValue = myDirection[row][col];
488  double diff = mitkValue - itkValue;
489  // if you decrease this value, you can see that there might be QUITE high inaccuracy!!!
490  if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT!
491  {
492  std::cout << "Had a difference of : " << diff;
493  std::cout << "Error: Casting altered Geometry!";
494  std::cout << "ITK Matrix:\n" << myDirection;
495  std::cout << "Mitk Matrix (With Spacing):\n" << matrix;
496  std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing();
497  MITK_TEST_CONDITION_REQUIRED(false == true, "");
498  return false;
499  }
500  }
501  }
502  }
503  }
504  }
505 
506  // Create 2D ITK Image from Scratch, cast to 2D MITK image, compare Geometries
507  Image2DType::Pointer image2DItk = Image2DType::New();
508  Image2DType::RegionType myRegion2D;
509  Image2DType::SizeType mySize2D;
510  Image2DType::IndexType myIndex2D;
511  Image2DType::SpacingType mySpacing2D;
512  Image2DType::DirectionType myDirection2D, rotMatrix;
513  mySpacing2D[0] = 31;
514  mySpacing2D[1] = 0.1;
515  myIndex2D[0] = -15;
516  myIndex2D[1] = 15;
517  mySize2D[0] = 10;
518  mySize2D[1] = 2;
519  myRegion2D.SetSize(mySize2D);
520  myRegion2D.SetIndex(myIndex2D);
521  image2DItk->SetSpacing(mySpacing2D);
522  image2DItk->SetRegions(myRegion2D);
523  image2DItk->Allocate();
524  image2DItk->FillBuffer(0);
525 
526  myDirection2D.SetIdentity();
527  rotMatrix.SetIdentity();
528 
529  // direction [row] [coloum]
530  MITK_TEST_OUTPUT(<< "Casting a rotated 2D ITK Image to a MITK Image and check if Geometry is still same");
531  for (double rotTheta = 0; rotTheta < (itk::Math::pi * 2); rotTheta += 0.1)
532  {
533  // Set Rotation
534  rotMatrix[0][0] = cos(rotTheta);
535  rotMatrix[0][1] = -sin(rotTheta);
536  rotMatrix[1][0] = sin(rotTheta);
537  rotMatrix[1][1] = cos(rotTheta);
538 
539  // Multiply matrizes
540  myDirection2D = myDirection2D * rotMatrix;
541  image2DItk->SetDirection(myDirection2D);
542  mitk::CastToMitkImage(image2DItk, mitkImage);
543  const mitk::AffineTransform3D::MatrixType &matrix =
544  mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix();
545 
546  // Compare MITK and ITK matrix
547  for (int row = 0; row < 3; row++)
548  {
549  for (int col = 0; col < 3; col++)
550  {
551  double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col];
552  if ((row == 2) && (col == row))
553  {
554  if (mitkValue != 1)
555  {
556  MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 "
557  "MUST be 0 and value for 2|2 must be 1");
558  return false;
559  }
560  }
561  else if ((row == 2) || (col == 2))
562  {
563  if (mitkValue != 0)
564  {
565  MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 "
566  "MUST be 0 and value for 2|2 must be 1");
567  return false;
568  }
569  }
570  else
571  {
572  double itkValue = myDirection2D[row][col];
573  double diff = mitkValue - itkValue;
574  // if you decrease this value, you can see that there might be QUITE high inaccuracy!!!
575  if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT!
576  {
577  std::cout << "Had a difference of : " << diff;
578  std::cout << "Error: Casting altered Geometry!";
579  std::cout << "ITK Matrix:\n" << myDirection2D;
580  std::cout << "Mitk Matrix (With Spacing):\n" << matrix;
581  std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing();
582  MITK_TEST_CONDITION_REQUIRED(false == true, "");
583  return false;
584  }
585  }
586  }
587  }
588  }
589 
590  // THIS WAS TESTED:
591  // 2D ITK -> 2D MITK,
592  // 3D ITK -> 3D MITK,
593 
594  // Still need to test: 2D MITK Image with ADDITIONAL INFORMATION IN MATRIX -> 2D ITK
595  // 1. Possibility: 3x3 MITK matrix can be converted without loss into 2x2 ITK matrix
596  // 2. Possibility: 3x3 MITK matrix can only be converted with loss into 2x2 ITK matrix
597  // .. before implementing this, we wait for further development in geometry classes (e.g. Geoemtry3D::SetRotation(..))
598 
599  return EXIT_SUCCESS;
600 }
601 
602 int mitkGeometry3DTest(int /*argc*/, char * /*argv*/ [])
603 {
605 
606  int result;
607 
608  MITK_TEST_CONDITION_REQUIRED((result = testItkImageIsCenterBased()) == EXIT_SUCCESS, "");
609 
610  MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = false");
611  MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(false)) == EXIT_SUCCESS, "");
612 
613  MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = true");
614  MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(true)) == EXIT_SUCCESS, "");
615 
616  MITK_TEST_OUTPUT(<< "Running test to see if Casting MITK to ITK and the other way around destroys geometry");
617  MITK_TEST_CONDITION_REQUIRED((result = testGeometryAfterCasting()) == EXIT_SUCCESS, "");
618 
619  MITK_TEST_END();
620 
621  return EXIT_SUCCESS;
622 }
bool testGetAxisVectorVariants(mitk::Geometry3D *geometry)
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)...
virtual bool GetImageGeometry() const
Is this an ImageGeometry?
Standard implementation of BaseGeometry.
void GetRotation(const mitk::BaseGeometry *geometry, TMatrixType &itkmatrix)
Vector3D GetAxisVector(unsigned int direction) const
Get vector along bounding-box in the specified direction in mm.
int mitkGeometry3DTest(int, char *[])
#define MITK_TEST_CONDITION_REQUIRED(COND, MSG)
Point3D GetCornerPoint(int id) const
Get the position of the corner number id (in world coordinates)
section GeneralTestsDeprecatedOldTestingStyle Deprecated macros All tests with MITK_TEST_BEGIN()
static Pointer New()
Constants for most interaction classes, due to the generic StateMachines.
int testIndexAndWorldConsistencyForVectors(mitk::Geometry3D *geometry3d)
static Matrix3D rotation
Point3D GetCenter() const
Get the center of the bounding-box in mm.
void FillVector3D(Tout &out, mitk::ScalarType x, mitk::ScalarType y, mitk::ScalarType z)
Definition: mitkArray.h:106
ScalarType GetExtentInMM(int direction) const
Get the extent of the bounding-box in the specified direction in mm.
int testItkImageIsCenterBased()
int testGeometry3D(bool imageGeometry)
#define MITK_TEST_OUTPUT(x)
Output some text.
#define MITK_TEST_CONDITION(COND, MSG)
void vtk2itk(const Tin &in, Tout &out)
bool testGetAxisVectorExtent(mitk::Geometry3D *geometry)
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
int testIndexAndWorldConsistencyForIndex(mitk::Geometry3D *geometry3d)
MITKNEWMODULE_EXPORT bool Equal(mitk::ExampleDataStructure *leftHandSide, mitk::ExampleDataStructure *rightHandSide, mitk::ScalarType eps, bool verbose)
Returns true if the example data structures are considered equal.
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
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
MITKCORE_EXPORT const ScalarType eps
int testIndexAndWorldConsistency(mitk::Geometry3D *geometry3d)
and MITK_TEST_END()
Operation, that holds everything necessary for an rotation operation on mitk::BaseData.
int testGeometryAfterCasting()
virtual const BoundingBoxType * GetBoundingBox()