diff --git a/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp b/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp index db37fc986a..be889b869c 100644 --- a/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp +++ b/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp @@ -1,277 +1,277 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkAbstractTransformGeometry.h" #include mitk::AbstractTransformGeometry::AbstractTransformGeometry() : m_Plane(NULL), m_FrameGeometry(NULL) { Initialize(); } mitk::AbstractTransformGeometry::AbstractTransformGeometry(const AbstractTransformGeometry& other) : Superclass(other) { if(other.m_ParametricBoundingBox.IsNotNull()) { this->SetParametricBounds(m_ParametricBoundingBox->GetBounds()); } this->SetPlane(other.m_Plane); this->SetFrameGeometry(other.m_FrameGeometry); } mitk::AbstractTransformGeometry::~AbstractTransformGeometry() { } void mitk::AbstractTransformGeometry::InternPostInitialize() { m_ItkVtkAbstractTransform = itk::VtkAbstractTransform::New(); } vtkAbstractTransform* mitk::AbstractTransformGeometry::GetVtkAbstractTransform() const { return m_ItkVtkAbstractTransform->GetVtkAbstractTransform(); } mitk::ScalarType mitk::AbstractTransformGeometry::GetParametricExtentInMM(int direction) const { if(m_Plane.IsNull()) { itkExceptionMacro(<<"m_Plane is NULL."); } return m_Plane->GetExtentInMM(direction); } const mitk::Transform3D* mitk::AbstractTransformGeometry::GetParametricTransform() const { return m_ItkVtkAbstractTransform; } bool mitk::AbstractTransformGeometry::Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const { - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); mitk::Point2D pt2d_mm; bool isInside; isInside = Map(pt3d_mm, pt2d_mm); Map(pt2d_mm, projectedPt3d_mm); return isInside; //Point3D pt3d_units; //pt3d_units = m_ItkVtkAbstractTransform->BackTransform(pt3d_mm); //pt3d_units[2] = 0; //projectedPt3d_mm = m_ItkVtkAbstractTransform->TransformPoint(pt3d_units); //return const_cast(m_BoundingBox.GetPointer())->IsInside(pt3d_units); } bool mitk::AbstractTransformGeometry::Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const { assert((m_ItkVtkAbstractTransform.IsNotNull()) && (m_Plane.IsNotNull())); Point3D pt3d_units; pt3d_units = m_ItkVtkAbstractTransform->BackTransform(pt3d_mm); return m_Plane->Map(pt3d_units, pt2d_mm); } void mitk::AbstractTransformGeometry::Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const { assert((m_ItkVtkAbstractTransform.IsNotNull()) && (m_Plane.IsNotNull())); m_Plane->Map(pt2d_mm, pt3d_mm); pt3d_mm = m_ItkVtkAbstractTransform->TransformPoint(pt3d_mm); } bool mitk::AbstractTransformGeometry::Project(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { itkExceptionMacro("not implemented yet - replace GetIndexToWorldTransform by m_ItkVtkAbstractTransform->GetInverseVtkAbstractTransform()"); - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; vec3d_units = GetIndexToWorldTransform()->GetInverseMatrix() * vec3d_mm; vec3d_units[2] = 0; projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units); Point3D pt3d_units; mitk::ScalarType temp[3]; unsigned int i, j; for (j = 0; j < 3; ++j) temp[j] = atPt3d_mm[j] - GetIndexToWorldTransform()->GetOffset()[j]; for (i = 0; i < 3; ++i) { pt3d_units[i] = 0.0; for (j = 0; j < 3; ++j) pt3d_units[i] += GetIndexToWorldTransform()->GetInverseMatrix()[i][j] * temp[j]; } - return const_cast(m_BoundingBox.GetPointer())->IsInside(pt3d_units); + return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool mitk::AbstractTransformGeometry::Project(const mitk::Vector3D &/*vec3d_mm*/, mitk::Vector3D &/*projectedVec3d_mm*/) const { MITK_WARN << "Need additional point! No standard value defined. Please use Project(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm). Unfortunatley this one is not implemented at the moment. Sorry :("; itkExceptionMacro("not implemented yet - replace GetIndexToWorldTransform by m_ItkVtkAbstractTransform->GetInverseVtkAbstractTransform()"); return false; } bool mitk::AbstractTransformGeometry::Map(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector2D &vec2d_mm) const { assert((m_ItkVtkAbstractTransform.IsNotNull()) && (m_Plane.IsNotNull())); ScalarType vtkpt[3], vtkvec[3]; itk2vtk(atPt3d_mm, vtkpt); itk2vtk(vec3d_mm, vtkvec); m_ItkVtkAbstractTransform->GetInverseVtkAbstractTransform()->TransformVectorAtPoint(vtkpt, vtkvec, vtkvec); mitk::Vector3D vec3d_units; vtk2itk(vtkvec, vec3d_units); return m_Plane->Map(atPt3d_mm, vec3d_units, vec2d_mm); } void mitk::AbstractTransformGeometry::Map(const mitk::Point2D & atPt2d_mm, const mitk::Vector2D &vec2d_mm, mitk::Vector3D &vec3d_mm) const { m_Plane->Map(atPt2d_mm, vec2d_mm, vec3d_mm); Point3D atPt3d_mm; Map(atPt2d_mm, atPt3d_mm); float vtkpt[3], vtkvec[3]; itk2vtk(atPt3d_mm, vtkpt); itk2vtk(vec3d_mm, vtkvec); m_ItkVtkAbstractTransform->GetVtkAbstractTransform()->TransformVectorAtPoint(vtkpt, vtkvec, vtkvec); vtk2itk(vtkvec, vec3d_mm); } void mitk::AbstractTransformGeometry::IndexToWorld(const mitk::Point2D &pt_units, mitk::Point2D &pt_mm) const { m_Plane->IndexToWorld(pt_units, pt_mm); } void mitk::AbstractTransformGeometry::WorldToIndex(const mitk::Point2D &pt_mm, mitk::Point2D &pt_units) const { m_Plane->WorldToIndex(pt_mm, pt_units); } void mitk::AbstractTransformGeometry::IndexToWorld(const mitk::Point2D & /*atPt2d_units*/, const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const { MITK_WARN<<"Warning! Call of the deprecated function AbstractTransformGeometry::IndexToWorld(point, vec, vec). Use AbstractTransformGeometry::IndexToWorld(vec, vec) instead!"; this->IndexToWorld(vec_units, vec_mm); } void mitk::AbstractTransformGeometry::IndexToWorld(const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const { m_Plane->IndexToWorld(vec_units, vec_mm); } void mitk::AbstractTransformGeometry::WorldToIndex(const mitk::Point2D & /*atPt2d_mm*/, const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const { MITK_WARN<<"Warning! Call of the deprecated function AbstractTransformGeometry::WorldToIndex(point, vec, vec). Use AbstractTransformGeometry::WorldToIndex(vec, vec) instead!"; this->WorldToIndex(vec_mm, vec_units); } void mitk::AbstractTransformGeometry::WorldToIndex(const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const { m_Plane->WorldToIndex(vec_mm, vec_units); } bool mitk::AbstractTransformGeometry::IsAbove(const mitk::Point3D& pt3d_mm) const { assert((m_ItkVtkAbstractTransform.IsNotNull()) && (m_Plane.IsNotNull())); Point3D pt3d_ParametricWorld; pt3d_ParametricWorld = m_ItkVtkAbstractTransform->BackTransform(pt3d_mm); Point3D pt3d_ParametricUnits; ((Geometry3D*)m_Plane)->WorldToIndex(pt3d_ParametricWorld, pt3d_ParametricUnits); return (pt3d_ParametricUnits[2] > m_ParametricBoundingBox->GetBounds()[4]); } void mitk::AbstractTransformGeometry::SetVtkAbstractTransform(vtkAbstractTransform* aVtkAbstractTransform) { m_ItkVtkAbstractTransform->SetVtkAbstractTransform(aVtkAbstractTransform); } void mitk::AbstractTransformGeometry::SetPlane(const mitk::PlaneGeometry* aPlane) { if(aPlane!=NULL) { m_Plane = static_cast(aPlane->Clone().GetPointer()); BoundingBox::BoundsArrayType b=m_Plane->GetBoundingBox()->GetBounds(); SetParametricBounds(b); CalculateFrameGeometry(); } else { if(m_Plane.IsNull()) return; m_Plane=NULL; } Modified(); } void mitk::AbstractTransformGeometry::CalculateFrameGeometry() { if((m_Plane.IsNull()) || (m_FrameGeometry.IsNotNull())) return; //@warning affine-transforms and bounding-box should be set by specific sub-classes! SetBounds(m_Plane->GetBoundingBox()->GetBounds()); } void mitk::AbstractTransformGeometry::SetFrameGeometry(const mitk::Geometry3D* frameGeometry) { if((frameGeometry != NULL) && (frameGeometry->IsValid())) { m_FrameGeometry = static_cast(frameGeometry->Clone().GetPointer()); SetIndexToWorldTransform(m_FrameGeometry->GetIndexToWorldTransform()); SetBounds(m_FrameGeometry->GetBounds()); } else { m_FrameGeometry = NULL; } } unsigned long mitk::AbstractTransformGeometry::GetMTime() const { if(Superclass::GetMTime()GetMTime()) return m_ItkVtkAbstractTransform->GetMTime(); return Superclass::GetMTime(); } void mitk::AbstractTransformGeometry::SetOversampling(mitk::ScalarType oversampling) { if(m_Plane.IsNull()) { itkExceptionMacro(<< "m_Plane is not set."); } mitk::BoundingBox::BoundsArrayType bounds = m_Plane->GetBounds(); bounds[1]*=oversampling; bounds[3]*=oversampling; bounds[5]*=oversampling; SetParametricBounds(bounds); } itk::LightObject::Pointer mitk::AbstractTransformGeometry::InternalClone() const { Self::Pointer newGeometry = new AbstractTransformGeometry(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } diff --git a/Core/Code/DataManagement/mitkBaseGeometry.cpp b/Core/Code/DataManagement/mitkBaseGeometry.cpp index 66228d6a62..a0ff6a9e1d 100644 --- a/Core/Code/DataManagement/mitkBaseGeometry.cpp +++ b/Core/Code/DataManagement/mitkBaseGeometry.cpp @@ -1,926 +1,930 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include "mitkBaseGeometry.h" #include "mitkvector.h" #include "mitkMatrixConvert.h" #include #include #include "mitkRotationOperation.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkApplyTransformMatrixOperation.h" #include "mitkPointOperation.h" #include "mitkInteractionConst.h" mitk::BaseGeometry::BaseGeometry(): Superclass(), mitk::OperationActor(), - m_Valid(true), m_FrameOfReferenceID(0), m_IndexToWorldTransformLastModified(0) + m_FrameOfReferenceID(0), m_IndexToWorldTransformLastModified(0) { FillVector3D(m_FloatSpacing, 1,1,1); m_VtkMatrix = vtkMatrix4x4::New(); m_VtkIndexToWorldTransform = vtkMatrixToLinearTransform::New(); m_VtkIndexToWorldTransform->SetInput(m_VtkMatrix); Initialize(); } mitk::BaseGeometry::BaseGeometry(const BaseGeometry& other): Superclass(), m_TimeBounds(other.m_TimeBounds), - m_Valid(other.m_Valid), m_FrameOfReferenceID(other.m_FrameOfReferenceID), m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified), m_Origin(other.m_Origin) + m_FrameOfReferenceID(other.m_FrameOfReferenceID), m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified), m_Origin(other.m_Origin) { // DEPRECATED(m_RotationQuaternion = other.m_RotationQuaternion); // AffineGeometryFrame SetBounds(other.GetBounds()); m_VtkMatrix = vtkMatrix4x4::New(); m_VtkMatrix->DeepCopy(other.m_VtkMatrix); FillVector3D(m_FloatSpacing,other.m_FloatSpacing[0],other.m_FloatSpacing[1],other.m_FloatSpacing[2]); m_VtkIndexToWorldTransform = vtkMatrixToLinearTransform::New(); m_VtkIndexToWorldTransform->DeepCopy(other.m_VtkIndexToWorldTransform); m_VtkIndexToWorldTransform->SetInput(m_VtkMatrix); other.InitializeGeometry(this); } mitk::BaseGeometry::~BaseGeometry() { m_VtkMatrix->Delete(); m_VtkIndexToWorldTransform->Delete(); } const mitk::Point3D& mitk::BaseGeometry::GetOrigin() const { return m_Origin; } void mitk::BaseGeometry::SetOrigin(const Point3D & origin) { if(origin!=GetOrigin()) { m_Origin = origin; m_IndexToWorldTransform->SetOffset(m_Origin.GetVectorFromOrigin()); Modified(); TransferItkToVtkTransform(); } } void mitk::BaseGeometry::TransferItkToVtkTransform() { // copy m_IndexToWorldTransform into m_VtkIndexToWorldTransform TransferItkTransformToVtkMatrix(m_IndexToWorldTransform.GetPointer(), m_VtkMatrix); m_VtkIndexToWorldTransform->Modified(); } void mitk::BaseGeometry::CopySpacingFromTransform(mitk::AffineTransform3D* transform, mitk::Vector3D& spacing, float floatSpacing[3]) { mitk::AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = transform->GetMatrix().GetVnlMatrix(); spacing[0]=vnlmatrix.get_column(0).magnitude(); spacing[1]=vnlmatrix.get_column(1).magnitude(); spacing[2]=vnlmatrix.get_column(2).magnitude(); floatSpacing[0]=spacing[0]; floatSpacing[1]=spacing[1]; floatSpacing[2]=spacing[2]; } void mitk::BaseGeometry::Initialize() { float b[6] = {0,1,0,1,0,1}; SetFloatBounds(b); if(m_IndexToWorldTransform.IsNull()) m_IndexToWorldTransform = TransformType::New(); else m_IndexToWorldTransform->SetIdentity(); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); m_VtkMatrix->Identity(); m_TimeBounds[0]=ScalarTypeNumericTraits::NonpositiveMin(); m_TimeBounds[1]=ScalarTypeNumericTraits::max(); m_FrameOfReferenceID = 0; this->InternPostInitialize(); } void mitk::BaseGeometry::SetFloatBounds(const float bounds[6]) { mitk::BoundingBox::BoundsArrayType b; const float *input = bounds; int i=0; for(mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6 ;++i) *it++ = (mitk::ScalarType)*input++; SetBounds(b); } void mitk::BaseGeometry::SetFloatBounds(const double bounds[6]) { mitk::BoundingBox::BoundsArrayType b; const double *input = bounds; int i=0; for(mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6 ;++i) *it++ = (mitk::ScalarType)*input++; SetBounds(b); } /** Initialize the geometry */ void mitk::BaseGeometry::InitializeGeometry(BaseGeometry* newGeometry) const { newGeometry->SetBounds(m_BoundingBox->GetBounds()); // we have to create a new transform!! newGeometry->SetTimeBounds(m_TimeBounds); newGeometry->SetFrameOfReferenceID(GetFrameOfReferenceID()); if(m_IndexToWorldTransform) { TransformType::Pointer indexToWorldTransform = TransformType::New(); indexToWorldTransform->SetCenter( m_IndexToWorldTransform->GetCenter() ); indexToWorldTransform->SetMatrix( m_IndexToWorldTransform->GetMatrix() ); indexToWorldTransform->SetOffset( m_IndexToWorldTransform->GetOffset() ); newGeometry->SetIndexToWorldTransform(indexToWorldTransform); } this->InternPostInitializeGeometry(newGeometry); } /** Set the bounds */ void mitk::BaseGeometry::SetBounds(const BoundsArrayType& bounds) { InternPreSetBounds(bounds); m_BoundingBox = BoundingBoxType::New(); BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New(); BoundingBoxType::PointType p; BoundingBoxType::PointIdentifier pointid; for(pointid=0; pointid<2;++pointid) { unsigned int i; - for(i=0; iInsertElement(pointid, p); } m_BoundingBox->SetPoints(pointscontainer); m_BoundingBox->ComputeBoundingBox(); this->Modified(); } void mitk::BaseGeometry::InternPreSetBounds(const BoundsArrayType& bounds){}; void mitk::BaseGeometry::SetIndexToWorldTransform(mitk::AffineTransform3D* transform) { InternPreSetIndexToWorldTransform(transform); if(m_IndexToWorldTransform.GetPointer() != transform) { m_IndexToWorldTransform = transform; CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); TransferItkToVtkTransform(); Modified(); } InternPostSetIndexToWorldTransform(transform); } void mitk::BaseGeometry::InternPreSetIndexToWorldTransform(mitk::AffineTransform3D* transform) {} void mitk::BaseGeometry::InternPostSetIndexToWorldTransform(mitk::AffineTransform3D* transform) {} const mitk::BaseGeometry::BoundsArrayType mitk::BaseGeometry::GetBounds() const { assert(m_BoundingBox.IsNotNull()); return m_BoundingBox->GetBounds(); } bool mitk::BaseGeometry::IsValid() const -{ - bool isValid = m_Valid; - isValid = isValid && this->InternPostIsValid(); - - return isValid; -} - -bool mitk::BaseGeometry::InternPostIsValid() const { return true; } const float* mitk::BaseGeometry::GetFloatSpacing() const { return m_FloatSpacing; } bool mitk::Equal(const BaseGeometry::TransformType *leftHandSide, const BaseGeometry::TransformType *rightHandSide, ScalarType eps, bool verbose ) { //Compare IndexToWorldTransform Matrix if( !mitk::MatrixEqualElementWise( leftHandSide->GetMatrix(), rightHandSide->GetMatrix() ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry::TransformType )] Index to World Transformation matrix differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetMatrix() << " : leftHandSide is " << leftHandSide->GetMatrix() << " and tolerance is " << eps; } return false; } return true; } bool mitk::Equal( const mitk::BaseGeometry::BoundingBoxType *leftHandSide, const mitk::BaseGeometry::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) { bool result = true; if( rightHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry::BoundingBoxType )] rightHandSide NULL."; return false; } if( leftHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry::BoundingBoxType )] leftHandSide NULL."; return false; } BaseGeometry::BoundsArrayType rightBounds = rightHandSide->GetBounds(); BaseGeometry::BoundsArrayType leftBounds = leftHandSide->GetBounds(); BaseGeometry::BoundsArrayType::Iterator itLeft = leftBounds.Begin(); for( BaseGeometry::BoundsArrayType::Iterator itRight = rightBounds.Begin(); itRight != rightBounds.End(); ++itRight) { if(( !mitk::Equal( *itLeft, *itRight, eps )) ) { if(verbose) { MITK_INFO << "[( BaseGeometry::BoundingBoxType )] bounds are not equal."; MITK_INFO << "rightHandSide is " << setprecision(12) << *itRight << " : leftHandSide is " << *itLeft << " and tolerance is " << eps; } result = false; } itLeft++; } return result; } bool mitk::Equal(const mitk::BaseGeometry *leftHandSide, const mitk::BaseGeometry *rightHandSide, ScalarType eps, bool verbose) { bool result = true; if( rightHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry )] rightHandSide NULL."; return false; } if( leftHandSide == NULL) { if(verbose) MITK_INFO << "[( BaseGeometry )] leftHandSide NULL."; return false; } //Compare spacings if( !mitk::Equal( leftHandSide->GetSpacing(), rightHandSide->GetSpacing(), eps ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Spacing differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetSpacing() << " : leftHandSide is " << leftHandSide->GetSpacing() << " and tolerance is " << eps; } result = false; } //Compare Origins if( !mitk::Equal( leftHandSide->GetOrigin(), rightHandSide->GetOrigin(), eps ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Origin differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetOrigin() << " : leftHandSide is " << leftHandSide->GetOrigin() << " and tolerance is " << eps; } result = false; } //Compare Axis and Extents for( unsigned int i=0; i<3; ++i) { if( !mitk::Equal( leftHandSide->GetAxisVector(i), rightHandSide->GetAxisVector(i), eps)) { if(verbose) { MITK_INFO << "[( BaseGeometry )] AxisVector #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetAxisVector(i) << " : leftHandSide is " << leftHandSide->GetAxisVector(i) << " and tolerance is " << eps; } result = false; } if( !mitk::Equal( leftHandSide->GetExtent(i), rightHandSide->GetExtent(i), eps) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Extent #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetExtent(i) << " : leftHandSide is " << leftHandSide->GetExtent(i) << " and tolerance is " << eps; } result = false; } } //Compare BoundingBoxes if( !mitk::Equal( leftHandSide->GetBoundingBox(), rightHandSide->GetBoundingBox(), eps, verbose) ) { result = false; } //Compare IndexToWorldTransform Matrix if( !mitk::Equal( leftHandSide->GetIndexToWorldTransform(), rightHandSide->GetIndexToWorldTransform(), eps, verbose) ) { result = false; } return result; } -void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D& aSpacing) +void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D& aSpacing, bool enforceSetSpacing ) { InternPreSetSpacing(aSpacing); - InternSetSpacing(aSpacing); + InternSetSpacing(aSpacing, enforceSetSpacing); } void mitk::BaseGeometry::InternPreSetSpacing(const mitk::Vector3D& aSpacing) {} -void mitk::BaseGeometry::InternSetSpacing(const mitk::Vector3D& aSpacing){ - if(mitk::Equal(m_Spacing, aSpacing) == false) +void mitk::BaseGeometry::InternSetSpacing(const mitk::Vector3D& aSpacing, bool enforceSetSpacing){ + if(mitk::Equal(m_Spacing, aSpacing) == false || enforceSetSpacing) { assert(aSpacing[0]>0 && aSpacing[1]>0 && aSpacing[2]>0); m_Spacing = aSpacing; AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = m_IndexToWorldTransform->GetMatrix().GetVnlMatrix(); mitk::VnlVector col; col = vnlmatrix.get_column(0); col.normalize(); col*=aSpacing[0]; vnlmatrix.set_column(0, col); col = vnlmatrix.get_column(1); col.normalize(); col*=aSpacing[1]; vnlmatrix.set_column(1, col); col = vnlmatrix.get_column(2); col.normalize(); col*=aSpacing[2]; vnlmatrix.set_column(2, col); Matrix3D matrix; matrix = vnlmatrix; AffineTransform3D::Pointer transform = AffineTransform3D::New(); transform->SetMatrix(matrix); transform->SetOffset(m_IndexToWorldTransform->GetOffset()); SetIndexToWorldTransform(transform.GetPointer()); itk2vtk(m_Spacing, m_FloatSpacing); } } mitk::Vector3D mitk::BaseGeometry::GetAxisVector(unsigned int direction) const { Vector3D frontToBack; frontToBack.SetVnlVector(m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction)); frontToBack *= GetExtent(direction); return frontToBack; } mitk::ScalarType mitk::BaseGeometry::GetExtent(unsigned int direction) const { assert(m_BoundingBox.IsNotNull()); - if (direction>=NDimensions) + if (direction>=m_NDimensions) mitkThrow() << "Direction is too big. This geometry is for 3D Data"; BoundsArrayType bounds = m_BoundingBox->GetBounds(); return bounds[direction*2+1]-bounds[direction*2]; } bool mitk::BaseGeometry::Is2DConvertable() { bool isConvertableWithoutLoss = true; do { if (this->GetSpacing()[2] != 1) { isConvertableWithoutLoss = false; break; } if (this->GetOrigin()[2] != 0) { isConvertableWithoutLoss = false; break; } mitk::Vector3D col0, col1, col2; col0.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0)); col1.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1)); col2.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2)); if ((col0[2] != 0) || (col1[2] != 0) || (col2[0] != 0) || (col2[1] != 0) || (col2[2] != 1)) { isConvertableWithoutLoss = false; break; } } while (0); return isConvertableWithoutLoss; } mitk::Point3D mitk::BaseGeometry::GetCenter() const { assert(m_BoundingBox.IsNotNull()); return m_IndexToWorldTransform->TransformPoint(m_BoundingBox->GetCenter()); } double mitk::BaseGeometry::GetDiagonalLength2() const { Vector3D diagonalvector = GetCornerPoint()-GetCornerPoint(false, false, false); return diagonalvector.GetSquaredNorm(); } //##Documentation //## @brief Get the length of the diagonal of the bounding-box in mm //## double mitk::BaseGeometry::GetDiagonalLength() const { return sqrt(GetDiagonalLength2()); } mitk::Point3D mitk::BaseGeometry::GetCornerPoint(int id) const { assert(id >= 0); assert(m_BoundingBox.IsNotNull()); BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; switch(id) { case 0: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[4]); break; case 1: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[5]); break; case 2: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[4]); break; case 3: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[5]); break; case 4: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[4]); break; case 5: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[5]); break; case 6: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[4]); break; case 7: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[5]); break; default: { itkExceptionMacro(<<"A cube only has 8 corners. These are labeled 0-7."); } } return m_IndexToWorldTransform->TransformPoint(cornerpoint); } mitk::Point3D mitk::BaseGeometry::GetCornerPoint(bool xFront, bool yFront, bool zFront) const { assert(m_BoundingBox.IsNotNull()); BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; cornerpoint[0] = (xFront ? bounds[0] : bounds[1]); cornerpoint[1] = (yFront ? bounds[2] : bounds[3]); cornerpoint[2] = (zFront ? bounds[4] : bounds[5]); return m_IndexToWorldTransform->TransformPoint(cornerpoint); } mitk::ScalarType mitk::BaseGeometry::GetExtentInMM(int direction) const { return m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction).magnitude()*GetExtent(direction); } void mitk::BaseGeometry::SetExtentInMM(int direction, ScalarType extentInMM) { ScalarType len = GetExtentInMM(direction); if(fabs(len - extentInMM)>=mitk::eps) { AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = m_IndexToWorldTransform->GetMatrix().GetVnlMatrix(); if(len>extentInMM) vnlmatrix.set_column(direction, vnlmatrix.get_column(direction)/len*extentInMM); else vnlmatrix.set_column(direction, vnlmatrix.get_column(direction)*extentInMM/len); Matrix3D matrix; matrix = vnlmatrix; m_IndexToWorldTransform->SetMatrix(matrix); Modified(); } InternPostSetExtentInMM(direction,extentInMM); } void mitk::BaseGeometry::InternPostSetExtentInMM(int direction, ScalarType extentInMM){}; bool mitk::BaseGeometry::IsInside(const mitk::Point3D& p) const { mitk::Point3D index; WorldToIndex(p, index); return IsIndexInside(index); } bool mitk::BaseGeometry::IsIndexInside(const mitk::Point3D& index) const { bool inside = false; //if it is an image geometry, we need to convert the index to discrete values //this is done by applying the rounding function also used in WorldToIndex (see line 323) inside = m_BoundingBox->IsInside(index); return inside; } //##Documentation //## @brief Convenience method for working with ITK indices template bool mitk::BaseGeometry::IsIndexInside(const itk::Index &index) const { int i, dim=index.GetIndexDimension(); Point3D pt_index; pt_index.Fill(0); for ( i = 0; i < dim; ++i ) { pt_index[i] = index[i]; } return IsIndexInside(pt_index); } void mitk::BaseGeometry::WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const { BackTransform(pt_mm, pt_units); } void mitk::BaseGeometry::WorldToIndex( const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const { BackTransform( vec_mm, vec_units); } void mitk::BaseGeometry::BackTransform(const mitk::Vector3D& in, mitk::Vector3D& out) const { // Get WorldToIndex transform if (m_IndexToWorldTransformLastModified != m_IndexToWorldTransform->GetMTime()) { m_InvertedTransform = TransformType::New(); if (!m_IndexToWorldTransform->GetInverse( m_InvertedTransform.GetPointer() )) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed." ); } m_IndexToWorldTransformLastModified = m_IndexToWorldTransform->GetMTime(); } // Check for valid matrix inversion const TransformType::MatrixType& inverse = m_InvertedTransform->GetMatrix(); if(inverse.GetVnlMatrix().has_nans()) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed. Matrix was: " << std::endl << m_IndexToWorldTransform->GetMatrix() << "Suggested inverted matrix is:" << std::endl << inverse ); } // Transform vector for (unsigned int i = 0; i < 3; i++) { out[i] = 0.0; for (unsigned int j = 0; j < 3; j++) { out[i] += inverse[i][j]*in[j]; } } } void mitk::BaseGeometry::BackTransform(const mitk::Point3D &in, mitk::Point3D& out) const { ScalarType temp[3]; unsigned int i, j; const TransformType::OffsetType& offset = m_IndexToWorldTransform->GetOffset(); // Remove offset for (j = 0; j < 3; j++) { temp[j] = in[j] - offset[j]; } // Get WorldToIndex transform if (m_IndexToWorldTransformLastModified != m_IndexToWorldTransform->GetMTime()) { m_InvertedTransform = TransformType::New(); if (!m_IndexToWorldTransform->GetInverse( m_InvertedTransform.GetPointer() )) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed." ); } m_IndexToWorldTransformLastModified = m_IndexToWorldTransform->GetMTime(); } // Check for valid matrix inversion const TransformType::MatrixType& inverse = m_InvertedTransform->GetMatrix(); if(inverse.GetVnlMatrix().has_nans()) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed. Matrix was: " << std::endl << m_IndexToWorldTransform->GetMatrix() << "Suggested inverted matrix is:" << std::endl << inverse ); } // Transform point for (i = 0; i < 3; i++) { out[i] = 0.0; for (j = 0; j < 3; j++) { out[i] += inverse[i][j]*temp[j]; } } } mitk::VnlVector mitk::BaseGeometry::GetOriginVnl() const { return const_cast(this)->m_Origin.GetVnlVector(); } vtkLinearTransform* mitk::BaseGeometry::GetVtkTransform() const { return (vtkLinearTransform*)m_VtkIndexToWorldTransform; } void mitk::BaseGeometry::SetIdentity() { m_IndexToWorldTransform->SetIdentity(); m_Origin.Fill(0); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); //NEW!!! Modified(); TransferItkToVtkTransform(); } void mitk::BaseGeometry::TransferVtkToItkTransform() { TransferVtkMatrixToItkTransform(m_VtkMatrix, m_IndexToWorldTransform.GetPointer()); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); } void mitk::BaseGeometry::Compose( const mitk::BaseGeometry::TransformType * other, bool pre ) { m_IndexToWorldTransform->Compose(other, pre); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); Modified(); TransferItkToVtkTransform(); } void mitk::BaseGeometry::Compose( const vtkMatrix4x4 * vtkmatrix, bool pre ) { mitk::BaseGeometry::TransformType::Pointer itkTransform = mitk::BaseGeometry::TransformType::New(); TransferVtkMatrixToItkTransform(vtkmatrix, itkTransform.GetPointer()); Compose(itkTransform, pre); } void mitk::BaseGeometry::Translate(const Vector3D & vector) { if((vector[0] != 0) || (vector[1] != 0) || (vector[2] != 0)) { this->SetOrigin(m_Origin + vector); } } void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const { pt_mm = m_IndexToWorldTransform->TransformPoint(pt_units); } void mitk::BaseGeometry::IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const { vec_mm = m_IndexToWorldTransform->TransformVector(vec_units); } #include void mitk::BaseGeometry::ExecuteOperation(Operation* operation) { vtkTransform *vtktransform = vtkTransform::New(); vtktransform->SetMatrix(m_VtkMatrix); switch (operation->GetOperationType()) { case OpNOTHING: break; case OpMOVE: { mitk::PointOperation *pointOp = dynamic_cast(operation); if (pointOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } mitk::Point3D newPos = pointOp->GetPoint(); ScalarType data[3]; vtktransform->GetPosition(data); vtktransform->PostMultiply(); vtktransform->Translate(newPos[0], newPos[1], newPos[2]); vtktransform->PreMultiply(); break; } case OpSCALE: { mitk::PointOperation *pointOp = dynamic_cast(operation); if (pointOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } mitk::Point3D newScale = pointOp->GetPoint(); ScalarType data[3]; /* calculate new scale: newscale = oldscale * (oldscale + scaletoadd)/oldscale */ data[0] = 1 + (newScale[0] / GetMatrixColumn(0).magnitude()); data[1] = 1 + (newScale[1] / GetMatrixColumn(1).magnitude()); data[2] = 1 + (newScale[2] / GetMatrixColumn(2).magnitude()); mitk::Point3D center = const_cast(m_BoundingBox.GetPointer())->GetCenter(); ScalarType pos[3]; vtktransform->GetPosition(pos); vtktransform->PostMultiply(); vtktransform->Translate(-pos[0], -pos[1], -pos[2]); vtktransform->Translate(-center[0], -center[1], -center[2]); vtktransform->PreMultiply(); vtktransform->Scale(data[0], data[1], data[2]); vtktransform->PostMultiply(); vtktransform->Translate(+center[0], +center[1], +center[2]); vtktransform->Translate(pos[0], pos[1], pos[2]); vtktransform->PreMultiply(); break; } case OpROTATE: { mitk::RotationOperation *rotateOp = dynamic_cast(operation); if (rotateOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } Vector3D rotationVector = rotateOp->GetVectorOfRotation(); Point3D center = rotateOp->GetCenterOfRotation(); ScalarType angle = rotateOp->GetAngleOfRotation(); vtktransform->PostMultiply(); vtktransform->Translate(-center[0], -center[1], -center[2]); vtktransform->RotateWXYZ(angle, rotationVector[0], rotationVector[1], rotationVector[2]); vtktransform->Translate(center[0], center[1], center[2]); vtktransform->PreMultiply(); break; } case OpRESTOREPLANEPOSITION: { //Copy necessary to avoid vtk warning vtkMatrix4x4* matrix = vtkMatrix4x4::New(); TransferItkTransformToVtkMatrix(dynamic_cast(operation)->GetTransform().GetPointer(), matrix); vtktransform->SetMatrix(matrix); break; } case OpAPPLYTRANSFORMMATRIX: { ApplyTransformMatrixOperation *applyMatrixOp = dynamic_cast< ApplyTransformMatrixOperation* >( operation ); vtktransform->SetMatrix(applyMatrixOp->GetMatrix()); break; } default: vtktransform->Delete(); return; } m_VtkMatrix->DeepCopy(vtktransform->GetMatrix()); TransferVtkToItkTransform(); Modified(); vtktransform->Delete(); } mitk::VnlVector mitk::BaseGeometry::GetMatrixColumn(unsigned int direction) const { return m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction); } mitk::BoundingBox::Pointer mitk::BaseGeometry::CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D* transform) const { mitk::BoundingBox::PointsContainer::Pointer pointscontainer=mitk::BoundingBox::PointsContainer::New(); mitk::BoundingBox::PointIdentifier pointid=0; unsigned char i; if(transform!=NULL) { mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New(); transform->GetInverse(inverse); for(i=0; i<8; ++i) pointscontainer->InsertElement( pointid++, inverse->TransformPoint( GetCornerPoint(i) )); } else { for(i=0; i<8; ++i) pointscontainer->InsertElement( pointid++, GetCornerPoint(i) ); } mitk::BoundingBox::Pointer result = mitk::BoundingBox::New(); result->SetPoints(pointscontainer); result->ComputeBoundingBox(); return result; } void mitk::BaseGeometry::SetTimeBounds(const TimeBounds& timebounds) { if(m_TimeBounds != timebounds) { m_TimeBounds = timebounds; Modified(); } InternPostSetTimeBounds(timebounds); } void mitk::BaseGeometry::InternPostSetTimeBounds(const TimeBounds& timebounds) {} const std::string mitk::BaseGeometry::GetTransformAsString( TransformType* transformType ) { std::ostringstream out; out << '['; for( int i=0; i<3; ++i ) { out << '['; for( int j=0; j<3; ++j ) out << transformType->GetMatrix().GetVnlMatrix().get(i, j) << ' '; out << ']'; } out << "]["; for( int i=0; i<3; ++i ) out << transformType->GetOffset()[i] << ' '; out << "]\0"; return out.str(); } void mitk::BaseGeometry::SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4* vtkmatrix) { m_VtkMatrix->DeepCopy(vtkmatrix); TransferVtkToItkTransform(); } void mitk::BaseGeometry::WorldToIndex(const mitk::Point3D & /*atPt3d_mm*/, const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const { MITK_WARN<<"Warning! Call of the deprecated function BaseGeometry::WorldToIndex(point, vec, vec). Use BaseGeometry::WorldToIndex(vec, vec) instead!"; //BackTransform(atPt3d_mm, vec_mm, vec_units); this->WorldToIndex(vec_mm, vec_units); } void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D &/*atPt3d_units*/, const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const { MITK_WARN<<"Warning! Call of the deprecated function BaseGeometry::IndexToWorld(point, vec, vec). Use BaseGeometry::IndexToWorld(vec, vec) instead!"; //vec_mm = m_IndexToWorldTransform->TransformVector(vec_units); this->IndexToWorld(vec_units, vec_mm); } void mitk::BaseGeometry::BackTransform(const mitk::Point3D &/*at*/, const mitk::Vector3D &in, mitk::Vector3D& out) const { MITK_INFO<<"Warning! Call of the deprecated function BaseGeometry::BackTransform(point, vec, vec). Use BaseGeometry::BackTransform(vec, vec) instead!"; //// Get WorldToIndex transform //if (m_IndexToWorldTransformLastModified != m_IndexToWorldTransform->GetMTime()) //{ // m_InvertedTransform = TransformType::New(); // if (!m_IndexToWorldTransform->GetInverse( m_InvertedTransform.GetPointer() )) // { // itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed." ); // } // m_IndexToWorldTransformLastModified = m_IndexToWorldTransform->GetMTime(); //} //// Check for valid matrix inversion //const TransformType::MatrixType& inverse = m_InvertedTransform->GetMatrix(); //if(inverse.GetVnlMatrix().has_nans()) //{ // itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed. Matrix was: " << std::endl // << m_IndexToWorldTransform->GetMatrix() << "Suggested inverted matrix is:" << std::endl // << inverse ); //} //// Transform vector //for (unsigned int i = 0; i < 3; i++) //{ // out[i] = 0.0; // for (unsigned int j = 0; j < 3; j++) // { // out[i] += inverse[i][j]*in[j]; // } //} this->BackTransform(in, out); } + +vtkMatrix4x4* mitk::BaseGeometry::GetVtkMatrix(){ + return m_VtkMatrix; +} + +bool mitk::BaseGeometry::IsBoundingBoxNull() const{ + return m_BoundingBox.IsNull(); +} + +bool mitk::BaseGeometry::IsIndexToWorldTransformNull() const{ + return m_IndexToWorldTransform.IsNull(); +} diff --git a/Core/Code/DataManagement/mitkBaseGeometry.h b/Core/Code/DataManagement/mitkBaseGeometry.h index 77e8ad4883..33874a3964 100644 --- a/Core/Code/DataManagement/mitkBaseGeometry.h +++ b/Core/Code/DataManagement/mitkBaseGeometry.h @@ -1,585 +1,584 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef BaseGeometry_H_HEADER_INCLUDED #define BaseGeometry_H_HEADER_INCLUDED #include #include #include "mitkoperationactor.h" #include #include "mitkvector.h" #include #include #include "itkScalableAffineTransform.h" #include class vtkMatrix4x4; class vtkMatrixToLinearTransform; class vtkLinearTransform; namespace mitk { //##Documentation //## @brief Standard 3D-BoundingBox typedef //## //## Standard 3D-BoundingBox typedef to get rid of template arguments (3D, type). typedef itk::BoundingBox BoundingBox; //##Documentation //## @brief Standard typedef for time-bounds typedef itk::FixedArray TimeBounds; typedef itk::FixedArray FixedArrayType; //##Documentation //## @brief BaseGeometry xxxxxxxxxxxxxx //## //## xxxxxxxxxxx //## //## Rule: everything is in mm (ms) if not stated otherwise. //## @ingroup Geometry class MITK_CORE_EXPORT BaseGeometry : public itk::Object, public OperationActor { public: mitkClassMacro(BaseGeometry, itk::Object); - //void testXYZ(){}; //xxx - // ********************************** TypeDef ********************************** typedef itk::QuaternionRigidTransform< ScalarType > QuaternionTransformType; typedef QuaternionTransformType::VnlQuaternionType VnlQuaternionType; typedef itk::ScalableAffineTransform TransformType; typedef itk::BoundingBox BoundingBoxType; typedef BoundingBoxType::BoundsArrayType BoundsArrayType; typedef BoundingBoxType::Pointer BoundingBoxPointer; // ********************************** Origin, Spacing ********************************** //##Documentation //## @brief Get the origin, e.g. the upper-left corner of the plane const Point3D& GetOrigin() const; //##Documentation //## @brief Set the origin, i.e. the upper-left corner of the plane //## void SetOrigin(const Point3D& origin); //##Documentation //## @brief Get the spacing (size of a pixel). //## itkGetConstReferenceMacro(Spacing, mitk::Vector3D); //##Documentation //## @brief Get the spacing as a float[3] array. const float* GetFloatSpacing() const; //##Documentation //## @brief Set the spacing (m_Spacing) - void SetSpacing(const mitk::Vector3D& aSpacing); + void SetSpacing(const mitk::Vector3D& aSpacing, bool enforceSetSpacing = false); //##Documentation //## @brief Get the origin as VnlVector //## //## \sa GetOrigin VnlVector GetOriginVnl() const; // ********************************** other functions ********************************** //##Documentation //## @brief Get the DICOM FrameOfReferenceID referring to the //## used world coordinate system itkGetConstMacro(FrameOfReferenceID, unsigned int); //##Documentation //## @brief Set the DICOM FrameOfReferenceID referring to the //## used world coordinate system itkSetMacro(FrameOfReferenceID, unsigned int); + itkGetConstMacro(IndexToWorldTransformLastModified, unsigned long); + //##Documentation //## @brief Is this Geometry3D in a state that is valid? - bool IsValid() const; + virtual bool IsValid() const; // ********************************** Initialize ********************************** //##Documentation //## @brief Initialize the Geometry3D void Initialize(); void InitializeGeometry(Self * newGeometry) const; static void CopySpacingFromTransform(mitk::AffineTransform3D* transform, mitk::Vector3D& spacing, float floatSpacing[3]); // ********************************** Transformations Set/Get ********************************** // a bit of a misuse, but we want only doxygen to see the following: #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get the transformation used to convert from index //## to world coordinates itkGetObjectMacro(IndexToWorldTransform, AffineTransform3D); #endif //## @brief Set the transformation used to convert from index //## to world coordinates void SetIndexToWorldTransform(mitk::AffineTransform3D* transform); //##Documentation //## @brief Convenience method for setting the ITK transform //## (m_IndexToWorldTransform) via an vtkMatrix4x4 //## \sa SetIndexToWorldTransform virtual void SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4* vtkmatrix); /** Set/Get the IndexToWorldTransform */ itkGetConstObjectMacro(IndexToWorldTransform, AffineTransform3D); itkGetObjectMacro(IndexToWorldTransform, AffineTransform3D); + vtkMatrix4x4* GetVtkMatrix(); + //##Documentation //## @brief Get the m_IndexToWorldTransform as a vtkLinearTransform vtkLinearTransform* GetVtkTransform() const; //##Documentation //## @brief Set the transform to identity and origin to 0 //## virtual void SetIdentity(); // ********************************** Transformations ********************************** //##Documentation //## @brief Copy the ITK transform //## (m_IndexToWorldTransform) to the VTK transform //## \sa SetIndexToWorldTransform void TransferItkToVtkTransform(); //##Documentation //## @brief Copy the VTK transform //## to the ITK transform (m_IndexToWorldTransform) //## \sa SetIndexToWorldTransform void TransferVtkToItkTransform(); //##Documentation //## @brief Compose new IndexToWorldTransform with a given transform. //## //## This method composes m_IndexToWorldTransform with another transform, //## modifying self to be the composition of self and other. //## If the argument pre is true, then other is precomposed with self; //## that is, the resulting transformation consists of first applying //## other to the source, followed by self. If pre is false or omitted, //## then other is post-composed with self; that is the resulting //## transformation consists of first applying self to the source, //## followed by other. void Compose( const BaseGeometry::TransformType * other, bool pre = 0 ); //##Documentation //## @brief Compose new IndexToWorldTransform with a given vtkMatrix4x4. //## //## Converts the vtkMatrix4x4 into a itk-transform and calls the previous method. void Compose( const vtkMatrix4x4 * vtkmatrix, bool pre = 0 ); //##Documentation //## @brief Translate the origin by a vector //## void Translate(const Vector3D& vector); //##Documentation //##@brief executes affine operations (translate, rotate, scale) virtual void ExecuteOperation(Operation* operation); //##Documentation //## @brief Convert world coordinates (in mm) of a \em point to (continuous!) index coordinates //## \warning If you need (discrete) integer index coordinates (e.g., for iterating easily over an image), //## use WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index). //## For further information about coordinates types, please see the Geometry documentation void WorldToIndex(const mitk::Point3D& pt_mm, mitk::Point3D& pt_units) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em vector //## \a vec_mm to (continuous!) index coordinates. //## For further information about coordinates types, please see the Geometry documentation void WorldToIndex(const mitk::Vector3D& vec_mm, mitk::Vector3D& vec_units) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em point to (discrete!) index coordinates. //## This method rounds to integer indices! //## For further information about coordinates types, please see the Geometry documentation template void WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index) const { typedef itk::Index IndexType; mitk::Point3D pt_units; this->WorldToIndex(pt_mm, pt_units); int i, dim=index.GetIndexDimension(); if(dim>3) { index.Fill(0); dim=3; } for(i=0;i( pt_units[i] ); } } //##Documentation //## @brief Convert (continuous or discrete) index coordinates of a \em vector //## \a vec_units to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation void IndexToWorld(const mitk::Vector3D& vec_units, mitk::Vector3D& vec_mm) const; //##Documentation //## @brief Convert (continuous or discrete) index coordinates of a \em point to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation void IndexToWorld(const mitk::Point3D& pt_units, mitk::Point3D& pt_mm) const; //##Documentation //## @brief Convert (discrete) index coordinates of a \em point to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation template void IndexToWorld(const itk::Index &index, mitk::Point3D& pt_mm ) const { mitk::Point3D pt_units; pt_units.Fill(0); int i, dim=index.GetIndexDimension(); if(dim>3) { dim=3; } for(i=0;i void ItkPhysicalPointToWorld(const itk::Point& itkPhysicalPoint, mitk::Point3D& pt_mm) const { mitk::vtk2itk(itkPhysicalPoint, pt_mm); } //##Documentation //## @brief Deprecated for use with ITK version 3.10 or newer. //## Convert world coordinates (in mm) of a \em point to //## ITK physical coordinates (in mm, but without a possible rotation) //## //## This method is useful if you have want to access an mitk::Image //## via an itk::Image. ITK v3.8 and older did not support rotated (tilted) //## images, i.e., ITK images are always parallel to the coordinate axes. //## When accessing a (possibly rotated) mitk::Image via an itk::Image //## the rotational part of the transformation in the Geometry3D is //## simply discarded; in other word: only the origin and spacing is //## used by ITK, not the complete matrix available in MITK. //## With WorldToItkPhysicalPoint you can convert an MITK world //## coordinate (including the rotation) into a coordinate that //## can be used with the ITK image as a ITK physical coordinate //## (excluding the rotation). template void WorldToItkPhysicalPoint(const mitk::Point3D& pt_mm, itk::Point& itkPhysicalPoint) const { mitk::vtk2itk(pt_mm, itkPhysicalPoint); } // ********************************** BoundingBox ********************************** /** Get the bounding box */ itkGetConstObjectMacro(BoundingBox, BoundingBoxType); //##Documentation //## @brief Get the time bounds (in ms) itkGetConstReferenceMacro(TimeBounds, TimeBounds); // a bit of a misuse, but we want only doxygen to see the following: #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get bounding box (in index/unit coordinates) itkGetConstObjectMacro(BoundingBox, BoundingBoxType); //##Documentation //## @brief Get bounding box (in index/unit coordinates) as a BoundsArrayType const BoundsArrayType GetBounds() const; #endif const BoundsArrayType GetBounds() const; //##Documentation //## \brief Set the bounding box (in index/unit coordinates) //## //## Only possible via the BoundsArray to make clear that a //## copy of the bounding-box is stored, not a reference to it. void SetBounds(const BoundsArrayType& bounds); //##Documentation //## @brief Set the bounding box (in index/unit coordinates) via a float array void SetFloatBounds(const float bounds[6]); //##Documentation //## @brief Set the bounding box (in index/unit coordinates) via a double array void SetFloatBounds(const double bounds[6]); //##Documentation //## @brief Get a VnlVector along bounding-box in the specified //## @a direction, length is spacing //## //## \sa GetAxisVector VnlVector GetMatrixColumn(unsigned int direction) const; //##Documentation //## @brief Calculates a bounding-box around the geometry relative //## to a coordinate system defined by a transform //## mitk::BoundingBox::Pointer CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D* transform) const; //##Documentation //## @brief Set the time bounds (in ms) void SetTimeBounds(const TimeBounds& timebounds); // ********************************** Geometry ********************************** #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get the extent of the bounding box (in index/unit coordinates) //## //## To access the extent in mm use GetExtentInMM ScalarType GetExtent(unsigned int direction) const; #endif /** Get the extent of the bounding box */ ScalarType GetExtent(unsigned int direction) const; //##Documentation //## @brief Get the extent of the bounding-box in the specified @a direction in mm //## //## Equals length of GetAxisVector(direction). ScalarType GetExtentInMM(int direction) const; //##Documentation //## @brief Get vector along bounding-box in the specified @a direction in mm //## //## The length of the vector is the size of the bounding-box in the //## specified @a direction in mm //## \sa GetMatrixColumn Vector3D GetAxisVector(unsigned int direction) const; //##Documentation //## @brief 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 virtual bool Is2DConvertable(); //##Documentation //## @brief Get the center of the bounding-box in mm //## Point3D GetCenter() const; //##Documentation //## @brief Get the squared length of the diagonal of the bounding-box in mm //## double GetDiagonalLength2() const; //##Documentation //## @brief Get the length of the diagonal of the bounding-box in mm //## double GetDiagonalLength() const; //##Documentation //## @brief Get the position of the corner number \a id (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. virtual Point3D GetCornerPoint(int id) const; //##Documentation //## @brief Get the position of a corner (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. virtual Point3D GetCornerPoint(bool xFront=true, bool yFront=true, bool zFront=true) const; //##Documentation //## @brief Set the extent of the bounding-box in the specified @a direction in mm //## //## @note This changes the matrix in the transform, @a not the bounds, which are given in units! void SetExtentInMM(int direction, ScalarType extentInMM); //##Documentation //## @brief Test whether the point \a p (world coordinates in mm) is //## inside the bounding box virtual bool IsInside(const mitk::Point3D& p) const; //##Documentation //## @brief Test whether the point \a p ((continous!)index coordinates in units) is //## inside the bounding box virtual bool IsIndexInside(const mitk::Point3D& index) const; //##Documentation //## @brief Convenience method for working with ITK indices template bool IsIndexInside(const itk::Index &index) const; protected: // ********************************** Constructor ********************************** BaseGeometry(); BaseGeometry(const BaseGeometry& other); virtual ~BaseGeometry(); - itkGetConstMacro(IndexToWorldTransformLastModified, unsigned long); - void BackTransform(const mitk::Point3D& in, mitk::Point3D& out) const; //Without redundant parameter Point3D void BackTransform(const mitk::Vector3D& in, mitk::Vector3D& out) const; //##Documentation //## @brief Deprecated void BackTransform(const mitk::Point3D& at, const mitk::Vector3D& in, mitk::Vector3D& out) const; static const std::string GetTransformAsString( TransformType* transformType ); - //Internal Functions - virtual bool InternPostIsValid() const; - virtual void InternPostInitialize() {}; virtual void InternPostInitializeGeometry(Self * newGeometry) const{}; virtual void InternPreSetBounds(const BoundsArrayType& bounds); virtual void InternPostSetExtentInMM(int direction, ScalarType extentInMM); virtual void InternPostSetTimeBounds(const TimeBounds& timebounds); virtual void InternPreSetIndexToWorldTransform(mitk::AffineTransform3D* transform); virtual void InternPostSetIndexToWorldTransform(mitk::AffineTransform3D* transform); virtual void InternPreSetSpacing(const mitk::Vector3D& aSpacing); - void InternSetSpacing(const mitk::Vector3D& aSpacing); + void InternSetSpacing(const mitk::Vector3D& aSpacing, bool enforceSetSpacing = false); + + itkGetConstMacro(NDimensions, unsigned int); + + bool IsBoundingBoxNull() const; + + bool IsIndexToWorldTransformNull() const; + private: // ********************************** Variables ********************************** + mitk::Vector3D m_Spacing; + AffineTransform3D::Pointer m_IndexToWorldTransform; + mutable BoundingBoxPointer m_BoundingBox; + vtkMatrixToLinearTransform* m_VtkIndexToWorldTransform; vtkMatrix4x4* m_VtkMatrix; - bool m_Valid; - unsigned int m_FrameOfReferenceID; mutable mitk::TimeBounds m_TimeBounds; - mutable BoundingBoxPointer m_BoundingBox; - //##Documentation //## @brief Origin, i.e. upper-left corner of the plane //## Point3D m_Origin; - //##Documentation - //## @brief Spacing of the data. Only significant if the geometry describes - //## an Image (m_ImageGeometry==true). - mitk::Vector3D m_Spacing; - - static const unsigned int NDimensions = 3; + static const unsigned int m_NDimensions = 3; - mutable TransformType::Pointer m_InvertedTransform; //this was private + mutable TransformType::Pointer m_InvertedTransform; - mutable unsigned long m_IndexToWorldTransformLastModified; //this was private + mutable unsigned long m_IndexToWorldTransformLastModified; float m_FloatSpacing[3]; //this was private - // DEPRECATED(VnlQuaternionType m_RotationQuaternion); //this was private + // DEPRECATED(VnlQuaternionType m_RotationQuaternion); }; // ********************************** Equal Functions ********************************** // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two bounding boxes (BoundingBoxType) for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the bounds (elementwise). * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal( const mitk::BaseGeometry::BoundingBoxType *leftHandSide, const mitk::BaseGeometry::BoundingBoxType *rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two geometries for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the spacing, origin, axisvectors, extents, the matrix of the * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag. * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * If you want to use different tolarance values for different parts of the geometry, feel free to use * the other comparison methods and write your own implementation of Equal. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::BaseGeometry* leftHandSide, const mitk::BaseGeometry* rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo /** * @brief Equal A function comparing two transforms (TransformType) for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the IndexToWorldTransform (elementwise). * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::BaseGeometry::TransformType *leftHandSide, const mitk::BaseGeometry::TransformType *rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo } // namespace mitk #endif /* BaseGeometry_H_HEADER_INCLUDED */ diff --git a/Core/Code/DataManagement/mitkDisplayGeometry.cpp b/Core/Code/DataManagement/mitkDisplayGeometry.cpp index 087299af7d..99b8e331ec 100644 --- a/Core/Code/DataManagement/mitkDisplayGeometry.cpp +++ b/Core/Code/DataManagement/mitkDisplayGeometry.cpp @@ -1,613 +1,613 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDisplayGeometry.h" itk::LightObject::Pointer mitk::DisplayGeometry::InternalClone() const { // itkExceptionMacro(<<"calling mitk::DisplayGeometry::Clone does not make much sense."); DisplayGeometry* returnValue = const_cast(this); return returnValue; } -bool mitk::DisplayGeometry::InternPostIsValid() const +bool mitk::DisplayGeometry::IsValid() const { return m_WorldGeometry.IsNotNull() && m_WorldGeometry->IsValid(); } unsigned long mitk::DisplayGeometry::GetMTime() const { if((m_WorldGeometry.IsNotNull()) && (Geometry2D::GetMTime() < m_WorldGeometry->GetMTime())) { Modified(); } return Geometry2D::GetMTime(); } const mitk::TimeBounds& mitk::DisplayGeometry::GetTimeBounds() const { if(m_WorldGeometry.IsNull()) { - return m_TimeBounds; + return this->GetTimeBounds(); } return m_WorldGeometry->GetTimeBounds(); } // size definition methods void mitk::DisplayGeometry::SetWorldGeometry(const Geometry2D* aWorldGeometry) { m_WorldGeometry = aWorldGeometry; Modified(); } bool mitk::DisplayGeometry::SetOriginInMM(const Vector2D& origin_mm) { m_OriginInMM = origin_mm; WorldToDisplay(m_OriginInMM, m_OriginInDisplayUnits); Modified(); return !this->RefitVisibleRect(); } mitk::Vector2D mitk::DisplayGeometry::GetOriginInMM() const { return m_OriginInMM; } mitk::Vector2D mitk::DisplayGeometry::GetOriginInDisplayUnits() const { return m_OriginInDisplayUnits; } void mitk::DisplayGeometry::SetSizeInDisplayUnits(unsigned int width, unsigned int height, bool keepDisplayedRegion) { Vector2D oldSizeInMM( m_SizeInMM ); Point2D oldCenterInMM; if(keepDisplayedRegion) { Point2D centerInDisplayUnits; centerInDisplayUnits[0] = m_SizeInDisplayUnits[0]*0.5; centerInDisplayUnits[1] = m_SizeInDisplayUnits[1]*0.5; DisplayToWorld(centerInDisplayUnits, oldCenterInMM); } m_SizeInDisplayUnits[0]=width; m_SizeInDisplayUnits[1]=height; if(m_SizeInDisplayUnits[0] <= 0) m_SizeInDisplayUnits[0] = 1; if(m_SizeInDisplayUnits[1] <= 0) m_SizeInDisplayUnits[1] = 1; DisplayToWorld(m_SizeInDisplayUnits, m_SizeInMM); if(keepDisplayedRegion) { Point2D positionOfOldCenterInCurrentDisplayUnits; WorldToDisplay(oldCenterInMM, positionOfOldCenterInCurrentDisplayUnits); Point2D currentNewCenterInDisplayUnits; currentNewCenterInDisplayUnits[0] = m_SizeInDisplayUnits[0]*0.5; currentNewCenterInDisplayUnits[1] = m_SizeInDisplayUnits[1]*0.5; Vector2D shift; shift=positionOfOldCenterInCurrentDisplayUnits.GetVectorFromOrigin()-currentNewCenterInDisplayUnits; MoveBy(shift); Zoom(m_SizeInMM.GetNorm()/oldSizeInMM.GetNorm(), currentNewCenterInDisplayUnits); } Modified(); } mitk::Vector2D mitk::DisplayGeometry::GetSizeInDisplayUnits() const { return m_SizeInDisplayUnits; } mitk::Vector2D mitk::DisplayGeometry::GetSizeInMM() const { return m_SizeInMM; } unsigned int mitk::DisplayGeometry::GetDisplayWidth() const { assert(m_SizeInDisplayUnits[0] >= 0); return (unsigned int)m_SizeInDisplayUnits[0]; } unsigned int mitk::DisplayGeometry::GetDisplayHeight() const { assert(m_SizeInDisplayUnits[1] >= 0); return (unsigned int)m_SizeInDisplayUnits[1]; } // zooming, panning, restriction of both void mitk::DisplayGeometry::SetConstrainZoomingAndPanning(bool constrain) { m_ConstrainZoomingAndPanning = constrain; if (m_ConstrainZoomingAndPanning) { this->RefitVisibleRect(); } } bool mitk::DisplayGeometry::GetConstrainZommingAndPanning() const { return m_ConstrainZoomingAndPanning; } bool mitk::DisplayGeometry::SetScaleFactor(ScalarType mmPerDisplayUnit) { if(mmPerDisplayUnit<0.0001) { mmPerDisplayUnit=0.0001; } m_ScaleFactorMMPerDisplayUnit = mmPerDisplayUnit; assert(m_ScaleFactorMMPerDisplayUnit < ScalarTypeNumericTraits::infinity()); DisplayToWorld(m_SizeInDisplayUnits, m_SizeInMM); return !this->RefitVisibleRect(); } mitk::ScalarType mitk::DisplayGeometry::GetScaleFactorMMPerDisplayUnit() const { return m_ScaleFactorMMPerDisplayUnit; } // Zooms with a factor (1.0=identity) around the specified center in display units bool mitk::DisplayGeometry::Zoom(ScalarType factor, const Point2D& centerInDisplayUnits) { assert(factor > 0); if ( SetScaleFactor(m_ScaleFactorMMPerDisplayUnit/factor) ) { return SetOriginInMM(m_OriginInMM-centerInDisplayUnits.GetVectorFromOrigin()*(1-factor)*m_ScaleFactorMMPerDisplayUnit); } else { return false; } } // Zooms with a factor (1.0=identity) around the specified center, but tries (if its within view contraints) to match the center in display units with the center in world coordinates. bool mitk::DisplayGeometry::ZoomWithFixedWorldCoordinates(ScalarType factor, const Point2D& focusDisplayUnits, const Point2D& focusUnitsInMM ) { assert(factor > 0); SetScaleFactor(m_ScaleFactorMMPerDisplayUnit/factor); SetOriginInMM(focusUnitsInMM.GetVectorFromOrigin()-focusDisplayUnits.GetVectorFromOrigin()*m_ScaleFactorMMPerDisplayUnit); return true; } bool mitk::DisplayGeometry::MoveBy(const Vector2D& shiftInDisplayUnits) { SetOriginInMM(m_OriginInMM+shiftInDisplayUnits*m_ScaleFactorMMPerDisplayUnit); Modified(); return !this->RefitVisibleRect(); } void mitk::DisplayGeometry::Fit() { if((m_WorldGeometry.IsNull()) || (m_WorldGeometry->IsValid() == false)) return; /// \FIXME: try to remove all the casts int width=(int)m_SizeInDisplayUnits[0]; int height=(int)m_SizeInDisplayUnits[1]; ScalarType w = width; ScalarType h = height; const ScalarType& widthInMM = m_WorldGeometry->GetParametricExtentInMM(0); const ScalarType& heightInMM = m_WorldGeometry->GetParametricExtentInMM(1); ScalarType aspRatio=((ScalarType)widthInMM)/heightInMM; ScalarType x = (ScalarType)w/widthInMM; ScalarType y = (ScalarType)h/heightInMM; if (x > y) { w = (int) (aspRatio*h); } else { h = (int) (w/aspRatio); } if(w>0) { SetScaleFactor(widthInMM/w); } Vector2D origin_display; origin_display[0]=-(width-w)/2.0; origin_display[1]=-(height-h)/2.0; SetOriginInMM(origin_display*m_ScaleFactorMMPerDisplayUnit); this->RefitVisibleRect(); Modified(); } // conversion methods void mitk::DisplayGeometry::DisplayToWorld(const Point2D &pt_display, Point2D &pt_mm) const { pt_mm[0]=m_ScaleFactorMMPerDisplayUnit*pt_display[0]+m_OriginInMM[0]; pt_mm[1]=m_ScaleFactorMMPerDisplayUnit*pt_display[1]+m_OriginInMM[1]; } void mitk::DisplayGeometry::WorldToDisplay(const Point2D &pt_mm, Point2D &pt_display) const { pt_display[0]=(pt_mm[0]-m_OriginInMM[0])*(1.0/m_ScaleFactorMMPerDisplayUnit); pt_display[1]=(pt_mm[1]-m_OriginInMM[1])*(1.0/m_ScaleFactorMMPerDisplayUnit); } void mitk::DisplayGeometry::DisplayToWorld(const Vector2D &vec_display, Vector2D &vec_mm) const { vec_mm=vec_display*m_ScaleFactorMMPerDisplayUnit; } void mitk::DisplayGeometry::WorldToDisplay(const Vector2D &vec_mm, Vector2D &vec_display) const { vec_display=vec_mm*(1.0/m_ScaleFactorMMPerDisplayUnit); } void mitk::DisplayGeometry::ULDisplayToMM(const Point2D &pt_ULdisplay, Point2D &pt_mm) const { ULDisplayToDisplay(pt_ULdisplay, pt_mm); DisplayToWorld(pt_mm, pt_mm); } void mitk::DisplayGeometry::MMToULDisplay(const Point2D &pt_mm, Point2D &pt_ULdisplay) const { WorldToDisplay(pt_mm, pt_ULdisplay); DisplayToULDisplay(pt_ULdisplay, pt_ULdisplay); } void mitk::DisplayGeometry::ULDisplayToMM(const Vector2D &vec_ULdisplay, Vector2D &vec_mm) const { ULDisplayToDisplay(vec_ULdisplay, vec_mm); DisplayToWorld(vec_mm, vec_mm); } void mitk::DisplayGeometry::MMToULDisplay(const Vector2D &vec_mm, Vector2D &vec_ULdisplay) const { WorldToDisplay(vec_mm, vec_ULdisplay); DisplayToULDisplay(vec_ULdisplay, vec_ULdisplay); } void mitk::DisplayGeometry::ULDisplayToDisplay(const Point2D &pt_ULdisplay, Point2D &pt_display) const { pt_display[0]=pt_ULdisplay[0]; pt_display[1]=GetDisplayHeight()-pt_ULdisplay[1]; } void mitk::DisplayGeometry::DisplayToULDisplay(const Point2D &pt_display, Point2D &pt_ULdisplay) const { ULDisplayToDisplay(pt_display, pt_ULdisplay); } void mitk::DisplayGeometry::ULDisplayToDisplay(const Vector2D &vec_ULdisplay, Vector2D &vec_display) const { vec_display[0]= vec_ULdisplay[0]; vec_display[1]=-vec_ULdisplay[1]; } void mitk::DisplayGeometry::DisplayToULDisplay(const Vector2D &vec_display, Vector2D &vec_ULdisplay) const { ULDisplayToDisplay(vec_display, vec_ULdisplay); } bool mitk::DisplayGeometry::Project(const Point3D &pt3d_mm, Point3D &projectedPt3d_mm) const { if(m_WorldGeometry.IsNotNull()) { return m_WorldGeometry->Project(pt3d_mm, projectedPt3d_mm); } else { return false; } } bool mitk::DisplayGeometry::Project(const Point3D & atPt3d_mm, const Vector3D &vec3d_mm, Vector3D &projectedVec3d_mm) const { if(m_WorldGeometry.IsNotNull()) { return m_WorldGeometry->Project(atPt3d_mm, vec3d_mm, projectedVec3d_mm); } else { return false; } } bool mitk::DisplayGeometry::Project(const Vector3D &vec3d_mm, Vector3D &projectedVec3d_mm) const { if(m_WorldGeometry.IsNotNull()) { return m_WorldGeometry->Project(vec3d_mm, projectedVec3d_mm); } else { return false; } } bool mitk::DisplayGeometry::Map(const Point3D &pt3d_mm, Point2D &pt2d_mm) const { if(m_WorldGeometry.IsNotNull()) { return m_WorldGeometry->Map(pt3d_mm, pt2d_mm); } else { return false; } } void mitk::DisplayGeometry::Map(const Point2D &pt2d_mm, Point3D &pt3d_mm) const { if(m_WorldGeometry.IsNull()) return; m_WorldGeometry->Map(pt2d_mm, pt3d_mm); } bool mitk::DisplayGeometry::Map(const Point3D & atPt3d_mm, const Vector3D &vec3d_mm, Vector2D &vec2d_mm) const { if(m_WorldGeometry.IsNotNull()) { return m_WorldGeometry->Map(atPt3d_mm, vec3d_mm, vec2d_mm); } else { return false; } } void mitk::DisplayGeometry::Map(const Point2D & atPt2d_mm, const Vector2D &vec2d_mm, Vector3D &vec3d_mm) const { if(m_WorldGeometry.IsNull()) return; m_WorldGeometry->Map(atPt2d_mm, vec2d_mm, vec3d_mm); } // protected methods mitk::DisplayGeometry::DisplayGeometry() :m_ScaleFactorMMPerDisplayUnit(1.0) ,m_WorldGeometry(NULL) ,m_ConstrainZoomingAndPanning(true) ,m_MaxWorldViewPercentage(1.0) ,m_MinWorldViewPercentage(0.1) { m_OriginInMM.Fill(0.0); m_OriginInDisplayUnits.Fill(0.0); m_SizeInMM.Fill(1.0); m_SizeInDisplayUnits.Fill(10.0); } mitk::DisplayGeometry::~DisplayGeometry() { } bool mitk::DisplayGeometry::RefitVisibleRect() { // do nothing if not asked to if (!m_ConstrainZoomingAndPanning) return false; // don't allow recursion (need to be fixed, singleton) static bool inRecalculate = false; if (inRecalculate) return false; inRecalculate = true; // rename some basic measures of the current viewport and world geometry (MM = milimeters Px = Pixels = display units) float displayXMM = m_OriginInMM[0]; float displayYMM = m_OriginInMM[1]; float displayWidthPx = m_SizeInDisplayUnits[0]; float displayHeightPx = m_SizeInDisplayUnits[1]; float displayWidthMM = m_SizeInDisplayUnits[0] * m_ScaleFactorMMPerDisplayUnit; float displayHeightMM = m_SizeInDisplayUnits[1] * m_ScaleFactorMMPerDisplayUnit; float worldWidthMM = m_WorldGeometry->GetParametricExtentInMM(0); float worldHeightMM = m_WorldGeometry->GetParametricExtentInMM(1); // reserve variables for the correction logic to save a corrected origin and zoom factor Vector2D newOrigin = m_OriginInMM; bool correctPanning = false; float newScaleFactor = m_ScaleFactorMMPerDisplayUnit; bool correctZooming = false; // start of the correction logic // zoom to big means: // at a given percentage of the world's width/height should be visible. Otherwise // the whole screen could show only one pixel // // zoom to small means: // zooming out should be limited at the point where the smaller of the world's sides is completely visible bool zoomXtooSmall = displayWidthPx * m_ScaleFactorMMPerDisplayUnit > m_MaxWorldViewPercentage * worldWidthMM; bool zoomXtooBig = displayWidthPx * m_ScaleFactorMMPerDisplayUnit < m_MinWorldViewPercentage * worldWidthMM; bool zoomYtooSmall = displayHeightPx * m_ScaleFactorMMPerDisplayUnit > m_MaxWorldViewPercentage * worldHeightMM; bool zoomYtooBig = displayHeightPx * m_ScaleFactorMMPerDisplayUnit < m_MinWorldViewPercentage * worldHeightMM; // constrain zooming in both direction if ( zoomXtooBig && zoomYtooBig) { double fx = worldWidthMM * m_MinWorldViewPercentage / displayWidthPx; double fy = worldHeightMM * m_MinWorldViewPercentage / displayHeightPx; newScaleFactor = fx < fy ? fx : fy; correctZooming = true; } // constrain zooming in x direction else if ( zoomXtooBig ) { newScaleFactor = worldWidthMM * m_MinWorldViewPercentage / displayWidthPx; correctZooming = true; } // constrain zooming in y direction else if ( zoomYtooBig ) { newScaleFactor = worldHeightMM * m_MinWorldViewPercentage / displayHeightPx; correctZooming = true; } // constrain zooming out // we stop zooming out at these situations: // // *** display // --- image // // ********************** // * * x side maxed out // * * // *--------------------* // *| |* // *| |* // *--------------------* // * * // * * // * * // ********************** // // ********************** // * |------| * y side maxed out // * | | * // * | | * // * | | * // * | | * // * | | * // * | | * // * | | * // * |------| * // ********************** // // In both situations we center the not-maxed out direction // if ( zoomXtooSmall && zoomYtooSmall ) { // determine and set the bigger scale factor float fx = worldWidthMM * m_MaxWorldViewPercentage / displayWidthPx; float fy = worldHeightMM * m_MaxWorldViewPercentage / displayHeightPx; newScaleFactor = fx > fy ? fx : fy; correctZooming = true; } // actually execute correction if (correctZooming) { SetScaleFactor(newScaleFactor); } displayWidthMM = m_SizeInDisplayUnits[0] * m_ScaleFactorMMPerDisplayUnit; displayHeightMM = m_SizeInDisplayUnits[1] * m_ScaleFactorMMPerDisplayUnit; // constrain panning if(worldWidthMM center x newOrigin[0] = (worldWidthMM - displayWidthMM) / 2.0; correctPanning = true; } else { // make sure left display border inside our world if (displayXMM < 0) { newOrigin[0] = 0; correctPanning = true; } // make sure right display border inside our world else if (displayXMM + displayWidthMM > worldWidthMM) { newOrigin[0] = worldWidthMM - displayWidthMM; correctPanning = true; } } if (worldHeightMM center y newOrigin[1] = (worldHeightMM - displayHeightMM) / 2.0; correctPanning = true; } else { // make sure top display border inside our world if (displayYMM + displayHeightMM > worldHeightMM) { newOrigin[1] = worldHeightMM - displayHeightMM; correctPanning = true; } // make sure bottom display border inside our world else if (displayYMM < 0) { newOrigin[1] = 0; correctPanning = true; } } if (correctPanning) { SetOriginInMM( newOrigin ); } inRecalculate = false; if ( correctPanning || correctZooming ) { Modified(); } // return true if any correction has been made return correctPanning || correctZooming; } void mitk::DisplayGeometry::PrintSelf(std::ostream& os, itk::Indent indent) const { if(m_WorldGeometry.IsNull()) { os << indent << " WorldGeometry: " << "NULL" << std::endl; } else { m_WorldGeometry->Print(os, indent); os << indent << " OriginInMM: " << m_OriginInMM << std::endl; os << indent << " OriginInDisplayUnits: " << m_OriginInDisplayUnits << std::endl; os << indent << " SizeInMM: " << m_SizeInMM << std::endl; os << indent << " SizeInDisplayUnits: " << m_SizeInDisplayUnits << std::endl; os << indent << " ScaleFactorMMPerDisplayUni: " << m_ScaleFactorMMPerDisplayUnit << std::endl; } Superclass::PrintSelf(os,indent); } diff --git a/Core/Code/DataManagement/mitkDisplayGeometry.h b/Core/Code/DataManagement/mitkDisplayGeometry.h index 1c9275c77b..fdae96e680 100644 --- a/Core/Code/DataManagement/mitkDisplayGeometry.h +++ b/Core/Code/DataManagement/mitkDisplayGeometry.h @@ -1,227 +1,227 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDisplayGeometry_h #define mitkDisplayGeometry_h #include "mitkGeometry2D.h" namespace mitk { /** \brief Describes the geometry on the display/screen for 2D display. The main purpose of this class is to convert between display coordinates (in display-units) and world coordinates (in mm). DisplayGeometry depends on the size of the display area (widget width and height, m_SizeInDisplayUnits) and on a Geometry2D (m_WoldGeometry). It represents a recangular view on this world-geometry. E.g., you can tell the DisplayGeometry to fit the world-geometry in the display area by calling Fit(). Provides methods for zooming and panning. Zooming and panning can be restricted within reasonable bounds by setting the ConstrainZoomingAndPanning flag. In these cases you can re-define what bounds you accept as "reasonable" by calling \warning \em Units refers to the units of the underlying world-geometry. Take care, whether these are really the units you want to convert to. E.g., when you want to convert a point \a pt_display (which is 2D) given in display coordinates into a point in units of a BaseData-object @a datum (the requested point is 3D!), use \code displaygeometry->DisplayToWorld(pt_display, pt2d_mm); displaygeometry->Map(pt2d_mm, pt3d_mm); datum->GetGeometry()->WorldToIndex(pt3d_mm, pt3d_datum_units); \endcode Even, if you want to convert the 2D point \a pt_display into a 2D point in units on a certain 2D geometry \a certaingeometry, it is safer to use \code displaygeometry->DisplayToWorld(pt_display, pt_mm); certaingeometry->WorldToIndex(pt_mm, pt_certain_geometry_units); \endcode unless you can be sure that the underlying geometry of \a displaygeometry is really the \a certaingeometry. \ingroup Geometry */ class MITK_CORE_EXPORT DisplayGeometry : public Geometry2D { public: mitkClassMacro(DisplayGeometry,Geometry2D); /// Method for creation through the object factory. itkNewMacro(Self); /// \brief duplicates the geometry, NOT useful for this sub-class virtual itk::LightObject::Pointer InternalClone() const; /// \return this objects modified time. virtual unsigned long GetMTime() const; virtual const TimeBounds& GetTimeBounds() const; // size definition methods virtual void SetWorldGeometry(const Geometry2D* aWorldGeometry); itkGetConstObjectMacro(WorldGeometry, Geometry2D); /// \return if new origin was within accepted limits virtual bool SetOriginInMM(const Vector2D& origin_mm); virtual Vector2D GetOriginInMM() const; virtual Vector2D GetOriginInDisplayUnits() const; /** \brief Set the size of the display in display units. This method must be called every time the display is resized (normally, the GUI-toolkit informs about resizing). \param keepDisplayedRegion: if \a true (the default), the displayed contents is zoomed/shrinked so that the displayed region is (approximately) the same as before: The point at the center will be kept at the center and the length of the diagonal of the displayed region \em in \em units will also be kept. When the aspect ration changes, the displayed region includes the old displayed region, but cannot be exaclty the same. */ virtual void SetSizeInDisplayUnits(unsigned int width, unsigned int height, bool keepDisplayedRegion=true); virtual Vector2D GetSizeInDisplayUnits() const; virtual Vector2D GetSizeInMM() const; unsigned int GetDisplayWidth() const; unsigned int GetDisplayHeight() const; // zooming, panning, restriction of both virtual void SetConstrainZoomingAndPanning(bool constrain); virtual bool GetConstrainZommingAndPanning() const; /// what percentage of the world should be visible at maximum zoom out (default 1.0, i.e. 100% of width or height) itkGetMacro(MaxWorldViewPercentage, float); itkSetMacro(MaxWorldViewPercentage, float); /// what percentage of the world should be visible at maximum zoom in (default 0.1, i.e. 10% of width or height) itkGetMacro(MinWorldViewPercentage, float); itkSetMacro(MinWorldViewPercentage, float); virtual bool SetScaleFactor(ScalarType mmPerDisplayUnit); ScalarType GetScaleFactorMMPerDisplayUnit() const; /** * \brief Zooms with a factor (1.0=identity) to/from the specified center in display units * \return true if zoom request was within accepted limits */ virtual bool Zoom(ScalarType factor, const Point2D& centerInDisplayUnits); /** * \brief Zooms with a factor (1.0=identity) to/from the specified center, trying to preserve the center of zoom in world coordiantes * * Same zoom as mentioned above but tries (if it's within view contraints) to match the center in display units with the center in world coordinates. * * \return true if zoom request was within accepted limits */ virtual bool ZoomWithFixedWorldCoordinates(ScalarType factor, const Point2D& focusDisplayUnits, const Point2D& focusUnitsInMM ); // \return true if move request was within accepted limits virtual bool MoveBy(const Vector2D& shiftInDisplayUnits); // \brief align display with world, make world completely visible virtual void Fit(); // conversion methods virtual void DisplayToWorld(const Point2D &pt_display, Point2D &pt_mm) const; virtual void WorldToDisplay(const Point2D &pt_mm, Point2D &pt_display) const; virtual void DisplayToWorld(const Vector2D &vec_display, Vector2D &vec_mm) const; virtual void WorldToDisplay(const Vector2D &vec_mm, Vector2D &vec_display) const; virtual void ULDisplayToMM(const Point2D &pt_ULdisplay, Point2D &pt_mm) const; virtual void MMToULDisplay(const Point2D &pt_mm, Point2D &pt_ULdisplay) const; virtual void ULDisplayToMM(const Vector2D &vec_ULdisplay, Vector2D &vec_mm) const; virtual void MMToULDisplay(const Vector2D &vec_mm, Vector2D &vec_ULdisplay) const; virtual void ULDisplayToDisplay(const Point2D &pt_ULdisplay, Point2D &pt_display) const; virtual void DisplayToULDisplay(const Point2D &pt_display, Point2D &pt_ULdisplay) const; virtual void ULDisplayToDisplay(const Vector2D &vec_ULdisplay, Vector2D &vec_display) const; virtual void DisplayToULDisplay(const Vector2D &vec_display, Vector2D &vec_ULdisplay) const; /** * \brief projects the given point onto current 2D world geometry plane */ virtual bool Project(const Point3D &pt3d_mm, Point3D &projectedPt3d_mm) const; /** * \brief projects the given vector onto current 2D world geometry plane. * \warning DEPRECATED, please use Project(const Vector3D &vec3d_mm, Vector3D &projectedVec3d_mm) instead */ virtual bool Project(const Point3D & atPt3d_mm, const Vector3D &vec3d_mm, Vector3D &projectedVec3d_mm) const; /** * \brief projects the given vector onto current 2D world geometry plane */ virtual bool Project(const Vector3D &vec3d_mm, Vector3D &projectedVec3d_mm) const; virtual bool Map(const Point3D &pt3d_mm, Point2D &pt2d_mm) const; virtual void Map(const Point2D &pt2d_mm, Point3D &pt3d_mm) const; virtual bool Map(const Point3D & atPt3d_mm, const Vector3D &vec3d_mm, Vector2D &vec2d_mm) const; virtual void Map(const Point2D & atPt2d_mm, const Vector2D &vec2d_mm, Vector3D &vec3d_mm) const; + virtual bool IsValid() const; + protected: DisplayGeometry(); virtual ~DisplayGeometry(); /** \brief Called after zooming/panning to restrict these operations to sensible measures. \return true if a correction in either zooming or panning was made Enforces a couple of constraints on the relation of the current viewport and the current world geometry. The basic logic in this lengthy method is:
  1. Make display region big enough (in case of too large zoom factors)
  2. Make display region small enough (so that the image cannot be scaled into a single screen pixel
  3. Correct panning for each border (left, right, bottom, top)
The little more complicated implementation is illustrated in the code itself. */ virtual bool RefitVisibleRect(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; Vector2D m_OriginInMM; Vector2D m_OriginInDisplayUnits; ScalarType m_ScaleFactorMMPerDisplayUnit; Vector2D m_SizeInMM; Vector2D m_SizeInDisplayUnits; Geometry2D::ConstPointer m_WorldGeometry; bool m_ConstrainZoomingAndPanning; float m_MaxWorldViewPercentage; float m_MinWorldViewPercentage; - - virtual bool InternPostIsValid() const; }; } // namespace #endif // include guard diff --git a/Core/Code/DataManagement/mitkGeometry2D.cpp b/Core/Code/DataManagement/mitkGeometry2D.cpp index 67d7668663..52f633b0c9 100644 --- a/Core/Code/DataManagement/mitkGeometry2D.cpp +++ b/Core/Code/DataManagement/mitkGeometry2D.cpp @@ -1,264 +1,264 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkGeometry2D.h" #include mitk::Geometry2D::Geometry2D() : m_ScaleFactorMMPerUnitX( 1.0 ), m_ScaleFactorMMPerUnitY( 1.0 ), m_ReferenceGeometry( NULL ) { } mitk::Geometry2D::Geometry2D(const Geometry2D& other) : Geometry3D(other), m_ScaleFactorMMPerUnitX( other.m_ScaleFactorMMPerUnitX), m_ScaleFactorMMPerUnitY( other.m_ScaleFactorMMPerUnitY), m_ReferenceGeometry( other.m_ReferenceGeometry ) { } mitk::Geometry2D::~Geometry2D() { } void mitk::Geometry2D::InternPostSetIndexToWorldTransform( mitk::AffineTransform3D* transform) { m_ScaleFactorMMPerUnitX=GetExtentInMM(0)/GetExtent(0); m_ScaleFactorMMPerUnitY=GetExtentInMM(1)/GetExtent(1); assert(m_ScaleFactorMMPerUnitXIsBoundingBoxNull()==false); Point3D pt3d_units; BackTransform(pt3d_mm, pt3d_units); pt2d_mm[0]=pt3d_units[0]*m_ScaleFactorMMPerUnitX; pt2d_mm[1]=pt3d_units[1]*m_ScaleFactorMMPerUnitY; pt3d_units[2]=0; - return const_cast(m_BoundingBox.GetPointer())->IsInside(pt3d_units); + return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } void mitk::Geometry2D::Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const { Point3D pt3d_units; pt3d_units[0]=pt2d_mm[0]/m_ScaleFactorMMPerUnitX; pt3d_units[1]=pt2d_mm[1]/m_ScaleFactorMMPerUnitY; pt3d_units[2]=0; pt3d_mm = GetParametricTransform()->TransformPoint(pt3d_units); } void mitk::Geometry2D::IndexToWorld( const mitk::Point2D &/*pt_units*/, mitk::Point2D &/*pt_mm*/) const { itkExceptionMacro(<< "No general transform possible (only affine) ==> no general" \ " IndexToWorld(const mitk::Point2D &pt_mm, mitk::Point2D &pt_units)" \ " possible. Has to be implemented in sub-class."); } void mitk::Geometry2D::WorldToIndex( const mitk::Point2D &/*pt_mm*/, mitk::Point2D &/*pt_units*/) const { itkExceptionMacro(<< "No general back transform possible (only affine) ==> no general" \ " WorldToIndex(const mitk::Point2D &pt_mm, mitk::Point2D &pt_units)" \ " possible. Has to be implemented in sub-class."); } void mitk::Geometry2D::IndexToWorld(const mitk::Point2D &/*atPt2d_units*/, const mitk::Vector2D &/*vec_units*/, mitk::Vector2D &/*vec_mm*/) const { itkExceptionMacro(<< "No general transform possible (only affine) ==> no general" \ " IndexToWorld(const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units)" \ " possible. Has to be implemented in sub-class."); } void mitk::Geometry2D::WorldToIndex(const mitk::Point2D &/*atPt2d_mm*/, const mitk::Vector2D &/*vec_mm*/, mitk::Vector2D &/*vec_units*/) const { itkExceptionMacro(<< "No general back transform possible (only affine) ==> no general" \ " WorldToIndex(const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units)" \ " possible. Has to be implemented in sub-class."); } void mitk::Geometry2D::SetSizeInUnits(mitk::ScalarType width, mitk::ScalarType height) { ScalarType bounds[6]={0, width, 0, height, 0, 1}; ScalarType extent, newextentInMM; if(GetExtent(0)>0) { extent = GetExtent(0); if(width>extent) newextentInMM = GetExtentInMM(0)/width*extent; else newextentInMM = GetExtentInMM(0)*extent/width; SetExtentInMM(0, newextentInMM); } if(GetExtent(1)>0) { extent = GetExtent(1); if(width>extent) newextentInMM = GetExtentInMM(1)/height*extent; else newextentInMM = GetExtentInMM(1)*extent/height; SetExtentInMM(1, newextentInMM); } SetBounds(bounds); } bool mitk::Geometry2D::Project( const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const { - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); Point3D pt3d_units; BackTransform(pt3d_mm, pt3d_units); pt3d_units[2] = 0; projectedPt3d_mm = GetParametricTransform()->TransformPoint(pt3d_units); - return const_cast(m_BoundingBox.GetPointer())->IsInside(pt3d_units); + return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool mitk::Geometry2D::Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; BackTransform(vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetParametricTransform()->TransformVector(vec3d_units); return true; } bool mitk::Geometry2D::Project(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { MITK_WARN << "Deprecated function! Call Project(vec3D,vec3D) instead."; - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; BackTransform(atPt3d_mm, vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetParametricTransform()->TransformVector(vec3d_units); Point3D pt3d_units; BackTransform(atPt3d_mm, pt3d_units); - return const_cast(m_BoundingBox.GetPointer())->IsInside(pt3d_units); + return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool mitk::Geometry2D::Map(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector2D &vec2d_mm) const { Point2D pt2d_mm_start, pt2d_mm_end; Point3D pt3d_mm_end; bool inside=Map(atPt3d_mm, pt2d_mm_start); pt3d_mm_end = atPt3d_mm+vec3d_mm; inside&=Map(pt3d_mm_end, pt2d_mm_end); vec2d_mm=pt2d_mm_end-pt2d_mm_start; return inside; } void mitk::Geometry2D::Map(const mitk::Point2D &/*atPt2d_mm*/, const mitk::Vector2D &/*vec2d_mm*/, mitk::Vector3D &/*vec3d_mm*/) const { //@todo implement parallel to the other Map method! assert(false); } mitk::ScalarType mitk::Geometry2D::SignedDistance(const mitk::Point3D& pt3d_mm) const { Point3D projectedPoint; Project(pt3d_mm, projectedPoint); Vector3D direction = pt3d_mm-projectedPoint; ScalarType distance = direction.GetNorm(); if(IsAbove(pt3d_mm) == false) distance*=-1.0; return distance; } bool mitk::Geometry2D::IsAbove(const mitk::Point3D& pt3d_mm) const { Point3D pt3d_units; Geometry3D::WorldToIndex(pt3d_mm, pt3d_units); - return (pt3d_units[2] > m_BoundingBox->GetBounds()[4]); + return (pt3d_units[2] > this->GetBoundingBox()->GetBounds()[4]); } itk::LightObject::Pointer mitk::Geometry2D::InternalClone() const { Self::Pointer newGeometry = new Geometry2D(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void mitk::Geometry2D::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf(os,indent); os << indent << " ScaleFactorMMPerUnitX: " << m_ScaleFactorMMPerUnitX << std::endl; os << indent << " ScaleFactorMMPerUnitY: " << m_ScaleFactorMMPerUnitY << std::endl; } void mitk::Geometry2D::SetReferenceGeometry( mitk::Geometry3D *geometry ) { m_ReferenceGeometry = geometry; } mitk::Geometry3D * mitk::Geometry2D::GetReferenceGeometry() const { return m_ReferenceGeometry; } bool mitk::Geometry2D::HasReferenceGeometry() const { return ( m_ReferenceGeometry != NULL ); } diff --git a/Core/Code/DataManagement/mitkGeometry3D.cpp b/Core/Code/DataManagement/mitkGeometry3D.cpp index aca3f37727..3f64899f8b 100644 --- a/Core/Code/DataManagement/mitkGeometry3D.cpp +++ b/Core/Code/DataManagement/mitkGeometry3D.cpp @@ -1,328 +1,328 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include "mitkGeometry3D.h" #include "mitkRotationOperation.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkApplyTransformMatrixOperation.h" #include "mitkPointOperation.h" #include "mitkInteractionConst.h" #include #include #include "mitkMatrixConvert.h" // Standard constructor for the New() macro. Sets the geometry to 3 dimensions mitk::Geometry3D::Geometry3D() : m_ImageGeometry(false) { } mitk::Geometry3D::Geometry3D(const Geometry3D& other) : BaseGeometry(other), m_ImageGeometry(other.m_ImageGeometry), m_ParametricBoundingBox(other.m_ParametricBoundingBox) { if (other.m_ParametricBoundingBox.IsNotNull()) { m_ParametricBoundingBox = other.m_ParametricBoundingBox->DeepCopy(); } } mitk::Geometry3D::~Geometry3D() { } void mitk::Geometry3D::SetParametricBounds(const BoundingBox::BoundsArrayType& bounds) { m_ParametricBoundingBox = BoundingBoxType::New(); BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New(); BoundingBoxType::PointType p; BoundingBoxType::PointIdentifier pointid; for(pointid=0; pointid<2;++pointid) { unsigned int i; - for(i=0; iInsertElement(pointid, p); } m_ParametricBoundingBox->SetPoints(pointscontainer); m_ParametricBoundingBox->ComputeBoundingBox(); this->Modified(); } itk::LightObject::Pointer mitk::Geometry3D::InternalClone() const { Self::Pointer newGeometry = new Self(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void mitk::Geometry3D::InternPostInitializeGeometry(Geometry3D * newGeometry) const { newGeometry->m_ImageGeometry = m_ImageGeometry; } void mitk::Geometry3D::PrintSelf(std::ostream& os, itk::Indent indent) const { os << indent << " IndexToWorldTransform: "; - if(m_IndexToWorldTransform.IsNull()) + if(this->IsIndexToWorldTransformNull()) os << "NULL" << std::endl; else { // from itk::MatrixOffsetTransformBase unsigned int i, j; os << std::endl; os << indent << "Matrix: " << std::endl; for (i = 0; i < 3; i++) { os << indent.GetNextIndent(); for (j = 0; j < 3; j++) { - os << m_IndexToWorldTransform->GetMatrix()[i][j] << " "; + os << this->GetIndexToWorldTransform()->GetMatrix()[i][j] << " "; } os << std::endl; } - os << indent << "Offset: " << m_IndexToWorldTransform->GetOffset() << std::endl; - os << indent << "Center: " << m_IndexToWorldTransform->GetCenter() << std::endl; - os << indent << "Translation: " << m_IndexToWorldTransform->GetTranslation() << std::endl; + os << indent << "Offset: " << this->GetIndexToWorldTransform()->GetOffset() << std::endl; + os << indent << "Center: " << this->GetIndexToWorldTransform()->GetCenter() << std::endl; + os << indent << "Translation: " << this->GetIndexToWorldTransform()->GetTranslation() << std::endl; os << indent << "Inverse: " << std::endl; for (i = 0; i < 3; i++) { os << indent.GetNextIndent(); for (j = 0; j < 3; j++) { - os << m_IndexToWorldTransform->GetInverseMatrix()[i][j] << " "; + os << this->GetIndexToWorldTransform()->GetInverseMatrix()[i][j] << " "; } os << std::endl; } // from itk::ScalableAffineTransform os << indent << "Scale : "; for (i = 0; i < 3; i++) { - os << m_IndexToWorldTransform->GetScale()[i] << " "; + os << this->GetIndexToWorldTransform()->GetScale()[i] << " "; } os << std::endl; } os << indent << " BoundingBox: "; - if(m_BoundingBox.IsNull()) + if(this->IsBoundingBoxNull()) os << "NULL" << std::endl; else { os << indent << "( "; for (unsigned int i=0; i<3; i++) { - os << m_BoundingBox->GetBounds()[2*i] << "," << m_BoundingBox->GetBounds()[2*i+1] << " "; + os << this->GetBoundingBox()->GetBounds()[2*i] << "," << this->GetBoundingBox()->GetBounds()[2*i+1] << " "; } os << " )" << std::endl; } - os << indent << " Origin: " << m_Origin << std::endl; + os << indent << " Origin: " << this->GetOrigin() << std::endl; os << indent << " ImageGeometry: " << m_ImageGeometry << std::endl; - os << indent << " Spacing: " << m_Spacing << std::endl; - os << indent << " TimeBounds: " << m_TimeBounds << std::endl; + os << indent << " Spacing: " << this->GetSpacing() << std::endl; + os << indent << " TimeBounds: " << this->GetTimeBounds() << std::endl; } mitk::Point3D mitk::Geometry3D::GetCornerPoint(int id) const { assert(id >= 0); - assert(m_BoundingBox.IsNotNull()); + assert(this->IsBoundingBoxNull()==false); - BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); + BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds(); Point3D cornerpoint; switch(id) { case 0: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[4]); break; case 1: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[5]); break; case 2: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[4]); break; case 3: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[5]); break; case 4: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[4]); break; case 5: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[5]); break; case 6: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[4]); break; case 7: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[5]); break; default: { itkExceptionMacro(<<"A cube only has 8 corners. These are labeled 0-7."); } } if(m_ImageGeometry) { // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the // bounding box. The bounding box itself is no image, so it is corner-based FillVector3D(cornerpoint, cornerpoint[0]-0.5, cornerpoint[1]-0.5, cornerpoint[2]-0.5); } - return m_IndexToWorldTransform->TransformPoint(cornerpoint); + return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint); } mitk::Point3D mitk::Geometry3D::GetCornerPoint(bool xFront, bool yFront, bool zFront) const { - assert(m_BoundingBox.IsNotNull()); - BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); + assert(this->IsBoundingBoxNull()==false); + BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds(); Point3D cornerpoint; cornerpoint[0] = (xFront ? bounds[0] : bounds[1]); cornerpoint[1] = (yFront ? bounds[2] : bounds[3]); cornerpoint[2] = (zFront ? bounds[4] : bounds[5]); if(m_ImageGeometry) { // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the // bounding box. The bounding box itself is no image, so it is corner-based FillVector3D(cornerpoint, cornerpoint[0]-0.5, cornerpoint[1]-0.5, cornerpoint[2]-0.5); } - return m_IndexToWorldTransform->TransformPoint(cornerpoint); + return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint); } void mitk::Geometry3D::ChangeImageGeometryConsideringOriginOffset( const bool isAnImageGeometry ) { // If Geometry is switched to ImageGeometry, you have to put an offset to the origin, because // imageGeometries origins are pixel-center-based // ... and remove the offset, if you switch an imageGeometry back to a normal geometry // For more information please see the Geometry documentation page if(m_ImageGeometry == isAnImageGeometry) return; const BoundingBox::BoundsArrayType& boundsarray = this->GetBoundingBox()->GetBounds(); Point3D originIndex; FillVector3D(originIndex, boundsarray[0], boundsarray[2], boundsarray[4]); if(isAnImageGeometry == true) FillVector3D( originIndex, originIndex[0] + 0.5, originIndex[1] + 0.5, originIndex[2] + 0.5 ); else FillVector3D( originIndex, originIndex[0] - 0.5, originIndex[1] - 0.5, originIndex[2] - 0.5 ); Point3D originWorld; originWorld = GetIndexToWorldTransform() ->TransformPoint( originIndex ); // instead could as well call IndexToWorld(originIndex,originWorld); SetOrigin(originWorld); this->SetImageGeometry(isAnImageGeometry); } bool mitk::Equal(const mitk::Geometry3D *leftHandSide, const mitk::Geometry3D *rightHandSide, ScalarType eps, bool verbose) { bool result = true; if( rightHandSide == NULL ) { if(verbose) MITK_INFO << "[( Geometry3D )] rightHandSide NULL."; return false; } if( leftHandSide == NULL) { if(verbose) MITK_INFO << "[( Geometry3D )] leftHandSide NULL."; return false; } //Compare spacings if( !mitk::Equal( leftHandSide->GetSpacing(), rightHandSide->GetSpacing(), eps ) ) { if(verbose) { MITK_INFO << "[( Geometry3D )] Spacing differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetSpacing() << " : leftHandSide is " << leftHandSide->GetSpacing() << " and tolerance is " << eps; } result = false; } //Compare Origins if( !mitk::Equal( leftHandSide->GetOrigin(), rightHandSide->GetOrigin(), eps ) ) { if(verbose) { MITK_INFO << "[( Geometry3D )] Origin differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetOrigin() << " : leftHandSide is " << leftHandSide->GetOrigin() << " and tolerance is " << eps; } result = false; } //Compare Axis and Extents for( unsigned int i=0; i<3; ++i) { if( !mitk::Equal( leftHandSide->GetAxisVector(i), rightHandSide->GetAxisVector(i), eps)) { if(verbose) { MITK_INFO << "[( Geometry3D )] AxisVector #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetAxisVector(i) << " : leftHandSide is " << leftHandSide->GetAxisVector(i) << " and tolerance is " << eps; } result = false; } if( !mitk::Equal( leftHandSide->GetExtent(i), rightHandSide->GetExtent(i), eps) ) { if(verbose) { MITK_INFO << "[( Geometry3D )] Extent #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetExtent(i) << " : leftHandSide is " << leftHandSide->GetExtent(i) << " and tolerance is " << eps; } result = false; } } //Compare ImageGeometry Flag if( rightHandSide->GetImageGeometry() != leftHandSide->GetImageGeometry() ) { if(verbose) { MITK_INFO << "[( Geometry3D )] GetImageGeometry is different."; MITK_INFO << "rightHandSide is " << rightHandSide->GetImageGeometry() << " : leftHandSide is " << leftHandSide->GetImageGeometry(); } result = false; } //Compare BoundingBoxes if( !mitk::Equal( leftHandSide->GetBoundingBox(), rightHandSide->GetBoundingBox(), eps, verbose) ) { result = false; } //Compare IndexToWorldTransform Matrix if( !mitk::Equal( leftHandSide->GetIndexToWorldTransform(), rightHandSide->GetIndexToWorldTransform(), eps, verbose) ) { result = false; } return result; } void mitk::Geometry3D::InternPostInitialize() { m_ImageGeometry = false; } diff --git a/Core/Code/DataManagement/mitkGeometry3D.h b/Core/Code/DataManagement/mitkGeometry3D.h index 11c3d47b15..b32085ad3f 100644 --- a/Core/Code/DataManagement/mitkGeometry3D.h +++ b/Core/Code/DataManagement/mitkGeometry3D.h @@ -1,299 +1,298 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD #define GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD #include #include #include "itkScalableAffineTransform.h" #include #include "mitkBaseGeometry.h" class vtkLinearTransform; namespace mitk { //##Documentation //## @brief Standard typedef for time-bounds typedef itk::FixedArray TimeBounds; typedef itk::FixedArray FixedArrayType; //##Documentation //## @brief Standard 3D-BoundingBox typedef //## //## Standard 3D-BoundingBox typedef to get rid of template arguments (3D, type). typedef itk::BoundingBox BoundingBox; //##Documentation //## @brief Describes the geometry of a data object //## //## At least, it can return the bounding box of the data object. //## //## The class holds //## \li a bounding box which is axes-parallel in intrinsic coordinates //## (often integer indices of pixels), to be accessed by //## GetBoundingBox() //## \li a transform to convert intrinsic coordinates into a //## world-coordinate system with coordinates in millimeters //## and milliseconds (all are floating point values), to //## be accessed by GetIndexToWorldTransform() //## \li a life span, i.e. a bounding box in time in ms (with //## start and end time), to be accessed by GetTimeBounds(). //## The default is minus infinity to plus infinity. //## //## Geometry3D and its sub-classes allow converting between //## intrinsic coordinates (called index or unit coordinates) //## and world-coordinates (called world or mm coordinates), //## e.g. WorldToIndex. //## In case you need integer index coordinates, provide an //## mitk::Index3D (or itk::Index) as target variable to //## WorldToIndex, otherwise you will get a continuous index //## (floating point values). //## //## An important sub-class is SlicedGeometry3D, which descibes //## data objects consisting of slices, e.g., objects of type Image. //## Conversions between world coordinates (in mm) and unit coordinates //## (e.g., pixels in the case of an Image) can be performed. //## //## For more information on related classes, see \ref Geometry. //## //## Geometry3D instances referring to an Image need a slightly //## different definition of corners, see SetImageGeometry. This //## is usualy automatically called by Image. //## //## Geometry3D have to be initialized in the method GenerateOutputInformation() //## of BaseProcess (or CopyInformation/ UpdateOutputInformation of BaseData, //## if possible, e.g., by analyzing pic tags in Image) subclasses. See also //## itk::ProcessObject::GenerateOutputInformation(), //## itk::DataObject::CopyInformation() and //## itk::DataObject::UpdateOutputInformation(). //## //## Rule: everything is in mm (ms) if not stated otherwise. //## @ingroup Geometry class MITK_CORE_EXPORT Geometry3D : public BaseGeometry { public: mitkClassMacro(Geometry3D, mitk::BaseGeometry); - //void testXYZ(){ int a=1; }; //xxxxxxxxxxxxxxxx typedef itk::QuaternionRigidTransform< ScalarType > QuaternionTransformType; typedef QuaternionTransformType::VnlQuaternionType VnlQuaternionType; /** Method for creation through the object factory. */ itkNewMacro(Self); mitkNewMacro1Param(Self,Self); //##Documentation //## @brief When switching from an Image Geometry to a normal Geometry (and the other way around), you have to change the origin as well (See Geometry Documentation)! This function will change the "isImageGeometry" bool flag and changes the origin respectively. virtual void ChangeImageGeometryConsideringOriginOffset( const bool isAnImageGeometry ); //##Documentation //## @brief Get the position of the corner number \a id (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. Point3D GetCornerPoint(int id) const; //##Documentation //## @brief Get the position of a corner (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. Point3D GetCornerPoint(bool xFront=true, bool yFront=true, bool zFront=true) const; //##Documentation //## @brief Is this an ImageGeometry? //## //## For more information, see SetImageGeometry itkGetConstMacro(ImageGeometry, bool); //##Documentation //## @brief Define that this Geometry3D is refering to an Image //## //## A geometry referring to an Image needs a slightly different //## definition of the position of the corners (see GetCornerPoint). //## The position of a voxel is defined by the position of its center. //## If we would use the origin (position of the (center of) the first //## voxel) as a corner and display this point, it would seem to be //## \em not at the corner but a bit within the image. Even worse for //## the opposite corner of the image: here the corner would appear //## outside the image (by half of the voxel diameter). Thus, we have //## to correct for this and to be able to do that, we need to know //## that the Geometry3D is referring to an Image. itkSetMacro(ImageGeometry, bool); itkBooleanMacro(ImageGeometry); //##Documentation //## @brief Test whether the point \a p ((continous!)index coordinates in units) is //## inside the bounding box bool IsIndexInside(const mitk::Point3D& index) const { bool inside = false; //if it is an image geometry, we need to convert the index to discrete values //this is done by applying the rounding function also used in WorldToIndex (see line 323) if (m_ImageGeometry) { mitk::Point3D discretIndex; discretIndex[0]=itk::Math::RoundHalfIntegerUp( index[0] ); discretIndex[1]=itk::Math::RoundHalfIntegerUp( index[1] ); discretIndex[2]=itk::Math::RoundHalfIntegerUp( index[2] ); - inside = m_BoundingBox->IsInside(discretIndex); + inside = this->GetBoundingBox()->IsInside(discretIndex); //we have to check if the index is at the upper border of each dimension, // because the boundingbox is not centerbased if (inside) { - const BoundingBox::BoundsArrayType& bounds = m_BoundingBox->GetBounds(); + const BoundingBox::BoundsArrayType& bounds = this->GetBoundingBox()->GetBounds(); if((discretIndex[0] == bounds[1]) || (discretIndex[1] == bounds[3]) || (discretIndex[2] == bounds[5])) inside = false; } } else - inside = m_BoundingBox->IsInside(index); + inside = this->GetBoundingBox()->IsInside(index); return inside; } //##Documentation //## @brief Convenience method for working with ITK indices template bool IsIndexInside(const itk::Index &index) const { int i, dim=index.GetIndexDimension(); Point3D pt_index; pt_index.Fill(0); for ( i = 0; i < dim; ++i ) { pt_index[i] = index[i]; } return IsIndexInside(pt_index); } //##Documentation //## @brief clones the geometry //## //## Overwrite in all sub-classes. //## Normally looks like: //## \code //## Self::Pointer newGeometry = new Self(*this); //## newGeometry->UnRegister(); //## return newGeometry.GetPointer(); //## \endcode virtual itk::LightObject::Pointer InternalClone() const; //Umzug: //##Documentation //## @brief Get the parametric bounding-box //## //## See AbstractTransformGeometry for an example usage of this. itkGetConstObjectMacro(ParametricBoundingBox, BoundingBox); //##Documentation //## @brief Get the parametric bounds //## //## See AbstractTransformGeometry for an example usage of this. const BoundingBox::BoundsArrayType& GetParametricBounds() const { assert(m_ParametricBoundingBox.IsNotNull()); return m_ParametricBoundingBox->GetBounds(); } //##Documentation //## @brief Get the parametric extent //## //## See AbstractTransformGeometry for an example usage of this. mitk::ScalarType GetParametricExtent(int direction) const { if (direction < 0 || direction>=3) mitkThrow() << "Invalid direction. Must be between either 0, 1 or 2. "; assert(m_ParametricBoundingBox.IsNotNull()); BoundingBoxType::BoundsArrayType bounds = m_ParametricBoundingBox->GetBounds(); return bounds[direction*2+1]-bounds[direction*2]; } //##Documentation //## @brief Get the parametric extent in mm //## //## See AbstractTransformGeometry for an example usage of this. virtual mitk::ScalarType GetParametricExtentInMM(int direction) const { return GetExtentInMM(direction); } //##Documentation //## @brief Get the parametric transform //## //## See AbstractTransformGeometry for an example usage of this. virtual const Transform3D* GetParametricTransform() const { - return m_IndexToWorldTransform; + return this->GetIndexToWorldTransform(); } protected: Geometry3D(); Geometry3D(const Geometry3D& other); virtual ~Geometry3D(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; //##Documentation //## @brief Set the parametric bounds //## //## Protected in this class, made public in some sub-classes, e.g., //## ExternAbstractTransformGeometry. virtual void SetParametricBounds(const BoundingBox::BoundsArrayType& bounds); mutable mitk::BoundingBox::Pointer m_ParametricBoundingBox; bool m_ImageGeometry; virtual void InternPostInitialize(); virtual void InternPostInitializeGeometry(Geometry3D* newGeometry) const; static const std::string INDEX_TO_OBJECT_TRANSFORM; static const std::string OBJECT_TO_NODE_TRANSFORM; static const std::string INDEX_TO_NODE_TRANSFORM; static const std::string INDEX_TO_WORLD_TRANSFORM; }; // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two geometries for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the spacing, origin, axisvectors, extents, the matrix of the * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag. * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * If you want to use different tolarance values for different parts of the geometry, feel free to use * the other comparison methods and write your own implementation of Equal. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::Geometry3D* leftHandSide, const mitk::Geometry3D* rightHandSide, ScalarType eps, bool verbose); } // namespace mitk #endif /* GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD */ diff --git a/Core/Code/DataManagement/mitkMatrixConvert.h b/Core/Code/DataManagement/mitkMatrixConvert.h index e7ff867d3d..b549a10fbc 100644 --- a/Core/Code/DataManagement/mitkMatrixConvert.h +++ b/Core/Code/DataManagement/mitkMatrixConvert.h @@ -1,151 +1,152 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKMATRIXCONVERT_H_HEADER_INCLUDED_C1EBD0AD #define MITKMATRIXCONVERT_H_HEADER_INCLUDED_C1EBD0AD #include "mitkGeometry3D.h" +#include "mitkBaseGeometry.h" #include "mitkItkMatrixHack.h" #include namespace mitk { template void TransferVtkMatrixToItkTransform(const vtkMatrix4x4* vtkmatrix, TTransformType * itkTransform) { if(itkTransform==NULL) return; typename TTransformType::MatrixType::InternalMatrixType& vnlMatrix = const_cast(itkTransform->GetMatrix().GetVnlMatrix()); for ( int i=0; i < 3; ++i) for( int j=0; j < 3; ++j ) vnlMatrix[i][j] = vtkmatrix->GetElement( i, j ); // *This* ensures m_MatrixMTime.Modified(), which is therewith not equal to // m_InverseMatrixMTime, thus a new inverse will be calculated (when // requested). static_cast*>(itkTransform)->MatrixChanged(); typename TTransformType::OffsetType offset; offset[0] = vtkmatrix->GetElement( 0, 3 ); offset[1] = vtkmatrix->GetElement( 1, 3 ); offset[2] = vtkmatrix->GetElement( 2, 3 ); itkTransform->SetOffset( offset ); } template void TransferItkTransformToVtkMatrix(const TTransformType * itkTransform, vtkMatrix4x4* vtkmatrix) { int i,j; for(i=0;i<3;++i) for(j=0;j<3;++j) vtkmatrix->SetElement(i, j, itkTransform->GetMatrix().GetVnlMatrix().get(i, j)); for(i=0;i<3;++i) vtkmatrix->SetElement(i, 3, itkTransform->GetOffset()[i]); for(i=0;i<3;++i) vtkmatrix->SetElement(3, i, 0.0); vtkmatrix->SetElement(3, 3, 1); } template void ConvertItkTransform(const TTransformType1* sourceTransform, TTransformType2* destTransform) { if((sourceTransform==NULL) || (destTransform==NULL)) return; // transfer offset const typename TTransformType1::OutputVectorType& sourceOffset = sourceTransform->GetOffset(); typename TTransformType2::OutputVectorType offset; offset[0] = sourceOffset[0]; offset[1] = sourceOffset[1]; offset[2] = sourceOffset[2]; destTransform->SetOffset( offset ); typename TTransformType1::MatrixType::InternalMatrixType& sourceVnlMatrix = const_cast(sourceTransform->GetMatrix().GetVnlMatrix()); //transfer matrix typename TTransformType2::MatrixType::InternalMatrixType& destVnlMatrix = const_cast(destTransform->GetMatrix().GetVnlMatrix()); for ( int i=0; i < 3; ++i) for( int j=0; j < 3; ++j ) destVnlMatrix[i][j] = sourceVnlMatrix[i][j]; // *This* ensures m_MatrixMTime.Modified(), which is therewith not equal to // m_InverseMatrixMTime, thus a new inverse will be calculated (when // requested). static_cast*>(destTransform)->MatrixChanged(); } template - void GetRotation(const mitk::Geometry3D * geometry, TMatrixType& itkmatrix) + void GetRotation(const mitk::BaseGeometry * geometry, TMatrixType& itkmatrix) { const mitk::Vector3D& spacing = geometry->GetSpacing(); - typename mitk::Geometry3D::TransformType::MatrixType::InternalMatrixType& geometryVnlMatrix = - const_cast(geometry->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix()); + typename mitk::BaseGeometry::TransformType::MatrixType::InternalMatrixType& geometryVnlMatrix = + const_cast(geometry->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix()); typename TMatrixType::InternalMatrixType& outputVnlMatrix = const_cast(itkmatrix.GetVnlMatrix()); for ( int i=0; i < 3; ++i) for( int j=0; j < 3; ++j ) outputVnlMatrix [i][j] = geometryVnlMatrix [i][j] / spacing[j]; } template void GetWorldToItkPhysicalTransform(const mitk::Geometry3D * geometry, TTransformType* itkTransform) { if(itkTransform==NULL) return; // get rotation matrix and offset from Geometry and transfer in TTransformType types typename TTransformType::MatrixType rotationMatrix; GetRotation(geometry, rotationMatrix); const typename mitk::Geometry3D::TransformType::OffsetType& geometryOffset = geometry->GetIndexToWorldTransform()->GetOffset(); vnl_vector vnlOffset(3); vnlOffset[0] = geometryOffset[0]; vnlOffset[1] = geometryOffset[1]; vnlOffset[2] = geometryOffset[2]; // do calculations typename TTransformType::MatrixType::InternalMatrixType inverseRotationVnlMatrix = rotationMatrix.GetTranspose(); vnlOffset -= inverseRotationVnlMatrix*vnlOffset; typename TTransformType::OutputVectorType offset;//vnl_vector offset; offset[0] = vnlOffset[0]; offset[1] = vnlOffset[1]; offset[2] = vnlOffset[2]; itkTransform->SetOffset( offset ); // copy in destination itkTransform typename TTransformType::MatrixType::InternalMatrixType& destVnlMatrix = const_cast(itkTransform->GetMatrix().GetVnlMatrix()); for ( int i=0; i < 3; ++i) for( int j=0; j < 3; ++j ) destVnlMatrix[i][j] = inverseRotationVnlMatrix[i][j]; // *This* ensures m_MatrixMTime.Modified(), which is therewith not equal to // m_InverseMatrixMTime, thus a new inverse will be calculated (when // requested). static_cast*>(itkTransform)->MatrixChanged(); } } #endif /* MITKMATRIXCONVERT_H_HEADER_INCLUDED_C1EBD0AD */ diff --git a/Core/Code/DataManagement/mitkPlaneGeometry.cpp b/Core/Code/DataManagement/mitkPlaneGeometry.cpp index c4caad38f1..08905e6425 100644 --- a/Core/Code/DataManagement/mitkPlaneGeometry.cpp +++ b/Core/Code/DataManagement/mitkPlaneGeometry.cpp @@ -1,727 +1,727 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPlaneGeometry.h" #include "mitkPlaneOperation.h" #include "mitkInteractionConst.h" #include "mitkLine.h" #include #include namespace mitk { PlaneGeometry::PlaneGeometry() { Initialize(); } PlaneGeometry::~PlaneGeometry() { } void PlaneGeometry::EnsurePerpendicularNormal(mitk::AffineTransform3D *transform) { //ensure row(2) of transform to be perpendicular to plane, keep length. VnlVector normal = vnl_cross_3d( transform->GetMatrix().GetVnlMatrix().get_column(0), transform->GetMatrix().GetVnlMatrix().get_column(1) ); normal.normalize(); ScalarType len = transform->GetMatrix() .GetVnlMatrix().get_column(2).two_norm(); if (len==0) len = 1; normal*=len; Matrix3D matrix = transform->GetMatrix(); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); } void PlaneGeometry::InternPreSetIndexToWorldTransform(mitk::AffineTransform3D *transform) { EnsurePerpendicularNormal(transform); } void PlaneGeometry::InternPreSetBounds(const BoundingBox::BoundsArrayType &bounds) { //currently the unit rectangle must be starting at the origin [0,0] assert(bounds[0]==0); assert(bounds[2]==0); //the unit rectangle must be two-dimensional assert(bounds[1]>0); assert(bounds[3]>0); } void PlaneGeometry::IndexToWorld( const Point2D &pt_units, Point2D &pt_mm ) const { pt_mm[0]=m_ScaleFactorMMPerUnitX*pt_units[0]; pt_mm[1]=m_ScaleFactorMMPerUnitY*pt_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D &pt_mm, Point2D &pt_units ) const { pt_units[0]=pt_mm[0]*(1.0/m_ScaleFactorMMPerUnitX); pt_units[1]=pt_mm[1]*(1.0/m_ScaleFactorMMPerUnitY); } void PlaneGeometry::IndexToWorld( const Point2D & /*atPt2d_units*/, const Vector2D &vec_units, Vector2D &vec_mm) const { MITK_WARN<<"Warning! Call of the deprecated function PlaneGeometry::IndexToWorld(point, vec, vec). Use PlaneGeometry::IndexToWorld(vec, vec) instead!"; this->IndexToWorld(vec_units, vec_mm); } void PlaneGeometry::IndexToWorld(const Vector2D &vec_units, Vector2D &vec_mm) const { vec_mm[0] = m_ScaleFactorMMPerUnitX * vec_units[0]; vec_mm[1] = m_ScaleFactorMMPerUnitY * vec_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D & /*atPt2d_mm*/, const Vector2D &vec_mm, Vector2D &vec_units) const { MITK_WARN<<"Warning! Call of the deprecated function PlaneGeometry::WorldToIndex(point, vec, vec). Use PlaneGeometry::WorldToIndex(vec, vec) instead!"; this->WorldToIndex(vec_mm, vec_units); } void PlaneGeometry::WorldToIndex( const Vector2D &vec_mm, Vector2D &vec_units) const { vec_units[0] = vec_mm[0] * ( 1.0 / m_ScaleFactorMMPerUnitX ); vec_units[1] = vec_mm[1] * ( 1.0 / m_ScaleFactorMMPerUnitY ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const Vector3D & spacing, PlaneGeometry::PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { AffineTransform3D::Pointer transform; transform = AffineTransform3D::New(); AffineTransform3D::MatrixType matrix; AffineTransform3D::MatrixType::InternalMatrixType &vnlmatrix = matrix.GetVnlMatrix(); vnlmatrix.set_identity(); vnlmatrix(0,0) = spacing[0]; vnlmatrix(1,1) = spacing[1]; vnlmatrix(2,2) = spacing[2]; transform->SetIdentity(); transform->SetMatrix(matrix); InitializeStandardPlane(width, height, transform.GetPointer(), planeorientation, zPosition, frontside, rotated); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const AffineTransform3D* transform, PlaneGeometry::PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { Superclass::Initialize(); //construct standard view Point3D origin; VnlVector rightDV(3), bottomDV(3); origin.Fill(0); int normalDirection; switch(planeorientation) { case Axial: if(frontside) { if(rotated==false) { FillVector3D(origin, 0, 0, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } else { FillVector3D(origin, width, height, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } else { if(rotated==false) { FillVector3D(origin, width, 0, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } else { FillVector3D(origin, 0, height, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } normalDirection = 2; break; case Frontal: if(frontside) { if(rotated==false) { FillVector3D(origin, 0, zPosition, 0); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, width, zPosition, height); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) { FillVector3D(origin, width, zPosition, 0); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, 0, zPosition, height); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 1; break; case Sagittal: if(frontside) { if(rotated==false) { FillVector3D(origin, zPosition, 0, 0); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, zPosition, width, height); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) { FillVector3D(origin, zPosition, width, 0); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, zPosition, 0, height); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 0; break; default: itkExceptionMacro("unknown PlaneOrientation"); } if ( transform != NULL ) { origin = transform->TransformPoint( origin ); rightDV = transform->TransformVector( rightDV ); bottomDV = transform->TransformVector( bottomDV ); } ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); if ( transform == NULL ) { this->SetMatrixByVectors( rightDV, bottomDV ); } else { this->SetMatrixByVectors( rightDV, bottomDV, transform->GetMatrix().GetVnlMatrix() .get_column(normalDirection).magnitude() ); } this->SetOrigin(origin); } void PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { this->SetReferenceGeometry( const_cast< Geometry3D * >( geometry3D ) ); ScalarType width, height; const BoundingBox::BoundsArrayType& boundsarray = geometry3D->GetBoundingBox()->GetBounds(); Vector3D originVector; FillVector3D(originVector, boundsarray[0], boundsarray[2], boundsarray[4]); if(geometry3D->GetImageGeometry()) { FillVector3D( originVector, originVector[0] - 0.5, originVector[1] - 0.5, originVector[2] - 0.5 ); } switch(planeorientation) { case Axial: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(1); break; case Frontal: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(2); break; case Sagittal: width = geometry3D->GetExtent(1); height = geometry3D->GetExtent(2); break; default: itkExceptionMacro("unknown PlaneOrientation"); } InitializeStandardPlane( width, height, geometry3D->GetIndexToWorldTransform(), planeorientation, zPosition, frontside, rotated ); ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); Point3D origin; originVector = geometry3D->GetIndexToWorldTransform() ->TransformVector( originVector ); origin = GetOrigin() + originVector; SetOrigin(origin); } void PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated ) { ScalarType zPosition; switch(planeorientation) { case Axial: zPosition = (top ? 0.5 : geometry3D->GetExtent(2)-1+0.5); break; case Frontal: zPosition = (top ? 0.5 : geometry3D->GetExtent(1)-1+0.5); break; case Sagittal: zPosition = (top ? 0.5 : geometry3D->GetExtent(0)-1+0.5); break; default: itkExceptionMacro("unknown PlaneOrientation"); } InitializeStandardPlane( geometry3D, planeorientation, zPosition, frontside, rotated ); } void PlaneGeometry::InitializeStandardPlane( const Vector3D &rightVector, const Vector3D &downVector, const Vector3D *spacing ) { InitializeStandardPlane( rightVector.GetVnlVector(), downVector.GetVnlVector(), spacing ); } void PlaneGeometry::InitializeStandardPlane( const VnlVector& rightVector, const VnlVector &downVector, const Vector3D *spacing ) { ScalarType width = rightVector.magnitude(); ScalarType height = downVector.magnitude(); InitializeStandardPlane( width, height, rightVector, downVector, spacing ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const Vector3D &rightVector, const Vector3D &downVector, const Vector3D *spacing ) { InitializeStandardPlane( width, height, rightVector.GetVnlVector(), downVector.GetVnlVector(), spacing ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const VnlVector &rightVector, const VnlVector &downVector, const Vector3D *spacing ) { assert(width > 0); assert(height > 0); VnlVector rightDV = rightVector; rightDV.normalize(); VnlVector downDV = downVector; downDV.normalize(); VnlVector normal = vnl_cross_3d(rightVector, downVector); normal.normalize(); if(spacing!=NULL) { rightDV *= (*spacing)[0]; downDV *= (*spacing)[1]; normal *= (*spacing)[2]; } AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightDV); matrix.GetVnlMatrix().set_column(1, downDV); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); - transform->SetOffset(m_IndexToWorldTransform->GetOffset()); + transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset()); ScalarType bounds[6] = { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); this->SetIndexToWorldTransform( transform ); } void PlaneGeometry::InitializePlane( const Point3D &origin, const Vector3D &normal ) { VnlVector rightVectorVnl(3), downVectorVnl; if( Equal( normal[1], 0.0f ) == false ) { FillVector3D( rightVectorVnl, 1.0f, -normal[0]/normal[1], 0.0f ); rightVectorVnl.normalize(); } else { FillVector3D( rightVectorVnl, 0.0f, 1.0f, 0.0f ); } downVectorVnl = vnl_cross_3d( normal.GetVnlVector(), rightVectorVnl ); downVectorVnl.normalize(); InitializeStandardPlane( rightVectorVnl, downVectorVnl ); SetOrigin(origin); } void PlaneGeometry::SetMatrixByVectors( const VnlVector &rightVector, const VnlVector &downVector, ScalarType thickness ) { VnlVector normal = vnl_cross_3d(rightVector, downVector); normal.normalize(); normal *= thickness; AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightVector); matrix.GetVnlMatrix().set_column(1, downVector); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); - transform->SetOffset(m_IndexToWorldTransform->GetOffset()); + transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset()); SetIndexToWorldTransform(transform); } Vector3D PlaneGeometry::GetNormal() const { Vector3D frontToBack; - frontToBack.SetVnlVector( m_IndexToWorldTransform + frontToBack.SetVnlVector( this->GetIndexToWorldTransform() ->GetMatrix().GetVnlMatrix().get_column(2) ); return frontToBack; } VnlVector PlaneGeometry::GetNormalVnl() const { - return m_IndexToWorldTransform + return this->GetIndexToWorldTransform() ->GetMatrix().GetVnlMatrix().get_column(2); } ScalarType PlaneGeometry::DistanceFromPlane( const Point3D &pt3d_mm ) const { return fabs(SignedDistance( pt3d_mm )); } ScalarType PlaneGeometry::SignedDistance( const Point3D &pt3d_mm ) const { return SignedDistanceFromPlane(pt3d_mm); } bool PlaneGeometry::IsAbove( const Point3D &pt3d_mm ) const { return SignedDistanceFromPlane(pt3d_mm) > 0; } bool PlaneGeometry::IntersectionLine( const PlaneGeometry* plane, Line3D& crossline ) const { Vector3D normal = this->GetNormal(); normal.Normalize(); Vector3D planeNormal = plane->GetNormal(); planeNormal.Normalize(); Vector3D direction = itk::CrossProduct( normal, planeNormal ); if ( direction.GetSquaredNorm() < eps ) return false; crossline.SetDirection( direction ); double N1dN2 = normal * planeNormal; double determinant = 1.0 - N1dN2 * N1dN2; Vector3D origin = this->GetOrigin().GetVectorFromOrigin(); Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin(); double d1 = normal * origin; double d2 = planeNormal * planeOrigin; double c1 = ( d1 - d2 * N1dN2 ) / determinant; double c2 = ( d2 - d1 * N1dN2 ) / determinant; Vector3D p = normal * c1 + planeNormal * c2; crossline.GetPoint().GetVnlVector() = p.GetVnlVector(); return true; } unsigned int PlaneGeometry::IntersectWithPlane2D( const PlaneGeometry* plane, Point2D& lineFrom, Point2D &lineTo ) const { Line3D crossline; if ( this->IntersectionLine( plane, crossline ) == false ) return 0; Point2D point2; Vector2D direction2; this->Map( crossline.GetPoint(), point2 ); this->Map( crossline.GetPoint(), crossline.GetDirection(), direction2 ); return Line3D::RectangleLineIntersection( 0, 0, GetExtentInMM(0), GetExtentInMM(1), point2, direction2, lineFrom, lineTo ); } double PlaneGeometry::Angle( const PlaneGeometry *plane ) const { return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2)); } double PlaneGeometry::Angle( const Line3D &line ) const { return vnl_math::pi_over_2 - angle( line.GetDirection().GetVnlVector(), GetMatrixColumn(2) ); } bool PlaneGeometry::IntersectionPoint( const Line3D &line, Point3D &intersectionPoint ) const { Vector3D planeNormal = this->GetNormal(); planeNormal.Normalize(); Vector3D lineDirection = line.GetDirection(); lineDirection.Normalize(); double t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; intersectionPoint = line.GetPoint() + lineDirection * t; return true; } bool PlaneGeometry::IntersectionPointParam( const Line3D &line, double &t ) const { Vector3D planeNormal = this->GetNormal(); Vector3D lineDirection = line.GetDirection(); t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; return true; } bool PlaneGeometry::IsParallel( const PlaneGeometry *plane ) const { return ( (Angle(plane) < 10.0 * mitk::sqrteps ) || ( Angle(plane) > ( vnl_math::pi - 10.0 * sqrteps ) ) ) ; } bool PlaneGeometry::IsOnPlane( const Point3D &point ) const { return Distance(point) < eps; } bool PlaneGeometry::IsOnPlane( const Line3D &line ) const { return ( (Distance( line.GetPoint() ) < eps) && (Distance( line.GetPoint2() ) < eps) ); } bool PlaneGeometry::IsOnPlane( const PlaneGeometry *plane ) const { return ( IsParallel( plane ) && (Distance( plane->GetOrigin() ) < eps) ); } Point3D PlaneGeometry::ProjectPointOntoPlane( const Point3D& pt ) const { ScalarType len = this->GetNormalVnl().two_norm(); return pt - this->GetNormal() * this->SignedDistanceFromPlane( pt ) / len; } itk::LightObject::Pointer PlaneGeometry::InternalClone() const { Self::Pointer newGeometry = new PlaneGeometry(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void PlaneGeometry::ExecuteOperation( Operation *operation ) { vtkTransform *transform = vtkTransform::New(); - transform->SetMatrix( m_VtkMatrix ); + transform->SetMatrix( this->GetVtkMatrix()); switch ( operation->GetOperationType() ) { case OpORIENT: { mitk::PlaneOperation *planeOp = dynamic_cast< mitk::PlaneOperation * >( operation ); if ( planeOp == NULL ) { return; } Point3D center = planeOp->GetPoint(); Vector3D orientationVector = planeOp->GetNormal(); Vector3D defaultVector; FillVector3D( defaultVector, 0.0, 0.0, 1.0 ); Vector3D rotationAxis = itk::CrossProduct( orientationVector, defaultVector ); //double rotationAngle = acos( orientationVector[2] / orientationVector.GetNorm() ); double rotationAngle = atan2( (double) rotationAxis.GetNorm(), (double) (orientationVector * defaultVector) ); rotationAngle *= 180.0 / vnl_math::pi; transform->PostMultiply(); transform->Identity(); transform->Translate( center[0], center[1], center[2] ); transform->RotateWXYZ( rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2] ); transform->Translate( -center[0], -center[1], -center[2] ); break; } case OpRESTOREPLANEPOSITION: { RestorePlanePositionOperation *op = dynamic_cast< mitk::RestorePlanePositionOperation* >(operation); if(op == NULL) { return; } AffineTransform3D::Pointer transform2 = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(0)); matrix.GetVnlMatrix().set_column(1, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(1)); matrix.GetVnlMatrix().set_column(2, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(2)); transform2->SetMatrix(matrix); Vector3D offset = op->GetTransform()->GetOffset(); transform2->SetOffset(offset); this->SetIndexToWorldTransform(transform2); ScalarType bounds[6] = {0, op->GetWidth(), 0, op->GetHeight(), 0 ,1 }; this->SetBounds(bounds); TransferItkToVtkTransform(); this->Modified(); transform->Delete(); return; } default: Superclass::ExecuteOperation( operation ); transform->Delete(); return; } - m_VtkMatrix->DeepCopy(transform->GetMatrix()); + this->GetVtkMatrix()->DeepCopy(transform->GetMatrix()); this->TransferVtkToItkTransform(); this->Modified(); transform->Delete(); } void PlaneGeometry::PrintSelf( std::ostream& os, itk::Indent indent ) const { Superclass::PrintSelf(os,indent); os << indent << " Normal: " << GetNormal() << std::endl; } } // namespace diff --git a/Core/Code/DataManagement/mitkSlicedGeometry3D.cpp b/Core/Code/DataManagement/mitkSlicedGeometry3D.cpp index 30d20fb3ef..58623d3f1a 100644 --- a/Core/Code/DataManagement/mitkSlicedGeometry3D.cpp +++ b/Core/Code/DataManagement/mitkSlicedGeometry3D.cpp @@ -1,1033 +1,1031 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkSlicedGeometry3D.h" #include "mitkPlaneGeometry.h" #include "mitkRotationOperation.h" #include "mitkPlaneOperation.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkApplyTransformMatrixOperation.h" #include "mitkInteractionConst.h" #include "mitkSliceNavigationController.h" const mitk::ScalarType PI = 3.14159265359; mitk::SlicedGeometry3D::SlicedGeometry3D() : m_EvenlySpaced( true ), m_Slices( 0 ), m_ReferenceGeometry( NULL ), m_SliceNavigationController( NULL ) { m_DirectionVector.Fill(0); this->InitializeSlicedGeometry( m_Slices ); } mitk::SlicedGeometry3D::SlicedGeometry3D(const SlicedGeometry3D& other) : Superclass(other), m_EvenlySpaced( other.m_EvenlySpaced ), m_Slices( other.m_Slices ), m_ReferenceGeometry( other.m_ReferenceGeometry ), m_SliceNavigationController( other.m_SliceNavigationController ) { m_DirectionVector.Fill(0); SetSpacing( other.GetSpacing() ); SetDirectionVector( other.GetDirectionVector() ); if ( m_EvenlySpaced ) { Geometry2D::Pointer geometry = other.m_Geometry2Ds[0]->Clone(); Geometry2D* geometry2D = dynamic_cast(geometry.GetPointer()); assert(geometry2D!=NULL); SetGeometry2D(geometry2D, 0); } else { unsigned int s; for ( s = 0; s < other.m_Slices; ++s ) { if ( other.m_Geometry2Ds[s].IsNull() ) { assert(other.m_EvenlySpaced); m_Geometry2Ds[s] = NULL; } else { Geometry2D* geometry2D = other.m_Geometry2Ds[s]->Clone(); assert(geometry2D!=NULL); SetGeometry2D(geometry2D, s); } } } } mitk::SlicedGeometry3D::~SlicedGeometry3D() { } mitk::Geometry2D * mitk::SlicedGeometry3D::GetGeometry2D( int s ) const { mitk::Geometry2D::Pointer geometry2D = NULL; if ( this->IsValidSlice(s) ) { geometry2D = m_Geometry2Ds[s]; // If (a) m_EvenlySpaced==true, (b) we don't have a Geometry2D stored // for the requested slice, and (c) the first slice (s=0) // is a PlaneGeometry instance, then we calculate the geometry of the // requested as the plane of the first slice shifted by m_Spacing[2]*s // in the direction of m_DirectionVector. if ( (m_EvenlySpaced) && (geometry2D.IsNull()) ) { PlaneGeometry *firstSlice = dynamic_cast< PlaneGeometry * > ( m_Geometry2Ds[0].GetPointer() ); if ( firstSlice != NULL ) { if ( (m_DirectionVector[0] == 0.0) && (m_DirectionVector[1] == 0.0) && (m_DirectionVector[2] == 0.0) ) { m_DirectionVector = firstSlice->GetNormal(); m_DirectionVector.Normalize(); } Vector3D direction; - direction = m_DirectionVector * m_Spacing[2]; + direction = m_DirectionVector * this->GetSpacing()[2]; mitk::PlaneGeometry::Pointer requestedslice; requestedslice = static_cast< mitk::PlaneGeometry * >( firstSlice->Clone().GetPointer() ); requestedslice->SetOrigin( requestedslice->GetOrigin() + direction * s ); geometry2D = requestedslice; m_Geometry2Ds[s] = geometry2D; } } return geometry2D; } else { return NULL; } } const mitk::BoundingBox * mitk::SlicedGeometry3D::GetBoundingBox() const { - assert(m_BoundingBox.IsNotNull()); - return m_BoundingBox.GetPointer(); + assert(this->IsBoundingBoxNull()==false); + return Superclass::GetBoundingBox(); } bool mitk::SlicedGeometry3D::SetGeometry2D( mitk::Geometry2D *geometry2D, int s ) { if ( this->IsValidSlice(s) ) { m_Geometry2Ds[s] = geometry2D; m_Geometry2Ds[s]->SetReferenceGeometry( m_ReferenceGeometry ); return true; } return false; } void mitk::SlicedGeometry3D::InitializeSlicedGeometry( unsigned int slices ) { Superclass::Initialize(); m_Slices = slices; Geometry2D::Pointer gnull = NULL; m_Geometry2Ds.assign( m_Slices, gnull ); Vector3D spacing; spacing.Fill( 1.0 ); this->SetSpacing( spacing ); m_DirectionVector.Fill( 0 ); } void mitk::SlicedGeometry3D::InitializeEvenlySpaced( mitk::Geometry2D* geometry2D, unsigned int slices, bool flipped ) { assert( geometry2D != NULL ); this->InitializeEvenlySpaced( geometry2D, geometry2D->GetExtentInMM(2)/geometry2D->GetExtent(2), slices, flipped ); } void mitk::SlicedGeometry3D::InitializeEvenlySpaced( mitk::Geometry2D* geometry2D, mitk::ScalarType zSpacing, unsigned int slices, bool flipped ) { assert( geometry2D != NULL ); assert( geometry2D->GetExtent(0) > 0 ); assert( geometry2D->GetExtent(1) > 0 ); geometry2D->Register(); Superclass::Initialize(); m_Slices = slices; BoundingBox::BoundsArrayType bounds = geometry2D->GetBounds(); bounds[4] = 0; bounds[5] = slices; // clear and reserve Geometry2D::Pointer gnull = NULL; m_Geometry2Ds.assign( m_Slices, gnull ); Vector3D directionVector = geometry2D->GetAxisVector(2); directionVector.Normalize(); directionVector *= zSpacing; if ( flipped == false ) { // Normally we should use the following four lines to create a copy of // the transform contrained in geometry2D, because it may not be changed // by us. But we know that SetSpacing creates a new transform without // changing the old (coming from geometry2D), so we can use the fifth // line instead. We check this at (**). // // AffineTransform3D::Pointer transform = AffineTransform3D::New(); // transform->SetMatrix(geometry2D->GetIndexToWorldTransform()->GetMatrix()); // transform->SetOffset(geometry2D->GetIndexToWorldTransform()->GetOffset()); // SetIndexToWorldTransform(transform); - m_IndexToWorldTransform = const_cast< AffineTransform3D * >( - geometry2D->GetIndexToWorldTransform() ); + this->SetIndexToWorldTransform( const_cast< AffineTransform3D * >( + geometry2D->GetIndexToWorldTransform() )); } else { directionVector *= -1.0; - m_IndexToWorldTransform = AffineTransform3D::New(); - m_IndexToWorldTransform->SetMatrix( + this->SetIndexToWorldTransform( AffineTransform3D::New()); + this->GetIndexToWorldTransform()->SetMatrix( geometry2D->GetIndexToWorldTransform()->GetMatrix() ); AffineTransform3D::OutputVectorType scaleVector; FillVector3D(scaleVector, 1.0, 1.0, -1.0); - m_IndexToWorldTransform->Scale(scaleVector, true); - m_IndexToWorldTransform->SetOffset( + this->GetIndexToWorldTransform()->Scale(scaleVector, true); + this->GetIndexToWorldTransform()->SetOffset( geometry2D->GetIndexToWorldTransform()->GetOffset() ); } mitk::Vector3D spacing; FillVector3D( spacing, geometry2D->GetExtentInMM(0) / bounds[1], geometry2D->GetExtentInMM(1) / bounds[3], zSpacing ); - // Ensure that spacing differs from m_Spacing to make SetSpacing change the - // matrix. - m_Spacing[2] = zSpacing - 1; - this->SetDirectionVector( directionVector ); this->SetBounds( bounds ); this->SetGeometry2D( geometry2D, 0 ); - this->SetSpacing( spacing ); + this->SetSpacing( spacing ,true); this->SetEvenlySpaced(); this->SetTimeBounds( geometry2D->GetTimeBounds() ); - assert(m_IndexToWorldTransform.GetPointer() + assert(this->GetIndexToWorldTransform() != geometry2D->GetIndexToWorldTransform()); // (**) see above. this->SetFrameOfReferenceID( geometry2D->GetFrameOfReferenceID() ); this->SetImageGeometry( geometry2D->GetImageGeometry() ); geometry2D->UnRegister(); } void mitk::SlicedGeometry3D::InitializePlanes( const mitk::Geometry3D *geometry3D, mitk::PlaneGeometry::PlaneOrientation planeorientation, bool top, bool frontside, bool rotated ) { m_ReferenceGeometry = const_cast< Geometry3D * >( geometry3D ); PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); planeGeometry->InitializeStandardPlane( geometry3D, top, planeorientation, frontside, rotated ); ScalarType viewSpacing = 1; unsigned int slices = 1; switch ( planeorientation ) { case PlaneGeometry::Axial: viewSpacing = geometry3D->GetSpacing()[2]; slices = (unsigned int) geometry3D->GetExtent( 2 ); break; case PlaneGeometry::Frontal: viewSpacing = geometry3D->GetSpacing()[1]; slices = (unsigned int) geometry3D->GetExtent( 1 ); break; case PlaneGeometry::Sagittal: viewSpacing = geometry3D->GetSpacing()[0]; slices = (unsigned int) geometry3D->GetExtent( 0 ); break; default: itkExceptionMacro("unknown PlaneOrientation"); } mitk::Vector3D normal = this->AdjustNormal( planeGeometry->GetNormal() ); ScalarType directedExtent = std::abs( m_ReferenceGeometry->GetExtentInMM( 0 ) * normal[0] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 1 ) * normal[1] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 2 ) * normal[2] ); if ( directedExtent >= viewSpacing ) { slices = static_cast< int >(directedExtent / viewSpacing + 0.5); } else { slices = 1; } bool flipped = (top == false); if ( frontside == false ) { flipped = !flipped; } if ( planeorientation == PlaneGeometry::Frontal ) { flipped = !flipped; } this->InitializeEvenlySpaced( planeGeometry, viewSpacing, slices, flipped ); } void mitk::SlicedGeometry3D ::ReinitializePlanes( const Point3D ¢er, const Point3D &referencePoint ) { // Need a reference frame to align the rotated planes if ( !m_ReferenceGeometry ) { return; } // Get first plane of plane stack PlaneGeometry *firstPlane = dynamic_cast< PlaneGeometry * >( m_Geometry2Ds[0].GetPointer() ); // If plane stack is empty, exit if ( firstPlane == NULL ) { return; } // Calculate the "directed" spacing when taking the plane (defined by its axes // vectors and normal) as the reference coordinate frame. // // This is done by calculating the radius of the ellipsoid defined by the // original volume spacing axes, in the direction of the respective axis of the // reference frame. mitk::Vector3D axis0 = firstPlane->GetAxisVector(0); mitk::Vector3D axis1 = firstPlane->GetAxisVector(1); mitk::Vector3D normal = firstPlane->GetNormal(); normal.Normalize(); Vector3D spacing; spacing[0] = this->CalculateSpacing( axis0 ); spacing[1] = this->CalculateSpacing( axis1 ); spacing[2] = this->CalculateSpacing( normal ); Superclass::SetSpacing( spacing ); // Now we need to calculate the number of slices in the plane's normal // direction, so that the entire volume is covered. This is done by first // calculating the dot product between the volume diagonal (the maximum // distance inside the volume) and the normal, and dividing this value by // the directed spacing calculated above. ScalarType directedExtent = std::abs( m_ReferenceGeometry->GetExtentInMM( 0 ) * normal[0] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 1 ) * normal[1] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 2 ) * normal[2] ); if ( directedExtent >= spacing[2] ) { m_Slices = static_cast< unsigned int >(directedExtent / spacing[2] + 0.5); } else { m_Slices = 1; } // The origin of our "first plane" needs to be adapted to this new extent. // To achieve this, we first calculate the current distance to the volume's // center, and then shift the origin in the direction of the normal by the // difference between this distance and half of the new extent. double centerOfRotationDistance = firstPlane->SignedDistanceFromPlane( center ); if ( centerOfRotationDistance > 0 ) { firstPlane->SetOrigin( firstPlane->GetOrigin() + normal * (centerOfRotationDistance - directedExtent / 2.0) ); m_DirectionVector = normal; } else { firstPlane->SetOrigin( firstPlane->GetOrigin() + normal * (directedExtent / 2.0 + centerOfRotationDistance) ); m_DirectionVector = -normal; } // Now we adjust this distance according with respect to the given reference // point: we need to make sure that the point is touched by one slice of the // new slice stack. double referencePointDistance = firstPlane->SignedDistanceFromPlane( referencePoint ); int referencePointSlice = static_cast< int >( referencePointDistance / spacing[2]); double alignmentValue = referencePointDistance / spacing[2] - referencePointSlice; firstPlane->SetOrigin( firstPlane->GetOrigin() + normal * alignmentValue * spacing[2] ); // Finally, we can clear the previous geometry stack and initialize it with // our re-initialized "first plane". m_Geometry2Ds.assign( m_Slices, Geometry2D::Pointer( NULL ) ); if ( m_Slices > 0 ) { m_Geometry2Ds[0] = firstPlane; } // Reinitialize SNC with new number of slices m_SliceNavigationController->GetSlice()->SetSteps( m_Slices ); this->Modified(); } double mitk::SlicedGeometry3D::CalculateSpacing( const mitk::Vector3D &d ) const { // Need the spacing of the underlying dataset / geometry if ( !m_ReferenceGeometry ) { return 1.0; } const mitk::Vector3D &spacing = m_ReferenceGeometry->GetSpacing(); return SlicedGeometry3D::CalculateSpacing( spacing, d ); } double mitk::SlicedGeometry3D::CalculateSpacing( const mitk::Vector3D spacing, const mitk::Vector3D &d ) { // The following can be derived from the ellipsoid equation // // 1 = x^2/a^2 + y^2/b^2 + z^2/c^2 // // where (a,b,c) = spacing of original volume (ellipsoid radii) // and (x,y,z) = scaled coordinates of vector d (according to ellipsoid) // double scaling = d[0]*d[0] / (spacing[0] * spacing[0]) + d[1]*d[1] / (spacing[1] * spacing[1]) + d[2]*d[2] / (spacing[2] * spacing[2]); scaling = sqrt( scaling ); return ( sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ) / scaling ); } mitk::Vector3D mitk::SlicedGeometry3D::AdjustNormal( const mitk::Vector3D &normal ) const { TransformType::Pointer inverse = TransformType::New(); m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse( inverse ); Vector3D transformedNormal = inverse->TransformVector( normal ); transformedNormal.Normalize(); return transformedNormal; } void mitk::SlicedGeometry3D::SetImageGeometry( const bool isAnImageGeometry ) { Superclass::SetImageGeometry( isAnImageGeometry ); mitk::Geometry3D* geometry; unsigned int s; for ( s = 0; s < m_Slices; ++s ) { geometry = m_Geometry2Ds[s]; if ( geometry!=NULL ) { geometry->SetImageGeometry( isAnImageGeometry ); } } } void mitk::SlicedGeometry3D::ChangeImageGeometryConsideringOriginOffset( const bool isAnImageGeometry ) { mitk::Geometry3D* geometry; unsigned int s; for ( s = 0; s < m_Slices; ++s ) { geometry = m_Geometry2Ds[s]; if ( geometry!=NULL ) { geometry->ChangeImageGeometryConsideringOriginOffset( isAnImageGeometry ); } } Superclass::ChangeImageGeometryConsideringOriginOffset( isAnImageGeometry ); } bool mitk::SlicedGeometry3D::IsValidSlice( int s ) const { return ((s >= 0) && (s < (int)m_Slices)); } void mitk::SlicedGeometry3D::SetReferenceGeometry( Geometry3D *referenceGeometry ) { m_ReferenceGeometry = referenceGeometry; std::vector::iterator it; for ( it = m_Geometry2Ds.begin(); it != m_Geometry2Ds.end(); ++it ) { (*it)->SetReferenceGeometry( referenceGeometry ); } } void mitk::SlicedGeometry3D::InternPreSetSpacing( const mitk::Vector3D &aSpacing ) { bool hasEvenlySpacedPlaneGeometry = false; mitk::Point3D origin; mitk::Vector3D rightDV, bottomDV; BoundingBox::BoundsArrayType bounds; assert(aSpacing[0]>0 && aSpacing[1]>0 && aSpacing[2]>0); // In case of evenly-spaced data: re-initialize instances of Geometry2D, // since the spacing influences them if ((m_EvenlySpaced) && (m_Geometry2Ds.size() > 0)) { mitk::Geometry2D::ConstPointer firstGeometry = m_Geometry2Ds[0].GetPointer(); const PlaneGeometry *planeGeometry = dynamic_cast< const PlaneGeometry * >( firstGeometry.GetPointer() ); if (planeGeometry != NULL ) { this->WorldToIndex( planeGeometry->GetOrigin(), origin ); this->WorldToIndex( planeGeometry->GetAxisVector(0), rightDV ); this->WorldToIndex( planeGeometry->GetAxisVector(1), bottomDV ); bounds = planeGeometry->GetBounds(); hasEvenlySpacedPlaneGeometry = true; } } InternSetSpacing(aSpacing); mitk::Geometry2D::Pointer firstGeometry; // In case of evenly-spaced data: re-initialize instances of Geometry2D, // since the spacing influences them if ( hasEvenlySpacedPlaneGeometry ) { //create planeGeometry according to new spacing this->IndexToWorld( origin, origin ); this->IndexToWorld( rightDV, rightDV ); this->IndexToWorld( bottomDV, bottomDV ); mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); planeGeometry->SetImageGeometry( this->GetImageGeometry() ); planeGeometry->SetReferenceGeometry( m_ReferenceGeometry ); + + //Store spacing, as Initialize... needs a pointer + mitk::Vector3D lokalSpacing = this->GetSpacing(); planeGeometry->InitializeStandardPlane( - rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &m_Spacing ); + rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &lokalSpacing ); planeGeometry->SetOrigin(origin); planeGeometry->SetBounds(bounds); firstGeometry = planeGeometry; } else if ( (m_EvenlySpaced) && (m_Geometry2Ds.size() > 0) ) { firstGeometry = m_Geometry2Ds[0].GetPointer(); } //clear and reserve Geometry2D::Pointer gnull=NULL; m_Geometry2Ds.assign(m_Slices, gnull); if ( m_Slices > 0 ) { m_Geometry2Ds[0] = firstGeometry; } this->Modified(); } void mitk::SlicedGeometry3D ::SetSliceNavigationController( SliceNavigationController *snc ) { m_SliceNavigationController = snc; } mitk::SliceNavigationController * mitk::SlicedGeometry3D::GetSliceNavigationController() { return m_SliceNavigationController; } void mitk::SlicedGeometry3D::SetEvenlySpaced(bool on) { if(m_EvenlySpaced!=on) { m_EvenlySpaced=on; this->Modified(); } } void mitk::SlicedGeometry3D ::SetDirectionVector( const mitk::Vector3D& directionVector ) { Vector3D newDir = directionVector; newDir.Normalize(); if ( newDir != m_DirectionVector ) { m_DirectionVector = newDir; this->Modified(); } } void mitk::SlicedGeometry3D::InternPostSetTimeBounds( const mitk::TimeBounds& timebounds ) { unsigned int s; for ( s = 0; s < m_Slices; ++s ) { if(m_Geometry2Ds[s].IsNotNull()) { m_Geometry2Ds[s]->SetTimeBounds( timebounds ); } } - m_TimeBounds = timebounds; } itk::LightObject::Pointer mitk::SlicedGeometry3D::InternalClone() const { Self::Pointer newGeometry = new SlicedGeometry3D(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void mitk::SlicedGeometry3D::PrintSelf( std::ostream& os, itk::Indent indent ) const { Superclass::PrintSelf(os,indent); os << indent << " EvenlySpaced: " << m_EvenlySpaced << std::endl; if ( m_EvenlySpaced ) { os << indent << " DirectionVector: " << m_DirectionVector << std::endl; } os << indent << " Slices: " << m_Slices << std::endl; os << std::endl; os << indent << " GetGeometry2D(0): "; if ( this->GetGeometry2D(0) == NULL ) { os << "NULL" << std::endl; } else { this->GetGeometry2D(0)->Print(os, indent); } } void mitk::SlicedGeometry3D::ExecuteOperation(Operation* operation) { switch ( operation->GetOperationType() ) { case OpNOTHING: break; case OpROTATE: if ( m_EvenlySpaced ) { // Need a reference frame to align the rotation if ( m_ReferenceGeometry ) { // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // Save first slice Geometry2D::Pointer geometry2D = m_Geometry2Ds[0]; RotationOperation *rotOp = dynamic_cast< RotationOperation * >( operation ); // Generate a RotationOperation using the dataset center instead of // the supplied rotation center. This is necessary so that the rotated // zero-plane does not shift away. The supplied center is instead used // to adjust the slice stack afterwards. Point3D center = m_ReferenceGeometry->GetCenter(); RotationOperation centeredRotation( rotOp->GetOperationType(), center, rotOp->GetVectorOfRotation(), rotOp->GetAngleOfRotation() ); // Rotate first slice geometry2D->ExecuteOperation( ¢eredRotation ); // Clear the slice stack and adjust it according to the center of // the dataset and the supplied rotation center (see documentation of // ReinitializePlanes) this->ReinitializePlanes( center, rotOp->GetCenterOfRotation() ); geometry2D->SetSpacing(this->GetSpacing()); if ( m_SliceNavigationController ) { m_SliceNavigationController->SelectSliceByPoint( rotOp->GetCenterOfRotation() ); m_SliceNavigationController->AdjustSliceStepperRange(); } Geometry3D::ExecuteOperation( ¢eredRotation ); } else { // we also have to consider the case, that there is no reference geometry available. if ( m_Geometry2Ds.size() > 0 ) { // Reach through to all slices in my container for (std::vector::iterator iter = m_Geometry2Ds.begin(); iter != m_Geometry2Ds.end(); ++iter) { // Test for empty slices, which can happen if evenly spaced geometry if ((*iter).IsNotNull()) { (*iter)->ExecuteOperation(operation); } } // rotate overall geometry RotationOperation *rotOp = dynamic_cast< RotationOperation * >( operation ); Geometry3D::ExecuteOperation( rotOp); } } } else { // Reach through to all slices for (std::vector::iterator iter = m_Geometry2Ds.begin(); iter != m_Geometry2Ds.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpORIENT: if ( m_EvenlySpaced ) { // get operation data PlaneOperation *planeOp = dynamic_cast< PlaneOperation * >( operation ); // Get first slice Geometry2D::Pointer geometry2D = m_Geometry2Ds[0]; PlaneGeometry *planeGeometry = dynamic_cast< PlaneGeometry * >( geometry2D.GetPointer() ); // Need a PlaneGeometry, a PlaneOperation and a reference frame to // carry out the re-orientation. If not all avaialble, stop here if ( !m_ReferenceGeometry || !planeGeometry || !planeOp ) { break; } // General Behavior: // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // // 1st Step: Reorient Normal Vector of first plane // Point3D center = planeOp->GetPoint(); //m_ReferenceGeometry->GetCenter(); mitk::Vector3D currentNormal = planeGeometry->GetNormal(); mitk::Vector3D newNormal; if (planeOp->AreAxisDefined()) { // If planeOp was defined by one centerpoint and two axis vectors newNormal = CrossProduct(planeOp->GetAxisVec0(), planeOp->GetAxisVec1()); } else { // If planeOp was defined by one centerpoint and one normal vector newNormal = planeOp->GetNormal(); } // Get Rotation axis und angle currentNormal.Normalize(); newNormal.Normalize(); ScalarType rotationAngle = angle(currentNormal.GetVnlVector(),newNormal.GetVnlVector()); rotationAngle *= 180.0 / vnl_math::pi; // from rad to deg Vector3D rotationAxis = itk::CrossProduct( currentNormal, newNormal ); if (std::abs(rotationAngle-180) < mitk::eps ) { // current Normal and desired normal are not linear independent!!(e.g 1,0,0 and -1,0,0). // Rotation Axis should be ANY vector that is 90� to current Normal mitk::Vector3D helpNormal; helpNormal = currentNormal; helpNormal[0] += 1; helpNormal[1] -= 1; helpNormal[2] += 1; helpNormal.Normalize(); rotationAxis = itk::CrossProduct( helpNormal, currentNormal ); } RotationOperation centeredRotation( mitk::OpROTATE, center, rotationAxis, rotationAngle ); // Rotate first slice geometry2D->ExecuteOperation( ¢eredRotation ); // Reinitialize planes and select slice, if my rotations are all done. if (!planeOp->AreAxisDefined()) { // Clear the slice stack and adjust it according to the center of // rotation and plane position (see documentation of ReinitializePlanes) this->ReinitializePlanes( center, planeOp->GetPoint() ); if ( m_SliceNavigationController ) { m_SliceNavigationController->SelectSliceByPoint( planeOp->GetPoint() ); m_SliceNavigationController->AdjustSliceStepperRange(); } } // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry) Geometry3D::ExecuteOperation( ¢eredRotation ); // // 2nd step. If axis vectors were defined, rotate the plane around its normal to fit these // if (planeOp->AreAxisDefined()) { mitk::Vector3D vecAxixNew = planeOp->GetAxisVec0(); vecAxixNew.Normalize(); mitk::Vector3D VecAxisCurr = geometry2D->GetAxisVector(0); VecAxisCurr.Normalize(); ScalarType rotationAngle = angle(VecAxisCurr.GetVnlVector(),vecAxixNew.GetVnlVector()); rotationAngle = rotationAngle * 180 / PI; // Rad to Deg // we rotate around the normal of the plane, but we do not know, if we need to rotate clockwise // or anti-clockwise. So we rotate around the crossproduct of old and new Axisvector. // Since both axis vectors lie in the plane, the crossproduct is the planes normal or the negative planes normal rotationAxis = itk::CrossProduct( VecAxisCurr, vecAxixNew ); if (std::abs(rotationAngle-180) < mitk::eps ) { // current axisVec and desired axisVec are not linear independent!!(e.g 1,0,0 and -1,0,0). // Rotation Axis can be just plane Normal. (have to rotate by 180�) rotationAxis = newNormal; } // Perfom Rotation mitk::RotationOperation op(mitk::OpROTATE, center, rotationAxis, rotationAngle); geometry2D->ExecuteOperation( &op ); // Apply changes on first slice to whole slice stack this->ReinitializePlanes( center, planeOp->GetPoint() ); if ( m_SliceNavigationController ) { m_SliceNavigationController->SelectSliceByPoint( planeOp->GetPoint() ); m_SliceNavigationController->AdjustSliceStepperRange(); } // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry) Geometry3D::ExecuteOperation( &op ); } } else { // Reach through to all slices for (std::vector::iterator iter = m_Geometry2Ds.begin(); iter != m_Geometry2Ds.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpRESTOREPLANEPOSITION: if ( m_EvenlySpaced ) { // Save first slice Geometry2D::Pointer geometry2D = m_Geometry2Ds[0]; PlaneGeometry* planeGeometry = dynamic_cast< PlaneGeometry * >( geometry2D.GetPointer() ); RestorePlanePositionOperation *restorePlaneOp = dynamic_cast< RestorePlanePositionOperation* >( operation ); // Need a PlaneGeometry, a PlaneOperation and a reference frame to // carry out the re-orientation if ( m_ReferenceGeometry && planeGeometry && restorePlaneOp ) { // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // Rotate first slice geometry2D->ExecuteOperation( restorePlaneOp ); m_DirectionVector = restorePlaneOp->GetDirectionVector(); double centerOfRotationDistance = planeGeometry->SignedDistanceFromPlane( m_ReferenceGeometry->GetCenter() ); if ( centerOfRotationDistance > 0 ) { m_DirectionVector = m_DirectionVector; } else { m_DirectionVector = -m_DirectionVector; } Vector3D spacing = restorePlaneOp->GetSpacing(); Superclass::SetSpacing( spacing ); // /*Now we need to calculate the number of slices in the plane's normal // direction, so that the entire volume is covered. This is done by first // calculating the dot product between the volume diagonal (the maximum // distance inside the volume) and the normal, and dividing this value by // the directed spacing calculated above.*/ ScalarType directedExtent = std::abs( m_ReferenceGeometry->GetExtentInMM( 0 ) * m_DirectionVector[0] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 1 ) * m_DirectionVector[1] ) + std::abs( m_ReferenceGeometry->GetExtentInMM( 2 ) * m_DirectionVector[2] ); if ( directedExtent >= spacing[2] ) { m_Slices = static_cast< unsigned int >(directedExtent / spacing[2] + 0.5); } else { m_Slices = 1; } m_Geometry2Ds.assign( m_Slices, Geometry2D::Pointer( NULL ) ); if ( m_Slices > 0 ) { m_Geometry2Ds[0] = geometry2D; } m_SliceNavigationController->GetSlice()->SetSteps( m_Slices ); this->Modified(); //End Reinitialization if ( m_SliceNavigationController ) { m_SliceNavigationController->GetSlice()->SetPos( restorePlaneOp->GetPos() ); m_SliceNavigationController->AdjustSliceStepperRange(); } Geometry3D::ExecuteOperation(restorePlaneOp); } } else { // Reach through to all slices for (std::vector::iterator iter = m_Geometry2Ds.begin(); iter != m_Geometry2Ds.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpAPPLYTRANSFORMMATRIX: // Clear all generated geometries and then transform only the first slice. // The other slices will be re-generated on demand // Save first slice Geometry2D::Pointer geometry2D = m_Geometry2Ds[0]; ApplyTransformMatrixOperation *applyMatrixOp = dynamic_cast< ApplyTransformMatrixOperation* >( operation ); // Apply transformation to first plane geometry2D->ExecuteOperation( applyMatrixOp ); // Generate a ApplyTransformMatrixOperation using the dataset center instead of // the supplied rotation center. The supplied center is instead used to adjust the // slice stack afterwards (see OpROTATE). Point3D center = m_ReferenceGeometry->GetCenter(); // Clear the slice stack and adjust it according to the center of // the dataset and the supplied rotation center (see documentation of // ReinitializePlanes) this->ReinitializePlanes( center, applyMatrixOp->GetReferencePoint() ); Geometry3D::ExecuteOperation( applyMatrixOp ); break; } this->Modified(); } diff --git a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp index 7e67da9e9c..d1120307c2 100644 --- a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp +++ b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp @@ -1,100 +1,100 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkThinPlateSplineCurvedGeometry.h" #include #include mitk::ThinPlateSplineCurvedGeometry::ThinPlateSplineCurvedGeometry() { m_InterpolatingAbstractTransform = m_ThinPlateSplineTransform = vtkThinPlateSplineTransform::New(); m_VtkTargetLandmarks = vtkPoints::New(); m_VtkProjectedLandmarks = vtkPoints::New(); m_ThinPlateSplineTransform->SetInverseIterations(5000); } mitk::ThinPlateSplineCurvedGeometry::ThinPlateSplineCurvedGeometry(const ThinPlateSplineCurvedGeometry& other ) : Superclass(other) { this->SetSigma(other.GetSigma()); } mitk::ThinPlateSplineCurvedGeometry::~ThinPlateSplineCurvedGeometry() { // don't need to delete m_ThinPlateSplineTransform, because it is // the same as m_InterpolatingAbstractTransform, which will be deleted // by the superclass. if(m_VtkTargetLandmarks!=NULL) m_VtkTargetLandmarks->Delete(); if(m_VtkProjectedLandmarks!=NULL) m_VtkProjectedLandmarks->Delete(); } -bool mitk::ThinPlateSplineCurvedGeometry::InternPostIsValid() const +bool mitk::ThinPlateSplineCurvedGeometry::IsValid() const { return m_TargetLandmarks.IsNotNull() && (m_TargetLandmarks->Size() >= 3) && m_LandmarkProjector.IsNotNull(); } void mitk::ThinPlateSplineCurvedGeometry::SetSigma(double sigma) { m_ThinPlateSplineTransform->SetSigma(sigma); } double mitk::ThinPlateSplineCurvedGeometry::GetSigma() const { return m_ThinPlateSplineTransform->GetSigma(); } void mitk::ThinPlateSplineCurvedGeometry::ComputeGeometry() { Superclass::ComputeGeometry(); const mitk::PointSet::DataType::PointsContainer *finalTargetLandmarks, *projectedTargetLandmarks; finalTargetLandmarks = m_LandmarkProjector->GetFinalTargetLandmarks(); projectedTargetLandmarks = m_LandmarkProjector->GetProjectedLandmarks(); mitk::PointSet::DataType::PointsContainer::ConstIterator targetIt, projectedIt; targetIt = finalTargetLandmarks->Begin(); projectedIt = projectedTargetLandmarks->Begin(); //initialize Thin-Plate-Spline m_VtkTargetLandmarks->Reset(); m_VtkProjectedLandmarks->Reset(); vtkIdType id; int size=finalTargetLandmarks->Size(); for(id=0; id < size; ++id, ++targetIt, ++projectedIt) { const mitk::PointSet::PointType& target = targetIt->Value(); m_VtkTargetLandmarks->InsertPoint(id, target[0], target[1], target[2]); const mitk::PointSet::PointType& projected = projectedIt->Value(); m_VtkProjectedLandmarks->InsertPoint(id, projected[0], projected[1], projected[2]); } m_VtkTargetLandmarks->Modified(); m_VtkProjectedLandmarks->Modified(); m_ThinPlateSplineTransform->SetSourceLandmarks(m_VtkProjectedLandmarks); m_ThinPlateSplineTransform->SetTargetLandmarks(m_VtkTargetLandmarks); } itk::LightObject::Pointer mitk::ThinPlateSplineCurvedGeometry::InternalClone() const { mitk::Geometry3D::Pointer newGeometry = new Self(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } diff --git a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.h b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.h index 4c64df35ae..78c3c3e53e 100644 --- a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.h +++ b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.h @@ -1,64 +1,64 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKTHINPLATESPLINECURVEDGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #define MITKTHINPLATESPLINECURVEDGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #include "mitkLandmarkProjectorBasedCurvedGeometry.h" class vtkPoints; class vtkThinPlateSplineTransform; namespace mitk { -//##Documentation -//## @brief Thin-plate-spline-based landmark-based curved geometry -//## -//## @ingroup Geometry -class MITK_CORE_EXPORT ThinPlateSplineCurvedGeometry : public LandmarkProjectorBasedCurvedGeometry -{ -public: - mitkClassMacro(ThinPlateSplineCurvedGeometry, LandmarkProjectorBasedCurvedGeometry); + //##Documentation + //## @brief Thin-plate-spline-based landmark-based curved geometry + //## + //## @ingroup Geometry + class MITK_CORE_EXPORT ThinPlateSplineCurvedGeometry : public LandmarkProjectorBasedCurvedGeometry + { + public: + mitkClassMacro(ThinPlateSplineCurvedGeometry, LandmarkProjectorBasedCurvedGeometry); - itkNewMacro(Self); + itkNewMacro(Self); - virtual void ComputeGeometry(); + virtual void ComputeGeometry(); - virtual itk::LightObject::Pointer InternalClone() const; + virtual itk::LightObject::Pointer InternalClone() const; - vtkThinPlateSplineTransform* GetThinPlateSplineTransform() const - { - return m_ThinPlateSplineTransform; - } + vtkThinPlateSplineTransform* GetThinPlateSplineTransform() const + { + return m_ThinPlateSplineTransform; + } - virtual void SetSigma(double sigma); - virtual double GetSigma() const; + virtual void SetSigma(double sigma); + virtual double GetSigma() const; -protected: - ThinPlateSplineCurvedGeometry(); - ThinPlateSplineCurvedGeometry(const ThinPlateSplineCurvedGeometry& other ); + virtual bool IsValid() const; - virtual ~ThinPlateSplineCurvedGeometry(); + protected: + ThinPlateSplineCurvedGeometry(); + ThinPlateSplineCurvedGeometry(const ThinPlateSplineCurvedGeometry& other ); - vtkThinPlateSplineTransform* m_ThinPlateSplineTransform; + virtual ~ThinPlateSplineCurvedGeometry(); - vtkPoints* m_VtkTargetLandmarks; - vtkPoints* m_VtkProjectedLandmarks; + vtkThinPlateSplineTransform* m_ThinPlateSplineTransform; - virtual bool InternPostIsValid() const; -}; + vtkPoints* m_VtkTargetLandmarks; + vtkPoints* m_VtkProjectedLandmarks; + }; } // namespace mitk #endif /* MITKTHINPLATESPLINECURVEDGEOMETRY_H_HEADER_INCLUDED_C1C68A2C */ diff --git a/Core/Code/Testing/mitkBaseGeometryTest.cpp b/Core/Code/Testing/mitkBaseGeometryTest.cpp index 322b2a5723..8dd907acb1 100644 --- a/Core/Code/Testing/mitkBaseGeometryTest.cpp +++ b/Core/Code/Testing/mitkBaseGeometryTest.cpp @@ -1,513 +1,1103 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include #include #include #include #include "mitkoperationactor.h" #include #include "mitkvector.h" #include #include #include "itkScalableAffineTransform.h" #include #include #include +#include "mitkRotationOperation.h" +#include "mitkInteractionConst.h" +#include +#include + class vtkMatrix4x4; class vtkMatrixToLinearTransform; class vtkLinearTransform; typedef itk::BoundingBox BoundingBox; typedef itk::BoundingBox BoundingBoxType; typedef BoundingBoxType::BoundsArrayType BoundsArrayType; typedef BoundingBoxType::Pointer BoundingBoxPointer; // Dummy instance of abstract base class class DummyTestClass : public mitk::BaseGeometry { public: DummyTestClass(){}; DummyTestClass(const DummyTestClass& other) : BaseGeometry(other){}; ~DummyTestClass(){}; mitkClassMacro(DummyTestClass, mitk::BaseGeometry); itkNewMacro(Self); mitkNewMacro1Param(Self,Self); - itkGetConstMacro(IndexToWorldTransformLastModified, unsigned long); - itkGetConstMacro(VtkMatrix, vtkMatrix4x4*); - - void DummyTestClass::ExecuteOperation(mitk::Operation* operation){}; itk::LightObject::Pointer DummyTestClass::InternalClone() const { Self::Pointer newGeometry = new Self(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } }; class mitkBaseGeometryTestSuite : public mitk::TestFixture { // List of Tests CPPUNIT_TEST_SUITE(mitkBaseGeometryTestSuite); //Constructor MITK_TEST(TestConstructors); //Set MITK_TEST(TestSetOrigin); MITK_TEST(TestSetBounds); MITK_TEST(TestSetFloatBounds); MITK_TEST(TestSetFloatBoundsDouble); MITK_TEST(TestSetFrameOfReferenceID); MITK_TEST(TestSetIndexToWorldTransform); MITK_TEST(TestSetSpacing); MITK_TEST(TestTransferItkToVtkTransform); MITK_TEST(TestSetIndexToWorldTransformByVtkMatrix); MITK_TEST(TestSetIdentity); //Equal MITK_TEST(Equal_CloneAndOriginal_ReturnsTrue); MITK_TEST(Equal_DifferentOrigin_ReturnsFalse); MITK_TEST(Equal_DifferentIndexToWorldTransform_ReturnsFalse); MITK_TEST(Equal_DifferentSpacing_ReturnsFalse); MITK_TEST(Equal_InputIsNull_ReturnsFalse); MITK_TEST(Equal_DifferentBoundingBox_ReturnsFalse); + //other Functions MITK_TEST(TestComposeTransform); MITK_TEST(TestComposeVtkMatrix); + MITK_TEST(TestTranslate); + MITK_TEST(TestIndexToWorld); + MITK_TEST(TestExecuteOperation); + MITK_TEST(TestCalculateBoundingBoxRelToTransform); + MITK_TEST(TestSetTimeBounds); + MITK_TEST(TestIs2DConvertable); + MITK_TEST(TestGetCornerPoint); + MITK_TEST(TestExtentInMM); + MITK_TEST(TestGetAxisVector); + MITK_TEST(TestGetCenter); CPPUNIT_TEST_SUITE_END(); // Used Variables private: mitk::Point3D aPoint; float aFloatSpacing[3]; mitk::Vector3D aSpacing; mitk::AffineTransform3D::Pointer aTransform; BoundingBoxPointer aBoundingBox; mitk::AffineTransform3D::MatrixType aMatrix; mitk::Point3D anotherPoint; mitk::Vector3D anotherSpacing; BoundingBoxPointer anotherBoundingBox; BoundingBoxPointer aThirdBoundingBox; mitk::AffineTransform3D::Pointer anotherTransform; mitk::AffineTransform3D::Pointer aThirdTransform; mitk::AffineTransform3D::MatrixType anotherMatrix; mitk::AffineTransform3D::MatrixType aThirdMatrix; DummyTestClass::Pointer aDummyGeometry; DummyTestClass::Pointer anotherDummyGeometry; public: // Set up for variables void setUp() { mitk::FillVector3D(aFloatSpacing, 1,1,1); mitk::FillVector3D(aSpacing, 1,1,1); mitk::FillVector3D(aPoint, 0,0,0); //Transform aTransform = mitk::AffineTransform3D::New(); aTransform->SetIdentity(); aMatrix.SetIdentity(); anotherTransform = mitk::AffineTransform3D::New(); anotherMatrix.SetIdentity(); anotherMatrix(1,1) = 2; anotherTransform->SetMatrix( anotherMatrix ); aThirdTransform = mitk::AffineTransform3D::New(); aThirdMatrix.SetIdentity(); aThirdMatrix(1,1) = 7; aThirdTransform->SetMatrix( aThirdMatrix ); //Bounding Box float bounds[6] = {0,1,0,1,0,1}; mitk::BoundingBox::BoundsArrayType b; const float *input = bounds; int j=0; for(mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); j < 6 ;++j) *it++ = (mitk::ScalarType)*input++; aBoundingBox = BoundingBoxType::New(); BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New(); BoundingBoxType::PointType p; BoundingBoxType::PointIdentifier pointid; for(pointid=0; pointid<2;++pointid) { unsigned int i; for(i=0; i<3; ++i) { p[i] = bounds[2*i+pointid]; } pointscontainer->InsertElement(pointid, p); } aBoundingBox->SetPoints(pointscontainer); aBoundingBox->ComputeBoundingBox(); anotherBoundingBox = BoundingBoxType::New(); p[0]=11; p[1]=12; p[2]=13; pointscontainer->InsertElement(1, p); anotherBoundingBox->SetPoints(pointscontainer); anotherBoundingBox->ComputeBoundingBox(); aThirdBoundingBox = BoundingBoxType::New(); p[0]=22; p[1]=23; p[2]=24; pointscontainer->InsertElement(1, p); aThirdBoundingBox->SetPoints(pointscontainer); aThirdBoundingBox->ComputeBoundingBox(); mitk::FillVector3D(anotherPoint, 2,3,4); - mitk::FillVector3D(anotherSpacing, 5,6,7); + mitk::FillVector3D(anotherSpacing, 5,6.5,7); aDummyGeometry = DummyTestClass::New(); aDummyGeometry->Initialize(); anotherDummyGeometry = aDummyGeometry->Clone(); } void tearDown() { aDummyGeometry = NULL; anotherDummyGeometry = NULL; } // Test functions void TestSetOrigin() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetOrigin(anotherPoint); CPPUNIT_ASSERT(anotherPoint==dummy->GetOrigin()); //undo changes, new and changed object need to be the same! dummy->SetOrigin(aPoint); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetFloatBounds(){ float bounds[6] = {0,11,0,12,0,13}; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFloatBounds(bounds); CPPUNIT_ASSERT(mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, true)); //Wrong bounds, test needs to fail bounds[1]=7; dummy->SetFloatBounds(bounds); CPPUNIT_ASSERT((mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, false))==false); //undo changes, new and changed object need to be the same! float originalBounds[6] = {0,1,0,1,0,1}; dummy->SetFloatBounds(originalBounds); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetBounds(){ DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetBounds(anotherBoundingBox->GetBounds()); CPPUNIT_ASSERT(mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, true)); //Test needs to fail now dummy->SetBounds(aThirdBoundingBox->GetBounds()); CPPUNIT_ASSERT(mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, false)==false); //undo changes, new and changed object need to be the same! dummy->SetBounds(aBoundingBox->GetBounds()); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetFloatBoundsDouble(){ double bounds[6] = {0,11,0,12,0,13}; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFloatBounds(bounds); CPPUNIT_ASSERT(mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, true)); //Test needs to fail now bounds[3]=7; dummy->SetFloatBounds(bounds); CPPUNIT_ASSERT(mitk::Equal( dummy->GetBoundingBox(), anotherBoundingBox, mitk::eps, false)==false); //undo changes, new and changed object need to be the same! double originalBounds[6] = {0,1,0,1,0,1}; dummy->SetFloatBounds(originalBounds); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetFrameOfReferenceID() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFrameOfReferenceID(5); CPPUNIT_ASSERT(dummy->GetFrameOfReferenceID()==5); //undo changes, new and changed object need to be the same! dummy->SetFrameOfReferenceID(0); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetIndexToWorldTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); CPPUNIT_ASSERT(mitk::Equal(anotherTransform,dummy->GetIndexToWorldTransform(),mitk::eps,true)); //Test needs to fail now dummy->SetIndexToWorldTransform(aThirdTransform); CPPUNIT_ASSERT(mitk::Equal(anotherTransform,dummy->GetIndexToWorldTransform(),mitk::eps,false)==false); //undo changes, new and changed object need to be the same! dummy->SetIndexToWorldTransform(aTransform); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetIndexToWorldTransformByVtkMatrix() { vtkMatrix4x4* vtkmatrix; vtkmatrix = vtkMatrix4x4::New(); vtkmatrix->Identity(); vtkmatrix->SetElement(1,1,2); DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); CPPUNIT_ASSERT(mitk::Equal(anotherTransform,dummy->GetIndexToWorldTransform(),mitk::eps,true)); //test needs to fail now vtkmatrix->SetElement(1,1,7); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); CPPUNIT_ASSERT(mitk::Equal(anotherTransform,dummy->GetIndexToWorldTransform(),mitk::eps,false)==false); //undo changes, new and changed object need to be the same! vtkmatrix->SetElement(1,1,1); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetIdentity() { DummyTestClass::Pointer dummy = DummyTestClass::New(); //Change IndextoWorldTransform and Origin dummy->SetIndexToWorldTransform(anotherTransform); dummy->SetOrigin(anotherPoint); //Set Identity should reset ITWT and Origin dummy->SetIdentity(); CPPUNIT_ASSERT(mitk::Equal(aTransform,dummy->GetIndexToWorldTransform(),mitk::eps,true)); CPPUNIT_ASSERT(aPoint==dummy->GetOrigin()); CPPUNIT_ASSERT(aSpacing==dummy->GetSpacing()); //new and changed object need to be the same! DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestSetSpacing() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetSpacing(anotherSpacing); CPPUNIT_ASSERT(anotherSpacing==dummy->GetSpacing()); //undo changes, new and changed object need to be the same! dummy->SetSpacing(aSpacing); DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); } void TestTransferItkToVtkTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); //calls TransferItkToVtkTransform mitk::AffineTransform3D::Pointer dummyTransform = dummy->GetIndexToWorldTransform(); CPPUNIT_ASSERT(mitk::MatrixEqualElementWise( anotherMatrix, dummyTransform->GetMatrix() )); } void TestConstructors() { //test standard constructor DummyTestClass::Pointer dummy1 = DummyTestClass::New(); bool test = dummy1->IsValid(); CPPUNIT_ASSERT(test == true); CPPUNIT_ASSERT(dummy1->GetFrameOfReferenceID() == 0); CPPUNIT_ASSERT(dummy1->GetIndexToWorldTransformLastModified() == 0); const float *dummy1FloatSpacing = dummy1->GetFloatSpacing(); CPPUNIT_ASSERT(dummy1FloatSpacing[0]==aFloatSpacing[0]); CPPUNIT_ASSERT(dummy1FloatSpacing[1]==aFloatSpacing[1]); CPPUNIT_ASSERT(dummy1FloatSpacing[2]==aFloatSpacing[2]); CPPUNIT_ASSERT(dummy1->GetSpacing() == aSpacing); CPPUNIT_ASSERT(dummy1->GetOrigin()==aPoint); CPPUNIT_ASSERT(mitk::Equal( dummy1->GetIndexToWorldTransform(), aTransform, mitk::eps, true)); CPPUNIT_ASSERT(mitk::Equal( dummy1->GetBoundingBox(), aBoundingBox, mitk::eps, true)); DummyTestClass::Pointer dummy2 = DummyTestClass::New(); dummy2->SetOrigin(anotherPoint); float bounds[6] = {0,11,0,12,0,13}; dummy2->SetFloatBounds(bounds); dummy2->SetIndexToWorldTransform(anotherTransform); dummy2->SetSpacing(anotherSpacing); DummyTestClass::Pointer dummy3; dummy3 = DummyTestClass::New(*dummy2); CPPUNIT_ASSERT(mitk::Equal(dummy3,dummy2,mitk::eps,true)); } //Equal Tests void Equal_CloneAndOriginal_ReturnsTrue() { CPPUNIT_ASSERT( mitk::Equal(aDummyGeometry, anotherDummyGeometry, mitk::eps,true)); } void Equal_DifferentOrigin_ReturnsFalse() { anotherDummyGeometry->SetOrigin(anotherPoint); CPPUNIT_ASSERT( mitk::Equal(aDummyGeometry, anotherDummyGeometry, mitk::eps,false)==false); } void Equal_DifferentIndexToWorldTransform_ReturnsFalse() { anotherDummyGeometry->SetIndexToWorldTransform(anotherTransform); CPPUNIT_ASSERT( mitk::Equal(aDummyGeometry, anotherDummyGeometry, mitk::eps,false)==false); } void Equal_DifferentSpacing_ReturnsFalse() { anotherDummyGeometry->SetSpacing(anotherSpacing); CPPUNIT_ASSERT( mitk::Equal(aDummyGeometry, anotherDummyGeometry, mitk::eps,false)==false); } void Equal_InputIsNull_ReturnsFalse() { DummyTestClass::Pointer geometryNull = NULL; CPPUNIT_ASSERT( mitk::Equal(geometryNull, anotherDummyGeometry, mitk::eps,false)==false); } void Equal_DifferentBoundingBox_ReturnsFalse() { //create different bounds to make the comparison false mitk::ScalarType bounds[ ] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; anotherDummyGeometry->SetBounds(bounds); CPPUNIT_ASSERT( mitk::Equal(aDummyGeometry, anotherDummyGeometry, mitk::eps,false)==false); } void TestComposeTransform(){ //Create Transformations to set and compare mitk::AffineTransform3D::Pointer transform1; transform1 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix1; matrix1.SetIdentity(); matrix1(1,1) = 2; transform1->SetMatrix( matrix1 ); //Spacing = 2 mitk::AffineTransform3D::Pointer transform2; transform2 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix2; matrix2.SetIdentity(); matrix2(1,1) = 2; transform2->SetMatrix( matrix2 ); //Spacing = 2 mitk::AffineTransform3D::Pointer transform3; transform3 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix3; matrix3.SetIdentity(); matrix3(1,1) = 4; transform3->SetMatrix( matrix3 ); //Spacing = 4 mitk::AffineTransform3D::Pointer transform4; transform4 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix4; matrix4.SetIdentity(); matrix4(1,1) = 0.25; transform4->SetMatrix( matrix4 ); //Spacing = 0.25 + //Vector to compare spacing + mitk::Vector3D expectedSpacing; + expectedSpacing.Fill(1.0); + expectedSpacing[1] = 4; + DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(transform1); //Spacing = 2 dummy->Compose(transform2); //Spacing = 4 - + CPPUNIT_ASSERT(mitk::Equal(dummy->GetSpacing(), expectedSpacing)); CPPUNIT_ASSERT(mitk::Equal(transform3,dummy->GetIndexToWorldTransform(),mitk::eps,true)); // 4=4 //undo changes, new and changed object need to be the same! dummy->Compose(transform4); //Spacing = 1 DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); // 1=1 } void TestComposeVtkMatrix(){ //Create Transformations to set and compare mitk::AffineTransform3D::Pointer transform1; transform1 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix1; matrix1.SetIdentity(); matrix1(1,1) = 2; transform1->SetMatrix( matrix1 ); //Spacing = 2 vtkMatrix4x4* vtkmatrix2; vtkmatrix2 = vtkMatrix4x4::New(); vtkmatrix2->Identity(); vtkmatrix2->SetElement(1,1,2); //Spacing = 2 mitk::AffineTransform3D::Pointer transform3; transform3 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix3; matrix3.SetIdentity(); matrix3(1,1) = 4; transform3->SetMatrix( matrix3 ); //Spacing = 4 vtkMatrix4x4* vtkmatrix4; vtkmatrix4 = vtkMatrix4x4::New(); vtkmatrix4->Identity(); vtkmatrix4->SetElement(1,1,0.25); //Spacing = 0.25 + //Vector to compare spacing + mitk::Vector3D expectedSpacing; + expectedSpacing.Fill(1.0); + expectedSpacing[1] = 4; + DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(transform1); //Spacing = 2 dummy->Compose(vtkmatrix2); //Spacing = 4 CPPUNIT_ASSERT(mitk::Equal(transform3,dummy->GetIndexToWorldTransform(),mitk::eps,true)); // 4=4 + CPPUNIT_ASSERT(mitk::Equal(dummy->GetSpacing(), expectedSpacing)); //undo changes, new and changed object need to be the same! dummy->Compose(vtkmatrix4); //Spacing = 1 DummyTestClass::Pointer newDummy = DummyTestClass::New(); CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); // 1=1 } + + void TestTranslate(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetOrigin(anotherPoint); + CPPUNIT_ASSERT(anotherPoint==dummy->GetOrigin()); + + //use some random values for translation + mitk::Vector3D translationVector; + translationVector.SetElement(0, 17.5f); + translationVector.SetElement(1, -32.3f); + translationVector.SetElement(2, 4.0f); + + //compute ground truth + mitk::Point3D tmpResult = anotherPoint + translationVector; + dummy->Translate(translationVector); + CPPUNIT_ASSERT( mitk::Equal( dummy->GetOrigin(), tmpResult )); + + //undo changes + translationVector*=-1; + dummy->Translate(translationVector); + CPPUNIT_ASSERT( mitk::Equal( dummy->GetOrigin(), anotherPoint )); + + //undo changes, new and changed object need to be the same! + translationVector.SetElement(0, -1 * anotherPoint[0]); + translationVector.SetElement(1, -1 * anotherPoint[1]); + translationVector.SetElement(2, -1 * anotherPoint[2]); + dummy->Translate(translationVector); + + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + // a part of the test requires axis-parallel coordinates + int testIndexAndWorldConsistency(DummyTestClass* dummyGeometry) + { + //Testing consistency of index and world coordinate systems + mitk::Point3D origin = dummyGeometry->GetOrigin(); + mitk::Point3D dummyPoint; + + //Testing index->world->index conversion consistency + dummyGeometry->WorldToIndex(origin, dummyPoint); + dummyGeometry->IndexToWorld(dummyPoint, dummyPoint); + CPPUNIT_ASSERT(dummyPoint == origin); + + //Testing WorldToIndex(origin, mitk::Point3D)==(0,0,0) + mitk::Point3D globalOrigin; + mitk::FillVector3D(globalOrigin, 0,0,0); + + mitk::Point3D originContinuousIndex; + dummyGeometry->WorldToIndex(origin, originContinuousIndex); + CPPUNIT_ASSERT(originContinuousIndex == globalOrigin); + + //Testing WorldToIndex(origin, itk::Index)==(0,0,0) + itk::Index<3> itkindex; + dummyGeometry->WorldToIndex(origin, itkindex); + itk::Index<3> globalOriginIndex; + mitk::vtk2itk(globalOrigin, globalOriginIndex); + CPPUNIT_ASSERT(itkindex == globalOriginIndex); + + //Testing WorldToIndex(origin-0.5*spacing, itk::Index)==(0,0,0) + mitk::Vector3D halfSpacingStep = dummyGeometry->GetSpacing()*0.5; + mitk::Matrix3D rotation; + mitk::Point3D originOffCenter = origin-halfSpacingStep; + dummyGeometry->WorldToIndex(originOffCenter, itkindex); + CPPUNIT_ASSERT(itkindex == globalOriginIndex); + + //Testing WorldToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0) + originOffCenter = origin+halfSpacingStep; + originOffCenter -= 0.0001; + dummyGeometry->WorldToIndex( originOffCenter, itkindex); + CPPUNIT_ASSERT(itkindex == globalOriginIndex); + + //Testing WorldToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); + originOffCenter = origin+halfSpacingStep; + itk::Index<3> global111; + mitk::FillVector3D(global111, 1,1,1); + dummyGeometry->WorldToIndex( originOffCenter, itkindex); + CPPUNIT_ASSERT(itkindex == global111); + + //Testing WorldToIndex(GetCenter())==BoundingBox.GetCenter + mitk::Point3D center = dummyGeometry->GetCenter(); + mitk::Point3D centerContIndex; + dummyGeometry->WorldToIndex(center, centerContIndex); + mitk::BoundingBox::ConstPointer boundingBox = dummyGeometry->GetBoundingBox(); + mitk::BoundingBox::PointType centerBounds = boundingBox->GetCenter(); + CPPUNIT_ASSERT(mitk::Equal(centerContIndex,centerBounds)); + + //Testing GetCenter()==IndexToWorld(BoundingBox.GetCenter) + center = dummyGeometry->GetCenter(); + mitk::Point3D centerBoundsInWorldCoords; + dummyGeometry->IndexToWorld(centerBounds, centerBoundsInWorldCoords); + CPPUNIT_ASSERT(mitk::Equal(center,centerBoundsInWorldCoords)); + + //Test using random point, + //Testing consistency of index and world coordinate systems + mitk::Point3D point; + mitk::FillVector3D(point,3.5,-2,4.6); + + //Testing index->world->index conversion consistency + dummyGeometry->WorldToIndex(point, dummyPoint); + dummyGeometry->IndexToWorld(dummyPoint, dummyPoint); + CPPUNIT_ASSERT(dummyPoint == point); + + return EXIT_SUCCESS; + } + + int testIndexAndWorldConsistencyForVectors(DummyTestClass* dummyGeometry) + { + //Testing consistency of index and world coordinate systems for vectors + mitk::Vector3D xAxisMM = dummyGeometry->GetAxisVector(0); + mitk::Vector3D xAxisContinuousIndex; + + mitk::Point3D p, pIndex, origin; + origin = dummyGeometry->GetOrigin(); + p[0] = xAxisMM[0]+origin[0]; + p[1] = xAxisMM[1]+origin[1]; + p[2] = xAxisMM[2]+origin[2]; + + dummyGeometry->WorldToIndex(p,pIndex); + + dummyGeometry->WorldToIndex(xAxisMM,xAxisContinuousIndex); + CPPUNIT_ASSERT(xAxisContinuousIndex[0] == pIndex[0]); + CPPUNIT_ASSERT(xAxisContinuousIndex[1] == pIndex[1]); + CPPUNIT_ASSERT(xAxisContinuousIndex[2] == pIndex[2]); + + dummyGeometry->IndexToWorld(xAxisContinuousIndex,xAxisContinuousIndex); + + dummyGeometry->IndexToWorld(pIndex,p); + + CPPUNIT_ASSERT(xAxisContinuousIndex == xAxisMM); + CPPUNIT_ASSERT(xAxisContinuousIndex[0] == p[0]-origin[0]); + CPPUNIT_ASSERT(xAxisContinuousIndex[1] == p[1]-origin[1]); + CPPUNIT_ASSERT(xAxisContinuousIndex[2] == p[2]-origin[2]); + + //Test consictency for random vector + mitk::Vector3D vector; + mitk::FillVector3D(vector, 2.5,-3.2,8.1); + mitk::Vector3D vectorContinuousIndex; + + p[0] = vector[0]+origin[0]; + p[1] = vector[1]+origin[1]; + p[2] = vector[2]+origin[2]; + + dummyGeometry->WorldToIndex(p,pIndex); + + dummyGeometry->WorldToIndex(vector,vectorContinuousIndex); + CPPUNIT_ASSERT(vectorContinuousIndex[0] == pIndex[0]); + CPPUNIT_ASSERT(vectorContinuousIndex[1] == pIndex[1]); + CPPUNIT_ASSERT(vectorContinuousIndex[2] == pIndex[2]); + + dummyGeometry->IndexToWorld(vectorContinuousIndex,vectorContinuousIndex); + + dummyGeometry->IndexToWorld(pIndex,p); + + CPPUNIT_ASSERT(vectorContinuousIndex == vector); + CPPUNIT_ASSERT(vectorContinuousIndex[0] == p[0]-origin[0]); + CPPUNIT_ASSERT(vectorContinuousIndex[1] == p[1]-origin[1]); + CPPUNIT_ASSERT(vectorContinuousIndex[2] == p[2]-origin[2]); + + return EXIT_SUCCESS; + } + + int testIndexAndWorldConsistencyForIndex(DummyTestClass* dummyGeometry) + { + //Testing consistency of index and world coordinate systems + + // creating testing data + itk::Index<4> itkIndex4, itkIndex4b; + itk::Index<3> itkIndex3, itkIndex3b; + itk::Index<2> itkIndex2, itkIndex2b; + mitk::Index3D mitkIndex, mitkIndexb; + + itkIndex4[0] = itkIndex4[1] = itkIndex4[2] = itkIndex4[3] = 4; + itkIndex3[0] = itkIndex3[1] = itkIndex3[2] = 6; + itkIndex2[0] = itkIndex2[1] = 2; + mitkIndex[0] = mitkIndex[1] = mitkIndex[2] = 13; + + // check for constistency + mitk::Point3D point; + dummyGeometry->IndexToWorld(itkIndex2,point); + dummyGeometry->WorldToIndex(point,itkIndex2b); + + CPPUNIT_ASSERT( + ((itkIndex2b[0] == itkIndex2[0]) && + (itkIndex2b[1] == itkIndex2[1]))); + //Testing itk::index<2> for IndexToWorld/WorldToIndex consistency + + dummyGeometry->IndexToWorld(itkIndex3,point); + dummyGeometry->WorldToIndex(point,itkIndex3b); + + CPPUNIT_ASSERT( + ((itkIndex3b[0] == itkIndex3[0]) && + (itkIndex3b[1] == itkIndex3[1]) && + (itkIndex3b[2] == itkIndex3[2]))); + //Testing itk::index<3> for IndexToWorld/WorldToIndex consistency + + dummyGeometry->IndexToWorld(itkIndex4,point); + dummyGeometry->WorldToIndex(point,itkIndex4b); + + CPPUNIT_ASSERT( + ((itkIndex4b[0] == itkIndex4[0]) && + (itkIndex4b[1] == itkIndex4[1]) && + (itkIndex4b[2] == itkIndex4[2]) && + (itkIndex4b[3] == 0))); + //Testing itk::index<3> for IndexToWorld/WorldToIndex consistency + + dummyGeometry->IndexToWorld(mitkIndex,point); + dummyGeometry->WorldToIndex(point,mitkIndexb); + + CPPUNIT_ASSERT( + ((mitkIndexb[0] == mitkIndex[0]) && + (mitkIndexb[1] == mitkIndex[1]) && + (mitkIndexb[2] == mitkIndex[2]))); + //Testing mitk::Index for IndexToWorld/WorldToIndex consistency + + return EXIT_SUCCESS; + } + + void TestIndexToWorld(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + + testIndexAndWorldConsistency(dummy); + testIndexAndWorldConsistencyForVectors(dummy); + testIndexAndWorldConsistencyForIndex(dummy); + + //Geometry must not have changed + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + + //Test with other geometries + dummy->SetOrigin(anotherPoint); + testIndexAndWorldConsistency(dummy); + testIndexAndWorldConsistencyForVectors(dummy); + testIndexAndWorldConsistencyForIndex(dummy); + + dummy->SetIndexToWorldTransform(anotherTransform); + testIndexAndWorldConsistency(dummy); + testIndexAndWorldConsistencyForVectors(dummy); + testIndexAndWorldConsistencyForIndex(dummy); + + dummy->SetOrigin(anotherPoint); + testIndexAndWorldConsistency(dummy); + testIndexAndWorldConsistencyForVectors(dummy); + testIndexAndWorldConsistencyForIndex(dummy); + + dummy->SetSpacing(anotherSpacing); + testIndexAndWorldConsistency(dummy); + testIndexAndWorldConsistencyForVectors(dummy); + testIndexAndWorldConsistencyForIndex(dummy); + } + + void TestExecuteOperation(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + + //Do same Operations with new Dummy and compare + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + + //Test operation Nothing + mitk::Operation* opN = new mitk::Operation(mitk::OpNOTHING); + dummy->ExecuteOperation(opN); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + + //Test operation Move + mitk::PointOperation* opP = new mitk::PointOperation(mitk::OpMOVE,anotherPoint); + dummy->ExecuteOperation(opP); + CPPUNIT_ASSERT(anotherPoint==dummy->GetOrigin()); + newDummy->SetOrigin(anotherPoint); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + + //Test operation Scale, Scale sets spacing to scale+1 + mitk::Point3D spacing; + spacing[0]=anotherSpacing[0]-1.; + spacing[1]=anotherSpacing[1]-1.; + spacing[2]=anotherSpacing[2]-1.; + + mitk::PointOperation* opS = new mitk::PointOperation(mitk::OpSCALE,spacing); + dummy->ExecuteOperation(opS); + CPPUNIT_ASSERT(anotherSpacing==dummy->GetSpacing()); + newDummy->SetSpacing(anotherSpacing); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + + //change Geometry to test more cases + dummy->SetIndexToWorldTransform(anotherTransform); + dummy->SetSpacing(anotherSpacing); + + //Testing a rotation of the geometry + double angle = 35.0; + mitk::Vector3D rotationVector; mitk::FillVector3D( rotationVector, 1, 0, 0 ); + mitk::Point3D center = dummy->GetCenter(); + mitk::RotationOperation* opR = new mitk::RotationOperation( mitk::OpROTATE, center, rotationVector, angle ); + dummy->ExecuteOperation(opR); + + mitk::Matrix3D rotation; + mitk::GetRotation(dummy, rotation); + mitk::Vector3D voxelStep=rotation*anotherSpacing; + mitk::Vector3D voxelStepIndex; + dummy->WorldToIndex(voxelStep, voxelStepIndex); + mitk::Vector3D expectedVoxelStepIndex; + expectedVoxelStepIndex.Fill(1); + CPPUNIT_ASSERT(mitk::Equal(voxelStepIndex,expectedVoxelStepIndex)); + + delete opR, opN, opS, opP; + } + + void TestCalculateBoundingBoxRelToTransform(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetExtentInMM(0,15); + dummy->SetExtentInMM(1,20); + dummy->SetExtentInMM(2,8); + + mitk::BoundingBox::Pointer dummyBoundingBox = dummy->CalculateBoundingBoxRelativeToTransform(anotherTransform); + + mitk::BoundingBox::PointsContainer::Pointer pointscontainer=mitk::BoundingBox::PointsContainer::New(); + mitk::BoundingBox::PointIdentifier pointid=0; + unsigned char i; + mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New(); + anotherTransform->GetInverse(inverse); + for(i=0; i<8; ++i) + pointscontainer->InsertElement( pointid++, inverse->TransformPoint( dummy->GetCornerPoint(i) )); + mitk::BoundingBox::Pointer result = mitk::BoundingBox::New(); + result->SetPoints(pointscontainer); + result->ComputeBoundingBox(); + + CPPUNIT_ASSERT(mitk::Equal(result,dummyBoundingBox,mitk::eps,true)); + + //dummy still needs to be unchanged, except for extend + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + newDummy->SetExtentInMM(0,15); + newDummy->SetExtentInMM(1,20); + newDummy->SetExtentInMM(2,8); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + void TestSetTimeBounds(){ + mitk::TimeBounds timeBounds; + timeBounds[0] = 1; + timeBounds[1] = 9; + + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetTimeBounds(timeBounds); + mitk::TimeBounds timeBounds2 = dummy->GetTimeBounds(); + + CPPUNIT_ASSERT(timeBounds[0]==timeBounds2[0]); + CPPUNIT_ASSERT(timeBounds[1]==timeBounds2[1]); + + //undo changes, new and changed object need to be the same! + timeBounds[0]=mitk::ScalarTypeNumericTraits::NonpositiveMin(); + timeBounds[1]=mitk::ScalarTypeNumericTraits::max(); + + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + void TestIs2DConvertable(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + + //new initialized geometry is 2D convertable + CPPUNIT_ASSERT(dummy->Is2DConvertable()); + + //Wrong Spacing needs to fail + dummy->SetSpacing(anotherSpacing); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + //undo + dummy->SetSpacing(aSpacing); + CPPUNIT_ASSERT(dummy->Is2DConvertable()); + + //Wrong Origin needs to fail + dummy->SetOrigin(anotherPoint); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + //undo + dummy->SetOrigin(aPoint); + CPPUNIT_ASSERT(dummy->Is2DConvertable()); + + //third dimension must not be transformed + mitk::AffineTransform3D::Pointer dummyTransform = mitk::AffineTransform3D::New(); + mitk::AffineTransform3D::MatrixType dummyMatrix; + dummyMatrix.SetIdentity(); + dummyTransform->SetMatrix( dummyMatrix ); + dummy->SetIndexToWorldTransform(dummyTransform); + + //identity matrix is 2DConvertable + CPPUNIT_ASSERT(dummy->Is2DConvertable()); + + dummyMatrix(0,2) = 3; + dummyTransform->SetMatrix( dummyMatrix ); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + + dummyMatrix.SetIdentity(); + dummyMatrix(1,2) = 0.4; + dummyTransform->SetMatrix( dummyMatrix ); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + + dummyMatrix.SetIdentity(); + dummyMatrix(2,2) = 3; + dummyTransform->SetMatrix( dummyMatrix ); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + + dummyMatrix.SetIdentity(); + dummyMatrix(2,1) = 3; + dummyTransform->SetMatrix( dummyMatrix ); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + + dummyMatrix.SetIdentity(); + dummyMatrix(2,0) = 3; + dummyTransform->SetMatrix( dummyMatrix ); + CPPUNIT_ASSERT(dummy->Is2DConvertable()==false); + + //undo changes, new and changed object need to be the same! + dummyMatrix.SetIdentity(); + dummyTransform->SetMatrix( dummyMatrix ); + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + void TestGetCornerPoint(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetIndexToWorldTransform(anotherTransform); + double bounds[6] = {0,11,0,12,0,13}; + dummy->SetFloatBounds(bounds); + + mitk::Point3D corner, refCorner; + + //Corner 0 + mitk::FillVector3D(refCorner,bounds[0],bounds[2],bounds[4]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(0); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(true,true,true); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 1 + mitk::FillVector3D(refCorner,bounds[0],bounds[2],bounds[5]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(1); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(true,true,false); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 2 + mitk::FillVector3D(refCorner,bounds[0],bounds[3],bounds[4]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(2); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(true,false,true); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 3 + mitk::FillVector3D(refCorner,bounds[0],bounds[3],bounds[5]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(3); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(true,false,false); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 4 + mitk::FillVector3D(refCorner,bounds[1],bounds[2],bounds[4]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(4); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(false,true,true); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 5 + mitk::FillVector3D(refCorner,bounds[1],bounds[2],bounds[5]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(5); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(false,true,false); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 6 + mitk::FillVector3D(refCorner,bounds[1],bounds[3],bounds[4]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(6); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(false,false,true); + CPPUNIT_ASSERT(refCorner==corner); + + //Corner 7 + mitk::FillVector3D(refCorner,bounds[1],bounds[3],bounds[5]); + refCorner = anotherTransform->TransformPoint(refCorner); + corner=dummy->GetCornerPoint(7); + CPPUNIT_ASSERT(refCorner==corner); + corner=dummy->GetCornerPoint(false,false,false); + CPPUNIT_ASSERT(refCorner==corner); + + //Wrong Corner needs to fail + // CPPUNIT_ASSERT_THROW(dummy->GetCornerPoint(-1),itk::ExceptionObject); + + //dummy geometry must not have changed! + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + newDummy->SetIndexToWorldTransform(anotherTransform); + newDummy->SetFloatBounds(bounds); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + void TestExtentInMM() + { + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetExtentInMM(0,50); + CPPUNIT_ASSERT(50==dummy->GetExtentInMM(0)); + //Vnl Matrix has changed. The next line only works because the spacing is 1! + CPPUNIT_ASSERT(50==dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0).magnitude()); + + //Smaller extent than original + dummy->SetExtentInMM(0,5); + CPPUNIT_ASSERT(5==dummy->GetExtentInMM(0)); + CPPUNIT_ASSERT(5==dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0).magnitude()); + + dummy->SetExtentInMM(1,4); + CPPUNIT_ASSERT(4==dummy->GetExtentInMM(1)); + CPPUNIT_ASSERT(4==dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1).magnitude()); + + dummy->SetExtentInMM(2,2.5); + CPPUNIT_ASSERT(2.5==dummy->GetExtentInMM(2)); + CPPUNIT_ASSERT(2.5==dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2).magnitude()); + } + + void TestGetAxisVector(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetIndexToWorldTransform(anotherTransform); + double bounds[6] = {0,11,0,12,0,13}; + dummy->SetFloatBounds(bounds); + + mitk::Vector3D vector; + mitk::FillVector3D(vector,bounds[1],0,0); + dummy->IndexToWorld(vector,vector); + CPPUNIT_ASSERT(dummy->GetAxisVector(0)==vector); + + mitk::FillVector3D(vector,0,bounds[3],0); + dummy->IndexToWorld(vector,vector); + CPPUNIT_ASSERT(dummy->GetAxisVector(1)==vector); + + mitk::FillVector3D(vector,0,0,bounds[5]); + dummy->IndexToWorld(vector,vector); + CPPUNIT_ASSERT(dummy->GetAxisVector(2)==vector); + } + + void TestGetCenter(){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + dummy->SetIndexToWorldTransform(anotherTransform); + double bounds[6] = {0,11,2,12,1,13}; + dummy->SetFloatBounds(bounds); + + mitk::Point3D refCenter; + for( int i=0;i<3;i++) + refCenter.SetElement(i,( bounds[2 * i] + bounds[2 * i + 1] ) / 2.0); + + dummy->IndexToWorld(refCenter,refCenter); + + CPPUNIT_ASSERT(dummy->GetCenter()==refCenter); + } + + /* + + void (){ + DummyTestClass::Pointer dummy = DummyTestClass::New(); + + CPPUNIT_ASSERT(); + + //undo changes, new and changed object need to be the same! + + DummyTestClass::Pointer newDummy = DummyTestClass::New(); + CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); + } + + */ };//end class mitkBaseGeometryTestSuite MITK_TEST_SUITE_REGISTRATION(mitkBaseGeometry)