diff --git a/Modules/Core/include/mitkBaseGeometry.h b/Modules/Core/include/mitkBaseGeometry.h
index bc32fdf6e2..2f90ab6e04 100644
--- a/Modules/Core/include/mitkBaseGeometry.h
+++ b/Modules/Core/include/mitkBaseGeometry.h
@@ -1,764 +1,775 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef BaseGeometry_H_HEADER_INCLUDED
 #define BaseGeometry_H_HEADER_INCLUDED
 
 #include "mitkOperationActor.h"
 #include <MitkCoreExports.h>
 #include <mitkCommon.h>
 
 #include "itkScalableAffineTransform.h"
 #include "mitkNumericTypes.h"
 #include <itkAffineGeometryFrame.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;
 
   typedef itk::AffineGeometryFrame<ScalarType, 3> AffineGeometryFrame3D;
 
   //##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 descibes
   //## data objects consisting of slices, e.g., objects of type Image.
   //## Conversions between world coordinates (in mm) and unit coordinates
   //## (e.g., pixels in the case of an Image) can be performed.
   //##
   //## For more information on related classes, see \ref Geometry.
   //##
   //## BaseGeometry instances referring to an Image need a slightly
   //## different definition of corners, see SetImageGeometry. This
   //## is usualy 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 sublass 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 = 0);
 
     //##Documentation
     //## @brief Compose new IndexToWorldTransform with a given vtkMatrix4x4.
     //##
     //## Converts the vtkMatrix4x4 into a itk-transform and calls the previous method.
     void Compose(const vtkMatrix4x4 *vtkmatrix, bool pre = 0);
 
     //##Documentation
     //## @brief Translate the origin by a vector
     //##
     void Translate(const Vector3D &vector);
 
     //##Documentation
     //##@brief executes affine operations (translate, rotate, scale)
     virtual void ExecuteOperation(Operation *operation) 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 ((continous!)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 refering to an Image
       //##
       //## A geometry referring to an Image needs a slightly different
       //## definition of the position of the corners (see GetCornerPoint).
       //## The position of a voxel is defined by the position of its center.
       //## If we would use the origin (position of the (center of) the first
       //## voxel) as a corner and display this point, it would seem to be
       //## \em not at the corner but a bit within the image. Even worse for
       //## the opposite corner of the image: here the corner would appear
       //## outside the image (by half of the voxel diameter). Thus, we have
       //## to correct for this and to be able to do that, we need to know
       //## that the 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);
     virtual ~BaseGeometry();
 
     //##Documentation
     //## @brief clones the geometry
     //##
     //## Overwrite in all sub-classes.
     //## Normally looks like:
     //## \code
     //##  Self::Pointer newGeometry = new Self(*this);
     //##  newGeometry->UnRegister();
     //##  return newGeometry.GetPointer();
     //## \endcode
     virtual itk::LightObject::Pointer InternalClone() const override = 0;
 
     virtual 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 beeing identical.
   * @warning This method is deprecated and will not be available in the future. Use the \a bool mitk::Equal(const
   * mitk::mitk::BaseGeometry& g1, const mitk::BaseGeometry& g2) instead.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the spacing, origin, axisvectors, extents, the matrix of the
   * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag.
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * If you want to use different tolarance values for different parts of the geometry, feel free to use
   * the other comparison methods and write your own implementation of Equal.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   DEPRECATED(MITKCORE_EXPORT bool Equal(
     const mitk::BaseGeometry *leftHandSide, const mitk::BaseGeometry *rightHandSide, ScalarType eps, bool verbose));
 
   /**
   * @brief Equal A function comparing two geometries for beeing identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the spacing, origin, axisvectors, extents, the matrix of the
   * IndexToWorldTransform (elementwise), the bounding (elementwise) and the ImageGeometry flag.
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * If you want to use different tolarance values for different parts of the geometry, feel free to use
   * the other comparison methods and write your own implementation of Equal.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry &leftHandSide,
                              const mitk::BaseGeometry &rightHandSide,
                              ScalarType eps,
                              bool verbose);
 
   /**
   * @brief Equal A function comparing two transforms (TransformType) for beeing identical.
   * @warning This method is deprecated and will not be available in the future. Use the \a bool mitk::Equal(const
   * mitk::mitk::BaseGeometry::TransformType& t1, const mitk::BaseGeometry::TransformType& t2) instead.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the IndexToWorldTransform (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   DEPRECATED(MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry::TransformType *leftHandSide,
                                         const mitk::BaseGeometry::TransformType *rightHandSide,
                                         ScalarType eps,
                                         bool verbose));
 
   /**
   * @brief Equal A function comparing two transforms (TransformType) for beeing identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the IndexToWorldTransform (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   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 beeing identical.
   * @warning This method is deprecated and will not be available in the future. Use the \a bool mitk::Equal(const
   * mitk::mitk::BaseGeometry::BoundingBoxType& b1, const mitk::BaseGeometry::BoundingBoxType& b2) instead.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the bounds (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   DEPRECATED(MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry::BoundingBoxType *leftHandSide,
                                         const mitk::BaseGeometry::BoundingBoxType *rightHandSide,
                                         ScalarType eps,
                                         bool verbose));
 
   /**
   * @brief Equal A function comparing two bounding boxes (BoundingBoxType) for beeing identical.
   *
   * @ingroup MITKTestingAPI
   *
   * The function compares the bounds (elementwise).
   *
   * The parameter eps is a tolarence value for all methods which are internally used for comparion.
   * @param rightHandSide Compare this against leftHandSide.
   * @param leftHandSide Compare this against rightHandSide.
   * @param eps Tolarence for comparison. You can use mitk::eps in most cases.
   * @param verbose Flag indicating if the user wants detailed console output or not.
   * @return True, if all comparison are true. False in any other case.
   */
   MITKCORE_EXPORT bool Equal(const mitk::BaseGeometry::BoundingBoxType &leftHandSide,
                              const mitk::BaseGeometry::BoundingBoxType &rightHandSide,
                              ScalarType eps,
                              bool verbose);
 } // namespace mitk
 
 #endif /* BaseGeometry_H_HEADER_INCLUDED */
