Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkPlaneGeometry.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 "mitkPlaneGeometry.h"
14 #include "mitkInteractionConst.h"
15 #include "mitkLine.h"
16 #include "mitkPlaneOperation.h"
17 
18 #include <itkSpatialOrientationAdapter.h>
19 
20 #include <vtkTransform.h>
21 
22 #include <vnl/vnl_cross.h>
23 
24 namespace mitk
25 {
26  PlaneGeometry::PlaneGeometry() : Superclass(), m_ReferenceGeometry(nullptr) { Initialize(); }
30  {
31  }
32 
34  {
35  bool rotation = true;
36 
37  auto matrix = transform->GetMatrix().GetVnlMatrix();
38  matrix.normalize_columns();
39 
40  auto det = vnl_determinant(matrix);
41  if (fabs(det-1.0) > epsilon)
42  {
43  MITK_WARN << "Invalid rotation matrix! Determinant != 1 (" << det << ")";
44  rotation = false;
45  }
46 
47  vnl_matrix_fixed<double, 3, 3> id; id.set_identity();
48  auto should_be_id = matrix*matrix.transpose();
49  should_be_id -= id;
50  auto max = should_be_id.absolute_value_max();
51  if (max > epsilon)
52  {
53  MITK_WARN << "Invalid rotation matrix! R*R^T != ID. Max value: " << max << " (should be 0)";
54  rotation = false;
55  }
56 
57  return rotation;
58  }
59 
61  {
62  CheckRotationMatrix(transform);
63  }
64 
66  {
67  // error: unused parameter 'bounds'
68  // this happens in release mode, where the assert macro is defined empty
69  // hence we "use" the parameter:
70  (void)bounds;
71 
72  // currently the unit rectangle must be starting at the origin [0,0]
73  assert(bounds[0] == 0);
74  assert(bounds[2] == 0);
75  // the unit rectangle must be two-dimensional
76  assert(bounds[1] > 0);
77  assert(bounds[3] > 0);
78  }
79 
80  void PlaneGeometry::IndexToWorld(const Point2D &pt_units, Point2D &pt_mm) const
81  {
82  pt_mm[0] = GetExtentInMM(0) / GetExtent(0) * pt_units[0];
83  pt_mm[1] = GetExtentInMM(1) / GetExtent(1) * pt_units[1];
84  }
85 
86  void PlaneGeometry::WorldToIndex(const Point2D &pt_mm, Point2D &pt_units) const
87  {
88  pt_units[0] = pt_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0)));
89  pt_units[1] = pt_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(1)));
90  }
91 
92  void PlaneGeometry::IndexToWorld(const Point2D & /*atPt2d_units*/, const Vector2D &vec_units, Vector2D &vec_mm) const
93  {
94  MITK_WARN << "Warning! Call of the deprecated function PlaneGeometry::IndexToWorld(point, vec, vec). Use "
95  "PlaneGeometry::IndexToWorld(vec, vec) instead!";
96  this->IndexToWorld(vec_units, vec_mm);
97  }
98 
99  void PlaneGeometry::IndexToWorld(const Vector2D &vec_units, Vector2D &vec_mm) const
100  {
101  vec_mm[0] = (GetExtentInMM(0) / GetExtent(0)) * vec_units[0];
102  vec_mm[1] = (GetExtentInMM(1) / GetExtent(1)) * vec_units[1];
103  }
104 
105  void PlaneGeometry::WorldToIndex(const Point2D & /*atPt2d_mm*/, const Vector2D &vec_mm, Vector2D &vec_units) const
106  {
107  MITK_WARN << "Warning! Call of the deprecated function PlaneGeometry::WorldToIndex(point, vec, vec). Use "
108  "PlaneGeometry::WorldToIndex(vec, vec) instead!";
109  this->WorldToIndex(vec_mm, vec_units);
110  }
111 
112  void PlaneGeometry::WorldToIndex(const Vector2D &vec_mm, Vector2D &vec_units) const
113  {
114  vec_units[0] = vec_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0)));
115  vec_units[1] = vec_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(1)));
116  }
117 
119  ScalarType height,
120  const Vector3D &spacing,
121  PlaneGeometry::PlaneOrientation planeorientation,
122  ScalarType zPosition,
123  bool frontside,
124  bool rotated,
125  bool top)
126  {
127  AffineTransform3D::Pointer transform;
128 
129  transform = AffineTransform3D::New();
130  AffineTransform3D::MatrixType matrix;
131  AffineTransform3D::MatrixType::InternalMatrixType &vnlmatrix = matrix.GetVnlMatrix();
132 
133  vnlmatrix.set_identity();
134  vnlmatrix(0, 0) = spacing[0];
135  vnlmatrix(1, 1) = spacing[1];
136  vnlmatrix(2, 2) = spacing[2];
137  transform->SetIdentity();
138  transform->SetMatrix(matrix);
139 
140  InitializeStandardPlane(width, height, transform.GetPointer(), planeorientation, zPosition, frontside, rotated, top);
141  }
142 
144  mitk::ScalarType height,
145  const AffineTransform3D *transform /* = nullptr */,
146  PlaneGeometry::PlaneOrientation planeorientation /* = Axial */,
147  mitk::ScalarType zPosition /* = 0 */,
148  bool frontside /* = true */,
149  bool rotated /* = false */,
150  bool top /* = true */)
151  {
153 
155 
156  // We define at the moment "frontside" as: axial from above,
157  // coronal from front (nose), saggital from right.
158  // TODO: Double check with medicals doctors or radiologists [ ].
159 
160  // We define the orientation in patient's view, e.g. LAI is in a axial cut
161  // (parallel to the triangle ear-ear-nose):
162  // first axis: To the left ear of the patient
163  // seecond axis: To the nose of the patient
164  // third axis: To the legs of the patient.
165 
166  // Options are: L/R left/right; A/P anterior/posterior; I/S inferior/superior
167  // (AKA caudal/cranial).
168  // We note on all cases in the following switch block r.h. for right handed
169  // or l.h. for left handed to describe the different cases.
170  // However, which system is chosen is defined at the end of the switch block.
171 
172  // CAVE / be careful: the vectors right and bottom are relative to the plane
173  // and do NOT describe e.g. the right side of the patient.
174 
175  Point3D origin;
177  VnlVector rightDV(3), bottomDV(3);
179  origin.Fill(0);
186  int normalDirection;
187 
188  switch (planeorientation) // Switch through our limited choice of standard planes:
189  {
190  case None:
193  case Axial:
194  if (frontside) // Radiologist's view from below. A cut along the triangle ear-ear-nose.
195  {
196  if (rotated == false)
208  {
209  FillVector3D(origin, 0, 0, zPosition);
210  FillVector3D(rightDV, 1, 0, 0);
211  FillVector3D(bottomDV, 0, 1, 0);
212  }
213  else // Origin rotated to the bottom-right corner, x=[-1; 0; 0], y=[0; -1; 0], z=[0; 0; 1],
214  // origin=[w,h,zpos]: RPI (r.h.)
215  { // Caveat emptor: Still using top-left as origin of index coordinate system!
216  FillVector3D(origin, width, height, zPosition);
217  FillVector3D(rightDV, -1, 0, 0);
218  FillVector3D(bottomDV, 0, -1, 0);
219  }
220  }
221  else // 'Backside, not frontside.' Neuro-Surgeons's view from above patient.
222  {
223  if (rotated == false) // x=[-1; 0; 0], y=[0; 1; 0], z=[0; 0; 1], origin=[w,0,zpos]: RAS (r.h.)
224  {
225  FillVector3D(origin, width, 0, zPosition);
226  FillVector3D(rightDV, -1, 0, 0);
227  FillVector3D(bottomDV, 0, 1, 0);
228  }
229  else // Origin in the bottom-left corner, x=[1; 0; 0], y=[0; -1; 0], z=[0; 0; 1],
230  // origin=[0,h,zpos]: LPS (r.h.)
231  {
232  FillVector3D(origin, 0, height, zPosition);
233  FillVector3D(rightDV, 1, 0, 0);
234  FillVector3D(bottomDV, 0, -1, 0);
235  }
236  }
237  normalDirection = 2; // That is S=Superior=z=third_axis=middlefinger in righthanded LPS-system.
238  break;
239 
240  // Frontal is known as Coronal in mitk. Plane cuts through patient's ear-ear-heel-heel:
241  case Frontal:
242  if (frontside)
243  {
244  if (rotated == false) // x=[1; 0; 0], y=[0; 0; 1], z=[0; 1; 0], origin=[0,zpos,0]: LAI (r.h.)
245  {
246  FillVector3D(origin, 0, zPosition, 0);
247  FillVector3D(rightDV, 1, 0, 0);
248  FillVector3D(bottomDV, 0, 0, 1);
249  }
250  else // x=[-1;0;0], y=[0;0;-1], z=[0;1;0], origin=[w,zpos,h]: RAS (r.h.)
251  {
252  FillVector3D(origin, width, zPosition, height);
253  FillVector3D(rightDV, -1, 0, 0);
254  FillVector3D(bottomDV, 0, 0, -1);
255  }
256  }
257  else
258  {
259  if (rotated == false) // x=[-1;0;0], y=[0;0;1], z=[0;1;0], origin=[w,zpos,0]: RPI (r.h.)
260  {
261  FillVector3D(origin, width, zPosition, 0);
262  FillVector3D(rightDV, -1, 0, 0);
263  FillVector3D(bottomDV, 0, 0, 1);
264  }
265  else // x=[1;0;0], y=[0;1;0], z=[0;0;-1], origin=[0,zpos,h]: LPS (r.h.)
266  {
267  FillVector3D(origin, 0, zPosition, height);
268  FillVector3D(rightDV, 1, 0, 0);
269  FillVector3D(bottomDV, 0, 0, -1);
270  }
271  }
272  normalDirection = 1; // Normal vector = posterior direction.
273  break;
274 
275  case Sagittal: // Sagittal=Medial plane, the symmetry-plane mirroring your face.
276  if (frontside)
277  {
278  if (rotated == false) // x=[0;1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,0,0]: LAI (r.h.)
279  {
280  FillVector3D(origin, zPosition, 0, 0);
281  FillVector3D(rightDV, 0, 1, 0);
282  FillVector3D(bottomDV, 0, 0, 1);
283  }
284  else // x=[0;-1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,w,h]: LPS (r.h.)
285  {
286  FillVector3D(origin, zPosition, width, height);
287  FillVector3D(rightDV, 0, -1, 0);
288  FillVector3D(bottomDV, 0, 0, -1);
289  }
290  }
291  else
292  {
293  if (rotated == false) // x=[0;-1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,w,0]: RPI (r.h.)
294  {
295  FillVector3D(origin, zPosition, width, 0);
296  FillVector3D(rightDV, 0, -1, 0);
297  FillVector3D(bottomDV, 0, 0, 1);
298  }
299  else // x=[0;1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,0,h]: RAS (r.h.)
300  {
301  FillVector3D(origin, zPosition, 0, height);
302  FillVector3D(rightDV, 0, 1, 0);
303  FillVector3D(bottomDV, 0, 0, -1);
304  }
305  }
306  normalDirection = 0; // Normal vector = Lateral direction: Left in a LPS-system.
307  break;
308 
309  default:
310  itkExceptionMacro("unknown PlaneOrientation");
311  }
312 
313  VnlVector normal(3);
314  FillVector3D(normal, 0, 0, 0);
315  normal[normalDirection] = top ? 1 : -1;
316 
317  if ( transform != nullptr )
318  {
319  origin = transform->TransformPoint( origin );
320  rightDV = transform->TransformVector( rightDV );
321  bottomDV = transform->TransformVector( bottomDV );
322  normal = transform->TransformVector( normal );
323  }
324 
325  ScalarType bounds[6] = {0, width, 0, height, 0, 1};
326  this->SetBounds(bounds);
327 
328  AffineTransform3D::Pointer planeTransform = AffineTransform3D::New();
329  Matrix3D matrix;
330  matrix.GetVnlMatrix().set_column(0, rightDV);
331  matrix.GetVnlMatrix().set_column(1, bottomDV);
332  matrix.GetVnlMatrix().set_column(2, normal);
333  planeTransform->SetMatrix(matrix);
334  planeTransform->SetOffset(this->GetIndexToWorldTransform()->GetOffset());
335  this->SetIndexToWorldTransform(planeTransform);
336 
337  this->SetOrigin(origin);
338  }
339 
340  std::vector< int > PlaneGeometry::CalculateDominantAxes(mitk::AffineTransform3D::MatrixType::InternalMatrixType& rotation_matrix)
341  {
342  std::vector< int > axes;
343 
344  bool dominant_axis_error = false;
345  for (int i = 0; i < 3; ++i)
346  {
347  int dominantAxis = itk::Function::Max3(
348  rotation_matrix[0][i],
349  rotation_matrix[1][i],
350  rotation_matrix[2][i]
351  );
352 
353  for (int j=0; j<i; ++j)
354  if (axes[j] == dominantAxis)
355  {
356  dominant_axis_error = true;
357  break;
358  }
359  if (dominant_axis_error)
360  break;
361 
362  axes.push_back(dominantAxis);
363  }
364 
365  if (dominant_axis_error)
366  {
367  MITK_DEBUG << "Error during dominant axis calculation. Using default.";
368  MITK_DEBUG << "This is either caused by an imperfect rotation matrix or if the rotation is axactly 45° around one or more axis.";
369  axes.clear();
370  for (int i = 0; i < 3; ++i)
371  axes.push_back(i);
372  }
373 
374  return axes;
375  }
376 
378  PlaneOrientation planeorientation,
379  ScalarType zPosition,
380  bool frontside,
381  bool rotated,
382  bool top)
383  {
384  this->SetReferenceGeometry(geometry3D);
385 
386  ScalarType width, height;
387 
388  // Inspired by:
389  // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
390 
391  mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
392 
393  matrix.GetVnlMatrix().normalize_columns();
394  mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetTranspose();
395 
397  auto axes = CalculateDominantAxes(inverseMatrix);
401  int directions[3];
402  ScalarType extents[3];
403  ScalarType spacings[3];
404  for (int i=0; i<3; ++i)
405  {
406  int dominantAxis = axes.at(i);
407  directions[i] = itk::Function::Sign(inverseMatrix[dominantAxis][i]);
408  extents[i] = geometry3D->GetExtent(dominantAxis);
409  spacings[i] = geometry3D->GetSpacing()[dominantAxis];
410  }
411 
412  // matrix(column) = inverseTransformMatrix(row) * flippedAxes * spacing
413  matrix[0][0] = inverseMatrix[axes[0]][0] * directions[0] * spacings[0];
414  matrix[1][0] = inverseMatrix[axes[0]][1] * directions[0] * spacings[0];
415  matrix[2][0] = inverseMatrix[axes[0]][2] * directions[0] * spacings[0];
416  matrix[0][1] = inverseMatrix[axes[1]][0] * directions[1] * spacings[1];
417  matrix[1][1] = inverseMatrix[axes[1]][1] * directions[1] * spacings[1];
418  matrix[2][1] = inverseMatrix[axes[1]][2] * directions[1] * spacings[1];
419  matrix[0][2] = inverseMatrix[axes[2]][0] * directions[2] * spacings[2];
420  matrix[1][2] = inverseMatrix[axes[2]][1] * directions[2] * spacings[2];
421  matrix[2][2] = inverseMatrix[axes[2]][2] * directions[2] * spacings[2];
422 
426  Point3D worldOrigin = geometry3D->GetOrigin();
427  for (int i = 0; i < 3; ++i)
428  {
430  double offset = directions[i] > 0 ? 0.0 : extents[i];
431 
432  if (geometry3D->GetImageGeometry())
433  {
434  offset += directions[i] * 0.5;
435  }
436 
437  for (int j = 0; j < 3; ++j)
438  {
439  worldOrigin[j] -= offset * matrix[j][i];
440  }
441  }
442 
443  switch(planeorientation)
444  {
445  case None:
448  case Axial:
449  width = extents[0];
450  height = extents[1];
451  break;
452  case Frontal:
453  width = extents[0];
454  height = extents[2];
455  break;
456  case Sagittal:
457  width = extents[1];
458  height = extents[2];
459  break;
460  default:
461  itkExceptionMacro("unknown PlaneOrientation");
462  }
463 
464  ScalarType bounds[6]= { 0, width, 0, height, 0, 1 };
465  this->SetBounds( bounds );
466 
467  AffineTransform3D::Pointer transform = AffineTransform3D::New();
468  transform->SetMatrix(matrix);
469  transform->SetOffset(worldOrigin.GetVectorFromOrigin());
470 
472  width, height, transform, planeorientation, zPosition, frontside, rotated, top);
473  }
474 
476  const BaseGeometry *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated)
477  {
479  int worldAxis;
480  switch(planeorientation)
481  {
482  case None:
485  case Axial:
486  worldAxis = 2;
487  break;
488  case Frontal:
489  worldAxis = 1;
490  break;
491  case Sagittal:
492  worldAxis = 0;
493  break;
494  default:
495  itkExceptionMacro("unknown PlaneOrientation");
496  }
497 
498  // Inspired by:
499  // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
500 
501  mitk::AffineTransform3D::ConstPointer affineTransform = geometry3D->GetIndexToWorldTransform();
502  mitk::AffineTransform3D::MatrixType matrix = affineTransform->GetMatrix();
503  matrix.GetVnlMatrix().normalize_columns();
504  mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
505 
507  int dominantAxis = CalculateDominantAxes(inverseMatrix).at(worldAxis);
508 
509  ScalarType zPosition = top ? 0.5 : geometry3D->GetExtent(dominantAxis) - 0.5;
510 
511  InitializeStandardPlane(geometry3D, planeorientation, zPosition, frontside, rotated, top);
512  }
513 
515  const Vector3D &downVector,
516  const Vector3D *spacing)
517  {
518  InitializeStandardPlane(rightVector.GetVnlVector(), downVector.GetVnlVector(), spacing);
519  }
520 
522  const VnlVector &downVector,
523  const Vector3D *spacing)
524  {
525  ScalarType width = rightVector.two_norm();
526  ScalarType height = downVector.two_norm();
527 
528  InitializeStandardPlane(width, height, rightVector, downVector, spacing);
529  }
530 
532  ScalarType height,
533  const Vector3D &rightVector,
534  const Vector3D &downVector,
535  const Vector3D *spacing)
536  {
537  InitializeStandardPlane(width, height, rightVector.GetVnlVector(), downVector.GetVnlVector(), spacing);
538  }
539 
541  ScalarType height,
542  const VnlVector &rightVector,
543  const VnlVector &downVector,
544  const Vector3D *spacing)
545  {
546  assert(width > 0);
547  assert(height > 0);
548 
549  VnlVector rightDV = rightVector;
550  rightDV.normalize();
551  VnlVector downDV = downVector;
552  downDV.normalize();
553  VnlVector normal = vnl_cross_3d(rightVector, downVector);
554  normal.normalize();
555  // Crossproduct vnl_cross_3d is always righthanded, but that is okay here
556  // because in this method we create a new IndexToWorldTransform and
557  // spacing with 1 or 3 negative components could still make it lefthanded.
558 
559  if (spacing != nullptr)
560  {
561  rightDV *= (*spacing)[0];
562  downDV *= (*spacing)[1];
563  normal *= (*spacing)[2];
564  }
565 
566  AffineTransform3D::Pointer transform = AffineTransform3D::New();
567  Matrix3D matrix;
568  matrix.GetVnlMatrix().set_column(0, rightDV);
569  matrix.GetVnlMatrix().set_column(1, downDV);
570  matrix.GetVnlMatrix().set_column(2, normal);
571  transform->SetMatrix(matrix);
572  transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset());
573 
574  ScalarType bounds[6] = {0, width, 0, height, 0, 1};
575  this->SetBounds(bounds);
576 
577  this->SetIndexToWorldTransform(transform);
578  }
579 
580  void PlaneGeometry::InitializePlane(const Point3D &origin, const Vector3D &normal)
581  {
582  VnlVector rightVectorVnl(3), downVectorVnl;
583 
584  if (Equal(normal[1], 0.0f) == false)
585  {
586  FillVector3D(rightVectorVnl, 1.0f, -normal[0] / normal[1], 0.0f);
587  rightVectorVnl.normalize();
588  }
589  else
590  {
591  FillVector3D(rightVectorVnl, 0.0f, 1.0f, 0.0f);
592  }
593  downVectorVnl = vnl_cross_3d(normal.GetVnlVector(), rightVectorVnl);
594  downVectorVnl.normalize();
595  // Crossproduct vnl_cross_3d is always righthanded.
596 
597  InitializeStandardPlane(rightVectorVnl, downVectorVnl);
598 
599  SetOrigin(origin);
600  }
601 
603  const VnlVector &downVector,
604  ScalarType thickness /* = 1.0 */)
605  {
606  VnlVector normal = vnl_cross_3d(rightVector, downVector);
607  normal.normalize();
608  normal *= thickness;
609  // Crossproduct vnl_cross_3d is always righthanded, but that is okay here
610  // because in this method we create a new IndexToWorldTransform and
611  // a negative thickness could still make it lefthanded.
612 
613  AffineTransform3D::Pointer transform = AffineTransform3D::New();
614  Matrix3D matrix;
615  matrix.GetVnlMatrix().set_column(0, rightVector);
616  matrix.GetVnlMatrix().set_column(1, downVector);
617  matrix.GetVnlMatrix().set_column(2, normal);
618  transform->SetMatrix(matrix);
619  transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset());
620  SetIndexToWorldTransform(transform);
621  }
622 
624  {
625  Vector3D frontToBack;
626  frontToBack.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2));
627 
628  return frontToBack;
629  }
630 
632  {
633  return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2);
634  }
635 
636  ScalarType PlaneGeometry::DistanceFromPlane(const Point3D &pt3d_mm) const { return fabs(SignedDistance(pt3d_mm)); }
637  ScalarType PlaneGeometry::SignedDistance(const Point3D &pt3d_mm) const { return SignedDistanceFromPlane(pt3d_mm); }
638  bool PlaneGeometry::IsAbove(const Point3D &pt3d_mm, bool considerBoundingBox) const
639  {
640  if (considerBoundingBox)
641  {
642  Point3D pt3d_units;
643  BaseGeometry::WorldToIndex(pt3d_mm, pt3d_units);
644  return (pt3d_units[2] > this->GetBoundingBox()->GetBounds()[4]);
645  }
646  else
647  return SignedDistanceFromPlane(pt3d_mm) > 0;
648  }
649 
650  bool PlaneGeometry::IntersectionLine(const PlaneGeometry *plane, Line3D &crossline) const
651  {
652  Vector3D normal = this->GetNormal();
653  normal.Normalize();
654 
655  Vector3D planeNormal = plane->GetNormal();
656  planeNormal.Normalize();
657 
658  Vector3D direction = itk::CrossProduct(normal, planeNormal);
659 
660  if (direction.GetSquaredNorm() < eps)
661  return false;
662 
663  crossline.SetDirection(direction);
664 
665  double N1dN2 = normal * planeNormal;
666  double determinant = 1.0 - N1dN2 * N1dN2;
667 
668  Vector3D origin = this->GetOrigin().GetVectorFromOrigin();
669  Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin();
670 
671  double d1 = normal * origin;
672  double d2 = planeNormal * planeOrigin;
673 
674  double c1 = (d1 - d2 * N1dN2) / determinant;
675  double c2 = (d2 - d1 * N1dN2) / determinant;
676 
677  Vector3D p = normal * c1 + planeNormal * c2;
678  crossline.GetPoint().GetVnlVector() = p.GetVnlVector();
679 
680  return true;
681  }
682 
683  unsigned int PlaneGeometry::IntersectWithPlane2D(const PlaneGeometry *plane, Point2D &lineFrom, Point2D &lineTo) const
684  {
685  Line3D crossline;
686  if (this->IntersectionLine(plane, crossline) == false)
687  return 0;
688 
689  Point2D point2;
690  Vector2D direction2;
691 
692  this->Map(crossline.GetPoint(), point2);
693  this->Map(crossline.GetPoint(), crossline.GetDirection(), direction2);
694 
696  0, 0, GetExtentInMM(0), GetExtentInMM(1), point2, direction2, lineFrom, lineTo);
697  }
698 
699  double PlaneGeometry::Angle(const PlaneGeometry *plane) const
700  {
701  return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2));
702  }
703 
704  double PlaneGeometry::Angle(const Line3D &line) const
705  {
706  return vnl_math::pi_over_2 - angle(line.GetDirection().GetVnlVector(), GetMatrixColumn(2));
707  }
708 
709  bool PlaneGeometry::IntersectionPoint(const Line3D &line, Point3D &intersectionPoint) const
710  {
711  Vector3D planeNormal = this->GetNormal();
712  planeNormal.Normalize();
713 
714  Vector3D lineDirection = line.GetDirection();
715  lineDirection.Normalize();
716 
717  double t = planeNormal * lineDirection;
718  if (fabs(t) < eps)
719  {
720  return false;
721  }
722 
723  Vector3D diff;
724  diff = this->GetOrigin() - line.GetPoint();
725  t = (planeNormal * diff) / t;
726 
727  intersectionPoint = line.GetPoint() + lineDirection * t;
728  return true;
729  }
730 
731  bool PlaneGeometry::IntersectionPointParam(const Line3D &line, double &t) const
732  {
733  Vector3D planeNormal = this->GetNormal();
734 
735  Vector3D lineDirection = line.GetDirection();
736 
737  t = planeNormal * lineDirection;
738 
739  if (fabs(t) < eps)
740  {
741  return false;
742  }
743 
744  Vector3D diff;
745  diff = this->GetOrigin() - line.GetPoint();
746  t = (planeNormal * diff) / t;
747  return true;
748  }
749 
750  bool PlaneGeometry::IsParallel(const PlaneGeometry *plane) const
751  {
752  return ((Angle(plane) < 10.0 * mitk::sqrteps) || (Angle(plane) > (vnl_math::pi - 10.0 * sqrteps)));
753  }
754 
755  bool PlaneGeometry::IsOnPlane(const Point3D &point) const { return Distance(point) < eps; }
757  {
758  return ((Distance(line.GetPoint()) < eps) && (Distance(line.GetPoint2()) < eps));
759  }
760 
761  bool PlaneGeometry::IsOnPlane(const PlaneGeometry *plane) const
762  {
763  return (IsParallel(plane) && (Distance(plane->GetOrigin()) < eps));
764  }
765 
767  {
768  ScalarType len = this->GetNormalVnl().two_norm();
769  return pt - this->GetNormal() * this->SignedDistanceFromPlane(pt) / len;
770  }
771 
772  itk::LightObject::Pointer PlaneGeometry::InternalClone() const
773  {
774  Self::Pointer newGeometry = new PlaneGeometry(*this);
775  newGeometry->UnRegister();
776  return newGeometry.GetPointer();
777  }
778 
780  {
781  vtkTransform *transform = vtkTransform::New();
782  transform->SetMatrix(this->GetVtkMatrix());
783 
784  switch (operation->GetOperationType())
785  {
786  case OpORIENT:
787  {
788  auto *planeOp = dynamic_cast<mitk::PlaneOperation *>(operation);
789  if (planeOp == nullptr)
790  {
791  return;
792  }
793 
794  Point3D center = planeOp->GetPoint();
795 
796  Vector3D orientationVector = planeOp->GetNormal();
797  Vector3D defaultVector;
798  FillVector3D(defaultVector, 0.0, 0.0, 1.0);
799 
800  Vector3D rotationAxis = itk::CrossProduct(orientationVector, defaultVector);
801  // double rotationAngle = acos( orientationVector[2] / orientationVector.GetNorm() );
802 
803  double rotationAngle = atan2((double)rotationAxis.GetNorm(), (double)(orientationVector * defaultVector));
804  rotationAngle *= 180.0 / vnl_math::pi;
805 
806  transform->PostMultiply();
807  transform->Identity();
808  transform->Translate(center[0], center[1], center[2]);
809  transform->RotateWXYZ(rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2]);
810  transform->Translate(-center[0], -center[1], -center[2]);
811  break;
812  }
814  {
815  auto *op = dynamic_cast<mitk::RestorePlanePositionOperation *>(operation);
816  if (op == nullptr)
817  {
818  return;
819  }
820 
821  AffineTransform3D::Pointer transform2 = AffineTransform3D::New();
822  Matrix3D matrix;
823  matrix.GetVnlMatrix().set_column(0, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(0));
824  matrix.GetVnlMatrix().set_column(1, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(1));
825  matrix.GetVnlMatrix().set_column(2, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(2));
826  transform2->SetMatrix(matrix);
827  Vector3D offset = op->GetTransform()->GetOffset();
828  transform2->SetOffset(offset);
829 
830  this->SetIndexToWorldTransform(transform2);
831  ScalarType bounds[6] = {0, op->GetWidth(), 0, op->GetHeight(), 0, 1};
832  this->SetBounds(bounds);
833  this->Modified();
834  transform->Delete();
835  return;
836  }
837  default:
838  Superclass::ExecuteOperation(operation);
839  transform->Delete();
840  return;
841  }
842 
843  this->SetVtkMatrixDeepCopy(transform);
844  this->Modified();
845  transform->Delete();
846  }
847 
848  void PlaneGeometry::PrintSelf(std::ostream &os, itk::Indent indent) const
849  {
850  Superclass::PrintSelf(os, indent);
851  os << indent << " ScaleFactorMMPerUnitX: " << GetExtentInMM(0) / GetExtent(0) << std::endl;
852  os << indent << " ScaleFactorMMPerUnitY: " << GetExtentInMM(1) / GetExtent(1) << std::endl;
853  os << indent << " Normal: " << GetNormal() << std::endl;
854  }
855 
856  bool PlaneGeometry::Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const
857  {
858  assert(this->IsBoundingBoxNull() == false);
859 
860  Point3D pt3d_units;
861  Superclass::WorldToIndex(pt3d_mm, pt3d_units);
862  pt2d_mm[0] = pt3d_units[0] * GetExtentInMM(0) / GetExtent(0);
863  pt2d_mm[1] = pt3d_units[1] * GetExtentInMM(1) / GetExtent(1);
864  pt3d_units[2] = 0;
865  return this->GetBoundingBox()->IsInside(pt3d_units);
866  }
867 
868  void PlaneGeometry::Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const
869  {
870  // pt2d_mm is measured from the origin of the world geometry (at leats it called form BaseRendere::Mouse...Event)
871  Point3D pt3d_units;
872  pt3d_units[0] = pt2d_mm[0] / (GetExtentInMM(0) / GetExtent(0));
873  pt3d_units[1] = pt2d_mm[1] / (GetExtentInMM(1) / GetExtent(1));
874  pt3d_units[2] = 0;
875  // pt3d_units is a continuos index. We divided it with the Scale Factor (= spacing in x and y) to convert it from mm
876  // to index units.
877  //
878  pt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units);
879  // now we convert the 3d index to a 3D world point in mm. We could have used IndexToWorld as well as
880  // GetITW->Transform...
881  }
882 
884  {
885  ScalarType bounds[6] = {0, width, 0, height, 0, 1};
886  ScalarType extent, newextentInMM;
887  if (GetExtent(0) > 0)
888  {
889  extent = GetExtent(0);
890  if (width > extent)
891  newextentInMM = GetExtentInMM(0) / width * extent;
892  else
893  newextentInMM = GetExtentInMM(0) * extent / width;
894  SetExtentInMM(0, newextentInMM);
895  }
896  if (GetExtent(1) > 0)
897  {
898  extent = GetExtent(1);
899  if (width > extent)
900  newextentInMM = GetExtentInMM(1) / height * extent;
901  else
902  newextentInMM = GetExtentInMM(1) * extent / height;
903  SetExtentInMM(1, newextentInMM);
904  }
905  SetBounds(bounds);
906  }
907 
908  bool PlaneGeometry::Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const
909  {
910  assert(this->IsBoundingBoxNull() == false);
911 
912  Point3D pt3d_units;
913  Superclass::WorldToIndex(pt3d_mm, pt3d_units);
914  pt3d_units[2] = 0;
915  projectedPt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units);
916  return this->GetBoundingBox()->IsInside(pt3d_units);
917  }
918 
919  bool PlaneGeometry::Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const
920  {
921  assert(this->IsBoundingBoxNull() == false);
922 
923  Vector3D vec3d_units;
924  Superclass::WorldToIndex(vec3d_mm, vec3d_units);
925  vec3d_units[2] = 0;
926  projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units);
927  return true;
928  }
929 
930  bool PlaneGeometry::Project(const mitk::Point3D &atPt3d_mm,
931  const mitk::Vector3D &vec3d_mm,
932  mitk::Vector3D &projectedVec3d_mm) const
933  {
934  MITK_WARN << "Deprecated function! Call Project(vec3D,vec3D) instead.";
935  assert(this->IsBoundingBoxNull() == false);
936 
937  Vector3D vec3d_units;
938  Superclass::WorldToIndex(atPt3d_mm, vec3d_mm, vec3d_units);
939  vec3d_units[2] = 0;
940  projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units);
941 
942  Point3D pt3d_units;
943  Superclass::WorldToIndex(atPt3d_mm, pt3d_units);
944  return this->GetBoundingBox()->IsInside(pt3d_units);
945  }
946 
947  bool PlaneGeometry::Map(const mitk::Point3D &atPt3d_mm,
948  const mitk::Vector3D &vec3d_mm,
949  mitk::Vector2D &vec2d_mm) const
950  {
951  Point2D pt2d_mm_start, pt2d_mm_end;
952  Point3D pt3d_mm_end;
953  bool inside = Map(atPt3d_mm, pt2d_mm_start);
954  pt3d_mm_end = atPt3d_mm + vec3d_mm;
955  inside &= Map(pt3d_mm_end, pt2d_mm_end);
956  vec2d_mm = pt2d_mm_end - pt2d_mm_start;
957  return inside;
958  }
959 
960  void PlaneGeometry::Map(const mitk::Point2D & /*atPt2d_mm*/,
961  const mitk::Vector2D & /*vec2d_mm*/,
962  mitk::Vector3D & /*vec3d_mm*/) const
963  {
964  //@todo implement parallel to the other Map method!
965  assert(false);
966  }
967 
970  bool PlaneGeometry::HasReferenceGeometry() const { return (m_ReferenceGeometry != nullptr); }
971 } // namespace
VnlVector GetNormalVnl() const
Normal of the plane as VnlVector.
virtual ScalarType SignedDistance(const Point3D &pt3d_mm) const
Descibes a line.
Definition: mitkLine.h:28
static char * line
Definition: svm.cpp:2870
virtual bool GetImageGeometry() const
Is this an ImageGeometry?
static std::vector< int > CalculateDominantAxes(mitk::AffineTransform3D::MatrixType::InternalMatrixType &rotation_matrix)
vtkMatrix4x4 * GetVtkMatrix()
void ExecuteOperation(Operation *operation) override
void SetDirection(const itk::Vector< TCoordRep, NPointDimension > &direction)
Set the direction vector of the line.
Definition: mitkLine.h:76
void SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
bool IsBoundingBoxNull() const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Base class of all Operation-classes.
Definition: mitkOperation.h:29
void PrintSelf(std::ostream &os, itk::Indent indent) const override
vnl_vector< ScalarType > VnlVector
Definition: mitkVector.h:134
double ScalarType
static bool CheckRotationMatrix(AffineTransform3D *transform, double epsilon=mitk::eps)
Check if matrix is a rotation matrix:
bool IntersectionPointParam(const Line3D &line, double &t) const
Calculate line parameter of intersection point between the plane and a line.
unsigned int IntersectWithPlane2D(const PlaneGeometry *plane, Point2D &lineFrom, Point2D &lineTo) const
Calculate two points where another plane intersects the border of this plane.
bool IsOnPlane(const Point3D &point) const
Returns whether the point is on the plane (bounding-box not considered)
#define MITK_DEBUG
Definition: mitkLogMacros.h:22
bool IntersectionPoint(const Line3D &line, Point3D &intersectionPoint) const
Calculate intersection point between the plane and a line.
ScalarType GetExtent(unsigned int direction) const
Set the time bounds (in ms)
DataCollection - Class to facilitate loading/accessing structured data.
void Initialize()
Initialize the BaseGeometry.
Constants for most interaction classes, due to the generic StateMachines.
static Matrix3D rotation
virtual bool Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const
Project a 3D point given in mm (pt3d_mm) onto the 2D geometry. The result is a 2D point in mm (pt2d_m...
virtual bool Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const
Project a 3D point given in mm (pt3d_mm) onto the 2D geometry. The result is a 3D point in mm (projec...
void CheckIndexToWorldTransform(mitk::AffineTransform3D *transform) override
CheckIndexToWorldTransform.
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.
const BaseGeometry * GetReferenceGeometry() const
Get the geometrical frame of reference for this PlaneGeometry.
virtual void IndexToWorld(const Point2D &pt_units, Point2D &pt_mm) const
const itk::Point< TCoordRep, NPointDimension > & GetPoint() const
Get start point of the line.
Definition: mitkLine.h:49
const mitk::BaseGeometry * m_ReferenceGeometry
bool HasReferenceGeometry() const
static Vector3D offset
bool IntersectionLine(const PlaneGeometry *plane, Line3D &crossline) const
Calculate the intersecting line of two planes.
const itk::Vector< TCoordRep, NPointDimension > & GetDirection() const
Get the direction vector of the line.
Definition: mitkLine.h:70
void SetVtkMatrixDeepCopy(vtkTransform *vtktransform)
ScalarType DistanceFromPlane(const Point3D &pt3d_mm) const
Distance of the point from the plane (bounding-box not considered)
#define MITK_WARN
Definition: mitkLogMacros.h:19
virtual void SetSizeInUnits(mitk::ScalarType width, mitk::ScalarType height)
Set the width and height of this 2D-geometry in units by calling SetBounds. This does not change the ...
virtual void WorldToIndex(const Point2D &pt_mm, Point2D &pt_units) const
double Angle(const PlaneGeometry *plane) const
Calculate the angle between two planes.
virtual bool IsAbove(const Point3D &pt3d_mm, bool considerBoundingBox=false) const
Calculates, whether a point is below or above the plane. There are two different calculation methods...
void SetOrigin(const Point3D &origin)
Set the origin, i.e. the upper-left corner of the plane.
ScalarType SignedDistanceFromPlane(const Point3D &pt3d_mm) const
Signed distance of the point from the plane (bounding-box not considered)
void Modified() const override
Overload of function Modified() to prohibit several calls of Modified() using the ModifiedLock class...
Operation for setting a plane (defined by its origin and normal)
void SetReferenceGeometry(const mitk::BaseGeometry *geometry)
Set the geometrical frame of reference in which this PlaneGeometry is placed.
itk::AffineGeometryFrame< ScalarType, 3 >::TransformType AffineTransform3D
static T max(T x, T y)
Definition: svm.cpp:56
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.
Point3D ProjectPointOntoPlane(const Point3D &pt) const
Returns the lot from the point to the plane.
bool IsParallel(const PlaneGeometry *plane) const
Returns whether the plane is parallel to another plane.
MITKCORE_EXPORT const ScalarType sqrteps
Vector3D GetNormal() const
Normal of the plane.
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 SetBounds(const BoundsArrayType &bounds)
Set the bounding box (in index/unit coordinates)
itk::Point< TCoordRep, NPointDimension > GetPoint2() const
Get end point of the line.
Definition: mitkLine.h:117
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
itk::LightObject::Pointer InternalClone() const override
clones the geometry
void SetExtentInMM(int direction, ScalarType extentInMM)
Set the extent of the bounding-box in the specified direction in mm.
static int RectangleLineIntersection(TCoordRep x1, TCoordRep y1, TCoordRep x2, TCoordRep y2, itk::Point< TCoordRep, 2 > p, itk::Vector< TCoordRep, 2 > d, itk::Point< TCoordRep, 2 > &s1, itk::Point< TCoordRep, 2 > &s2)
Calculates the intersection points of a straight line in 2D with a rectangle.
Definition: mitkLine.h:260
void SetMatrixByVectors(const VnlVector &rightVector, const VnlVector &downVector, ScalarType thickness=1.0)
Initialize plane by right-/down-vector.
void CheckBounds(const BoundsArrayType &bounds) override
CheckBounds.
virtual void InitializePlane(const Point3D &origin, const Vector3D &normal)
Initialize plane by origin and normal (size is 1.0 mm in all directions, direction of right-/down-vec...
MITKCORE_EXPORT const ScalarType eps
ScalarType Distance(const Point3D &pt3d_mm) const
Distance of the point from the geometry (bounding-box not considered)
virtual void InitializeStandardPlane(const BaseGeometry *geometry3D, PlaneOrientation planeorientation=Axial, ScalarType zPosition=0, bool frontside=true, bool rotated=false, bool top=true)
Initialize a plane with orientation planeorientation (default: axial) with respect to BaseGeometry (d...
Describes a two-dimensional, rectangular plane.
OperationType GetOperationType()
const BoundsArrayType GetBounds() const
BaseGeometry Describes the geometry of a data object.
void ExecuteOperation(Operation *operation) override
executes affine operations (translate, rotate, scale)
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.
VnlVector GetMatrixColumn(unsigned int direction) const
Get a VnlVector along bounding-box in the specified direction, length is spacing. ...
BoundingBoxType::BoundsArrayType BoundsArrayType
virtual const BoundingBoxType * GetBoundingBox()