diff --git a/Modules/Core/include/mitkPlaneGeometry.h b/Modules/Core/include/mitkPlaneGeometry.h
index 4e296c981b..c287d67541 100644
--- a/Modules/Core/include/mitkPlaneGeometry.h
+++ b/Modules/Core/include/mitkPlaneGeometry.h
@@ -1,608 +1,611 @@
 /*===================================================================
 
 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 "mitkBaseGeometry.h"
 #include "mitkRestorePlanePositionOperation.h"
 #include <MitkCoreExports.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,
                                          bool top = true);
 
     /**
     * \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,
                                          bool top = true);
 
     /**
     * \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,
                                          bool top = true);
 
     /**
     * \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
-    *
+    * \brief Check if matrix is a rotation matrix:
+    * - determinant is 1?
+    * - R*R^T is ID?
+    * Output warning otherwise.
     */
-    static void EnsurePerpendicularNormal(AffineTransform3D *transform);
+    static bool CheckRotationMatrix(AffineTransform3D *transform, double epsilon=mitk::eps);
 
     /**
     * \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;
 
     itk::LightObject::Pointer InternalClone() const override;
 
     /** Implements operation to re-orient the plane */
     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(const mitk::BaseGeometry *geometry);
 
     /**
     * \brief Get the geometrical frame of reference for this PlaneGeometry.
     */
     const BaseGeometry *GetReferenceGeometry() const;
     bool HasReferenceGeometry() const;
 
+    static std::vector< int > CalculateDominantAxes(mitk::AffineTransform3D::MatrixType::InternalMatrixType& rotation_matrix);
+
   protected:
     PlaneGeometry();
 
     PlaneGeometry(const PlaneGeometry &other);
 
     ~PlaneGeometry() override;
 
     void PrintSelf(std::ostream &os, itk::Indent indent) const override;
 
     const 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();};
     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!)
     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);.
     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 b0868e59b5..f1d1afbc86 100644
--- a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
+++ b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
@@ -1,984 +1,975 @@
 /*===================================================================
 
 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 "mitkInteractionConst.h"
 #include "mitkLine.h"
 #include "mitkPlaneOperation.h"
 
 #include <itkSpatialOrientationAdapter.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)
   {
   }
 
-  void PlaneGeometry::EnsurePerpendicularNormal(mitk::AffineTransform3D *transform)
+  bool PlaneGeometry::CheckRotationMatrix(mitk::AffineTransform3D *transform, double epsilon)
   {
-    /** \brief ensure column(2) of indexToWorldTransform-matrix to be perpendicular to plane, keep length and
-     * handedness.
-     */
+    bool rotation = true;
 
-    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.
+    auto matrix = transform->GetMatrix().GetVnlMatrix();
+    matrix.normalize_columns();
 
-    ScalarType len = transform->GetMatrix().GetVnlMatrix().get_column(2).two_norm();
-    if (len == 0)
+    auto det = vnl_determinant(matrix);
+    if (fabs(det-1.0) > epsilon)
     {
-      len = 1;
+      MITK_WARN << "Matrix is not a rotation matrix! Determinant≠1 (" << det << ")";
+      rotation = false;
     }
-    normal *= len;
 
-    // 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)
+    vnl_matrix_fixed<double, 3, 3> id; id.set_identity();
+    auto should_be_id = matrix*matrix.transpose();
+    should_be_id -= id;
+    auto max = should_be_id.absolute_value_max();
+    if (max > epsilon)
     {
-      MITK_DEBUG << "EnsurePerpendicularNormal(): Lefthanded geometry preserved, rh-normal: [ " << normal << " ],";
-      normal *= (-1.0);
-      MITK_DEBUG << "lh-normal: [ " << normal << " ], original vector zed is: [ " << zed << " ]";
+      MITK_WARN << "Matrix is not a rotation matrix! R*R^T≠ID. Max value: " << max << " (should be 0)";
+      rotation = false;
     }
 
