diff --git a/Modules/Core/include/mitkBaseGeometry.h b/Modules/Core/include/mitkBaseGeometry.h
index 833189965d..4d00557d8d 100644
--- a/Modules/Core/include/mitkBaseGeometry.h
+++ b/Modules/Core/include/mitkBaseGeometry.h
@@ -1,752 +1,763 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #ifndef mitkBaseGeometry_h
 #define mitkBaseGeometry_h
 
 #include "mitkOperationActor.h"
 #include <MitkCoreExports.h>
 #include <mitkCommon.h>
 
 #include "itkScalableAffineTransform.h"
 #include "mitkNumericTypes.h"
 #include <itkBoundingBox.h>
 #include <itkIndex.h>
 #include <itkQuaternionRigidTransform.h>
 #include <mitkAffineTransform3D.h>
 
 #include <mitkGeometryTransformHolder.h>
 #include <vtkTransform.h>
 
 class vtkMatrix4x4;
 class vtkMatrixToLinearTransform;
 class vtkLinearTransform;
 
 namespace mitk
 {
   //##Documentation
   //## @brief Standard 3D-BoundingBox typedef
   //##
   //## Standard 3D-BoundingBox typedef to get rid of template arguments (3D, type).
   typedef itk::BoundingBox<unsigned long, 3, ScalarType> BoundingBox;
 
   //##Documentation
   //## @brief Standard typedef for time-bounds
   typedef itk::FixedArray<ScalarType, 2> TimeBounds;
   typedef itk::FixedArray<ScalarType, 3> FixedArrayType;
 
   //##Documentation
   //## @brief BaseGeometry Describes the geometry of a data object
   //##
   //## The class holds
   //## \li a bounding box which is axes-parallel in intrinsic coordinates
   //## (often integer indices of pixels), to be accessed by
   //## GetBoundingBox()
   //## \li a transform to convert intrinsic coordinates into a
   //## world-coordinate system with coordinates in millimeters
   //## and milliseconds (all are floating point values), to
   //## be accessed by GetIndexToWorldTransform()
   //## \li an origin and spacing to define the geometry
   //##
   //## BaseGeometry and its sub-classes allow converting between
   //## intrinsic coordinates (called index or unit coordinates)
   //## and world-coordinates (called world or mm coordinates),
   //## e.g. WorldToIndex.
   //## In case you need integer index coordinates, provide an
   //## mitk::Index3D (or itk::Index) as target variable to
   //## WorldToIndex, otherwise you will get a continuous index
   //## (floating point values).
   //##
   //## An important sub-class is SlicedGeometry3D, which describes
   //## data objects consisting of slices, e.g., objects of type Image.
   //## Conversions between world coordinates (in mm) and unit coordinates
   //## (e.g., pixels in the case of an Image) can be performed.
   //##
   //## For more information on related classes, see \ref Geometry.
   //##
   //## BaseGeometry instances referring to an Image need a slightly
   //## different definition of corners, see SetImageGeometry. This
   //## is usually automatically called by Image.
   //##
   //## BaseGeometry have to be initialized in the method GenerateOutputInformation()
   //## of BaseProcess (or CopyInformation/ UpdateOutputInformation of BaseData,
   //## if possible, e.g., by analyzing pic tags in Image) subclasses. See also
   //## itk::ProcessObject::GenerateOutputInformation(),
   //## itk::DataObject::CopyInformation() and
   //## itk::DataObject::UpdateOutputInformation().
   //##
   //## At least, it can return the bounding box of the data object.
   //##
   //## The BaseGeometry class is an abstract class. The most simple implementation
   //## is the subclass Geometry3D.
   //##
   //## Rule: everything is in mm (ms) if not stated otherwise.
   //## @ingroup Geometry
   class MITKCORE_EXPORT BaseGeometry : public itk::Object, public OperationActor
   {
   public:
     mitkClassMacroItkParent(BaseGeometry, itk::Object);
     itkCloneMacro(Self);
 
       // ********************************** TypeDef **********************************
 
       typedef GeometryTransformHolder::TransformType TransformType;
     typedef itk::BoundingBox<unsigned long, 3, ScalarType> BoundingBoxType;
     typedef BoundingBoxType::BoundsArrayType BoundsArrayType;
     typedef BoundingBoxType::Pointer BoundingBoxPointer;
 
     // ********************************** Origin, Spacing **********************************
 
     //##Documentation
     //## @brief Get the origin, e.g. the upper-left corner of the plane
     const Point3D GetOrigin() const;
 
     //##Documentation
     //## @brief Set the origin, i.e. the upper-left corner of the plane
     //##
     void SetOrigin(const Point3D &origin);
 
     //##Documentation
     //## @brief Get the spacing (size of a pixel).
     //##
     const mitk::Vector3D GetSpacing() const;
 
     //##Documentation
     //## @brief Set the spacing (m_Spacing).
     //##
     //##The spacing is also changed in the IndexToWorldTransform.
     void SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing = false);
 
     //##Documentation
     //## @brief Get the origin as VnlVector
     //##
     //## \sa GetOrigin
     VnlVector GetOriginVnl() const;
 
     // ********************************** other functions **********************************
 
     //##Documentation
     //## @brief Get the DICOM FrameOfReferenceID referring to the
     //## used world coordinate system
     itkGetConstMacro(FrameOfReferenceID, unsigned int);
     //##Documentation
     //## @brief Set the DICOM FrameOfReferenceID referring to the
     //## used world coordinate system
     itkSetMacro(FrameOfReferenceID, unsigned int);
 
     itkGetConstMacro(IndexToWorldTransformLastModified, unsigned long);
 
     //##Documentation
     //## @brief Overload of function Modified() to prohibit several calls of Modified() using the ModifiedLock class.
     //##
     //## For the use of Modified(), see class ModifiedLock.
     void Modified() const override;
 
     friend class ModifiedLock;
 
     //##Documentation
     //## @brief Is this BaseGeometry in a state that is valid?
     //##
     //## This function returns always true in the BaseGeometry class. Other implementations are possible in subclasses.
     virtual bool IsValid() const;
 
     // ********************************** Initialize **********************************
 
     //##Documentation
     //## @brief Initialize the BaseGeometry
     void Initialize();
 
     void InitializeGeometry(Self *newGeometry) const;
 
     // ********************************** Transformations Set/Get **********************************
 
     //##Documentation
     //## @brief Get the transformation used to convert from index
     //## to world coordinates
     mitk::AffineTransform3D *GetIndexToWorldTransform();
 
     //##Documentation
     //## @brief Get the transformation used to convert from index
     //## to world coordinates
     const mitk::AffineTransform3D *GetIndexToWorldTransform() const;
 
     //## @brief Set the transformation used to convert from index
     //## to world coordinates. The spacing of the new transform is
     //## copied to m_spacing.
     void SetIndexToWorldTransform(mitk::AffineTransform3D *transform);
 
     //##Documentation
     //## @brief Convenience method for setting the ITK transform
     //## (m_IndexToWorldTransform) via an vtkMatrix4x4.The spacing of
     //## the new transform is copied to m_spacing.
     //## \sa SetIndexToWorldTransform
     void SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4 *vtkmatrix);
 
     //## @brief Set the transformation used to convert from index
     //## to world coordinates.This function keeps the original spacing.
     void SetIndexToWorldTransformWithoutChangingSpacing(mitk::AffineTransform3D *transform);
 
     //##Documentation
     //## @brief Convenience method for setting the ITK transform
     //## (m_IndexToWorldTransform) via an vtkMatrix4x4. This function keeps the original spacing.
     //## \sa SetIndexToWorldTransform
     void SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkMatrix4x4 *vtkmatrix);
 
     //## Get the Vtk Matrix which describes the transform.
     vtkMatrix4x4 *GetVtkMatrix();
 
     //##Documentation
     //## @brief Get the m_IndexToWorldTransform as a vtkLinearTransform
     vtkLinearTransform *GetVtkTransform() const;
 
     //##Documentation
     //## @brief Set the transform to identity, the spacing to 1 and origin to 0
     //##
     void SetIdentity();
 
     // ********************************** Transformations **********************************
 
     //##Documentation
     //## @brief Compose new IndexToWorldTransform with a given transform.
     //##
     //## This method composes m_IndexToWorldTransform with another transform,
     //## modifying self to be the composition of self and other.
     //## If the argument pre is true, then other is precomposed with self;
     //## that is, the resulting transformation consists of first applying
     //## other to the source, followed by self. If pre is false or omitted,
     //## then other is post-composed with self; that is the resulting
     //## transformation consists of first applying self to the source,
     //## followed by other.
     //## This method also changes m_spacing.
     void Compose(const TransformType *other, bool pre = false);
 
     //##Documentation
     //## @brief Compose new IndexToWorldTransform with a given vtkMatrix4x4.
     //##
     //## Converts the vtkMatrix4x4 into a itk-transform and calls the previous method.
     void Compose(const vtkMatrix4x4 *vtkmatrix, bool pre = false);
 
     //##Documentation
     //## @brief Translate the origin by a vector
     //##
     void Translate(const Vector3D &vector);
 
     //##Documentation
     //##@brief executes affine operations (translate, rotate, scale)
     void ExecuteOperation(Operation *operation) override;
 
     //##Documentation
     //## @brief Convert world coordinates (in mm) of a \em point to (continuous!) index coordinates
     //## \warning If you need (discrete) integer index coordinates (e.g., for iterating easily over an image),
     //## use WorldToIndex(const mitk::Point3D& pt_mm, itk::Index<VIndexDimension> &index).
     //## For further information about coordinates types, please see the Geometry documentation
     void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const;
 
     //##Documentation
     //## @brief Convert world coordinates (in mm) of a \em vector
     //## \a vec_mm to (continuous!) index coordinates.
     //## For further information about coordinates types, please see the Geometry documentation
     void WorldToIndex(const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const;
 
     //##Documentation
     //## @brief Convert world coordinates (in mm) of a \em point to (discrete!) index coordinates.
     //## This method rounds to integer indices!
     //## For further information about coordinates types, please see the Geometry documentation
     template <unsigned int VIndexDimension>
     void WorldToIndex(const mitk::Point3D &pt_mm, itk::Index<VIndexDimension> &index) const
     {
       typedef itk::Index<VIndexDimension> IndexType;
       mitk::Point3D pt_units;
       this->WorldToIndex(pt_mm, pt_units);
       int i, dim = index.GetIndexDimension();
       if (dim > 3)
       {
         index.Fill(0);
         dim = 3;
       }
       for (i = 0; i < dim; ++i)
       {
         index[i] = itk::Math::RoundHalfIntegerUp<typename IndexType::IndexValueType>(pt_units[i]);
       }
     }
 
     //##Documentation
     //## @brief Convert (continuous or discrete) index coordinates of a \em vector
     //## \a vec_units to world coordinates (in mm)
     //## For further information about coordinates types, please see the Geometry documentation
     void IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const;
 
     //##Documentation
     //## @brief Convert (continuous or discrete) index coordinates of a \em point to world coordinates (in mm)
     //## For further information about coordinates types, please see the Geometry documentation
     void IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const;
 
     //##Documentation
     //## @brief Convert (discrete) index coordinates of a \em point to world coordinates (in mm)
     //## For further information about coordinates types, please see the Geometry documentation
     template <unsigned int VIndexDimension>
     void IndexToWorld(const itk::Index<VIndexDimension> &index, mitk::Point3D &pt_mm) const
     {
       mitk::Point3D pt_units;
       pt_units.Fill(0);
       int i, dim = index.GetIndexDimension();
       if (dim > 3)
       {
         dim = 3;
       }
       for (i = 0; i < dim; ++i)
       {
         pt_units[i] = index[i];
       }
 
       IndexToWorld(pt_units, pt_mm);
     }
 
     //##Documentation
     //## @brief Convert (continuous or discrete) index coordinates of a \em vector
     //## \a vec_units to world coordinates (in mm)
     //## @deprecated First parameter (Point3D) is not used. If possible, please use void IndexToWorld(const
     // mitk::Vector3D& vec_units, mitk::Vector3D& vec_mm) const.
     //## For further information about coordinates types, please see the Geometry documentation
     void IndexToWorld(const mitk::Point3D &atPt3d_units, const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const;
 
     //##Documentation
     //## @brief Convert world coordinates (in mm) of a \em vector
     //## \a vec_mm to (continuous!) index coordinates.
     //## @deprecated First parameter (Point3D) is not used. If possible, please use void WorldToIndex(const
     // mitk::Vector3D& vec_mm, mitk::Vector3D& vec_units) const.
     //## For further information about coordinates types, please see the Geometry documentation
     void WorldToIndex(const mitk::Point3D &atPt3d_mm, const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const;
 
     //##Documentation
     //## @brief Deprecated for use with ITK version 3.10 or newer.
     //## Convert ITK physical coordinates of a \em point (in mm,
     //## but without a rotation) into MITK world coordinates (in mm)
     //##
     //## For more information, see WorldToItkPhysicalPoint.
     template <class TCoordRep>
     void ItkPhysicalPointToWorld(const itk::Point<TCoordRep, 3> &itkPhysicalPoint, mitk::Point3D &pt_mm) const
     {
       mitk::vtk2itk(itkPhysicalPoint, pt_mm);
     }
 
     //##Documentation
     //## @brief Deprecated for use with ITK version 3.10 or newer.
     //## Convert world coordinates (in mm) of a \em point to
     //## ITK physical coordinates (in mm, but without a possible rotation)
     //##
     //## This method is useful if you have want to access an mitk::Image
     //## via an itk::Image. ITK v3.8 and older did not support rotated (tilted)
     //## images, i.e., ITK images are always parallel to the coordinate axes.
     //## When accessing a (possibly rotated) mitk::Image via an itk::Image
     //## the rotational part of the transformation in the BaseGeometry is
     //## simply discarded; in other word: only the origin and spacing is
     //## used by ITK, not the complete matrix available in MITK.
     //## With WorldToItkPhysicalPoint you can convert an MITK world
     //## coordinate (including the rotation) into a coordinate that
     //## can be used with the ITK image as a ITK physical coordinate
     //## (excluding the rotation).
     template <class TCoordRep>
     void WorldToItkPhysicalPoint(const mitk::Point3D &pt_mm, itk::Point<TCoordRep, 3> &itkPhysicalPoint) const
     {
       mitk::vtk2itk(pt_mm, itkPhysicalPoint);
     }
 
     // ********************************** BoundingBox **********************************
 
     /** Get the bounding box */
     itkGetConstObjectMacro(BoundingBox, BoundingBoxType);
 
 // a bit of a misuse, but we want only doxygen to see the following:
 #ifdef DOXYGEN_SKIP
     //##Documentation
     //## @brief Get bounding box (in index/unit coordinates)
     itkGetConstObjectMacro(BoundingBox, BoundingBoxType);
     //##Documentation
     //## @brief Get bounding box (in index/unit coordinates) as a BoundsArrayType
     const BoundsArrayType GetBounds() const;
 #endif
     const BoundsArrayType GetBounds() const;
 
     //##Documentation
     //## \brief Set the bounding box (in index/unit coordinates)
     //##
     //## Only possible via the BoundsArray to make clear that a
     //## copy of the bounding-box is stored, not a reference to it.
     void SetBounds(const BoundsArrayType &bounds);
 
     //##Documentation
     //## @brief Set the bounding box (in index/unit coordinates) via a float array
     void SetFloatBounds(const float bounds[6]);
     //##Documentation
     //## @brief Set the bounding box (in index/unit coordinates) via a double array
     void SetFloatBounds(const double bounds[6]);
 
     //##Documentation
     //## @brief Get a VnlVector along bounding-box in the specified
     //## @a direction, length is spacing
     //##
     //## \sa GetAxisVector
     VnlVector GetMatrixColumn(unsigned int direction) const;
 
     //##Documentation
     //## @brief Calculates a bounding-box around the geometry relative
     //## to a coordinate system defined by a transform
     //##
     mitk::BoundingBox::Pointer CalculateBoundingBoxRelativeToTransform(const mitk::AffineTransform3D *transform) const;
 
 //##Documentation
 //## @brief Set the time bounds (in ms)
 // void SetTimeBounds(const TimeBounds& timebounds);
 
 // ********************************** Geometry **********************************
 
 #ifdef DOXYGEN_SKIP
     //##Documentation
     //## @brief Get the extent of the bounding box (in index/unit coordinates)
     //##
     //## To access the extent in mm use GetExtentInMM
     ScalarType GetExtent(unsigned int direction) const;
 #endif
 
     /** Get the extent of the bounding box */
     ScalarType GetExtent(unsigned int direction) const;
 
     //##Documentation
     //## @brief Get the extent of the bounding-box in the specified @a direction in mm
     //##
     //## Equals length of GetAxisVector(direction).
     ScalarType GetExtentInMM(int direction) const;
 
     //##Documentation
     //## @brief Get vector along bounding-box in the specified @a direction in mm
     //##
     //## The length of the vector is the size of the bounding-box in the
     //## specified @a direction in mm
     //## \sa GetMatrixColumn
     Vector3D GetAxisVector(unsigned int direction) const;
 
     //##Documentation
     //## @brief Checks, if the given geometry can be converted to 2D without information loss
     //## e.g. when a 2D image is saved, the matrix is usually cropped to 2x2, and when you load it back to MITK
     //## it will be filled with standard values. This function checks, if information would be lost during this
     //## procedure
     virtual bool Is2DConvertable();
 
     //##Documentation
     //## @brief Get the center of the bounding-box in mm
     //##
     Point3D GetCenter() const;
 
     //##Documentation
     //## @brief Get the squared length of the diagonal of the bounding-box in mm
     //##
     double GetDiagonalLength2() const;
 
     //##Documentation
     //## @brief Get the length of the diagonal of the bounding-box in mm
     //##
     double GetDiagonalLength() const;
 
     //##Documentation
     //## @brief Get the position of the corner number \a id (in world coordinates)
     //##
     //## See SetImageGeometry for how a corner is defined on images.
     Point3D GetCornerPoint(int id) const;
 
     //##Documentation
     //## @brief Get the position of a corner (in world coordinates)
     //##
     //## See SetImageGeometry for how a corner is defined on images.
     Point3D GetCornerPoint(bool xFront = true, bool yFront = true, bool zFront = true) const;
 
     //##Documentation
     //## @brief Set the extent of the bounding-box in the specified @a direction in mm
     //##
     //## @note This changes the matrix in the transform, @a not the bounds, which are given in units!
     void SetExtentInMM(int direction, ScalarType extentInMM);
 
     //##Documentation
     //## @brief Test whether the point \a p (world coordinates in mm) is
     //## inside the bounding box
     bool IsInside(const mitk::Point3D &p) const;
 
     //##Documentation
     //## @brief Test whether the point \a p ((continuous!)index coordinates in units) is
     //## inside the bounding box
     bool IsIndexInside(const mitk::Point3D &index) const;
 
     //##Documentation
     //## @brief Convenience method for working with ITK indices
     template <unsigned int VIndexDimension>
     bool IsIndexInside(const itk::Index<VIndexDimension> &index) const
     {
       int i, dim = index.GetIndexDimension();
       Point3D pt_index;
       pt_index.Fill(0);
       for (i = 0; i < dim; ++i)
       {
         pt_index[i] = index[i];
       }
       return IsIndexInside(pt_index);
     }
 
     // ********************************* Image Geometry ********************************
     //##Documentation
     //## @brief When switching from an Image Geometry to a normal Geometry (and the other way around), you have to
     //change
     // the origin as well (See Geometry Documentation)! This function will change the "isImageGeometry" bool flag and
     // changes the origin respectively.
     virtual void ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry);
 
     //##Documentation
     //## @brief Is this an ImageGeometry?
     //##
     //## For more information, see SetImageGeometry
     itkGetConstMacro(ImageGeometry, bool)
       //##Documentation
       //## @brief Define that this BaseGeometry is referring to an Image
       //##
       //## A geometry referring to an Image needs a slightly different
       //## definition of the position of the corners (see GetCornerPoint).
       //## The position of a voxel is defined by the position of its center.
       //## If we would use the origin (position of the (center of) the first
       //## voxel) as a corner and display this point, it would seem to be
       //## \em not at the corner but a bit within the image. Even worse for
       //## the opposite corner of the image: here the corner would appear
       //## outside the image (by half of the voxel diameter). Thus, we have
       //## to correct for this and to be able to do that, we need to know
       //## that the BaseGeometry is referring to an Image.
       itkSetMacro(ImageGeometry, bool);
       itkBooleanMacro(ImageGeometry);
 
         const GeometryTransformHolder *GetGeometryTransformHolder() const;
 
+    //##Documentation
+    //## @brief One to one mapping of axes to world orientations.
+    //##
+    //## The result is stored in the output argument that must be an array of three int values.
+    //## The elements of the array will be the axis indices that correspond to the sagittal,
+    //## coronal and axial orientations, in this order. It is guaranteed that each axis will
+    //## be mapped to different orientations.
+    //##
+    //## @param axes Output argument that will store the axis indices for each orientation.
+    void MapAxesToOrientations(int axes[]) const;
+
   protected:
     // ********************************** Constructor **********************************
     BaseGeometry();
     BaseGeometry(const BaseGeometry &other);
     ~BaseGeometry() override;
 
     itk::LightObject::Pointer InternalClone() const override = 0;
 
     void PrintSelf(std::ostream &os, itk::Indent indent) const override;
 
     static const std::string GetTransformAsString(TransformType *transformType);
 
     itkGetConstMacro(NDimensions, unsigned int);
 
       bool IsBoundingBoxNull() const;
 
     bool IsIndexToWorldTransformNull() const;
 
     void SetVtkMatrixDeepCopy(vtkTransform *vtktransform);
 
     void _SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing = false);
 
     //##Documentation
     //## @brief PreSetSpacing
     //##
     //## These virtual function allows a different beahiour in subclasses.
     //## Do implement them in every subclass of BaseGeometry. If not needed, use
     //## {Superclass::PreSetSpacing();};
     virtual void PreSetSpacing(const mitk::Vector3D & /*aSpacing*/){};
 
     //##Documentation
     //## @brief CheckBounds
     //##
     //## This function is called in SetBounds. Assertions can be implemented in this function (see PlaneGeometry.cpp).
     //## If you implement this function in a subclass, make sure, that all classes were your class inherits from
     //## have an implementation of CheckBounds
     //## (e.g. inheritance BaseGeometry <- A <- B. Implementation of CheckBounds in class B needs implementation in A as
     // well!)
     virtual void CheckBounds(const BoundsArrayType & /*bounds*/){};
 
     //##Documentation
     //## @brief CheckIndexToWorldTransform
     //##
     //## This function is called in SetIndexToWorldTransform. Assertions can be implemented in this function (see
     // PlaneGeometry.cpp).
     //## In Subclasses of BaseGeometry, implement own conditions or call Superclass::CheckBounds(bounds);.
     virtual void CheckIndexToWorldTransform(mitk::AffineTransform3D * /*transform*/){};
 
   private:
     GeometryTransformHolder *m_GeometryTransform;
 
     void InitializeGeometryTransformHolder(const BaseGeometry *otherGeometry);
 
     //##Documentation
     //## @brief Bounding Box, which is axes-parallel in intrinsic coordinates
     //## (often integer indices of pixels)
     BoundingBoxPointer m_BoundingBox;
 
     unsigned int m_FrameOfReferenceID;
 
     // mitk::TimeBounds m_TimeBounds;
 
     static const unsigned int m_NDimensions = 3;
 
     mutable TransformType::Pointer m_InvertedTransform;
 
     mutable unsigned long m_IndexToWorldTransformLastModified;
 
     bool m_ImageGeometry;
 
     //##Documentation
     //## @brief ModifiedLockFlag is used to prohibit the call of Modified()
     //##
     //## For the use of this Flag, see class ModifiedLock. This flag should only be set
     //## by the ModifiedLock class!
     bool m_ModifiedLockFlag;
 
     //##Documentation
     //## @brief ModifiedcalledFlag is used to collect calls of Modified().
     //##
     //## For the use of this Flag, see class ModifiedLock. This flag should only be set
     //## by the Modified() function!
     mutable bool m_ModifiedCalledFlag;
   };
 
   // ********************************** Equal Functions **********************************
   //
   // Static compare functions mainly for testing
   //
   /**
   * @brief Equal A function comparing two geometries for being identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the spacing, origin, axisvectors, extents, the matrix of the
   * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag.
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparison.
   * If you want to use different tolerance values for different parts of the geometry, feel free to use
   * the other comparison methods and write your own implementation of Equal.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param coordinateEps Tolerance for comparison of all spatial aspects (spacing, origin and grid alignment).
   * You can use mitk::eps in most cases.
   * @param directionEps Tolerance for comparison of all directional aspects (axis). You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry& leftHandSide,
     const mitk::BaseGeometry& rightHandSide,
     ScalarType coordinateEps,
     ScalarType directionEps,
     bool verbose = false);
 
   /**
   * @brief Equal A function comparing two geometries for being identical.
   *
   * @ingroup MITKTestingAPI
   *
   * This is an overloaded version that uses a single tolerance for spatial and directional aspects. For more details,
   * see the other overloaded version.
   *
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry &leftHandSide,
                              const mitk::BaseGeometry &rightHandSide,
                              ScalarType eps = mitk::eps,
                              bool verbose = false);
 
   /**
   * @brief Equal A function comparing two transforms (TransformType) for being identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the IndexToWorldTransform (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparison.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry::TransformType &leftHandSide,
                              const mitk::BaseGeometry::TransformType &rightHandSide,
                              ScalarType eps,
                              bool verbose);
 
   /**
   * @brief Equal A function comparing two bounding boxes (BoundingBoxType) for being identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the bounds (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparison.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry::BoundingBoxType &leftHandSide,
                              const mitk::BaseGeometry::BoundingBoxType &rightHandSide,
                              ScalarType eps,
                              bool verbose);
 
   /**
   * @brief A function checks if a test geometry is a sub geometry of
   * a given reference geometry.
   *
   * Sub geometry means that both geometries have the same voxel grid (same spacing, same axes,
   * origin is on voxel grid), but the bounding box of the checked geometry is contained or equal
   * to the bounding box of the reference geometry.\n
   * By this definition equal geometries are always sub geometries of each other.
   *
   * The function checks the spacing, origin, axis vectors, extents, the matrix of the
   * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag.
   *
   * The parameter eps is a tolerance value for all methods which are internally used for comparison.
   * @param testGeo Geometry that should be checked if it is a sub geometry of referenceGeo.
   * @param referenceGeo Geometry that should contain testedGeometry as sub geometry.
   * @param coordinateEps Tolerance for comparison of all spatial aspects (spacing, origin and grid alignment).
   * You can use mitk::eps in most cases.
   * @param directionEps Tolerance for comparison of all directional aspects (axis). You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparisons are true. False otherwise.
   */
   MITKCORE_EXPORT bool IsSubGeometry(const mitk::BaseGeometry& testGeo,
     const mitk::BaseGeometry& referenceGeo,
     ScalarType coordinateEps,
     ScalarType directionEps,
     bool verbose = false);
 
   /**
   * @brief A function checks if a test geometry is a sub geometry of
   * a given reference geometry.
   *
   * This is a overloaded version that uses a single tolerance for spatial and directional aspects. For more details,
   * see the other overloaded version.
   *
   * @param testGeo Geometry that should be checked if it is a sub geometry of referenceGeo.
   * @param referenceGeo Geometry that should contain testedGeometry as sub geometry.
   * @param eps Tolarence for comparison (both spatial and directional). You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False otherwise.
   */
   MITKCORE_EXPORT bool IsSubGeometry(const mitk::BaseGeometry& testGeo,
     const mitk::BaseGeometry& referenceGeo,
     ScalarType eps = mitk::eps,
     bool verbose = false);
 
 } // namespace mitk
 
 #endif
