Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkSlicedGeometry3D.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 <itkSpatialOrientationAdapter.h>
14 
15 #include "mitkSlicedGeometry3D.h"
18 #include "mitkInteractionConst.h"
19 #include "mitkPlaneGeometry.h"
20 #include "mitkPlaneOperation.h"
22 #include "mitkRotationOperation.h"
24 
25 const mitk::ScalarType PI = 3.14159265359;
26 
28  : m_EvenlySpaced(true), m_Slices(0), m_ReferenceGeometry(nullptr), m_SliceNavigationController(nullptr)
29 {
30  m_DirectionVector.Fill(0);
32 }
33 
35  : Superclass(other),
37  m_Slices(other.m_Slices),
40 {
41  m_DirectionVector.Fill(0);
42  SetSpacing(other.GetSpacing());
44 
45  if (m_EvenlySpaced)
46  {
47  assert(!other.m_PlaneGeometries.empty() && "This may happen when you use one of the old Initialize methods, which had a bool parameter that is implicitly casted to the number of slices now.");
48  PlaneGeometry::Pointer geometry = other.m_PlaneGeometries[0]->Clone();
49  assert(geometry.IsNotNull());
50  SetPlaneGeometry(geometry, 0);
51  }
52  else
53  {
54  unsigned int s;
55  for (s = 0; s < other.m_Slices; ++s)
56  {
57  if (other.m_PlaneGeometries[s].IsNull())
58  {
59  assert(other.m_EvenlySpaced);
60  m_PlaneGeometries[s] = nullptr;
61  }
62  else
63  {
64  PlaneGeometry *geometry2D = other.m_PlaneGeometries[s]->Clone();
65  assert(geometry2D != nullptr);
66  SetPlaneGeometry(geometry2D, s);
67  }
68  }
69  }
70 }
71 
73 {
74 }
75 
77 {
78  mitk::PlaneGeometry::Pointer geometry2D = nullptr;
79 
80  if (this->IsValidSlice(s))
81  {
82  geometry2D = m_PlaneGeometries[s];
83 
84  // If (a) m_EvenlySpaced==true, (b) we don't have a PlaneGeometry stored
85  // for the requested slice, and (c) the first slice (s=0)
86  // is a PlaneGeometry instance, then we calculate the geometry of the
87  // requested as the plane of the first slice shifted by m_Spacing[2]*s
88  // in the direction of m_DirectionVector.
89  if ((m_EvenlySpaced) && (geometry2D.IsNull()))
90  {
91  PlaneGeometry *firstSlice = m_PlaneGeometries[0];
92 
93  if (firstSlice != nullptr &&
94  dynamic_cast<AbstractTransformGeometry *>(m_PlaneGeometries[0].GetPointer()) == nullptr)
95  {
96  if ((m_DirectionVector[0] == 0.0) && (m_DirectionVector[1] == 0.0) && (m_DirectionVector[2] == 0.0))
97  {
98  m_DirectionVector = firstSlice->GetNormal();
99  m_DirectionVector.Normalize();
100  }
101 
102  Vector3D direction;
103  direction = m_DirectionVector * this->GetSpacing()[2];
104 
105  mitk::PlaneGeometry::Pointer requestedslice;
106  requestedslice = static_cast<mitk::PlaneGeometry *>(firstSlice->Clone().GetPointer());
107 
108  requestedslice->SetOrigin(requestedslice->GetOrigin() + direction * s);
109 
110  geometry2D = requestedslice;
111  m_PlaneGeometries[s] = geometry2D;
112  }
113  }
114  return geometry2D;
115  }
116  else
117  {
118  return nullptr;
119  }
120 }
121 
123 {
124  assert(this->IsBoundingBoxNull() == false);
126 }
127 
129 {
130  if (this->IsValidSlice(s))
131  {
132  m_PlaneGeometries[s] = geometry2D;
133  m_PlaneGeometries[s]->SetReferenceGeometry(m_ReferenceGeometry);
134  return true;
135  }
136  return false;
137 }
138 
140 {
142  m_Slices = slices;
143 
144  PlaneGeometry::Pointer gnull = nullptr;
145  m_PlaneGeometries.assign(m_Slices, gnull);
146 
147  Vector3D spacing;
148  spacing.Fill(1.0);
149  this->SetSpacing(spacing);
150 
151  m_DirectionVector.Fill(0);
152 }
153 
155 {
156  assert(geometry2D != nullptr);
157  this->InitializeEvenlySpaced(geometry2D, geometry2D->GetExtentInMM(2) / geometry2D->GetExtent(2), slices);
158 }
159 
161  mitk::ScalarType zSpacing,
162  unsigned int slices)
163 {
164  assert(geometry2D != nullptr);
165  assert(geometry2D->GetExtent(0) > 0);
166  assert(geometry2D->GetExtent(1) > 0);
167 
168  geometry2D->Register();
169 
171  m_Slices = slices;
172 
173  BoundingBox::BoundsArrayType bounds = geometry2D->GetBounds();
174  bounds[4] = 0;
175  bounds[5] = slices;
176 
177  // clear and reserve
178  PlaneGeometry::Pointer gnull = nullptr;
179  m_PlaneGeometries.assign(m_Slices, gnull);
180 
181  Vector3D directionVector = geometry2D->GetAxisVector(2);
182  directionVector.Normalize();
183  directionVector *= zSpacing;
184 
185  // Normally we should use the following four lines to create a copy of
186  // the transform contrained in geometry2D, because it may not be changed
187  // by us. But we know that SetSpacing creates a new transform without
188  // changing the old (coming from geometry2D), so we can use the fifth
189  // line instead. We check this at (**).
190  //
191  // AffineTransform3D::Pointer transform = AffineTransform3D::New();
192  // transform->SetMatrix(geometry2D->GetIndexToWorldTransform()->GetMatrix());
193  // transform->SetOffset(geometry2D->GetIndexToWorldTransform()->GetOffset());
194  // SetIndexToWorldTransform(transform);
195 
197 
198  mitk::Vector3D spacing;
199  FillVector3D(spacing, geometry2D->GetExtentInMM(0) / bounds[1], geometry2D->GetExtentInMM(1) / bounds[3], zSpacing);
200 
201  this->SetDirectionVector(directionVector);
202  this->SetBounds(bounds);
203  this->SetPlaneGeometry(geometry2D, 0);
204  this->SetSpacing(spacing, true);
205  this->SetEvenlySpaced();
206 
207  // this->SetTimeBounds( geometry2D->GetTimeBounds() );
208 
209  assert(this->GetIndexToWorldTransform() != geometry2D->GetIndexToWorldTransform()); // (**) see above.
210 
211  this->SetFrameOfReferenceID(geometry2D->GetFrameOfReferenceID());
212  this->SetImageGeometry(geometry2D->GetImageGeometry());
213 
214  geometry2D->UnRegister();
215 }
216 
218  mitk::PlaneGeometry::PlaneOrientation planeorientation,
219  bool top,
220  bool frontside,
221  bool rotated)
222 {
223  m_ReferenceGeometry = geometry3D;
224 
226  planeGeometry->InitializeStandardPlane(geometry3D, top, planeorientation, frontside, rotated);
227 
228  int worldAxis =
229  planeorientation == PlaneGeometry::Sagittal ? 0 :
230  planeorientation == PlaneGeometry::Frontal ? 1 : 2;
231 
232  // Inspired by:
233  // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
234 
235  mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
236  matrix.GetVnlMatrix().normalize_columns();
237  mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetTranspose();
238 
239  int dominantAxis = planeGeometry->CalculateDominantAxes(inverseMatrix).at(worldAxis);
240  ScalarType viewSpacing = geometry3D->GetSpacing()[dominantAxis];
241 
245  auto slices = static_cast<unsigned int>(geometry3D->GetExtent(dominantAxis) + 0.5);
246  if ( slices == 0 && geometry3D->GetExtent(dominantAxis) > 0) {
247  // require at least one slice if there is _some_ extent
248  slices = 1;
249  }
250 
251 #ifndef NDEBUG
252  int upDirection = itk::Function::Sign(inverseMatrix[dominantAxis][worldAxis]);
253 
257  Vector3D worldPlaneNormal = inverseMatrix.get_row(dominantAxis) * (upDirection * viewSpacing);
258 
260  Vector3D standardPlaneNormal = planeGeometry->GetNormal();
261 
265  assert((standardPlaneNormal - (top ? 1.0 : -1.0) * worldPlaneNormal).GetSquaredNorm() < 0.000001);
266 #endif
267 
268  this->InitializeEvenlySpaced(planeGeometry, viewSpacing, slices);
269 
270 #ifndef NDEBUG
271  Vector3D zAxisVector = this->GetAxisVector(2);
274  Vector3D upscaledStandardPlaneNormal = standardPlaneNormal;
275  upscaledStandardPlaneNormal *= slices;
276  assert((zAxisVector - upscaledStandardPlaneNormal).GetSquaredNorm() < 0.000001);
277 
282 // ScalarType det = vnl_det(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix());
283 // MITK_DEBUG << "world axis: " << worldAxis << (det > 0 ? " ; right-handed" : " ; left-handed");
284 #endif
285 }
286 
287 void mitk::SlicedGeometry3D::ReinitializePlanes(const Point3D &center, const Point3D &referencePoint)
288 {
289  // Need a reference frame to align the rotated planes
290  if (!m_ReferenceGeometry)
291  {
292  return;
293  }
294 
295  // Get first plane of plane stack
296  PlaneGeometry *firstPlane = m_PlaneGeometries[0];
297 
298  // If plane stack is empty, exit
299  if (!firstPlane || dynamic_cast<AbstractTransformGeometry *>(firstPlane))
300  {
301  return;
302  }
303 
304  // Calculate the "directed" spacing when taking the plane (defined by its axes
305  // vectors and normal) as the reference coordinate frame.
306  //
307  // This is done by calculating the radius of the ellipsoid defined by the
308  // original volume spacing axes, in the direction of the respective axis of the
309  // reference frame.
310  mitk::Vector3D axis0 = firstPlane->GetAxisVector(0);
311  mitk::Vector3D axis1 = firstPlane->GetAxisVector(1);
312  mitk::Vector3D normal = firstPlane->GetNormal();
313  normal.Normalize();
314 
315  Vector3D spacing;
316  spacing[0] = this->CalculateSpacing(axis0);
317  spacing[1] = this->CalculateSpacing(axis1);
318  spacing[2] = this->CalculateSpacing(normal);
319 
320  Superclass::SetSpacing(spacing);
321 
322  // Now we need to calculate the number of slices in the plane's normal
323  // direction, so that the entire volume is covered. This is done by first
324  // calculating the dot product between the volume diagonal (the maximum
325  // distance inside the volume) and the normal, and dividing this value by
326  // the directed spacing calculated above.
327  ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * normal[0]) +
328  std::abs(m_ReferenceGeometry->GetExtentInMM(1) * normal[1]) +
329  std::abs(m_ReferenceGeometry->GetExtentInMM(2) * normal[2]);
330 
331  if (directedExtent >= spacing[2])
332  {
333  m_Slices = static_cast<unsigned int>(directedExtent / spacing[2] + 0.5);
334  }
335  else
336  {
337  m_Slices = 1;
338  }
339 
340  // The origin of our "first plane" needs to be adapted to this new extent.
341  // To achieve this, we first calculate the current distance to the volume's
342  // center, and then shift the origin in the direction of the normal by the
343  // difference between this distance and half of the new extent.
344  double centerOfRotationDistance = firstPlane->SignedDistanceFromPlane(center);
345 
346  if (centerOfRotationDistance > 0)
347  {
348  firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (centerOfRotationDistance - directedExtent / 2.0));
349  m_DirectionVector = normal;
350  }
351  else
352  {
353  firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (directedExtent / 2.0 + centerOfRotationDistance));
354  m_DirectionVector = -normal;
355  }
356 
357  // Now we adjust this distance according with respect to the given reference
358  // point: we need to make sure that the point is touched by one slice of the
359  // new slice stack.
360  double referencePointDistance = firstPlane->SignedDistanceFromPlane(referencePoint);
361 
362  auto referencePointSlice = static_cast<int>(referencePointDistance / spacing[2]);
363 
364  double alignmentValue = referencePointDistance / spacing[2] - referencePointSlice;
365 
366  firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * alignmentValue * spacing[2]);
367 
368  // Finally, we can clear the previous geometry stack and initialize it with
369  // our re-initialized "first plane".
371 
372  if (m_Slices > 0)
373  {
374  m_PlaneGeometries[0] = firstPlane;
375  }
376 
377  // Reinitialize SNC with new number of slices
379 
380  this->Modified();
381 }
382 
384 {
385  // Need the spacing of the underlying dataset / geometry
386  if (!m_ReferenceGeometry)
387  {
388  return 1.0;
389  }
390 
391  const mitk::Vector3D &spacing = m_ReferenceGeometry->GetSpacing();
392  return SlicedGeometry3D::CalculateSpacing(spacing, d);
393 }
394 
396 {
397  // The following can be derived from the ellipsoid equation
398  //
399  // 1 = x^2/a^2 + y^2/b^2 + z^2/c^2
400  //
401  // where (a,b,c) = spacing of original volume (ellipsoid radii)
402  // and (x,y,z) = scaled coordinates of vector d (according to ellipsoid)
403  //
404  double scaling = d[0] * d[0] / (spacing[0] * spacing[0]) + d[1] * d[1] / (spacing[1] * spacing[1]) +
405  d[2] * d[2] / (spacing[2] * spacing[2]);
406 
407  scaling = sqrt(scaling);
408 
409  return (sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]) / scaling);
410 }
411 
413 {
414  TransformType::Pointer inverse = TransformType::New();
415  m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse(inverse);
416 
417  Vector3D transformedNormal = inverse->TransformVector(normal);
418 
419  transformedNormal.Normalize();
420  return transformedNormal;
421 }
422 
423 void mitk::SlicedGeometry3D::SetImageGeometry(const bool isAnImageGeometry)
424 {
425  Superclass::SetImageGeometry(isAnImageGeometry);
426 
427  unsigned int s;
428  for (s = 0; s < m_Slices; ++s)
429  {
430  mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
431  if (geometry != nullptr)
432  {
433  geometry->SetImageGeometry(isAnImageGeometry);
434  }
435  }
436 }
437 
439 {
440  unsigned int s;
441  for (s = 0; s < m_Slices; ++s)
442  {
443  mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
444  if (geometry != nullptr)
445  {
446  geometry->ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry);
447  }
448  }
449 
451 }
452 
454 {
455  return ((s >= 0) && (s < (int)m_Slices));
456 }
457 
459 {
460  return m_ReferenceGeometry;
461 }
462 
464 {
465  m_ReferenceGeometry = referenceGeometry;
466 
467  std::vector<PlaneGeometry::Pointer>::iterator it;
468 
469  for (it = m_PlaneGeometries.begin(); it != m_PlaneGeometries.end(); ++it)
470  {
471  (*it)->SetReferenceGeometry(referenceGeometry);
472  }
473 }
474 
476 {
477  return ( m_ReferenceGeometry != nullptr );
478 }
479 
481 {
482  bool hasEvenlySpacedPlaneGeometry = false;
483  mitk::Point3D origin;
484  mitk::Vector3D rightDV, bottomDV;
486 
487  // Check for valid spacing
488  if (!(aSpacing[0] > 0 && aSpacing[1] > 0 && aSpacing[2] > 0))
489  {
490  mitkThrow() << "You try to set a spacing with at least one element equal or "
491  "smaller to \"0\". This might lead to a crash during rendering. Please double"
492  " check your data!";
493  }
494 
495  // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
496  // since the spacing influences them
497  if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
498  {
499  const PlaneGeometry *planeGeometry = m_PlaneGeometries[0];
500 
501  if (planeGeometry && !dynamic_cast<const AbstractTransformGeometry *>(planeGeometry))
502  {
503  this->WorldToIndex(planeGeometry->GetOrigin(), origin);
504  this->WorldToIndex(planeGeometry->GetAxisVector(0), rightDV);
505  this->WorldToIndex(planeGeometry->GetAxisVector(1), bottomDV);
506 
507  bounds = planeGeometry->GetBounds();
508  hasEvenlySpacedPlaneGeometry = true;
509  }
510  }
511 
512  BaseGeometry::_SetSpacing(aSpacing);
513 
514  mitk::PlaneGeometry::Pointer firstGeometry;
515 
516  // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
517  // since the spacing influences them
518  if (hasEvenlySpacedPlaneGeometry)
519  {
520  // create planeGeometry according to new spacing
521  this->IndexToWorld(origin, origin);
522  this->IndexToWorld(rightDV, rightDV);
523  this->IndexToWorld(bottomDV, bottomDV);
524 
526  planeGeometry->SetImageGeometry(this->GetImageGeometry());
527 
528  planeGeometry->SetReferenceGeometry(m_ReferenceGeometry);
529 
530  // Store spacing, as Initialize... needs a pointer
531  mitk::Vector3D lokalSpacing = this->GetSpacing();
532  planeGeometry->InitializeStandardPlane(rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &lokalSpacing);
533  planeGeometry->SetOrigin(origin);
534  planeGeometry->SetBounds(bounds);
535 
536  firstGeometry = planeGeometry;
537  }
538  else if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
539  {
540  firstGeometry = m_PlaneGeometries[0].GetPointer();
541  }
542 
543  // clear and reserve
544  PlaneGeometry::Pointer gnull = nullptr;
545  m_PlaneGeometries.assign(m_Slices, gnull);
546 
547  if (m_Slices > 0)
548  {
549  m_PlaneGeometries[0] = firstGeometry;
550  }
551 
552  this->Modified();
553 }
554 
556 {
558 }
559 
561 {
563 }
564 
566 {
567  if (m_EvenlySpaced != on)
568  {
569  m_EvenlySpaced = on;
570  this->Modified();
571  }
572 }
573 
575 {
576  Vector3D newDir = directionVector;
577  newDir.Normalize();
578  if (newDir != m_DirectionVector)
579  {
580  m_DirectionVector = newDir;
581  this->Modified();
582  }
583 }
584 
585 // void
586 // mitk::SlicedGeometry3D::SetTimeBounds( const mitk::TimeBounds& timebounds )
587 //{
588 // Superclass::SetTimeBounds( timebounds );
589 //
590 // unsigned int s;
591 // for ( s = 0; s < m_Slices; ++s )
592 // {
593 // if(m_Geometry2Ds[s].IsNotNull())
594 // {
595 // m_Geometry2Ds[s]->SetTimeBounds( timebounds );
596 // }
597 // }
598 // m_TimeBounds = timebounds;
599 //}
600 
601 itk::LightObject::Pointer mitk::SlicedGeometry3D::InternalClone() const
602 {
603  Self::Pointer newGeometry = new SlicedGeometry3D(*this);
604  newGeometry->UnRegister();
605  return newGeometry.GetPointer();
606 }
607 
608 void mitk::SlicedGeometry3D::PrintSelf(std::ostream &os, itk::Indent indent) const
609 {
610  Superclass::PrintSelf(os, indent);
611  os << indent << " EvenlySpaced: " << m_EvenlySpaced << std::endl;
612  if (m_EvenlySpaced)
613  {
614  os << indent << " DirectionVector: " << m_DirectionVector << std::endl;
615  }
616  os << indent << " Slices: " << m_Slices << std::endl;
617 
618  os << std::endl;
619  os << indent << " GetPlaneGeometry(0): ";
620  if (this->GetPlaneGeometry(0) == nullptr)
621  {
622  os << "nullptr" << std::endl;
623  }
624  else
625  {
626  this->GetPlaneGeometry(0)->Print(os, indent);
627  }
628 }
629 
631 {
632  PlaneGeometry::Pointer geometry2D;
633  ApplyTransformMatrixOperation *applyMatrixOp;
634  Point3D center;
635 
636  switch (operation->GetOperationType())
637  {
638  case OpNOTHING:
639  break;
640 
641  case OpROTATE:
642  if (m_EvenlySpaced)
643  {
644  // Need a reference frame to align the rotation
646  {
647  // Clear all generated geometries and then rotate only the first slice.
648  // The other slices will be re-generated on demand
649 
650  // Save first slice
652 
653  auto *rotOp = dynamic_cast<RotationOperation *>(operation);
654 
655  // Generate a RotationOperation using the dataset center instead of
656  // the supplied rotation center. This is necessary so that the rotated
657  // zero-plane does not shift away. The supplied center is instead used
658  // to adjust the slice stack afterwards.
660 
661  RotationOperation centeredRotation(
662  rotOp->GetOperationType(), center, rotOp->GetVectorOfRotation(), rotOp->GetAngleOfRotation());
663 
664  // Rotate first slice
665  geometry2D->ExecuteOperation(&centeredRotation);
666 
667  // Clear the slice stack and adjust it according to the center of
668  // the dataset and the supplied rotation center (see documentation of
669  // ReinitializePlanes)
670  this->ReinitializePlanes(center, rotOp->GetCenterOfRotation());
671 
672  geometry2D->SetSpacing(this->GetSpacing());
673 
675  {
676  m_SliceNavigationController->SelectSliceByPoint(rotOp->GetCenterOfRotation());
678  }
679 
680  BaseGeometry::ExecuteOperation(&centeredRotation);
681  }
682  else
683  {
684  // we also have to consider the case, that there is no reference geometry available.
685  if (m_PlaneGeometries.size() > 0)
686  {
687  // Reach through to all slices in my container
688  for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
689  {
690  // Test for empty slices, which can happen if evenly spaced geometry
691  if ((*iter).IsNotNull())
692  {
693  (*iter)->ExecuteOperation(operation);
694  }
695  }
696 
697  // rotate overall geometry
698  auto *rotOp = dynamic_cast<RotationOperation *>(operation);
700  }
701  }
702  }
703  else
704  {
705  // Reach through to all slices
706  for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
707  {
708  (*iter)->ExecuteOperation(operation);
709  }
710  }
711  break;
712 
713  case OpORIENT:
714  if (m_EvenlySpaced)
715  {
716  // get operation data
717  auto *planeOp = dynamic_cast<PlaneOperation *>(operation);
718 
719  // Get first slice
720  PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
721 
722  // Need a PlaneGeometry, a PlaneOperation and a reference frame to
723  // carry out the re-orientation. If not all avaialble, stop here
724  if (!m_ReferenceGeometry ||
725  (!planeGeometry || dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer())) || !planeOp)
726  {
727  break;
728  }
729 
730  // General Behavior:
731  // Clear all generated geometries and then rotate only the first slice.
732  // The other slices will be re-generated on demand
733 
734  //
735  // 1st Step: Reorient Normal Vector of first plane
736  //
737  Point3D center = planeOp->GetPoint(); // m_ReferenceGeometry->GetCenter();
738  mitk::Vector3D currentNormal = planeGeometry->GetNormal();
739  mitk::Vector3D newNormal;
740  if (planeOp->AreAxisDefined())
741  {
742  // If planeOp was defined by one centerpoint and two axis vectors
743  newNormal = CrossProduct(planeOp->GetAxisVec0(), planeOp->GetAxisVec1());
744  }
745  else
746  {
747  // If planeOp was defined by one centerpoint and one normal vector
748  newNormal = planeOp->GetNormal();
749  }
750 
751  // Get Rotation axis und angle
752  currentNormal.Normalize();
753  newNormal.Normalize();
754  ScalarType rotationAngle = angle(currentNormal.GetVnlVector(), newNormal.GetVnlVector());
755 
756  rotationAngle *= 180.0 / vnl_math::pi; // from rad to deg
757  Vector3D rotationAxis = itk::CrossProduct(currentNormal, newNormal);
758  if (std::abs(rotationAngle - 180) < mitk::eps)
759  {
760  // current Normal and desired normal are not linear independent!!(e.g 1,0,0 and -1,0,0).
761  // Rotation Axis should be ANY vector that is 90� to current Normal
762  mitk::Vector3D helpNormal;
763  helpNormal = currentNormal;
764  helpNormal[0] += 1;
765  helpNormal[1] -= 1;
766  helpNormal[2] += 1;
767  helpNormal.Normalize();
768  rotationAxis = itk::CrossProduct(helpNormal, currentNormal);
769  }
770 
771  RotationOperation centeredRotation(mitk::OpROTATE, center, rotationAxis, rotationAngle);
772 
773  // Rotate first slice
774  planeGeometry->ExecuteOperation(&centeredRotation);
775 
776  // Reinitialize planes and select slice, if my rotations are all done.
777  if (!planeOp->AreAxisDefined())
778  {
779  // Clear the slice stack and adjust it according to the center of
780  // rotation and plane position (see documentation of ReinitializePlanes)
781  this->ReinitializePlanes(center, planeOp->GetPoint());
782  planeGeometry->SetSpacing(this->GetSpacing());
783 
785  {
786  m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
788  }
789  }
790 
791  // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
792  BaseGeometry::ExecuteOperation(&centeredRotation);
793 
794  //
795  // 2nd step. If axis vectors were defined, rotate the plane around its normal to fit these
796  //
797 
798  if (planeOp->AreAxisDefined())
799  {
800  mitk::Vector3D vecAxixNew = planeOp->GetAxisVec0();
801  vecAxixNew.Normalize();
802  mitk::Vector3D VecAxisCurr = planeGeometry->GetAxisVector(0);
803  VecAxisCurr.Normalize();
804 
805  ScalarType rotationAngle = angle(VecAxisCurr.GetVnlVector(), vecAxixNew.GetVnlVector());
806  rotationAngle = rotationAngle * 180 / PI; // Rad to Deg
807 
808  // we rotate around the normal of the plane, but we do not know, if we need to rotate clockwise
809  // or anti-clockwise. So we rotate around the crossproduct of old and new Axisvector.
810  // Since both axis vectors lie in the plane, the crossproduct is the planes normal or the negative planes
811  // normal
812 
813  rotationAxis = itk::CrossProduct(VecAxisCurr, vecAxixNew);
814  if (std::abs(rotationAngle - 180) < mitk::eps)
815  {
816  // current axisVec and desired axisVec are not linear independent!!(e.g 1,0,0 and -1,0,0).
817  // Rotation Axis can be just plane Normal. (have to rotate by 180�)
818  rotationAxis = newNormal;
819  }
820 
821  // Perfom Rotation
822  mitk::RotationOperation op(mitk::OpROTATE, center, rotationAxis, rotationAngle);
823  planeGeometry->ExecuteOperation(&op);
824 
825  // Apply changes on first slice to whole slice stack
826  this->ReinitializePlanes(center, planeOp->GetPoint());
827  planeGeometry->SetSpacing(this->GetSpacing());
828 
830  {
831  m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
833  }
834 
835  // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
837  }
838  }
839  else
840  {
841  // Reach through to all slices
842  for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
843  {
844  (*iter)->ExecuteOperation(operation);
845  }
846  }
847  break;
848 
850  if (m_EvenlySpaced)
851  {
852  // Save first slice
853  PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
854 
855  auto *restorePlaneOp = dynamic_cast<RestorePlanePositionOperation *>(operation);
856 
857  // Need a PlaneGeometry, a PlaneOperation and a reference frame to
858  // carry out the re-orientation
859  if (m_ReferenceGeometry &&
860  (planeGeometry && dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer()) == nullptr) &&
861  restorePlaneOp)
862  {
863  // Clear all generated geometries and then rotate only the first slice.
864  // The other slices will be re-generated on demand
865 
866  // Rotate first slice
867  planeGeometry->ExecuteOperation(restorePlaneOp);
868 
869  m_DirectionVector = restorePlaneOp->GetDirectionVector();
870 
871  double centerOfRotationDistance = planeGeometry->SignedDistanceFromPlane(m_ReferenceGeometry->GetCenter());
872 
873  if (centerOfRotationDistance <= 0)
874  {
876  }
877 
878  Vector3D spacing = restorePlaneOp->GetSpacing();
879 
880  Superclass::SetSpacing(spacing);
881 
882  // /*Now we need to calculate the number of slices in the plane's normal
883  // direction, so that the entire volume is covered. This is done by first
884  // calculating the dot product between the volume diagonal (the maximum
885  // distance inside the volume) and the normal, and dividing this value by
886  // the directed spacing calculated above.*/
887  ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * m_DirectionVector[0]) +
890 
891  if (directedExtent >= spacing[2])
892  {
893  m_Slices = static_cast<unsigned int>(directedExtent / spacing[2] + 0.5);
894  }
895  else
896  {
897  m_Slices = 1;
898  }
899 
901 
902  if (m_Slices > 0)
903  {
904  m_PlaneGeometries[0] = planeGeometry;
905  }
906 
908 
909  this->Modified();
910 
911  // End Reinitialization
912 
914  {
915  m_SliceNavigationController->GetSlice()->SetPos(restorePlaneOp->GetPos());
917  }
918  BaseGeometry::ExecuteOperation(restorePlaneOp);
919  }
920  }
921  else
922  {
923  // Reach through to all slices
924  for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
925  {
926  (*iter)->ExecuteOperation(operation);
927  }
928  }
929  break;
930 
932 
933  // Clear all generated geometries and then transform only the first slice.
934  // The other slices will be re-generated on demand
935 
936  // Save first slice
937  geometry2D = m_PlaneGeometries[0];
938 
939  applyMatrixOp = dynamic_cast<ApplyTransformMatrixOperation *>(operation);
940 
941  // Apply transformation to first plane
942  geometry2D->ExecuteOperation(applyMatrixOp);
943 
944  // Generate a ApplyTransformMatrixOperation using the dataset center instead of
945  // the supplied rotation center. The supplied center is instead used to adjust the
946  // slice stack afterwards (see OpROTATE).
947  center = m_ReferenceGeometry->GetCenter();
948 
949  // Clear the slice stack and adjust it according to the center of
950  // the dataset and the supplied rotation center (see documentation of
951  // ReinitializePlanes)
952  this->ReinitializePlanes(center, applyMatrixOp->GetReferencePoint());
953 
954  BaseGeometry::ExecuteOperation(applyMatrixOp);
955  break;
956 
957  default: // let handle by base class if we don't do anything
959  }
960 
961  this->Modified();
962 }
virtual void InitializePlanes(const mitk::BaseGeometry *geometry3D, mitk::PlaneGeometry::PlaneOrientation planeorientation, bool top=true, bool frontside=true, bool rotated=false)
Completely initialize this instance as evenly-spaced plane slices parallel to a side of the provided ...
Stepper * GetSlice()
Get the Stepper through the slices.
mitk::Vector3D AdjustNormal(const mitk::Vector3D &normal) const
itk::BoundingBox< unsigned long, 3, ScalarType > BoundingBox
Standard 3D-BoundingBox typedef.
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?
mitk::SliceNavigationController * m_SliceNavigationController
void SelectSliceByPoint(const mitk::Point3D &point)
Positions the SNC according to the specified point.
void SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
Set the spacing (m_Spacing).
virtual void InitializeSlicedGeometry(unsigned int slices)
Tell this instance how many PlaneGeometries it shall manage. Bounding box and the PlaneGeometries mus...
void SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
bool IsBoundingBoxNull() const
Vector3D GetAxisVector(unsigned int direction) const
Get vector along bounding-box in the specified direction in mm.
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Base class of all Operation-classes.
Definition: mitkOperation.h:29
double ScalarType
void _SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
virtual void ReinitializePlanes(const Point3D &center, const Point3D &referencePoint)
virtual void SetEvenlySpaced(bool on=true)
ScalarType GetExtent(unsigned int direction) const
Set the time bounds (in ms)
mitk::SliceNavigationController * GetSliceNavigationController()
Pointer Clone() const
void Initialize()
Initialize the BaseGeometry.
virtual bool SetPlaneGeometry(mitk::PlaneGeometry *geometry2D, int s)
Set PlaneGeometry of slice s.
Constants for most interaction classes, due to the generic StateMachines.
virtual const mitk::Vector3D & GetDirectionVector() const
Point3D GetCenter() const
Get the center of the bounding-box in mm.
itk::LightObject::Pointer InternalClone() const override
clones the geometry
void PreSetSpacing(const mitk::Vector3D &aSpacing) override
PreSetSpacing.
virtual void InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, unsigned int slices)
Completely initialize this instance as evenly-spaced with slices parallel to the provided PlaneGeomet...
virtual void SetSliceNavigationController(mitk::SliceNavigationController *snc)
Set the SliceNavigationController corresponding to this sliced geometry.
std::vector< PlaneGeometry::Pointer > m_PlaneGeometries
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.
Controls the selection of the slice the associated BaseRenderer will display.
const mitk::BoundingBox * GetBoundingBox() const override
const mitk::ScalarType PI
virtual mitk::PlaneGeometry * GetPlaneGeometry(int s) const
Returns the PlaneGeometry of the slice (s).
void SetImageGeometry(const bool isAnImageGeometry) override
Define that this BaseGeometry is refering to an Image.
void PrintSelf(std::ostream &os, itk::Indent indent) const override
void ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry) override
When switching from an Image Geometry to a normal Geometry (and the other way around), you have to.
#define mitkThrow()
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)
virtual bool IsValidSlice(int s=0) const
Check whether a slice exists.
const mitk::BaseGeometry * m_ReferenceGeometry
virtual void SetPos(unsigned int pos)
Definition: mitkStepper.h:55
virtual unsigned int GetFrameOfReferenceID() const
Get the DICOM FrameOfReferenceID referring to the used world coordinate system.
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.
Describes the geometry of a data object consisting of slices.
void AdjustSliceStepperRange()
Adjusts the numerical range of the slice stepper according to the current geometry orientation of thi...
Vector3D GetNormal() const
Normal of the plane.
virtual void ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry)
When switching from an Image Geometry to a normal Geometry (and the other way around), you have to.
void SetBounds(const BoundsArrayType &bounds)
Set the bounding box (in index/unit coordinates)
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
static double CalculateSpacing(const mitk::Vector3D &spacing, const mitk::Vector3D &d)
static Pointer New()
virtual const BaseGeometry * GetReferenceGeometry() const
virtual void SetFrameOfReferenceID(unsigned int _arg)
Set the DICOM FrameOfReferenceID referring to the used world coordinate system.
MITKCORE_EXPORT const ScalarType eps
Describes a two-dimensional, rectangular plane.
OperationType GetOperationType()
Operation, that holds everything necessary for an rotation operation on mitk::BaseData.
virtual void SetDirectionVector(const mitk::Vector3D &directionVector)
Set/Get the vector between slices for the evenly-spaced case (m_EvenlySpaced==true).
virtual void SetReferenceGeometry(const BaseGeometry *referenceGeometry)
void ExecuteOperation(Operation *operation) override
executes affine operations (translate, rotate, scale)
const BoundsArrayType GetBounds() const
virtual void SetImageGeometry(bool _arg)
Define that this BaseGeometry is refering to an Image.
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.
BoundingBoxType::BoundsArrayType BoundsArrayType
virtual void SetSteps(unsigned int _arg)
virtual const BoundingBoxType * GetBoundingBox()