Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkBaseGeometry.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 <iomanip>
14 #include <sstream>
15 #include <bitset>
16 
17 #include <vtkMatrix4x4.h>
18 #include <vtkMatrixToLinearTransform.h>
19 
21 #include "mitkBaseGeometry.h"
23 #include "mitkInteractionConst.h"
24 #include "mitkMatrixConvert.h"
25 #include "mitkModifiedLock.h"
26 #include "mitkPointOperation.h"
28 #include "mitkRotationOperation.h"
29 #include "mitkScaleOperation.h"
30 #include "mitkVector.h"
31 #include "mitkMatrix.h"
32 
34  : Superclass(),
36  m_FrameOfReferenceID(0),
37  m_IndexToWorldTransformLastModified(0),
38  m_ImageGeometry(false),
39  m_ModifiedLockFlag(false),
40  m_ModifiedCalledFlag(false)
41 {
42  m_GeometryTransform = new GeometryTransformHolder();
43  Initialize();
44 }
45 
47  : Superclass(),
49  m_FrameOfReferenceID(other.m_FrameOfReferenceID),
50  m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified),
51  m_ImageGeometry(other.m_ImageGeometry),
52  m_ModifiedLockFlag(false),
53  m_ModifiedCalledFlag(false)
54 {
55  m_GeometryTransform = new GeometryTransformHolder(*other.GetGeometryTransformHolder());
56  other.InitializeGeometry(this);
57 }
58 
60 {
61  delete m_GeometryTransform;
62 }
63 
64 void mitk::BaseGeometry::SetVtkMatrixDeepCopy(vtkTransform *vtktransform)
65 {
66  m_GeometryTransform->SetVtkMatrixDeepCopy(vtktransform);
67 }
68 
70 {
71  return m_GeometryTransform->GetOrigin();
72 }
73 
75 {
76  mitk::ModifiedLock lock(this);
77 
78  if (origin != GetOrigin())
79  {
80  m_GeometryTransform->SetOrigin(origin);
81  Modified();
82  }
83 }
84 
86 {
87  return m_GeometryTransform->GetSpacing();
88 }
89 
91 {
92  float b[6] = {0, 1, 0, 1, 0, 1};
93  SetFloatBounds(b);
94 
95  m_GeometryTransform->Initialize();
96 
97  m_FrameOfReferenceID = 0;
98 
99  m_ImageGeometry = false;
100 }
101 
102 void mitk::BaseGeometry::SetFloatBounds(const float bounds[6])
103 {
105  const float *input = bounds;
106  int i = 0;
107  for (mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6; ++i)
108  *it++ = (mitk::ScalarType)*input++;
109  SetBounds(b);
110 }
111 
112 void mitk::BaseGeometry::SetFloatBounds(const double bounds[6])
113 {
115  const double *input = bounds;
116  int i = 0;
117  for (mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6; ++i)
118  *it++ = (mitk::ScalarType)*input++;
119  SetBounds(b);
120 }
121 
124 {
125  newGeometry->SetBounds(m_BoundingBox->GetBounds());
126 
128 
129  newGeometry->InitializeGeometryTransformHolder(this);
130 
131  newGeometry->m_ImageGeometry = m_ImageGeometry;
132 }
133 
134 void mitk::BaseGeometry::InitializeGeometryTransformHolder(const BaseGeometry *otherGeometry)
135 {
136  this->m_GeometryTransform->Initialize(otherGeometry->GetGeometryTransformHolder());
137 }
138 
141 {
142  mitk::ModifiedLock lock(this);
143 
144  this->CheckBounds(bounds);
145 
146  m_BoundingBox = BoundingBoxType::New();
147 
148  BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New();
149  BoundingBoxType::PointType p;
150  BoundingBoxType::PointIdentifier pointid;
151 
152  for (pointid = 0; pointid < 2; ++pointid)
153  {
154  unsigned int i;
155  for (i = 0; i < m_NDimensions; ++i)
156  {
157  p[i] = bounds[2 * i + pointid];
158  }
159  pointscontainer->InsertElement(pointid, p);
160  }
161 
162  m_BoundingBox->SetPoints(pointscontainer);
163  m_BoundingBox->ComputeBoundingBox();
164  this->Modified();
165 }
166 
168 {
169  mitk::ModifiedLock lock(this);
170 
171  CheckIndexToWorldTransform(transform);
172 
173  m_GeometryTransform->SetIndexToWorldTransform(transform);
174  Modified();
175 }
176 
178 {
179  // security check
180  mitk::Vector3D originalSpacing = this->GetSpacing();
181 
182  mitk::ModifiedLock lock(this);
183 
184  CheckIndexToWorldTransform(transform);
185 
186  m_GeometryTransform->SetIndexToWorldTransformWithoutChangingSpacing(transform);
187  Modified();
188 
189  // Security check. Spacig must not have changed
190  if (!mitk::Equal(originalSpacing, this->GetSpacing()))
191  {
192  MITK_WARN << "Spacing has changed in a method, where the spacing must not change.";
193  assert(false);
194  }
195 }
196 
198 {
199  assert(m_BoundingBox.IsNotNull());
200  return m_BoundingBox->GetBounds();
201 }
202 
204 {
205  return true;
206 }
207 
208 void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing)
209 {
210  PreSetSpacing(aSpacing);
211  _SetSpacing(aSpacing, enforceSetSpacing);
212 }
213 
214 void mitk::BaseGeometry::_SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing)
215 {
216  m_GeometryTransform->SetSpacing(aSpacing, enforceSetSpacing);
217 }
218 
220 {
221  Vector3D frontToBack;
222  frontToBack.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction));
223  frontToBack *= GetExtent(direction);
224  return frontToBack;
225 }
226 
227 mitk::ScalarType mitk::BaseGeometry::GetExtent(unsigned int direction) const
228 {
229  assert(m_BoundingBox.IsNotNull());
230  if (direction >= m_NDimensions)
231  mitkThrow() << "Direction is too big. This geometry is for 3D Data";
232  BoundsArrayType bounds = m_BoundingBox->GetBounds();
233  return bounds[direction * 2 + 1] - bounds[direction * 2];
234 }
235 
237 {
238  bool isConvertableWithoutLoss = true;
239  do
240  {
241  if (this->GetSpacing()[2] != 1)
242  {
243  isConvertableWithoutLoss = false;
244  break;
245  }
246  if (this->GetOrigin()[2] != 0)
247  {
248  isConvertableWithoutLoss = false;
249  break;
250  }
251  mitk::Vector3D col0, col1, col2;
252  col0.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0));
253  col1.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1));
254  col2.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2));
255 
256  if ((col0[2] != 0) || (col1[2] != 0) || (col2[0] != 0) || (col2[1] != 0) || (col2[2] != 1))
257  {
258  isConvertableWithoutLoss = false;
259  break;
260  }
261  } while (false);
262 
263  return isConvertableWithoutLoss;
264 }
265 
267 {
268  assert(m_BoundingBox.IsNotNull());
269  Point3D c = m_BoundingBox->GetCenter();
270  if (m_ImageGeometry)
271  {
272  // Get Center returns the middel of min and max pixel index. In corner based images, this is the right position.
273  // In center based images (imageGeometry == true), the index needs to be shifted back.
274  c[0] -= 0.5;
275  c[1] -= 0.5;
276  c[2] -= 0.5;
277  }
278  this->IndexToWorld(c, c);
279  return c;
280 }
281 
283 {
284  Vector3D diagonalvector = GetCornerPoint() - GetCornerPoint(false, false, false);
285  return diagonalvector.GetSquaredNorm();
286 }
287 
289 {
290  return sqrt(GetDiagonalLength2());
291 }
292 
294 {
295  assert(id >= 0);
296  assert(this->IsBoundingBoxNull() == false);
297 
298  BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds();
299 
300  Point3D cornerpoint;
301  switch (id)
302  {
303  case 0:
304  FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[4]);
305  break;
306  case 1:
307  FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[5]);
308  break;
309  case 2:
310  FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[4]);
311  break;
312  case 3:
313  FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[5]);
314  break;
315  case 4:
316  FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[4]);
317  break;
318  case 5:
319  FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[5]);
320  break;
321  case 6:
322  FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[4]);
323  break;
324  case 7:
325  FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[5]);
326  break;
327  default:
328  {
329  itkExceptionMacro(<< "A cube only has 8 corners. These are labeled 0-7.");
330  }
331  }
332  if (m_ImageGeometry)
333  {
334  // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the
335  // bounding box. The bounding box itself is no image, so it is corner-based
336  FillVector3D(cornerpoint, cornerpoint[0] - 0.5, cornerpoint[1] - 0.5, cornerpoint[2] - 0.5);
337  }
338  return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint);
339 }
340 
341 mitk::Point3D mitk::BaseGeometry::GetCornerPoint(bool xFront, bool yFront, bool zFront) const
342 {
343  assert(this->IsBoundingBoxNull() == false);
344  BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds();
345 
346  Point3D cornerpoint;
347  cornerpoint[0] = (xFront ? bounds[0] : bounds[1]);
348  cornerpoint[1] = (yFront ? bounds[2] : bounds[3]);
349  cornerpoint[2] = (zFront ? bounds[4] : bounds[5]);
350  if (m_ImageGeometry)
351  {
352  // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the
353  // bounding box. The bounding box itself is no image, so it is corner-based
354  FillVector3D(cornerpoint, cornerpoint[0] - 0.5, cornerpoint[1] - 0.5, cornerpoint[2] - 0.5);
355  }
356 
357  return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint);
358 }
359 
361 {
362  return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction).magnitude() *
363  GetExtent(direction);
364 }
365 
366 void mitk::BaseGeometry::SetExtentInMM(int direction, ScalarType extentInMM)
367 {
368  mitk::ModifiedLock lock(this);
369 
370  ScalarType len = GetExtentInMM(direction);
371  if (fabs(len - extentInMM) >= mitk::eps)
372  {
373  AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix;
374  vnlmatrix = m_GeometryTransform->GetVnlMatrix();
375  if (len > extentInMM)
376  vnlmatrix.set_column(direction, vnlmatrix.get_column(direction) / len * extentInMM);
377  else
378  vnlmatrix.set_column(direction, vnlmatrix.get_column(direction) * extentInMM / len);
379  Matrix3D matrix;
380  matrix = vnlmatrix;
381  m_GeometryTransform->SetMatrix(matrix);
382 
383  Modified();
384  }
385 }
386 
388 {
389  mitk::Point3D index;
390  WorldToIndex(p, index);
391  return IsIndexInside(index);
392 }
393 
395 {
396  bool inside = false;
397  // if it is an image geometry, we need to convert the index to discrete values
398  // this is done by applying the rounding function also used in WorldToIndex (see line 323)
399  if (m_ImageGeometry)
400  {
401  mitk::Point3D discretIndex;
402  discretIndex[0] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[0]);
403  discretIndex[1] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[1]);
404  discretIndex[2] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[2]);
405 
406  inside = this->GetBoundingBox()->IsInside(discretIndex);
407  // we have to check if the index is at the upper border of each dimension,
408  // because the boundingbox is not centerbased
409  if (inside)
410  {
411  const BoundingBox::BoundsArrayType &bounds = this->GetBoundingBox()->GetBounds();
412  if ((discretIndex[0] == bounds[1]) || (discretIndex[1] == bounds[3]) || (discretIndex[2] == bounds[5]))
413  inside = false;
414  }
415  }
416  else
417  inside = this->GetBoundingBox()->IsInside(index);
418 
419  return inside;
420 }
421 
423 {
424  mitk::Vector3D tempIn, tempOut;
425  const TransformType::OffsetType &offset = this->GetIndexToWorldTransform()->GetOffset();
426  tempIn = pt_mm.GetVectorFromOrigin() - offset;
427 
428  WorldToIndex(tempIn, tempOut);
429 
430  pt_units = tempOut;
431 }
432 
434 {
435  // Get WorldToIndex transform
436  if (m_IndexToWorldTransformLastModified != this->GetIndexToWorldTransform()->GetMTime())
437  {
438  if (!m_InvertedTransform)
439  {
440  m_InvertedTransform = TransformType::New();
441  }
442  if (!this->GetIndexToWorldTransform()->GetInverse(m_InvertedTransform.GetPointer()))
443  {
444  itkExceptionMacro("Internal ITK matrix inversion error, cannot proceed.");
445  }
446  m_IndexToWorldTransformLastModified = this->GetIndexToWorldTransform()->GetMTime();
447  }
448 
449  // Check for valid matrix inversion
450  const TransformType::MatrixType &inverse = m_InvertedTransform->GetMatrix();
451  if (inverse.GetVnlMatrix().has_nans())
452  {
453  itkExceptionMacro("Internal ITK matrix inversion error, cannot proceed. Matrix was: "
454  << std::endl
455  << this->GetIndexToWorldTransform()->GetMatrix()
456  << "Suggested inverted matrix is:"
457  << std::endl
458  << inverse);
459  }
460 
461  vec_units = inverse * vec_mm;
462 }
463 
465  const mitk::Vector3D &vec_mm,
466  mitk::Vector3D &vec_units) const
467 {
468  MITK_WARN << "Warning! Call of the deprecated function BaseGeometry::WorldToIndex(point, vec, vec). Use "
469  "BaseGeometry::WorldToIndex(vec, vec) instead!";
470  this->WorldToIndex(vec_mm, vec_units);
471 }
472 
474 {
475  return GetOrigin().GetVnlVector();
476 }
477 
478 vtkLinearTransform *mitk::BaseGeometry::GetVtkTransform() const
479 {
480  return m_GeometryTransform->GetVtkTransform();
481 }
482 
484 {
485  mitk::ModifiedLock lock(this);
486 
487  m_GeometryTransform->SetIdentity();
488  Modified();
489 }
490 
492 {
493  mitk::ModifiedLock lock(this);
494  m_GeometryTransform->Compose(other, pre);
495  Modified();
496 }
497 
498 void mitk::BaseGeometry::Compose(const vtkMatrix4x4 *vtkmatrix, bool pre)
499 {
500  mitk::BaseGeometry::TransformType::Pointer itkTransform = mitk::BaseGeometry::TransformType::New();
501  TransferVtkMatrixToItkTransform(vtkmatrix, itkTransform.GetPointer());
502  Compose(itkTransform, pre);
503 }
504 
506 {
507  if ((vector[0] != 0) || (vector[1] != 0) || (vector[2] != 0))
508  {
509  this->SetOrigin(this->GetOrigin() + vector);
510  }
511 }
512 
514 {
515  pt_mm = this->GetIndexToWorldTransform()->TransformPoint(pt_units);
516 }
517 
519 {
520  vec_mm = this->GetIndexToWorldTransform()->TransformVector(vec_units);
521 }
522 
524 {
525  mitk::ModifiedLock lock(this);
526 
527  vtkTransform *vtktransform = vtkTransform::New();
528  vtktransform->SetMatrix(this->GetVtkMatrix());
529  switch (operation->GetOperationType())
530  {
531  case OpNOTHING:
532  break;
533  case OpMOVE:
534  {
535  auto *pointOp = dynamic_cast<mitk::PointOperation *>(operation);
536  if (pointOp == nullptr)
537  {
538  MITK_ERROR << "Point move operation is null!";
539  return;
540  }
541  mitk::Point3D newPos = pointOp->GetPoint();
542  ScalarType data[3];
543  vtktransform->GetPosition(data);
544  vtktransform->PostMultiply();
545  vtktransform->Translate(newPos[0], newPos[1], newPos[2]);
546  vtktransform->PreMultiply();
547  break;
548  }
549  case OpSCALE:
550  {
551  auto *scaleOp = dynamic_cast<mitk::ScaleOperation *>(operation);
552  if (scaleOp == nullptr)
553  {
554  MITK_ERROR << "Scale operation is null!";
555  return;
556  }
557  mitk::Point3D newScale = scaleOp->GetScaleFactor();
558  ScalarType scalefactor[3];
559 
560  scalefactor[0] = 1 + (newScale[0] / GetMatrixColumn(0).magnitude());
561  scalefactor[1] = 1 + (newScale[1] / GetMatrixColumn(1).magnitude());
562  scalefactor[2] = 1 + (newScale[2] / GetMatrixColumn(2).magnitude());
563 
564  mitk::Point3D anchor = scaleOp->GetScaleAnchorPoint();
565 
566  vtktransform->PostMultiply();
567  vtktransform->Translate(-anchor[0], -anchor[1], -anchor[2]);
568  vtktransform->Scale(scalefactor[0], scalefactor[1], scalefactor[2]);
569  vtktransform->Translate(anchor[0], anchor[1], anchor[2]);
570  break;
571  }
572  case OpROTATE:
573  {
574  auto *rotateOp = dynamic_cast<mitk::RotationOperation *>(operation);
575  if (rotateOp == nullptr)
576  {
577  MITK_ERROR << "Rotation operation is null!";
578  return;
579  }
580  Vector3D rotationVector = rotateOp->GetVectorOfRotation();
581  Point3D center = rotateOp->GetCenterOfRotation();
582  ScalarType angle = rotateOp->GetAngleOfRotation();
583  vtktransform->PostMultiply();
584  vtktransform->Translate(-center[0], -center[1], -center[2]);
585  vtktransform->RotateWXYZ(angle, rotationVector[0], rotationVector[1], rotationVector[2]);
586  vtktransform->Translate(center[0], center[1], center[2]);
587  vtktransform->PreMultiply();
588  break;
589  }
591  {
592  // Copy necessary to avoid vtk warning
593  vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
595  dynamic_cast<mitk::RestorePlanePositionOperation *>(operation)->GetTransform().GetPointer(), matrix);
596  vtktransform->SetMatrix(matrix);
597  matrix->Delete();
598  break;
599  }
601  {
602  auto *applyMatrixOp = dynamic_cast<ApplyTransformMatrixOperation *>(operation);
603  vtktransform->SetMatrix(applyMatrixOp->GetMatrix());
604  break;
605  }
606  default:
607  vtktransform->Delete();
608  return;
609  }
610  this->SetVtkMatrixDeepCopy(vtktransform);
611  Modified();
612  vtktransform->Delete();
613 }
614 
616 {
617  return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction);
618 }
619 
621  const mitk::AffineTransform3D *transform) const
622 {
623  mitk::BoundingBox::PointsContainer::Pointer pointscontainer = mitk::BoundingBox::PointsContainer::New();
624 
625  mitk::BoundingBox::PointIdentifier pointid = 0;
626 
627  unsigned char i;
628  if (transform != nullptr)
629  {
630  mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New();
631  transform->GetInverse(inverse);
632  for (i = 0; i < 8; ++i)
633  pointscontainer->InsertElement(pointid++, inverse->TransformPoint(GetCornerPoint(i)));
634  }
635  else
636  {
637  for (i = 0; i < 8; ++i)
638  pointscontainer->InsertElement(pointid++, GetCornerPoint(i));
639  }
640 
641  mitk::BoundingBox::Pointer result = mitk::BoundingBox::New();
642  result->SetPoints(pointscontainer);
643  result->ComputeBoundingBox();
644 
645  return result;
646 }
647 
649 {
650  std::ostringstream out;
651 
652  out << '[';
653 
654  for (int i = 0; i < 3; ++i)
655  {
656  out << '[';
657  for (int j = 0; j < 3; ++j)
658  out << transformType->GetMatrix().GetVnlMatrix().get(i, j) << ' ';
659  out << ']';
660  }
661 
662  out << "][";
663 
664  for (int i = 0; i < 3; ++i)
665  out << transformType->GetOffset()[i] << ' ';
666 
667  out << "]\0";
668 
669  return out.str();
670 }
671 
673 {
674  m_GeometryTransform->SetIndexToWorldTransformByVtkMatrix(vtkmatrix);
675 }
676 
678 {
679  m_GeometryTransform->SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkmatrix);
680 }
681 
682 void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D & /*atPt3d_units*/,
683  const mitk::Vector3D &vec_units,
684  mitk::Vector3D &vec_mm) const
685 {
686  MITK_WARN << "Warning! Call of the deprecated function BaseGeometry::IndexToWorld(point, vec, vec). Use "
687  "BaseGeometry::IndexToWorld(vec, vec) instead!";
688  // vec_mm = m_IndexToWorldTransform->TransformVector(vec_units);
689  this->IndexToWorld(vec_units, vec_mm);
690 }
691 
693 {
694  return m_GeometryTransform->GetVtkMatrix();
695 }
696 
698 {
699  return m_BoundingBox.IsNull();
700 }
701 
703 {
704  return m_GeometryTransform->IsIndexToWorldTransformNull();
705 }
706 
708 {
709  // If Geometry is switched to ImageGeometry, you have to put an offset to the origin, because
710  // imageGeometries origins are pixel-center-based
711  // ... and remove the offset, if you switch an imageGeometry back to a normal geometry
712  // For more information please see the Geometry documentation page
713 
714  if (m_ImageGeometry == isAnImageGeometry)
715  return;
716 
717  const BoundingBox::BoundsArrayType &boundsarray = this->GetBoundingBox()->GetBounds();
718 
719  Point3D originIndex;
720  FillVector3D(originIndex, boundsarray[0], boundsarray[2], boundsarray[4]);
721 
722  if (isAnImageGeometry == true)
723  FillVector3D(originIndex, originIndex[0] + 0.5, originIndex[1] + 0.5, originIndex[2] + 0.5);
724  else
725  FillVector3D(originIndex, originIndex[0] - 0.5, originIndex[1] - 0.5, originIndex[2] - 0.5);
726 
727  Point3D originWorld;
728 
729  originWorld = GetIndexToWorldTransform()->TransformPoint(originIndex);
730  // instead could as well call IndexToWorld(originIndex,originWorld);
731 
732  SetOrigin(originWorld);
733 
734  this->SetImageGeometry(isAnImageGeometry);
735 }
736 
737 void mitk::BaseGeometry::PrintSelf(std::ostream &os, itk::Indent indent) const
738 {
739  os << indent << " IndexToWorldTransform: ";
740  if (this->IsIndexToWorldTransformNull())
741  os << "nullptr" << std::endl;
742  else
743  {
744  // from itk::MatrixOffsetTransformBase
745  unsigned int i, j;
746  os << std::endl;
747  os << indent << "Matrix: " << std::endl;
748  for (i = 0; i < 3; i++)
749  {
750  os << indent.GetNextIndent();
751  for (j = 0; j < 3; j++)
752  {
753  os << this->GetIndexToWorldTransform()->GetMatrix()[i][j] << " ";
754  }
755  os << std::endl;
756  }
757 
758  os << indent << "Offset: " << this->GetIndexToWorldTransform()->GetOffset() << std::endl;
759  os << indent << "Center: " << this->GetIndexToWorldTransform()->GetCenter() << std::endl;
760  os << indent << "Translation: " << this->GetIndexToWorldTransform()->GetTranslation() << std::endl;
761 
762  os << indent << "Inverse: " << std::endl;
763  for (i = 0; i < 3; i++)
764  {
765  os << indent.GetNextIndent();
766  for (j = 0; j < 3; j++)
767  {
768  os << this->GetIndexToWorldTransform()->GetInverseMatrix()[i][j] << " ";
769  }
770  os << std::endl;
771  }
772 
773  // from itk::ScalableAffineTransform
774  os << indent << "Scale : ";
775  for (i = 0; i < 3; i++)
776  {
777  os << this->GetIndexToWorldTransform()->GetScale()[i] << " ";
778  }
779  os << std::endl;
780  }
781 
782  os << indent << " BoundingBox: ";
783  if (this->IsBoundingBoxNull())
784  os << "nullptr" << std::endl;
785  else
786  {
787  os << indent << "( ";
788  for (unsigned int i = 0; i < 3; i++)
789  {
790  os << this->GetBoundingBox()->GetBounds()[2 * i] << "," << this->GetBoundingBox()->GetBounds()[2 * i + 1] << " ";
791  }
792  os << " )" << std::endl;
793  }
794 
795  os << indent << " Origin: " << this->GetOrigin() << std::endl;
796  os << indent << " ImageGeometry: " << this->GetImageGeometry() << std::endl;
797  os << indent << " Spacing: " << this->GetSpacing() << std::endl;
798 }
799 
801 {
802  if (!m_ModifiedLockFlag)
803  Superclass::Modified();
804  else
805  m_ModifiedCalledFlag = true;
806 }
807 
809 {
810  return m_GeometryTransform->GetIndexToWorldTransform();
811 }
812 
814 {
815  return m_GeometryTransform->GetIndexToWorldTransform();
816 }
817 
819 {
820  return m_GeometryTransform;
821 }
822 
824  const mitk::BaseGeometry::BoundingBoxType *rightHandSide,
825  ScalarType eps,
826  bool verbose)
827 {
828  if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
829  {
830  MITK_ERROR << "mitk::Equal( const mitk::Geometry3D::BoundingBoxType *leftHandSide, const "
831  "mitk::Geometry3D::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) does not with nullptr "
832  "pointer input.";
833  return false;
834  }
835  return Equal(*leftHandSide, *rightHandSide, eps, verbose);
836 }
837 
839  const mitk::BaseGeometry::BoundingBoxType &rightHandSide,
840  ScalarType eps,
841  bool verbose)
842 {
843  bool result = true;
844 
845  BaseGeometry::BoundsArrayType rightBounds = rightHandSide.GetBounds();
846  BaseGeometry::BoundsArrayType leftBounds = leftHandSide.GetBounds();
847  BaseGeometry::BoundsArrayType::Iterator itLeft = leftBounds.Begin();
848  for (BaseGeometry::BoundsArrayType::Iterator itRight = rightBounds.Begin(); itRight != rightBounds.End(); ++itRight)
849  {
850  if ((!mitk::Equal(*itLeft, *itRight, eps)))
851  {
852  if (verbose)
853  {
854  MITK_INFO << "[( Geometry3D::BoundingBoxType )] bounds are not equal.";
855  MITK_INFO << "rightHandSide is " << setprecision(12) << *itRight << " : leftHandSide is " << *itLeft
856  << " and tolerance is " << eps;
857  }
858  result = false;
859  }
860  itLeft++;
861  }
862  return result;
863 }
864 
865 bool mitk::Equal(const mitk::BaseGeometry *leftHandSide,
866  const mitk::BaseGeometry *rightHandSide,
867  ScalarType eps,
868  bool verbose)
869 {
870  if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
871  {
872  MITK_ERROR << "mitk::Equal(const mitk::Geometry3D *leftHandSide, const mitk::Geometry3D *rightHandSide, ScalarType "
873  "eps, bool verbose) does not with nullptr pointer input.";
874  return false;
875  }
876  return Equal(*leftHandSide, *rightHandSide, eps, verbose);
877 }
878 
879 bool mitk::Equal(const mitk::BaseGeometry &leftHandSide,
880  const mitk::BaseGeometry &rightHandSide,
881  ScalarType eps,
882  bool verbose)
883 {
884  bool result = true;
885 
886  // Compare spacings
887  if (!mitk::Equal(leftHandSide.GetSpacing(), rightHandSide.GetSpacing(), eps))
888  {
889  if (verbose)
890  {
891  MITK_INFO << "[( Geometry3D )] Spacing differs.";
892  MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetSpacing() << " : leftHandSide is "
893  << leftHandSide.GetSpacing() << " and tolerance is " << eps;
894  }
895  result = false;
896  }
897 
898  // Compare Origins
899  if (!mitk::Equal(leftHandSide.GetOrigin(), rightHandSide.GetOrigin(), eps))
900  {
901  if (verbose)
902  {
903  MITK_INFO << "[( Geometry3D )] Origin differs.";
904  MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetOrigin() << " : leftHandSide is "
905  << leftHandSide.GetOrigin() << " and tolerance is " << eps;
906  }
907  result = false;
908  }
909 
910  // Compare Axis and Extents
911  for (unsigned int i = 0; i < 3; ++i)
912  {
913  if (!mitk::Equal(leftHandSide.GetAxisVector(i), rightHandSide.GetAxisVector(i), eps))
914  {
915  if (verbose)
916  {
917  MITK_INFO << "[( Geometry3D )] AxisVector #" << i << " differ";
918  MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetAxisVector(i) << " : leftHandSide is "
919  << leftHandSide.GetAxisVector(i) << " and tolerance is " << eps;
920  }
921  result = false;
922  }
923 
924  if (!mitk::Equal(leftHandSide.GetExtent(i), rightHandSide.GetExtent(i), eps))
925  {
926  if (verbose)
927  {
928  MITK_INFO << "[( Geometry3D )] Extent #" << i << " differ";
929  MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetExtent(i) << " : leftHandSide is "
930  << leftHandSide.GetExtent(i) << " and tolerance is " << eps;
931  }
932  result = false;
933  }
934  }
935 
936  // Compare ImageGeometry Flag
937  if (rightHandSide.GetImageGeometry() != leftHandSide.GetImageGeometry())
938  {
939  if (verbose)
940  {
941  MITK_INFO << "[( Geometry3D )] GetImageGeometry is different.";
942  MITK_INFO << "rightHandSide is " << rightHandSide.GetImageGeometry() << " : leftHandSide is "
943  << leftHandSide.GetImageGeometry();
944  }
945  result = false;
946  }
947 
948  // Compare FrameOfReference ID
949  if (rightHandSide.GetFrameOfReferenceID() != leftHandSide.GetFrameOfReferenceID())
950  {
951  if (verbose)
952  {
953  MITK_INFO << "[( Geometry3D )] GetFrameOfReferenceID is different.";
954  MITK_INFO << "rightHandSide is " << rightHandSide.GetFrameOfReferenceID() << " : leftHandSide is "
955  << leftHandSide.GetFrameOfReferenceID();
956  }
957  result = false;
958  }
959 
960  // Compare BoundingBoxes
961  if (!mitk::Equal(*leftHandSide.GetBoundingBox(), *rightHandSide.GetBoundingBox(), eps, verbose))
962  {
963  result = false;
964  }
965 
966  // Compare IndexToWorldTransform Matrix
967  if (!mitk::Equal(*leftHandSide.GetIndexToWorldTransform(), *rightHandSide.GetIndexToWorldTransform(), eps, verbose))
968  {
969  result = false;
970  }
971  return result;
972 }
973 
975  const mitk::BaseGeometry::TransformType *rightHandSide,
976  ScalarType eps,
977  bool verbose)
978 {
979  if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
980  {
981  MITK_ERROR << "mitk::Equal(const Geometry3D::TransformType *leftHandSide, const Geometry3D::TransformType "
982  "*rightHandSide, ScalarType eps, bool verbose ) does not with nullptr pointer input.";
983  return false;
984  }
985  return Equal(*leftHandSide, *rightHandSide, eps, verbose);
986 }
987 
989  const mitk::BaseGeometry::TransformType &rightHandSide,
990  ScalarType eps,
991  bool verbose)
992 {
993  // Compare IndexToWorldTransform Matrix
994  if (!mitk::MatrixEqualElementWise(leftHandSide.GetMatrix(), rightHandSide.GetMatrix(), eps))
995  {
996  if (verbose)
997  {
998  MITK_INFO << "[( Geometry3D::TransformType )] Index to World Transformation matrix differs.";
999  MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetMatrix() << " : leftHandSide is "
1000  << leftHandSide.GetMatrix() << " and tolerance is " << eps;
1001  }
1002  return false;
1003  }
1004  return true;
1005 }
1006 
1008  const mitk::BaseGeometry& referenceGeo,
1009  ScalarType eps,
1010  bool verbose)
1011 {
1012  bool result = true;
1013 
1014  // Compare spacings (must be equal)
1015  const auto testedSpacing = testGeo.GetSpacing();
1016  if (!mitk::Equal(testedSpacing, referenceGeo.GetSpacing(), eps))
1017  {
1018  if (verbose)
1019  {
1020  MITK_INFO << "[( Geometry3D )] Spacing differs.";
1021  MITK_INFO << "testedGeometry is " << setprecision(12) << testedSpacing << " : referenceGeometry is "
1022  << referenceGeo.GetSpacing() << " and tolerance is " << eps;
1023  }
1024  result = false;
1025  }
1026 
1027  // Compare ImageGeometry Flag (must be equal)
1028  if (referenceGeo.GetImageGeometry() != testGeo.GetImageGeometry())
1029  {
1030  if (verbose)
1031  {
1032  MITK_INFO << "[( Geometry3D )] GetImageGeometry is different.";
1033  MITK_INFO << "referenceGeo is " << referenceGeo.GetImageGeometry() << " : testGeo is "
1034  << testGeo.GetImageGeometry();
1035  }
1036  result = false;
1037  }
1038 
1039  // Compare IndexToWorldTransform Matrix (must be equal -> same axis directions)
1040  if (!Equal(*(testGeo.GetIndexToWorldTransform()), *(referenceGeo.GetIndexToWorldTransform()), eps, verbose))
1041  {
1042  result = false;
1043  }
1044 
1045  // Compare BoundingBoxes (must be <=, thus corners of the tests geom must be inside the reference)
1046  for (int i = 0; i<8; ++i)
1047  {
1048  auto testCorner = testGeo.GetCornerPoint(i);
1049  if (testGeo.GetImageGeometry())
1050  {
1051  //if it is an image geometry the upper bounds indicate the first index that is not inside
1052  //the geometry. Thus we will get the lower corener of the "first" voxel outside,
1053  //but we need the corner of the "last" voxel inside the test geometry.
1054  //To achive that we supstract one spacing from each index elment that is constructed
1055  //by an upper bound value (see implementation of BaseGeometry::GetCorner().
1056  std::bitset<sizeof(int)> bs(i);
1057  if (bs.test(0))
1058  {
1059  testCorner[2] -= testedSpacing[2];
1060  }
1061  if (bs.test(1))
1062  {
1063  testCorner[1] -= testedSpacing[1];
1064  }
1065  if (bs.test(2))
1066  {
1067  testCorner[0] -= testedSpacing[0];
1068  }
1069  }
1070 
1071  if (!referenceGeo.IsInside(testCorner))
1072  {
1073  if (verbose)
1074  {
1075  MITK_INFO << "[( Geometry3D )] corner point is not inside. ";
1076  MITK_INFO << "referenceGeo is " << setprecision(12) << referenceGeo << " : tested corner is "
1077  << testGeo.GetCornerPoint(i);
1078  }
1079  result = false;
1080  }
1081  }
1082 
1083  // check grid of test geometry is on the grid of the reference geometry. This is important as the
1084  // boundingbox is only checked for containing the tested geometry, but if a corner (one is enough
1085  // as we know that axis and spacing are equal) of the tested geometry is on the grid it is really
1086  // a sub geometry (as they have the same spacing and axis).
1087  auto cornerOffset = testGeo.GetCornerPoint(0) - referenceGeo.GetCornerPoint(0);
1088  for (unsigned int i = 0; i < 3; ++i)
1089  {
1090  auto pixelCountContinous = cornerOffset[i] / testedSpacing[i];
1091  auto pixelCount = std::round(pixelCountContinous);
1092  if (std::abs(pixelCount - pixelCountContinous) > eps)
1093  {
1094  if (verbose)
1095  {
1096  MITK_INFO << "[( Geometry3D )] Tested geometry is not on the grid of the reference geometry. ";
1097  MITK_INFO << "referenceGeo is " << setprecision(15) << referenceGeo << " : tested corner offset in pixels is "
1098  << pixelCountContinous << " for axis "<<i;
1099  }
1100  result = false;
1101  }
1102  }
1103 
1104  return result;
1105 }
void TransferVtkMatrixToItkTransform(const vtkMatrix4x4 *vtkmatrix, TTransformType *itkTransform)
void SetVtkMatrixDeepCopy(vtkTransform *vtktransform)
itk::BoundingBox< unsigned long, 3, ScalarType > BoundingBoxType
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?
vtkMatrix4x4 * GetVtkMatrix()
void SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkMatrix4x4 *vtkmatrix)
Convenience method for setting the ITK transform (m_IndexToWorldTransform) via an vtkMatrix4x4...
void SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
Set the spacing (m_Spacing).
void SetIndexToWorldTransformWithoutChangingSpacing(mitk::AffineTransform3D *transform)
MITKCORE_EXPORT bool IsSubGeometry(const mitk::BaseGeometry &testGeo, const mitk::BaseGeometry &referenceGeo, ScalarType eps, bool verbose)
A function checks if a test geometry is a sub geometry of a given reference geometry.
void SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
bool MatrixEqualElementWise(const vnl_matrix_fixed< TCoordRep, NRows, NCols > &matrix1, const vnl_matrix_fixed< TCoordRep, NRows, NCols > &matrix2, mitk::ScalarType epsilon=mitk::eps)
Check for element-wise matrix equality with a user defined accuracy.
Definition: mitkMatrix.h:140
bool IsBoundingBoxNull() const
BoundingBoxType::BoundsArrayType BoundsArrayType
#define MITK_INFO
Definition: mitkLogMacros.h:18
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
vnl_vector< ScalarType > VnlVector
Definition: mitkVector.h:134
#define MITK_ERROR
Definition: mitkLogMacros.h:20
bool IsIndexInside(const mitk::Point3D &index) const
Test whether the point p ((continous!)index coordinates in units) is inside the bounding box...
double ScalarType
void SetIdentity()
Set the transform to identity, the spacing to 1 and origin to 0.
vtkLinearTransform * GetVtkTransform() const
Get the m_IndexToWorldTransform as a vtkLinearTransform.
bool IsInside(const mitk::Point3D &p) const
Test whether the point p (world coordinates in mm) is inside the bounding box.
void _SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
double GetDiagonalLength2() const
Get the squared length of the diagonal of the bounding-box in mm.
GeometryTransformHolder::TransformType TransformType
ModifiedLock manages the calls of Modified() functions.
ScalarType GetExtent(unsigned int direction) const
Set the time bounds (in ms)
DataCollection - Class to facilitate loading/accessing structured data.
Point3D GetCornerPoint(int id) const
Get the position of the corner number id (in world coordinates)
void TransferItkTransformToVtkMatrix(const TTransformType *itkTransform, vtkMatrix4x4 *vtkmatrix)
AffineTransform3D::MatrixType::InternalMatrixType GetVnlMatrix()
void Initialize()
Initialize the BaseGeometry.
Constants for most interaction classes, due to the generic StateMachines.
void SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
Set the spacing (m_Spacing).
void SetIndexToWorldTransformWithoutChangingSpacing(mitk::AffineTransform3D *transform)
Point3D GetCenter() const
Get the center of the bounding-box in mm.
void Translate(const Vector3D &vector)
Translate the origin by a vector.
virtual void CheckIndexToWorldTransform(mitk::AffineTransform3D *)
CheckIndexToWorldTransform.
abstract class, that can be used by Undo to undo an operation.
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 GeometryTransformHolder * GetGeometryTransformHolder() const
void Compose(const TransformType *other, bool pre=false)
Compose new IndexToWorldTransform with a given transform.
virtual void CheckBounds(const BoundsArrayType &)
CheckBounds.
void InitializeGeometry(Self *newGeometry) const
virtual void PreSetSpacing(const mitk::Vector3D &)
PreSetSpacing.
static Vector3D offset
static const std::string GetTransformAsString(TransformType *transformType)
void SetVtkMatrixDeepCopy(vtkTransform *vtktransform)
#define MITK_WARN
Definition: mitkLogMacros.h:19
void Compose(const TransformType *other, bool pre=false)
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
#define mitkThrow()
bool verbose(false)
void SetOrigin(const Point3D &origin)
Set the origin, i.e. the upper-left corner of the plane.
vtkLinearTransform * GetVtkTransform() const
Get the m_IndexToWorldTransform as a vtkLinearTransform.
void Modified() const override
Overload of function Modified() to prohibit several calls of Modified() using the ModifiedLock class...
Operation that handles all actions on one Point.
bool IsIndexToWorldTransformNull() const
itk::AffineGeometryFrame< ScalarType, 3 >::TransformType AffineTransform3D
virtual unsigned int GetFrameOfReferenceID() const
Get the DICOM FrameOfReferenceID referring to the used world coordinate system.
void SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
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.
virtual bool IsValid() const
Is this BaseGeometry in a state that is valid?
virtual void ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry)
When switching from an Image Geometry to a normal Geometry (and the other way around), you have to.
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)
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
The ScaleOperation is an operation to scale any mitk::BaseGeometry.
void SetExtentInMM(int direction, ScalarType extentInMM)
Set the extent of the bounding-box in the specified direction in mm.
double GetDiagonalLength() const
Get the length of the diagonal of the bounding-box in mm.
void SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4 *vtkmatrix)
Convenience method for setting the ITK transform (m_IndexToWorldTransform) via an vtkMatrix4x4...
VnlVector GetOriginVnl() const
Get the origin as VnlVector.
virtual void SetFrameOfReferenceID(unsigned int _arg)
Set the DICOM FrameOfReferenceID referring to the used world coordinate system.
mitk::BoundingBox::Pointer CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D *transform) const
Calculates a bounding-box around the geometry relative to a coordinate system defined by a transform...
MITKCORE_EXPORT const ScalarType eps
void SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkMatrix4x4 *vtkmatrix)
Convenience method for setting the ITK transform (m_IndexToWorldTransform) via an vtkMatrix4x4...
virtual void SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4 *vtkmatrix)
Convenience method for setting the ITK transform (m_IndexToWorldTransform) via an vtkMatrix4x4...
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.
OperationType GetOperationType()
virtual bool Is2DConvertable()
Checks, if the given geometry can be converted to 2D without information loss e.g. when a 2D image is saved, the matrix is usually cropped to 2x2, and when you load it back to MITK it will be filled with standard values. This function checks, if information would be lost during this procedure.
Operation, that holds everything necessary for an rotation operation on mitk::BaseData.
void SetOrigin(const Point3D &origin)
Set the origin, i.e. the upper-left corner of the plane.
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
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.
VnlVector GetMatrixColumn(unsigned int direction) const
Get a VnlVector along bounding-box in the specified direction, length is spacing. ...
BoundingBoxType::BoundsArrayType BoundsArrayType
void SetFloatBounds(const float bounds[6])
Set the bounding box (in index/unit coordinates) via a float array.
virtual const BoundingBoxType * GetBoundingBox()