diff --git a/Modules/Core/include/mitkPlaneGeometry.h b/Modules/Core/include/mitkPlaneGeometry.h
index 0fec6a0794..81a0c3c0e1 100644
--- a/Modules/Core/include/mitkPlaneGeometry.h
+++ b/Modules/Core/include/mitkPlaneGeometry.h
@@ -1,608 +1,606 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #ifndef mitkPlaneGeometry_h
 #define mitkPlaneGeometry_h
 
 #include <MitkCoreExports.h>
 
 #include <mitkAnatomicalPlanes.h>
 #include <mitkBaseGeometry.h>
 #include <mitkRestorePlanePositionOperation.h>
 
 #include <vnl/vnl_cross.h>
 
 namespace mitk
 {
   template <class TCoordRep, unsigned int NPointDimension>
   class Line;
   typedef Line<ScalarType, 3> Line3D;
 
   /**
    * \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
    */
 
   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);
 
     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 AnatomicalPlane
     * (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,
                                          AnatomicalPlane planeorientation = AnatomicalPlane::Axial,
                                          ScalarType zPosition = 0,
                                          bool frontside = true,
                                          bool rotated = false,
                                          bool top = true);
 
     /**
     * \brief Initialize a plane with orientation \a AnatomicalPlane
     * (default: axial) with respect to \a BaseGeometry (default: identity).
     * Spacing also taken from \a BaseGeometry.
     *
     * \param geometry3D
     * \param top if \a true, create plane at top, otherwise at bottom
     * (for AnatomicalPlane Axial, for other plane locations respectively)
     * \param planeorientation
     * \param frontside
     * \param rotated
     */
     virtual void InitializeStandardPlane(const BaseGeometry *geometry3D,
                                          bool top,
                                          AnatomicalPlane planeorientation = AnatomicalPlane::Axial,
                                          bool frontside = true,
                                          bool rotated = false);
 
     /**
     * \brief Initialize a plane with orientation \a AnatomicalPlane
     * (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,
                                          AnatomicalPlane planeorientation = AnatomicalPlane::Axial,
                                          ScalarType zPosition = 0,
                                          bool frontside = true,
                                          bool rotated = false,
                                          bool top = true);
 
     /**
     * \brief Initialize plane with orientation \a AnatomicalPlane
     * (default: axial) given width, height and spacing.
     *
     */
     virtual void InitializeStandardPlane(ScalarType width,
                                          ScalarType height,
                                          const Vector3D &spacing,
                                          AnatomicalPlane planeorientation = AnatomicalPlane::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 Check if matrix is a rotation matrix:
     * - determinant is 1?
     * - R*R^T is ID?
     * Output warning otherwise.
     */
     static bool CheckRotationMatrix(AffineTransform3D *transform, double epsilon=1e-6);
 
     /**
     * \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 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 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 if 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
diff --git a/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp b/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
index 55ea6e4b07..6a3a8352e7 100644
--- a/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
+++ b/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
@@ -1,1101 +1,1185 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #include <iomanip>
 #include <sstream>
 #include <bitset>
 
 #include <vtkMatrix4x4.h>
 #include <vtkMatrixToLinearTransform.h>
 
 #include "mitkApplyTransformMatrixOperation.h"
 #include "mitkBaseGeometry.h"
 #include "mitkGeometryTransformHolder.h"
 #include "mitkInteractionConst.h"
 #include "mitkMatrixConvert.h"
 #include "mitkModifiedLock.h"
 #include "mitkPointOperation.h"
 #include "mitkRestorePlanePositionOperation.h"
 #include "mitkRotationOperation.h"
 #include "mitkScaleOperation.h"
 #include "mitkVector.h"
 #include "mitkMatrix.h"
 
 mitk::BaseGeometry::BaseGeometry()
   : Superclass(),
     mitk::OperationActor(),
     m_FrameOfReferenceID(0),
     m_IndexToWorldTransformLastModified(0),
     m_ImageGeometry(false),
     m_ModifiedLockFlag(false),
     m_ModifiedCalledFlag(false)
 {
   m_GeometryTransform = new GeometryTransformHolder();
   Initialize();
 }
 
 mitk::BaseGeometry::BaseGeometry(const BaseGeometry &other)
   : Superclass(),
     mitk::OperationActor(),
     m_FrameOfReferenceID(other.m_FrameOfReferenceID),
     m_IndexToWorldTransformLastModified(other.m_IndexToWorldTransformLastModified),
     m_ImageGeometry(other.m_ImageGeometry),
     m_ModifiedLockFlag(false),
     m_ModifiedCalledFlag(false)
 {
   m_GeometryTransform = new GeometryTransformHolder(*other.GetGeometryTransformHolder());
   other.InitializeGeometry(this);
 }
 
 mitk::BaseGeometry::~BaseGeometry()
 {
   delete m_GeometryTransform;
 }
 
 void mitk::BaseGeometry::SetVtkMatrixDeepCopy(vtkTransform *vtktransform)
 {
   m_GeometryTransform->SetVtkMatrixDeepCopy(vtktransform);
 }
 
 const mitk::Point3D mitk::BaseGeometry::GetOrigin() const
 {
   return m_GeometryTransform->GetOrigin();
 }
 
 void mitk::BaseGeometry::SetOrigin(const Point3D &origin)
 {
   mitk::ModifiedLock lock(this);
 
   if (origin != GetOrigin())
   {
     m_GeometryTransform->SetOrigin(origin);
     Modified();
   }
 }
 
 const mitk::Vector3D mitk::BaseGeometry::GetSpacing() const
 {
   return m_GeometryTransform->GetSpacing();
 }
 
 void mitk::BaseGeometry::Initialize()
 {
   float b[6] = {0, 1, 0, 1, 0, 1};
   SetFloatBounds(b);
 
   m_GeometryTransform->Initialize();
 
   m_FrameOfReferenceID = 0;
 
   m_ImageGeometry = false;
 }
 
 void mitk::BaseGeometry::SetFloatBounds(const float bounds[6])
 {
   mitk::BoundingBox::BoundsArrayType b;
   const float *input = bounds;
   int i = 0;
   for (mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6; ++i)
     *it++ = (mitk::ScalarType)*input++;
   SetBounds(b);
 }
 
 void mitk::BaseGeometry::SetFloatBounds(const double bounds[6])
 {
   mitk::BoundingBox::BoundsArrayType b;
   const double *input = bounds;
   int i = 0;
   for (mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); i < 6; ++i)
     *it++ = (mitk::ScalarType)*input++;
   SetBounds(b);
 }
 
 /** Initialize the geometry */
 void mitk::BaseGeometry::InitializeGeometry(BaseGeometry *newGeometry) const
 {
   newGeometry->SetBounds(m_BoundingBox->GetBounds());
 
   newGeometry->SetFrameOfReferenceID(GetFrameOfReferenceID());
 
   newGeometry->InitializeGeometryTransformHolder(this);
 
   newGeometry->m_ImageGeometry = m_ImageGeometry;
 }
 
 void mitk::BaseGeometry::InitializeGeometryTransformHolder(const BaseGeometry *otherGeometry)
 {
   this->m_GeometryTransform->Initialize(otherGeometry->GetGeometryTransformHolder());
 }
 
 /** Set the bounds */
 void mitk::BaseGeometry::SetBounds(const BoundsArrayType &bounds)
 {
   mitk::ModifiedLock lock(this);
 
   this->CheckBounds(bounds);
 
   m_BoundingBox = BoundingBoxType::New();
 
   BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New();
   BoundingBoxType::PointType p;
   BoundingBoxType::PointIdentifier pointid;
 
   for (pointid = 0; pointid < 2; ++pointid)
   {
     unsigned int i;
     for (i = 0; i < m_NDimensions; ++i)
     {
       p[i] = bounds[2 * i + pointid];
     }
     pointscontainer->InsertElement(pointid, p);
   }
 
   m_BoundingBox->SetPoints(pointscontainer);
   m_BoundingBox->ComputeBoundingBox();
   this->Modified();
 }
 
 void mitk::BaseGeometry::SetIndexToWorldTransform(mitk::AffineTransform3D *transform)
 {
   mitk::ModifiedLock lock(this);
 
   CheckIndexToWorldTransform(transform);
 
   m_GeometryTransform->SetIndexToWorldTransform(transform);
   Modified();
 }
 
 void mitk::BaseGeometry::SetIndexToWorldTransformWithoutChangingSpacing(mitk::AffineTransform3D *transform)
 {
   // security check
   mitk::Vector3D originalSpacing = this->GetSpacing();
 
   mitk::ModifiedLock lock(this);
 
   CheckIndexToWorldTransform(transform);
 
   m_GeometryTransform->SetIndexToWorldTransformWithoutChangingSpacing(transform);
   Modified();
 
   // Security check. Spacig must not have changed
   if (!mitk::Equal(originalSpacing, this->GetSpacing()))
   {
     MITK_WARN << "Spacing has changed in a method, where the spacing must not change.";
     assert(false);
   }
 }
 
 const mitk::BaseGeometry::BoundsArrayType mitk::BaseGeometry::GetBounds() const
 {
   assert(m_BoundingBox.IsNotNull());
   return m_BoundingBox->GetBounds();
 }
 
 bool mitk::BaseGeometry::IsValid() const
 {
   return true;
 }
 
 void mitk::BaseGeometry::SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing)
 {
   PreSetSpacing(aSpacing);
   _SetSpacing(aSpacing, enforceSetSpacing);
 }
 
 void mitk::BaseGeometry::_SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing)
 {
   m_GeometryTransform->SetSpacing(aSpacing, enforceSetSpacing);
 }
 
 mitk::Vector3D mitk::BaseGeometry::GetAxisVector(unsigned int direction) const
 {
   Vector3D frontToBack;
   frontToBack.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction).as_ref());
   frontToBack *= GetExtent(direction);
   return frontToBack;
 }
 
 mitk::ScalarType mitk::BaseGeometry::GetExtent(unsigned int direction) const
 {
   assert(m_BoundingBox.IsNotNull());
   if (direction >= m_NDimensions)
     mitkThrow() << "Direction is too big. This geometry is for 3D Data";
   BoundsArrayType bounds = m_BoundingBox->GetBounds();
   return bounds[direction * 2 + 1] - bounds[direction * 2];
 }
 
 bool mitk::BaseGeometry::Is2DConvertable()
 {
   bool isConvertableWithoutLoss = true;
   do
   {
     if (this->GetSpacing()[2] != 1)
     {
       isConvertableWithoutLoss = false;
       break;
     }
     if (this->GetOrigin()[2] != 0)
     {
       isConvertableWithoutLoss = false;
       break;
     }
     mitk::Vector3D col0, col1, col2;
     col0.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0).as_ref());
     col1.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1).as_ref());
     col2.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2).as_ref());
 
     if ((col0[2] != 0) || (col1[2] != 0) || (col2[0] != 0) || (col2[1] != 0) || (col2[2] != 1))
     {
       isConvertableWithoutLoss = false;
       break;
     }
   } while (false);
 
   return isConvertableWithoutLoss;
 }
 
 mitk::Point3D mitk::BaseGeometry::GetCenter() const
 {
   assert(m_BoundingBox.IsNotNull());
   Point3D c = m_BoundingBox->GetCenter();
   if (m_ImageGeometry)
   {
     // Get Center returns the middel of min and max pixel index. In corner based images, this is the right position.
     // In center based images (imageGeometry == true), the index needs to be shifted back.
     c[0] -= 0.5;
     c[1] -= 0.5;
     c[2] -= 0.5;
   }
   this->IndexToWorld(c, c);
   return c;
 }
 
 double mitk::BaseGeometry::GetDiagonalLength2() const
 {
   Vector3D diagonalvector = GetCornerPoint() - GetCornerPoint(false, false, false);
   return diagonalvector.GetSquaredNorm();
 }
 
 double mitk::BaseGeometry::GetDiagonalLength() const
 {
   return sqrt(GetDiagonalLength2());
 }
 
 mitk::Point3D mitk::BaseGeometry::GetCornerPoint(int id) const
 {
   assert(id >= 0);
   assert(this->IsBoundingBoxNull() == false);
 
   BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds();
 
   Point3D cornerpoint;
   switch (id)
   {
     case 0:
       FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[4]);
       break;
     case 1:
       FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[5]);
       break;
     case 2:
       FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[4]);
       break;
     case 3:
       FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[5]);
       break;
     case 4:
       FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[4]);
       break;
     case 5:
       FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[5]);
       break;
     case 6:
       FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[4]);
       break;
     case 7:
       FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[5]);
       break;
     default:
     {
       itkExceptionMacro(<< "A cube only has 8 corners. These are labeled 0-7.");
     }
   }
   if (m_ImageGeometry)
   {
     // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the
     // bounding box. The bounding box itself is no image, so it is corner-based
     FillVector3D(cornerpoint, cornerpoint[0] - 0.5, cornerpoint[1] - 0.5, cornerpoint[2] - 0.5);
   }
   return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint);
 }
 
 mitk::Point3D mitk::BaseGeometry::GetCornerPoint(bool xFront, bool yFront, bool zFront) const
 {
   assert(this->IsBoundingBoxNull() == false);
   BoundingBox::BoundsArrayType bounds = this->GetBoundingBox()->GetBounds();
 
   Point3D cornerpoint;
   cornerpoint[0] = (xFront ? bounds[0] : bounds[1]);
   cornerpoint[1] = (yFront ? bounds[2] : bounds[3]);
   cornerpoint[2] = (zFront ? bounds[4] : bounds[5]);
   if (m_ImageGeometry)
   {
     // Here i have to adjust the 0.5 offset manually, because the cornerpoint is the corner of the
     // bounding box. The bounding box itself is no image, so it is corner-based
     FillVector3D(cornerpoint, cornerpoint[0] - 0.5, cornerpoint[1] - 0.5, cornerpoint[2] - 0.5);
   }
 
   return this->GetIndexToWorldTransform()->TransformPoint(cornerpoint);
 }
 
 mitk::ScalarType mitk::BaseGeometry::GetExtentInMM(int direction) const
 {
   return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction).magnitude() *
          GetExtent(direction);
 }
 
 void mitk::BaseGeometry::SetExtentInMM(int direction, ScalarType extentInMM)
 {
   mitk::ModifiedLock lock(this);
 
   ScalarType len = GetExtentInMM(direction);
   if (fabs(len - extentInMM) >= mitk::eps)
   {
     AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix;
     vnlmatrix = m_GeometryTransform->GetVnlMatrix();
     if (len > extentInMM)
       vnlmatrix.set_column(direction, vnlmatrix.get_column(direction) / len * extentInMM);
     else
       vnlmatrix.set_column(direction, vnlmatrix.get_column(direction) * extentInMM / len);
     Matrix3D matrix;
     matrix = vnlmatrix;
     m_GeometryTransform->SetMatrix(matrix);
 
     Modified();
   }
 }
 
 bool mitk::BaseGeometry::IsInside(const mitk::Point3D &p) const
 {
   mitk::Point3D index;
   WorldToIndex(p, index);
   return IsIndexInside(index);
 }
 
 bool mitk::BaseGeometry::IsIndexInside(const mitk::Point3D &index) const
 {
   bool inside = false;
   // if it is an image geometry, we need to convert the index to discrete values
   // this is done by applying the rounding function also used in WorldToIndex (see line 323)
   if (m_ImageGeometry)
   {
     mitk::Point3D discretIndex;
     discretIndex[0] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[0]);
     discretIndex[1] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[1]);
     discretIndex[2] = itk::Math::RoundHalfIntegerUp<mitk::ScalarType>(index[2]);
 
     inside = this->GetBoundingBox()->IsInside(discretIndex);
     // we have to check if the index is at the upper border of each dimension,
     // because the boundingbox is not centerbased
     if (inside)
     {
       const BoundingBox::BoundsArrayType &bounds = this->GetBoundingBox()->GetBounds();
       if ((discretIndex[0] == bounds[1]) || (discretIndex[1] == bounds[3]) || (discretIndex[2] == bounds[5]))
         inside = false;
     }
   }
   else
     inside = this->GetBoundingBox()->IsInside(index);
 
   return inside;
 }
 
 void mitk::BaseGeometry::WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
 {
   mitk::Vector3D tempIn, tempOut;
   const TransformType::OffsetType &offset = this->GetIndexToWorldTransform()->GetOffset();
   tempIn = pt_mm.GetVectorFromOrigin() - offset;
 
   WorldToIndex(tempIn, tempOut);
 
   pt_units = Point3D(tempOut);
 }
 
 void mitk::BaseGeometry::WorldToIndex(const mitk::Vector3D &vec_mm, mitk::Vector3D &vec_units) const
 {
   // Get WorldToIndex transform
   if (m_IndexToWorldTransformLastModified != this->GetIndexToWorldTransform()->GetMTime())
   {
     if (!m_InvertedTransform)
     {
       m_InvertedTransform = TransformType::New();
     }
     if (!this->GetIndexToWorldTransform()->GetInverse(m_InvertedTransform.GetPointer()))
     {
       itkExceptionMacro("Internal ITK matrix inversion error, cannot proceed.");
     }
     m_IndexToWorldTransformLastModified = this->GetIndexToWorldTransform()->GetMTime();
   }
 
   // Check for valid matrix inversion
   const TransformType::MatrixType &inverse = m_InvertedTransform->GetMatrix();
   if (inverse.GetVnlMatrix().has_nans())
   {
     itkExceptionMacro("Internal ITK matrix inversion error, cannot proceed. Matrix was: "
                       << std::endl
                       << this->GetIndexToWorldTransform()->GetMatrix()
                       << "Suggested inverted matrix is:"
                       << std::endl
                       << inverse);
   }
 
   vec_units = inverse * vec_mm;
 }
 
 void mitk::BaseGeometry::WorldToIndex(const mitk::Point3D & /*atPt3d_mm*/,
                                       const mitk::Vector3D &vec_mm,
                                       mitk::Vector3D &vec_units) const
 {
   MITK_WARN << "Warning! Call of the deprecated function BaseGeometry::WorldToIndex(point, vec, vec). Use "
                "BaseGeometry::WorldToIndex(vec, vec) instead!";
   this->WorldToIndex(vec_mm, vec_units);
 }
 
 mitk::VnlVector mitk::BaseGeometry::GetOriginVnl() const
 {
   return GetOrigin().GetVnlVector();
 }
 
 vtkLinearTransform *mitk::BaseGeometry::GetVtkTransform() const
 {
   return m_GeometryTransform->GetVtkTransform();
 }
 
 void mitk::BaseGeometry::SetIdentity()
 {
   mitk::ModifiedLock lock(this);
 
   m_GeometryTransform->SetIdentity();
   Modified();
 }
 
 void mitk::BaseGeometry::Compose(const mitk::BaseGeometry::TransformType *other, bool pre)
 {
   mitk::ModifiedLock lock(this);
   m_GeometryTransform->Compose(other, pre);
   Modified();
 }
 
 void mitk::BaseGeometry::Compose(const vtkMatrix4x4 *vtkmatrix, bool pre)
 {
   mitk::BaseGeometry::TransformType::Pointer itkTransform = mitk::BaseGeometry::TransformType::New();
   TransferVtkMatrixToItkTransform(vtkmatrix, itkTransform.GetPointer());
   Compose(itkTransform, pre);
 }
 
 void mitk::BaseGeometry::Translate(const Vector3D &vector)
 {
   if ((vector[0] != 0) || (vector[1] != 0) || (vector[2] != 0))
   {
     this->SetOrigin(this->GetOrigin() + vector);
   }
 }
 
 void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D &pt_units, mitk::Point3D &pt_mm) const
 {
   pt_mm = this->GetIndexToWorldTransform()->TransformPoint(pt_units);
 }
 
 void mitk::BaseGeometry::IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const
 {
   vec_mm = this->GetIndexToWorldTransform()->TransformVector(vec_units);
 }
 
 void mitk::BaseGeometry::ExecuteOperation(Operation *operation)
 {
   mitk::ModifiedLock lock(this);
 
   vtkTransform *vtktransform = vtkTransform::New();
   vtktransform->SetMatrix(this->GetVtkMatrix());
   switch (operation->GetOperationType())
   {
     case OpNOTHING:
       break;
     case OpMOVE:
     {
       auto *pointOp = dynamic_cast<mitk::PointOperation *>(operation);
       if (pointOp == nullptr)
       {
         MITK_ERROR << "Point move operation is null!";
         return;
       }
       mitk::Point3D newPos = pointOp->GetPoint();
       ScalarType data[3];
       vtktransform->GetPosition(data);
       vtktransform->PostMultiply();
       vtktransform->Translate(newPos[0], newPos[1], newPos[2]);
       vtktransform->PreMultiply();
       break;
     }
     case OpSCALE:
     {
       auto *scaleOp = dynamic_cast<mitk::ScaleOperation *>(operation);
       if (scaleOp == nullptr)
       {
         MITK_ERROR << "Scale operation is null!";
         return;
       }
       mitk::Point3D newScale = scaleOp->GetScaleFactor();
       ScalarType scalefactor[3];
 
       scalefactor[0] = 1 + (newScale[0] / GetMatrixColumn(0).magnitude());
       scalefactor[1] = 1 + (newScale[1] / GetMatrixColumn(1).magnitude());
       scalefactor[2] = 1 + (newScale[2] / GetMatrixColumn(2).magnitude());
 
       mitk::Point3D anchor = scaleOp->GetScaleAnchorPoint();
 
       vtktransform->PostMultiply();
       vtktransform->Translate(-anchor[0], -anchor[1], -anchor[2]);
       vtktransform->Scale(scalefactor[0], scalefactor[1], scalefactor[2]);
       vtktransform->Translate(anchor[0], anchor[1], anchor[2]);
       break;
     }
     case OpROTATE:
     {
       auto *rotateOp = dynamic_cast<mitk::RotationOperation *>(operation);
       if (rotateOp == nullptr)
       {
         MITK_ERROR << "Rotation operation is null!";
         return;
       }
       Vector3D rotationVector = rotateOp->GetVectorOfRotation();
       Point3D center = rotateOp->GetCenterOfRotation();
       ScalarType angle = rotateOp->GetAngleOfRotation();
       vtktransform->PostMultiply();
       vtktransform->Translate(-center[0], -center[1], -center[2]);
       vtktransform->RotateWXYZ(angle, rotationVector[0], rotationVector[1], rotationVector[2]);
       vtktransform->Translate(center[0], center[1], center[2]);
       vtktransform->PreMultiply();
       break;
     }
     case OpRESTOREPLANEPOSITION:
     {
       // Copy necessary to avoid vtk warning
       vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
       TransferItkTransformToVtkMatrix(
         dynamic_cast<mitk::RestorePlanePositionOperation *>(operation)->GetTransform().GetPointer(), matrix);
       vtktransform->SetMatrix(matrix);
       matrix->Delete();
       break;
     }
     case OpAPPLYTRANSFORMMATRIX:
     {
       auto *applyMatrixOp = dynamic_cast<ApplyTransformMatrixOperation *>(operation);
       vtktransform->SetMatrix(applyMatrixOp->GetMatrix());
       break;
     }
     default:
       vtktransform->Delete();
       return;
   }
   this->SetVtkMatrixDeepCopy(vtktransform);
   Modified();
   vtktransform->Delete();
 }
 
 mitk::VnlVector mitk::BaseGeometry::GetMatrixColumn(unsigned int direction) const
 {
   return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(direction).as_ref();
 }
 
 mitk::BoundingBox::Pointer mitk::BaseGeometry::CalculateBoundingBoxRelativeToTransform(
   const mitk::AffineTransform3D *transform) const
 {
   mitk::BoundingBox::PointsContainer::Pointer pointscontainer = mitk::BoundingBox::PointsContainer::New();
 
   mitk::BoundingBox::PointIdentifier pointid = 0;
 
   unsigned char i;
   if (transform != nullptr)
   {
     mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New();
     transform->GetInverse(inverse);
     for (i = 0; i < 8; ++i)
       pointscontainer->InsertElement(pointid++, inverse->TransformPoint(GetCornerPoint(i)));
   }
   else
   {
     for (i = 0; i < 8; ++i)
       pointscontainer->InsertElement(pointid++, GetCornerPoint(i));
   }
 
   mitk::BoundingBox::Pointer result = mitk::BoundingBox::New();
   result->SetPoints(pointscontainer);
   result->ComputeBoundingBox();
 
   return result;
 }
 
 const std::string mitk::BaseGeometry::GetTransformAsString(TransformType *transformType)
 {
   std::ostringstream out;
 
   out << '[';
 
   for (int i = 0; i < 3; ++i)
   {
     out << '[';
     for (int j = 0; j < 3; ++j)
       out << transformType->GetMatrix().GetVnlMatrix().get(i, j) << ' ';
     out << ']';
   }
 
   out << "][";
 
   for (int i = 0; i < 3; ++i)
     out << transformType->GetOffset()[i] << ' ';
 
   out << "]\0";
 
   return out.str();
 }
 
 void mitk::BaseGeometry::SetIndexToWorldTransformByVtkMatrix(vtkMatrix4x4 *vtkmatrix)
 {
   m_GeometryTransform->SetIndexToWorldTransformByVtkMatrix(vtkmatrix);
 }
 
 void mitk::BaseGeometry::SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkMatrix4x4 *vtkmatrix)
 {
   m_GeometryTransform->SetIndexToWorldTransformByVtkMatrixWithoutChangingSpacing(vtkmatrix);
 }
 
 void mitk::BaseGeometry::IndexToWorld(const mitk::Point3D & /*atPt3d_units*/,
                                       const mitk::Vector3D &vec_units,
                                       mitk::Vector3D &vec_mm) const
 {
   MITK_WARN << "Warning! Call of the deprecated function BaseGeometry::IndexToWorld(point, vec, vec). Use "
                "BaseGeometry::IndexToWorld(vec, vec) instead!";
   // vec_mm = m_IndexToWorldTransform->TransformVector(vec_units);
   this->IndexToWorld(vec_units, vec_mm);
 }
 
 vtkMatrix4x4 *mitk::BaseGeometry::GetVtkMatrix()
 {
   return m_GeometryTransform->GetVtkMatrix();
 }
 
 bool mitk::BaseGeometry::IsBoundingBoxNull() const
 {
   return m_BoundingBox.IsNull();
 }
 
 bool mitk::BaseGeometry::IsIndexToWorldTransformNull() const
 {
   return m_GeometryTransform->IsIndexToWorldTransformNull();
 }
 
 void mitk::BaseGeometry::ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry)
 {
   // If Geometry is switched to ImageGeometry, you have to put an offset to the origin, because
   // imageGeometries origins are pixel-center-based
   // ... and remove the offset, if you switch an imageGeometry back to a normal geometry
   // For more information please see the Geometry documentation page
 
   if (m_ImageGeometry == isAnImageGeometry)
     return;
 
   const BoundingBox::BoundsArrayType &boundsarray = this->GetBoundingBox()->GetBounds();
 
   Point3D originIndex;
   FillVector3D(originIndex, boundsarray[0], boundsarray[2], boundsarray[4]);
 
   if (isAnImageGeometry == true)
     FillVector3D(originIndex, originIndex[0] + 0.5, originIndex[1] + 0.5, originIndex[2] + 0.5);
   else
     FillVector3D(originIndex, originIndex[0] - 0.5, originIndex[1] - 0.5, originIndex[2] - 0.5);
 
   Point3D originWorld;
 
   originWorld = GetIndexToWorldTransform()->TransformPoint(originIndex);
   // instead could as well call  IndexToWorld(originIndex,originWorld);
 
   SetOrigin(originWorld);
 
   this->SetImageGeometry(isAnImageGeometry);
 }
 
 void mitk::BaseGeometry::PrintSelf(std::ostream &os, itk::Indent indent) const
 {
   os << indent << " IndexToWorldTransform: ";
   if (this->IsIndexToWorldTransformNull())
     os << "nullptr" << std::endl;
   else
   {
     // from itk::MatrixOffsetTransformBase
     unsigned int i, j;
     os << std::endl;
     os << indent << "Matrix: " << std::endl;
     for (i = 0; i < 3; i++)
     {
       os << indent.GetNextIndent();
       for (j = 0; j < 3; j++)
       {
         os << this->GetIndexToWorldTransform()->GetMatrix()[i][j] << " ";
       }
       os << std::endl;
     }
 
     os << indent << "Offset: " << this->GetIndexToWorldTransform()->GetOffset() << std::endl;
     os << indent << "Center: " << this->GetIndexToWorldTransform()->GetCenter() << std::endl;
     os << indent << "Translation: " << this->GetIndexToWorldTransform()->GetTranslation() << std::endl;
 
     auto inverse = mitk::AffineTransform3D::New();
     if (this->GetIndexToWorldTransform()->GetInverse(inverse))
     {
       os << indent << "Inverse: " << std::endl;
       for (i = 0; i < 3; i++)
       {
         os << indent.GetNextIndent();
         for (j = 0; j < 3; j++)
         {
           os << inverse->GetMatrix()[i][j] << " ";
         }
         os << std::endl;
       }
     }
 
     // from itk::ScalableAffineTransform
     os << indent << "Scale : ";
     for (i = 0; i < 3; i++)
     {
       os << this->GetIndexToWorldTransform()->GetScale()[i] << " ";
     }
     os << std::endl;
   }
 
   os << indent << " BoundingBox: ";
   if (this->IsBoundingBoxNull())
     os << "nullptr" << std::endl;
   else
   {
     os << indent << "( ";
     for (unsigned int i = 0; i < 3; i++)
     {
       os << this->GetBoundingBox()->GetBounds()[2 * i] << "," << this->GetBoundingBox()->GetBounds()[2 * i + 1] << " ";
     }
     os << " )" << std::endl;
   }
 
   os << indent << " Origin: " << this->GetOrigin() << std::endl;
   os << indent << " ImageGeometry: " << this->GetImageGeometry() << std::endl;
   os << indent << " Spacing: " << this->GetSpacing() << std::endl;
 }
 
 void mitk::BaseGeometry::Modified() const
 {
   if (!m_ModifiedLockFlag)
     Superclass::Modified();
   else
     m_ModifiedCalledFlag = true;
 }
 
 mitk::AffineTransform3D *mitk::BaseGeometry::GetIndexToWorldTransform()
 {
   return m_GeometryTransform->GetIndexToWorldTransform();
 }
 
 const mitk::AffineTransform3D *mitk::BaseGeometry::GetIndexToWorldTransform() const
 {
   return m_GeometryTransform->GetIndexToWorldTransform();
 }
 
 const mitk::GeometryTransformHolder *mitk::BaseGeometry::GetGeometryTransformHolder() const
 {
   return m_GeometryTransform;
 }
 
