diff --git a/Core/Code/DataManagement/mitkBaseGeometry.cpp b/Core/Code/DataManagement/mitkBaseGeometry.cpp index 4d7d129416..b4c5572fdf 100644 --- a/Core/Code/DataManagement/mitkBaseGeometry.cpp +++ b/Core/Code/DataManagement/mitkBaseGeometry.cpp @@ -1,900 +1,904 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include "mitkBaseGeometry.h" #include "mitkvector.h" #include "mitkMatrixConvert.h" #include #include #include "mitkRotationOperation.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkApplyTransformMatrixOperation.h" #include "mitkPointOperation.h" #include "mitkInteractionConst.h" mitk::BaseGeometry::BaseGeometry(): Superclass(), mitk::OperationActor(), m_Valid(true), m_FrameOfReferenceID(0), m_IndexToWorldTransformLastModified(0) { FillVector3D(m_FloatSpacing, 1,1,1); m_VtkMatrix = vtkMatrix4x4::New(); m_VtkIndexToWorldTransform = vtkMatrixToLinearTransform::New(); m_VtkIndexToWorldTransform->SetInput(m_VtkMatrix); Initialize(); } mitk::BaseGeometry::BaseGeometry(const BaseGeometry& other): Superclass(), m_TimeBounds(other.m_TimeBounds), m_Valid(other.m_Valid), m_FrameOfReferenceID(other.m_FrameOfReferenceID), m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified), m_Origin(other.m_Origin) { // DEPRECATED(m_RotationQuaternion = other.m_RotationQuaternion); // AffineGeometryFrame SetBounds(other.GetBounds()); m_VtkMatrix = vtkMatrix4x4::New(); m_VtkMatrix->DeepCopy(other.m_VtkMatrix); FillVector3D(m_FloatSpacing,other.m_FloatSpacing[0],other.m_FloatSpacing[1],other.m_FloatSpacing[2]); m_VtkIndexToWorldTransform = vtkMatrixToLinearTransform::New(); m_VtkIndexToWorldTransform->DeepCopy(other.m_VtkIndexToWorldTransform); m_VtkIndexToWorldTransform->SetInput(m_VtkMatrix); other.InitializeGeometry(this); } mitk::BaseGeometry::~BaseGeometry() { m_VtkMatrix->Delete(); m_VtkIndexToWorldTransform->Delete(); } const mitk::Point3D& mitk::BaseGeometry::GetOrigin() const { return m_Origin; } void mitk::BaseGeometry::SetOrigin(const Point3D & origin) { if(origin!=GetOrigin()) { m_Origin = origin; m_IndexToWorldTransform->SetOffset(m_Origin.GetVectorFromOrigin()); Modified(); TransferItkToVtkTransform(); } } void mitk::BaseGeometry::TransferItkToVtkTransform() { // copy m_IndexToWorldTransform into m_VtkIndexToWorldTransform TransferItkTransformToVtkMatrix(m_IndexToWorldTransform.GetPointer(), m_VtkMatrix); m_VtkIndexToWorldTransform->Modified(); } void mitk::BaseGeometry::CopySpacingFromTransform(mitk::AffineTransform3D* transform, mitk::Vector3D& spacing, float floatSpacing[3]) { mitk::AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = transform->GetMatrix().GetVnlMatrix(); spacing[0]=vnlmatrix.get_column(0).magnitude(); spacing[1]=vnlmatrix.get_column(1).magnitude(); spacing[2]=vnlmatrix.get_column(2).magnitude(); floatSpacing[0]=spacing[0]; floatSpacing[1]=spacing[1]; floatSpacing[2]=spacing[2]; } void mitk::BaseGeometry::Initialize() { float b[6] = {0,1,0,1,0,1}; SetFloatBounds(b); if(m_IndexToWorldTransform.IsNull()) m_IndexToWorldTransform = TransformType::New(); else m_IndexToWorldTransform->SetIdentity(); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); m_VtkMatrix->Identity(); m_TimeBounds[0]=ScalarTypeNumericTraits::NonpositiveMin(); m_TimeBounds[1]=ScalarTypeNumericTraits::max(); m_FrameOfReferenceID = 0; this->InternPostInitialize(); } void mitk::BaseGeometry::SetFloatBounds(const float bounds[6]) { mitk::BoundingBox::BoundsArrayType b; const float *input = bounds; int i=0; for(mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6 ;++i) *it++ = (mitk::ScalarType)*input++; SetBounds(b); } void mitk::BaseGeometry::SetFloatBounds(const double bounds[6]) { mitk::BoundingBox::BoundsArrayType b; const double *input = bounds; int i=0; for(mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6 ;++i) *it++ = (mitk::ScalarType)*input++; SetBounds(b); } /** Initialize the geometry */ void mitk::BaseGeometry::InitializeGeometry(BaseGeometry* newGeometry) const { newGeometry->SetBounds(m_BoundingBox->GetBounds()); // we have to create a new transform!! newGeometry->SetTimeBounds(m_TimeBounds); newGeometry->SetFrameOfReferenceID(GetFrameOfReferenceID()); if(m_IndexToWorldTransform) { TransformType::Pointer indexToWorldTransform = TransformType::New(); indexToWorldTransform->SetCenter( m_IndexToWorldTransform->GetCenter() ); indexToWorldTransform->SetMatrix( m_IndexToWorldTransform->GetMatrix() ); indexToWorldTransform->SetOffset( m_IndexToWorldTransform->GetOffset() ); newGeometry->SetIndexToWorldTransform(indexToWorldTransform); } this->InternPostInitializeGeometry(newGeometry); } /** Set the bounds */ void mitk::BaseGeometry::SetBounds(const BoundsArrayType& bounds) { + InternPreSetBounds(bounds); + m_BoundingBox = BoundingBoxType::New(); BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New(); BoundingBoxType::PointType p; BoundingBoxType::PointIdentifier pointid; for(pointid=0; pointid<2;++pointid) { unsigned int i; for(i=0; iInsertElement(pointid, p); } m_BoundingBox->SetPoints(pointscontainer); m_BoundingBox->ComputeBoundingBox(); this->Modified(); } +void mitk::BaseGeometry::InternPreSetBounds(const BoundsArrayType& bounds){}; + void mitk::BaseGeometry::SetIndexToWorldTransform(mitk::AffineTransform3D* transform) { if(m_IndexToWorldTransform.GetPointer() != transform) { m_IndexToWorldTransform = transform; CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); TransferItkToVtkTransform(); Modified(); } } const mitk::BaseGeometry::BoundsArrayType mitk::BaseGeometry::GetBounds() const { assert(m_BoundingBox.IsNotNull()); return m_BoundingBox->GetBounds(); } bool mitk::BaseGeometry::IsValid() const { bool isValid = m_Valid; isValid = isValid && this->InternPostIsValid(); return isValid; } bool mitk::BaseGeometry::InternPostIsValid() const { return true; } const float* mitk::BaseGeometry::GetFloatSpacing() const { return m_FloatSpacing; } bool mitk::Equal(const BaseGeometry::TransformType *leftHandSide, const BaseGeometry::TransformType *rightHandSide, ScalarType eps, bool verbose ) { //Compare IndexToWorldTransform Matrix if( !mitk::MatrixEqualElementWise( leftHandSide->GetMatrix(), rightHandSide->GetMatrix() ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry::TransformType )] Index to World Transformation matrix differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetMatrix() << " : leftHandSide is " << leftHandSide->GetMatrix() << " and tolerance is " << eps; } return false; } return true; } bool mitk::Equal( const mitk::BaseGeometry::BoundingBoxType *leftHandSide, const mitk::BaseGeometry::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) { bool result = true; if( rightHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry::BoundingBoxType )] rightHandSide NULL."; return false; } if( leftHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry::BoundingBoxType )] leftHandSide NULL."; return false; } BaseGeometry::BoundsArrayType rightBounds = rightHandSide->GetBounds(); BaseGeometry::BoundsArrayType leftBounds = leftHandSide->GetBounds(); BaseGeometry::BoundsArrayType::Iterator itLeft = leftBounds.Begin(); for( BaseGeometry::BoundsArrayType::Iterator itRight = rightBounds.Begin(); itRight != rightBounds.End(); ++itRight) { if(( !mitk::Equal( *itLeft, *itRight, eps )) ) { if(verbose) { MITK_INFO << "[( BaseGeometry::BoundingBoxType )] bounds are not equal."; MITK_INFO << "rightHandSide is " << setprecision(12) << *itRight << " : leftHandSide is " << *itLeft << " and tolerance is " << eps; } result = false; } itLeft++; } return result; } bool mitk::Equal(const mitk::BaseGeometry *leftHandSide, const mitk::BaseGeometry *rightHandSide, ScalarType eps, bool verbose) { bool result = true; if( rightHandSide == NULL ) { if(verbose) MITK_INFO << "[( BaseGeometry )] rightHandSide NULL."; return false; } if( leftHandSide == NULL) { if(verbose) MITK_INFO << "[( BaseGeometry )] leftHandSide NULL."; return false; } //Compare spacings if( !mitk::Equal( leftHandSide->GetSpacing(), rightHandSide->GetSpacing(), eps ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Spacing differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetSpacing() << " : leftHandSide is " << leftHandSide->GetSpacing() << " and tolerance is " << eps; } result = false; } //Compare Origins if( !mitk::Equal( leftHandSide->GetOrigin(), rightHandSide->GetOrigin(), eps ) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Origin differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetOrigin() << " : leftHandSide is " << leftHandSide->GetOrigin() << " and tolerance is " << eps; } result = false; } //Compare Axis and Extents for( unsigned int i=0; i<3; ++i) { if( !mitk::Equal( leftHandSide->GetAxisVector(i), rightHandSide->GetAxisVector(i), eps)) { if(verbose) { MITK_INFO << "[( BaseGeometry )] AxisVector #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetAxisVector(i) << " : leftHandSide is " << leftHandSide->GetAxisVector(i) << " and tolerance is " << eps; } result = false; } if( !mitk::Equal( leftHandSide->GetExtent(i), rightHandSide->GetExtent(i), eps) ) { if(verbose) { MITK_INFO << "[( BaseGeometry )] Extent #" << i << " differ"; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide->GetExtent(i) << " : leftHandSide is " << leftHandSide->GetExtent(i) << " and tolerance is " << eps; } result = false; } } //Compare BoundingBoxes if( !mitk::Equal( leftHandSide->GetBoundingBox(), rightHandSide->GetBoundingBox(), eps, verbose) ) { result = false; } //Compare IndexToWorldTransform Matrix if( !mitk::Equal( leftHandSide->GetIndexToWorldTransform(), rightHandSide->GetIndexToWorldTransform(), eps, verbose) ) { result = false; } return result; } void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D& aSpacing) { if(mitk::Equal(m_Spacing, aSpacing) == false) { assert(aSpacing[0]>0 && aSpacing[1]>0 && aSpacing[2]>0); m_Spacing = aSpacing; AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = m_IndexToWorldTransform->GetMatrix().GetVnlMatrix(); mitk::VnlVector col; col = vnlmatrix.get_column(0); col.normalize(); col*=aSpacing[0]; vnlmatrix.set_column(0, col); col = vnlmatrix.get_column(1); col.normalize(); col*=aSpacing[1]; vnlmatrix.set_column(1, col); col = vnlmatrix.get_column(2); col.normalize(); col*=aSpacing[2]; vnlmatrix.set_column(2, col); Matrix3D matrix; matrix = vnlmatrix; AffineTransform3D::Pointer transform = AffineTransform3D::New(); transform->SetMatrix(matrix); transform->SetOffset(m_IndexToWorldTransform->GetOffset()); SetIndexToWorldTransform(transform.GetPointer()); itk2vtk(m_Spacing, m_FloatSpacing); } } mitk::Vector3D mitk::BaseGeometry::GetAxisVector(unsigned int direction) const { Vector3D frontToBack; frontToBack.SetVnlVector(m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction)); frontToBack *= GetExtent(direction); return frontToBack; } mitk::ScalarType mitk::BaseGeometry::GetExtent(unsigned int direction) const { assert(m_BoundingBox.IsNotNull()); if (direction>=NDimensions) mitkThrow() << "Direction is too big. This geometry is for 3D Data"; BoundsArrayType bounds = m_BoundingBox->GetBounds(); return bounds[direction*2+1]-bounds[direction*2]; } bool mitk::BaseGeometry::Is2DConvertable() { bool isConvertableWithoutLoss = true; do { if (this->GetSpacing()[2] != 1) { isConvertableWithoutLoss = false; break; } if (this->GetOrigin()[2] != 0) { isConvertableWithoutLoss = false; break; } mitk::Vector3D col0, col1, col2; col0.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0)); col1.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1)); col2.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2)); if ((col0[2] != 0) || (col1[2] != 0) || (col2[0] != 0) || (col2[1] != 0) || (col2[2] != 1)) { isConvertableWithoutLoss = false; break; } } while (0); return isConvertableWithoutLoss; } mitk::Point3D mitk::BaseGeometry::GetCenter() const { assert(m_BoundingBox.IsNotNull()); return m_IndexToWorldTransform->TransformPoint(m_BoundingBox->GetCenter()); } double mitk::BaseGeometry::GetDiagonalLength2() const { Vector3D diagonalvector = GetCornerPoint()-GetCornerPoint(false, false, false); return diagonalvector.GetSquaredNorm(); } //##Documentation //## @brief Get the length of the diagonal of the bounding-box in mm //## double mitk::BaseGeometry::GetDiagonalLength() const { return sqrt(GetDiagonalLength2()); } mitk::Point3D mitk::BaseGeometry::GetCornerPoint(int id) const { assert(id >= 0); assert(m_BoundingBox.IsNotNull()); BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; switch(id) { case 0: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[4]); break; case 1: FillVector3D(cornerpoint, bounds[0],bounds[2],bounds[5]); break; case 2: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[4]); break; case 3: FillVector3D(cornerpoint, bounds[0],bounds[3],bounds[5]); break; case 4: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[4]); break; case 5: FillVector3D(cornerpoint, bounds[1],bounds[2],bounds[5]); break; case 6: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[4]); break; case 7: FillVector3D(cornerpoint, bounds[1],bounds[3],bounds[5]); break; default: { itkExceptionMacro(<<"A cube only has 8 corners. These are labeled 0-7."); } } return m_IndexToWorldTransform->TransformPoint(cornerpoint); } mitk::Point3D mitk::BaseGeometry::GetCornerPoint(bool xFront, bool yFront, bool zFront) const { assert(m_BoundingBox.IsNotNull()); BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; cornerpoint[0] = (xFront ? bounds[0] : bounds[1]); cornerpoint[1] = (yFront ? bounds[2] : bounds[3]); cornerpoint[2] = (zFront ? bounds[4] : bounds[5]); return m_IndexToWorldTransform->TransformPoint(cornerpoint); } mitk::ScalarType mitk::BaseGeometry::GetExtentInMM(int direction) const { return m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction).magnitude()*GetExtent(direction); } void mitk::BaseGeometry::SetExtentInMM(int direction, ScalarType extentInMM) { ScalarType len = GetExtentInMM(direction); if(fabs(len - extentInMM)>=mitk::eps) { AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = m_IndexToWorldTransform->GetMatrix().GetVnlMatrix(); if(len>extentInMM) vnlmatrix.set_column(direction, vnlmatrix.get_column(direction)/len*extentInMM); else vnlmatrix.set_column(direction, vnlmatrix.get_column(direction)*extentInMM/len); Matrix3D matrix; matrix = vnlmatrix; m_IndexToWorldTransform->SetMatrix(matrix); Modified(); } } bool mitk::BaseGeometry::IsInside(const mitk::Point3D& p) const { mitk::Point3D index; WorldToIndex(p, index); return IsIndexInside(index); } bool mitk::BaseGeometry::IsIndexInside(const mitk::Point3D& index) const { bool inside = false; //if it is an image geometry, we need to convert the index to discrete values //this is done by applying the rounding function also used in WorldToIndex (see line 323) inside = m_BoundingBox->IsInside(index); return inside; } //##Documentation //## @brief Convenience method for working with ITK indices template bool mitk::BaseGeometry::IsIndexInside(const itk::Index &index) const { int i, dim=index.GetIndexDimension(); Point3D pt_index; pt_index.Fill(0); for ( i = 0; i < dim; ++i ) { pt_index[i] = index[i]; } return IsIndexInside(pt_index); } void mitk::BaseGeometry::WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const { BackTransform(pt_mm, pt_units); } void mitk::BaseGeometry::WorldToIndex( const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const { BackTransform( vec_mm, vec_units); } void mitk::BaseGeometry::BackTransform(const mitk::Vector3D& in, mitk::Vector3D& out) const { // Get WorldToIndex transform if (m_IndexToWorldTransformLastModified != m_IndexToWorldTransform->GetMTime()) { m_InvertedTransform = TransformType::New(); if (!m_IndexToWorldTransform->GetInverse( m_InvertedTransform.GetPointer() )) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed." ); } m_IndexToWorldTransformLastModified = m_IndexToWorldTransform->GetMTime(); } // Check for valid matrix inversion const TransformType::MatrixType& inverse = m_InvertedTransform->GetMatrix(); if(inverse.GetVnlMatrix().has_nans()) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed. Matrix was: " << std::endl << m_IndexToWorldTransform->GetMatrix() << "Suggested inverted matrix is:" << std::endl << inverse ); } // Transform vector for (unsigned int i = 0; i < 3; i++) { out[i] = 0.0; for (unsigned int j = 0; j < 3; j++) { out[i] += inverse[i][j]*in[j]; } } } void mitk::BaseGeometry::BackTransform(const mitk::Point3D &in, mitk::Point3D& out) const { ScalarType temp[3]; unsigned int i, j; const TransformType::OffsetType& offset = m_IndexToWorldTransform->GetOffset(); // Remove offset for (j = 0; j < 3; j++) { temp[j] = in[j] - offset[j]; } // Get WorldToIndex transform if (m_IndexToWorldTransformLastModified != m_IndexToWorldTransform->GetMTime()) { m_InvertedTransform = TransformType::New(); if (!m_IndexToWorldTransform->GetInverse( m_InvertedTransform.GetPointer() )) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed." ); } m_IndexToWorldTransformLastModified = m_IndexToWorldTransform->GetMTime(); } // Check for valid matrix inversion const TransformType::MatrixType& inverse = m_InvertedTransform->GetMatrix(); if(inverse.GetVnlMatrix().has_nans()) { itkExceptionMacro( "Internal ITK matrix inversion error, cannot proceed. Matrix was: " << std::endl << m_IndexToWorldTransform->GetMatrix() << "Suggested inverted matrix is:" << std::endl << inverse ); } // Transform point for (i = 0; i < 3; i++) { out[i] = 0.0; for (j = 0; j < 3; j++) { out[i] += inverse[i][j]*temp[j]; } } } mitk::VnlVector mitk::BaseGeometry::GetOriginVnl() const { return const_cast(this)->m_Origin.GetVnlVector(); } vtkLinearTransform* mitk::BaseGeometry::GetVtkTransform() const { return (vtkLinearTransform*)m_VtkIndexToWorldTransform; } void mitk::BaseGeometry::SetIdentity() { m_IndexToWorldTransform->SetIdentity(); m_Origin.Fill(0); Modified(); TransferItkToVtkTransform(); } void mitk::BaseGeometry::TransferVtkToItkTransform() { TransferVtkMatrixToItkTransform(m_VtkMatrix, m_IndexToWorldTransform.GetPointer()); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); } void mitk::BaseGeometry::Compose( const mitk::BaseGeometry::TransformType * other, bool pre ) { m_IndexToWorldTransform->Compose(other, pre); CopySpacingFromTransform(m_IndexToWorldTransform, m_Spacing, m_FloatSpacing); vtk2itk(m_IndexToWorldTransform->GetOffset(), m_Origin); Modified(); TransferItkToVtkTransform(); } void mitk::BaseGeometry::Compose( const vtkMatrix4x4 * vtkmatrix, bool pre ) { mitk::BaseGeometry::TransformType::Pointer itkTransform = mitk::BaseGeometry::TransformType::New(); TransferVtkMatrixToItkTransform(vtkmatrix, itkTransform.GetPointer()); Compose(itkTransform, pre); } void mitk::BaseGeometry::Translate(const Vector3D & vector) { if((vector[0] != 0) || (vector[1] != 0) || (vector[2] != 0)) { this->SetOrigin(m_Origin + vector); } } void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const { pt_mm = m_IndexToWorldTransform->TransformPoint(pt_units); } void mitk::BaseGeometry::IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const { vec_mm = m_IndexToWorldTransform->TransformVector(vec_units); } #include void mitk::BaseGeometry::ExecuteOperation(Operation* operation) { vtkTransform *vtktransform = vtkTransform::New(); vtktransform->SetMatrix(m_VtkMatrix); switch (operation->GetOperationType()) { case OpNOTHING: break; case OpMOVE: { mitk::PointOperation *pointOp = dynamic_cast(operation); if (pointOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } mitk::Point3D newPos = pointOp->GetPoint(); ScalarType data[3]; vtktransform->GetPosition(data); vtktransform->PostMultiply(); vtktransform->Translate(newPos[0], newPos[1], newPos[2]); vtktransform->PreMultiply(); break; } case OpSCALE: { mitk::PointOperation *pointOp = dynamic_cast(operation); if (pointOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } mitk::Point3D newScale = pointOp->GetPoint(); ScalarType data[3]; /* calculate new scale: newscale = oldscale * (oldscale + scaletoadd)/oldscale */ data[0] = 1 + (newScale[0] / GetMatrixColumn(0).magnitude()); data[1] = 1 + (newScale[1] / GetMatrixColumn(1).magnitude()); data[2] = 1 + (newScale[2] / GetMatrixColumn(2).magnitude()); mitk::Point3D center = const_cast(m_BoundingBox.GetPointer())->GetCenter(); ScalarType pos[3]; vtktransform->GetPosition(pos); vtktransform->PostMultiply(); vtktransform->Translate(-pos[0], -pos[1], -pos[2]); vtktransform->Translate(-center[0], -center[1], -center[2]); vtktransform->PreMultiply(); vtktransform->Scale(data[0], data[1], data[2]); vtktransform->PostMultiply(); vtktransform->Translate(+center[0], +center[1], +center[2]); vtktransform->Translate(pos[0], pos[1], pos[2]); vtktransform->PreMultiply(); break; } case OpROTATE: { mitk::RotationOperation *rotateOp = dynamic_cast(operation); if (rotateOp == NULL) { //mitk::StatusBar::GetInstance()->DisplayText("received wrong type of operation!See mitkAffineInteractor.cpp", 10000); return; } Vector3D rotationVector = rotateOp->GetVectorOfRotation(); Point3D center = rotateOp->GetCenterOfRotation(); ScalarType angle = rotateOp->GetAngleOfRotation(); vtktransform->PostMultiply(); vtktransform->Translate(-center[0], -center[1], -center[2]); vtktransform->RotateWXYZ(angle, rotationVector[0], rotationVector[1], rotationVector[2]); vtktransform->Translate(center[0], center[1], center[2]); vtktransform->PreMultiply(); break; } case OpRESTOREPLANEPOSITION: { //Copy necessary to avoid vtk warning vtkMatrix4x4* matrix = vtkMatrix4x4::New(); TransferItkTransformToVtkMatrix(dynamic_cast(operation)->GetTransform().GetPointer(), matrix); vtktransform->SetMatrix(matrix); break; } case OpAPPLYTRANSFORMMATRIX: { ApplyTransformMatrixOperation *applyMatrixOp = dynamic_cast< ApplyTransformMatrixOperation* >( operation ); vtktransform->SetMatrix(applyMatrixOp->GetMatrix()); break; } default: vtktransform->Delete(); return; } m_VtkMatrix->DeepCopy(vtktransform->GetMatrix()); TransferVtkToItkTransform(); Modified(); vtktransform->Delete(); } mitk::VnlVector mitk::BaseGeometry::GetMatrixColumn(unsigned int direction) const { return m_IndexToWorldTransform->GetMatrix().GetVnlMatrix().get_column(direction); } mitk::BoundingBox::Pointer mitk::BaseGeometry::CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D* transform) const { mitk::BoundingBox::PointsContainer::Pointer pointscontainer=mitk::BoundingBox::PointsContainer::New(); mitk::BoundingBox::PointIdentifier pointid=0; unsigned char i; if(transform!=NULL) { mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New(); transform->GetInverse(inverse); for(i=0; i<8; ++i) pointscontainer->InsertElement( pointid++, inverse->TransformPoint( GetCornerPoint(i) )); } else { for(i=0; i<8; ++i) pointscontainer->InsertElement( pointid++, GetCornerPoint(i) ); } mitk::BoundingBox::Pointer result = mitk::BoundingBox::New(); result->SetPoints(pointscontainer); result->ComputeBoundingBox(); return result; } void mitk::BaseGeometry::SetTimeBounds(const TimeBounds& timebounds) { if(m_TimeBounds != timebounds) { m_TimeBounds = timebounds; Modified(); } } 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); } diff --git a/Core/Code/DataManagement/mitkBaseGeometry.h b/Core/Code/DataManagement/mitkBaseGeometry.h index 9c57919293..a30bdeef2d 100644 --- a/Core/Code/DataManagement/mitkBaseGeometry.h +++ b/Core/Code/DataManagement/mitkBaseGeometry.h @@ -1,571 +1,575 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef BaseGeometry_H_HEADER_INCLUDED #define BaseGeometry_H_HEADER_INCLUDED #include #include #include "mitkoperationactor.h" #include #include "mitkvector.h" #include #include #include "itkScalableAffineTransform.h" #include class vtkMatrix4x4; class vtkMatrixToLinearTransform; class vtkLinearTransform; namespace mitk { //##Documentation //## @brief Standard 3D-BoundingBox typedef //## //## Standard 3D-BoundingBox typedef to get rid of template arguments (3D, type). typedef itk::BoundingBox BoundingBox; //##Documentation //## @brief Standard typedef for time-bounds typedef itk::FixedArray TimeBounds; typedef itk::FixedArray FixedArrayType; //##Documentation //## @brief BaseGeometry xxxxxxxxxxxxxx //## //## xxxxxxxxxxx //## //## Rule: everything is in mm (ms) if not stated otherwise. //## @ingroup Geometry class MITK_CORE_EXPORT BaseGeometry : public itk::Object, public OperationActor { public: mitkClassMacro(BaseGeometry, itk::Object); + //void testXYZ(){}; //xxx + // ********************************** TypeDef ********************************** typedef itk::QuaternionRigidTransform< ScalarType > QuaternionTransformType; typedef QuaternionTransformType::VnlQuaternionType VnlQuaternionType; typedef itk::ScalableAffineTransform TransformType; typedef itk::BoundingBox BoundingBoxType; typedef BoundingBoxType::BoundsArrayType BoundsArrayType; typedef BoundingBoxType::Pointer BoundingBoxPointer; // ********************************** Origin, Spacing ********************************** //##Documentation //## @brief Get the origin, e.g. the upper-left corner of the plane const Point3D& GetOrigin() const; //##Documentation //## @brief Set the origin, i.e. the upper-left corner of the plane //## void SetOrigin(const Point3D& origin); //##Documentation //## @brief Get the spacing (size of a pixel). //## itkGetConstReferenceMacro(Spacing, mitk::Vector3D); //##Documentation //## @brief Get the spacing as a float[3] array. const float* GetFloatSpacing() const; //##Documentation //## @brief Set the spacing (m_Spacing) virtual void SetSpacing(const mitk::Vector3D& aSpacing); //##Documentation //## @brief Get the origin as VnlVector //## //## \sa GetOrigin VnlVector GetOriginVnl() const; // ********************************** other functions ********************************** //##Documentation //## @brief Get the DICOM FrameOfReferenceID referring to the //## used world coordinate system itkGetConstMacro(FrameOfReferenceID, unsigned int); //##Documentation //## @brief Set the DICOM FrameOfReferenceID referring to the //## used world coordinate system itkSetMacro(FrameOfReferenceID, unsigned int); //##Documentation //## @brief Is this Geometry3D in a state that is valid? bool IsValid() const; // ********************************** Initialize ********************************** //##Documentation //## @brief Initialize the Geometry3D void Initialize(); void InitializeGeometry(Self * newGeometry) const; static void CopySpacingFromTransform(mitk::AffineTransform3D* transform, mitk::Vector3D& spacing, float floatSpacing[3]); // ********************************** Transformations Set/Get ********************************** // a bit of a misuse, but we want only doxygen to see the following: #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get the transformation used to convert from index //## to world coordinates itkGetObjectMacro(IndexToWorldTransform, AffineTransform3D); #endif //## @brief Set the transformation used to convert from index //## to world coordinates virtual void SetIndexToWorldTransform(mitk::AffineTransform3D* transform); //##Documentation //## @brief Convenience method for setting the ITK transform //## (m_IndexToWorldTransform) via an vtkMatrix4x4 //## \sa SetIndexToWorldTransform virtual void SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4* vtkmatrix); /** Set/Get the IndexToWorldTransform */ itkGetConstObjectMacro(IndexToWorldTransform, AffineTransform3D); itkGetObjectMacro(IndexToWorldTransform, AffineTransform3D); //##Documentation //## @brief Get the m_IndexToWorldTransform as a vtkLinearTransform vtkLinearTransform* GetVtkTransform() const; //##Documentation //## @brief Set the transform to identity and origin to 0 //## virtual void SetIdentity(); // ********************************** Transformations ********************************** //##Documentation //## @brief Copy the ITK transform //## (m_IndexToWorldTransform) to the VTK transform //## \sa SetIndexToWorldTransform void TransferItkToVtkTransform(); //##Documentation //## @brief Copy the VTK transform //## to the ITK transform (m_IndexToWorldTransform) //## \sa SetIndexToWorldTransform void TransferVtkToItkTransform(); //##Documentation //## @brief Compose new IndexToWorldTransform with a given transform. //## //## This method composes m_IndexToWorldTransform with another transform, //## modifying self to be the composition of self and other. //## If the argument pre is true, then other is precomposed with self; //## that is, the resulting transformation consists of first applying //## other to the source, followed by self. If pre is false or omitted, //## then other is post-composed with self; that is the resulting //## transformation consists of first applying self to the source, //## followed by other. void Compose( const BaseGeometry::TransformType * other, bool pre = 0 ); //##Documentation //## @brief Compose new IndexToWorldTransform with a given vtkMatrix4x4. //## //## Converts the vtkMatrix4x4 into a itk-transform and calls the previous method. void Compose( const vtkMatrix4x4 * vtkmatrix, bool pre = 0 ); //##Documentation //## @brief Translate the origin by a vector //## void Translate(const Vector3D& vector); //##Documentation //##@brief executes affine operations (translate, rotate, scale) virtual void ExecuteOperation(Operation* operation); //##Documentation //## @brief Convert world coordinates (in mm) of a \em point to (continuous!) index coordinates //## \warning If you need (discrete) integer index coordinates (e.g., for iterating easily over an image), //## use WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index). //## For further information about coordinates types, please see the Geometry documentation void WorldToIndex(const mitk::Point3D& pt_mm, mitk::Point3D& pt_units) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em vector //## \a vec_mm to (continuous!) index coordinates. //## For further information about coordinates types, please see the Geometry documentation void WorldToIndex(const mitk::Vector3D& vec_mm, mitk::Vector3D& vec_units) const; //##Documentation //## @brief Convert world coordinates (in mm) of a \em point to (discrete!) index coordinates. //## This method rounds to integer indices! //## For further information about coordinates types, please see the Geometry documentation template void WorldToIndex(const mitk::Point3D& pt_mm, itk::Index &index) const { typedef itk::Index IndexType; mitk::Point3D pt_units; this->WorldToIndex(pt_mm, pt_units); int i, dim=index.GetIndexDimension(); if(dim>3) { index.Fill(0); dim=3; } for(i=0;i( pt_units[i] ); } } //##Documentation //## @brief Convert (continuous or discrete) index coordinates of a \em vector //## \a vec_units to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation void IndexToWorld(const mitk::Vector3D& vec_units, mitk::Vector3D& vec_mm) const; //##Documentation //## @brief Convert (continuous or discrete) index coordinates of a \em point to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation void IndexToWorld(const mitk::Point3D& pt_units, mitk::Point3D& pt_mm) const; //##Documentation //## @brief Convert (discrete) index coordinates of a \em point to world coordinates (in mm) //## For further information about coordinates types, please see the Geometry documentation template void IndexToWorld(const itk::Index &index, mitk::Point3D& pt_mm ) const { mitk::Point3D pt_units; pt_units.Fill(0); int i, dim=index.GetIndexDimension(); if(dim>3) { dim=3; } for(i=0;i void ItkPhysicalPointToWorld(const itk::Point& itkPhysicalPoint, mitk::Point3D& pt_mm) const { mitk::vtk2itk(itkPhysicalPoint, pt_mm); } //##Documentation //## @brief Deprecated for use with ITK version 3.10 or newer. //## Convert world coordinates (in mm) of a \em point to //## ITK physical coordinates (in mm, but without a possible rotation) //## //## This method is useful if you have want to access an mitk::Image //## via an itk::Image. ITK v3.8 and older did not support rotated (tilted) //## images, i.e., ITK images are always parallel to the coordinate axes. //## When accessing a (possibly rotated) mitk::Image via an itk::Image //## the rotational part of the transformation in the Geometry3D is //## simply discarded; in other word: only the origin and spacing is //## used by ITK, not the complete matrix available in MITK. //## With WorldToItkPhysicalPoint you can convert an MITK world //## coordinate (including the rotation) into a coordinate that //## can be used with the ITK image as a ITK physical coordinate //## (excluding the rotation). template void WorldToItkPhysicalPoint(const mitk::Point3D& pt_mm, itk::Point& itkPhysicalPoint) const { mitk::vtk2itk(pt_mm, itkPhysicalPoint); } // ********************************** BoundingBox ********************************** /** Get the bounding box */ itkGetConstObjectMacro(BoundingBox, BoundingBoxType); //##Documentation //## @brief Get the time bounds (in ms) itkGetConstReferenceMacro(TimeBounds, TimeBounds); // a bit of a misuse, but we want only doxygen to see the following: #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get bounding box (in index/unit coordinates) itkGetConstObjectMacro(BoundingBox, BoundingBoxType); //##Documentation //## @brief Get bounding box (in index/unit coordinates) as a BoundsArrayType const BoundsArrayType GetBounds() const; #endif const BoundsArrayType GetBounds() const; //##Documentation //## \brief Set the bounding box (in index/unit coordinates) //## //## Only possible via the BoundsArray to make clear that a //## copy of the bounding-box is stored, not a reference to it. - virtual void SetBounds(const BoundsArrayType& bounds); + void SetBounds(const BoundsArrayType& bounds); //##Documentation //## @brief Set the bounding box (in index/unit coordinates) via a float array void SetFloatBounds(const float bounds[6]); //##Documentation //## @brief Set the bounding box (in index/unit coordinates) via a double array void SetFloatBounds(const double bounds[6]); //##Documentation //## @brief Get a VnlVector along bounding-box in the specified //## @a direction, length is spacing //## //## \sa GetAxisVector VnlVector GetMatrixColumn(unsigned int direction) const; //##Documentation //## @brief Calculates a bounding-box around the geometry relative //## to a coordinate system defined by a transform //## mitk::BoundingBox::Pointer CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D* transform) const; //##Documentation //## @brief Set the time bounds (in ms) virtual void SetTimeBounds(const TimeBounds& timebounds); // ********************************** Geometry ********************************** #ifdef DOXYGEN_SKIP //##Documentation //## @brief Get the extent of the bounding box (in index/unit coordinates) //## //## To access the extent in mm use GetExtentInMM ScalarType GetExtent(unsigned int direction) const; #endif /** Get the extent of the bounding box */ ScalarType GetExtent(unsigned int direction) const; //##Documentation //## @brief Get the extent of the bounding-box in the specified @a direction in mm //## //## Equals length of GetAxisVector(direction). ScalarType GetExtentInMM(int direction) const; //##Documentation //## @brief Get vector along bounding-box in the specified @a direction in mm //## //## The length of the vector is the size of the bounding-box in the //## specified @a direction in mm //## \sa GetMatrixColumn Vector3D GetAxisVector(unsigned int direction) const; //##Documentation //## @brief Checks, if the given geometry can be converted to 2D without information loss //## e.g. when a 2D image is saved, the matrix is usually cropped to 2x2, and when you load it back to MITK //## it will be filled with standard values. This function checks, if information would be lost during this //## procedure virtual bool Is2DConvertable(); //##Documentation //## @brief Get the center of the bounding-box in mm //## Point3D GetCenter() const; //##Documentation //## @brief Get the squared length of the diagonal of the bounding-box in mm //## double GetDiagonalLength2() const; //##Documentation //## @brief Get the length of the diagonal of the bounding-box in mm //## double GetDiagonalLength() const; //##Documentation //## @brief Get the position of the corner number \a id (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. virtual Point3D GetCornerPoint(int id) const; //##Documentation //## @brief Get the position of a corner (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. virtual Point3D GetCornerPoint(bool xFront=true, bool yFront=true, bool zFront=true) const; //##Documentation //## @brief Set the extent of the bounding-box in the specified @a direction in mm //## //## @note This changes the matrix in the transform, @a not the bounds, which are given in units! virtual void SetExtentInMM(int direction, ScalarType extentInMM); //##Documentation //## @brief Test whether the point \a p (world coordinates in mm) is //## inside the bounding box virtual bool IsInside(const mitk::Point3D& p) const; //##Documentation //## @brief Test whether the point \a p ((continous!)index coordinates in units) is //## inside the bounding box virtual bool IsIndexInside(const mitk::Point3D& index) const; //##Documentation //## @brief Convenience method for working with ITK indices template bool IsIndexInside(const itk::Index &index) const; protected: // ********************************** Constructor ********************************** BaseGeometry(); BaseGeometry(const BaseGeometry& other); virtual ~BaseGeometry(); itkGetConstMacro(IndexToWorldTransformLastModified, unsigned long); void BackTransform(const mitk::Point3D& in, mitk::Point3D& out) const; //Without redundant parameter Point3D void BackTransform(const mitk::Vector3D& in, mitk::Vector3D& out) const; //##Documentation //## @brief Deprecated void BackTransform(const mitk::Point3D& at, const mitk::Vector3D& in, mitk::Vector3D& out) const; static const std::string GetTransformAsString( TransformType* transformType ); //Internal Functions virtual bool InternPostIsValid() const; virtual void InternPostInitialize() {}; virtual void InternPostInitializeGeometry(Self * newGeometry) const{}; + virtual void InternPreSetBounds(const BoundsArrayType& bounds); + // ********************************** Variables ********************************** AffineTransform3D::Pointer m_IndexToWorldTransform; vtkMatrixToLinearTransform* m_VtkIndexToWorldTransform; vtkMatrix4x4* m_VtkMatrix; bool m_Valid; unsigned int m_FrameOfReferenceID; mutable mitk::TimeBounds m_TimeBounds; mutable BoundingBoxPointer m_BoundingBox; //##Documentation //## @brief Origin, i.e. upper-left corner of the plane //## Point3D m_Origin; //##Documentation //## @brief Spacing of the data. Only significant if the geometry describes //## an Image (m_ImageGeometry==true). mitk::Vector3D m_Spacing; static const unsigned int NDimensions = 3; mutable TransformType::Pointer m_InvertedTransform; //this was private mutable unsigned long m_IndexToWorldTransformLastModified; //this was private float m_FloatSpacing[3]; //this was private // DEPRECATED(VnlQuaternionType m_RotationQuaternion); //this was private }; // ********************************** Equal Functions ********************************** // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two bounding boxes (BoundingBoxType) for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the bounds (elementwise). * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal( const mitk::BaseGeometry::BoundingBoxType *leftHandSide, const mitk::BaseGeometry::BoundingBoxType *rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two geometries for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the spacing, origin, axisvectors, extents, the matrix of the * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag. * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * If you want to use different tolarance values for different parts of the geometry, feel free to use * the other comparison methods and write your own implementation of Equal. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::BaseGeometry* leftHandSide, const mitk::BaseGeometry* rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo /** * @brief Equal A function comparing two transforms (TransformType) for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the IndexToWorldTransform (elementwise). * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::BaseGeometry::TransformType *leftHandSide, const mitk::BaseGeometry::TransformType *rightHandSide, mitk::ScalarType eps, bool verbose); //ToDo } // namespace mitk #endif /* BaseGeometry_H_HEADER_INCLUDED */ diff --git a/Core/Code/DataManagement/mitkGeometry2D.h b/Core/Code/DataManagement/mitkGeometry2D.h index 45bb00953a..10b62fcf89 100644 --- a/Core/Code/DataManagement/mitkGeometry2D.h +++ b/Core/Code/DataManagement/mitkGeometry2D.h @@ -1,273 +1,269 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ - #ifndef GEOMETRY2D_H_HEADER_INCLUDED_C1F4D8E0 #define GEOMETRY2D_H_HEADER_INCLUDED_C1F4D8E0 #include #include "mitkGeometry3D.h" namespace mitk { - /** * \brief Describes the geometry of a two-dimensional object * * Describes a two-dimensional manifold, i.e., to put it simply, * an object that can be described using a 2D coordinate-system. * * Geometry2D can map points between 3D world coordinates * (in mm) and the described 2D coordinate-system (in mm) by first projecting * the 3D point onto the 2D manifold and then calculating the 2D-coordinates * (in mm). These 2D-mm-coordinates can be further converted into * 2D-unit-coordinates (e.g., pixels), giving a parameter representation of * the object with parameter values inside a rectangle * (e.g., [0,0]..[width, height]), which is the bounding box (bounding range * in z-direction always [0]..[1]). * * A Geometry2D describes the 2D representation within a 3D object and is * therefore itself a Geometry3D (derived from Geometry3D). For example, * a single CT-image (slice) is 2D in the sense that you can access the * pixels using 2D-coordinates, but is also 3D, as the pixels are really * voxels, thus have an extension (thickness) in the 3rd dimension. * * Most often, instances of Geometry2D will be used to descibe a plane, * which is represented by the sub-class PlaneGeometry, but curved * surfaces are also possible. * * Optionally, a reference Geometry3D can be specified, which usually would * be the geometry associated with the underlying dataset. This is currently * used for calculating the intersection of inclined / rotated planes * (represented as Geometry2D) with the bounding box of the associated * Geometry3D. * * \warning The Geometry2Ds are not necessarily up-to-date and not even * initialized. As described in the previous paragraph, one of the * Generate-/Copy-/UpdateOutputInformation methods have to initialize it. * mitk::BaseData::GetGeometry2D() makes sure, that the Geometry2D is * up-to-date before returning it (by setting the update extent appropriately * and calling UpdateOutputInformation). * * Rule: everything is in mm (or ms for temporal information) if not * stated otherwise. * \ingroup Geometry */ class MITK_CORE_EXPORT Geometry2D : public mitk::Geometry3D { public: mitkClassMacro(Geometry2D, mitk::Geometry3D); itkNewMacro(Self); /** * \brief Project a 3D point given in mm (\a pt3d_mm) onto the 2D * geometry. The result is a 2D point in mm (\a pt2d_mm). * * The result is a 2D point in mm (\a pt2d_mm) relative to the upper-left * corner of the geometry. To convert this point into units (e.g., pixels * in case of an image), use WorldToIndex. * \return true projection was possible * \sa Project(const mitk::Point3D &pt3d_mm, mitk::Point3D * &projectedPt3d_mm) */ virtual bool Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const; /** * \brief Converts a 2D point given in mm (\a pt2d_mm) relative to the * upper-left corner of the geometry into the corresponding * world-coordinate (a 3D point in mm, \a pt3d_mm). * * To convert a 2D point given in units (e.g., pixels in case of an * image) into a 2D point given in mm (as required by this method), use * IndexToWorld. */ virtual void Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const; /** * \brief Convert a 2D point given in units (e.g., pixels in case of an * image) into a 2D point given in mm */ virtual void IndexToWorld( const mitk::Point2D &pt_units, mitk::Point2D &pt_mm) const; /** * \brief Convert a 2D point given in mm into a 2D point given in mm * (e.g., pixels in case of an image) */ virtual void WorldToIndex( const mitk::Point2D &pt_mm, mitk::Point2D &pt_units) const; /** * \brief Convert a 2D vector given in units (e.g., pixels in case of an * image) into a 2D vector given in mm * \warning strange: in contrast to vtkTransform the class itk::Transform * does not have the parameter, \em where the vector that is to be * transformed is located. This method here should also need this * information for general transforms. */ virtual void IndexToWorld( const mitk::Point2D &atPt2d_units, const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const; /** * \brief Convert a 2D vector given in mm into a 2D point vector in mm * (e.g., pixels in case of an image) * \warning strange: in contrast to vtkTransform the class itk::Transform * does not have the parameter, \em where the vector that is to be * transformed is located. This method here should also need this * information for general transforms. */ virtual void WorldToIndex( const mitk::Point2D &atPt2d_mm, const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const; /** * \brief Set the width and height of this 2D-geometry in units by calling * SetBounds. This does \a not change the extent in mm! * * For an image, this is the number of pixels in x-/y-direction. * \note In contrast to calling SetBounds directly, this does \a not change * the extent in mm! */ virtual void SetSizeInUnits(mitk::ScalarType width, mitk::ScalarType height); /** * \brief Project a 3D point given in mm (\a pt3d_mm) onto the 2D * geometry. The result is a 3D point in mm (\a projectedPt3d_mm). * * \return true projection was possible */ virtual bool Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const; /** * \brief Project a 3D vector given in mm (\a vec3d_mm) onto the 2D * geometry. The result is a 2D vector in mm (\a vec2d_mm). * * The result is a 2D vector in mm (\a vec2d_mm) relative to the * upper-left * corner of the geometry. To convert this point into units (e.g., pixels * in case of an image), use WorldToIndex. * \return true projection was possible * \sa Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D * &projectedVec3d_mm) */ virtual bool Map(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector2D &vec2d_mm) const; /** * \brief Converts a 2D vector given in mm (\a vec2d_mm) relative to the * upper-left corner of the geometry into the corresponding * world-coordinate (a 3D vector in mm, \a vec3d_mm). * * To convert a 2D vector given in units (e.g., pixels in case of an * image) into a 2D vector given in mm (as required by this method), use * IndexToWorld. */ virtual void Map(const mitk::Point2D & atPt2d_mm, const mitk::Vector2D &vec2d_mm, mitk::Vector3D &vec3d_mm) const; /** * \brief Project a 3D vector given in mm (\a vec3d_mm) onto the 2D * geometry. The result is a 3D vector in mm (\a projectedVec3d_mm). * * DEPRECATED. Use Project(vector,vector) instead * * \return true projection was possible */ virtual bool Project(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const; /** * \brief Project a 3D vector given in mm (\a vec3d_mm) onto the 2D * geometry. The result is a 3D vector in mm (\a projectedVec3d_mm). * * \return true projection was possible */ virtual bool Project( const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const; /** * \brief Distance of the point from the geometry * (bounding-box \em not considered) * */ inline ScalarType Distance(const Point3D& pt3d_mm) const { return fabs(SignedDistance(pt3d_mm)); } /** * \brief Signed distance of the point from the geometry * (bounding-box \em not considered) * */ virtual ScalarType SignedDistance(const Point3D& pt3d_mm) const; /** * \brief Test if the point is above the geometry * (bounding-box \em not considered) * */ virtual bool IsAbove(const Point3D& pt3d_mm) const; virtual void SetIndexToWorldTransform(mitk::AffineTransform3D* transform); virtual void SetExtentInMM(int direction, ScalarType extentInMM); virtual itk::LightObject::Pointer InternalClone() const; /** * \brief Set the geometrical frame of reference in which this Geometry2D * is placed. * * This would usually be the Geometry3D of the underlying dataset, but * setting it is optional. */ void SetReferenceGeometry( mitk::Geometry3D *geometry ); /** * \brief Get the geometrical frame of reference for this Geometry2D. */ Geometry3D *GetReferenceGeometry() const; bool HasReferenceGeometry() const; - protected: Geometry2D(); Geometry2D(const Geometry2D& other); virtual ~Geometry2D(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; /** * \brief factor to convert x-coordinates from mm to units and vice versa * */ mutable mitk::ScalarType m_ScaleFactorMMPerUnitX; /** * \brief factor to convert y-coordinates from mm to units and vice versa * */ mutable mitk::ScalarType m_ScaleFactorMMPerUnitY; mitk::Geometry3D *m_ReferenceGeometry; }; - } // namespace mitk #endif /* GEOMETRY2D_H_HEADER_INCLUDED_C1F4D8E0 */ diff --git a/Core/Code/DataManagement/mitkGeometry3D.h b/Core/Code/DataManagement/mitkGeometry3D.h index b39c8ad254..11c3d47b15 100644 --- a/Core/Code/DataManagement/mitkGeometry3D.h +++ b/Core/Code/DataManagement/mitkGeometry3D.h @@ -1,298 +1,299 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD #define GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD #include #include #include "itkScalableAffineTransform.h" #include #include "mitkBaseGeometry.h" class vtkLinearTransform; namespace mitk { //##Documentation //## @brief Standard typedef for time-bounds typedef itk::FixedArray TimeBounds; typedef itk::FixedArray FixedArrayType; //##Documentation //## @brief Standard 3D-BoundingBox typedef //## //## Standard 3D-BoundingBox typedef to get rid of template arguments (3D, type). typedef itk::BoundingBox BoundingBox; //##Documentation //## @brief Describes the geometry of a data object //## //## At least, it can return the bounding box of the data object. //## //## The class holds //## \li a bounding box which is axes-parallel in intrinsic coordinates //## (often integer indices of pixels), to be accessed by //## GetBoundingBox() //## \li a transform to convert intrinsic coordinates into a //## world-coordinate system with coordinates in millimeters //## and milliseconds (all are floating point values), to //## be accessed by GetIndexToWorldTransform() //## \li a life span, i.e. a bounding box in time in ms (with //## start and end time), to be accessed by GetTimeBounds(). //## The default is minus infinity to plus infinity. //## //## Geometry3D and its sub-classes allow converting between //## intrinsic coordinates (called index or unit coordinates) //## and world-coordinates (called world or mm coordinates), //## e.g. WorldToIndex. //## In case you need integer index coordinates, provide an //## mitk::Index3D (or itk::Index) as target variable to //## WorldToIndex, otherwise you will get a continuous index //## (floating point values). //## //## An important sub-class is SlicedGeometry3D, which descibes //## data objects consisting of slices, e.g., objects of type Image. //## Conversions between world coordinates (in mm) and unit coordinates //## (e.g., pixels in the case of an Image) can be performed. //## //## For more information on related classes, see \ref Geometry. //## //## Geometry3D instances referring to an Image need a slightly //## different definition of corners, see SetImageGeometry. This //## is usualy automatically called by Image. //## //## Geometry3D have to be initialized in the method GenerateOutputInformation() //## of BaseProcess (or CopyInformation/ UpdateOutputInformation of BaseData, //## if possible, e.g., by analyzing pic tags in Image) subclasses. See also //## itk::ProcessObject::GenerateOutputInformation(), //## itk::DataObject::CopyInformation() and //## itk::DataObject::UpdateOutputInformation(). //## //## Rule: everything is in mm (ms) if not stated otherwise. //## @ingroup Geometry class MITK_CORE_EXPORT Geometry3D : public BaseGeometry { public: mitkClassMacro(Geometry3D, mitk::BaseGeometry); + //void testXYZ(){ int a=1; }; //xxxxxxxxxxxxxxxx typedef itk::QuaternionRigidTransform< ScalarType > QuaternionTransformType; typedef QuaternionTransformType::VnlQuaternionType VnlQuaternionType; /** Method for creation through the object factory. */ itkNewMacro(Self); mitkNewMacro1Param(Self,Self); //##Documentation //## @brief When switching from an Image Geometry to a normal Geometry (and the other way around), you have to change the origin as well (See Geometry Documentation)! This function will change the "isImageGeometry" bool flag and changes the origin respectively. virtual void ChangeImageGeometryConsideringOriginOffset( const bool isAnImageGeometry ); //##Documentation //## @brief Get the position of the corner number \a id (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. Point3D GetCornerPoint(int id) const; //##Documentation //## @brief Get the position of a corner (in world coordinates) //## //## See SetImageGeometry for how a corner is defined on images. Point3D GetCornerPoint(bool xFront=true, bool yFront=true, bool zFront=true) const; //##Documentation //## @brief Is this an ImageGeometry? //## //## For more information, see SetImageGeometry itkGetConstMacro(ImageGeometry, bool); //##Documentation //## @brief Define that this Geometry3D is refering to an Image //## //## A geometry referring to an Image needs a slightly different //## definition of the position of the corners (see GetCornerPoint). //## The position of a voxel is defined by the position of its center. //## If we would use the origin (position of the (center of) the first //## voxel) as a corner and display this point, it would seem to be //## \em not at the corner but a bit within the image. Even worse for //## the opposite corner of the image: here the corner would appear //## outside the image (by half of the voxel diameter). Thus, we have //## to correct for this and to be able to do that, we need to know //## that the Geometry3D is referring to an Image. itkSetMacro(ImageGeometry, bool); itkBooleanMacro(ImageGeometry); //##Documentation //## @brief Test whether the point \a p ((continous!)index coordinates in units) is //## inside the bounding box bool IsIndexInside(const mitk::Point3D& index) const { bool inside = false; //if it is an image geometry, we need to convert the index to discrete values //this is done by applying the rounding function also used in WorldToIndex (see line 323) if (m_ImageGeometry) { mitk::Point3D discretIndex; discretIndex[0]=itk::Math::RoundHalfIntegerUp( index[0] ); discretIndex[1]=itk::Math::RoundHalfIntegerUp( index[1] ); discretIndex[2]=itk::Math::RoundHalfIntegerUp( index[2] ); inside = m_BoundingBox->IsInside(discretIndex); //we have to check if the index is at the upper border of each dimension, // because the boundingbox is not centerbased if (inside) { const BoundingBox::BoundsArrayType& bounds = m_BoundingBox->GetBounds(); if((discretIndex[0] == bounds[1]) || (discretIndex[1] == bounds[3]) || (discretIndex[2] == bounds[5])) inside = false; } } else inside = m_BoundingBox->IsInside(index); return inside; } //##Documentation //## @brief Convenience method for working with ITK indices template bool IsIndexInside(const itk::Index &index) const { int i, dim=index.GetIndexDimension(); Point3D pt_index; pt_index.Fill(0); for ( i = 0; i < dim; ++i ) { pt_index[i] = index[i]; } return IsIndexInside(pt_index); } //##Documentation //## @brief clones the geometry //## //## Overwrite in all sub-classes. //## Normally looks like: //## \code //## Self::Pointer newGeometry = new Self(*this); //## newGeometry->UnRegister(); //## return newGeometry.GetPointer(); //## \endcode virtual itk::LightObject::Pointer InternalClone() const; //Umzug: //##Documentation //## @brief Get the parametric bounding-box //## //## See AbstractTransformGeometry for an example usage of this. itkGetConstObjectMacro(ParametricBoundingBox, BoundingBox); //##Documentation //## @brief Get the parametric bounds //## //## See AbstractTransformGeometry for an example usage of this. const BoundingBox::BoundsArrayType& GetParametricBounds() const { assert(m_ParametricBoundingBox.IsNotNull()); return m_ParametricBoundingBox->GetBounds(); } //##Documentation //## @brief Get the parametric extent //## //## See AbstractTransformGeometry for an example usage of this. mitk::ScalarType GetParametricExtent(int direction) const { if (direction < 0 || direction>=3) mitkThrow() << "Invalid direction. Must be between either 0, 1 or 2. "; assert(m_ParametricBoundingBox.IsNotNull()); BoundingBoxType::BoundsArrayType bounds = m_ParametricBoundingBox->GetBounds(); return bounds[direction*2+1]-bounds[direction*2]; } //##Documentation //## @brief Get the parametric extent in mm //## //## See AbstractTransformGeometry for an example usage of this. virtual mitk::ScalarType GetParametricExtentInMM(int direction) const { return GetExtentInMM(direction); } //##Documentation //## @brief Get the parametric transform //## //## See AbstractTransformGeometry for an example usage of this. virtual const Transform3D* GetParametricTransform() const { return m_IndexToWorldTransform; } protected: Geometry3D(); Geometry3D(const Geometry3D& other); virtual ~Geometry3D(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; //##Documentation //## @brief Set the parametric bounds //## //## Protected in this class, made public in some sub-classes, e.g., //## ExternAbstractTransformGeometry. virtual void SetParametricBounds(const BoundingBox::BoundsArrayType& bounds); mutable mitk::BoundingBox::Pointer m_ParametricBoundingBox; bool m_ImageGeometry; virtual void InternPostInitialize(); virtual void InternPostInitializeGeometry(Geometry3D* newGeometry) const; static const std::string INDEX_TO_OBJECT_TRANSFORM; static const std::string OBJECT_TO_NODE_TRANSFORM; static const std::string INDEX_TO_NODE_TRANSFORM; static const std::string INDEX_TO_WORLD_TRANSFORM; }; // // Static compare functions mainly for testing // /** * @brief Equal A function comparing two geometries for beeing identical. * * @ingroup MITKTestingAPI * * The function compares the spacing, origin, axisvectors, extents, the matrix of the * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag. * * The parameter eps is a tolarence value for all methods which are internally used for comparion. * If you want to use different tolarance values for different parts of the geometry, feel free to use * the other comparison methods and write your own implementation of Equal. * @param rightHandSide Compare this against leftHandSide. * @param leftHandSide Compare this against rightHandSide. * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return True, if all comparison are true. False in any other case. */ MITK_CORE_EXPORT bool Equal(const mitk::Geometry3D* leftHandSide, const mitk::Geometry3D* rightHandSide, ScalarType eps, bool verbose); } // namespace mitk #endif /* GEOMETRY3D_H_HEADER_INCLUDED_C1EBD0AD */ diff --git a/Core/Code/DataManagement/mitkPlaneGeometry.cpp b/Core/Code/DataManagement/mitkPlaneGeometry.cpp index 364fe98271..560b029def 100644 --- a/Core/Code/DataManagement/mitkPlaneGeometry.cpp +++ b/Core/Code/DataManagement/mitkPlaneGeometry.cpp @@ -1,777 +1,735 @@ /*=================================================================== 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 { + mitk::PlaneGeometry::PlaneGeometry() + { + Initialize(); + } -mitk::PlaneGeometry::PlaneGeometry() -{ - Initialize(); -} - - -mitk::PlaneGeometry::~PlaneGeometry() -{ -} - - -void -PlaneGeometry::Initialize() -{ - Superclass::Initialize(); -} - - -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::SetIndexToWorldTransform(mitk::AffineTransform3D *transform) -{ - EnsurePerpendicularNormal(transform); - - Superclass::SetIndexToWorldTransform(transform); -} - + mitk::PlaneGeometry::~PlaneGeometry() + { + } -void -PlaneGeometry::SetBounds(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::Initialize() + { + Superclass::Initialize(); + } - Superclass::SetBounds(bounds); -} + 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::SetIndexToWorldTransform(mitk::AffineTransform3D *transform) + { + EnsurePerpendicularNormal(transform); -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]; -} + Superclass::SetIndexToWorldTransform(transform); + } + void + PlaneGeometry::InternPreSetBounds(const BoundingBox::BoundsArrayType &bounds) + { + //currently the unit rectangle must be starting at the origin [0,0] + assert(bounds[0]==0); + assert(bounds[2]==0); + //the unit rectangle must be two-dimensional + assert(bounds[1]>0); + assert(bounds[3]>0); + } -void PlaneGeometry::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 &pt_units, Point2D &pt_mm ) const + { + pt_mm[0]=m_ScaleFactorMMPerUnitX*pt_units[0]; + pt_mm[1]=m_ScaleFactorMMPerUnitY*pt_units[1]; + } -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::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 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::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::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::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 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::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(); + 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); + } - //construct standard view - Point3D origin; - VnlVector rightDV(3), bottomDV(3); - origin.Fill(0); - int normalDirection; - switch(planeorientation) + 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 ); - } + } + 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 ); + 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() + if ( transform == NULL ) + { + this->SetMatrixByVectors( rightDV, bottomDV ); + } + else + { + this->SetMatrixByVectors( + rightDV, bottomDV, + transform->GetMatrix().GetVnlMatrix() .get_column(normalDirection).magnitude() - ); - } - - this->SetOrigin(origin); -} + ); + } + this->SetOrigin(origin); + } -void -PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, - PlaneOrientation planeorientation, ScalarType zPosition, - bool frontside, bool rotated ) -{ - this->SetReferenceGeometry( const_cast< Geometry3D * >( geometry3D ) ); + void + PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, + PlaneOrientation planeorientation, ScalarType zPosition, + bool frontside, bool rotated ) + { + this->SetReferenceGeometry( const_cast< Geometry3D * >( geometry3D ) ); - ScalarType width, height; + ScalarType width, height; - const BoundingBox::BoundsArrayType& boundsarray = - geometry3D->GetBoundingBox()->GetBounds(); + const BoundingBox::BoundsArrayType& boundsarray = + geometry3D->GetBoundingBox()->GetBounds(); - Vector3D originVector; - FillVector3D(originVector, boundsarray[0], boundsarray[2], boundsarray[4]); + 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) - { + 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 ); + InitializeStandardPlane( width, height, + geometry3D->GetIndexToWorldTransform(), + planeorientation, zPosition, frontside, rotated ); - origin = GetOrigin() + originVector; - SetOrigin(origin); -} + ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; + this->SetBounds( bounds ); + Point3D origin; + originVector = geometry3D->GetIndexToWorldTransform() + ->TransformVector( originVector ); -void -PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, - bool top, PlaneOrientation planeorientation, bool frontside, bool rotated ) -{ - ScalarType zPosition; + origin = GetOrigin() + originVector; + SetOrigin(origin); + } - switch(planeorientation) + void + PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, + bool top, PlaneOrientation planeorientation, bool frontside, bool rotated ) { + ScalarType zPosition; + + switch(planeorientation) + { case Axial: zPosition = (top ? 0.5 : geometry3D->GetExtent(2)-1+0.5); break; case Frontal: zPosition = (top ? 0.5 : geometry3D->GetExtent(1)-1+0.5); break; case Sagittal: zPosition = (top ? 0.5 : geometry3D->GetExtent(0)-1+0.5); break; default: itkExceptionMacro("unknown PlaneOrientation"); - } - - InitializeStandardPlane( geometry3D, planeorientation, - zPosition, frontside, rotated ); -} - - -void -PlaneGeometry::InitializeStandardPlane( const Vector3D &rightVector, - const Vector3D &downVector, const Vector3D *spacing ) -{ - InitializeStandardPlane( rightVector.GetVnlVector(), - downVector.GetVnlVector(), spacing ); -} - - -void -PlaneGeometry::InitializeStandardPlane( const VnlVector& rightVector, - const VnlVector &downVector, const Vector3D *spacing ) -{ - ScalarType width = rightVector.magnitude(); - ScalarType height = downVector.magnitude(); - - InitializeStandardPlane( width, height, rightVector, downVector, spacing ); -} - - -void -PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, - ScalarType height, const Vector3D &rightVector, const Vector3D &downVector, - const Vector3D *spacing ) -{ - InitializeStandardPlane( - width, height, - rightVector.GetVnlVector(), downVector.GetVnlVector(), - spacing ); -} - - -void -PlaneGeometry::InitializeStandardPlane( - mitk::ScalarType width, ScalarType height, - const VnlVector &rightVector, const VnlVector &downVector, - const Vector3D *spacing ) -{ - assert(width > 0); - assert(height > 0); + } - VnlVector rightDV = rightVector; rightDV.normalize(); - VnlVector downDV = downVector; downDV.normalize(); - VnlVector normal = vnl_cross_3d(rightVector, downVector); - normal.normalize(); + InitializeStandardPlane( geometry3D, planeorientation, + zPosition, frontside, rotated ); + } - if(spacing!=NULL) + void + PlaneGeometry::InitializeStandardPlane( const Vector3D &rightVector, + const Vector3D &downVector, const Vector3D *spacing ) { - rightDV *= (*spacing)[0]; - downDV *= (*spacing)[1]; - normal *= (*spacing)[2]; + InitializeStandardPlane( rightVector.GetVnlVector(), + downVector.GetVnlVector(), spacing ); } - AffineTransform3D::Pointer transform = AffineTransform3D::New(); - Matrix3D matrix; - matrix.GetVnlMatrix().set_column(0, rightDV); - matrix.GetVnlMatrix().set_column(1, downDV); - matrix.GetVnlMatrix().set_column(2, normal); - transform->SetMatrix(matrix); - transform->SetOffset(m_IndexToWorldTransform->GetOffset()); - - 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 ) + void + PlaneGeometry::InitializeStandardPlane( const VnlVector& rightVector, + const VnlVector &downVector, const Vector3D *spacing ) { - FillVector3D( rightVectorVnl, 1.0f, -normal[0]/normal[1], 0.0f ); - rightVectorVnl.normalize(); + ScalarType width = rightVector.magnitude(); + ScalarType height = downVector.magnitude(); + + InitializeStandardPlane( width, height, rightVector, downVector, spacing ); } - else + + void + PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, + ScalarType height, const Vector3D &rightVector, const Vector3D &downVector, + const Vector3D *spacing ) { - FillVector3D( rightVectorVnl, 0.0f, 1.0f, 0.0f ); + InitializeStandardPlane( + width, height, + rightVector.GetVnlVector(), downVector.GetVnlVector(), + spacing ); } - downVectorVnl = vnl_cross_3d( normal.GetVnlVector(), rightVectorVnl ); - downVectorVnl.normalize(); - - InitializeStandardPlane( rightVectorVnl, downVectorVnl ); - - SetOrigin(origin); -} - - -void -PlaneGeometry::SetMatrixByVectors( const VnlVector &rightVector, - const VnlVector &downVector, ScalarType thickness ) -{ - VnlVector normal = vnl_cross_3d(rightVector, downVector); - normal.normalize(); - normal *= thickness; - - AffineTransform3D::Pointer transform = AffineTransform3D::New(); - Matrix3D matrix; - matrix.GetVnlMatrix().set_column(0, rightVector); - matrix.GetVnlMatrix().set_column(1, downVector); - matrix.GetVnlMatrix().set_column(2, normal); - transform->SetMatrix(matrix); - transform->SetOffset(m_IndexToWorldTransform->GetOffset()); - SetIndexToWorldTransform(transform); -} - - -Vector3D -PlaneGeometry::GetNormal() const -{ - Vector3D frontToBack; - frontToBack.SetVnlVector( m_IndexToWorldTransform - ->GetMatrix().GetVnlMatrix().get_column(2) ); - - return frontToBack; -} + void + PlaneGeometry::InitializeStandardPlane( + mitk::ScalarType width, ScalarType height, + const VnlVector &rightVector, const VnlVector &downVector, + const Vector3D *spacing ) + { + assert(width > 0); + assert(height > 0); -VnlVector -PlaneGeometry::GetNormalVnl() const -{ - return m_IndexToWorldTransform - ->GetMatrix().GetVnlMatrix().get_column(2); -} + 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]; + } -ScalarType -PlaneGeometry::DistanceFromPlane( const Point3D &pt3d_mm ) const -{ - return fabs(SignedDistance( pt3d_mm )); -} + AffineTransform3D::Pointer transform = AffineTransform3D::New(); + Matrix3D matrix; + matrix.GetVnlMatrix().set_column(0, rightDV); + matrix.GetVnlMatrix().set_column(1, downDV); + matrix.GetVnlMatrix().set_column(2, normal); + transform->SetMatrix(matrix); + transform->SetOffset(m_IndexToWorldTransform->GetOffset()); + ScalarType bounds[6] = { 0, width, 0, height, 0, 1 }; + this->SetBounds( bounds ); -ScalarType -PlaneGeometry::SignedDistance( const Point3D &pt3d_mm ) const -{ - return SignedDistanceFromPlane(pt3d_mm); -} + this->SetIndexToWorldTransform( transform ); + } + void + PlaneGeometry::InitializePlane( const Point3D &origin, + const Vector3D &normal ) + { + VnlVector rightVectorVnl(3), downVectorVnl; -bool -PlaneGeometry::IsAbove( const Point3D &pt3d_mm ) const -{ - return SignedDistanceFromPlane(pt3d_mm) > 0; -} + 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 ); -bool -PlaneGeometry::IntersectionLine( - const PlaneGeometry* plane, Line3D& crossline ) const -{ - Vector3D normal = this->GetNormal(); - normal.Normalize(); + SetOrigin(origin); + } - Vector3D planeNormal = plane->GetNormal(); - planeNormal.Normalize(); + void + PlaneGeometry::SetMatrixByVectors( const VnlVector &rightVector, + const VnlVector &downVector, ScalarType thickness ) + { + VnlVector normal = vnl_cross_3d(rightVector, downVector); + normal.normalize(); + normal *= thickness; + + AffineTransform3D::Pointer transform = AffineTransform3D::New(); + Matrix3D matrix; + matrix.GetVnlMatrix().set_column(0, rightVector); + matrix.GetVnlMatrix().set_column(1, downVector); + matrix.GetVnlMatrix().set_column(2, normal); + transform->SetMatrix(matrix); + transform->SetOffset(m_IndexToWorldTransform->GetOffset()); + SetIndexToWorldTransform(transform); + } - Vector3D direction = itk::CrossProduct( normal, planeNormal ); + Vector3D + PlaneGeometry::GetNormal() const + { + Vector3D frontToBack; + frontToBack.SetVnlVector( m_IndexToWorldTransform + ->GetMatrix().GetVnlMatrix().get_column(2) ); - if ( direction.GetSquaredNorm() < eps ) - return false; + return frontToBack; + } - crossline.SetDirection( direction ); + VnlVector + PlaneGeometry::GetNormalVnl() const + { + return m_IndexToWorldTransform + ->GetMatrix().GetVnlMatrix().get_column(2); + } - double N1dN2 = normal * planeNormal; - double determinant = 1.0 - N1dN2 * N1dN2; + ScalarType + PlaneGeometry::DistanceFromPlane( const Point3D &pt3d_mm ) const + { + return fabs(SignedDistance( pt3d_mm )); + } - Vector3D origin = this->GetOrigin().GetVectorFromOrigin(); - Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin(); + ScalarType + PlaneGeometry::SignedDistance( const Point3D &pt3d_mm ) const + { + return SignedDistanceFromPlane(pt3d_mm); + } - double d1 = normal * origin; - double d2 = planeNormal * planeOrigin; + bool + PlaneGeometry::IsAbove( const Point3D &pt3d_mm ) const + { + return SignedDistanceFromPlane(pt3d_mm) > 0; + } - double c1 = ( d1 - d2 * N1dN2 ) / determinant; - double c2 = ( d2 - d1 * N1dN2 ) / determinant; + bool + PlaneGeometry::IntersectionLine( + const PlaneGeometry* plane, Line3D& crossline ) const + { + Vector3D normal = this->GetNormal(); + normal.Normalize(); - Vector3D p = normal * c1 + planeNormal * c2; - crossline.GetPoint().GetVnlVector() = p.GetVnlVector(); + Vector3D planeNormal = plane->GetNormal(); + planeNormal.Normalize(); - return true; -} + Vector3D direction = itk::CrossProduct( normal, planeNormal ); + if ( direction.GetSquaredNorm() < eps ) + return false; -unsigned int -PlaneGeometry::IntersectWithPlane2D( - const PlaneGeometry* plane, Point2D& lineFrom, Point2D &lineTo ) const -{ - Line3D crossline; - if ( this->IntersectionLine( plane, crossline ) == false ) - return 0; + crossline.SetDirection( direction ); - Point2D point2; - Vector2D direction2; + double N1dN2 = normal * planeNormal; + double determinant = 1.0 - N1dN2 * N1dN2; - this->Map( crossline.GetPoint(), point2 ); - this->Map( crossline.GetPoint(), crossline.GetDirection(), direction2 ); + Vector3D origin = this->GetOrigin().GetVectorFromOrigin(); + Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin(); - return - Line3D::RectangleLineIntersection( - 0, 0, GetExtentInMM(0), GetExtentInMM(1), - point2, direction2, lineFrom, lineTo ); -} + double d1 = normal * origin; + double d2 = planeNormal * planeOrigin; + double c1 = ( d1 - d2 * N1dN2 ) / determinant; + double c2 = ( d2 - d1 * N1dN2 ) / determinant; -double PlaneGeometry::Angle( const PlaneGeometry *plane ) const -{ - return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2)); -} + Vector3D p = normal * c1 + planeNormal * c2; + crossline.GetPoint().GetVnlVector() = p.GetVnlVector(); + return true; + } -double PlaneGeometry::Angle( const Line3D &line ) const -{ - return vnl_math::pi_over_2 - - angle( line.GetDirection().GetVnlVector(), GetMatrixColumn(2) ); -} + 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; -bool PlaneGeometry::IntersectionPoint( - const Line3D &line, Point3D &intersectionPoint ) const -{ - Vector3D planeNormal = this->GetNormal(); - planeNormal.Normalize(); + this->Map( crossline.GetPoint(), point2 ); + this->Map( crossline.GetPoint(), crossline.GetDirection(), direction2 ); - Vector3D lineDirection = line.GetDirection(); - lineDirection.Normalize(); + return + Line3D::RectangleLineIntersection( + 0, 0, GetExtentInMM(0), GetExtentInMM(1), + point2, direction2, lineFrom, lineTo ); + } - double t = planeNormal * lineDirection; - if ( fabs( t ) < eps ) + double PlaneGeometry::Angle( const PlaneGeometry *plane ) const { - return false; + return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2)); } - Vector3D diff; - diff = this->GetOrigin() - line.GetPoint(); - t = ( planeNormal * diff ) / t; - - intersectionPoint = line.GetPoint() + lineDirection * t; - return true; -} + 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(); -bool -PlaneGeometry::IntersectionPointParam( const Line3D &line, double &t ) const -{ - Vector3D planeNormal = this->GetNormal(); + Vector3D lineDirection = line.GetDirection(); + lineDirection.Normalize(); - Vector3D lineDirection = line.GetDirection(); + double t = planeNormal * lineDirection; + if ( fabs( t ) < eps ) + { + return false; + } - t = planeNormal * lineDirection; + Vector3D diff; + diff = this->GetOrigin() - line.GetPoint(); + t = ( planeNormal * diff ) / t; - if ( fabs( t ) < eps ) - { - return false; + intersectionPoint = line.GetPoint() + lineDirection * t; + return true; } - 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::IntersectionPointParam( const Line3D &line, double &t ) const + { + Vector3D planeNormal = this->GetNormal(); -bool -PlaneGeometry::IsOnPlane( const Point3D &point ) const -{ - return Distance(point) < eps; -} + Vector3D lineDirection = line.GetDirection(); + t = planeNormal * lineDirection; -bool -PlaneGeometry::IsOnPlane( const Line3D &line ) const -{ - return ( (Distance( line.GetPoint() ) < eps) - && (Distance( line.GetPoint2() ) < eps) ); -} + if ( fabs( t ) < eps ) + { + return false; + } + Vector3D diff; + diff = this->GetOrigin() - line.GetPoint(); + t = ( planeNormal * diff ) / t; + return true; + } -bool -PlaneGeometry::IsOnPlane( const PlaneGeometry *plane ) const -{ - return ( IsParallel( plane ) && (Distance( plane->GetOrigin() ) < eps) ); -} + 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; + } -Point3D -PlaneGeometry::ProjectPointOntoPlane( const Point3D& pt ) const -{ - ScalarType len = this->GetNormalVnl().two_norm(); - return pt - this->GetNormal() * this->SignedDistanceFromPlane( pt ) / len; -} + 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) ); + } -itk::LightObject::Pointer -PlaneGeometry::InternalClone() const -{ - Self::Pointer newGeometry = new PlaneGeometry(*this); - newGeometry->UnRegister(); - return newGeometry.GetPointer(); -} + Point3D + PlaneGeometry::ProjectPointOntoPlane( const Point3D& pt ) const + { + ScalarType len = this->GetNormalVnl().two_norm(); + return pt - this->GetNormal() * this->SignedDistanceFromPlane( pt ) / len; + } -void -PlaneGeometry::ExecuteOperation( Operation *operation ) -{ - vtkTransform *transform = vtkTransform::New(); - transform->SetMatrix( m_VtkMatrix ); + itk::LightObject::Pointer + PlaneGeometry::InternalClone() const + { + Self::Pointer newGeometry = new PlaneGeometry(*this); + newGeometry->UnRegister(); + return newGeometry.GetPointer(); + } - switch ( operation->GetOperationType() ) + void + PlaneGeometry::ExecuteOperation( Operation *operation ) { - case OpORIENT: + vtkTransform *transform = vtkTransform::New(); + transform->SetMatrix( m_VtkMatrix ); + + switch ( operation->GetOperationType() ) { - mitk::PlaneOperation *planeOp = dynamic_cast< mitk::PlaneOperation * >( operation ); - if ( planeOp == NULL ) + case OpORIENT: { - return; - } + mitk::PlaneOperation *planeOp = dynamic_cast< mitk::PlaneOperation * >( operation ); + if ( planeOp == NULL ) + { + return; + } - Point3D center = planeOp->GetPoint(); + Point3D center = planeOp->GetPoint(); - Vector3D orientationVector = planeOp->GetNormal(); - Vector3D defaultVector; - FillVector3D( defaultVector, 0.0, 0.0, 1.0 ); + 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() ); + 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; + 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) + 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; } - - 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(); + default: + Superclass::ExecuteOperation( operation ); transform->Delete(); return; - } - default: - Superclass::ExecuteOperation( operation ); + + m_VtkMatrix->DeepCopy(transform->GetMatrix()); + this->TransferVtkToItkTransform(); + this->Modified(); transform->Delete(); - return; } - m_VtkMatrix->DeepCopy(transform->GetMatrix()); - this->TransferVtkToItkTransform(); - this->Modified(); - transform->Delete(); -} - -void PlaneGeometry::PrintSelf( std::ostream& os, itk::Indent indent ) const -{ - Superclass::PrintSelf(os,indent); - os << indent << " Normal: " << GetNormal() << std::endl; -} - - + void PlaneGeometry::PrintSelf( std::ostream& os, itk::Indent indent ) const + { + Superclass::PrintSelf(os,indent); + os << indent << " Normal: " << GetNormal() << std::endl; + } } // namespace diff --git a/Core/Code/DataManagement/mitkPlaneGeometry.h b/Core/Code/DataManagement/mitkPlaneGeometry.h index 82fd504f1b..b95ec89c37 100644 --- a/Core/Code/DataManagement/mitkPlaneGeometry.h +++ b/Core/Code/DataManagement/mitkPlaneGeometry.h @@ -1,434 +1,394 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ - #ifndef PLANEGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #define PLANEGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #include #include "mitkGeometry2D.h" #include "mitkRestorePlanePositionOperation.h" #include namespace mitk { + template < class TCoordRep, unsigned int NPointDimension > class Line; + typedef Line Line3D; -template < class TCoordRep, unsigned int NPointDimension > class Line; -typedef Line Line3D; - - - -/** - * \brief Describes a two-dimensional, rectangular plane - * - * \ingroup Geometry - */ -class MITK_CORE_EXPORT PlaneGeometry : public Geometry2D -{ -public: - mitkClassMacro(PlaneGeometry,Geometry2D); + /** + * \brief Describes a two-dimensional, rectangular plane + * + * \ingroup Geometry + */ + class MITK_CORE_EXPORT PlaneGeometry : public Geometry2D + { + public: + mitkClassMacro(PlaneGeometry,Geometry2D); - /** Method for creation through the object factory. */ - itkNewMacro(Self); + /** Method for creation through the object factory. */ + itkNewMacro(Self); - enum PlaneOrientation - { + enum PlaneOrientation + { #ifdef _MSC_VER - Transversal, // deprecated + Transversal, // deprecated #endif - Axial = 0, - Sagittal, - Frontal - }; + Axial = 0, + Sagittal, + Frontal + }; #ifdef __GNUC__ - __attribute__ ((deprecated)) static const PlaneOrientation Transversal = PlaneOrientation(Axial); + __attribute__ ((deprecated)) static const PlaneOrientation Transversal = PlaneOrientation(Axial); #endif - virtual void IndexToWorld(const Point2D &pt_units, Point2D &pt_mm) const; - - virtual void WorldToIndex(const Point2D &pt_mm, Point2D &pt_units) const; - - - //##Documentation - //## @brief Convert (continuous or discrete) index coordinates of a \em vector - //## \a vec_units to world coordinates (in mm) - //## @deprecated First parameter (Point2D) is not used. If possible, please use void IndexToWorld(const mitk::Vector2D& vec_units, mitk::Vector2D& vec_mm) const. - //## For further information about coordinates types, please see the Geometry documentation - virtual void IndexToWorld(const mitk::Point2D &atPt2d_untis, const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const; - - //##Documentation - //## @brief Convert (continuous or discrete) index coordinates of a \em vector - //## \a vec_units to world coordinates (in mm) - //## For further information about coordinates types, please see the Geometry documentation - virtual void IndexToWorld(const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const; - - //##Documentation - //## @brief Convert world coordinates (in mm) of a \em vector - //## \a vec_mm to (continuous!) index coordinates. - //## @deprecated First parameter (Point2D) is not used. If possible, please use void WorldToIndex(const mitk::Vector2D& vec_mm, mitk::Vector2D& vec_units) const. - //## For further information about coordinates types, please see the Geometry documentation - virtual void WorldToIndex(const mitk::Point2D &atPt2d_mm, const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const; - - //##Documentation - //## @brief Convert world coordinates (in mm) of a \em vector - //## \a vec_mm to (continuous!) index coordinates. - //## For further information about coordinates types, please see the Geometry documentation - virtual void WorldToIndex(const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const; - - - virtual void Initialize(); - - /** - * \brief Initialize a plane with orientation \a planeorientation - * (default: axial) with respect to \a geometry3D (default: identity). - * Spacing also taken from \a geometry3D. - * - * \warning A former version of this method created a geometry with unit - * spacing. For unit spacing use - * - * \code - * // for in-plane unit spacing: - * thisgeometry->SetSizeInUnits(thisgeometry->GetExtentInMM(0), - * thisgeometry->GetExtentInMM(1)); - * // additionally, for unit spacing in normal direction (former version - * // did not do this): - * thisgeometry->SetExtentInMM(2, 1.0); - * \endcode - */ - virtual void InitializeStandardPlane( const Geometry3D* geometry3D, - PlaneOrientation planeorientation = Axial, ScalarType zPosition = 0, - bool frontside=true, bool rotated=false ); - - - /** - * \brief Initialize a plane with orientation \a planeorientation - * (default: axial) with respect to \a geometry3D (default: identity). - * Spacing also taken from \a geometry3D. - * - * \param top if \a true, create plane at top, otherwise at bottom - * (for PlaneOrientation Axial, for other plane locations respectively) - */ - virtual void InitializeStandardPlane( const Geometry3D* geometry3D, bool top, - PlaneOrientation planeorientation = Axial, - bool frontside=true, bool rotated=false ); - - - /** - * \brief Initialize a plane with orientation \a planeorientation - * (default: axial) with respect to \a transform (default: identity) - * given width and height in units. - * - */ - virtual void InitializeStandardPlane( ScalarType width, ScalarType height, - const AffineTransform3D* transform = NULL, - PlaneOrientation planeorientation = Axial, - ScalarType zPosition = 0, bool frontside=true, bool rotated=false ); - - - /** - * \brief Initialize plane with orientation \a planeorientation - * (default: axial) given width, height and spacing. - * - */ - virtual void InitializeStandardPlane( ScalarType width, ScalarType height, - const Vector3D & spacing, PlaneOrientation planeorientation = Axial, - ScalarType zPosition = 0, bool frontside = true, bool rotated = false ); - - /** - * \brief Initialize plane by width and height in pixels, right-/down-vector - * (itk) to describe orientation in world-space (vectors will be normalized) - * and spacing (default: 1.0 mm in all directions). - * - * The vectors are normalized and multiplied by the respective spacing before - * they are set in the matrix. - */ - virtual void InitializeStandardPlane( ScalarType width, ScalarType height, - const Vector3D& rightVector, const Vector3D& downVector, - const Vector3D *spacing = NULL ); - - - /** - * \brief Initialize plane by width and height in pixels, - * right-/down-vector (vnl) to describe orientation in world-space (vectors - * will be normalized) and spacing (default: 1.0 mm in all directions). - * - * The vectors are normalized and multiplied by the respective spacing - * before they are set in the matrix. - */ - virtual void InitializeStandardPlane( ScalarType width, ScalarType height, - const VnlVector& rightVector, const VnlVector& downVector, - const Vector3D * spacing = NULL ); - - - /** - * \brief Initialize plane by right-/down-vector (itk) and spacing - * (default: 1.0 mm in all directions). - * - * The length of the right-/-down-vector is used as width/height in units, - * respectively. Then, the vectors are normalized and multiplied by the - * respective spacing before they are set in the matrix. - */ - virtual void InitializeStandardPlane( const Vector3D& rightVector, - const Vector3D& downVector, const Vector3D * spacing = NULL ); - - - /** - * \brief Initialize plane by right-/down-vector (vnl) and spacing - * (default: 1.0 mm in all directions). - * - * The length of the right-/-down-vector is used as width/height in units, - * respectively. Then, the vectors are normalized and multiplied by the - * respective spacing before they are set in the matrix. - */ - virtual void InitializeStandardPlane( const VnlVector& rightVector, - const VnlVector& downVector, const Vector3D * spacing = NULL ); - - /** - * \brief Initialize plane by origin and normal (size is 1.0 mm in - * all directions, direction of right-/down-vector valid but - * undefined). - * - */ - virtual void InitializePlane( const Point3D& origin, const Vector3D& normal); - - /** - * \brief Initialize plane by right-/down-vector. - * - * \warning The vectors are set into the matrix as they are, - * \em without normalization! - */ - void SetMatrixByVectors( const VnlVector& rightVector, - const VnlVector& downVector, ScalarType thickness=1.0 ); - - - /** - * \brief Change \a transform so that the third column of the - * transform-martix is perpendicular to the first two columns - * - */ - static void EnsurePerpendicularNormal( AffineTransform3D* transform ); - - - /** - * \brief Normal of the plane - * - */ - Vector3D GetNormal() const; - - - /** - * \brief Normal of the plane as VnlVector - * - */ - VnlVector GetNormalVnl() const; - - - virtual ScalarType SignedDistance( const Point3D& pt3d_mm ) const; - - - virtual bool IsAbove( const Point3D& pt3d_mm ) const; - - - /** - * \brief Distance of the point from the plane - * (bounding-box \em not considered) - * - */ - ScalarType DistanceFromPlane( const Point3D& pt3d_mm ) const ; - - - /** - * \brief Signed distance of the point from the plane - * (bounding-box \em not considered) - * - * > 0 : point is in the direction of the direction vector. - */ - inline ScalarType SignedDistanceFromPlane( const Point3D& pt3d_mm ) const - { - ScalarType len = GetNormalVnl().two_norm(); - - if( len == 0 ) - return 0; - - return (pt3d_mm-GetOrigin())*GetNormal() / len; - } - + virtual void IndexToWorld(const Point2D &pt_units, Point2D &pt_mm) const; + + virtual void WorldToIndex(const Point2D &pt_mm, Point2D &pt_units) const; + + //##Documentation + //## @brief Convert (continuous or discrete) index coordinates of a \em vector + //## \a vec_units to world coordinates (in mm) + //## @deprecated First parameter (Point2D) is not used. If possible, please use void IndexToWorld(const mitk::Vector2D& vec_units, mitk::Vector2D& vec_mm) const. + //## For further information about coordinates types, please see the Geometry documentation + virtual void IndexToWorld(const mitk::Point2D &atPt2d_untis, const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const; + + //##Documentation + //## @brief Convert (continuous or discrete) index coordinates of a \em vector + //## \a vec_units to world coordinates (in mm) + //## For further information about coordinates types, please see the Geometry documentation + virtual void IndexToWorld(const mitk::Vector2D &vec_units, mitk::Vector2D &vec_mm) const; + + //##Documentation + //## @brief Convert world coordinates (in mm) of a \em vector + //## \a vec_mm to (continuous!) index coordinates. + //## @deprecated First parameter (Point2D) is not used. If possible, please use void WorldToIndex(const mitk::Vector2D& vec_mm, mitk::Vector2D& vec_units) const. + //## For further information about coordinates types, please see the Geometry documentation + virtual void WorldToIndex(const mitk::Point2D &atPt2d_mm, const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const; + + //##Documentation + //## @brief Convert world coordinates (in mm) of a \em vector + //## \a vec_mm to (continuous!) index coordinates. + //## For further information about coordinates types, please see the Geometry documentation + virtual void WorldToIndex(const mitk::Vector2D &vec_mm, mitk::Vector2D &vec_units) const; + + virtual void Initialize(); + + /** + * \brief Initialize a plane with orientation \a planeorientation + * (default: axial) with respect to \a geometry3D (default: identity). + * Spacing also taken from \a geometry3D. + * + * \warning A former version of this method created a geometry with unit + * spacing. For unit spacing use + * + * \code + * // for in-plane unit spacing: + * thisgeometry->SetSizeInUnits(thisgeometry->GetExtentInMM(0), + * thisgeometry->GetExtentInMM(1)); + * // additionally, for unit spacing in normal direction (former version + * // did not do this): + * thisgeometry->SetExtentInMM(2, 1.0); + * \endcode + */ + virtual void InitializeStandardPlane( const Geometry3D* geometry3D, + PlaneOrientation planeorientation = Axial, ScalarType zPosition = 0, + bool frontside=true, bool rotated=false ); + + /** + * \brief Initialize a plane with orientation \a planeorientation + * (default: axial) with respect to \a geometry3D (default: identity). + * Spacing also taken from \a geometry3D. + * + * \param top if \a true, create plane at top, otherwise at bottom + * (for PlaneOrientation Axial, for other plane locations respectively) + */ + virtual void InitializeStandardPlane( const Geometry3D* geometry3D, bool top, + PlaneOrientation planeorientation = Axial, + bool frontside=true, bool rotated=false ); + + /** + * \brief Initialize a plane with orientation \a planeorientation + * (default: axial) with respect to \a transform (default: identity) + * given width and height in units. + * + */ + virtual void InitializeStandardPlane( ScalarType width, ScalarType height, + const AffineTransform3D* transform = NULL, + PlaneOrientation planeorientation = Axial, + ScalarType zPosition = 0, bool frontside=true, bool rotated=false ); + + /** + * \brief Initialize plane with orientation \a planeorientation + * (default: axial) given width, height and spacing. + * + */ + virtual void InitializeStandardPlane( ScalarType width, ScalarType height, + const Vector3D & spacing, PlaneOrientation planeorientation = Axial, + ScalarType zPosition = 0, bool frontside = true, bool rotated = false ); + + /** + * \brief Initialize plane by width and height in pixels, right-/down-vector + * (itk) to describe orientation in world-space (vectors will be normalized) + * and spacing (default: 1.0 mm in all directions). + * + * The vectors are normalized and multiplied by the respective spacing before + * they are set in the matrix. + */ + virtual void InitializeStandardPlane( ScalarType width, ScalarType height, + const Vector3D& rightVector, const Vector3D& downVector, + const Vector3D *spacing = NULL ); + + /** + * \brief Initialize plane by width and height in pixels, + * right-/down-vector (vnl) to describe orientation in world-space (vectors + * will be normalized) and spacing (default: 1.0 mm in all directions). + * + * The vectors are normalized and multiplied by the respective spacing + * before they are set in the matrix. + */ + virtual void InitializeStandardPlane( ScalarType width, ScalarType height, + const VnlVector& rightVector, const VnlVector& downVector, + const Vector3D * spacing = NULL ); + + /** + * \brief Initialize plane by right-/down-vector (itk) and spacing + * (default: 1.0 mm in all directions). + * + * The length of the right-/-down-vector is used as width/height in units, + * respectively. Then, the vectors are normalized and multiplied by the + * respective spacing before they are set in the matrix. + */ + virtual void InitializeStandardPlane( const Vector3D& rightVector, + const Vector3D& downVector, const Vector3D * spacing = NULL ); + + /** + * \brief Initialize plane by right-/down-vector (vnl) and spacing + * (default: 1.0 mm in all directions). + * + * The length of the right-/-down-vector is used as width/height in units, + * respectively. Then, the vectors are normalized and multiplied by the + * respective spacing before they are set in the matrix. + */ + virtual void InitializeStandardPlane( const VnlVector& rightVector, + const VnlVector& downVector, const Vector3D * spacing = NULL ); + + /** + * \brief Initialize plane by origin and normal (size is 1.0 mm in + * all directions, direction of right-/down-vector valid but + * undefined). + * + */ + virtual void InitializePlane( const Point3D& origin, const Vector3D& normal); + + /** + * \brief Initialize plane by right-/down-vector. + * + * \warning The vectors are set into the matrix as they are, + * \em without normalization! + */ + void SetMatrixByVectors( const VnlVector& rightVector, + const VnlVector& downVector, ScalarType thickness=1.0 ); + + /** + * \brief Change \a transform so that the third column of the + * transform-martix is perpendicular to the first two columns + * + */ + static void EnsurePerpendicularNormal( AffineTransform3D* transform ); + + /** + * \brief Normal of the plane + * + */ + Vector3D GetNormal() const; + + /** + * \brief Normal of the plane as VnlVector + * + */ + VnlVector GetNormalVnl() const; + + virtual ScalarType SignedDistance( const Point3D& pt3d_mm ) const; + + virtual bool IsAbove( const Point3D& pt3d_mm ) const; + + /** + * \brief Distance of the point from the plane + * (bounding-box \em not considered) + * + */ + ScalarType DistanceFromPlane( const Point3D& pt3d_mm ) const ; + + /** + * \brief Signed distance of the point from the plane + * (bounding-box \em not considered) + * + * > 0 : point is in the direction of the direction vector. + */ + inline ScalarType SignedDistanceFromPlane( const Point3D& pt3d_mm ) const + { + ScalarType len = GetNormalVnl().two_norm(); - /** - * \brief Distance of the plane from another plane - * (bounding-box \em not considered) - * - * Result is 0 if planes are not parallel. - */ - ScalarType DistanceFromPlane(const PlaneGeometry* plane) const - { - return fabs(SignedDistanceFromPlane(plane)); - } + if( len == 0 ) + return 0; + return (pt3d_mm-GetOrigin())*GetNormal() / len; + } - /** - * \brief Signed distance of the plane from another plane - * (bounding-box \em not considered) - * - * Result is 0 if planes are not parallel. - */ - inline ScalarType SignedDistanceFromPlane( const PlaneGeometry *plane ) const - { - if(IsParallel(plane)) + /** + * \brief Distance of the plane from another plane + * (bounding-box \em not considered) + * + * Result is 0 if planes are not parallel. + */ + ScalarType DistanceFromPlane(const PlaneGeometry* plane) const { - return SignedDistance(plane->GetOrigin()); + return fabs(SignedDistanceFromPlane(plane)); } - return 0; - } - - - /** - * \brief Calculate the intersecting line of two planes - * - * \return \a true planes are intersecting - * \return \a false planes do not intersect - */ - bool IntersectionLine( const PlaneGeometry *plane, Line3D &crossline ) const; - - - /** - * \brief Calculate two points where another plane intersects the border of this plane - * - * \return number of intersection points (0..2). First interection point (if existing) - * is returned in \a lineFrom, second in \a lineTo. - */ - unsigned int IntersectWithPlane2D(const PlaneGeometry *plane, - Point2D &lineFrom, Point2D &lineTo ) const ; - - - /** - * \brief Calculate the angle between two planes - * - * \return angle in radiants - */ - double Angle( const PlaneGeometry *plane ) const; - - - /** - * \brief Calculate the angle between the plane and a line - * - * \return angle in radiants - */ - double Angle( const Line3D &line ) const; - - - /** - * \brief Calculate intersection point between the plane and a line - * - * \param intersectionPoint intersection point - * \return \a true if \em unique intersection exists, i.e., if line - * is \em not on or parallel to the plane - */ - bool IntersectionPoint( const Line3D &line, - Point3D &intersectionPoint ) const; - - - /** - * \brief Calculate line parameter of intersection point between the - * plane and a line - * - * \param t parameter of line: intersection point is - * line.GetPoint()+t*line.GetDirection() - * \return \a true if \em unique intersection exists, i.e., if line - * is \em not on or parallel to the plane - */ - bool IntersectionPointParam( const Line3D &line, double &t ) const; - - - /** - * \brief Returns whether the plane is parallel to another plane - * - * @return true iff the normal vectors both point to the same or exactly oposit direction - */ - bool IsParallel( const PlaneGeometry *plane ) const; - - /** - * \brief Returns whether the point is on the plane - * (bounding-box \em not considered) - */ - bool IsOnPlane( const Point3D &point ) const; - - - /** - * \brief Returns whether the line is on the plane - * (bounding-box \em not considered) - */ - bool IsOnPlane( const Line3D &line ) const; - - - /** - * \brief Returns whether the plane is on the plane - * (bounding-box \em not considered) - * - * @return true iff the normal vector of the planes point to the same or the exactly oposit direction and - * the distance of the planes is < eps - * - */ - bool IsOnPlane( const PlaneGeometry *plane ) const; - - - /** - * \brief Returns the lot from the point to the plane - */ - Point3D ProjectPointOntoPlane( const Point3D &pt ) const; - - - virtual void SetIndexToWorldTransform( AffineTransform3D *transform); - - - virtual void SetBounds( const BoundingBox::BoundsArrayType &bounds ); - - - virtual itk::LightObject::Pointer InternalClone() const; - - - /** Implements operation to re-orient the plane */ - virtual void ExecuteOperation( Operation *operation ); - - -protected: - PlaneGeometry(); - - virtual ~PlaneGeometry(); - - virtual void PrintSelf( std::ostream &os, itk::Indent indent ) const; - -private: - /** - * \brief Compares plane with another plane: \a true if IsOnPlane - * (bounding-box \em not considered) - */ - virtual bool operator==( const PlaneGeometry * ) const { return false; }; - - /** - * \brief Compares plane with another plane: \a false if IsOnPlane - * (bounding-box \em not considered) - */ - virtual bool operator!=( const PlaneGeometry * ) const { return false; }; - -}; + /** + * \brief Signed distance of the plane from another plane + * (bounding-box \em not considered) + * + * Result is 0 if planes are not parallel. + */ + inline ScalarType SignedDistanceFromPlane( const PlaneGeometry *plane ) const + { + if(IsParallel(plane)) + { + return SignedDistance(plane->GetOrigin()); + } + return 0; + } + /** + * \brief Calculate the intersecting line of two planes + * + * \return \a true planes are intersecting + * \return \a false planes do not intersect + */ + bool IntersectionLine( const PlaneGeometry *plane, Line3D &crossline ) const; + + /** + * \brief Calculate two points where another plane intersects the border of this plane + * + * \return number of intersection points (0..2). First interection point (if existing) + * is returned in \a lineFrom, second in \a lineTo. + */ + unsigned int IntersectWithPlane2D(const PlaneGeometry *plane, + Point2D &lineFrom, Point2D &lineTo ) const ; + + /** + * \brief Calculate the angle between two planes + * + * \return angle in radiants + */ + double Angle( const PlaneGeometry *plane ) const; + + /** + * \brief Calculate the angle between the plane and a line + * + * \return angle in radiants + */ + double Angle( const Line3D &line ) const; + + /** + * \brief Calculate intersection point between the plane and a line + * + * \param intersectionPoint intersection point + * \return \a true if \em unique intersection exists, i.e., if line + * is \em not on or parallel to the plane + */ + bool IntersectionPoint( const Line3D &line, + Point3D &intersectionPoint ) const; + + /** + * \brief Calculate line parameter of intersection point between the + * plane and a line + * + * \param t parameter of line: intersection point is + * line.GetPoint()+t*line.GetDirection() + * \return \a true if \em unique intersection exists, i.e., if line + * is \em not on or parallel to the plane + */ + bool IntersectionPointParam( const Line3D &line, double &t ) const; + + /** + * \brief Returns whether the plane is parallel to another plane + * + * @return true iff the normal vectors both point to the same or exactly oposit direction + */ + bool IsParallel( const PlaneGeometry *plane ) const; + + /** + * \brief Returns whether the point is on the plane + * (bounding-box \em not considered) + */ + bool IsOnPlane( const Point3D &point ) const; + + /** + * \brief Returns whether the line is on the plane + * (bounding-box \em not considered) + */ + bool IsOnPlane( const Line3D &line ) const; + + /** + * \brief Returns whether the plane is on the plane + * (bounding-box \em not considered) + * + * @return true iff the normal vector of the planes point to the same or the exactly oposit direction and + * the distance of the planes is < eps + * + */ + bool IsOnPlane( const PlaneGeometry *plane ) const; + + /** + * \brief Returns the lot from the point to the plane + */ + Point3D ProjectPointOntoPlane( const Point3D &pt ) const; + + virtual void SetIndexToWorldTransform( AffineTransform3D *transform); + + virtual itk::LightObject::Pointer InternalClone() const; + + /** Implements operation to re-orient the plane */ + virtual void ExecuteOperation( Operation *operation ); + + protected: + PlaneGeometry(); + + virtual ~PlaneGeometry(); + + virtual void PrintSelf( std::ostream &os, itk::Indent indent ) const; + + virtual void InternPreSetBounds( const BoundingBox::BoundsArrayType &bounds ); + + private: + /** + * \brief Compares plane with another plane: \a true if IsOnPlane + * (bounding-box \em not considered) + */ + virtual bool operator==( const PlaneGeometry * ) const { return false; }; + + /** + * \brief Compares plane with another plane: \a false if IsOnPlane + * (bounding-box \em not considered) + */ + virtual bool operator!=( const PlaneGeometry * ) const { return false; }; + }; } // namespace mitk - #endif /* PLANEGEOMETRY_H_HEADER_INCLUDED_C1C68A2C */