-    // 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.
-    }
+    return rotation;
   }
 
   void PlaneGeometry::CheckIndexToWorldTransform(mitk::AffineTransform3D *transform)
   {
-    EnsurePerpendicularNormal(transform);
+    CheckRotationMatrix(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,
                                               bool top)
   {
     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, top);
   }
 
   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 */,
                                               bool top /* = true */)
   {
     Superclass::Initialize();
 
     /// construct standard view.
 
     // 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):
     // 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.
 
     // 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;
     /** 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 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---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.)
           {    // 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.)
           {
             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;
 
       // 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;
 
       default:
         itkExceptionMacro("unknown PlaneOrientation");
     }
 
     VnlVector normal(3);
     FillVector3D(normal, 0, 0, 0);
     normal[normalDirection] = top ? 1 : -1;
 
     if ( transform != nullptr )
     {
       origin = transform->TransformPoint( origin );
       rightDV = transform->TransformVector( rightDV );
       bottomDV = transform->TransformVector( bottomDV );
       normal = transform->TransformVector( normal );
     }
 
     ScalarType bounds[6] = {0, width, 0, height, 0, 1};
     this->SetBounds(bounds);
 
     AffineTransform3D::Pointer planeTransform = AffineTransform3D::New();
     Matrix3D matrix;
     matrix.GetVnlMatrix().set_column(0, rightDV);
     matrix.GetVnlMatrix().set_column(1, bottomDV);
     matrix.GetVnlMatrix().set_column(2, normal);
     planeTransform->SetMatrix(matrix);
     planeTransform->SetOffset(this->GetIndexToWorldTransform()->GetOffset());
     this->SetIndexToWorldTransform(planeTransform);
 
     this->SetOrigin(origin);
   }
 
+  std::vector< int > PlaneGeometry::CalculateDominantAxes(mitk::AffineTransform3D::MatrixType::InternalMatrixType& rotation_matrix)
+  {
+    std::vector< int > axes;
+
+    bool dominant_axis_error = false;
+    for (int i = 0; i < 3; ++i)
+    {
+      int dominantAxis = itk::Function::Max3(
+          rotation_matrix[0][i],
+          rotation_matrix[1][i],
+          rotation_matrix[2][i]
+      );
+
+      for (int j=0; j<i; ++j)
+        if (axes[j] == dominantAxis)
+        {
+          dominant_axis_error = true;
+          break;
+        }
+      if (dominant_axis_error)
+        break;
+
+      axes.push_back(dominantAxis);
+    }
+
+    if (dominant_axis_error)
+    {
+      MITK_DEBUG << "Error during dominant axis calculation. Using default.";
+      MITK_DEBUG << "This is either caused by an imperfect rotation matrix or if the rotation is axactly 45° around one or more axis.";
+      axes.clear();
+      for (int i = 0; i < 3; ++i)
+        axes.push_back(i);
+    }
+
+    return axes;
+  }
+
   void PlaneGeometry::InitializeStandardPlane(const BaseGeometry *geometry3D,
                                               PlaneOrientation planeorientation,
                                               ScalarType zPosition,
                                               bool frontside,
                                               bool rotated,
                                               bool top)
   {
     this->SetReferenceGeometry(geometry3D);
 
     ScalarType width, height;
 
     // Inspired by:
     // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
 
     mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
 
     matrix.GetVnlMatrix().normalize_columns();
-    mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
+    mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetTranspose();
 
     /// The index of the sagittal, coronal and axial axes in the reference geometry.
-    int axes[3];
+    auto axes = CalculateDominantAxes(inverseMatrix);
     /// The direction of the sagittal, coronal and axial axes in the reference geometry.
     /// +1 means that the direction is straight, i.e. greater index translates to greater
     /// world coordinate. -1 means that the direction is inverted.
     int directions[3];
     ScalarType extents[3];
     ScalarType spacings[3];
-    for (int i = 0; i < 3; ++i)
+    for (int i=0; i<3; ++i)
     {
-      int dominantAxis = itk::Function::Max3(
-          inverseMatrix[0][i],
-          inverseMatrix[1][i],
-          inverseMatrix[2][i]
-      );
-      axes[i] = dominantAxis;
+      int dominantAxis = axes.at(i);
       directions[i] = itk::Function::Sign(inverseMatrix[dominantAxis][i]);
       extents[i] = geometry3D->GetExtent(dominantAxis);
       spacings[i] = geometry3D->GetSpacing()[dominantAxis];
     }
 
     // matrix(column) = inverseTransformMatrix(row) * flippedAxes * spacing
     matrix[0][0] = inverseMatrix[axes[0]][0] * directions[0] * spacings[0];
     matrix[1][0] = inverseMatrix[axes[0]][1] * directions[0] * spacings[0];
     matrix[2][0] = inverseMatrix[axes[0]][2] * directions[0] * spacings[0];
     matrix[0][1] = inverseMatrix[axes[1]][0] * directions[1] * spacings[1];
     matrix[1][1] = inverseMatrix[axes[1]][1] * directions[1] * spacings[1];
     matrix[2][1] = inverseMatrix[axes[1]][2] * directions[1] * spacings[1];
     matrix[0][2] = inverseMatrix[axes[2]][0] * directions[2] * spacings[2];
     matrix[1][2] = inverseMatrix[axes[2]][1] * directions[2] * spacings[2];
     matrix[2][2] = inverseMatrix[axes[2]][2] * directions[2] * spacings[2];
 
     /// The "world origin" is the corner with the lowest physical coordinates.
     /// We use it as a reference point so that we get the correct anatomical
     /// orientations.
     Point3D worldOrigin = geometry3D->GetOrigin();
     for (int i = 0; i < 3; ++i)
     {
       /// The distance of the plane origin from the world origin in voxels.
       double offset = directions[i] > 0 ? 0.0 : extents[i];
 
       if (geometry3D->GetImageGeometry())
       {
         offset += directions[i] * 0.5;
       }
 
       for (int j = 0; j < 3; ++j)
       {
         worldOrigin[j] -= offset * matrix[j][i];
       }
     }
 
     switch(planeorientation)
     {
     case None:
     /** Orientation 'None' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case Axial:
       width  = extents[0];
       height = extents[1];
       break;
     case Frontal:
       width  = extents[0];
       height = extents[2];
       break;
     case Sagittal:
       width  = extents[1];
       height = extents[2];
       break;
     default:
       itkExceptionMacro("unknown PlaneOrientation");
     }
 
     ScalarType bounds[6]= { 0, width, 0, height, 0, 1 };
     this->SetBounds( bounds );
 
     AffineTransform3D::Pointer transform = AffineTransform3D::New();
     transform->SetMatrix(matrix);
     transform->SetOffset(worldOrigin.GetVectorFromOrigin());
 
     InitializeStandardPlane(
       width, height, transform, planeorientation, zPosition, frontside, rotated, top);
   }
 
   void PlaneGeometry::InitializeStandardPlane(
     const BaseGeometry *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated)
   {
     /// The index of the sagittal, coronal and axial axes in world coordinate system.
     int worldAxis;
     switch(planeorientation)
     {
     case None:
     /** Orientation 'None' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case Axial:
       worldAxis = 2;
       break;
     case Frontal:
       worldAxis = 1;
       break;
     case Sagittal:
       worldAxis = 0;
       break;
     default:
       itkExceptionMacro("unknown PlaneOrientation");
     }
 
     // Inspired by:
     // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
 
     mitk::AffineTransform3D::ConstPointer affineTransform = geometry3D->GetIndexToWorldTransform();
     mitk::AffineTransform3D::MatrixType matrix = affineTransform->GetMatrix();
     matrix.GetVnlMatrix().normalize_columns();
     mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
 
     /// The index of the sagittal, coronal and axial axes in the reference geometry.
-    int dominantAxis = itk::Function::Max3(
-        inverseMatrix[0][worldAxis],
-        inverseMatrix[1][worldAxis],
-        inverseMatrix[2][worldAxis]);
+    int dominantAxis = CalculateDominantAxes(inverseMatrix).at(worldAxis);
 
     ScalarType zPosition = top ? 0.5 : geometry3D->GetExtent(dominantAxis) - 0.5;
 
     InitializeStandardPlane(geometry3D, planeorientation, zPosition, frontside, rotated, top);
   }
 
   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.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:
       {
         auto *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:
       {
         auto *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 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 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 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(const mitk::BaseGeometry *geometry) { m_ReferenceGeometry = geometry; }
   const mitk::BaseGeometry *PlaneGeometry::GetReferenceGeometry() const { return m_ReferenceGeometry; }
   bool PlaneGeometry::HasReferenceGeometry() const { return (m_ReferenceGeometry != nullptr); }
 } // namespace
diff --git a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
index 9a3a8e1e4e..4e4cd940a1 100644
--- a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
+++ b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
@@ -1,970 +1,966 @@
 /*===================================================================
 
 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 <itkSpatialOrientationAdapter.h>
 
 #include "mitkSlicedGeometry3D.h"
 #include "mitkAbstractTransformGeometry.h"
 #include "mitkApplyTransformMatrixOperation.h"
 #include "mitkInteractionConst.h"
 #include "mitkPlaneGeometry.h"
 #include "mitkPlaneOperation.h"
 #include "mitkRestorePlanePositionOperation.h"
 #include "mitkRotationOperation.h"
 #include "mitkSliceNavigationController.h"
 
 const mitk::ScalarType PI = 3.14159265359;
 
 mitk::SlicedGeometry3D::SlicedGeometry3D()
   : m_EvenlySpaced(true), m_Slices(0), m_ReferenceGeometry(nullptr), m_SliceNavigationController(nullptr)
 {
   m_DirectionVector.Fill(0);
   this->InitializeSlicedGeometry(m_Slices);
 }
 
 mitk::SlicedGeometry3D::SlicedGeometry3D(const SlicedGeometry3D &other)
   : Superclass(other),
     m_EvenlySpaced(other.m_EvenlySpaced),
     m_Slices(other.m_Slices),
     m_ReferenceGeometry(other.m_ReferenceGeometry),
     m_SliceNavigationController(other.m_SliceNavigationController)
 {
   m_DirectionVector.Fill(0);
   SetSpacing(other.GetSpacing());
   SetDirectionVector(other.GetDirectionVector());
 
   if (m_EvenlySpaced)
   {
     assert(!other.m_PlaneGeometries.empty() && "This may happen when you use one of the old Initialize methods, which had a bool parameter that is implicitly casted to the number of slices now.");
     PlaneGeometry::Pointer geometry = other.m_PlaneGeometries[0]->Clone();
     assert(geometry.IsNotNull());
     SetPlaneGeometry(geometry, 0);
   }
   else
   {
     unsigned int s;
     for (s = 0; s < other.m_Slices; ++s)
     {
       if (other.m_PlaneGeometries[s].IsNull())
       {
         assert(other.m_EvenlySpaced);
         m_PlaneGeometries[s] = nullptr;
       }
       else
       {
         PlaneGeometry *geometry2D = other.m_PlaneGeometries[s]->Clone();
         assert(geometry2D != nullptr);
         SetPlaneGeometry(geometry2D, s);
       }
     }
   }
 }
 
 mitk::SlicedGeometry3D::~SlicedGeometry3D()
 {
 }
 
 mitk::PlaneGeometry *mitk::SlicedGeometry3D::GetPlaneGeometry(int s) const
 {
   mitk::PlaneGeometry::Pointer geometry2D = nullptr;
 
   if (this->IsValidSlice(s))
   {
     geometry2D = m_PlaneGeometries[s];
 
     // If (a) m_EvenlySpaced==true, (b) we don't have a PlaneGeometry stored
     // for the requested slice, and (c) the first slice (s=0)
     // is a PlaneGeometry instance, then we calculate the geometry of the
     // requested as the plane of the first slice shifted by m_Spacing[2]*s
     // in the direction of m_DirectionVector.
     if ((m_EvenlySpaced) && (geometry2D.IsNull()))
     {
       PlaneGeometry *firstSlice = m_PlaneGeometries[0];
 
       if (firstSlice != nullptr &&
           dynamic_cast<AbstractTransformGeometry *>(m_PlaneGeometries[0].GetPointer()) == nullptr)
       {
         if ((m_DirectionVector[0] == 0.0) && (m_DirectionVector[1] == 0.0) && (m_DirectionVector[2] == 0.0))
         {
           m_DirectionVector = firstSlice->GetNormal();
           m_DirectionVector.Normalize();
         }
 
         Vector3D direction;
         direction = m_DirectionVector * this->GetSpacing()[2];
 
         mitk::PlaneGeometry::Pointer requestedslice;
         requestedslice = static_cast<mitk::PlaneGeometry *>(firstSlice->Clone().GetPointer());
 
         requestedslice->SetOrigin(requestedslice->GetOrigin() + direction * s);
 
         geometry2D = requestedslice;
         m_PlaneGeometries[s] = geometry2D;
       }
     }
     return geometry2D;
   }
   else
   {
     return nullptr;
   }
 }
 
 const mitk::BoundingBox *mitk::SlicedGeometry3D::GetBoundingBox() const
 {
   assert(this->IsBoundingBoxNull() == false);
   return Superclass::GetBoundingBox();
 }
 
 bool mitk::SlicedGeometry3D::SetPlaneGeometry(mitk::PlaneGeometry *geometry2D, int s)
 {
   if (this->IsValidSlice(s))
   {
     m_PlaneGeometries[s] = geometry2D;
     m_PlaneGeometries[s]->SetReferenceGeometry(m_ReferenceGeometry);
     return true;
   }
   return false;
 }
 
 void mitk::SlicedGeometry3D::InitializeSlicedGeometry(unsigned int slices)
 {
   Superclass::Initialize();
   m_Slices = slices;
 
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   Vector3D spacing;
   spacing.Fill(1.0);
   this->SetSpacing(spacing);
 
   m_DirectionVector.Fill(0);
 }
 
 void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, unsigned int slices)
 {
   assert(geometry2D != nullptr);
   this->InitializeEvenlySpaced(geometry2D, geometry2D->GetExtentInMM(2) / geometry2D->GetExtent(2), slices);
 }
 
 void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D,
                                                     mitk::ScalarType zSpacing,
                                                     unsigned int slices)
 {
   assert(geometry2D != nullptr);
   assert(geometry2D->GetExtent(0) > 0);
   assert(geometry2D->GetExtent(1) > 0);
 
   geometry2D->Register();
 
   Superclass::Initialize();
   m_Slices = slices;
 
   BoundingBox::BoundsArrayType bounds = geometry2D->GetBounds();
   bounds[4] = 0;
   bounds[5] = slices;
 
   // clear and reserve
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   Vector3D directionVector = geometry2D->GetAxisVector(2);
   directionVector.Normalize();
   directionVector *= zSpacing;
 
   // Normally we should use the following four lines to create a copy of
   // the transform contrained in geometry2D, because it may not be changed
   // by us. But we know that SetSpacing creates a new transform without
   // changing the old (coming from geometry2D), so we can use the fifth
   // line instead. We check this at (**).
   //
   // AffineTransform3D::Pointer transform = AffineTransform3D::New();
   // transform->SetMatrix(geometry2D->GetIndexToWorldTransform()->GetMatrix());
   // transform->SetOffset(geometry2D->GetIndexToWorldTransform()->GetOffset());
   // SetIndexToWorldTransform(transform);
 
   this->SetIndexToWorldTransform(geometry2D->GetIndexToWorldTransform());
 
   mitk::Vector3D spacing;
   FillVector3D(spacing, geometry2D->GetExtentInMM(0) / bounds[1], geometry2D->GetExtentInMM(1) / bounds[3], zSpacing);
 
   this->SetDirectionVector(directionVector);
   this->SetBounds(bounds);
   this->SetPlaneGeometry(geometry2D, 0);
   this->SetSpacing(spacing, true);
   this->SetEvenlySpaced();
 
   // this->SetTimeBounds( geometry2D->GetTimeBounds() );
 
   assert(this->GetIndexToWorldTransform() != geometry2D->GetIndexToWorldTransform()); // (**) see above.
 
   this->SetFrameOfReferenceID(geometry2D->GetFrameOfReferenceID());
   this->SetImageGeometry(geometry2D->GetImageGeometry());
 
   geometry2D->UnRegister();
 }
 
 void mitk::SlicedGeometry3D::InitializePlanes(const mitk::BaseGeometry *geometry3D,
                                               mitk::PlaneGeometry::PlaneOrientation planeorientation,
                                               bool top,
                                               bool frontside,
                                               bool rotated)
 {
   m_ReferenceGeometry = geometry3D;
 
   PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New();
   planeGeometry->InitializeStandardPlane(geometry3D, top, planeorientation, frontside, rotated);
 
   int worldAxis =
       planeorientation == PlaneGeometry::Sagittal ? 0 :
       planeorientation == PlaneGeometry::Frontal  ? 1 : 2;
 
   // Inspired by:
   // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
 
   mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
   matrix.GetVnlMatrix().normalize_columns();
-  mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
-
-  int dominantAxis = itk::Function::Max3(
-      inverseMatrix[0][worldAxis],
-      inverseMatrix[1][worldAxis],
-      inverseMatrix[2][worldAxis]);
+  mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetTranspose();
 
+  int dominantAxis = planeGeometry->CalculateDominantAxes(inverseMatrix).at(worldAxis);
   ScalarType viewSpacing = geometry3D->GetSpacing()[dominantAxis];
 
   /// Although the double value returned by GetExtent() holds a round number,
   /// you need to add 0.5 to safely convert it to unsigned it. I have seen a
   /// case when the result was less by one without this.
   auto slices = static_cast<unsigned int>(geometry3D->GetExtent(dominantAxis) + 0.5);
   if ( slices == 0 && geometry3D->GetExtent(dominantAxis) > 0) {
       // require at least one slice if there is _some_ extent
       slices = 1;
   }
 
 #ifndef NDEBUG
   int upDirection = itk::Function::Sign(inverseMatrix[dominantAxis][worldAxis]);
 
   /// The normal vector of an imaginary plane that points from the world origin (bottom left back
   /// corner or the world, with the lowest physical coordinates) towards the inside of the volume,
   /// along the renderer axis. Length is the slice thickness.
   Vector3D worldPlaneNormal = inverseMatrix.get_row(dominantAxis) * (upDirection * viewSpacing);
 
   /// The normal of the standard plane geometry just created.
   Vector3D standardPlaneNormal = planeGeometry->GetNormal();
 
   /// The standard plane must be parallel to the 'world plane'. The normal of the standard plane
   /// must point against the world plane if and only if 'top' is 'false'. The length of the
   /// standard plane normal must be equal to the slice thickness.
   assert((standardPlaneNormal - (top ? 1.0 : -1.0) * worldPlaneNormal).GetSquaredNorm() < 0.000001);
 #endif
 
   this->InitializeEvenlySpaced(planeGeometry, viewSpacing, slices);
 
 #ifndef NDEBUG
   /// The standard plane normal and the z axis vector of the sliced geometry must point in
   /// the same direction.
   Vector3D zAxisVector = this->GetAxisVector(2);
   Vector3D upscaledStandardPlaneNormal = standardPlaneNormal;
   upscaledStandardPlaneNormal *= slices;
   assert((zAxisVector - upscaledStandardPlaneNormal).GetSquaredNorm() < 0.000001);
 
   /// You can use this test is to check the handedness of the coordinate system of the current
   /// geometry. In principle, you can use either left- or right-handed coordinate systems, but
   /// you normally want it to be consistent, that is the handedness should be the same across
   /// the renderers of the same viewer.
 //  ScalarType det = vnl_det(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix());
 //  MITK_DEBUG << "world axis: " << worldAxis << (det > 0 ? " ; right-handed" : " ; left-handed");
 #endif
 }
 
 void mitk::SlicedGeometry3D::ReinitializePlanes(const Point3D &center, const Point3D &referencePoint)
 {
   // Need a reference frame to align the rotated planes
   if (!m_ReferenceGeometry)
   {
     return;
   }
 
   // Get first plane of plane stack
   PlaneGeometry *firstPlane = m_PlaneGeometries[0];
 
   // If plane stack is empty, exit
   if (!firstPlane || dynamic_cast<AbstractTransformGeometry *>(firstPlane))
   {
     return;
   }
 
   // Calculate the "directed" spacing when taking the plane (defined by its axes
   // vectors and normal) as the reference coordinate frame.
   //
   // This is done by calculating the radius of the ellipsoid defined by the
   // original volume spacing axes, in the direction of the respective axis of the
   // reference frame.
   mitk::Vector3D axis0 = firstPlane->GetAxisVector(0);
   mitk::Vector3D axis1 = firstPlane->GetAxisVector(1);
   mitk::Vector3D normal = firstPlane->GetNormal();
   normal.Normalize();
 
   Vector3D spacing;
   spacing[0] = this->CalculateSpacing(axis0);
   spacing[1] = this->CalculateSpacing(axis1);
   spacing[2] = this->CalculateSpacing(normal);
 
   Superclass::SetSpacing(spacing);
 
   // Now we need to calculate the number of slices in the plane's normal
   // direction, so that the entire volume is covered. This is done by first
   // calculating the dot product between the volume diagonal (the maximum
   // distance inside the volume) and the normal, and dividing this value by
   // the directed spacing calculated above.
   ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * normal[0]) +
                               std::abs(m_ReferenceGeometry->GetExtentInMM(1) * normal[1]) +
                               std::abs(m_ReferenceGeometry->GetExtentInMM(2) * normal[2]);
 
   if (directedExtent >= spacing[2])
   {
     m_Slices = static_cast<unsigned int>(directedExtent / spacing[2] + 0.5);
   }
   else
   {
     m_Slices = 1;
   }
 
   // The origin of our "first plane" needs to be adapted to this new extent.
   // To achieve this, we first calculate the current distance to the volume's
   // center, and then shift the origin in the direction of the normal by the
   // difference between this distance and half of the new extent.
   double centerOfRotationDistance = firstPlane->SignedDistanceFromPlane(center);
 
   if (centerOfRotationDistance > 0)
   {
     firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (centerOfRotationDistance - directedExtent / 2.0));
     m_DirectionVector = normal;
   }
   else
   {
     firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (directedExtent / 2.0 + centerOfRotationDistance));
     m_DirectionVector = -normal;
   }
 
   // Now we adjust this distance according with respect to the given reference
   // point: we need to make sure that the point is touched by one slice of the
   // new slice stack.
   double referencePointDistance = firstPlane->SignedDistanceFromPlane(referencePoint);
 
   auto referencePointSlice = static_cast<int>(referencePointDistance / spacing[2]);
 
   double alignmentValue = referencePointDistance / spacing[2] - referencePointSlice;
 
   firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * alignmentValue * spacing[2]);
 
   // Finally, we can clear the previous geometry stack and initialize it with
   // our re-initialized "first plane".
   m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr));
 
   if (m_Slices > 0)
   {
     m_PlaneGeometries[0] = firstPlane;
   }
 
   // Reinitialize SNC with new number of slices
   m_SliceNavigationController->GetSlice()->SetSteps(m_Slices);
 
   this->Modified();
 }
 
 double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &d) const
 {
   // Need the spacing of the underlying dataset / geometry
   if (!m_ReferenceGeometry)
   {
     return 1.0;
   }
 
   const mitk::Vector3D &spacing = m_ReferenceGeometry->GetSpacing();
   return SlicedGeometry3D::CalculateSpacing(spacing, d);
 }
 
 double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &spacing, const mitk::Vector3D &d)
 {
   // The following can be derived from the ellipsoid equation
   //
   //   1 = x^2/a^2 + y^2/b^2 + z^2/c^2
   //
   // where (a,b,c) = spacing of original volume (ellipsoid radii)
   // and   (x,y,z) = scaled coordinates of vector d (according to ellipsoid)
   //
   double scaling = d[0] * d[0] / (spacing[0] * spacing[0]) + d[1] * d[1] / (spacing[1] * spacing[1]) +
                    d[2] * d[2] / (spacing[2] * spacing[2]);
 
   scaling = sqrt(scaling);
 
   return (sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]) / scaling);
 }
 
 mitk::Vector3D mitk::SlicedGeometry3D::AdjustNormal(const mitk::Vector3D &normal) const
 {
   TransformType::Pointer inverse = TransformType::New();
   m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse(inverse);
 
   Vector3D transformedNormal = inverse->TransformVector(normal);
 
   transformedNormal.Normalize();
   return transformedNormal;
 }
 
 void mitk::SlicedGeometry3D::SetImageGeometry(const bool isAnImageGeometry)
 {
   Superclass::SetImageGeometry(isAnImageGeometry);
 
   unsigned int s;
   for (s = 0; s < m_Slices; ++s)
   {
     mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
     if (geometry != nullptr)
     {
       geometry->SetImageGeometry(isAnImageGeometry);
     }
   }
 }
 
 void mitk::SlicedGeometry3D::ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry)
 {
   unsigned int s;
   for (s = 0; s < m_Slices; ++s)
   {
     mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
     if (geometry != nullptr)
     {
       geometry->ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry);
     }
   }
 
   Superclass::ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry);
 }
 
 bool mitk::SlicedGeometry3D::IsValidSlice(int s) const
 {
   return ((s >= 0) && (s < (int)m_Slices));
 }
 
 const mitk::BaseGeometry *mitk::SlicedGeometry3D::GetReferenceGeometry() const
 {
   return m_ReferenceGeometry;
 }
 
 void mitk::SlicedGeometry3D::SetReferenceGeometry(const BaseGeometry *referenceGeometry)
 {
   m_ReferenceGeometry = referenceGeometry;
 
   std::vector<PlaneGeometry::Pointer>::iterator it;
 
   for (it = m_PlaneGeometries.begin(); it != m_PlaneGeometries.end(); ++it)
   {
     (*it)->SetReferenceGeometry(referenceGeometry);
   }
 }
 
 bool mitk::SlicedGeometry3D::HasReferenceGeometry() const
 {
   return ( m_ReferenceGeometry != nullptr );
 }
 
 void mitk::SlicedGeometry3D::PreSetSpacing(const mitk::Vector3D &aSpacing)
 {
   bool hasEvenlySpacedPlaneGeometry = false;
   mitk::Point3D origin;
   mitk::Vector3D rightDV, bottomDV;
   BoundingBox::BoundsArrayType bounds;
 
   // Check for valid spacing
   if (!(aSpacing[0] > 0 && aSpacing[1] > 0 && aSpacing[2] > 0))
   {
     mitkThrow() << "You try to set a spacing with at least one element equal or "
                    "smaller to \"0\". This might lead to a crash during rendering. Please double"
                    " check your data!";
   }
 
   // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
   // since the spacing influences them
   if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
   {
     const PlaneGeometry *planeGeometry = m_PlaneGeometries[0];
 
     if (planeGeometry && !dynamic_cast<const AbstractTransformGeometry *>(planeGeometry))
     {
       this->WorldToIndex(planeGeometry->GetOrigin(), origin);
       this->WorldToIndex(planeGeometry->GetAxisVector(0), rightDV);
       this->WorldToIndex(planeGeometry->GetAxisVector(1), bottomDV);
 
       bounds = planeGeometry->GetBounds();
       hasEvenlySpacedPlaneGeometry = true;
     }
   }
 
   BaseGeometry::_SetSpacing(aSpacing);
 
   mitk::PlaneGeometry::Pointer firstGeometry;
 
   // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
   // since the spacing influences them
   if (hasEvenlySpacedPlaneGeometry)
   {
     // create planeGeometry according to new spacing
     this->IndexToWorld(origin, origin);
     this->IndexToWorld(rightDV, rightDV);
     this->IndexToWorld(bottomDV, bottomDV);
 
     mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New();
     planeGeometry->SetImageGeometry(this->GetImageGeometry());
 
     planeGeometry->SetReferenceGeometry(m_ReferenceGeometry);
 
     // Store spacing, as Initialize... needs a pointer
     mitk::Vector3D lokalSpacing = this->GetSpacing();
     planeGeometry->InitializeStandardPlane(rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &lokalSpacing);
     planeGeometry->SetOrigin(origin);
     planeGeometry->SetBounds(bounds);
 
     firstGeometry = planeGeometry;
   }
   else if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
   {
     firstGeometry = m_PlaneGeometries[0].GetPointer();
   }
 
   // clear and reserve
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   if (m_Slices > 0)
   {
     m_PlaneGeometries[0] = firstGeometry;
   }
 
   this->Modified();
 }
 
 void mitk::SlicedGeometry3D::SetSliceNavigationController(SliceNavigationController *snc)
 {
   m_SliceNavigationController = snc;
 }
 
 mitk::SliceNavigationController *mitk::SlicedGeometry3D::GetSliceNavigationController()
 {
   return m_SliceNavigationController;
 }
 
 void mitk::SlicedGeometry3D::SetEvenlySpaced(bool on)
 {
   if (m_EvenlySpaced != on)
   {
     m_EvenlySpaced = on;
     this->Modified();
   }
 }
 
 void mitk::SlicedGeometry3D::SetDirectionVector(const mitk::Vector3D &directionVector)
 {
   Vector3D newDir = directionVector;
   newDir.Normalize();
   if (newDir != m_DirectionVector)
   {
     m_DirectionVector = newDir;
     this->Modified();
   }
 }
 
 // void
 // mitk::SlicedGeometry3D::SetTimeBounds( const mitk::TimeBounds& timebounds )
 //{
 //  Superclass::SetTimeBounds( timebounds );
 //
 //  unsigned int s;
 //  for ( s = 0; s < m_Slices; ++s )
 //  {
 //    if(m_Geometry2Ds[s].IsNotNull())
 //    {
 //      m_Geometry2Ds[s]->SetTimeBounds( timebounds );
 //    }
 //  }
 //  m_TimeBounds = timebounds;
 //}
 
 itk::LightObject::Pointer mitk::SlicedGeometry3D::InternalClone() const
 {
   Self::Pointer newGeometry = new SlicedGeometry3D(*this);
   newGeometry->UnRegister();
   return newGeometry.GetPointer();
 }
 
 void mitk::SlicedGeometry3D::PrintSelf(std::ostream &os, itk::Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
   os << indent << " EvenlySpaced: " << m_EvenlySpaced << std::endl;
   if (m_EvenlySpaced)
   {
     os << indent << " DirectionVector: " << m_DirectionVector << std::endl;
   }
   os << indent << " Slices: " << m_Slices << std::endl;
 
   os << std::endl;
   os << indent << " GetPlaneGeometry(0): ";
   if (this->GetPlaneGeometry(0) == nullptr)
   {
     os << "nullptr" << std::endl;
   }
   else
   {
     this->GetPlaneGeometry(0)->Print(os, indent);
   }
 }
 
 void mitk::SlicedGeometry3D::ExecuteOperation(Operation *operation)
 {
   PlaneGeometry::Pointer geometry2D;
   ApplyTransformMatrixOperation *applyMatrixOp;
   Point3D center;
 
   switch (operation->GetOperationType())
   {
     case OpNOTHING:
       break;
 
     case OpROTATE:
       if (m_EvenlySpaced)
       {
         // Need a reference frame to align the rotation
         if (m_ReferenceGeometry)
         {
           // Clear all generated geometries and then rotate only the first slice.
           // The other slices will be re-generated on demand
 
           // Save first slice
           PlaneGeometry::Pointer geometry2D = m_PlaneGeometries[0];
 
           auto *rotOp = dynamic_cast<RotationOperation *>(operation);
 
           // Generate a RotationOperation using the dataset center instead of
           // the supplied rotation center. This is necessary so that the rotated
           // zero-plane does not shift away. The supplied center is instead used
           // to adjust the slice stack afterwards.
           Point3D center = m_ReferenceGeometry->GetCenter();
 
           RotationOperation centeredRotation(
             rotOp->GetOperationType(), center, rotOp->GetVectorOfRotation(), rotOp->GetAngleOfRotation());
 
           // Rotate first slice
           geometry2D->ExecuteOperation(&centeredRotation);
 
           // Clear the slice stack and adjust it according to the center of
           // the dataset and the supplied rotation center (see documentation of
           // ReinitializePlanes)
           this->ReinitializePlanes(center, rotOp->GetCenterOfRotation());
 
           geometry2D->SetSpacing(this->GetSpacing());
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->SelectSliceByPoint(rotOp->GetCenterOfRotation());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
 
           BaseGeometry::ExecuteOperation(&centeredRotation);
         }
         else
         {
           // we also have to consider the case, that there is no reference geometry available.
           if (m_PlaneGeometries.size() > 0)
           {
             // Reach through to all slices in my container
             for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
             {
               // Test for empty slices, which can happen if evenly spaced geometry
               if ((*iter).IsNotNull())
               {
                 (*iter)->ExecuteOperation(operation);
               }
             }
 
             // rotate overall geometry
             auto *rotOp = dynamic_cast<RotationOperation *>(operation);
             BaseGeometry::ExecuteOperation(rotOp);
           }
         }
       }
       else
       {
         // Reach through to all slices
         for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
         {
           (*iter)->ExecuteOperation(operation);
         }
       }
       break;
 
     case OpORIENT:
       if (m_EvenlySpaced)
       {
         // get operation data
         auto *planeOp = dynamic_cast<PlaneOperation *>(operation);
 
         // Get first slice
         PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
 
         // Need a PlaneGeometry, a PlaneOperation and a reference frame to
         // carry out the re-orientation. If not all avaialble, stop here
         if (!m_ReferenceGeometry ||
             (!planeGeometry || dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer())) || !planeOp)
         {
           break;
         }
 
         // General Behavior:
         // Clear all generated geometries and then rotate only the first slice.
         // The other slices will be re-generated on demand
 
         //
         // 1st Step: Reorient Normal Vector of first plane
         //
         Point3D center = planeOp->GetPoint(); // m_ReferenceGeometry->GetCenter();
         mitk::Vector3D currentNormal = planeGeometry->GetNormal();
         mitk::Vector3D newNormal;
         if (planeOp->AreAxisDefined())
         {
           // If planeOp was defined by one centerpoint and two axis vectors
           newNormal = CrossProduct(planeOp->GetAxisVec0(), planeOp->GetAxisVec1());
         }
         else
         {
           // If planeOp was defined by one centerpoint and one normal vector
           newNormal = planeOp->GetNormal();
         }
 
         // Get Rotation axis und angle
         currentNormal.Normalize();
         newNormal.Normalize();
         ScalarType rotationAngle = angle(currentNormal.GetVnlVector(), newNormal.GetVnlVector());
 
         rotationAngle *= 180.0 / vnl_math::pi; // from rad to deg
         Vector3D rotationAxis = itk::CrossProduct(currentNormal, newNormal);
         if (std::abs(rotationAngle - 180) < mitk::eps)
         {
           // current Normal and desired normal are not linear independent!!(e.g 1,0,0 and -1,0,0).
           // Rotation Axis should be ANY vector that is 90� to current Normal
           mitk::Vector3D helpNormal;
           helpNormal = currentNormal;
           helpNormal[0] += 1;
           helpNormal[1] -= 1;
           helpNormal[2] += 1;
           helpNormal.Normalize();
           rotationAxis = itk::CrossProduct(helpNormal, currentNormal);
         }
 
         RotationOperation centeredRotation(mitk::OpROTATE, center, rotationAxis, rotationAngle);
 
         // Rotate first slice
         planeGeometry->ExecuteOperation(&centeredRotation);
 
         // Reinitialize planes and select slice, if my rotations are all done.
         if (!planeOp->AreAxisDefined())
         {
           // Clear the slice stack and adjust it according to the center of
           // rotation and plane position (see documentation of ReinitializePlanes)
           this->ReinitializePlanes(center, planeOp->GetPoint());
           planeGeometry->SetSpacing(this->GetSpacing());
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
         }
 
         // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
         BaseGeometry::ExecuteOperation(&centeredRotation);
 
         //
         // 2nd step. If axis vectors were defined, rotate the plane around its normal to fit these
         //
 
         if (planeOp->AreAxisDefined())
         {
           mitk::Vector3D vecAxixNew = planeOp->GetAxisVec0();
           vecAxixNew.Normalize();
           mitk::Vector3D VecAxisCurr = planeGeometry->GetAxisVector(0);
           VecAxisCurr.Normalize();
 
           ScalarType rotationAngle = angle(VecAxisCurr.GetVnlVector(), vecAxixNew.GetVnlVector());
           rotationAngle = rotationAngle * 180 / PI; // Rad to Deg
 
           // we rotate around the normal of the plane, but we do not know, if we need to rotate clockwise
           // or anti-clockwise. So we rotate around the crossproduct of old and new Axisvector.
           // Since both axis vectors lie in the plane, the crossproduct is the planes normal or the negative planes
           // normal
 
           rotationAxis = itk::CrossProduct(VecAxisCurr, vecAxixNew);
           if (std::abs(rotationAngle - 180) < mitk::eps)
           {
             // current axisVec and desired axisVec are not linear independent!!(e.g 1,0,0 and -1,0,0).
             // Rotation Axis can be just plane Normal. (have to rotate by 180�)
             rotationAxis = newNormal;
           }
 
           // Perfom Rotation
           mitk::RotationOperation op(mitk::OpROTATE, center, rotationAxis, rotationAngle);
           planeGeometry->ExecuteOperation(&op);
 
           // Apply changes on first slice to whole slice stack
           this->ReinitializePlanes(center, planeOp->GetPoint());
           planeGeometry->SetSpacing(this->GetSpacing());
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
 
           // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
           BaseGeometry::ExecuteOperation(&op);
         }
       }
       else
       {
         // Reach through to all slices
         for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
         {
           (*iter)->ExecuteOperation(operation);
         }
       }
       break;
 
     case OpRESTOREPLANEPOSITION:
       if (m_EvenlySpaced)
       {
         // Save first slice
         PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
 
         auto *restorePlaneOp = dynamic_cast<RestorePlanePositionOperation *>(operation);
 
         // Need a PlaneGeometry, a PlaneOperation and a reference frame to
         // carry out the re-orientation
         if (m_ReferenceGeometry &&
             (planeGeometry && dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer()) == nullptr) &&
             restorePlaneOp)
         {
           // Clear all generated geometries and then rotate only the first slice.
           // The other slices will be re-generated on demand
 
           // Rotate first slice
           planeGeometry->ExecuteOperation(restorePlaneOp);
 
           m_DirectionVector = restorePlaneOp->GetDirectionVector();
 
           double centerOfRotationDistance = planeGeometry->SignedDistanceFromPlane(m_ReferenceGeometry->GetCenter());
 
           if (centerOfRotationDistance <= 0)
           {
             m_DirectionVector = -m_DirectionVector;
           }
 
           Vector3D spacing = restorePlaneOp->GetSpacing();
 
           Superclass::SetSpacing(spacing);
 
           // /*Now we need to calculate the number of slices in the plane's normal
           // direction, so that the entire volume is covered. This is done by first
           // calculating the dot product between the volume diagonal (the maximum
           // distance inside the volume) and the normal, and dividing this value by
           // the directed spacing calculated above.*/
           ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * m_DirectionVector[0]) +
                                       std::abs(m_ReferenceGeometry->GetExtentInMM(1) * m_DirectionVector[1]) +
                                       std::abs(m_ReferenceGeometry->GetExtentInMM(2) * m_DirectionVector[2]);
 
           if (directedExtent >= spacing[2])
           {
             m_Slices = static_cast<unsigned int>(directedExtent / spacing[2] + 0.5);
           }
           else
           {
             m_Slices = 1;
           }
 
           m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr));
 
           if (m_Slices > 0)
           {
             m_PlaneGeometries[0] = planeGeometry;
           }
 
           m_SliceNavigationController->GetSlice()->SetSteps(m_Slices);
 
           this->Modified();
 
           // End Reinitialization
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->GetSlice()->SetPos(restorePlaneOp->GetPos());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
           BaseGeometry::ExecuteOperation(restorePlaneOp);
         }
       }
       else
       {
         // Reach through to all slices
         for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
         {
           (*iter)->ExecuteOperation(operation);
         }
       }
       break;
 
     case OpAPPLYTRANSFORMMATRIX:
 
       // Clear all generated geometries and then transform only the first slice.
       // The other slices will be re-generated on demand
 
       // Save first slice
       geometry2D = m_PlaneGeometries[0];
 
       applyMatrixOp = dynamic_cast<ApplyTransformMatrixOperation *>(operation);
 
       // Apply transformation to first plane
       geometry2D->ExecuteOperation(applyMatrixOp);
 
       // Generate a ApplyTransformMatrixOperation using the dataset center instead of
       // the supplied rotation center. The supplied center is instead used to adjust the
       // slice stack afterwards (see OpROTATE).
       center = m_ReferenceGeometry->GetCenter();
 
       // Clear the slice stack and adjust it according to the center of
       // the dataset and the supplied rotation center (see documentation of
       // ReinitializePlanes)
       this->ReinitializePlanes(center, applyMatrixOp->GetReferencePoint());
 
       BaseGeometry::ExecuteOperation(applyMatrixOp);
       break;
 
     default: // let handle by base class if we don't do anything
       BaseGeometry::ExecuteOperation(operation);
   }
 
   this->Modified();
 }