+void mitk::BaseGeometry::MapAxesToOrientations(int axes[]) const
+{
+  auto affineTransform = this->GetIndexToWorldTransform();
+  auto matrix = affineTransform->GetMatrix();
+  matrix.GetVnlMatrix().normalize_columns();
+  auto inverseMatrix = matrix.GetInverse();
+
+  bool mapped[3] = {false, false, false};
+
+  // We need to allow an epsilon difference to ignore rounding.
+  const double eps = 0.0001;
+
+  for (int orientation = 0; orientation < 3; ++orientation)
+  {
+    auto absX = std::abs(inverseMatrix[0][orientation]);
+    auto absY = std::abs(inverseMatrix[1][orientation]);
+    auto absZ = std::abs(inverseMatrix[2][orientation]);
+
+    // First we check if there is a single maximum value. If there is, we found the axis
+    // that corresponds to the given orientation. If there is no single maximum value,
+    // we choose one from the the two or three axes that have the maximum value, but we
+    // need to make sure that we do not map the same axis to different orientations.
+    // Equal values are valid if the volume is rotated by exactly 45 degrees around one
+    // axis. If the volume is rotated by 45 degrees around two axes, you will get single
+    // maximum values at the same axes for two different orientations. In this case,
+    // the axis is mapped to one of the orientations, and for the other orientation we
+    // choose a different axis that has not been mapped yet, even if it is not a maximum.
+
+    if (absX > absY + eps)
+    {
+      if (absX > absZ + eps)
+      {
+        // x is the greatest
+        int axis = !mapped[0] ? 0 : !mapped[1] ? 1 : 2;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+      else
+      {
+        // z is the greatest OR x and z are equal and greater than y
+        int axis = !mapped[2] ? 2 : !mapped[0] ? 0 : 1;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+    }
+    else if (absY > absX + eps)
+    {
+      if (absY > absZ + eps)
+      {
+        // y is the greatest
+        int axis = !mapped[1] ? 1 : !mapped[2] ? 2 : 0;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+      else
+      {
+        // z is the greatest OR y and z are equal and greater than x
+        int axis = !mapped[2] ? 2 : !mapped[1] ? 1 : 0;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+    }
+    else
+    {
+      if (absZ > absX + eps)
+      {
+        // z is the greatest
+        int axis = !mapped[2] ? 2 : !mapped[0] ? 0 : 1;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+      else
+      {
+        // x and y are equal and greater than z OR x and y and z are equal
+        int axis = !mapped[0] ? 0 : !mapped[1] ? 1 : 2;
+        axes[orientation] = axis;
+        mapped[axis] = true;
+      }
+    }
+  }
+
+  assert(mapped[0] && mapped[1] && mapped[2]);
+}
+
 bool mitk::Equal(const mitk::BaseGeometry::BoundingBoxType &leftHandSide,
                  const mitk::BaseGeometry::BoundingBoxType &rightHandSide,
                  ScalarType eps,
                  bool verbose)
 {
   bool result = true;
 
   BaseGeometry::BoundsArrayType rightBounds = rightHandSide.GetBounds();
   BaseGeometry::BoundsArrayType leftBounds = leftHandSide.GetBounds();
   BaseGeometry::BoundsArrayType::Iterator itLeft = leftBounds.Begin();
   for (BaseGeometry::BoundsArrayType::Iterator itRight = rightBounds.Begin(); itRight != rightBounds.End(); ++itRight)
   {
     if ((!mitk::Equal(*itLeft, *itRight, eps)))
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D::BoundingBoxType )] bounds are not equal.";
         MITK_INFO << "rightHandSide is " << setprecision(12) << *itRight << " : leftHandSide is " << *itLeft
                   << " and tolerance is " << eps;
       }
       result = false;
     }
     itLeft++;
   }
   return result;
 }
 
 bool mitk::Equal(const mitk::BaseGeometry &leftHandSide,
                  const mitk::BaseGeometry &rightHandSide,
                  ScalarType coordinateEps,
                  ScalarType directionEps,
                  bool verbose)
 {
   bool result = true;
 
   // Compare spacings
   if (!mitk::Equal(leftHandSide.GetSpacing(), rightHandSide.GetSpacing(), coordinateEps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] Spacing differs.";
       MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetSpacing() << " : leftHandSide is "
                 << leftHandSide.GetSpacing() << " and tolerance is " << coordinateEps;
     }
     result = false;
   }
 
   // Compare Origins
   if (!mitk::Equal(leftHandSide.GetOrigin(), rightHandSide.GetOrigin(), coordinateEps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] Origin differs.";
       MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetOrigin() << " : leftHandSide is "
                 << leftHandSide.GetOrigin() << " and tolerance is " << coordinateEps;
     }
     result = false;
   }
 
   // Compare Axis and Extents
   for (unsigned int i = 0; i < 3; ++i)
   {
     if (!mitk::Equal(leftHandSide.GetAxisVector(i), rightHandSide.GetAxisVector(i), directionEps))
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] AxisVector #" << i << " differ";
         MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetAxisVector(i) << " : leftHandSide is "
                   << leftHandSide.GetAxisVector(i) << " and tolerance is " << directionEps;
       }
       result = false;
     }
 
     if (!mitk::Equal(leftHandSide.GetExtent(i), rightHandSide.GetExtent(i), coordinateEps))
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] Extent #" << i << " differ";
         MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetExtent(i) << " : leftHandSide is "
                   << leftHandSide.GetExtent(i) << " and tolerance is " << coordinateEps;
       }
       result = false;
     }
   }
 
   // Compare ImageGeometry Flag
   if (rightHandSide.GetImageGeometry() != leftHandSide.GetImageGeometry())
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] GetImageGeometry is different.";
       MITK_INFO << "rightHandSide is " << rightHandSide.GetImageGeometry() << " : leftHandSide is "
                 << leftHandSide.GetImageGeometry();
     }
     result = false;
   }
 
   // Compare FrameOfReference ID
   if (rightHandSide.GetFrameOfReferenceID() != leftHandSide.GetFrameOfReferenceID())
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] GetFrameOfReferenceID is different.";
       MITK_INFO << "rightHandSide is " << rightHandSide.GetFrameOfReferenceID() << " : leftHandSide is "
                 << leftHandSide.GetFrameOfReferenceID();
     }
     result = false;
   }
 
   // Compare BoundingBoxes
   if (!mitk::Equal(*leftHandSide.GetBoundingBox(), *rightHandSide.GetBoundingBox(), coordinateEps, verbose))
   {
     result = false;
   }
 
   // Compare IndexToWorldTransform Matrix
   if (!mitk::Equal(*leftHandSide.GetIndexToWorldTransform(), *rightHandSide.GetIndexToWorldTransform(), directionEps, verbose))
   {
     result = false;
   }
   return result;
 }
 
 bool mitk::Equal(const mitk::BaseGeometry& leftHandSide,
   const mitk::BaseGeometry& rightHandSide,
   ScalarType eps,
   bool verbose)
 {
   return Equal(leftHandSide, rightHandSide, eps, eps, verbose);
 }
 
 bool mitk::Equal(const mitk::BaseGeometry::TransformType &leftHandSide,
                  const mitk::BaseGeometry::TransformType &rightHandSide,
                  ScalarType eps,
                  bool verbose)
 {
   // Compare IndexToWorldTransform Matrix
   if (!mitk::MatrixEqualElementWise(leftHandSide.GetMatrix(), rightHandSide.GetMatrix(), eps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D::TransformType )] Index to World Transformation matrix differs.";
       MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetMatrix() << " : leftHandSide is "
                 << leftHandSide.GetMatrix() << " and tolerance is " << eps;
     }
     return false;
   }
   return true;
 }
 
 bool mitk::IsSubGeometry(const mitk::BaseGeometry& testGeo,
   const mitk::BaseGeometry& referenceGeo,
   ScalarType coordinateEps,
   ScalarType directionEps,
   bool verbose)
 {
   bool result = true;
 
   // Compare spacings (must be equal)
   const auto testedSpacing = testGeo.GetSpacing();
   if (!mitk::Equal(testedSpacing, referenceGeo.GetSpacing(), coordinateEps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] Spacing differs.";
       MITK_INFO << "testedGeometry is " << setprecision(12) << testedSpacing << " : referenceGeometry is "
         << referenceGeo.GetSpacing() << " and tolerance is " << coordinateEps;
     }
     result = false;
   }
 
   // Compare ImageGeometry Flag (must be equal)
   if (referenceGeo.GetImageGeometry() != testGeo.GetImageGeometry())
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] GetImageGeometry is different.";
       MITK_INFO << "referenceGeo is " << referenceGeo.GetImageGeometry() << " : testGeo is "
         << testGeo.GetImageGeometry();
     }
     result = false;
   }
 
   // Compare IndexToWorldTransform Matrix (must be equal -> same axis directions)
   if (!Equal(*(testGeo.GetIndexToWorldTransform()), *(referenceGeo.GetIndexToWorldTransform()), directionEps, verbose))
   {
     result = false;
   }
 
   //check if the geometry is within or equal to the bounds of reference geomentry.
   for (int i = 0; i<8; ++i)
   {
     auto testCorner = testGeo.GetCornerPoint(i);
     bool isInside = false;
     mitk::Point3D testCornerIndex;
     referenceGeo.WorldToIndex(testCorner, testCornerIndex);
 
     std::bitset<sizeof(int)> bs(i);
     //To regard the coordinateEps, we subtract or add it to the index elements
     //depending on whether it was constructed by a lower or an upper bound value
     //(see implementation of BaseGeometry::GetCorner()).
     if (bs.test(0))
     {
       testCornerIndex[2] -= coordinateEps;
     }
     else
     {
       testCornerIndex[2] += coordinateEps;
     }
 
     if (bs.test(1))
     {
       testCornerIndex[1] -= coordinateEps;
     }
     else
     {
       testCornerIndex[1] += coordinateEps;
     }
 
     if (bs.test(2))
     {
       testCornerIndex[0] -= coordinateEps;
     }
     else
     {
       testCornerIndex[0] += coordinateEps;
     }
 
     isInside = referenceGeo.IsIndexInside(testCornerIndex);
 
     if (!isInside)
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] corner point is not inside. ";
         MITK_INFO << "referenceGeo is " << setprecision(12) << referenceGeo << " : tested corner is "
           << testGeo.GetCornerPoint(i);
       }
       result = false;
     }
   }
 
   // check grid of test geometry is on the grid of the reference geometry. This is important as the
   // boundingbox is only checked for containing the tested geometry, but if a corner (one is enough
   // as we know that axis and spacing are equal, due to equal transfor (see above)) of the tested geometry
   // is on the grid it is really a sub geometry (as they have the same spacing and axis).
   auto cornerOffset = testGeo.GetCornerPoint(0) - referenceGeo.GetCornerPoint(0);
   mitk::Vector3D cornerIndexOffset;
   referenceGeo.WorldToIndex(cornerOffset, cornerIndexOffset);
   for (unsigned int i = 0; i < 3; ++i)
   {
     auto pixelCountContinous = cornerIndexOffset[i];
     auto pixelCount = std::round(pixelCountContinous);
     if (std::abs(pixelCount - pixelCountContinous) > coordinateEps)
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] Tested geometry is not on the grid of the reference geometry. ";
         MITK_INFO << "referenceGeo is " << setprecision(15) << referenceGeo << " : tested corner offset in pixels is "
           << pixelCountContinous << " for axis "<<i;
       }
       result = false;
     }
   }
 
   return result;
 }
 
 bool mitk::IsSubGeometry(const mitk::BaseGeometry& testGeo,
   const mitk::BaseGeometry& referenceGeo,
   ScalarType eps,
   bool verbose)
 {
   return IsSubGeometry(testGeo, referenceGeo, eps, eps, verbose);
 }