diff --git a/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp b/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
index 91aa7fbbf9..0e368b4cc5 100644
--- a/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
+++ b/Modules/Core/src/DataManagement/mitkBaseGeometry.cpp
@@ -1,1007 +1,1091 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <iomanip>
 #include <sstream>
 
 #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"
 
 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));
   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));
     col1.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1));
     col2.SetVnlVector(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2));
 
     if ((col0[2] != 0) || (col1[2] != 0) || (col2[0] != 0) || (col2[1] != 0) || (col2[2] != 1))
     {
       isConvertableWithoutLoss = false;
       break;
     }
   } while (0);
 
   return isConvertableWithoutLoss;
 }
 
 mitk::Point3D mitk::BaseGeometry::GetCenter() const
 {
   assert(m_BoundingBox.IsNotNull());
   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 = 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 const_cast<Self *>(this)->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:
     {
       mitk::PointOperation *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:
     {
       mitk::ScaleOperation *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:
     {
       mitk::RotationOperation *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:
     {
       ApplyTransformMatrixOperation *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);
 }
 
 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;
 
     os << indent << "Inverse: " << std::endl;
     for (i = 0; i < 3; i++)
     {
       os << indent.GetNextIndent();
       for (j = 0; j < 3; j++)
       {
         os << this->GetIndexToWorldTransform()->GetInverseMatrix()[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
+{
+  mitk::AffineTransform3D::ConstPointer affineTransform = this->GetIndexToWorldTransform();
+  mitk::AffineTransform3D::MatrixType matrix = affineTransform->GetMatrix();
+  matrix.GetVnlMatrix().normalize_columns();
+  mitk::AffineTransform3D::MatrixType::InternalMatrixType 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)
+  {
+    double absX = std::abs(inverseMatrix[0][orientation]);
+    double absY = std::abs(inverseMatrix[1][orientation]);
+    double 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)
 {
   if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
   {
     MITK_ERROR << "mitk::Equal( const mitk::Geometry3D::BoundingBoxType *leftHandSide, const "
                   "mitk::Geometry3D::BoundingBoxType *rightHandSide, ScalarType eps, bool verbose ) does not with nullptr "
                   "pointer input.";
     return false;
   }
   return Equal(*leftHandSide, *rightHandSide, eps, verbose);
 }
 
 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 eps,
                  bool verbose)
 {
   if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
   {
     MITK_ERROR << "mitk::Equal(const mitk::Geometry3D *leftHandSide, const mitk::Geometry3D *rightHandSide, ScalarType "
                   "eps, bool verbose) does not with nullptr pointer input.";
     return false;
   }
   return Equal(*leftHandSide, *rightHandSide, eps, verbose);
 }
 
 bool mitk::Equal(const mitk::BaseGeometry &leftHandSide,
                  const mitk::BaseGeometry &rightHandSide,
                  ScalarType eps,
                  bool verbose)
 {
   bool result = true;
 
   // Compare spacings
   if (!mitk::Equal(leftHandSide.GetSpacing(), rightHandSide.GetSpacing(), eps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] Spacing differs.";
       MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetSpacing() << " : leftHandSide is "
                 << leftHandSide.GetSpacing() << " and tolerance is " << eps;
     }
     result = false;
   }
 
   // Compare Origins
   if (!mitk::Equal(leftHandSide.GetOrigin(), rightHandSide.GetOrigin(), eps))
   {
     if (verbose)
     {
       MITK_INFO << "[( Geometry3D )] Origin differs.";
       MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetOrigin() << " : leftHandSide is "
                 << leftHandSide.GetOrigin() << " and tolerance is " << eps;
     }
     result = false;
   }
 
   // Compare Axis and Extents
   for (unsigned int i = 0; i < 3; ++i)
   {
     if (!mitk::Equal(leftHandSide.GetAxisVector(i), rightHandSide.GetAxisVector(i), eps))
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] AxisVector #" << i << " differ";
         MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetAxisVector(i) << " : leftHandSide is "
                   << leftHandSide.GetAxisVector(i) << " and tolerance is " << eps;
       }
       result = false;
     }
 
     if (!mitk::Equal(leftHandSide.GetExtent(i), rightHandSide.GetExtent(i), eps))
     {
       if (verbose)
       {
         MITK_INFO << "[( Geometry3D )] Extent #" << i << " differ";
         MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetExtent(i) << " : leftHandSide is "
                   << leftHandSide.GetExtent(i) << " and tolerance is " << eps;
       }
       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(), eps, verbose))
   {
     result = false;
   }
 
   // Compare IndexToWorldTransform Matrix
   if (!mitk::Equal(*leftHandSide.GetIndexToWorldTransform(), *rightHandSide.GetIndexToWorldTransform(), eps, verbose))
   {
     result = false;
   }
   return result;
 }
 
 bool mitk::Equal(const mitk::BaseGeometry::TransformType *leftHandSide,
                  const mitk::BaseGeometry::TransformType *rightHandSide,
                  ScalarType eps,
                  bool verbose)
 {
   if ((leftHandSide == nullptr) || (rightHandSide == nullptr))
   {
     MITK_ERROR << "mitk::Equal(const Geometry3D::TransformType *leftHandSide, const Geometry3D::TransformType "
                   "*rightHandSide, ScalarType eps, bool verbose ) does not with nullptr pointer input.";
     return false;
   }
   return Equal(*leftHandSide, *rightHandSide, 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;
 }
diff --git a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
index c9d64bffec..a27792ac80 100644
--- a/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
+++ b/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp
@@ -1,984 +1,968 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPlaneGeometry.h"
 #include "mitkInteractionConst.h"
 #include "mitkLine.h"
 #include "mitkPlaneOperation.h"
 
 #include <itkSpatialOrientationAdapter.h>
 
 #include <vtkTransform.h>
 
 #include <vnl/vnl_cross.h>
 
 namespace mitk
 {
   PlaneGeometry::PlaneGeometry() : Superclass(), m_ReferenceGeometry(nullptr) { Initialize(); }
   PlaneGeometry::~PlaneGeometry() {}
   PlaneGeometry::PlaneGeometry(const PlaneGeometry &other)
     : Superclass(other), m_ReferenceGeometry(other.m_ReferenceGeometry)
   {
   }
 
   void PlaneGeometry::EnsurePerpendicularNormal(mitk::AffineTransform3D *transform)
   {
     /** \brief ensure column(2) of indexToWorldTransform-matrix to be perpendicular to plane, keep length and
      * handedness.
      */
 
     VnlVector normal = vnl_cross_3d(transform->GetMatrix().GetVnlMatrix().get_column(0),
                                     transform->GetMatrix().GetVnlMatrix().get_column(1));
     normal.normalize(); // Now normal is a righthand normal unit vector, perpendicular to the plane.
 
     ScalarType len = transform->GetMatrix().GetVnlMatrix().get_column(2).two_norm();
     if (len == 0)
     {
       len = 1;
     }
     normal *= len;
 
     // Get the existing normal vector zed:
     vnl_vector_fixed<double, 3> zed = transform->GetMatrix().GetVnlMatrix().get_column(2);
 
     /** If det(matrix)<0, multiply normal vector by (-1) to keep geometry lefthanded. */
     if (vnl_determinant(transform->GetMatrix().GetVnlMatrix()) < 0)
     {
       MITK_DEBUG << "EnsurePerpendicularNormal(): Lefthanded geometry preserved, rh-normal: [ " << normal << " ],";
       normal *= (-1.0);
       MITK_DEBUG << "lh-normal: [ " << normal << " ], original vector zed is: [ " << zed << " ]";
     }
 
     // Now lets compare and only replace if necessary and only then warn the user:
 
     // float epsilon is precise enough here, but we need to respect numerical condition:
     // Higham, N., 2002, Accuracy and Stability of Numerical Algorithms,
     // SIAM, page 37, 2nd edition:
     double feps = std::numeric_limits<float>::epsilon();
     double zedsMagnitude = zed.two_norm();
     feps = feps * zedsMagnitude * 2;
 
     /** Check if normal (3. column) was perpendicular: If not, replace with calculated normal vector: */
     if (normal != zed)
     {
       vnl_vector_fixed<double, 3> parallel;
       for (unsigned int i = 0; i < 3; ++i)
       {
         parallel[i] = normal[i] / zed[i]; // Remember linear algebra: checking for parallelity.
       }
       // Checking if really not paralell i.e. non-colinear vectors by comparing these floating point numbers:
       if ((parallel[0] + feps < parallel[1] || feps + parallel[1] < parallel[0]) &&
           (parallel[0] + feps < parallel[2] || feps + parallel[2] < parallel[0]))
       {
         MITK_WARN
           << "EnsurePerpendicularNormal(): Plane geometry was _/askew/_, so here it gets rectified by substituting"
           << " the 3rd column of the indexToWorldMatrix with an appropriate normal vector: [ " << normal
           << " ], original vector zed was: [ " << zed << " ].";
 
         Matrix3D matrix = transform->GetMatrix();
         matrix.GetVnlMatrix().set_column(2, normal);
         transform->SetMatrix(matrix);
       }
     }
     else
     {
       // Nothing to do, 3rd column of indexToWorldTransformMatrix already was perfectly perpendicular.
     }
   }
 
   void PlaneGeometry::CheckIndexToWorldTransform(mitk::AffineTransform3D *transform)
   {
     EnsurePerpendicularNormal(transform);
   }
 
   void PlaneGeometry::CheckBounds(const BoundingBox::BoundsArrayType &bounds)
   {
     // error: unused parameter 'bounds'
     // this happens in release mode, where the assert macro is defined empty
     // hence we "use" the parameter:
     (void)bounds;
 
     // currently the unit rectangle must be starting at the origin [0,0]
     assert(bounds[0] == 0);
     assert(bounds[2] == 0);
     // the unit rectangle must be two-dimensional
     assert(bounds[1] > 0);
     assert(bounds[3] > 0);
   }
 
   void PlaneGeometry::IndexToWorld(const Point2D &pt_units, Point2D &pt_mm) const
   {
     pt_mm[0] = GetExtentInMM(0) / GetExtent(0) * pt_units[0];
     pt_mm[1] = GetExtentInMM(1) / GetExtent(1) * pt_units[1];
   }
 
   void PlaneGeometry::WorldToIndex(const Point2D &pt_mm, Point2D &pt_units) const
   {
     pt_units[0] = pt_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0)));
     pt_units[1] = pt_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(1)));
   }
 
   void PlaneGeometry::IndexToWorld(const Point2D & /*atPt2d_units*/, const Vector2D &vec_units, Vector2D &vec_mm) const
   {
     MITK_WARN << "Warning! Call of the deprecated function PlaneGeometry::IndexToWorld(point, vec, vec). Use "
                  "PlaneGeometry::IndexToWorld(vec, vec) instead!";
     this->IndexToWorld(vec_units, vec_mm);
   }
 
   void PlaneGeometry::IndexToWorld(const Vector2D &vec_units, Vector2D &vec_mm) const
   {
     vec_mm[0] = (GetExtentInMM(0) / GetExtent(0)) * vec_units[0];
     vec_mm[1] = (GetExtentInMM(1) / GetExtent(1)) * vec_units[1];
   }
 
   void PlaneGeometry::WorldToIndex(const Point2D & /*atPt2d_mm*/, const Vector2D &vec_mm, Vector2D &vec_units) const
   {
     MITK_WARN << "Warning! Call of the deprecated function PlaneGeometry::WorldToIndex(point, vec, vec). Use "
                  "PlaneGeometry::WorldToIndex(vec, vec) instead!";
     this->WorldToIndex(vec_mm, vec_units);
   }
 
   void PlaneGeometry::WorldToIndex(const Vector2D &vec_mm, Vector2D &vec_units) const
   {
     vec_units[0] = vec_mm[0] * (1.0 / (GetExtentInMM(0) / GetExtent(0)));
     vec_units[1] = vec_mm[1] * (1.0 / (GetExtentInMM(1) / GetExtent(1)));
   }
 
   void PlaneGeometry::InitializeStandardPlane(mitk::ScalarType width,
                                               ScalarType height,
                                               const Vector3D &spacing,
                                               PlaneGeometry::PlaneOrientation planeorientation,
                                               ScalarType zPosition,
                                               bool frontside,
                                               bool rotated,
                                               bool top)
   {
     AffineTransform3D::Pointer transform;
 
     transform = AffineTransform3D::New();
     AffineTransform3D::MatrixType matrix;
     AffineTransform3D::MatrixType::InternalMatrixType &vnlmatrix = matrix.GetVnlMatrix();
 
     vnlmatrix.set_identity();
     vnlmatrix(0, 0) = spacing[0];
     vnlmatrix(1, 1) = spacing[1];
     vnlmatrix(2, 2) = spacing[2];
     transform->SetIdentity();
     transform->SetMatrix(matrix);
 
     InitializeStandardPlane(width, height, transform.GetPointer(), planeorientation, zPosition, frontside, rotated, top);
   }
 
   void PlaneGeometry::InitializeStandardPlane(mitk::ScalarType width,
                                               mitk::ScalarType height,
                                               const AffineTransform3D *transform /* = nullptr */,
                                               PlaneGeometry::PlaneOrientation planeorientation /* = Axial */,
                                               mitk::ScalarType zPosition /* = 0 */,
                                               bool frontside /* = true */,
                                               bool rotated /* = false */,
                                               bool top /* = true */)
   {
     Superclass::Initialize();
 
     /// construct standard view.
 
     // We define at the moment "frontside" as: axial from above,
     // coronal from front (nose), saggital from right.
     // TODO: Double check with medicals doctors or radiologists [ ].
 
     // We define the orientation in patient's view, e.g. LAI is in a axial cut
     // (parallel to the triangle ear-ear-nose):
     // first axis: To the left ear of the patient
     // seecond axis: To the nose of the patient
     // third axis: To the legs of the patient.
 
     // Options are: L/R left/right; A/P anterior/posterior; I/S inferior/superior
     // (AKA caudal/cranial).
     // We note on all cases in the following switch block r.h. for right handed
     // or l.h. for left handed to describe the different cases.
     // However, which system is chosen is defined at the end of the switch block.
 
     // CAVE / be careful: the vectors right and bottom are relative to the plane
     // and do NOT describe e.g. the right side of the patient.
 
     Point3D origin;
     /** Bottom means downwards, DV means Direction Vector. Both relative to the image! */
     VnlVector rightDV(3), bottomDV(3);
     /** Origin of this plane is by default a zero vector and implicitly in the top-left corner: */
     origin.Fill(0);
     /** This is different to all definitions in MITK, except the QT mouse clicks.
     *   But it is like this here and we don't want to change a running system.
     *   Just be aware, that IN THIS FUNCTION we define the origin at the top left (e.g. your screen). */
 
     /** NormalDirection defines which axis (i.e. column index in the transform matrix)
     * is perpendicular to the plane: */
     int normalDirection;
 
     switch (planeorientation) // Switch through our limited choice of standard planes:
     {
       case None:
       /** Orientation 'None' shall be done like the axial plane orientation,
        *  for whatever reasons. */
       case Axial:
         if (frontside) // Radiologist's view from below. A cut along the triangle ear-ear-nose.
         {
           if (rotated == false)
           /** Origin in the top-left corner, x=[1; 0; 0], y=[0; 1; 0], z=[0; 0; 1],
           *   origin=[0,0,zpos]: LAI (r.h.)
           *
           *  0---rightDV---->                            |
           *  |                                           |
           *  |  Picture of a finite, rectangular plane   |
           *  |  ( insert LOLCAT-scan here ^_^ )          |
           *  |                                           |
           *  v  _________________________________________|
           *
           */
           {
             FillVector3D(origin, 0, 0, zPosition);
             FillVector3D(rightDV, 1, 0, 0);
             FillVector3D(bottomDV, 0, 1, 0);
           }
           else // Origin rotated to the bottom-right corner, x=[-1; 0; 0], y=[0; -1; 0], z=[0; 0; 1],
                // origin=[w,h,zpos]: RPI (r.h.)
           {    // Caveat emptor:  Still  using  top-left  as  origin  of  index  coordinate  system!
             FillVector3D(origin, width, height, zPosition);
             FillVector3D(rightDV, -1, 0, 0);
             FillVector3D(bottomDV, 0, -1, 0);
           }
         }
         else // 'Backside, not frontside.' Neuro-Surgeons's view from above patient.
         {
           if (rotated == false) // x=[-1; 0; 0], y=[0; 1; 0], z=[0; 0; 1], origin=[w,0,zpos]:  RAS (r.h.)
           {
             FillVector3D(origin, width, 0, zPosition);
             FillVector3D(rightDV, -1, 0, 0);
             FillVector3D(bottomDV, 0, 1, 0);
           }
           else // Origin in the bottom-left corner, x=[1; 0; 0], y=[0; -1; 0], z=[0; 0; 1],
                // origin=[0,h,zpos]:  LPS (r.h.)
           {
             FillVector3D(origin, 0, height, zPosition);
             FillVector3D(rightDV, 1, 0, 0);
             FillVector3D(bottomDV, 0, -1, 0);
           }
         }
         normalDirection = 2; // That is S=Superior=z=third_axis=middlefinger in righthanded LPS-system.
         break;
 
       // Frontal is known as Coronal in mitk. Plane cuts through patient's ear-ear-heel-heel:
       case Frontal:
         if (frontside)
         {
           if (rotated == false) // x=[1; 0; 0], y=[0; 0; 1], z=[0; 1; 0], origin=[0,zpos,0]: LAI (r.h.)
           {
             FillVector3D(origin, 0, zPosition, 0);
             FillVector3D(rightDV, 1, 0, 0);
             FillVector3D(bottomDV, 0, 0, 1);
           }
           else // x=[-1;0;0], y=[0;0;-1], z=[0;1;0], origin=[w,zpos,h]:  RAS  (r.h.)
           {
             FillVector3D(origin, width, zPosition, height);
             FillVector3D(rightDV, -1, 0, 0);
             FillVector3D(bottomDV, 0, 0, -1);
           }
         }
         else
         {
           if (rotated == false) //  x=[-1;0;0], y=[0;0;1], z=[0;1;0], origin=[w,zpos,0]: RPI (r.h.)
           {
             FillVector3D(origin, width, zPosition, 0);
             FillVector3D(rightDV, -1, 0, 0);
             FillVector3D(bottomDV, 0, 0, 1);
           }
           else //  x=[1;0;0], y=[0;1;0], z=[0;0;-1], origin=[0,zpos,h]: LPS (r.h.)
           {
             FillVector3D(origin, 0, zPosition, height);
             FillVector3D(rightDV, 1, 0, 0);
             FillVector3D(bottomDV, 0, 0, -1);
           }
         }
         normalDirection = 1; // Normal vector = posterior direction.
         break;
 
       case Sagittal: // Sagittal=Medial plane, the symmetry-plane mirroring your face.
         if (frontside)
         {
           if (rotated == false) //  x=[0;1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,0,0]:  LAI (r.h.)
           {
             FillVector3D(origin, zPosition, 0, 0);
             FillVector3D(rightDV, 0, 1, 0);
             FillVector3D(bottomDV, 0, 0, 1);
           }
           else //  x=[0;-1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,w,h]:  LPS (r.h.)
           {
             FillVector3D(origin, zPosition, width, height);
             FillVector3D(rightDV, 0, -1, 0);
             FillVector3D(bottomDV, 0, 0, -1);
           }
         }
         else
         {
           if (rotated == false) //  x=[0;-1;0], y=[0;0;1], z=[1;0;0], origin=[zpos,w,0]:  RPI (r.h.)
           {
             FillVector3D(origin, zPosition, width, 0);
             FillVector3D(rightDV, 0, -1, 0);
             FillVector3D(bottomDV, 0, 0, 1);
           }
           else //  x=[0;1;0], y=[0;0;-1], z=[1;0;0], origin=[zpos,0,h]:  RAS (r.h.)
           {
             FillVector3D(origin, zPosition, 0, height);
             FillVector3D(rightDV, 0, 1, 0);
             FillVector3D(bottomDV, 0, 0, -1);
           }
         }
         normalDirection = 0; // Normal vector = Lateral direction: Left in a LPS-system.
         break;
 
       default:
         itkExceptionMacro("unknown PlaneOrientation");
     }
 
     VnlVector normal(3);
     FillVector3D(normal, 0, 0, 0);
     normal[normalDirection] = top ? 1 : -1;
 
     if ( transform != nullptr )
     {
       origin = transform->TransformPoint( origin );
       rightDV = transform->TransformVector( rightDV );
       bottomDV = transform->TransformVector( bottomDV );
       normal = transform->TransformVector( normal );
     }
 
     ScalarType bounds[6] = {0, width, 0, height, 0, 1};
     this->SetBounds(bounds);
 
     AffineTransform3D::Pointer planeTransform = AffineTransform3D::New();
     Matrix3D matrix;
     matrix.GetVnlMatrix().set_column(0, rightDV);
     matrix.GetVnlMatrix().set_column(1, bottomDV);
     matrix.GetVnlMatrix().set_column(2, normal);
     planeTransform->SetMatrix(matrix);
     planeTransform->SetOffset(this->GetIndexToWorldTransform()->GetOffset());
     this->SetIndexToWorldTransform(planeTransform);
 
     this->SetOrigin(origin);
   }
 
   void PlaneGeometry::InitializeStandardPlane(const BaseGeometry *geometry3D,
                                               PlaneOrientation planeorientation,
                                               ScalarType zPosition,
                                               bool frontside,
                                               bool rotated,
                                               bool top)
   {
     this->SetReferenceGeometry(geometry3D);
 
     ScalarType width, height;
 
-    // Inspired by:
-    // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
-
     mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
 
     matrix.GetVnlMatrix().normalize_columns();
     mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
 
     /// The index of the sagittal, coronal and axial axes in the reference geometry.
     int axes[3];
+    geometry3D->MapAxesToOrientations(axes);
+
     /// 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 = itk::Function::Max3(
-          inverseMatrix[0][i],
-          inverseMatrix[1][i],
-          inverseMatrix[2][i]
-      );
-      axes[i] = dominantAxis;
-      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 None:
     /** Orientation 'None' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case Axial:
       width  = extents[0];
       height = extents[1];
       break;
     case Frontal:
       width  = extents[0];
       height = extents[2];
       break;
     case Sagittal:
       width  = extents[1];
       height = extents[2];
       break;
     default:
       itkExceptionMacro("unknown PlaneOrientation");
     }
 
     ScalarType bounds[6]= { 0, width, 0, height, 0, 1 };
     this->SetBounds( bounds );
 
     AffineTransform3D::Pointer transform = AffineTransform3D::New();
     transform->SetMatrix(matrix);
     transform->SetOffset(worldOrigin.GetVectorFromOrigin());
 
     InitializeStandardPlane(
       width, height, transform, planeorientation, zPosition, frontside, rotated, top);
   }
 
   void PlaneGeometry::InitializeStandardPlane(
     const BaseGeometry *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated)
   {
     /// The index of the sagittal, coronal and axial axes in world coordinate system.
     int worldAxis;
     switch(planeorientation)
     {
     case None:
     /** Orientation 'None' shall be done like the axial plane orientation,
      *  for whatever reasons. */
     case Axial:
       worldAxis = 2;
       break;
     case Frontal:
       worldAxis = 1;
       break;
     case Sagittal:
       worldAxis = 0;
       break;
     default:
       itkExceptionMacro("unknown PlaneOrientation");
     }
 
