diff --git a/Modules/Core/include/mitkPlaneGeometry.h b/Modules/Core/include/mitkPlaneGeometry.h index d622cf45ff..4c853e7ff7 100644 --- a/Modules/Core/include/mitkPlaneGeometry.h +++ b/Modules/Core/include/mitkPlaneGeometry.h @@ -1,585 +1,595 @@ /*=================================================================== 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. ===================================================================*/ /** * \brief Describes the geometry of a plane object * * Describes a two-dimensional manifold, i.e., to put it simply, * an object that can be described using a 2D coordinate-system. * * PlaneGeometry 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 PlaneGeometry describes the 2D representation within a 3D object (derived from BaseGeometry). 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. * * * Optionally, a reference BaseGeometry 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 PlaneGeometry) with the bounding box of the associated * BaseGeometry. * * \warning The PlaneGeometry 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::GetPlaneGeometry() makes sure, that the PlaneGeometry 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 */ #ifndef PLANEGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #define PLANEGEOMETRY_H_HEADER_INCLUDED_C1C68A2C #include #include "mitkBaseGeometry.h" #include "mitkRestorePlanePositionOperation.h" #include namespace mitk { template < class TCoordRep, unsigned int NPointDimension > class Line; typedef Line Line3D; class PlaneGeometry; /** \deprecatedSince{2014_10} This class is deprecated. Please use PlaneGeometry instead. */ DEPRECATED(typedef PlaneGeometry Geometry2D); /** * \brief Describes a two-dimensional, rectangular plane * * \ingroup Geometry */ class MITKCORE_EXPORT PlaneGeometry : public BaseGeometry { public: mitkClassMacro(PlaneGeometry, BaseGeometry); /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) enum PlaneOrientation { Axial, Sagittal, Frontal, // also known as "Coronal" in mitk. None //This defines the PlaneGeometry for the 3D renderWindow which //curiously also needs a PlaneGeometry. This should be reconsidered some time. }; 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; /** * \brief Initialize a plane with orientation \a planeorientation * (default: axial) with respect to \a BaseGeometry (default: identity). * Spacing also taken from \a BaseGeometry. * * \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 BaseGeometry* 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 BaseGeometry (default: identity). * Spacing also taken from \a BaseGeometry. * * \param top if \a true, create plane at top, otherwise at bottom * (for PlaneOrientation Axial, for other plane locations respectively) */ virtual void InitializeStandardPlane(const BaseGeometry* 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. * * \a Rotated means rotated by 180 degrees (1/2 rotation) within the plane. * Rotation by 90 degrees (1/4 rotation) is not implemented as of now. * * \a Frontside/Backside: * Viewed from below = frontside in the axial case; * (radiologist's view versus neuro-surgeon's view, see: * http://www.itk.org/Wiki/images/e/ed/DICOM-OrientationDiagram-Radiologist-vs-NeuroSurgeon.png ) * Viewed from front = frontside in the coronal case; * Viewed from left = frontside in the sagittal case. * * \a Cave/Caution: Currently only RPI, LAI, LPS and RAS in the three standard planes are covered, * i.e. 12 cases of 144: 3 standard planes * 48 coordinate orientations = 144 cases. */ virtual void InitializeStandardPlane(ScalarType width, ScalarType height, const AffineTransform3D* transform = nullptr, 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. + * + * This overloaded version of InitializeStandardPlane() creates only righthanded + * coordinate orientations, unless spacing contains 1 or 3 negative entries. + * */ virtual void InitializeStandardPlane(ScalarType width, ScalarType height, const Vector3D& rightVector, const Vector3D& downVector, const Vector3D *spacing = nullptr); /** * \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. + * + * This overloaded version of InitializeStandardPlane() creates only righthanded + * coordinate orientations, unless spacing contains 1 or 3 negative entries. + * */ virtual void InitializeStandardPlane(ScalarType width, ScalarType height, const VnlVector& rightVector, const VnlVector& downVector, const Vector3D * spacing = nullptr); /** * \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 = nullptr); /** * \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 = nullptr); /** * \brief Initialize plane by origin and normal (size is 1.0 mm in * all directions, direction of right-/down-vector valid but * undefined). - * + * \warning This function can only produce righthanded coordinate orientation, not lefthanded. */ 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! + * This function creates a righthanded IndexToWorldTransform, + * only a negative thickness could still make it lefthanded. */ 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; /** * \brief Calculates, whether a point is below or above the plane. There are two different *calculation methods, with or without consideration of the bounding box. */ virtual bool IsAbove(const Point3D& pt3d_mm, bool considerBoundingBox = false) 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; } /** * \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)); } /** * \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 itk::LightObject::Pointer InternalClone() const override; /** Implements operation to re-orient the plane */ virtual void ExecuteOperation(Operation *operation) override; /** * \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 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 Set the geometrical frame of reference in which this PlaneGeometry * is placed. * * This would usually be the BaseGeometry of the underlying dataset, but * setting it is optional. */ void SetReferenceGeometry(mitk::BaseGeometry *geometry); /** * \brief Get the geometrical frame of reference for this PlaneGeometry. */ BaseGeometry *GetReferenceGeometry() const; bool HasReferenceGeometry() const; protected: PlaneGeometry(); PlaneGeometry(const PlaneGeometry& other); virtual ~PlaneGeometry(); virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; mitk::BaseGeometry *m_ReferenceGeometry; //##Documentation //## @brief PreSetSpacing //## //## These virtual function allows a different beahiour in subclasses. //## Do implement them in every subclass of BaseGeometry. If not needed, use //## {Superclass::PreSetSpacing();}; virtual void PreSetSpacing(const mitk::Vector3D& aSpacing) override{ Superclass::PreSetSpacing(aSpacing); }; //##Documentation //## @brief CheckBounds //## //## This function is called in SetBounds. Assertions can be implemented in this function (see PlaneGeometry.cpp). //## If you implement this function in a subclass, make sure, that all classes were your class inherits from //## have an implementation of CheckBounds //## (e.g. inheritance BaseGeometry <- A <- B. Implementation of CheckBounds in class B needs implementation in A as well!) virtual void CheckBounds(const BoundsArrayType& bounds) override; //##Documentation //## @brief CheckIndexToWorldTransform //## //## This function is called in SetIndexToWorldTransform. Assertions can be implemented in this function (see PlaneGeometry.cpp). //## In Subclasses of BaseGeometry, implement own conditions or call Superclass::CheckBounds(bounds);. virtual void CheckIndexToWorldTransform(mitk::AffineTransform3D* transform) override; 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 */ diff --git a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp index c85be4ab72..c94b8ff526 100644 --- a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp +++ b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp @@ -1,916 +1,985 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPlaneGeometry.h" #include "mitkPlaneOperation.h" #include "mitkInteractionConst.h" #include "mitkLine.h" #include #include namespace mitk { PlaneGeometry::PlaneGeometry() : Superclass(), m_ReferenceGeometry( nullptr ) { Initialize(); } PlaneGeometry::~PlaneGeometry() { } - PlaneGeometry::PlaneGeometry(const PlaneGeometry& other) - : Superclass(other), - m_ReferenceGeometry( other.m_ReferenceGeometry ) + PlaneGeometry::PlaneGeometry( const PlaneGeometry& other ) + : Superclass( other ), m_ReferenceGeometry( other.m_ReferenceGeometry ) { } - void - PlaneGeometry::EnsurePerpendicularNormal(mitk::AffineTransform3D *transform) + void PlaneGeometry::EnsurePerpendicularNormal( mitk::AffineTransform3D *transform ) { - //ensure row(2) of transform to be perpendicular to plane, keep length. + /** \brief ensure column(2) of indexToWorldTransform-matrix to be perpendicular to plane, keep length and handedness. */ + VnlVector normal = vnl_cross_3d( transform->GetMatrix().GetVnlMatrix().get_column(0), transform->GetMatrix().GetVnlMatrix().get_column(1) ); + normal.normalize(); // Now normal is a righthand normal unit vector, perpendicular to the plane. - normal.normalize(); ScalarType len = transform->GetMatrix().GetVnlMatrix().get_column(2).two_norm(); - - if (len==0){ len = 1; } + if (len==0) { len = 1; } normal *= len; - Matrix3D matrix = transform->GetMatrix(); - matrix.GetVnlMatrix().set_column(2, normal); - transform->SetMatrix( matrix ); + + // Get the existing normal vector zed: + vnl_vector_fixed zed = transform->GetMatrix().GetVnlMatrix().get_column(2); + + /** If det(matrix)<0, multiply normal vector by (-1) to keep geometry lefthanded. */ + if( vnl_determinant( transform->GetMatrix().GetVnlMatrix()) < 0 ) + { + MITK_DEBUG << "EnsurePerpendicularNormal(): Lefthanded geometry preserved, rh-normal: [ " << normal << " ],"; + normal *= (-1.0); + MITK_DEBUG << "lh-normal: [ " << normal << " ], original vector zed is: [ " << zed << " ]"; + } + + + // Now lets compare and only replace if necessary and only then warn the user: + + // float epsilon is precise enough here, but we need to respect numerical condition: + // Higham, N., 2002, Accuracy and Stability of Numerical Algorithms, + // SIAM, page 37, 2nd edition: + double feps = std::numeric_limits::epsilon(); + double zedsMagnitude = zed.two_norm(); + feps = feps * zedsMagnitude * 2; + + /** Check if normal (3. column) was perpendicular: If not, replace with calculated normal vector: */ + if( normal != zed ) + { + vnl_vector_fixed parallel; + for ( unsigned int i = 0; i<3 ; ++i ) + { + parallel[i] = normal[i] / zed[i]; // Remember linear algebra: checking for parallelity. + } + // Checking if really not paralell i.e. non-colinear vectors by comparing these floating point numbers: + if( (parallel[0]+feps < parallel[1] || feps+parallel[1] < parallel[0]) + && (parallel[0]+feps < parallel[2] || feps+parallel[2] < parallel[0]) ) + { + MITK_WARN << "EnsurePerpendicularNormal(): Plane geometry was _/askew/_, so here it gets rectified by substituting" + << " the 3rd column of the indexToWorldMatrix with an appropriate normal vector: [ " << normal + << " ], original vector zed was: [ " << zed << " ]."; + + Matrix3D matrix = transform->GetMatrix(); + matrix.GetVnlMatrix().set_column( 2, normal ); + transform->SetMatrix( matrix ); + } + } + else + { + // Nothing to do, 3rd column of indexToWorldTransformMatrix already was perfectly perpendicular. + } } void PlaneGeometry::CheckIndexToWorldTransform(mitk::AffineTransform3D* transform) { EnsurePerpendicularNormal(transform); } void PlaneGeometry::CheckBounds(const BoundingBox::BoundsArrayType &bounds) { // error: unused parameter 'bounds' // this happens in release mode, where the assert macro is defined empty // hence we "use" the parameter: (void)bounds; //currently the unit rectangle must be starting at the origin [0,0] assert(bounds[0]==0); assert(bounds[2]==0); //the unit rectangle must be two-dimensional assert(bounds[1]>0); assert(bounds[3]>0); } void PlaneGeometry::IndexToWorld( const Point2D &pt_units, Point2D &pt_mm ) const { pt_mm[0] = GetExtentInMM(0) / GetExtent(0)*pt_units[0]; pt_mm[1] = GetExtentInMM(1) / GetExtent(1)*pt_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D &pt_mm, Point2D &pt_units ) const { pt_units[0] = pt_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0))); pt_units[1] = pt_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(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::IndexToWorld(const Vector2D &vec_units, Vector2D &vec_mm) const { vec_mm[0] = (GetExtentInMM(0) / GetExtent(0)) * vec_units[0]; vec_mm[1] = (GetExtentInMM(1) / GetExtent(1)) * vec_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D & /*atPt2d_mm*/, const Vector2D &vec_mm, Vector2D &vec_units) const { MITK_WARN<<"Warning! Call of the deprecated function PlaneGeometry::WorldToIndex(point, vec, vec). Use PlaneGeometry::WorldToIndex(vec, vec) instead!"; this->WorldToIndex(vec_mm, vec_units); } void PlaneGeometry::WorldToIndex( const Vector2D &vec_mm, Vector2D &vec_units) const { vec_units[0] = vec_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0))); vec_units[1] = vec_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(1))); } 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, mitk::ScalarType height, const AffineTransform3D* transform /* = nullptr */, PlaneGeometry::PlaneOrientation planeorientation /* = Axial */, mitk::ScalarType zPosition /* = 0 */, bool frontside /* = true */, bool rotated /* = false */ ) { Superclass::Initialize(); - //construct standard view + /// construct standard view. - // We define at the moment "frontside" as: axial from above, coronal from front (nose), saggital from right. - // TODO: Double check with medics [ ]. + // We define at the moment "frontside" as: axial from above, + // coronal from front (nose), saggital from right. + // TODO: Double check with medicals doctors or radiologists [ ]. - // We define the orientation in patient's view, e.g. LAI is in a axial cut (parallel to the triangle ear-ear-nose): + // We define the orientation in patient's view, e.g. LAI is in a axial cut + // (parallel to the triangle ear-ear-nose): // first axis: To the left ear of the patient // seecond axis: To the nose of the patient // third axis: To the legs of the patient. - // Options are: L/R left/right; A/P anterior/posterior; I/S inferior/superior (AKA caudal/cranial). - // We note on all cases in the following switch block r.h. for right handed or l.h. for left handed to describe the - // different cases. However, which system is chosen is defined at the end of the switch block. + // Options are: L/R left/right; A/P anterior/posterior; I/S inferior/superior + // (AKA caudal/cranial). + // We note on all cases in the following switch block r.h. for right handed + // or l.h. for left handed to describe the different cases. + // However, which system is chosen is defined at the end of the switch block. - // Cave/careful: the vectors right and bottom are relative to the plane + // CAVE / be careful: the vectors right and bottom are relative to the plane // and do NOT describe e.g. the right side of the patient. Point3D origin; - VnlVector rightDV(3), bottomDV(3); // Bottom means downwards, DV means Direction Vector. Both relative to the image! - origin.Fill(0); /// Origin of this plane is by default a zero vector and implicitly in the top-left corner: - /// This is different to all definitions in MITK, except the QT mouse clicks. - // But it is like this here and we don't want to change a running system. - // Just be aware, that IN THIS FUNCTION we define the origin at the top left (e.g. your screen). - - int normalDirection; // NormalDirection defines which axis (i.e. column index in the transform matrix) - // is perpendicular to the plane. + /** Bottom means downwards, DV means Direction Vector. Both relative to the image! */ + VnlVector rightDV(3), bottomDV(3); + /** Origin of this plane is by default a zero vector and implicitly in the top-left corner: */ + origin.Fill(0); + /** This is different to all definitions in MITK, except the QT mouse clicks. + * But it is like this here and we don't want to change a running system. + * Just be aware, that IN THIS FUNCTION we define the origin at the top left (e.g. your screen). */ + + /** NormalDirection defines which axis (i.e. column index in the transform matrix) + * is perpendicular to the plane: */ + int normalDirection; switch(planeorientation) // Switch through our limited choice of standard planes: { - case None: /// Orientation 'None' shall be done like the axial plane orientation, for whatever reasons. + case None: + /** Orientation 'None' shall be done like the axial plane orientation, + * for whatever reasons. */ case Axial: if(frontside) // Radiologist's view from below. A cut along the triangle ear-ear-nose. { - if(rotated==false) /// Origin in the top-left corner, x=[1; 0; 0], y=[0; 1; 0], z=[0; 0; 1], origin=[0,0,zpos]: LAI (r.h.) - { /** 0--right--> - * | Picture of a finite, rectangular plane... - * | (insert LOL/CAT-scan here). - * v */ + if(rotated==false) + /** Origin in the top-left corner, x=[1; 0; 0], y=[0; 1; 0], z=[0; 0; 1], + * origin=[0,0,zpos]: LAI (r.h.) + * + * 0---rightDV----> | + * | | + * | Picture of a finite, rectangular plane | + * | ( insert LOLCAT-scan here ^_^ ) | + * | | + * v _________________________________________| + * + */ + { FillVector3D(origin, 0, 0, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } - else // Origin rotated to the bottom-right corner, x=[-1; 0; 0], y=[0; -1; 0], z=[0; 0; 1], origin=[w,h,zpos]: RPI (r.h.) + else // Origin rotated to the bottom-right corner, x=[-1; 0; 0], y=[0; -1; 0], z=[0; 0; 1], + // origin=[w,h,zpos]: RPI (r.h.) { // Caveat emptor: Still using top-left as origin of index coordinate system! FillVector3D(origin, width, height, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } else // 'Backside, not frontside.' Neuro-Surgeons's view from above patient. { if(rotated==false) // x=[-1; 0; 0], y=[0; 1; 0], z=[0; 0; 1], origin=[w,0,zpos]: RAS (r.h.) { FillVector3D(origin, width, 0, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } - else // Origin in the bottom-left corner, x=[1; 0; 0], y=[0; -1; 0], z=[0; 0; 1], origin=[0,h,zpos]: LPS (r.h.) + else // Origin in the bottom-left corner, x=[1; 0; 0], y=[0; -1; 0], z=[0; 0; 1], + // origin=[0,h,zpos]: LPS (r.h.) { FillVector3D(origin, 0, height, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } normalDirection = 2; // That is S=Superior=z=third_axis=middlefinger in righthanded LPS-system. break; - case Frontal: // Frontal is known as Coronal in mitk. Plane cuts through patient's ear-ear-heel-heel + // Frontal is known as Coronal in mitk. Plane cuts through patient's ear-ear-heel-heel: + case Frontal: if(frontside) { if(rotated==false) // x=[1; 0; 0], y=[0; 0; 1], z=[0; 1; 0], origin=[0,zpos,0]: LAI (r.h.) { FillVector3D(origin, 0, zPosition, 0); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else // x=[-1;0;0], y=[0;0;-1], z=[0;1;0], origin=[w,zpos,h]: RAS (r.h.) { FillVector3D(origin, width, zPosition, height); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) // x=[-1;0;0], y=[0;0;1], z=[0;1;0], origin=[w,zpos,0]: RPI (r.h.) { FillVector3D(origin, width, zPosition, 0); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else // x=[1;0;0], y=[0;1;0], z=[0;0;-1], origin=[0,zpos,h]: LPS (r.h.) { FillVector3D(origin, 0, zPosition, height); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 1; // Normal vector = posterior direction. break; case Sagittal: // Sagittal=Medial plane, the symmetry-plane mirroring your face. if(frontside) { if(rotated==false) // x=[0;1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,0,0]: LAI (r.h.) { FillVector3D(origin, zPosition, 0, 0); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, 1); } else // x=[0;-1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,w,h]: LPS (r.h.) { FillVector3D(origin, zPosition, width, height); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) // x=[0;-1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,w,0]: RPI (r.h.) { FillVector3D(origin, zPosition, width, 0); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, 1); } else // x=[0;1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,0,h]: RAS (r.h.) { FillVector3D(origin, zPosition, 0, height); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 0; // Normal vector = Lateral direction: Left in a LPS-system. break; - // <2015-10-23, Esther Wild and Martin Hettich:> We found a problem: - // Till now this function only covers 4 cases: RPI, LAI, LPS, RAS; - // And we scroll the standard planes along each axis in these cases. --> 3*4 = 12. - // So far only 4 of 48 possible coordinate orientations and - // therefore 12 of 144 cases of standard planes are treated here. - // - // TODO: To handle r.h./l.h. systems correctly: check and correct normal vector and BoundingBox. - // TODO: Also see bugs.mitk.org #11477, #18662, #19347. - default: itkExceptionMacro("unknown PlaneOrientation"); } + /// Checking if lefthanded or righthanded: + mitk::ScalarType lhOrRhSign=(+1); if ( transform != nullptr ) { + lhOrRhSign = ( (0 < vnl_determinant( transform->GetMatrix().GetVnlMatrix() )) ? (+1) : (-1) ); + + MITK_DEBUG << "mitk::PlaneGeometry::InitializeStandardPlane(): lhOrRhSign, normalDirection, NameOfClass, ObjectName = " + << lhOrRhSign << ", " << normalDirection << ", " << transform->GetNameOfClass() + << ", " << transform->GetObjectName() << "."; + origin = transform->TransformPoint( origin ); rightDV = transform->TransformVector( rightDV ); bottomDV = transform->TransformVector( bottomDV ); - } - - ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; // TODO: check signs according to l.h./r.h. system. - this->SetBounds( bounds ); - if ( transform == nullptr ) - { - this->SetMatrixByVectors( rightDV, bottomDV ); // TODO: Spacing = 1 or -1 ? + /// Signing normal vector according to l.h./r.h. coordinate system: + this->SetMatrixByVectors( rightDV, + bottomDV, + transform->GetMatrix().GetVnlMatrix().get_column(normalDirection).two_norm() * lhOrRhSign ); } - else + else if ( transform == nullptr ) { - this->SetMatrixByVectors( rightDV, bottomDV, - transform->GetMatrix().GetVnlMatrix().get_column(normalDirection).magnitude() ); - // TODO: check signs according to l.h./r.h. system. + this->SetMatrixByVectors( rightDV, bottomDV ); } - this->SetOrigin(origin); + ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; + + this->SetBounds( bounds ); + + this->SetOrigin( origin ); } + void PlaneGeometry::InitializeStandardPlane( const BaseGeometry *geometry3D, PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { this->SetReferenceGeometry( const_cast< BaseGeometry * >( geometry3D ) ); ScalarType width, height; const BoundingBox::BoundsArrayType& boundsarray = geometry3D->GetBoundingBox()->GetBounds(); Vector3D originVector; FillVector3D(originVector, boundsarray[0], boundsarray[2], boundsarray[4]); if(geometry3D->GetImageGeometry()) { FillVector3D( originVector, originVector[0] - 0.5, originVector[1] - 0.5, originVector[2] - 0.5 ); } switch(planeorientation) { case None: case Axial: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(1); break; case Frontal: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(2); break; case Sagittal: width = geometry3D->GetExtent(1); height = geometry3D->GetExtent(2); break; default: itkExceptionMacro("unknown PlaneOrientation"); } InitializeStandardPlane( width, height, geometry3D->GetIndexToWorldTransform(), planeorientation, zPosition, frontside, rotated ); ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); Point3D origin; originVector = geometry3D->GetIndexToWorldTransform() ->TransformVector( originVector ); origin = GetOrigin() + originVector; SetOrigin(origin); } void PlaneGeometry::InitializeStandardPlane( const BaseGeometry *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated ) { ScalarType zPosition; switch(planeorientation) { case Axial: zPosition = (top ? 0.5 : geometry3D->GetExtent(2)-0.5); break; case Frontal: zPosition = (top ? 0.5 : geometry3D->GetExtent(1)-0.5); break; case Sagittal: zPosition = (top ? 0.5 : geometry3D->GetExtent(0)-0.5); break; case None: zPosition = (top ? 0 : geometry3D->GetExtent(2)-1.0); 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(); + ScalarType width = rightVector.two_norm(); + ScalarType height = downVector.two_norm(); 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(); + // Crossproduct vnl_cross_3d is always righthanded, but that is okay here + // because in this method we create a new IndexToWorldTransform and + // spacing with 1 or 3 negative components could still make it lefthanded. if(spacing!=nullptr) { rightDV *= (*spacing)[0]; downDV *= (*spacing)[1]; normal *= (*spacing)[2]; } AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightDV); matrix.GetVnlMatrix().set_column(1, downDV); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset()); ScalarType bounds[6] = { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); this->SetIndexToWorldTransform( transform ); } void PlaneGeometry::InitializePlane( const Point3D &origin, const Vector3D &normal ) { VnlVector rightVectorVnl(3), downVectorVnl; if( Equal( normal[1], 0.0f ) == false ) { FillVector3D( rightVectorVnl, 1.0f, -normal[0]/normal[1], 0.0f ); rightVectorVnl.normalize(); } else { FillVector3D( rightVectorVnl, 0.0f, 1.0f, 0.0f ); } downVectorVnl = vnl_cross_3d( normal.GetVnlVector(), rightVectorVnl ); downVectorVnl.normalize(); + // Crossproduct vnl_cross_3d is always righthanded. InitializeStandardPlane( rightVectorVnl, downVectorVnl ); SetOrigin(origin); } void PlaneGeometry::SetMatrixByVectors( const VnlVector &rightVector, const VnlVector &downVector, ScalarType thickness /* = 1.0 */ ) { VnlVector normal = vnl_cross_3d(rightVector, downVector); normal.normalize(); normal *= thickness; + // Crossproduct vnl_cross_3d is always righthanded, but that is okay here + // because in this method we create a new IndexToWorldTransform and + // a negative thickness could still make it lefthanded. AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightVector); matrix.GetVnlMatrix().set_column(1, downVector); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); transform->SetOffset(this->GetIndexToWorldTransform()->GetOffset()); SetIndexToWorldTransform(transform); } Vector3D PlaneGeometry::GetNormal() const { Vector3D frontToBack; frontToBack.SetVnlVector( this->GetIndexToWorldTransform() ->GetMatrix().GetVnlMatrix().get_column(2) ); return frontToBack; } VnlVector PlaneGeometry::GetNormalVnl() const { return this->GetIndexToWorldTransform() ->GetMatrix().GetVnlMatrix().get_column(2); } ScalarType PlaneGeometry::DistanceFromPlane( const Point3D &pt3d_mm ) const { return fabs(SignedDistance( pt3d_mm )); } ScalarType PlaneGeometry::SignedDistance( const Point3D &pt3d_mm ) const { return SignedDistanceFromPlane(pt3d_mm); } bool PlaneGeometry::IsAbove( const Point3D &pt3d_mm , bool considerBoundingBox) const { if(considerBoundingBox) { Point3D pt3d_units; BaseGeometry::WorldToIndex(pt3d_mm, pt3d_units); return (pt3d_units[2] > this->GetBoundingBox()->GetBounds()[4]); } else return SignedDistanceFromPlane(pt3d_mm) > 0; } bool PlaneGeometry::IntersectionLine( const PlaneGeometry* plane, Line3D& crossline ) const { Vector3D normal = this->GetNormal(); normal.Normalize(); Vector3D planeNormal = plane->GetNormal(); planeNormal.Normalize(); Vector3D direction = itk::CrossProduct( normal, planeNormal ); if ( direction.GetSquaredNorm() < eps ) return false; crossline.SetDirection( direction ); double N1dN2 = normal * planeNormal; double determinant = 1.0 - N1dN2 * N1dN2; Vector3D origin = this->GetOrigin().GetVectorFromOrigin(); Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin(); double d1 = normal * origin; double d2 = planeNormal * planeOrigin; double c1 = ( d1 - d2 * N1dN2 ) / determinant; double c2 = ( d2 - d1 * N1dN2 ) / determinant; Vector3D p = normal * c1 + planeNormal * c2; crossline.GetPoint().GetVnlVector() = p.GetVnlVector(); return true; } unsigned int PlaneGeometry::IntersectWithPlane2D( const PlaneGeometry* plane, Point2D& lineFrom, Point2D &lineTo ) const { Line3D crossline; if ( this->IntersectionLine( plane, crossline ) == false ) return 0; Point2D point2; Vector2D direction2; this->Map( crossline.GetPoint(), point2 ); this->Map( crossline.GetPoint(), crossline.GetDirection(), direction2 ); return Line3D::RectangleLineIntersection( 0, 0, GetExtentInMM(0), GetExtentInMM(1), point2, direction2, lineFrom, lineTo ); } double PlaneGeometry::Angle( const PlaneGeometry *plane ) const { return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2)); } double PlaneGeometry::Angle( const Line3D &line ) const { return vnl_math::pi_over_2 - angle( line.GetDirection().GetVnlVector(), GetMatrixColumn(2) ); } bool PlaneGeometry::IntersectionPoint( const Line3D &line, Point3D &intersectionPoint ) const { Vector3D planeNormal = this->GetNormal(); planeNormal.Normalize(); Vector3D lineDirection = line.GetDirection(); lineDirection.Normalize(); double t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; intersectionPoint = line.GetPoint() + lineDirection * t; return true; } bool PlaneGeometry::IntersectionPointParam( const Line3D &line, double &t ) const { Vector3D planeNormal = this->GetNormal(); Vector3D lineDirection = line.GetDirection(); t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; return true; } bool PlaneGeometry::IsParallel( const PlaneGeometry *plane ) const { return ( (Angle(plane) < 10.0 * mitk::sqrteps ) || ( Angle(plane) > ( vnl_math::pi - 10.0 * sqrteps ) ) ) ; } bool PlaneGeometry::IsOnPlane( const Point3D &point ) const { return Distance(point) < eps; } bool PlaneGeometry::IsOnPlane( const Line3D &line ) const { return ( (Distance( line.GetPoint() ) < eps) && (Distance( line.GetPoint2() ) < eps) ); } bool PlaneGeometry::IsOnPlane( const PlaneGeometry *plane ) const { return ( IsParallel( plane ) && (Distance( plane->GetOrigin() ) < eps) ); } Point3D PlaneGeometry::ProjectPointOntoPlane( const Point3D& pt ) const { ScalarType len = this->GetNormalVnl().two_norm(); return pt - this->GetNormal() * this->SignedDistanceFromPlane( pt ) / len; } itk::LightObject::Pointer PlaneGeometry::InternalClone() const { Self::Pointer newGeometry = new PlaneGeometry(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void PlaneGeometry::ExecuteOperation( Operation *operation ) { vtkTransform *transform = vtkTransform::New(); transform->SetMatrix( this->GetVtkMatrix()); switch ( operation->GetOperationType() ) { case OpORIENT: { mitk::PlaneOperation *planeOp = dynamic_cast< mitk::PlaneOperation * >( operation ); if ( planeOp == nullptr ) { return; } Point3D center = planeOp->GetPoint(); Vector3D orientationVector = planeOp->GetNormal(); Vector3D defaultVector; FillVector3D( defaultVector, 0.0, 0.0, 1.0 ); Vector3D rotationAxis = itk::CrossProduct( orientationVector, defaultVector ); //double rotationAngle = acos( orientationVector[2] / orientationVector.GetNorm() ); double rotationAngle = atan2( (double) rotationAxis.GetNorm(), (double) (orientationVector * defaultVector) ); rotationAngle *= 180.0 / vnl_math::pi; transform->PostMultiply(); transform->Identity(); transform->Translate( center[0], center[1], center[2] ); transform->RotateWXYZ( rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2] ); transform->Translate( -center[0], -center[1], -center[2] ); break; } case OpRESTOREPLANEPOSITION: { RestorePlanePositionOperation *op = dynamic_cast< mitk::RestorePlanePositionOperation* >(operation); if(op == nullptr) { 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); this->Modified(); transform->Delete(); return; } default: Superclass::ExecuteOperation( operation ); transform->Delete(); return; } this->SetVtkMatrixDeepCopy(transform); this->Modified(); transform->Delete(); } void PlaneGeometry::PrintSelf( std::ostream& os, itk::Indent indent ) const { Superclass::PrintSelf(os,indent); os << indent << " ScaleFactorMMPerUnitX: " << GetExtentInMM(0) / GetExtent(0) << std::endl; os << indent << " ScaleFactorMMPerUnitY: " << GetExtentInMM(1) / GetExtent(1) << std::endl; os << indent << " Normal: " << GetNormal() << std::endl; } bool PlaneGeometry::Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const { assert(this->IsBoundingBoxNull()==false); Point3D pt3d_units; Superclass::WorldToIndex(pt3d_mm, pt3d_units); pt2d_mm[0] = pt3d_units[0] * GetExtentInMM(0) / GetExtent(0); pt2d_mm[1] = pt3d_units[1] * GetExtentInMM(1) / GetExtent(1); pt3d_units[2]=0; return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } void PlaneGeometry::Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const { //pt2d_mm is measured from the origin of the world geometry (at leats it called form BaseRendere::Mouse...Event) Point3D pt3d_units; pt3d_units[0] = pt2d_mm[0] / (GetExtentInMM(0) / GetExtent(0)); pt3d_units[1] = pt2d_mm[1] / (GetExtentInMM(1) / GetExtent(1)); pt3d_units[2]=0; //pt3d_units is a continuos index. We divided it with the Scale Factor (= spacing in x and y) to convert it from mm to index units. // pt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units); //now we convert the 3d index to a 3D world point in mm. We could have used IndexToWorld as well as GetITW->Transform... } void PlaneGeometry::SetSizeInUnits(mitk::ScalarType width, mitk::ScalarType height) { ScalarType bounds[6]={0, width, 0, height, 0, 1}; ScalarType extent, newextentInMM; if(GetExtent(0)>0) { extent = GetExtent(0); if(width>extent) newextentInMM = GetExtentInMM(0)/width*extent; else newextentInMM = GetExtentInMM(0)*extent/width; SetExtentInMM(0, newextentInMM); } if(GetExtent(1)>0) { extent = GetExtent(1); if(width>extent) newextentInMM = GetExtentInMM(1)/height*extent; else newextentInMM = GetExtentInMM(1)*extent/height; SetExtentInMM(1, newextentInMM); } SetBounds(bounds); } bool PlaneGeometry::Project( const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const { assert(this->IsBoundingBoxNull()==false); Point3D pt3d_units; Superclass::WorldToIndex(pt3d_mm, pt3d_units); pt3d_units[2] = 0; projectedPt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units); return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool PlaneGeometry::Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; Superclass::WorldToIndex(vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units); return true; } bool PlaneGeometry::Project(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const { MITK_WARN << "Deprecated function! Call Project(vec3D,vec3D) instead."; assert(this->IsBoundingBoxNull()==false); Vector3D vec3d_units; Superclass::WorldToIndex(atPt3d_mm, vec3d_mm, vec3d_units); vec3d_units[2] = 0; projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units); Point3D pt3d_units; Superclass::WorldToIndex(atPt3d_mm, pt3d_units); return const_cast(this->GetBoundingBox())->IsInside(pt3d_units); } bool PlaneGeometry::Map(const mitk::Point3D & atPt3d_mm, const mitk::Vector3D &vec3d_mm, mitk::Vector2D &vec2d_mm) const { Point2D pt2d_mm_start, pt2d_mm_end; Point3D pt3d_mm_end; bool inside=Map(atPt3d_mm, pt2d_mm_start); pt3d_mm_end = atPt3d_mm+vec3d_mm; inside&=Map(pt3d_mm_end, pt2d_mm_end); vec2d_mm=pt2d_mm_end-pt2d_mm_start; return inside; } void PlaneGeometry::Map(const mitk::Point2D &/*atPt2d_mm*/, const mitk::Vector2D &/*vec2d_mm*/, mitk::Vector3D &/*vec3d_mm*/) const { //@todo implement parallel to the other Map method! assert(false); } void PlaneGeometry::SetReferenceGeometry( mitk::BaseGeometry *geometry ) { m_ReferenceGeometry = geometry; } mitk::BaseGeometry * PlaneGeometry::GetReferenceGeometry() const { return m_ReferenceGeometry; } bool PlaneGeometry::HasReferenceGeometry() const { return ( m_ReferenceGeometry != nullptr ); } } // namespace diff --git a/Modules/Core/test/mitkPlaneGeometryTest.cpp b/Modules/Core/test/mitkPlaneGeometryTest.cpp index e467070d98..5d31893fc1 100644 --- a/Modules/Core/test/mitkPlaneGeometryTest.cpp +++ b/Modules/Core/test/mitkPlaneGeometryTest.cpp @@ -1,845 +1,944 @@ /*=================================================================== 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 "mitkRotationOperation.h" #include "mitkInteractionConst.h" #include "mitkLine.h" #include "mitkGeometry3D.h" #include "mitkThinPlateSplineCurvedGeometry.h" #include "mitkSlicedGeometry3D.h" +#include "mitkAffineTransform3D.h" +#include "mitkBaseGeometry.h" #include #include #include #include #include #include static const mitk::ScalarType testEps = 1E-9; // the epsilon used in this test == at least float precision. class mitkPlaneGeometryTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPlaneGeometryTestSuite); MITK_TEST(TestInitializeStandardPlane); MITK_TEST(TestProjectPointOntoPlane); MITK_TEST(TestPlaneGeometryCloning); MITK_TEST(TestInheritance); MITK_TEST(TestSetExtendInMM); MITK_TEST(TestRotate); MITK_TEST(TestClone); MITK_TEST(TestPlaneComparison); MITK_TEST(TestAxialInitialization); MITK_TEST(TestFrontalInitialization); MITK_TEST(TestSaggitalInitialization); + MITK_TEST(TestLefthandedCoordinateSystem); // Currently commented out, see See bug 15990 // MITK_TEST(testPlaneGeometryInitializeOrder); MITK_TEST(TestIntersectionPoint); MITK_TEST(TestCase1210); CPPUNIT_TEST_SUITE_END(); private: // private test members that are initialized by setUp() mitk::PlaneGeometry::Pointer planegeometry; mitk::Point3D origin; mitk::Vector3D right, bottom, normal; mitk::ScalarType width, height; mitk::ScalarType widthInMM, heightInMM, thicknessInMM; public: void setUp() override { planegeometry = mitk::PlaneGeometry::New(); width = 100; widthInMM = width; height = 200; heightInMM = height; thicknessInMM = 1.0; mitk::FillVector3D(origin, 4.5, 7.3, 11.2); mitk::FillVector3D(right, widthInMM, 0, 0); mitk::FillVector3D(bottom, 0, heightInMM, 0); mitk::FillVector3D(normal, 0, 0, thicknessInMM); planegeometry->InitializeStandardPlane(right.GetVnlVector(), bottom.GetVnlVector()); planegeometry->SetOrigin(origin); } void tearDown() override { } // This test verifies inheritance behaviour, this test will fail if the behaviour changes in the future void TestInheritance() { mitk::PlaneGeometry::Pointer plane = mitk::PlaneGeometry::New(); mitk::Geometry3D::Pointer g3d = dynamic_cast < mitk::Geometry3D* > ( plane.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("Planegeometry should not be castable to Geometry 3D", g3d.IsNull()); mitk::BaseGeometry::Pointer base = dynamic_cast < mitk::BaseGeometry* > ( plane.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("Planegeometry should be castable to BaseGeometry", base.IsNotNull()); base = nullptr; g3d = mitk::Geometry3D::New(); base = dynamic_cast < mitk::BaseGeometry* > ( g3d.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("Geometry3D should be castable to BaseGeometry", base.IsNotNull()); g3d=nullptr; mitk::SlicedGeometry3D::Pointer sliced = mitk::SlicedGeometry3D::New(); g3d = dynamic_cast < mitk::Geometry3D* > ( sliced.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("SlicedGeometry3D should not be castable to Geometry3D", g3d.IsNull()); plane=nullptr; mitk::ThinPlateSplineCurvedGeometry::Pointer thin = mitk::ThinPlateSplineCurvedGeometry::New(); plane = dynamic_cast < mitk::PlaneGeometry* > ( thin.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("AbstractTransformGeometry should be castable to PlaneGeometry", plane.IsNotNull()); plane = mitk::PlaneGeometry::New(); mitk::AbstractTransformGeometry::Pointer atg = dynamic_cast < mitk::AbstractTransformGeometry* > ( plane.GetPointer() ); CPPUNIT_ASSERT_MESSAGE("PlaneGeometry should not be castable to AbstractTransofrmGeometry", atg.IsNull()); } + + void TestLefthandedCoordinateSystem() + { + /** + * @brief This method tests InitializeStandardPlane() and IndexToWorld() + * with a left-handed coordinate orientation or indexToWorldMatrix. + * + * Of course this test does not use standard Parameters, which are right-handed. + * See also discussion of bug #11477: http://bugs.mitk.org/show_bug.cgi?id=11477 + */ + planegeometry = mitk::PlaneGeometry::New(); + width = 100; + widthInMM = 5; + height = 200; + heightInMM = 3; + thicknessInMM = 1.0; + mitk::FillVector3D(right, widthInMM, 0, 0); + mitk::FillVector3D(bottom, 0, heightInMM, 0); + // This one negative sign results in lefthanded coordinate orientation and det(matrix) < 0. + mitk::FillVector3D(normal, 0, 0, -thicknessInMM); + + mitk::AffineTransform3D::Pointer transform = mitk::AffineTransform3D::New(); + + mitk::AffineTransform3D::MatrixType matrix; + mitk::AffineTransform3D::MatrixType::InternalMatrixType& vnl_matrix = matrix.GetVnlMatrix(); + + vnl_matrix.set_column(0, right); + vnl_matrix.set_column(1, bottom); + vnl_matrix.set_column(2, normal); + + // making sure that we didn't screw up this special test case or else fail deadly: + assert( vnl_determinant(vnl_matrix) < 0.0 ); + + transform->SetIdentity(); + transform->SetMatrix(matrix); + + planegeometry->InitializeStandardPlane( width, height, transform ); // Crux of the matter. + CPPUNIT_ASSERT_MESSAGE("Testing if IndexToWorldMatrix is correct after InitializeStandardPlane( width, height, transform ) ", + mitk::MatrixEqualElementWise( planegeometry->GetIndexToWorldTransform()->GetMatrix(), matrix ) ); + + mitk::Point3D p_index; + p_index[0] = 10.; + p_index[1] = 10.; + p_index[2] = 0.; + mitk::Point3D p_world; + mitk::Point3D p_expectedResult; + p_expectedResult[0] = 50.; + p_expectedResult[1] = 30.; + p_expectedResult[2] = 0.; + + ((mitk::BaseGeometry::Pointer) planegeometry)->IndexToWorld(p_index, p_world); // Crux of the matter. + CPPUNIT_ASSERT_MESSAGE( "Testing if IndexToWorld(a,b) function works correctly with lefthanded matrix ", + mitk::Equal(p_world, p_expectedResult, testEps ) ); + } + + // See bug 1210 // Test does not use standard Parameters void TestCase1210() { mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); mitk::Point3D origin; mitk::Vector3D right, down, spacing; mitk::FillVector3D(origin, 4.5, 7.3, 11.2); mitk::FillVector3D(right, - 1.015625, 1.015625, 1.1999969482421875 - ); + 1.015625, 1.015625, 1.1999969482421875 + ); mitk::FillVector3D(down, - 1.4012984643248170709237295832899161312802619418765e-45, 0, 0 - ); + 1.4012984643248170709237295832899161312802619418765e-45, 0, 0 + ); mitk::FillVector3D(spacing, - 0, 1.4713633875410579244699160624544119378442750389703e-43, 9.2806360452222355258639080851310540729807238879469e-32 - ); + 0, 1.4713633875410579244699160624544119378442750389703e-43, + 9.2806360452222355258639080851310540729807238879469e-32 + ); std::cout << "Testing InitializeStandardPlane(rightVector, downVector, spacing = NULL): "<InitializeStandardPlane(right, down, &spacing)); /* std::cout << "Testing width, height and thickness (in units): "; if((mitk::Equal(planegeometry->GetExtent(0),width)==false) || (mitk::Equal(planegeometry->GetExtent(1),height)==false) || (mitk::Equal(planegeometry->GetExtent(2),1)==false) ) { std::cout<<"[FAILED]"<GetExtentInMM(0),widthInMM)==false) || (mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM)==false) || (mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM)==false) ) { std::cout<<"[FAILED]"< 0. - * - */ + * @brief This method tests method IntersectionPoint + * + * See also bug #7151. (ref 2 this test: iggy) + * This test was written due to incorrect calculation of the intersection point + * between a given line and plane. This only occured when the pointdistance of + * the line was less than 1. + * Test Behavour: + * ============== + * we have a given line and a given plane. + * we let the line intersect the plane. + * when testing several positions on the line the resulting intersection point must be the same + * we test a position where the distance between the correspoinding points is < 0 and another position where the distance is > 0. + * + */ // Test does not use standard Parameters void TestIntersectionPoint() { //init plane with its parameter mitk::PlaneGeometry::Pointer myPlaneGeometry = mitk::PlaneGeometry::New(); mitk::Point3D origin; origin[0] = 0.0; origin[1] = 2.0; origin[2] = 0.0; mitk::Vector3D normal; normal[0] = 0.0; normal[1] = 1.0; normal[2] = 0.0; myPlaneGeometry->InitializePlane(origin,normal); //generate points and line for intersection testing //point distance of given line > 1 mitk::Point3D pointP1; pointP1[0] = 2.0; pointP1[1] = 1.0; pointP1[2] = 0.0; mitk::Point3D pointP2; pointP2[0] = 2.0; pointP2[1] = 4.0; pointP2[2] = 0.0; mitk::Vector3D lineDirection; lineDirection[0] = pointP2[0] - pointP1[0]; lineDirection[1] = pointP2[1] - pointP1[1]; lineDirection[2] = pointP2[2] - pointP1[2]; mitk::Line3D xingline( pointP1, lineDirection ); mitk::Point3D calcXingPoint; myPlaneGeometry->IntersectionPoint(xingline, calcXingPoint); //point distance of given line < 1 mitk::Point3D pointP3; pointP3[0] = 2.0; pointP3[1] = 2.2; pointP3[2] = 0.0; mitk::Point3D pointP4; pointP4[0] = 2.0; pointP4[1] = 1.7; pointP4[2] = 0.0; mitk::Vector3D lineDirection2; lineDirection2[0] = pointP4[0] - pointP3[0]; lineDirection2[1] = pointP4[1] - pointP3[1]; lineDirection2[2] = pointP4[2] - pointP3[2]; mitk::Line3D xingline2( pointP3, lineDirection2 ); mitk::Point3D calcXingPoint2; myPlaneGeometry->IntersectionPoint( xingline2, calcXingPoint2 ); //intersection points must be the same CPPUNIT_ASSERT_MESSAGE("Failed to calculate Intersection Point", calcXingPoint == calcXingPoint2); } /** - * @brief This method tests method ProjectPointOntoPlane. - * - * See also bug #3409. - */ + * @brief This method tests method ProjectPointOntoPlane. + * + * See also bug #3409. + */ // Test does not use standard Parameters void TestProjectPointOntoPlane() { mitk::PlaneGeometry::Pointer myPlaneGeometry = mitk::PlaneGeometry::New(); //create normal mitk::Vector3D normal; normal[0] = 0.0; normal[1] = 0.0; normal[2] = 1.0; //create origin mitk::Point3D origin; origin[0] = -27.582859; origin[1] = 50; origin[2] = 200.27742; //initialize plane geometry myPlaneGeometry->InitializePlane(origin,normal); //output to descripe the test std::cout << "Testing PlaneGeometry according to bug #3409" << std::endl; std::cout << "Our normal is: " << normal << std::endl; std::cout << "So ALL projected points should have exactly the same z-value!" << std::endl; //create a number of points mitk::Point3D myPoints[5]; myPoints[0][0] = -27.582859; myPoints[0][1] = 50.00; myPoints[0][2] = 200.27742; myPoints[1][0] = -26.58662; myPoints[1][1] = 50.00; myPoints[1][2] = 200.19026; myPoints[2][0] = -26.58662; myPoints[2][1] = 50.00; myPoints[2][2] = 200.33124; myPoints[3][0] = 104.58662; myPoints[3][1] = 452.12313; myPoints[3][2] = 866.41236; myPoints[4][0] = -207.58662; myPoints[4][1] = 312.00; myPoints[4][2] = -300.12346; //project points onto plane mitk::Point3D myProjectedPoints[5]; for ( unsigned int i = 0; i < 5; ++i ) { myProjectedPoints[i] = myPlaneGeometry->ProjectPointOntoPlane( myPoints[i] ); } //compare z-values with z-value of plane (should be equal) bool allPointsOnPlane = true; for (auto & myProjectedPoint : myProjectedPoints) { if ( fabs(myProjectedPoint[2] - origin[2]) > mitk::sqrteps ) { allPointsOnPlane = false; } } CPPUNIT_ASSERT_MESSAGE("All points lie not on the same plane", allPointsOnPlane); } void TestPlaneGeometryCloning() { mitk::PlaneGeometry::Pointer geometry2D = createPlaneGeometry(); try { mitk::PlaneGeometry::Pointer clone = geometry2D->Clone(); itk::Matrix matrix = clone->GetIndexToWorldTransform()->GetMatrix(); CPPUNIT_ASSERT_MESSAGE("Test if matrix element exists...", matrix[0][0] == 31); double origin = geometry2D->GetOrigin()[0]; CPPUNIT_ASSERT_MESSAGE("First Point of origin as expected...", mitk::Equal(origin, 8)); double spacing = geometry2D->GetSpacing()[0]; CPPUNIT_ASSERT_MESSAGE("First Point of spacing as expected...", mitk::Equal(spacing, 31)); } catch (...) { CPPUNIT_FAIL("Error during access on a member of cloned geometry"); } // direction [row] [coloum] MITK_TEST_OUTPUT( << "Casting a rotated 2D ITK Image to a MITK Image and check if Geometry is still same" ); } void TestPlaneGeometryInitializeOrder() { mitk::Vector3D mySpacing; mySpacing[0] = 31; mySpacing[1] = 0.1; mySpacing[2] = 5.4; mitk::Point3D myOrigin; myOrigin[0] = 8; myOrigin[1] = 9; myOrigin[2] = 10; mitk::AffineTransform3D::Pointer myTransform = mitk::AffineTransform3D::New(); itk::Matrix transMatrix; transMatrix.Fill(0); transMatrix[0][0] = 1; transMatrix[1][1] = 2; transMatrix[2][2] = 4; myTransform->SetMatrix(transMatrix); mitk::PlaneGeometry::Pointer geometry2D1 = mitk::PlaneGeometry::New(); geometry2D1->SetIndexToWorldTransform(myTransform); geometry2D1->SetSpacing(mySpacing); geometry2D1->SetOrigin(myOrigin); mitk::PlaneGeometry::Pointer geometry2D2 = mitk::PlaneGeometry::New(); geometry2D2->SetSpacing(mySpacing); geometry2D2->SetOrigin(myOrigin); geometry2D2->SetIndexToWorldTransform(myTransform); mitk::PlaneGeometry::Pointer geometry2D3 = mitk::PlaneGeometry::New(); geometry2D3->SetIndexToWorldTransform(myTransform); geometry2D3->SetSpacing(mySpacing); geometry2D3->SetOrigin(myOrigin); geometry2D3->SetIndexToWorldTransform(myTransform); - CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 1 matches that of Geometry 2.", mitk::Equal(geometry2D1->GetOrigin(), geometry2D2->GetOrigin())); - CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 1 match those of Geometry 3.", mitk::Equal(geometry2D1->GetOrigin(), geometry2D3->GetOrigin())); - CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 2 match those of Geometry 3.", mitk::Equal(geometry2D2->GetOrigin(), geometry2D3->GetOrigin())); - - CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 1 match those of Geometry 2.", mitk::Equal(geometry2D1->GetSpacing(), geometry2D2->GetSpacing())); - CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 1 match those of Geometry 3.", mitk::Equal(geometry2D1->GetSpacing(), geometry2D3->GetSpacing())); - CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 2 match those of Geometry 3.", mitk::Equal(geometry2D2->GetSpacing(), geometry2D3->GetSpacing())); - - CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 1 match those of Geometry 2.", compareMatrix(geometry2D1->GetIndexToWorldTransform()->GetMatrix(), geometry2D2->GetIndexToWorldTransform()->GetMatrix())); - CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 1 match those of Geometry 3.", compareMatrix(geometry2D1->GetIndexToWorldTransform()->GetMatrix(), geometry2D3->GetIndexToWorldTransform()->GetMatrix())); - CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 2 match those of Geometry 3.", compareMatrix(geometry2D2->GetIndexToWorldTransform()->GetMatrix(), geometry2D3->GetIndexToWorldTransform()->GetMatrix())); + CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 1 matches that of Geometry 2.", + mitk::Equal(geometry2D1->GetOrigin(), geometry2D2->GetOrigin())); + CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 1 match those of Geometry 3.", + mitk::Equal(geometry2D1->GetOrigin(), geometry2D3->GetOrigin())); + CPPUNIT_ASSERT_MESSAGE("Origin of Geometry 2 match those of Geometry 3.", + mitk::Equal(geometry2D2->GetOrigin(), geometry2D3->GetOrigin())); + + CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 1 match those of Geometry 2.", + mitk::Equal(geometry2D1->GetSpacing(), geometry2D2->GetSpacing())); + CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 1 match those of Geometry 3.", + mitk::Equal(geometry2D1->GetSpacing(), geometry2D3->GetSpacing())); + CPPUNIT_ASSERT_MESSAGE("Spacing of Geometry 2 match those of Geometry 3.", + mitk::Equal(geometry2D2->GetSpacing(), geometry2D3->GetSpacing())); + + CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 1 match those of Geometry 2.", + compareMatrix(geometry2D1->GetIndexToWorldTransform()->GetMatrix(), geometry2D2->GetIndexToWorldTransform()->GetMatrix())); + CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 1 match those of Geometry 3.", + compareMatrix(geometry2D1->GetIndexToWorldTransform()->GetMatrix(), geometry2D3->GetIndexToWorldTransform()->GetMatrix())); + CPPUNIT_ASSERT_MESSAGE("Transformation of Geometry 2 match those of Geometry 3.", + compareMatrix(geometry2D2->GetIndexToWorldTransform()->GetMatrix(), geometry2D3->GetIndexToWorldTransform()->GetMatrix())); } void TestInitializeStandardPlane() { - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: width", mitk::Equal(planegeometry->GetExtent(0),width, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: height", mitk::Equal(planegeometry->GetExtent(1),height, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: depth", mitk::Equal(planegeometry->GetExtent(2),1, testEps)); - - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: width in mm", mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: heght in mm", mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: depth in mm", mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps) ); - - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorRight", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorBottom", mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorNormal", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: width", + mitk::Equal(planegeometry->GetExtent(0),width, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: height", + mitk::Equal(planegeometry->GetExtent(1),height, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: depth", + mitk::Equal(planegeometry->GetExtent(2),1, testEps)); + + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: width in mm", + mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: heght in mm", + mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: depth in mm", + mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps) ); + + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorRight", + mitk::Equal(planegeometry->GetAxisVector(0), right, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorBottom", + mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with default Spacing: AxisVectorNormal", + mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps) ); mitk::Vector3D spacing; thicknessInMM = 1.5; normal.Normalize(); normal *= thicknessInMM; mitk::FillVector3D(spacing, 1.0, 1.0, thicknessInMM); planegeometry->InitializeStandardPlane(right.GetVnlVector(), bottom.GetVnlVector(), &spacing); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: width", mitk::Equal(planegeometry->GetExtent(0),width, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: height", mitk::Equal(planegeometry->GetExtent(1),height, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: depth", mitk::Equal(planegeometry->GetExtent(2),1, testEps)); - - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: width in mm", mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: height in mm", mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: depth in mm", mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps) ); - - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorRight", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorBottom", mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps) ); - CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorNormal", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: width", + mitk::Equal(planegeometry->GetExtent(0),width, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: height", + mitk::Equal(planegeometry->GetExtent(1),height, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: depth", + mitk::Equal(planegeometry->GetExtent(2),1, testEps)); + + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: width in mm", + mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: height in mm", + mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: depth in mm", + mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps) ); + + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorRight", + mitk::Equal(planegeometry->GetAxisVector(0), right, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorBottom", + mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps) ); + CPPUNIT_ASSERT_MESSAGE("Testing correct Standard Plane initialization with custom Spacing: AxisVectorNormal", + mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps) ); ; } void TestSetExtendInMM() { normal.Normalize(); normal *= thicknessInMM; planegeometry->SetExtentInMM(2, thicknessInMM); - CPPUNIT_ASSERT_MESSAGE("Testing SetExtentInMM(2, ...), querying by GetExtentInMM(2): ", mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing SetExtentInMM(2, ...), querying by GetAxisVector(2) and comparing to normal: ", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing SetExtentInMM(2, ...), querying by GetExtentInMM(2): ", + mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing SetExtentInMM(2, ...), querying by GetAxisVector(2) and comparing to normal: ", + mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); planegeometry->SetOrigin(origin); CPPUNIT_ASSERT_MESSAGE("Testing SetOrigin", mitk::Equal(planegeometry->GetOrigin(), origin, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() after SetOrigin: Right", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps) ); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() after SetOrigin: Bottom", mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps) ); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() after SetOrigin: Normal", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps) ); mappingTests2D(planegeometry, width, height, widthInMM, heightInMM, origin, right, bottom); } void TestRotate() { // Changing the IndexToWorldTransform to a rotated version by SetIndexToWorldTransform() (keep origin): mitk::AffineTransform3D::Pointer transform = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = planegeometry->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix(); mitk::VnlVector axis(3); mitk::FillVector3D(axis, 1.0, 1.0, 1.0); axis.normalize(); vnl_quaternion rotation(axis, 0.223); vnlmatrix = rotation.rotation_matrix_transpose()*vnlmatrix; mitk::Matrix3D matrix; matrix = vnlmatrix; transform->SetMatrix(matrix); transform->SetOffset(planegeometry->GetIndexToWorldTransform()->GetOffset()); right.SetVnlVector( rotation.rotation_matrix_transpose()*right.GetVnlVector() ); bottom.SetVnlVector(rotation.rotation_matrix_transpose()*bottom.GetVnlVector()); normal.SetVnlVector(rotation.rotation_matrix_transpose()*normal.GetVnlVector()); planegeometry->SetIndexToWorldTransform(transform); //The origin changed,because m_Origin=m_IndexToWorldTransform->GetOffset()+GetAxisVector(2)*0.5 //and the AxisVector changes due to the rotation. In other words: the rotation was done around //the corner of the box, not around the planes origin. Now change it to a rotation around //the origin, simply by re-setting the origin to the original one: planegeometry->SetOrigin(origin); CPPUNIT_ASSERT_MESSAGE("Testing whether SetIndexToWorldTransform kept origin: ", mitk::Equal(planegeometry->GetOrigin(), origin, testEps) ); mitk::Point2D point; point[0] = 4; point[1] = 3; mitk::Point2D dummy; planegeometry->WorldToIndex(point, dummy); planegeometry->IndexToWorld(dummy, dummy); - CPPUNIT_ASSERT_MESSAGE("Testing consistancy of index and world coordinates.", dummy == point); + CPPUNIT_ASSERT_MESSAGE("Testing consistency of index and world coordinates.", dummy == point); CPPUNIT_ASSERT_MESSAGE("Testing width of rotated version: ", mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height of rotated version: ", mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness of rotated version: ", mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of rotated version: right ", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of rotated version: bottom", mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of rotated version: normal", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(0).GetNorm(),planegeometry->GetExtentInMM(0), testEps)); + mitk::Equal(planegeometry->GetAxisVector(0).GetNorm(),planegeometry->GetExtentInMM(0), testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(1).GetNorm(),planegeometry->GetExtentInMM(1), testEps)); + mitk::Equal(planegeometry->GetAxisVector(1).GetNorm(),planegeometry->GetExtentInMM(1), testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(2).GetNorm(),planegeometry->GetExtentInMM(2), testEps)); + mitk::Equal(planegeometry->GetAxisVector(2).GetNorm(),planegeometry->GetExtentInMM(2), testEps)); mappingTests2D(planegeometry, width, height, widthInMM, heightInMM, origin, right, bottom); width *= 2; height *= 3; planegeometry->SetSizeInUnits(width, height); - CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: ", mitk::Equal(planegeometry->GetExtent(0),width, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: ", mitk::Equal(planegeometry->GetExtent(1),height, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: ", mitk::Equal(planegeometry->GetExtent(2),1, testEps)); - - CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of version with changed size in units: ", mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of version with changed size in units: ", mitk::Equal(planegeometry->GetExtentInMM(1), heightInMM, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of version with changed size in units: ", mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps)); - - CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: right ", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: bottom", mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps)); - CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: normal", mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: ", + mitk::Equal(planegeometry->GetExtent(0),width, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: ", + mitk::Equal(planegeometry->GetExtent(1),height, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing SetSizeInUnits() of rotated version: " +, mitk::Equal(planegeometry->GetExtent(2),1, testEps)); + + CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of version with changed size in units: ", + mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of version with changed size in units: ", + mitk::Equal(planegeometry->GetExtentInMM(1), heightInMM, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of version with changed size in units: ", + mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps)); + + CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: right ", + mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: bottom", + mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps)); + CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of version with changed size in units: normal", + mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(0).GetNorm(),planegeometry->GetExtentInMM(0), testEps)); + mitk::Equal(planegeometry->GetAxisVector(0).GetNorm(),planegeometry->GetExtentInMM(0), testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(1).GetNorm(),planegeometry->GetExtentInMM(1), testEps)); + mitk::Equal(planegeometry->GetAxisVector(1).GetNorm(),planegeometry->GetExtentInMM(1), testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ", - mitk::Equal(planegeometry->GetAxisVector(2).GetNorm(),planegeometry->GetExtentInMM(2), testEps)); + mitk::Equal(planegeometry->GetAxisVector(2).GetNorm(),planegeometry->GetExtentInMM(2), testEps)); mappingTests2D(planegeometry, width, height, widthInMM, heightInMM, origin, right, bottom); } void TestClone() { mitk::PlaneGeometry::Pointer clonedplanegeometry = dynamic_cast(planegeometry->Clone().GetPointer()); // Cave: Statement below is negated! CPPUNIT_ASSERT_MESSAGE("Testing Clone(): ", ! ((clonedplanegeometry.IsNull()) || (clonedplanegeometry->GetReferenceCount()!=1))); CPPUNIT_ASSERT_MESSAGE("Testing origin of cloned version: ", mitk::Equal(clonedplanegeometry->GetOrigin(), origin, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width (in units) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtent(0),width, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in units) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtent(1),height, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing extent (in units) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtent(2),1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtentInMM(1), heightInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of cloned version: ", mitk::Equal(clonedplanegeometry->GetExtentInMM(2), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of cloned version: right", mitk::Equal(clonedplanegeometry->GetAxisVector(0), right, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of cloned version: bottom", mitk::Equal(clonedplanegeometry->GetAxisVector(1), bottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of cloned version: normal", mitk::Equal(clonedplanegeometry->GetAxisVector(2), normal, testEps)); mappingTests2D(clonedplanegeometry, width, height, widthInMM, heightInMM, origin, right, bottom); } void TestSaggitalInitialization() { mitk::Point3D cornerpoint0 = planegeometry->GetCornerPoint(0); mitk::PlaneGeometry::Pointer clonedplanegeometry = dynamic_cast(planegeometry->Clone().GetPointer()); // Testing InitializeStandardPlane(clonedplanegeometry, planeorientation = Sagittal, zPosition = 0, frontside=true): planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::PlaneGeometry::Sagittal); mitk::Vector3D newright, newbottom, newnormal; mitk::ScalarType newthicknessInMM; newright = bottom; newthicknessInMM = widthInMM/width*1.0; // extent in normal direction is 1; newnormal = right; newnormal.Normalize(); newnormal *= newthicknessInMM; newbottom = normal; newbottom.Normalize(); newbottom *= thicknessInMM; CPPUNIT_ASSERT_MESSAGE("Testing GetCornerPoint(0) of sagitally initialized version:", mitk::Equal(planegeometry->GetCornerPoint(0), cornerpoint0, testEps)); //ok, corner was fine, so we can dare to believe the origin is ok. origin = planegeometry->GetOrigin(); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtent(0), height, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtent(1), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(0), heightInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagitally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagitally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagitally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagitally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(2), newnormal, testEps)); mappingTests2D(planegeometry, height, 1, heightInMM, thicknessInMM, origin, newright, newbottom); // set origin back to the one of the axial slice: origin = clonedplanegeometry->GetOrigin(); // Testing backside initialization: InitializeStandardPlane(clonedplanegeometry, planeorientation = Axial, zPosition = 0, frontside=false, rotated=true): planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::PlaneGeometry::Axial, 0, false, true); mitk::Point3D backsideorigin; backsideorigin=origin+clonedplanegeometry->GetAxisVector(1);//+clonedplanegeometry->GetAxisVector(2); CPPUNIT_ASSERT_MESSAGE("Testing origin of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetOrigin(), backsideorigin, testEps)); mitk::Point3D backsidecornerpoint0; backsidecornerpoint0 = cornerpoint0+clonedplanegeometry->GetAxisVector(1);//+clonedplanegeometry->GetAxisVector(2); CPPUNIT_ASSERT_MESSAGE("Testing GetCornerPoint(0) of sagitally initialized version: ", mitk::Equal(planegeometry->GetCornerPoint(0), backsidecornerpoint0, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of backsidedly, axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(0), width, testEps)); + mitk::Equal(planegeometry->GetExtent(0), width, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of backsidedly, axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(1), height, testEps)); + mitk::Equal(planegeometry->GetExtent(1), height, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of backsidedly, axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); + mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(1), heightInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetAxisVector(1), -bottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of backsidedly, axially initialized version: ", mitk::Equal(planegeometry->GetAxisVector(2), -normal, testEps)); mappingTests2D(planegeometry, width, height, widthInMM, heightInMM, backsideorigin, right, -bottom); } void TestFrontalInitialization() { mitk::Point3D cornerpoint0 = planegeometry->GetCornerPoint(0); mitk::PlaneGeometry::Pointer clonedplanegeometry = dynamic_cast(planegeometry->Clone().GetPointer()); //-------- mitk::Vector3D newright, newbottom, newnormal; mitk::ScalarType newthicknessInMM; // Testing InitializeStandardPlane(clonedplanegeometry, planeorientation = Frontal, zPosition = 0, frontside=true) planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::PlaneGeometry::Frontal); newright = right; newbottom = normal; newbottom.Normalize(); newbottom *= thicknessInMM; newthicknessInMM = heightInMM/height*1.0/*extent in normal direction is 1*/; newnormal = -bottom; newnormal.Normalize(); newnormal *= newthicknessInMM; CPPUNIT_ASSERT_MESSAGE("Testing GetCornerPoint(0) of frontally initialized version: ", mitk::Equal(planegeometry->GetCornerPoint(0), cornerpoint0, testEps)); //ok, corner was fine, so we can dare to believe the origin is ok. origin = planegeometry->GetOrigin(); CPPUNIT_ASSERT_MESSAGE("Testing width (in units) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(0), width, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in units) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(1), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness (in units) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(2), newnormal, testEps)); mappingTests2D(planegeometry, width, 1, widthInMM, thicknessInMM, origin, newright, newbottom); //Changing plane to in-plane unit spacing using SetSizeInUnits: planegeometry->SetSizeInUnits(planegeometry->GetExtentInMM(0), planegeometry->GetExtentInMM(1)); CPPUNIT_ASSERT_MESSAGE( "Testing origin of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetOrigin(), origin, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(2), newnormal, testEps)); mappingTests2D(planegeometry, widthInMM, thicknessInMM, widthInMM, thicknessInMM, origin, newright, newbottom); // Changing plane to unit spacing also in normal direction using SetExtentInMM(2, 1.0): planegeometry->SetExtentInMM(2, 1.0); newnormal.Normalize(); CPPUNIT_ASSERT_MESSAGE("Testing origin of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetOrigin(), origin, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetExtentInMM(2), 1.0, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, frontally initialized version: ", mitk::Equal(planegeometry->GetAxisVector(2), newnormal, testEps)); mappingTests2D(planegeometry, widthInMM, thicknessInMM, widthInMM, thicknessInMM, origin, newright, newbottom); } void TestAxialInitialization() { mitk::Point3D cornerpoint0 = planegeometry->GetCornerPoint(0); // Clone, move, rotate and test for 'IsParallel' and 'IsOnPlane' mitk::PlaneGeometry::Pointer clonedplanegeometry = dynamic_cast(planegeometry->Clone().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Testing Clone(): ", ! ((clonedplanegeometry.IsNull()) || (clonedplanegeometry->GetReferenceCount()!=1)) ); std::cout << "Testing InitializeStandardPlane(clonedplanegeometry, planeorientation = Axial, zPosition = 0, frontside=true): " <InitializeStandardPlane(clonedplanegeometry); CPPUNIT_ASSERT_MESSAGE("Testing origin of axially initialized version: ", mitk::Equal(planegeometry->GetOrigin(), origin)); - CPPUNIT_ASSERT_MESSAGE("Testing GetCornerPoint(0) of axially initialized version: ", mitk::Equal(planegeometry->GetCornerPoint(0), cornerpoint0)); + CPPUNIT_ASSERT_MESSAGE("Testing GetCornerPoint(0) of axially initialized version: ", + mitk::Equal(planegeometry->GetCornerPoint(0), cornerpoint0)); CPPUNIT_ASSERT_MESSAGE("Testing width (in units) of axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(0), width, testEps)); + mitk::Equal(planegeometry->GetExtent(0), width, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in units) of axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(1), height, testEps)); + mitk::Equal(planegeometry->GetExtent(1), height, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness (in units) of axially initialized version (should be same as in mm due to unit spacing, except for thickness, which is always 1): ", - mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); + mitk::Equal(planegeometry->GetExtent(2), 1, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of axially initialized version: ", - mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); + mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of axially initialized version: ", - mitk::Equal(planegeometry->GetExtentInMM(1), heightInMM, testEps)); + mitk::Equal(planegeometry->GetExtentInMM(1), heightInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of axially initialized version: ", - mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps)); + mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ", - mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); + mitk::Equal(planegeometry->GetAxisVector(0), right, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ", - mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps)); + mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps)); CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ", - mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); + mitk::Equal(planegeometry->GetAxisVector(2), normal, testEps)); mappingTests2D(planegeometry, width, height, widthInMM, heightInMM, origin, right, bottom); } void TestPlaneComparison() { // Clone, move, rotate and test for 'IsParallel' and 'IsOnPlane' mitk::PlaneGeometry::Pointer clonedplanegeometry2 = dynamic_cast(planegeometry->Clone().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Testing Clone(): ", ! ((clonedplanegeometry2.IsNull()) || (clonedplanegeometry2->GetReferenceCount()!=1)) ); CPPUNIT_ASSERT_MESSAGE("Testing wheter original and clone are at the same position", clonedplanegeometry2->IsOnPlane(planegeometry.GetPointer())); CPPUNIT_ASSERT_MESSAGE(" Asserting that origin is on the plane cloned plane:", clonedplanegeometry2->IsOnPlane(origin)); mitk::VnlVector newaxis(3); mitk::FillVector3D(newaxis, 1.0, 1.0, 1.0); newaxis.normalize(); vnl_quaternion rotation2(newaxis, 0.0); mitk::Vector3D clonednormal = clonedplanegeometry2->GetNormal(); mitk::Point3D clonedorigin = clonedplanegeometry2->GetOrigin(); auto planerot = new mitk::RotationOperation( mitk::OpROTATE, origin, clonedplanegeometry2->GetAxisVector( 0 ), 180.0 ); clonedplanegeometry2->ExecuteOperation( planerot ); CPPUNIT_ASSERT_MESSAGE(" Asserting that a flipped plane is still on the original plane: ", clonedplanegeometry2->IsOnPlane(planegeometry.GetPointer())); clonedorigin += clonednormal; clonedplanegeometry2->SetOrigin( clonedorigin ); CPPUNIT_ASSERT_MESSAGE("Testing if the translated (cloned, flipped) plane is parallel to its origin plane: ", clonedplanegeometry2->IsParallel(planegeometry)); delete planerot; planerot = new mitk::RotationOperation( mitk::OpROTATE, origin, clonedplanegeometry2->GetAxisVector( 0 ), 0.5 ); clonedplanegeometry2->ExecuteOperation( planerot ); CPPUNIT_ASSERT_MESSAGE("Testing if a non-paralell plane gets recognized as not paralell [rotation +0.5 degree] : ", ! clonedplanegeometry2->IsParallel(planegeometry)); delete planerot; planerot = new mitk::RotationOperation( mitk::OpROTATE, origin, clonedplanegeometry2->GetAxisVector( 0 ), -1.0 ); clonedplanegeometry2->ExecuteOperation( planerot ); CPPUNIT_ASSERT_MESSAGE("Testing if a non-paralell plane gets recognized as not paralell [rotation -0.5 degree] : ", ! clonedplanegeometry2->IsParallel(planegeometry)); delete planerot; planerot = new mitk::RotationOperation( mitk::OpROTATE, origin, clonedplanegeometry2->GetAxisVector( 0 ), 360.5 ); clonedplanegeometry2->ExecuteOperation( planerot ); CPPUNIT_ASSERT_MESSAGE("Testing if a non-paralell plane gets recognized as paralell [rotation 360 degree] : ", clonedplanegeometry2->IsParallel(planegeometry)); } private: // helper Methods for the Tests mitk::PlaneGeometry::Pointer createPlaneGeometry() { mitk::Vector3D mySpacing; mySpacing[0] = 31; mySpacing[1] = 0.1; mySpacing[2] = 5.4; mitk::Point3D myOrigin; myOrigin[0] = 8; myOrigin[1] = 9; myOrigin[2] = 10; mitk::AffineTransform3D::Pointer myTransform = mitk::AffineTransform3D::New(); itk::Matrix transMatrix; transMatrix.Fill(0); transMatrix[0][0] = 1; transMatrix[1][1] = 2; transMatrix[2][2] = 4; myTransform->SetMatrix(transMatrix); mitk::PlaneGeometry::Pointer geometry2D = mitk::PlaneGeometry::New(); geometry2D->SetIndexToWorldTransform(myTransform); geometry2D->SetSpacing(mySpacing); geometry2D->SetOrigin(myOrigin); return geometry2D; } bool compareMatrix(itk::Matrix left, itk::Matrix right) { bool equal = true; for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) equal &= mitk::Equal(left[i][j], right[i][j]); return equal; } /** - * This function tests for correct mapping and is called several times from other tests - **/ + * This function tests for correct mapping and is called several times from other tests + **/ void mappingTests2D(const mitk::PlaneGeometry* planegeometry, const mitk::ScalarType& width, const mitk::ScalarType& height, const mitk::ScalarType& widthInMM, const mitk::ScalarType& heightInMM, const mitk::Point3D& origin, const mitk::Vector3D& right, const mitk::Vector3D& bottom) { std::cout << "Testing mapping Map(pt2d_mm(x=widthInMM/2.3,y=heightInMM/2.5), pt3d_mm) and compare with expected: "; mitk::Point2D pt2d_mm; mitk::Point3D pt3d_mm, expected_pt3d_mm; pt2d_mm[0] = widthInMM/2.3; pt2d_mm[1] = heightInMM/2.5; expected_pt3d_mm = origin+right*(pt2d_mm[0]/right.GetNorm())+bottom*(pt2d_mm[1]/bottom.GetNorm()); planegeometry->Map(pt2d_mm, pt3d_mm); CPPUNIT_ASSERT_MESSAGE("Testing mapping Map(pt2d_mm(x=widthInMM/2.3,y=heightInMM/2.5), pt3d_mm) and compare with expected", mitk::Equal(pt3d_mm, expected_pt3d_mm, testEps)); std::cout << "Testing mapping Map(pt3d_mm, pt2d_mm) and compare with expected: "; mitk::Point2D testpt2d_mm; planegeometry->Map(pt3d_mm, testpt2d_mm); std::cout << std::setprecision(12) << "Expected pt2d_mm " << pt2d_mm << std::endl; std::cout << std::setprecision(12) << "Result testpt2d_mm " << testpt2d_mm << std::endl; std::cout << std::setprecision(12) << "10*mitk::eps " << 10*mitk::eps << std::endl; //This eps is temporarily set to 10*mitk::eps. See bug #15037 for details. CPPUNIT_ASSERT_MESSAGE("Testing mapping Map(pt3d_mm, pt2d_mm) and compare with expected", mitk::Equal(pt2d_mm, testpt2d_mm, 10*mitk::eps)); std::cout << "Testing IndexToWorld(pt2d_units, pt2d_mm) and compare with expected: "; mitk::Point2D pt2d_units; pt2d_units[0] = width/2.0; pt2d_units[1] = height/2.0; pt2d_mm[0] = widthInMM/2.0; pt2d_mm[1] = heightInMM/2.0; planegeometry->IndexToWorld(pt2d_units, testpt2d_mm); std::cout << std::setprecision(12) << "Expected pt2d_mm " << pt2d_mm << std::endl; std::cout << std::setprecision(12) << "Result testpt2d_mm " << testpt2d_mm << std::endl; std::cout << std::setprecision(12) << "10*mitk::eps " << 10*mitk::eps << std::endl; //This eps is temporarily set to 10*mitk::eps. See bug #15037 for details. CPPUNIT_ASSERT_MESSAGE("Testing IndexToWorld(pt2d_units, pt2d_mm) and compare with expected: ", mitk::Equal(pt2d_mm, testpt2d_mm, 10*mitk::eps)); std::cout << "Testing WorldToIndex(pt2d_mm, pt2d_units) and compare with expected: "; mitk::Point2D testpt2d_units; planegeometry->WorldToIndex(pt2d_mm, testpt2d_units); std::cout << std::setprecision(12) << "Expected pt2d_units " << pt2d_units << std::endl; std::cout << std::setprecision(12) << "Result testpt2d_units " << testpt2d_units << std::endl; std::cout << std::setprecision(12) << "10*mitk::eps " << 10*mitk::eps << std::endl; //This eps is temporarily set to 10*mitk::eps. See bug #15037 for details. CPPUNIT_ASSERT_MESSAGE("Testing WorldToIndex(pt2d_mm, pt2d_units) and compare with expected:", mitk::Equal(pt2d_units, testpt2d_units, 10*mitk::eps)); } }; -MITK_TEST_SUITE_REGISTRATION(mitkPlaneGeometry) \ No newline at end of file +MITK_TEST_SUITE_REGISTRATION(mitkPlaneGeometry)