diff --git a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
index 30d2f6829b..585d7f9c1a 100644
--- a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
+++ b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
@@ -1,975 +1,929 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #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)
   {
   }
 
   bool PlaneGeometry::CheckRotationMatrix(mitk::AffineTransform3D *transform, double epsilon)
   {
     bool rotation = true;
 
     auto matrix = transform->GetMatrix().GetVnlMatrix();
     matrix.normalize_columns();
 
     auto det = vnl_determinant(matrix);
     if (fabs(det-1.0) > epsilon)
     {
       MITK_WARN << "Invalid rotation matrix! Determinant != 1 (" << det << ")";
       rotation = false;
     }
 
     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_WARN << "Invalid rotation matrix! R*R^T != ID. Max value: " << max << " (should be 0)";
       rotation = false;
     }
 
     return rotation;
   }
 
   void PlaneGeometry::CheckIndexToWorldTransform(mitk::AffineTransform3D *transform)
   {
     this->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,
                                               AnatomicalPlane 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 */,
                                               AnatomicalPlane 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), sagittal 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 AnatomicalPlane::Original:
       /** Orientation 'Original' shall be done like the axial plane,
        *  for whatever reasons. */
       case AnatomicalPlane::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;
 
       case AnatomicalPlane::Coronal: // Coronal=Frontal plane; cuts through patient's ear-ear-heel-heel:
         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 AnatomicalPlane::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 AnatomicalPlane");
     }
 
     VnlVector normal(3);
     FillVector3D(normal, 0, 0, 0);
     normal[normalDirection] = top ? 1 : -1;
 
     if ( transform != nullptr )
     {
       origin = transform->TransformPoint( origin );
       rightDV = transform->TransformVector( rightDV ).as_ref();
       bottomDV = transform->TransformVector( bottomDV ).as_ref();
       normal = transform->TransformVector( normal ).as_ref();
     }
 
     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.as_ref());
     matrix.GetVnlMatrix().set_column(1, bottomDV.as_ref());
     matrix.GetVnlMatrix().set_column(2, normal.as_ref());
     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,
                                               AnatomicalPlane 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();
