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 <MitkCoreExports.h>
#include "mitkBaseGeometry.h"
#include "mitkRestorePlanePositionOperation.h"
#include <vnl/vnl_cross.h>
namespace mitk {
template < class TCoordRep, unsigned int NPointDimension > class Line;
typedef Line<ScalarType, 3> 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 <vtkTransform.h>
#include <vnl/vnl_cross.h>
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<double, 3> 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<float>::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<double, 3> 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<BoundingBox*>(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<BoundingBox*>(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<BoundingBox*>(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 <mitkTestingMacros.h>
#include <mitkTestFixture.h>
#include <vnl/vnl_quaternion.h>
#include <vnl/vnl_quaternion.txx>
#include <fstream>
#include <iomanip>
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): "<<std::endl;
CPPUNIT_ASSERT_NO_THROW(planegeometry->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]"<<std::endl;
return EXIT_FAILURE;
}
std::cout<<"[PASSED]"<<std::endl;
std::cout << "Testing width, height and thickness (in mm): ";
if((mitk::Equal(planegeometry->GetExtentInMM(0),widthInMM)==false) ||
(mitk::Equal(planegeometry->GetExtentInMM(1),heightInMM)==false) ||
(mitk::Equal(planegeometry->GetExtentInMM(2),thicknessInMM)==false)
)
{
std::cout<<"[FAILED]"<<std::endl;
return EXIT_FAILURE;
}
- */
+ */
}
/**
- * @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.
- *
- */
+ * @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<mitk::ScalarType,3,3> 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<mitk::ScalarType, 3,3> 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<mitk::ScalarType> 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<mitk::PlaneGeometry*>(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<mitk::PlaneGeometry*>(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<mitk::PlaneGeometry*>(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<mitk::PlaneGeometry*>(planegeometry->Clone().GetPointer());
CPPUNIT_ASSERT_MESSAGE("Testing Clone(): ", ! ((clonedplanegeometry.IsNull()) || (clonedplanegeometry->GetReferenceCount()!=1)) );
std::cout << "Testing InitializeStandardPlane(clonedplanegeometry, planeorientation = Axial, zPosition = 0, frontside=true): " <<std::endl;
planegeometry->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<mitk::PlaneGeometry*>(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<mitk::ScalarType> 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<mitk::ScalarType, 3,3> 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<mitk::ScalarType, 3,3> left, itk::Matrix<mitk::ScalarType, 3,3> 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)