-    // Inspired by:
-    // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
-
-    mitk::AffineTransform3D::ConstPointer affineTransform = geometry3D->GetIndexToWorldTransform();
-    mitk::AffineTransform3D::MatrixType matrix = affineTransform->GetMatrix();
-    matrix.GetVnlMatrix().normalize_columns();
-    mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
-
-    /// The index of the sagittal, coronal and axial axes in the reference geometry.
-    int dominantAxis = itk::Function::Max3(
-        inverseMatrix[0][worldAxis],
-        inverseMatrix[1][worldAxis],
-        inverseMatrix[2][worldAxis]);
+    int 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));
 
     return frontToBack;
   }
 
   VnlVector PlaneGeometry::GetNormalVnl() const
   {
     return this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2);
   }
 
   ScalarType PlaneGeometry::DistanceFromPlane(const Point3D &pt3d_mm) const { return fabs(SignedDistance(pt3d_mm)); }
   ScalarType PlaneGeometry::SignedDistance(const Point3D &pt3d_mm) const { return SignedDistanceFromPlane(pt3d_mm); }
   bool PlaneGeometry::IsAbove(const Point3D &pt3d_mm, bool considerBoundingBox) const
   {
     if (considerBoundingBox)
     {
       Point3D pt3d_units;
       BaseGeometry::WorldToIndex(pt3d_mm, pt3d_units);
       return (pt3d_units[2] > this->GetBoundingBox()->GetBounds()[4]);
     }
     else
       return SignedDistanceFromPlane(pt3d_mm) > 0;
   }
 
   bool PlaneGeometry::IntersectionLine(const PlaneGeometry *plane, Line3D &crossline) const
   {
     Vector3D normal = this->GetNormal();
     normal.Normalize();
 
     Vector3D planeNormal = plane->GetNormal();
     planeNormal.Normalize();
 
     Vector3D direction = itk::CrossProduct(normal, planeNormal);
 
     if (direction.GetSquaredNorm() < eps)
       return false;
 
     crossline.SetDirection(direction);
 
     double N1dN2 = normal * planeNormal;
     double determinant = 1.0 - N1dN2 * N1dN2;
 
     Vector3D origin = this->GetOrigin().GetVectorFromOrigin();
     Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin();
 
     double d1 = normal * origin;
     double d2 = planeNormal * planeOrigin;
 
     double c1 = (d1 - d2 * N1dN2) / determinant;
     double c2 = (d2 - d1 * N1dN2) / determinant;
 
     Vector3D p = normal * c1 + planeNormal * c2;
     crossline.GetPoint().GetVnlVector() = p.GetVnlVector();
 
     return true;
   }
 
   unsigned int PlaneGeometry::IntersectWithPlane2D(const PlaneGeometry *plane, Point2D &lineFrom, Point2D &lineTo) const
   {
     Line3D crossline;
     if (this->IntersectionLine(plane, crossline) == false)
       return 0;
 
     Point2D point2;
     Vector2D direction2;
 
     this->Map(crossline.GetPoint(), point2);
     this->Map(crossline.GetPoint(), crossline.GetDirection(), direction2);
 
     return Line3D::RectangleLineIntersection(
       0, 0, GetExtentInMM(0), GetExtentInMM(1), point2, direction2, lineFrom, lineTo);
   }
 
   double PlaneGeometry::Angle(const PlaneGeometry *plane) const
   {
     return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2));
   }
 
   double PlaneGeometry::Angle(const Line3D &line) const
   {
     return vnl_math::pi_over_2 - angle(line.GetDirection().GetVnlVector(), GetMatrixColumn(2));
   }
 
   bool PlaneGeometry::IntersectionPoint(const Line3D &line, Point3D &intersectionPoint) const
   {
     Vector3D planeNormal = this->GetNormal();
     planeNormal.Normalize();
 
     Vector3D lineDirection = line.GetDirection();
     lineDirection.Normalize();
 
     double t = planeNormal * lineDirection;
     if (fabs(t) < eps)
     {
       return false;
     }
 
     Vector3D diff;
     diff = this->GetOrigin() - line.GetPoint();
     t = (planeNormal * diff) / t;
 
     intersectionPoint = line.GetPoint() + lineDirection * t;
     return true;
   }
 
   bool PlaneGeometry::IntersectionPointParam(const Line3D &line, double &t) const
   {
     Vector3D planeNormal = this->GetNormal();
 
     Vector3D lineDirection = line.GetDirection();
 
     t = planeNormal * lineDirection;
 
     if (fabs(t) < eps)
     {
       return false;
     }
 
     Vector3D diff;
     diff = this->GetOrigin() - line.GetPoint();
     t = (planeNormal * diff) / t;
     return true;
   }
 
   bool PlaneGeometry::IsParallel(const PlaneGeometry *plane) const
   {
     return ((Angle(plane) < 10.0 * mitk::sqrteps) || (Angle(plane) > (vnl_math::pi - 10.0 * sqrteps)));
   }
 
   bool PlaneGeometry::IsOnPlane(const Point3D &point) const { return Distance(point) < eps; }
   bool PlaneGeometry::IsOnPlane(const Line3D &line) const
   {
     return ((Distance(line.GetPoint()) < eps) && (Distance(line.GetPoint2()) < eps));
   }
 
   bool PlaneGeometry::IsOnPlane(const PlaneGeometry *plane) const
   {
     return (IsParallel(plane) && (Distance(plane->GetOrigin()) < eps));
   }
 
   Point3D PlaneGeometry::ProjectPointOntoPlane(const Point3D &pt) const
   {
     ScalarType len = this->GetNormalVnl().two_norm();
     return pt - this->GetNormal() * this->SignedDistanceFromPlane(pt) / len;
   }
 
   itk::LightObject::Pointer PlaneGeometry::InternalClone() const
   {
     Self::Pointer newGeometry = new PlaneGeometry(*this);
     newGeometry->UnRegister();
     return newGeometry.GetPointer();
   }
 
   void PlaneGeometry::ExecuteOperation(Operation *operation)
   {
     vtkTransform *transform = vtkTransform::New();
     transform->SetMatrix(this->GetVtkMatrix());
 
     switch (operation->GetOperationType())
     {
       case OpORIENT:
       {
         mitk::PlaneOperation *planeOp = dynamic_cast<mitk::PlaneOperation *>(operation);
         if (planeOp == nullptr)
         {
           return;
         }
 
         Point3D center = planeOp->GetPoint();
 
         Vector3D orientationVector = planeOp->GetNormal();
         Vector3D defaultVector;
         FillVector3D(defaultVector, 0.0, 0.0, 1.0);
 
         Vector3D rotationAxis = itk::CrossProduct(orientationVector, defaultVector);
         // double rotationAngle = acos( orientationVector[2] / orientationVector.GetNorm() );
 
         double rotationAngle = atan2((double)rotationAxis.GetNorm(), (double)(orientationVector * defaultVector));
         rotationAngle *= 180.0 / vnl_math::pi;
 
         transform->PostMultiply();
         transform->Identity();
         transform->Translate(center[0], center[1], center[2]);
         transform->RotateWXYZ(rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2]);
         transform->Translate(-center[0], -center[1], -center[2]);
         break;
       }
       case OpRESTOREPLANEPOSITION:
       {
         RestorePlanePositionOperation *op = dynamic_cast<mitk::RestorePlanePositionOperation *>(operation);
         if (op == nullptr)
         {
           return;
         }
 
         AffineTransform3D::Pointer transform2 = AffineTransform3D::New();
         Matrix3D matrix;
         matrix.GetVnlMatrix().set_column(0, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(0));
         matrix.GetVnlMatrix().set_column(1, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(1));
         matrix.GetVnlMatrix().set_column(2, op->GetTransform()->GetMatrix().GetVnlMatrix().get_column(2));
         transform2->SetMatrix(matrix);
         Vector3D offset = op->GetTransform()->GetOffset();
         transform2->SetOffset(offset);
 
         this->SetIndexToWorldTransform(transform2);
         ScalarType bounds[6] = {0, op->GetWidth(), 0, op->GetHeight(), 0, 1};
         this->SetBounds(bounds);
         this->Modified();
         transform->Delete();
         return;
       }
       default:
         Superclass::ExecuteOperation(operation);
         transform->Delete();
         return;
     }
 
     this->SetVtkMatrixDeepCopy(transform);
     this->Modified();
     transform->Delete();
   }
 
   void PlaneGeometry::PrintSelf(std::ostream &os, itk::Indent indent) const
   {
     Superclass::PrintSelf(os, indent);
     os << indent << " ScaleFactorMMPerUnitX: " << GetExtentInMM(0) / GetExtent(0) << std::endl;
     os << indent << " ScaleFactorMMPerUnitY: " << GetExtentInMM(1) / GetExtent(1) << std::endl;
     os << indent << " Normal: " << GetNormal() << std::endl;
   }
 
   bool PlaneGeometry::Map(const mitk::Point3D &pt3d_mm, mitk::Point2D &pt2d_mm) const
   {
     assert(this->IsBoundingBoxNull() == false);
 
     Point3D pt3d_units;
     Superclass::WorldToIndex(pt3d_mm, pt3d_units);
     pt2d_mm[0] = pt3d_units[0] * GetExtentInMM(0) / GetExtent(0);
     pt2d_mm[1] = pt3d_units[1] * GetExtentInMM(1) / GetExtent(1);
     pt3d_units[2] = 0;
     return const_cast<BoundingBox *>(this->GetBoundingBox())->IsInside(pt3d_units);
   }
 
   void PlaneGeometry::Map(const mitk::Point2D &pt2d_mm, mitk::Point3D &pt3d_mm) const
   {
     // pt2d_mm is measured from the origin of the world geometry (at leats it called form BaseRendere::Mouse...Event)
     Point3D pt3d_units;
     pt3d_units[0] = pt2d_mm[0] / (GetExtentInMM(0) / GetExtent(0));
     pt3d_units[1] = pt2d_mm[1] / (GetExtentInMM(1) / GetExtent(1));
     pt3d_units[2] = 0;
     // pt3d_units is a continuos index. We divided it with the Scale Factor (= spacing in x and y) to convert it from mm
     // to index units.
     //
     pt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units);
     // now we convert the 3d index to a 3D world point in mm. We could have used IndexToWorld as well as
     // GetITW->Transform...
   }
 
   void PlaneGeometry::SetSizeInUnits(mitk::ScalarType width, mitk::ScalarType height)
   {
     ScalarType bounds[6] = {0, width, 0, height, 0, 1};
     ScalarType extent, newextentInMM;
     if (GetExtent(0) > 0)
     {
       extent = GetExtent(0);
       if (width > extent)
         newextentInMM = GetExtentInMM(0) / width * extent;
       else
         newextentInMM = GetExtentInMM(0) * extent / width;
       SetExtentInMM(0, newextentInMM);
     }
     if (GetExtent(1) > 0)
     {
       extent = GetExtent(1);
       if (width > extent)
         newextentInMM = GetExtentInMM(1) / height * extent;
       else
         newextentInMM = GetExtentInMM(1) * extent / height;
       SetExtentInMM(1, newextentInMM);
     }
     SetBounds(bounds);
   }
 
   bool PlaneGeometry::Project(const mitk::Point3D &pt3d_mm, mitk::Point3D &projectedPt3d_mm) const
   {
     assert(this->IsBoundingBoxNull() == false);
 
     Point3D pt3d_units;
     Superclass::WorldToIndex(pt3d_mm, pt3d_units);
     pt3d_units[2] = 0;
     projectedPt3d_mm = GetIndexToWorldTransform()->TransformPoint(pt3d_units);
     return const_cast<BoundingBox *>(this->GetBoundingBox())->IsInside(pt3d_units);
   }
 
   bool PlaneGeometry::Project(const mitk::Vector3D &vec3d_mm, mitk::Vector3D &projectedVec3d_mm) const
   {
     assert(this->IsBoundingBoxNull() == false);
 
     Vector3D vec3d_units;
     Superclass::WorldToIndex(vec3d_mm, vec3d_units);
     vec3d_units[2] = 0;
     projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units);
     return true;
   }
 
   bool PlaneGeometry::Project(const mitk::Point3D &atPt3d_mm,
                               const mitk::Vector3D &vec3d_mm,
                               mitk::Vector3D &projectedVec3d_mm) const
   {
     MITK_WARN << "Deprecated function! Call Project(vec3D,vec3D) instead.";
     assert(this->IsBoundingBoxNull() == false);
 
     Vector3D vec3d_units;
     Superclass::WorldToIndex(atPt3d_mm, vec3d_mm, vec3d_units);
     vec3d_units[2] = 0;
     projectedVec3d_mm = GetIndexToWorldTransform()->TransformVector(vec3d_units);
 
     Point3D pt3d_units;
     Superclass::WorldToIndex(atPt3d_mm, pt3d_units);
     return const_cast<BoundingBox *>(this->GetBoundingBox())->IsInside(pt3d_units);
   }
 
   bool PlaneGeometry::Map(const mitk::Point3D &atPt3d_mm,
                           const mitk::Vector3D &vec3d_mm,
                           mitk::Vector2D &vec2d_mm) const
   {
     Point2D pt2d_mm_start, pt2d_mm_end;
     Point3D pt3d_mm_end;
     bool inside = Map(atPt3d_mm, pt2d_mm_start);
     pt3d_mm_end = atPt3d_mm + vec3d_mm;
     inside &= Map(pt3d_mm_end, pt2d_mm_end);
     vec2d_mm = pt2d_mm_end - pt2d_mm_start;
     return inside;
   }
 
   void PlaneGeometry::Map(const mitk::Point2D & /*atPt2d_mm*/,
                           const mitk::Vector2D & /*vec2d_mm*/,
                           mitk::Vector3D & /*vec3d_mm*/) const
   {
     //@todo implement parallel to the other Map method!
     assert(false);
   }
 
   void PlaneGeometry::SetReferenceGeometry(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 1bca8b48ba..50af2033d3 100644
--- a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
+++ b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp
@@ -1,966 +1,958 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <itkSpatialOrientationAdapter.h>
 
 #include "mitkSlicedGeometry3D.h"
 #include "mitkAbstractTransformGeometry.h"
 #include "mitkApplyTransformMatrixOperation.h"
 #include "mitkInteractionConst.h"
 #include "mitkPlaneGeometry.h"
 #include "mitkPlaneOperation.h"
 #include "mitkRestorePlanePositionOperation.h"
 #include "mitkRotationOperation.h"
 #include "mitkSliceNavigationController.h"
 
 const mitk::ScalarType PI = 3.14159265359;
 
 mitk::SlicedGeometry3D::SlicedGeometry3D()
   : m_EvenlySpaced(true), m_Slices(0), m_ReferenceGeometry(nullptr), m_SliceNavigationController(nullptr)
 {
   m_DirectionVector.Fill(0);
   this->InitializeSlicedGeometry(m_Slices);
 }
 
 mitk::SlicedGeometry3D::SlicedGeometry3D(const SlicedGeometry3D &other)
   : Superclass(other),
     m_EvenlySpaced(other.m_EvenlySpaced),
     m_Slices(other.m_Slices),
     m_ReferenceGeometry(other.m_ReferenceGeometry),
     m_SliceNavigationController(other.m_SliceNavigationController)
 {
   m_DirectionVector.Fill(0);
   SetSpacing(other.GetSpacing());
   SetDirectionVector(other.GetDirectionVector());
 
   if (m_EvenlySpaced)
   {
     assert(!other.m_PlaneGeometries.empty() && "This may happen when you use one of the old Initialize methods, which had a bool parameter that is implicitly casted to the number of slices now.");
     PlaneGeometry::Pointer geometry = other.m_PlaneGeometries[0]->Clone();
     assert(geometry.IsNotNull());
     SetPlaneGeometry(geometry, 0);
   }
   else
   {
     unsigned int s;
     for (s = 0; s < other.m_Slices; ++s)
     {
       if (other.m_PlaneGeometries[s].IsNull())
       {
         assert(other.m_EvenlySpaced);
         m_PlaneGeometries[s] = nullptr;
       }
       else
       {
         PlaneGeometry *geometry2D = other.m_PlaneGeometries[s]->Clone();
         assert(geometry2D != nullptr);
         SetPlaneGeometry(geometry2D, s);
       }
     }
   }
 }
 
 mitk::SlicedGeometry3D::~SlicedGeometry3D()
 {
 }
 
 mitk::PlaneGeometry *mitk::SlicedGeometry3D::GetPlaneGeometry(int s) const
 {
   mitk::PlaneGeometry::Pointer geometry2D = nullptr;
 
   if (this->IsValidSlice(s))
   {
     geometry2D = m_PlaneGeometries[s];
 
     // If (a) m_EvenlySpaced==true, (b) we don't have a PlaneGeometry stored
     // for the requested slice, and (c) the first slice (s=0)
     // is a PlaneGeometry instance, then we calculate the geometry of the
     // requested as the plane of the first slice shifted by m_Spacing[2]*s
     // in the direction of m_DirectionVector.
     if ((m_EvenlySpaced) && (geometry2D.IsNull()))
     {
       PlaneGeometry *firstSlice = m_PlaneGeometries[0];
 
       if (firstSlice != nullptr &&
           dynamic_cast<AbstractTransformGeometry *>(m_PlaneGeometries[0].GetPointer()) == nullptr)
       {
         if ((m_DirectionVector[0] == 0.0) && (m_DirectionVector[1] == 0.0) && (m_DirectionVector[2] == 0.0))
         {
           m_DirectionVector = firstSlice->GetNormal();
           m_DirectionVector.Normalize();
         }
 
         Vector3D direction;
         direction = m_DirectionVector * this->GetSpacing()[2];
 
         mitk::PlaneGeometry::Pointer requestedslice;
         requestedslice = static_cast<mitk::PlaneGeometry *>(firstSlice->Clone().GetPointer());
 
         requestedslice->SetOrigin(requestedslice->GetOrigin() + direction * s);
 
         geometry2D = requestedslice;
         m_PlaneGeometries[s] = geometry2D;
       }
     }
     return geometry2D;
   }
   else
   {
     return nullptr;
   }
 }
 
 const mitk::BoundingBox *mitk::SlicedGeometry3D::GetBoundingBox() const
 {
   assert(this->IsBoundingBoxNull() == false);
   return Superclass::GetBoundingBox();
 }
 
 bool mitk::SlicedGeometry3D::SetPlaneGeometry(mitk::PlaneGeometry *geometry2D, int s)
 {
   if (this->IsValidSlice(s))
   {
     m_PlaneGeometries[s] = geometry2D;
     m_PlaneGeometries[s]->SetReferenceGeometry(m_ReferenceGeometry);
     return true;
   }
   return false;
 }
 
 void mitk::SlicedGeometry3D::InitializeSlicedGeometry(unsigned int slices)
 {
   Superclass::Initialize();
   m_Slices = slices;
 
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   Vector3D spacing;
   spacing.Fill(1.0);
   this->SetSpacing(spacing);
 
   m_DirectionVector.Fill(0);
 }
 
 void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, unsigned int slices)
 {
   assert(geometry2D != nullptr);
   this->InitializeEvenlySpaced(geometry2D, geometry2D->GetExtentInMM(2) / geometry2D->GetExtent(2), slices);
 }
 
 void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D,
                                                     mitk::ScalarType zSpacing,
                                                     unsigned int slices)
 {
   assert(geometry2D != nullptr);
   assert(geometry2D->GetExtent(0) > 0);
   assert(geometry2D->GetExtent(1) > 0);
 
   geometry2D->Register();
 
   Superclass::Initialize();
   m_Slices = slices;
 
   BoundingBox::BoundsArrayType bounds = geometry2D->GetBounds();
   bounds[4] = 0;
   bounds[5] = slices;
 
   // clear and reserve
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   Vector3D directionVector = geometry2D->GetAxisVector(2);
   directionVector.Normalize();
   directionVector *= zSpacing;
 
   // Normally we should use the following four lines to create a copy of
   // the transform contrained in geometry2D, because it may not be changed
   // by us. But we know that SetSpacing creates a new transform without
   // changing the old (coming from geometry2D), so we can use the fifth
   // line instead. We check this at (**).
   //
   // AffineTransform3D::Pointer transform = AffineTransform3D::New();
   // transform->SetMatrix(geometry2D->GetIndexToWorldTransform()->GetMatrix());
   // transform->SetOffset(geometry2D->GetIndexToWorldTransform()->GetOffset());
   // SetIndexToWorldTransform(transform);
 
   this->SetIndexToWorldTransform(const_cast<AffineTransform3D *>(geometry2D->GetIndexToWorldTransform()));
 
   mitk::Vector3D spacing;
   FillVector3D(spacing, geometry2D->GetExtentInMM(0) / bounds[1], geometry2D->GetExtentInMM(1) / bounds[3], zSpacing);
 
   this->SetDirectionVector(directionVector);
   this->SetBounds(bounds);
   this->SetPlaneGeometry(geometry2D, 0);
   this->SetSpacing(spacing, true);
   this->SetEvenlySpaced();
 
   // this->SetTimeBounds( geometry2D->GetTimeBounds() );
 
   assert(this->GetIndexToWorldTransform() != geometry2D->GetIndexToWorldTransform()); // (**) see above.
 
   this->SetFrameOfReferenceID(geometry2D->GetFrameOfReferenceID());
   this->SetImageGeometry(geometry2D->GetImageGeometry());
 
   geometry2D->UnRegister();
 }
 
 void mitk::SlicedGeometry3D::InitializePlanes(const mitk::BaseGeometry *geometry3D,
                                               mitk::PlaneGeometry::PlaneOrientation planeorientation,
                                               bool top,
                                               bool frontside,
                                               bool rotated)
 {
   m_ReferenceGeometry = geometry3D;
 
   PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New();
   planeGeometry->InitializeStandardPlane(geometry3D, top, planeorientation, frontside, rotated);
 
   int worldAxis =
       planeorientation == PlaneGeometry::Sagittal ? 0 :
       planeorientation == PlaneGeometry::Frontal  ? 1 : 2;
 
-  // Inspired by:
-  // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3
+  int axes[3];
+  geometry3D->MapAxesToOrientations(axes);
+  int axis = axes[worldAxis];
 
+  ScalarType viewSpacing = geometry3D->GetSpacing()[axis];
+  unsigned int slices = static_cast<unsigned int>(geometry3D->GetExtent(axis));
+
+#ifndef NDEBUG
   mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix();
   matrix.GetVnlMatrix().normalize_columns();
   mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse();
 
-  int dominantAxis = itk::Function::Max3(
-      inverseMatrix[0][worldAxis],
-      inverseMatrix[1][worldAxis],
-      inverseMatrix[2][worldAxis]);
-
-  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.
-  unsigned int slices = static_cast<unsigned int>(geometry3D->GetExtent(dominantAxis) + 0.5);
-
-#ifndef NDEBUG
-  int upDirection = itk::Function::Sign(inverseMatrix[dominantAxis][worldAxis]);
+  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);
 
   int referencePointSlice = static_cast<int>(referencePointDistance / spacing[2]);
 
   double alignmentValue = referencePointDistance / spacing[2] - referencePointSlice;
 
   firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * alignmentValue * spacing[2]);
 
   // Finally, we can clear the previous geometry stack and initialize it with
   // our re-initialized "first plane".
   m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr));
 
   if (m_Slices > 0)
   {
     m_PlaneGeometries[0] = firstPlane;
   }
 
   // Reinitialize SNC with new number of slices
   m_SliceNavigationController->GetSlice()->SetSteps(m_Slices);
 
   this->Modified();
 }
 
 double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &d) const
 {
   // Need the spacing of the underlying dataset / geometry
   if (!m_ReferenceGeometry)
   {
     return 1.0;
   }
 
   const mitk::Vector3D &spacing = m_ReferenceGeometry->GetSpacing();
   return SlicedGeometry3D::CalculateSpacing(spacing, d);
 }
 
 double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &spacing, const mitk::Vector3D &d)
 {
   // The following can be derived from the ellipsoid equation
   //
   //   1 = x^2/a^2 + y^2/b^2 + z^2/c^2
   //
   // where (a,b,c) = spacing of original volume (ellipsoid radii)
   // and   (x,y,z) = scaled coordinates of vector d (according to ellipsoid)
   //
   double scaling = d[0] * d[0] / (spacing[0] * spacing[0]) + d[1] * d[1] / (spacing[1] * spacing[1]) +
                    d[2] * d[2] / (spacing[2] * spacing[2]);
 
   scaling = sqrt(scaling);
 
   return (sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]) / scaling);
 }
 
 mitk::Vector3D mitk::SlicedGeometry3D::AdjustNormal(const mitk::Vector3D &normal) const
 {
   TransformType::Pointer inverse = TransformType::New();
   m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse(inverse);
 
   Vector3D transformedNormal = inverse->TransformVector(normal);
 
   transformedNormal.Normalize();
   return transformedNormal;
 }
 
 void mitk::SlicedGeometry3D::SetImageGeometry(const bool isAnImageGeometry)
 {
   Superclass::SetImageGeometry(isAnImageGeometry);
 
   unsigned int s;
   for (s = 0; s < m_Slices; ++s)
   {
     mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
     if (geometry != nullptr)
     {
       geometry->SetImageGeometry(isAnImageGeometry);
     }
   }
 }
 
 void mitk::SlicedGeometry3D::ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry)
 {
   unsigned int s;
   for (s = 0; s < m_Slices; ++s)
   {
     mitk::BaseGeometry *geometry = m_PlaneGeometries[s];
     if (geometry != nullptr)
     {
       geometry->ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry);
     }
   }
 
   Superclass::ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry);
 }
 
 bool mitk::SlicedGeometry3D::IsValidSlice(int s) const
 {
   return ((s >= 0) && (s < (int)m_Slices));
 }
 
 const mitk::BaseGeometry *mitk::SlicedGeometry3D::GetReferenceGeometry() const
 {
   return m_ReferenceGeometry;
 }
 
 void mitk::SlicedGeometry3D::SetReferenceGeometry(const BaseGeometry *referenceGeometry)
 {
   m_ReferenceGeometry = referenceGeometry;
 
   std::vector<PlaneGeometry::Pointer>::iterator it;
 
   for (it = m_PlaneGeometries.begin(); it != m_PlaneGeometries.end(); ++it)
   {
     (*it)->SetReferenceGeometry(referenceGeometry);
   }
 }
 
 bool mitk::SlicedGeometry3D::HasReferenceGeometry() const
 {
   return ( m_ReferenceGeometry != nullptr );
 }
 
 void mitk::SlicedGeometry3D::PreSetSpacing(const mitk::Vector3D &aSpacing)
 {
   bool hasEvenlySpacedPlaneGeometry = false;
   mitk::Point3D origin;
   mitk::Vector3D rightDV, bottomDV;
   BoundingBox::BoundsArrayType bounds;
 
   // Check for valid spacing
   if (!(aSpacing[0] > 0 && aSpacing[1] > 0 && aSpacing[2] > 0))
   {
     mitkThrow() << "You try to set a spacing with at least one element equal or "
                    "smaller to \"0\". This might lead to a crash during rendering. Please double"
                    " check your data!";
   }
 
   // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
   // since the spacing influences them
   if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
   {
     const PlaneGeometry *planeGeometry = m_PlaneGeometries[0];
 
     if (planeGeometry && !dynamic_cast<const AbstractTransformGeometry *>(planeGeometry))
     {
       this->WorldToIndex(planeGeometry->GetOrigin(), origin);
       this->WorldToIndex(planeGeometry->GetAxisVector(0), rightDV);
       this->WorldToIndex(planeGeometry->GetAxisVector(1), bottomDV);
 
       bounds = planeGeometry->GetBounds();
       hasEvenlySpacedPlaneGeometry = true;
     }
   }
 
   BaseGeometry::_SetSpacing(aSpacing);
 
   mitk::PlaneGeometry::Pointer firstGeometry;
 
   // In case of evenly-spaced data: re-initialize instances of PlaneGeometry,
   // since the spacing influences them
   if (hasEvenlySpacedPlaneGeometry)
   {
     // create planeGeometry according to new spacing
     this->IndexToWorld(origin, origin);
     this->IndexToWorld(rightDV, rightDV);
     this->IndexToWorld(bottomDV, bottomDV);
 
     mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New();
     planeGeometry->SetImageGeometry(this->GetImageGeometry());
 
     planeGeometry->SetReferenceGeometry(m_ReferenceGeometry);
 
     // Store spacing, as Initialize... needs a pointer
     mitk::Vector3D lokalSpacing = this->GetSpacing();
     planeGeometry->InitializeStandardPlane(rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &lokalSpacing);
     planeGeometry->SetOrigin(origin);
     planeGeometry->SetBounds(bounds);
 
     firstGeometry = planeGeometry;
   }
   else if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0))
   {
     firstGeometry = m_PlaneGeometries[0].GetPointer();
   }
 
   // clear and reserve
   PlaneGeometry::Pointer gnull = nullptr;
   m_PlaneGeometries.assign(m_Slices, gnull);
 
   if (m_Slices > 0)
   {
     m_PlaneGeometries[0] = firstGeometry;
   }
 
   this->Modified();
 }
 
 void mitk::SlicedGeometry3D::SetSliceNavigationController(SliceNavigationController *snc)
 {
   m_SliceNavigationController = snc;
 }
 
 mitk::SliceNavigationController *mitk::SlicedGeometry3D::GetSliceNavigationController()
 {
   return m_SliceNavigationController;
 }
 
 void mitk::SlicedGeometry3D::SetEvenlySpaced(bool on)
 {
   if (m_EvenlySpaced != on)
   {
     m_EvenlySpaced = on;
     this->Modified();
   }
 }
 
 void mitk::SlicedGeometry3D::SetDirectionVector(const mitk::Vector3D &directionVector)
 {
   Vector3D newDir = directionVector;
   newDir.Normalize();
   if (newDir != m_DirectionVector)
   {
     m_DirectionVector = newDir;
     this->Modified();
   }
 }
 
 // void
 // mitk::SlicedGeometry3D::SetTimeBounds( const mitk::TimeBounds& timebounds )
 //{
 //  Superclass::SetTimeBounds( timebounds );
 //
 //  unsigned int s;
 //  for ( s = 0; s < m_Slices; ++s )
 //  {
 //    if(m_Geometry2Ds[s].IsNotNull())
 //    {
 //      m_Geometry2Ds[s]->SetTimeBounds( timebounds );
 //    }
 //  }
 //  m_TimeBounds = timebounds;
 //}
 
 itk::LightObject::Pointer mitk::SlicedGeometry3D::InternalClone() const
 {
   Self::Pointer newGeometry = new SlicedGeometry3D(*this);
   newGeometry->UnRegister();
   return newGeometry.GetPointer();
 }
 
 void mitk::SlicedGeometry3D::PrintSelf(std::ostream &os, itk::Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
   os << indent << " EvenlySpaced: " << m_EvenlySpaced << std::endl;
   if (m_EvenlySpaced)
   {
     os << indent << " DirectionVector: " << m_DirectionVector << std::endl;
   }
   os << indent << " Slices: " << m_Slices << std::endl;
 
   os << std::endl;
   os << indent << " GetPlaneGeometry(0): ";
   if (this->GetPlaneGeometry(0) == nullptr)
   {
     os << "nullptr" << std::endl;
   }
   else
   {
     this->GetPlaneGeometry(0)->Print(os, indent);
   }
 }
 
 void mitk::SlicedGeometry3D::ExecuteOperation(Operation *operation)
 {
   PlaneGeometry::Pointer geometry2D;
   ApplyTransformMatrixOperation *applyMatrixOp;
   Point3D center;
 
   switch (operation->GetOperationType())
   {
     case OpNOTHING:
       break;
 
     case OpROTATE:
       if (m_EvenlySpaced)
       {
         // Need a reference frame to align the rotation
         if (m_ReferenceGeometry)
         {
           // Clear all generated geometries and then rotate only the first slice.
           // The other slices will be re-generated on demand
 
           // Save first slice
           PlaneGeometry::Pointer geometry2D = m_PlaneGeometries[0];
 
           RotationOperation *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
             RotationOperation *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
         PlaneOperation *planeOp = dynamic_cast<PlaneOperation *>(operation);
 
         // Get first slice
         PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
 
         // Need a PlaneGeometry, a PlaneOperation and a reference frame to
         // carry out the re-orientation. If not all avaialble, stop here
         if (!m_ReferenceGeometry ||
             (!planeGeometry || dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer())) || !planeOp)
         {
           break;
         }
 
         // General Behavior:
         // Clear all generated geometries and then rotate only the first slice.
         // The other slices will be re-generated on demand
 
         //
         // 1st Step: Reorient Normal Vector of first plane
         //
         Point3D center = planeOp->GetPoint(); // m_ReferenceGeometry->GetCenter();
         mitk::Vector3D currentNormal = planeGeometry->GetNormal();
         mitk::Vector3D newNormal;
         if (planeOp->AreAxisDefined())
         {
           // If planeOp was defined by one centerpoint and two axis vectors
           newNormal = CrossProduct(planeOp->GetAxisVec0(), planeOp->GetAxisVec1());
         }
         else
         {
           // If planeOp was defined by one centerpoint and one normal vector
           newNormal = planeOp->GetNormal();
         }
 
         // Get Rotation axis und angle
         currentNormal.Normalize();
         newNormal.Normalize();
         ScalarType rotationAngle = angle(currentNormal.GetVnlVector(), newNormal.GetVnlVector());
 
         rotationAngle *= 180.0 / vnl_math::pi; // from rad to deg
         Vector3D rotationAxis = itk::CrossProduct(currentNormal, newNormal);
         if (std::abs(rotationAngle - 180) < mitk::eps)
         {
           // current Normal and desired normal are not linear independent!!(e.g 1,0,0 and -1,0,0).
           // Rotation Axis should be ANY vector that is 90� to current Normal
           mitk::Vector3D helpNormal;
           helpNormal = currentNormal;
           helpNormal[0] += 1;
           helpNormal[1] -= 1;
           helpNormal[2] += 1;
           helpNormal.Normalize();
           rotationAxis = itk::CrossProduct(helpNormal, currentNormal);
         }
 
         RotationOperation centeredRotation(mitk::OpROTATE, center, rotationAxis, rotationAngle);
 
         // Rotate first slice
         planeGeometry->ExecuteOperation(&centeredRotation);
 
         // Reinitialize planes and select slice, if my rotations are all done.
         if (!planeOp->AreAxisDefined())
         {
           // Clear the slice stack and adjust it according to the center of
           // rotation and plane position (see documentation of ReinitializePlanes)
           this->ReinitializePlanes(center, planeOp->GetPoint());
           planeGeometry->SetSpacing(this->GetSpacing());
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
         }
 
         // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
         BaseGeometry::ExecuteOperation(&centeredRotation);
 
         //
         // 2nd step. If axis vectors were defined, rotate the plane around its normal to fit these
         //
 
         if (planeOp->AreAxisDefined())
         {
           mitk::Vector3D vecAxixNew = planeOp->GetAxisVec0();
           vecAxixNew.Normalize();
           mitk::Vector3D VecAxisCurr = planeGeometry->GetAxisVector(0);
           VecAxisCurr.Normalize();
 
           ScalarType rotationAngle = angle(VecAxisCurr.GetVnlVector(), vecAxixNew.GetVnlVector());
           rotationAngle = rotationAngle * 180 / PI; // Rad to Deg
 
           // we rotate around the normal of the plane, but we do not know, if we need to rotate clockwise
           // or anti-clockwise. So we rotate around the crossproduct of old and new Axisvector.
           // Since both axis vectors lie in the plane, the crossproduct is the planes normal or the negative planes
           // normal
 
           rotationAxis = itk::CrossProduct(VecAxisCurr, vecAxixNew);
           if (std::abs(rotationAngle - 180) < mitk::eps)
           {
             // current axisVec and desired axisVec are not linear independent!!(e.g 1,0,0 and -1,0,0).
             // Rotation Axis can be just plane Normal. (have to rotate by 180�)
             rotationAxis = newNormal;
           }
 
           // Perfom Rotation
           mitk::RotationOperation op(mitk::OpROTATE, center, rotationAxis, rotationAngle);
           planeGeometry->ExecuteOperation(&op);
 
           // Apply changes on first slice to whole slice stack
           this->ReinitializePlanes(center, planeOp->GetPoint());
           planeGeometry->SetSpacing(this->GetSpacing());
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
 
           // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry)
           BaseGeometry::ExecuteOperation(&op);
         }
       }
       else
       {
         // Reach through to all slices
         for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
         {
           (*iter)->ExecuteOperation(operation);
         }
       }
       break;
 
     case OpRESTOREPLANEPOSITION:
       if (m_EvenlySpaced)
       {
         // Save first slice
         PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0];
 
         RestorePlanePositionOperation *restorePlaneOp = dynamic_cast<RestorePlanePositionOperation *>(operation);
 
         // Need a PlaneGeometry, a PlaneOperation and a reference frame to
         // carry out the re-orientation
         if (m_ReferenceGeometry &&
             (planeGeometry && dynamic_cast<AbstractTransformGeometry *>(planeGeometry.GetPointer()) == nullptr) &&
             restorePlaneOp)
         {
           // Clear all generated geometries and then rotate only the first slice.
           // The other slices will be re-generated on demand
 
           // Rotate first slice
           planeGeometry->ExecuteOperation(restorePlaneOp);
 
           m_DirectionVector = restorePlaneOp->GetDirectionVector();
 
           double centerOfRotationDistance = planeGeometry->SignedDistanceFromPlane(m_ReferenceGeometry->GetCenter());
 
           if (centerOfRotationDistance <= 0)
           {
             m_DirectionVector = -m_DirectionVector;
           }
 
           Vector3D spacing = restorePlaneOp->GetSpacing();
 
           Superclass::SetSpacing(spacing);
 
           // /*Now we need to calculate the number of slices in the plane's normal
           // direction, so that the entire volume is covered. This is done by first
           // calculating the dot product between the volume diagonal (the maximum
           // distance inside the volume) and the normal, and dividing this value by
           // the directed spacing calculated above.*/
           ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * m_DirectionVector[0]) +
                                       std::abs(m_ReferenceGeometry->GetExtentInMM(1) * m_DirectionVector[1]) +
                                       std::abs(m_ReferenceGeometry->GetExtentInMM(2) * m_DirectionVector[2]);
 
           if (directedExtent >= spacing[2])
           {
             m_Slices = static_cast<unsigned int>(directedExtent / spacing[2] + 0.5);
           }
           else
           {
             m_Slices = 1;
           }
 
           m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr));
 
           if (m_Slices > 0)
           {
             m_PlaneGeometries[0] = planeGeometry;
           }
 
           m_SliceNavigationController->GetSlice()->SetSteps(m_Slices);
 
           this->Modified();
 
           // End Reinitialization
 
           if (m_SliceNavigationController)
           {
             m_SliceNavigationController->GetSlice()->SetPos(restorePlaneOp->GetPos());
             m_SliceNavigationController->AdjustSliceStepperRange();
           }
           BaseGeometry::ExecuteOperation(restorePlaneOp);
         }
       }
       else
       {
         // Reach through to all slices
         for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter)
         {
           (*iter)->ExecuteOperation(operation);
         }
       }
       break;
 
     case OpAPPLYTRANSFORMMATRIX:
 
       // Clear all generated geometries and then transform only the first slice.
       // The other slices will be re-generated on demand
 
       // Save first slice
       geometry2D = m_PlaneGeometries[0];
 
       applyMatrixOp = dynamic_cast<ApplyTransformMatrixOperation *>(operation);
 
       // Apply transformation to first plane
       geometry2D->ExecuteOperation(applyMatrixOp);
 
       // Generate a ApplyTransformMatrixOperation using the dataset center instead of
       // the supplied rotation center. The supplied center is instead used to adjust the
       // slice stack afterwards (see OpROTATE).
       center = m_ReferenceGeometry->GetCenter();
 
       // Clear the slice stack and adjust it according to the center of
       // the dataset and the supplied rotation center (see documentation of
       // ReinitializePlanes)
       this->ReinitializePlanes(center, applyMatrixOp->GetReferencePoint());
 
       BaseGeometry::ExecuteOperation(applyMatrixOp);
       break;
 
     default: // let handle by base class if we don't do anything
       BaseGeometry::ExecuteOperation(operation);
   }
 
   this->Modified();
 }