-
+    auto matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
     matrix.GetVnlMatrix().normalize_columns();
-    mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetTranspose();
+    auto inverseMatrix = matrix.GetInverse();
+
+    // The index of the sagittal, coronal and axial axes in the reference geometry.
+    int axes[3];
+    geometry3D->MapAxesToOrientations(axes);
 
-    /// The index of the sagittal, coronal and axial axes in the reference geometry.
-    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.
+    // 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)
     {
-      int dominantAxis = axes.at(i);
-      directions[i] = itk::Function::Sign(inverseMatrix[dominantAxis][i]);
-      extents[i] = geometry3D->GetExtent(dominantAxis);
-      spacings[i] = geometry3D->GetSpacing()[dominantAxis];
+      int axis = axes[i];
+      directions[i] = itk::Function::Sign(inverseMatrix[axis][i]);
+      extents[i] = geometry3D->GetExtent(axis);
+      spacings[i] = geometry3D->GetSpacing()[axis];
     }
 
     // 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 AnatomicalPlane::Original:
     /** Orientation 'Original' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case AnatomicalPlane::Axial:
       width  = extents[0];
       height = extents[1];
       break;
     case AnatomicalPlane::Coronal:
       width  = extents[0];
       height = extents[2];
       break;
     case AnatomicalPlane::Sagittal:
       width  = extents[1];
       height = extents[2];
       break;
     default:
       itkExceptionMacro("unknown AnatomicalPlane");
     }
 
     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,
                                               AnatomicalPlane planeorientation,
                                               bool frontside,
                                               bool rotated)
   {
     /// The index of the sagittal, coronal and axial axes in world coordinate system.
     int worldAxis;
     switch(planeorientation)
     {
     case AnatomicalPlane::Original:
     /** Orientation 'Original' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case AnatomicalPlane::Axial:
       worldAxis = 2;
       break;
     case AnatomicalPlane::Coronal:
       worldAxis = 1;
       break;
     case AnatomicalPlane::Sagittal:
       worldAxis = 0;
       break;
     default:
       itkExceptionMacro("unknown AnatomicalPlane");
     }
 
-    // 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 = CalculateDominantAxes(inverseMatrix).at(worldAxis);
+    int axes[3];
+    geometry3D->MapAxesToOrientations(axes);
+    int axis = axes[worldAxis];
 
-    ScalarType zPosition = top ? 0.5 : geometry3D->GetExtent(dominantAxis) - 0.5;
+    ScalarType zPosition = top ? 0.5 : geometry3D->GetExtent(axis) - 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).as_ref());
 
     return frontToBack;
   }
 
   VnlVector PlaneGeometry::GetNormalVnl() const
   {
     return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2).as_ref();
   }
 
   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()[0] = p.GetVnlVector()[0];
     crossline.GetPoint()[1] = p.GetVnlVector()[1];
     crossline.GetPoint()[2] = p.GetVnlVector()[2];
 
     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 least 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 continuous 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 6e23ef32d9..d03779a6a1 100644
--- a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
+++ b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
@@ -1,962 +1,960 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #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 contained 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::AnatomicalPlane orientation,
                                               bool top,
                                               bool frontside,
                                               bool rotated)
 {
   m_ReferenceGeometry = geometry3D;
 
   PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New();
   planeGeometry->InitializeStandardPlane(geometry3D, top, orientation, frontside, rotated);
 
   int worldAxis = 
     orientation == AnatomicalPlane::Sagittal ? 0 :
     orientation == AnatomicalPlane::Coronal  ? 1 : 2;
 
-  // Inspired by:
-  // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
+  int axes[3];
+  geometry3D->MapAxesToOrientations(axes);
+  int axis = axes[worldAxis];
 
-  mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
-  matrix.GetVnlMatrix().normalize_columns();
-  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;
+  ScalarType viewSpacing = geometry3D->GetSpacing()[axis];
+  unsigned int slices = static_cast<unsigned int>(geometry3D->GetExtent(axis) + 0.5);
+
+  if (slices == 0 && geometry3D->GetExtent(axis) > 0)
+  {
+    // require at least one slice if there is _some_ extent
+    slices = 1;
   }
 
 #ifndef NDEBUG
-  int upDirection = itk::Function::Sign(inverseMatrix[dominantAxis][worldAxis]);
+  auto matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
+  matrix.GetVnlMatrix().normalize_columns();
+  auto inverseMatrix = matrix.GetInverse();
+
+  int upDirection = itk::Function::Sign(inverseMatrix[axis][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);
+  Vector3D worldPlaneNormal = inverseMatrix.get_row(axis) * (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->GetStepper()->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 available, 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;
           }
 
           // Perform 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->GetStepper()->SetSteps(m_Slices);
 
           this->Modified();
 
           // End Reinitialization
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->GetStepper()->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 e64ebc6e26..737222fc05 100644
--- a/Modules/Core/test/mitkPlaneGeometryTest.cpp
+++ b/Modules/Core/test/mitkPlaneGeometryTest.cpp
@@ -1,1093 +1,1084 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #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(TestCoronalInitialization);
   MITK_TEST(TestSagittalInitialization);
   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(), 1e-8);
     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: https://phabricator.mitk.org/T11477
      */
     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 occurred when the pointdistance of
    * the line was less than 1.
    * Test Behaviour:
    * ===============
    * 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 corresponding 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 describe 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] [column]
     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 TestSagittalInitialization()
   {
     mitk::Point3D cornerpoint0 = planegeometry->GetCornerPoint(0);
     mitk::PlaneGeometry::Pointer clonedplanegeometry = planegeometry->Clone();
 
     // Testing InitializeStandardPlane(clonedplanegeometry, anatomicalplane = Sagittal, zPosition = 0, frontside=true):
     planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::AnatomicalPlane::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 sagittally 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 sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(0), height, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(1), 1, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in units) of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(2), 1, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(0), heightInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing width, height and thickness (in mm) of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagittally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of sagittally 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, anatomicalplane = Axial, zPosition
     // = 0, frontside=false, rotated=true):
     planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::AnatomicalPlane::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 sagittally 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 TestCoronalInitialization()
   {
     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, anatomicalplane = Coronal, zPosition = 0, frontside=true)
     planegeometry->InitializeStandardPlane(clonedplanegeometry, mitk::AnatomicalPlane::Coronal);
     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 coronally 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 coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(0), width, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing height (in units) of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(1), 1, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing thickness (in units) of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtent(2), 1, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing width (in mm) of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing height (in mm) of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing thickness (in mm) of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of coronally 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, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetOrigin(), origin, testEps));
 
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(0), widthInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(2), 1, testEps));
 
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(2), newthicknessInMM, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally 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, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetOrigin(), origin, testEps));
 
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(0), widthInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in units) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtent(2), 1, testEps));
 
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(0), widthInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(1), thicknessInMM, testEps));
     CPPUNIT_ASSERT_MESSAGE(
       "Testing width, height and thickness (in mm) of unit spaced, coronally initialized version: ",
       mitk::Equal(planegeometry->GetExtentInMM(2), 1.0, testEps));
 
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(0), newright, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally initialized version: ",
                            mitk::Equal(planegeometry->GetAxisVector(1), newbottom, testEps));
     CPPUNIT_ASSERT_MESSAGE("Testing GetAxisVector() of unit spaced, coronally 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, anatomicalplane = 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 whether 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-parallel plane gets recognized as not parallel  [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-parallel plane gets recognized as not parallel  [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-parallel plane gets recognized as parallel  [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)