diff --git a/Modules/Core/test/mitkPlaneGeometryTest.cpp b/Modules/Core/test/mitkPlaneGeometryTest.cpp
index ba4a70cbff..10ea8ad159 100644
--- a/Modules/Core/test/mitkPlaneGeometryTest.cpp
+++ b/Modules/Core/test/mitkPlaneGeometryTest.cpp
@@ -1,1078 +1,1097 @@
 /*===================================================================
 
 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 "mitkAffineTransform3D.h"
 #include "mitkBaseGeometry.h"
 #include "mitkGeometry3D.h"
 #include "mitkInteractionConst.h"
 #include "mitkLine.h"
 #include "mitkPlaneGeometry.h"
 #include "mitkRotationOperation.h"
 #include "mitkSlicedGeometry3D.h"
 #include "mitkThinPlateSplineCurvedGeometry.h"
 
 #include <mitkTestFixture.h>
 #include <mitkTestingMacros.h>
 
 #include <vnl/vnl_quaternion.h>
 #include <vnl/vnl_quaternion.hxx>
 
 #include <fstream>
 #include <iomanip>
+#include <mitkIOUtil.h>
+#include <mitkImage.h>
 
 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);
+  MITK_TEST(TestDominantAxesError);
+  MITK_TEST(TestCheckRotationMatrix);
 
   // 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, spacing;
   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);
     mitk::FillVector3D(spacing, 1.0, 1.0, thicknessInMM);
 
     planegeometry->InitializeStandardPlane(right, bottom);
     planegeometry->SetOrigin(origin);
     planegeometry->SetSpacing(spacing);
   }
 
   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());
 
     g3d = mitk::Geometry3D::New();
     base = dynamic_cast<mitk::BaseGeometry *>(g3d.GetPointer());
     CPPUNIT_ASSERT_MESSAGE("Geometry3D should be castable to BaseGeometry", base.IsNotNull());
 
     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());
 
     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 TestDominantAxesError()
+  {
+    auto image = mitk::IOUtil::Load<mitk::Image>(GetTestDataFilePath("NotQuiteARotationMatrix.nrrd"));
+    auto matrix = image->GetGeometry()->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().transpose();
+    std::vector< int > axes = mitk::PlaneGeometry::CalculateDominantAxes(matrix);
+    CPPUNIT_ASSERT_MESSAGE("Domiant axes cannot be determined in this dataset. Output should be default ordering.", axes.at(0)==0 && axes.at(1)==1 && axes.at(2)==2);
+  }
+
+  void TestCheckRotationMatrix()
+  {
+    auto image = mitk::IOUtil::Load<mitk::Image>(GetTestDataFilePath("NotQuiteARotationMatrix.nrrd"));
+    bool is_rotation = mitk::PlaneGeometry::CheckRotationMatrix(image->GetGeometry()->GetIndexToWorldTransform());
+    CPPUNIT_ASSERT_MESSAGE("Since the test data matrix is not quite a rotation matrix, this should be detected.", !is_rotation);
+  }
+
   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);
 
     mitk::FillVector3D(down, 1.4012984643248170709237295832899161312802619418765e-45, 0, 0);
     mitk::FillVector3D(spacing,
                        0,
                        1.4713633875410579244699160624544119378442750389703e-43,
                        9.2806360452222355258639080851310540729807238879469e-32);
 
     std::cout << "Testing InitializeStandardPlane(rightVector, downVector, spacing = nullptr): " << 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.
    *
    */
   // 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.
    */
   // 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()));
   }
 
   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));
 
     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));
 
     ;
   }
 
   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));
 
     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 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));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ",
       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));
 
     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 GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ",
       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));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing GetAxisVector(direction).GetNorm() != planegeometry->GetExtentInMM(direction) of rotated version: ",
       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 = planegeometry->Clone();
 
     // 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));
     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));
     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));
 
     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)); // T22254: Flipped sign
 
     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)); // T22254: Flipped sign
 
     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)); // T22254: Flipped sign
 
     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)); // T22254: Flipped sign
     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 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));
     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));
     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));
 
     CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of axially initialized version: ",
                            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));
     CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of axially initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(2), thicknessInMM, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(0), right, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(1), bottom, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of axially initialized version: ",
                            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
    **/
   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)