diff --git a/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp b/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp index db80da8ae6..5a79c7edc9 100644 --- a/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp +++ b/Core/Code/DataManagement/mitkAbstractTransformGeometry.cpp @@ -1,318 +1,318 @@ /*=================================================================== 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) +mitk::AbstractTransformGeometry::AbstractTransformGeometry() : Superclass(), m_Plane(NULL), m_FrameGeometry(NULL) { Initialize(); } mitk::AbstractTransformGeometry::AbstractTransformGeometry(const AbstractTransformGeometry& other) : Superclass(other), m_ParametricBoundingBox(other.m_ParametricBoundingBox) { if(other.m_ParametricBoundingBox.IsNotNull()) { m_ParametricBoundingBox = other.m_ParametricBoundingBox->DeepCopy(); this->SetParametricBounds(m_ParametricBoundingBox->GetBounds()); } this->SetPlane(other.m_Plane); this->SetFrameGeometry(other.m_FrameGeometry); } mitk::AbstractTransformGeometry::~AbstractTransformGeometry() { } void mitk::AbstractTransformGeometry::PostInitialize() { 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 itk::Transform* mitk::AbstractTransformGeometry::GetParametricTransform() const { return m_ItkVtkAbstractTransform; } bool mitk::AbstractTransformGeometry::Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const { 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(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(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, bool considerBoundingBox) const { assert((m_ItkVtkAbstractTransform.IsNotNull()) && (m_Plane.IsNotNull())); Point3D pt3d_ParametricWorld; pt3d_ParametricWorld = m_ItkVtkAbstractTransform->BackTransform(pt3d_mm); Point3D pt3d_ParametricUnits; ((BaseGeometry*)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::BaseGeometry* 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(); } void mitk::AbstractTransformGeometry::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(); } const mitk::BoundingBox::BoundsArrayType& mitk::AbstractTransformGeometry::GetParametricBounds() const { assert(m_ParametricBoundingBox.IsNotNull()); return m_ParametricBoundingBox->GetBounds(); } mitk::ScalarType mitk::AbstractTransformGeometry::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]; } diff --git a/Core/Code/DataManagement/mitkBaseGeometry.cpp b/Core/Code/DataManagement/mitkBaseGeometry.cpp index 4f10cb7607..fb51ca75e3 100644 --- a/Core/Code/DataManagement/mitkBaseGeometry.cpp +++ b/Core/Code/DataManagement/mitkBaseGeometry.cpp @@ -1,1095 +1,1086 @@ /*=================================================================== 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" #include "mitkModifiedLock.h" mitk::BaseGeometry::BaseGeometry(): Superclass(), mitk::OperationActor(), m_FrameOfReferenceID(0), m_IndexToWorldTransformLastModified(0), m_ImageGeometry(false), m_ModifiedLockFlag(false), m_ModifiedCalledFlag(false) { m_VtkMatrix = vtkMatrix4x4::New(); m_VtkIndexToWorldTransform = vtkMatrixToLinearTransform::New(); m_VtkIndexToWorldTransform->SetInput(m_VtkMatrix); Initialize(); } mitk::BaseGeometry::BaseGeometry(const BaseGeometry& other): Superclass(), mitk::OperationActor(), //m_TimeBounds(other.m_TimeBounds), m_FrameOfReferenceID(other.m_FrameOfReferenceID), m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified), m_Origin(other.m_Origin), m_ImageGeometry(other.m_ImageGeometry), m_ModifiedLockFlag(false), m_ModifiedCalledFlag(false) { - // DEPRECATED(m_RotationQuaternion = other.m_RotationQuaternion); - // AffineGeometryFrame - SetBounds(other.GetBounds()); m_VtkMatrix = vtkMatrix4x4::New(); - m_VtkMatrix->DeepCopy(other.m_VtkMatrix); - 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) { mitk::ModifiedLock lock(this); if(origin!=GetOrigin()) { m_Origin = origin; m_IndexToWorldTransform->SetOffset(m_Origin.GetVectorFromOrigin()); Modified(); TransferItkToVtkTransform(); } } void mitk::BaseGeometry::TransferItkToVtkTransform() { mitk::ModifiedLock lock(this); // copy m_IndexToWorldTransform into m_VtkIndexToWorldTransform TransferItkTransformToVtkMatrix(m_IndexToWorldTransform.GetPointer(), m_VtkMatrix); m_VtkIndexToWorldTransform->Modified(); } static void CopySpacingFromTransform(mitk::AffineTransform3D* transform, mitk::Vector3D& spacing) { 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(); } 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); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); m_VtkMatrix->Identity(); //m_TimeBounds[0]=ScalarTypeNumericTraits::NonpositiveMin(); m_TimeBounds[1]=ScalarTypeNumericTraits::max(); m_FrameOfReferenceID = 0; m_ImageGeometry = false; this->PostInitialize(); } void mitk::BaseGeometry::PostInitializeGeometry(BaseGeometry * newGeometry) const { newGeometry->m_ImageGeometry = m_ImageGeometry; } 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() ); + TransformType::Pointer indexToWorldTransform = m_IndexToWorldTransform->Clone(); newGeometry->SetIndexToWorldTransform(indexToWorldTransform); } this->PostInitializeGeometry(newGeometry); } void mitk::BaseGeometry::PostInitialize() { } /** Set the bounds */ void mitk::BaseGeometry::SetBounds(const BoundsArrayType& bounds) { mitk::ModifiedLock lock(this); PreSetBounds(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::PreSetBounds(const BoundsArrayType& /*bounds*/){}; void mitk::BaseGeometry::SetIndexToWorldTransform(mitk::AffineTransform3D* transform) { mitk::ModifiedLock lock(this); PreSetIndexToWorldTransform(transform); if(m_IndexToWorldTransform.GetPointer() != transform) { m_IndexToWorldTransform = transform; CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); TransferItkToVtkTransform(); Modified(); } PostSetIndexToWorldTransform(transform); } void mitk::BaseGeometry::PreSetIndexToWorldTransform(mitk::AffineTransform3D* /*transform*/) {} void mitk::BaseGeometry::PostSetIndexToWorldTransform(mitk::AffineTransform3D* /*transform*/) {} const mitk::BaseGeometry::BoundsArrayType mitk::BaseGeometry::GetBounds() const { assert(m_BoundingBox.IsNotNull()); return m_BoundingBox->GetBounds(); } bool mitk::BaseGeometry::IsValid() const { return true; } void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D& aSpacing, bool enforceSetSpacing ) { PreSetSpacing(aSpacing); _SetSpacing(aSpacing, enforceSetSpacing); } void mitk::BaseGeometry::PreSetSpacing(const mitk::Vector3D& /*aSpacing*/) {} void mitk::BaseGeometry::_SetSpacing(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()); } } 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>=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(this->IsBoundingBoxNull()==false); 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 this->GetIndexToWorldTransform()->TransformPoint(cornerpoint); } mitk::Point3D mitk::BaseGeometry::GetCornerPoint(bool xFront, bool yFront, bool zFront) const { 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 this->GetIndexToWorldTransform()->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) { mitk::ModifiedLock lock(this); 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(); } PostSetExtentInMM(direction,extentInMM); } void mitk::BaseGeometry::PostSetExtentInMM(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) 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 = 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 = this->GetBoundingBox()->GetBounds(); if((discretIndex[0] == bounds[1]) || (discretIndex[1] == bounds[3]) || (discretIndex[2] == bounds[5])) inside = false; } } else inside = this->GetBoundingBox()->IsInside(index); return inside; } 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() { mitk::ModifiedLock lock(this); m_IndexToWorldTransform->SetIdentity(); m_Origin.Fill(0); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing); Modified(); TransferItkToVtkTransform(); } void mitk::BaseGeometry::TransferVtkToItkTransform() { TransferVtkMatrixToItkTransform(m_VtkMatrix, m_IndexToWorldTransform.GetPointer()); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); } void mitk::BaseGeometry::Compose( const mitk::BaseGeometry::TransformType * other, bool pre ) { mitk::ModifiedLock lock(this); m_IndexToWorldTransform->Compose(other, pre); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing); 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) { mitk::ModifiedLock lock(this); 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) //{ // mitk::ModifiedLock lock(this); // // if(m_TimeBounds != timebounds) // { // m_TimeBounds = timebounds; // Modified(); // } // PostSetTimeBounds(timebounds); //} // //void mitk::BaseGeometry::PostSetTimeBounds(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(); } void mitk::BaseGeometry::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); } //itk::LightObject::Pointer mitk::BaseGeometry::InternalClone() const //{ // Self::Pointer newGeometry = new Self(*this); // newGeometry->UnRegister(); // return newGeometry.GetPointer(); //} void mitk::BaseGeometry::PrintSelf(std::ostream& os, itk::Indent indent) const { os << indent << " IndexToWorldTransform: "; 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 << this->GetIndexToWorldTransform()->GetMatrix()[i][j] << " "; } os << 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 << this->GetIndexToWorldTransform()->GetInverseMatrix()[i][j] << " "; } os << std::endl; } // from itk::ScalableAffineTransform os << indent << "Scale : "; for (i = 0; i < 3; i++) { os << this->GetIndexToWorldTransform()->GetScale()[i] << " "; } os << std::endl; } os << indent << " BoundingBox: "; if(this->IsBoundingBoxNull()) os << "NULL" << std::endl; else { os << indent << "( "; for (unsigned int i=0; i<3; i++) { os << this->GetBoundingBox()->GetBounds()[2*i] << "," << this->GetBoundingBox()->GetBounds()[2*i+1] << " "; } os << " )" << std::endl; } os << indent << " Origin: " << this->GetOrigin() << std::endl; os << indent << " ImageGeometry: " << this->GetImageGeometry() << std::endl; os << indent << " Spacing: " << this->GetSpacing() << std::endl; //os << indent << " TimeBounds: " << this->GetTimeBounds() << std::endl; } void mitk::BaseGeometry::Modified() const{ if(!m_ModifiedLockFlag) Superclass::Modified(); else m_ModifiedCalledFlag = true; } bool mitk::Equal( const mitk::BaseGeometry::BoundingBoxType *leftHandSide, const mitk::BaseGeometry::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) { if(( leftHandSide == NULL) || ( rightHandSide == NULL )) { MITK_ERROR << "mitk::Equal( const mitk::Geometry3D::BoundingBoxType *leftHandSide, const mitk::Geometry3D::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) does not with NULL pointer input."; return false; } return Equal( *leftHandSide, *rightHandSide, eps, verbose); } bool mitk::Equal( const mitk::BaseGeometry::BoundingBoxType& leftHandSide, const mitk::BaseGeometry::BoundingBoxType& rightHandSide, ScalarType eps, bool verbose ) { bool result = true; 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 << "[( Geometry3D::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) { if(( leftHandSide == NULL) || ( rightHandSide == NULL )) { MITK_ERROR << "mitk::Equal(const mitk::Geometry3D *leftHandSide, const mitk::Geometry3D *rightHandSide, ScalarType eps, bool verbose) does not with NULL pointer input."; return false; } return Equal( *leftHandSide, *rightHandSide, eps, verbose); } bool mitk::Equal(const mitk::BaseGeometry& leftHandSide, const mitk::BaseGeometry& rightHandSide, ScalarType eps, bool verbose) { bool result = true; //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; } bool mitk::Equal(const BaseGeometry::TransformType *leftHandSide, const BaseGeometry::TransformType *rightHandSide, ScalarType eps, bool verbose ) { if(( leftHandSide == NULL) || ( rightHandSide == NULL )) { MITK_ERROR << "mitk::Equal(const Geometry3D::TransformType *leftHandSide, const Geometry3D::TransformType *rightHandSide, ScalarType eps, bool verbose ) does not with NULL pointer input."; return false; } return Equal( *leftHandSide, *rightHandSide, eps, verbose); } 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 << "[( Geometry3D::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; } diff --git a/Core/Code/DataManagement/mitkDisplayGeometry.cpp b/Core/Code/DataManagement/mitkDisplayGeometry.cpp index ab93699403..7ad44c83e7 100644 --- a/Core/Code/DataManagement/mitkDisplayGeometry.cpp +++ b/Core/Code/DataManagement/mitkDisplayGeometry.cpp @@ -1,613 +1,614 @@ /*=================================================================== 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::IsValid() const { return m_WorldGeometry.IsNotNull() && m_WorldGeometry->IsValid(); } unsigned long mitk::DisplayGeometry::GetMTime() const { if((m_WorldGeometry.IsNotNull()) && (PlaneGeometry::GetMTime() < m_WorldGeometry->GetMTime())) { Modified(); } return PlaneGeometry::GetMTime(); } //const mitk::TimeBounds& mitk::DisplayGeometry::GetTimeBounds() const //{ // if(m_WorldGeometry.IsNull()) // { // return m_TimeBounds; // } // // return m_WorldGeometry->GetTimeBounds(); //} // size definition methods void mitk::DisplayGeometry::SetWorldGeometry(const PlaneGeometry* 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-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 < itk::NumericTraits::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->GetExtentInMM(0); const ScalarType& heightInMM = m_WorldGeometry->GetExtentInMM(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) + : PlaneGeometry() + ,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->GetExtentInMM(0); float worldHeightMM = m_WorldGeometry->GetExtentInMM(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/mitkLandmarkProjectorBasedCurvedGeometry.cpp b/Core/Code/DataManagement/mitkLandmarkProjectorBasedCurvedGeometry.cpp index 9e8ac72940..7aa5dc7ef1 100644 --- a/Core/Code/DataManagement/mitkLandmarkProjectorBasedCurvedGeometry.cpp +++ b/Core/Code/DataManagement/mitkLandmarkProjectorBasedCurvedGeometry.cpp @@ -1,83 +1,83 @@ /*=================================================================== 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 "mitkLandmarkProjectorBasedCurvedGeometry.h" #include mitk::LandmarkProjectorBasedCurvedGeometry::LandmarkProjectorBasedCurvedGeometry() - : m_LandmarkProjector(NULL), m_InterpolatingAbstractTransform(NULL), m_TargetLandmarks(NULL) + : AbstractTransformGeometry(), m_LandmarkProjector(NULL), m_InterpolatingAbstractTransform(NULL), m_TargetLandmarks(NULL) { } mitk::LandmarkProjectorBasedCurvedGeometry::LandmarkProjectorBasedCurvedGeometry(const mitk::LandmarkProjectorBasedCurvedGeometry& other) : Superclass(other) { SetTargetLandmarks(other.m_TargetLandmarks); this->SetLandmarkProjector(other.m_LandmarkProjector); this->ComputeGeometry(); } mitk::LandmarkProjectorBasedCurvedGeometry::~LandmarkProjectorBasedCurvedGeometry() { if(m_InterpolatingAbstractTransform!=NULL) m_InterpolatingAbstractTransform->Delete(); } void mitk::LandmarkProjectorBasedCurvedGeometry::SetLandmarkProjector(mitk::LandmarkProjector* aLandmarkProjector) { itkDebugMacro("setting LandmarkProjector to " << aLandmarkProjector ); if(m_LandmarkProjector != aLandmarkProjector) { m_LandmarkProjector = aLandmarkProjector; if(m_LandmarkProjector.IsNotNull()) { if(m_FrameGeometry.IsNotNull()) m_LandmarkProjector->SetFrameGeometry(m_FrameGeometry); if(m_InterpolatingAbstractTransform == NULL) { itkWarningMacro(<<"m_InterpolatingAbstractTransform not set."); } m_LandmarkProjector->SetInterpolatingAbstractTransform(GetInterpolatingAbstractTransform()); SetVtkAbstractTransform(m_LandmarkProjector->GetCompleteAbstractTransform()); } Modified(); } } void mitk::LandmarkProjectorBasedCurvedGeometry::SetFrameGeometry(const mitk::BaseGeometry* frameGeometry) { Superclass::SetFrameGeometry(frameGeometry); if(m_LandmarkProjector.IsNotNull()) m_LandmarkProjector->SetFrameGeometry(frameGeometry); } void mitk::LandmarkProjectorBasedCurvedGeometry::ComputeGeometry() { if(m_LandmarkProjector.IsNull()) { itkExceptionMacro(<< "m_LandmarkProjector is not set."); } m_LandmarkProjector->ProjectLandmarks(m_TargetLandmarks); SetPlane(m_LandmarkProjector->GetParameterPlane()); } itk::LightObject::Pointer mitk::LandmarkProjectorBasedCurvedGeometry::InternalClone() const { mitk::BaseGeometry::Pointer newGeometry = new LandmarkProjectorBasedCurvedGeometry(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } diff --git a/Core/Code/DataManagement/mitkPlaneGeometry.cpp b/Core/Code/DataManagement/mitkPlaneGeometry.cpp index 8ad0309caf..eca60654bd 100644 --- a/Core/Code/DataManagement/mitkPlaneGeometry.cpp +++ b/Core/Code/DataManagement/mitkPlaneGeometry.cpp @@ -1,916 +1,916 @@ /*=================================================================== 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() - : m_ScaleFactorMMPerUnitX( 1.0 ), + : Superclass(), m_ScaleFactorMMPerUnitX( 1.0 ), m_ScaleFactorMMPerUnitY( 1.0 ), m_ReferenceGeometry( NULL ) { Initialize(); } PlaneGeometry::~PlaneGeometry() { } PlaneGeometry::PlaneGeometry(const PlaneGeometry& other) : Superclass(other), m_ScaleFactorMMPerUnitX( other.m_ScaleFactorMMPerUnitX), m_ScaleFactorMMPerUnitY( other.m_ScaleFactorMMPerUnitY), m_ReferenceGeometry( other.m_ReferenceGeometry ) { } 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::PreSetIndexToWorldTransform(mitk::AffineTransform3D *transform) { EnsurePerpendicularNormal(transform); } void PlaneGeometry::PreSetBounds(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 BaseGeometry *geometry3D, PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { this->SetReferenceGeometry( const_cast< BaseGeometry * >( 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 BaseGeometry *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(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(this->GetIndexToWorldTransform()->GetOffset()); SetIndexToWorldTransform(transform); } Vector3D PlaneGeometry::GetNormal() const { Vector3D frontToBack; frontToBack.SetVnlVector( this->GetIndexToWorldTransform() ->GetMatrix().GetVnlMatrix().get_column(2) ); return frontToBack; } VnlVector PlaneGeometry::GetNormalVnl() const { 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); } //Function from Geometry2D // mitk::ScalarType // PlaneGeometry::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 PlaneGeometry::IsAbove( const Point3D &pt3d_mm , bool considerBoundingBox) const { if(considerBoundingBox) { Point3D pt3d_units; BaseGeometry::WorldToIndex(pt3d_mm, pt3d_units); return (pt3d_units[2] > this->GetBoundingBox()->GetBounds()[4]); } else 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( 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; } 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 << " ScaleFactorMMPerUnitX: " << m_ScaleFactorMMPerUnitX << std::endl; os << indent << " ScaleFactorMMPerUnitY: " << m_ScaleFactorMMPerUnitY << std::endl; os << indent << " Normal: " << GetNormal() << std::endl; } void PlaneGeometry::PostSetIndexToWorldTransform( mitk::AffineTransform3D* transform) { m_ScaleFactorMMPerUnitX=GetExtentInMM(0)/GetExtent(0); m_ScaleFactorMMPerUnitY=GetExtentInMM(1)/GetExtent(1); assert(m_ScaleFactorMMPerUnitX::infinity()); assert(m_ScaleFactorMMPerUnitY::infinity()); } void PlaneGeometry::PostSetExtentInMM(int direction, ScalarType extentInMM) { m_ScaleFactorMMPerUnitX=GetExtentInMM(0)/GetExtent(0); m_ScaleFactorMMPerUnitY=GetExtentInMM(1)/GetExtent(1); assert(m_ScaleFactorMMPerUnitX::infinity()); assert(m_ScaleFactorMMPerUnitY::infinity()); } bool PlaneGeometry::Map( const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const { assert(this->IsBoundingBoxNull()==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(this->GetBoundingBox())->IsInside(pt3d_units); } void PlaneGeometry::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 = GetIndexToWorldTransform()->TransformPoint(pt3d_units); } void PlaneGeometry::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 PlaneGeometry::Project( const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const { assert(this->IsBoundingBoxNull()==false); Point3D pt3d_units; BackTransform(pt3d_mm, pt3d_units); pt3d_units[2] = 0; projectedPt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units); return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool PlaneGeometry::Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; BackTransform(vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units); return true; } bool PlaneGeometry::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(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; BackTransform(atPt3d_mm, vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units); Point3D pt3d_units; BackTransform(atPt3d_mm, pt3d_units); return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool PlaneGeometry::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 PlaneGeometry::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); } void PlaneGeometry::SetReferenceGeometry( mitk::BaseGeometry *geometry ) { m_ReferenceGeometry = geometry; } mitk::BaseGeometry * PlaneGeometry::GetReferenceGeometry() const { return m_ReferenceGeometry; } bool PlaneGeometry::HasReferenceGeometry() const { return ( m_ReferenceGeometry != NULL ); } } // namespace diff --git a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp index 562f2d6745..33ddbac8d7 100644 --- a/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp +++ b/Core/Code/DataManagement/mitkThinPlateSplineCurvedGeometry.cpp @@ -1,100 +1,101 @@ /*=================================================================== 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() + : Superclass() { 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::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::BaseGeometry::Pointer newGeometry = new Self(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } diff --git a/Core/Code/Testing/mitkImageTest.cpp b/Core/Code/Testing/mitkImageTest.cpp index 978fee8647..545ceb6db8 100644 --- a/Core/Code/Testing/mitkImageTest.cpp +++ b/Core/Code/Testing/mitkImageTest.cpp @@ -1,526 +1,517 @@ /*=================================================================== 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. ===================================================================*/ // mitk includes #include #include #include #include "mitkItkImageFileReader.h" #include #include #include "mitkImageGenerator.h" #include "mitkImageReadAccessor.h" #include "mitkException.h" #include "mitkPixelTypeMultiplex.h" #include "mitkImagePixelReadAccessor.h" #include "mitkImageSliceSelector.h" // itk includes #include #include // stl includes #include // vtk includes #include // Checks if reference count is correct after using GetVtkImageData() bool ImageVtkDataReferenceCheck(const char* fname) { const std::string filename = std::string(fname); mitk::ItkImageFileReader::Pointer imageReader = mitk::ItkImageFileReader::New(); try { imageReader->SetFileName(filename); imageReader->Update(); } catch(...) { MITK_TEST_FAILED_MSG(<< "Could not read file for testing: " << filename); return false; } { mitk::Image::Pointer image = imageReader->GetOutput(); vtkImageData* vtk = image->GetVtkImageData(); if(vtk == NULL) return false; } return true; } template void TestRandomPixelAccess( const mitk::PixelType ptype, mitk::Image::Pointer image, mitk::Point3D & point, mitk::ScalarType & value ) { // generate a random point in world coordinates mitk::Point3D xMax, yMax, zMax, xMaxIndex, yMaxIndex, zMaxIndex; xMaxIndex.Fill(0.0f); yMaxIndex.Fill(0.0f); zMaxIndex.Fill(0.0f); xMaxIndex[0] = image->GetLargestPossibleRegion().GetSize()[0]; yMaxIndex[1] = image->GetLargestPossibleRegion().GetSize()[1]; zMaxIndex[2] = image->GetLargestPossibleRegion().GetSize()[2]; image->GetGeometry()->IndexToWorld(xMaxIndex, xMax); image->GetGeometry()->IndexToWorld(yMaxIndex, yMax); image->GetGeometry()->IndexToWorld(zMaxIndex, zMax); MITK_INFO << "Origin " << image->GetGeometry()->GetOrigin()[0] << " "<< image->GetGeometry()->GetOrigin()[1] << " "<< image->GetGeometry()->GetOrigin()[2] << ""; MITK_INFO << "MaxExtend " << xMax[0] << " "<< yMax[1] << " "<< zMax[2] << ""; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randomGenerator = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randomGenerator->Initialize( std::rand() ); // initialize with random value, to get sensible random points for the image point[0] = randomGenerator->GetUniformVariate( image->GetGeometry()->GetOrigin()[0], xMax[0]); point[1] = randomGenerator->GetUniformVariate( image->GetGeometry()->GetOrigin()[1], yMax[1]); point[2] = randomGenerator->GetUniformVariate( image->GetGeometry()->GetOrigin()[2], zMax[2]); MITK_INFO << "RandomPoint " << point[0] << " "<< point[1] << " "<< point[2] << ""; // test values and max/min mitk::ScalarType imageMin = image->GetStatistics()->GetScalarValueMin(); mitk::ScalarType imageMax = image->GetStatistics()->GetScalarValueMax(); // test accessing PixelValue with coordinate leading to a negative index const mitk::Point3D geom_origin = image->GetGeometry()->GetOrigin(); const mitk::Point3D geom_center = image->GetGeometry()->GetCenter(); // shift position from origin outside of the image ( in the opposite direction to [center-origin] vector which points in the inside) mitk::Point3D position = geom_origin + (geom_origin - geom_center); MITK_INFO << "Testing access outside of the image"; unsigned int dim = image->GetDimension(); if(dim == 3 || dim == 4){ mitk::ImagePixelReadAccessor imAccess3(image,image->GetVolumeData(0)); // Comparison ?>=0 not needed since all position[i] and timestep are unsigned int // (position[0]>=0 && position[1] >=0 && position[2]>=0 && timestep>=0) // bug-11978 : we still need to catch index with negative values if ( point[0] < 0 || point[1] < 0 || point[2] < 0 ) { MITK_WARN << "Given position ("<< point << ") is out of image range, returning 0." ; } else { value = static_cast(imAccess3.GetPixelByWorldCoordinates(point)); MITK_TEST_CONDITION( (value >= imageMin && value <= imageMax), "Value returned is between max/min"); } itk::Index<3> itkIndex; image->GetGeometry()->WorldToIndex(position, itkIndex); MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception); imAccess3.GetPixelByIndexSafe(itkIndex); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception); } MITK_INFO << imageMin << " "<< imageMax << " "<< value << ""; } class mitkImageTestClass { public: void SetClonedGeometry_None_ClonedEqualInput() { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateRandomImage(100, 100, 100, 1, 0.2, 0.3, 0.4); //----------------- // geometry information for image mitk::Point3D origin; mitk::Vector3D right, bottom; mitk::Vector3D spacing; mitk::FillVector3D(origin, 17.0, 19.92, 7.83); mitk::FillVector3D(right, 1.0, 2.0, 3.0); mitk::FillVector3D(bottom, 0.0, -3.0, 2.0); mitk::FillVector3D(spacing, 0.78, 0.91, 2.23); //InitializeStandardPlane(rightVector, downVector, spacing) mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); planegeometry->InitializeStandardPlane(100, 100, right, bottom, &spacing); planegeometry->SetOrigin(origin); planegeometry->ChangeImageGeometryConsideringOriginOffset(true); image->SetClonedGeometry(planegeometry); mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); itk::ScalableAffineTransform* frameNew = imageGeometry->GetIndexToWorldTransform(); itk::ScalableAffineTransform* frameOld = planegeometry->GetIndexToWorldTransform(); - bool matrixEqual = true; - for (int i = 0; i < 16; ++i) - { - double valueNew = *(frameNew->GetMatrix()[i]); - double valueOld = *(frameOld->GetMatrix()[i]); - //MITK_INFO << "Index: " << i << " Old: " << valueOld << " New: " << valueNew << " Difference:" << valueOld-valueNew<< std::endl; - matrixEqual = matrixEqual && mitk::Equal(valueNew, valueOld, mitk::eps); - } + bool matrixEqual = mitk::Equal(imageGeometry, planegeometry, mitk::eps, false); - // Disabled because this test fails on the dashboard. Does not fail on my machine. - // See Bug 6505 - // MITK_TEST_CONDITION(matrixEqual, "Matrix elements of cloned matrix equal original matrix"); + MITK_TEST_CONDITION(matrixEqual, "Matrix elements of cloned matrix equal original matrix"); } }; int mitkImageTest(int argc, char* argv[]) { MITK_TEST_BEGIN(mitkImageTest); mitkImageTestClass tester; tester.SetClonedGeometry_None_ClonedEqualInput(); //Create Image out of nowhere mitk::Image::Pointer imgMem = mitk::Image::New(); mitk::PixelType pt = mitk::MakeScalarPixelType(); unsigned int dim[]={100,100,20}; MITK_TEST_CONDITION_REQUIRED( imgMem.IsNotNull(), "An image was created. "); // Initialize image imgMem->Initialize( pt, 3, dim); MITK_TEST_CONDITION_REQUIRED( imgMem->IsInitialized(), "Image::IsInitialized() ?"); MITK_TEST_CONDITION_REQUIRED( imgMem->GetPixelType() == pt, "PixelType was set correctly."); int *p = NULL; int *p2 = NULL; try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int*)imgMemAcc.GetData(); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION( p != NULL, "GetData() returned not-NULL pointer."); // filling image const unsigned int size = dim[0]*dim[1]*dim[2]; for(unsigned int i=0; iGetSliceData(dim[2]/2)); p2 = (int*)imgMemAcc.GetData(); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED( p2 != NULL, "Valid slice data returned"); unsigned int xy_size = dim[0]*dim[1]; unsigned int start_mid_slice = (dim[2]/2)*xy_size; isEqual = true; for(unsigned int i=0; i(); imgMem->Initialize( pType , 3, dim); MITK_TEST_CONDITION_REQUIRED(imgMem->GetDimension()== 3, "Testing initialization parameter dimension!"); MITK_TEST_CONDITION_REQUIRED(imgMem->GetPixelType() == pType, "Testing initialization parameter pixeltype!"); MITK_TEST_CONDITION_REQUIRED(imgMem->GetDimension(0) == dim[0] && imgMem->GetDimension(1)== dim[1] && imgMem->GetDimension(2)== dim[2], "Testing initialization of dimensions!"); MITK_TEST_CONDITION( imgMem->IsInitialized(), "Image is initialized."); // Setting volume again: try { mitk::ImageReadAccessor imgMemAcc(imgMem); imgMem->SetVolume(imgMemAcc.GetData()); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } //----------------- // geometry information for image mitk::Point3D origin; mitk::Vector3D right, bottom; mitk::Vector3D spacing; mitk::FillVector3D(origin, 17.0, 19.92, 7.83); mitk::FillVector3D(right, 1.0, 2.0, 3.0); mitk::FillVector3D(bottom, 0.0, -3.0, 2.0); mitk::FillVector3D(spacing, 0.78, 0.91, 2.23); //InitializeStandardPlane(rightVector, downVector, spacing) mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); planegeometry->InitializeStandardPlane(100, 100, right, bottom, &spacing); planegeometry->SetOrigin(origin); // Testing Initialize(const mitk::PixelType& type, const mitk::Geometry3D& geometry, unsigned int slices) with PlaneGeometry and GetData(): "; imgMem->Initialize( mitk::MakePixelType(), *planegeometry); MITK_TEST_CONDITION_REQUIRED( imgMem->GetGeometry()->GetOrigin() == static_cast(planegeometry)->GetOrigin(), "Testing correct setting of geometry via initialize!"); try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int*)imgMemAcc.GetData(); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED( p!=NULL, "GetData() returned valid pointer."); // Testing Initialize(const mitk::PixelType& type, int sDim, const mitk::PlaneGeometry& geometry) and GetData(): "; imgMem->Initialize( mitk::MakePixelType() , 40, *planegeometry); try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int*)imgMemAcc.GetData(); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED( p != NULL, "GetData() returned valid pointer."); //----------------- // testing origin information and methods MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetGeometry()->GetOrigin(), origin), "Testing correctness of origin via GetGeometry()->GetOrigin(): "); // Setting origin via SetOrigin(origin): "; mitk::FillVector3D(origin, 37.0, 17.92, 27.83); imgMem->SetOrigin(origin); // Test origin MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetGeometry()->GetOrigin(), origin), "Testing correctness of changed origin via GetGeometry()->GetOrigin(): "); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetSlicedGeometry()->GetPlaneGeometry(0)->GetOrigin(), origin), "Testing correctness of changed origin via GetSlicedGeometry()->GetPlaneGeometry(0)->GetOrigin(): "); //----------------- // testing spacing information and methodsunsigned int dim[]={100,100,20}; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(imgMem->GetGeometry()->GetSpacing(), spacing), "Testing correct spacing from Geometry3D!"); mitk::FillVector3D(spacing, 7.0, 0.92, 1.83); imgMem->SetSpacing(spacing); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetGeometry()->GetSpacing(), spacing), "Testing correctness of changed spacing via GetGeometry()->GetSpacing(): "); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetSlicedGeometry()->GetPlaneGeometry(0)->GetSpacing(), spacing), "Testing correctness of changed spacing via GetSlicedGeometry()->GetPlaneGeometry(0)->GetSpacing(): "); mitk::Image::Pointer vecImg = mitk::Image::New(); try { mitk::ImageReadAccessor imgMemAcc(imgMem); vecImg->Initialize( imgMem->GetPixelType(), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/ ); vecImg->SetImportChannel(const_cast(imgMemAcc.GetData()), 0, mitk::Image::CopyMemory ); vecImg->SetImportChannel(const_cast(imgMemAcc.GetData()), 1, mitk::Image::CopyMemory ); mitk::ImageReadAccessor vecImgAcc(vecImg); mitk::ImageReadAccessor vecImgAcc0(vecImg, vecImg->GetChannelData(0)); mitk::ImageReadAccessor vecImgAcc1(vecImg, vecImg->GetChannelData(1)); MITK_TEST_CONDITION_REQUIRED(vecImgAcc0.GetData() != NULL && vecImgAcc1.GetData() != NULL, "Testing set and return of channel data!"); MITK_TEST_CONDITION_REQUIRED( vecImg->IsValidSlice(0,0,1) , ""); MITK_TEST_OUTPUT(<< " Testing whether CopyMemory worked"); MITK_TEST_CONDITION_REQUIRED(imgMemAcc.GetData() != vecImgAcc.GetData(), ""); MITK_TEST_OUTPUT(<< " Testing destruction after SetImportChannel"); vecImg = NULL; MITK_TEST_CONDITION_REQUIRED(vecImg.IsNull() , "testing destruction!"); } catch (mitk::Exception& e) { MITK_ERROR << e.what(); } //----------------- MITK_TEST_OUTPUT(<< "Testing initialization via vtkImageData"); MITK_TEST_OUTPUT(<< " Setting up vtkImageData"); vtkImageData* vtkimage = vtkImageData::New(); vtkimage->Initialize(); vtkimage->SetDimensions( 2, 3, 4); double vtkorigin[] = {-350,-358.203, -1363.5}; vtkimage->SetOrigin(vtkorigin); mitk::Point3D vtkoriginAsMitkPoint; mitk::vtk2itk(vtkorigin, vtkoriginAsMitkPoint); double vtkspacing[] = {1.367, 1.367, 2}; vtkimage->SetSpacing(vtkspacing); vtkimage->AllocateScalars(VTK_SHORT,1); std::cout<<"[PASSED]"<Initialize(vtkimage); MITK_TEST_CONDITION_REQUIRED(mitkByVtkImage->IsInitialized(), ""); vtkimage->Delete(); MITK_TEST_OUTPUT(<< " Testing whether spacing has been correctly initialized from vtkImageData"); mitk::Vector3D spacing2 = mitkByVtkImage->GetGeometry()->GetSpacing(); mitk::Vector3D vtkspacingAsMitkVector; mitk::vtk2itk(vtkspacing, vtkspacingAsMitkVector); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(spacing2,vtkspacingAsMitkVector), ""); MITK_TEST_OUTPUT(<< " Testing whether GetSlicedGeometry(0)->GetOrigin() has been correctly initialized from vtkImageData"); mitk::Point3D origin2 = mitkByVtkImage->GetSlicedGeometry(0)->GetOrigin(); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(origin2,vtkoriginAsMitkPoint), ""); MITK_TEST_OUTPUT(<< " Testing whether GetGeometry()->GetOrigin() has been correctly initialized from vtkImageData"); origin2 = mitkByVtkImage->GetGeometry()->GetOrigin(); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(origin2,vtkoriginAsMitkPoint), ""); // TODO test the following initializers on channel-incorporation // void mitk::Image::Initialize(const mitk::PixelType& type, unsigned int dimension, unsigned int *dimensions, unsigned int channels) // void mitk::Image::Initialize(const mitk::PixelType& type, int sDim, const mitk::PlaneGeometry& geometry2d, bool flipped, unsigned int channels, int tDim ) // void mitk::Image::Initialize(const mitk::Image* image) // void mitk::Image::Initialize(const mitkIpPicDescriptor* pic, int channels, int tDim, int sDim) //mitk::Image::Pointer vecImg = mitk::Image::New(); //vecImg->Initialize(PixelType(typeid(float), 6, itk::ImageIOBase::SYMMETRICSECONDRANKTENSOR), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/, false /*shiftBoundingBoxMinimumToZero*/ ); //vecImg->Initialize(PixelType(typeid(itk::Vector)), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/, false /*shiftBoundingBoxMinimumToZero*/ ); // testing access by index coordinates and by world coordinates MITK_TEST_CONDITION_REQUIRED(argc == 2, "Check if test image is accessible!"); const std::string filename = std::string(argv[1]); mitk::ItkImageFileReader::Pointer imageReader = mitk::ItkImageFileReader::New(); try { imageReader->SetFileName(filename); imageReader->Update(); } catch(...) { MITK_TEST_FAILED_MSG(<< "Could not read file for testing: " << filename); return 0; } mitk::Image::Pointer image = imageReader->GetOutput(); mitk::Point3D point; mitk::ScalarType value = -1.; mitkPixelTypeMultiplex3(TestRandomPixelAccess,image->GetImageDescriptor()->GetChannelTypeById(0),image,point,value) { // testing the clone method of mitk::Image mitk::Image::Pointer cloneImage = image->Clone(); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetDimension() == image->GetDimension(), "Clone (testing dimension)"); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetPixelType() == image->GetPixelType(), "Clone (testing pixel type)"); // After cloning an image the geometry of both images should be equal too MITK_TEST_CONDITION_REQUIRED(cloneImage->GetGeometry()->GetOrigin() == image->GetGeometry()->GetOrigin(), "Clone (testing origin)"); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetGeometry()->GetSpacing() == image->GetGeometry()->GetSpacing(), "Clone (testing spacing)"); MITK_TEST_CONDITION_REQUIRED(mitk::MatrixEqualElementWise(cloneImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix(), image->GetGeometry()->GetIndexToWorldTransform()->GetMatrix()), "Clone (testing transformation matrix)"); MITK_TEST_CONDITION_REQUIRED(mitk::MatrixEqualElementWise(cloneImage->GetTimeGeometry()->GetGeometryForTimeStep(cloneImage->GetDimension(3)-1)->GetIndexToWorldTransform()->GetMatrix(), cloneImage->GetTimeGeometry()->GetGeometryForTimeStep(image->GetDimension(3)-1)->GetIndexToWorldTransform()->GetMatrix()), "Clone(testing time sliced geometry)"); for (unsigned int i = 0u; i < cloneImage->GetDimension(); ++i) { MITK_TEST_CONDITION_REQUIRED(cloneImage->GetDimension(i) == image->GetDimension(i), "Clone (testing dimension " << i << ")"); } } //access via itk if(image->GetDimension()> 3) // CastToItk only works with 3d images so we need to check for 4d images { mitk::ImageTimeSelector::Pointer selector = mitk::ImageTimeSelector::New(); selector->SetTimeNr(0); selector->SetInput(image); selector->Update(); image = selector->GetOutput(); } if(image->GetDimension()==3) { typedef itk::Image ItkFloatImage3D; ItkFloatImage3D::Pointer itkimage; try { mitk::CastToItkImage(image, itkimage); MITK_TEST_CONDITION_REQUIRED(itkimage.IsNotNull(), "Test conversion to itk::Image!"); } catch (std::exception& e) { MITK_INFO << e.what(); } mitk::Point3D itkPhysicalPoint; image->GetGeometry()->WorldToItkPhysicalPoint(point, itkPhysicalPoint); MITK_INFO << "ITKPoint " << itkPhysicalPoint[0] << " "<< itkPhysicalPoint[1] << " "<< itkPhysicalPoint[2] << ""; mitk::Point3D backTransformedPoint; image->GetGeometry()->ItkPhysicalPointToWorld(itkPhysicalPoint, backTransformedPoint); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(point,backTransformedPoint), "Testing world->itk-physical->world consistency"); itk::Index<3> idx; bool status = itkimage->TransformPhysicalPointToIndex(itkPhysicalPoint, idx); MITK_INFO << "ITK Index " << idx[0] << " "<< idx[1] << " "<< idx[2] << ""; if(status && value != -1.) { float valByItk = itkimage->GetPixel(idx); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(valByItk, value), "Compare value of pixel returned by mitk in comparison to itk"); } else { MITK_WARN<< "Index is out buffered region!"; } } else { MITK_INFO << "Image does not contain three dimensions, some test cases are skipped!"; } // clone generated 3D image with one slice in z direction (cf. bug 11058) unsigned int* threeDdim = new unsigned int[3]; threeDdim[0] = 100; threeDdim[1] = 200; threeDdim[2] = 1; mitk::Image::Pointer threeDImage = mitk::Image::New(); threeDImage->Initialize(mitk::MakeScalarPixelType(), 3, threeDdim); mitk::Image::Pointer cloneThreeDImage = threeDImage->Clone(); // check that the clone image has the same dimensionality as the source image MITK_TEST_CONDITION_REQUIRED( cloneThreeDImage->GetDimension() == 3, "Testing if the clone image initializes with 3D!"); MITK_TEST_CONDITION_REQUIRED( ImageVtkDataReferenceCheck(argv[1]), "Checking reference count of Image after using GetVtkImageData()"); MITK_TEST_END(); }