# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1558018066 -7200
# Node ID 7719eb852dd547289704501608d30245bbf28463
# Parent  dbc1d8bfc68a51652a5047d9ac7cc85538456e1a
new class: VolumeImageGeometry

diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Toolbox/FiniteProjectiveCamera.cpp
--- a/Framework/Toolbox/FiniteProjectiveCamera.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Toolbox/FiniteProjectiveCamera.cpp	Thu May 16 16:47:46 2019 +0200
@@ -316,7 +316,7 @@
     LOG(WARNING) << "Output pixel format: " << Orthanc::EnumerationToString(target.GetFormat());
 
     std::auto_ptr<OrthancStone::ParallelSlices> slices(source.GetGeometry(projection));
-    const OrthancStone::Vector pixelSpacing = source.GetVoxelDimensions(projection);
+    const OrthancStone::Vector pixelSpacing = source.GetGeometry().GetVoxelDimensions(projection);
     const unsigned int targetWidth = target.GetWidth();
     const unsigned int targetHeight = target.GetHeight();
 
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Toolbox/OrientedBoundingBox.cpp
--- a/Framework/Toolbox/OrientedBoundingBox.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Toolbox/OrientedBoundingBox.cpp	Thu May 16 16:47:46 2019 +0200
@@ -37,8 +37,8 @@
       throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);      
     }
 
-    const CoordinateSystem3D& geometry = image.GetAxialGeometry();
-    Vector dim = image.GetVoxelDimensions(VolumeProjection_Axial);
+    const CoordinateSystem3D& geometry = image.GetGeometry().GetAxialGeometry();
+    Vector dim = image.GetGeometry().GetVoxelDimensions(VolumeProjection_Axial);
 
     u_ = geometry.GetAxisX();
     v_ = geometry.GetAxisY();
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Toolbox/ShearWarpProjectiveTransform.cpp
--- a/Framework/Toolbox/ShearWarpProjectiveTransform.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Toolbox/ShearWarpProjectiveTransform.cpp	Thu May 16 16:47:46 2019 +0200
@@ -380,8 +380,8 @@
     
     // Compute the "world" matrix that maps the source volume to the
     // (0,0,0)->(1,1,1) unit cube
-    Vector origin = source.GetCoordinates(0, 0, 0);
-    Vector ps = source.GetVoxelDimensions(VolumeProjection_Axial);
+    Vector origin = source.GetGeometry().GetCoordinates(0, 0, 0);
+    Vector ps = source.GetGeometry().GetVoxelDimensions(VolumeProjection_Axial);
     Matrix world = LinearAlgebra::Product(
       GeometryToolbox::CreateScalingMatrix(1.0 / ps[0], 1.0 / ps[1], 1.0 / ps[2]),
       GeometryToolbox::CreateTranslationMatrix(-origin[0], -origin[1], -origin[2]));
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Toolbox/VolumeImageGeometry.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/VolumeImageGeometry.cpp	Thu May 16 16:47:46 2019 +0200
@@ -0,0 +1,205 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2019 Osimis S.A., Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Affero General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#include "VolumeImageGeometry.h"
+
+#include "../Toolbox/GeometryToolbox.h"
+
+#include <Core/OrthancException.h>
+
+
+namespace OrthancStone
+{
+  void VolumeImageGeometry::Invalidate()
+  {
+    Vector p = (axialGeometry_.GetOrigin() +
+                static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal());
+        
+    coronalGeometry_ = CoordinateSystem3D(p,
+                                          axialGeometry_.GetAxisX(),
+                                          -axialGeometry_.GetNormal());
+    
+    sagittalGeometry_ = CoordinateSystem3D(p,
+                                           axialGeometry_.GetAxisY(),
+                                           -axialGeometry_.GetNormal());
+
+    Vector origin = (
+      axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],
+                                                -0.5 * voxelDimensions_[1]) -
+      0.5 * voxelDimensions_[2] * axialGeometry_.GetNormal());
+
+    Vector scaling;
+    
+    if (width_ == 0 ||
+        height_ == 0 ||
+        depth_ == 0)
+    {
+      LinearAlgebra::AssignVector(scaling, 1, 1, 1);
+    }
+    else
+    {
+      scaling = (
+        axialGeometry_.GetAxisX() * voxelDimensions_[0] * static_cast<double>(width_) +
+        axialGeometry_.GetAxisY() * voxelDimensions_[1] * static_cast<double>(height_) +
+        axialGeometry_.GetNormal() * voxelDimensions_[2] * static_cast<double>(depth_));
+    }
+
+    transform_ = LinearAlgebra::Product(
+      GeometryToolbox::CreateTranslationMatrix(origin[0], origin[1], origin[2]),
+      GeometryToolbox::CreateScalingMatrix(scaling[0], scaling[1], scaling[2]));
+
+    LinearAlgebra::InvertMatrix(transformInverse_, transform_);
+  }
+
+  
+  VolumeImageGeometry::VolumeImageGeometry() :
+    width_(0),
+    height_(0),
+    depth_(0)
+  {
+    LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1);
+    Invalidate();
+  }
+
+
+  void VolumeImageGeometry::SetSize(unsigned int width,
+                                    unsigned int height,
+                                    unsigned int depth)
+  {
+    width_ = width;
+    height_ = height;
+    depth_ = depth;
+    Invalidate();
+  }
+
+  
+  void VolumeImageGeometry::SetAxialGeometry(const CoordinateSystem3D& geometry)
+  {
+    axialGeometry_ = geometry;
+    Invalidate();
+  }
+
+
+  void VolumeImageGeometry::SetVoxelDimensions(double x,
+                                               double y,
+                                               double z)
+  {
+    if (x <= 0 ||
+        y <= 0 ||
+        z <= 0)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+    else
+    {
+      LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
+      Invalidate();
+    }
+  }
+
+
+  Vector VolumeImageGeometry::GetVoxelDimensions(VolumeProjection projection) const
+  {
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        return voxelDimensions_;
+
+      case VolumeProjection_Coronal:
+        return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
+
+      case VolumeProjection_Sagittal:
+        return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+  }
+
+
+  void VolumeImageGeometry::GetSliceSize(unsigned int& width,
+                                         unsigned int& height,
+                                         VolumeProjection projection)
+  {
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        width = width_;
+        height = height_;
+        break;
+
+      case VolumeProjection_Coronal:
+        width = width_;
+        height = depth_;
+        break;
+
+      case VolumeProjection_Sagittal:
+        width = height_;
+        height = depth_;
+        break;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+  }
+
+
+
+  Vector VolumeImageGeometry::GetCoordinates(float x,
+                                             float y,
+                                             float z) const
+  {
+    Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
+
+    assert(LinearAlgebra::IsNear(p[3], 1));  // Affine transform, no perspective effect
+
+    // Back to non-homogeneous coordinates
+    return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
+  }
+
+
+  bool VolumeImageGeometry::DetectProjection(VolumeProjection& projection,
+                                             const CoordinateSystem3D& plane) const
+  {
+    if (GeometryToolbox::IsParallel(plane.GetNormal(),
+                                    axialGeometry_.GetNormal()))
+    {
+      projection = VolumeProjection_Axial;
+      return true;
+    }
+    else if (GeometryToolbox::IsParallel(plane.GetNormal(),
+                                         coronalGeometry_.GetNormal()))
+    {
+      projection = VolumeProjection_Coronal;
+      return true;
+    }
+    else if (GeometryToolbox::IsParallel(plane.GetNormal(),
+                                         sagittalGeometry_.GetNormal()))
+    {
+      projection = VolumeProjection_Sagittal;
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+}
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Toolbox/VolumeImageGeometry.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/VolumeImageGeometry.h	Thu May 16 16:47:46 2019 +0200
@@ -0,0 +1,104 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, Belgium
+ * Copyright (C) 2017-2019 Osimis S.A., Belgium
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU Affero General Public License
+ * as published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ **/
+
+
+#pragma once
+
+#include "../StoneEnumerations.h"
+#include "CoordinateSystem3D.h"
+
+namespace OrthancStone
+{
+  class VolumeImageGeometry
+  {
+  private:
+    unsigned int           width_;
+    unsigned int           height_;
+    unsigned int           depth_;
+    CoordinateSystem3D     axialGeometry_;
+    CoordinateSystem3D     coronalGeometry_;
+    CoordinateSystem3D     sagittalGeometry_;
+    Vector                 voxelDimensions_;
+    Matrix                 transform_;
+    Matrix                 transformInverse_;
+
+    void Invalidate();
+
+  public:
+    VolumeImageGeometry();
+
+    unsigned int GetWidth() const
+    {
+      return width_;
+    }
+
+    unsigned int GetHeight() const
+    {
+      return height_;
+    }
+
+    unsigned int GetDepth() const
+    {
+      return depth_;
+    }
+
+    const CoordinateSystem3D& GetAxialGeometry() const
+    {
+      return axialGeometry_;
+    }
+
+    const CoordinateSystem3D& GetCoronalGeometry() const
+    {
+      return coronalGeometry_;
+    }
+
+    const CoordinateSystem3D& GetSagittalGeometry() const
+    {
+      return sagittalGeometry_;
+    }
+
+    void SetSize(unsigned int width,
+                 unsigned int height,
+                 unsigned int depth);
+
+    // Set the geometry of the first axial slice (i.e. the one whose
+    // depth == 0)
+    void SetAxialGeometry(const CoordinateSystem3D& geometry);
+
+    void SetVoxelDimensions(double x,
+                            double y,
+                            double z);
+
+    Vector GetVoxelDimensions(VolumeProjection projection) const;
+
+    void GetSliceSize(unsigned int& width,
+                      unsigned int& height,
+                      VolumeProjection projection);
+    
+    // Get the 3D position of a point in the volume, where x, y and z
+    // lie in the [0;1] range
+    Vector GetCoordinates(float x,
+                          float y,
+                          float z) const;
+
+    bool DetectProjection(VolumeProjection& projection,
+                          const CoordinateSystem3D& plane) const;
+  };
+}
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Volumes/ImageBuffer3D.cpp
--- a/Framework/Volumes/ImageBuffer3D.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Volumes/ImageBuffer3D.cpp	Thu May 16 16:47:46 2019 +0200
@@ -118,13 +118,11 @@
     computeRange_(computeRange),
     hasRange_(false)
   {
-    LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1);
+    geometry_.SetSize(width, height, depth);
 
     LOG(INFO) << "Created a 3D image of size " << width << "x" << height
               << "x" << depth << " in " << Orthanc::EnumerationToString(format)
               << " (" << (GetEstimatedMemorySize() / (1024ll * 1024ll)) << "MB)";
-
-    UpdateGeometry();
   }
 
 
@@ -134,123 +132,57 @@
   }
 
 
-  void ImageBuffer3D::SetAxialGeometry(const CoordinateSystem3D& geometry)
-  {
-    axialGeometry_ = geometry;
-    UpdateGeometry();
-  }
-
-
-  void ImageBuffer3D::SetVoxelDimensions(double x,
-                                         double y,
-                                         double z)
-  {
-    if (x <= 0 ||
-        y <= 0 ||
-        z <= 0)
-    {
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-    }
-    else
-    {
-      LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
-      UpdateGeometry();
-    }
-  }
-
-
-  Vector ImageBuffer3D::GetVoxelDimensions(VolumeProjection projection) const
-  {
-    switch (projection)
-    {
-    case VolumeProjection_Axial:
-      return voxelDimensions_;
-
-    case VolumeProjection_Coronal:
-      return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
-
-    case VolumeProjection_Sagittal:
-      return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
-
-    default:
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-    }
-  }
-
-
-  void ImageBuffer3D::GetSliceSize(unsigned int& width,
-                                   unsigned int& height,
-                                   VolumeProjection projection)
-  {
-    switch (projection)
-    {
-    case VolumeProjection_Axial:
-      width = width_;
-      height = height_;
-      break;
-
-    case VolumeProjection_Coronal:
-      width = width_;
-      height = depth_;
-      break;
-
-    case VolumeProjection_Sagittal:
-      width = height_;
-      height = depth_;
-      break;
-
-    default:
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
-    }
-  }
 
 
   ParallelSlices* ImageBuffer3D::GetGeometry(VolumeProjection projection) const
   {
+    const Vector dimensions = geometry_.GetVoxelDimensions(VolumeProjection_Axial);
+    const CoordinateSystem3D& axial = geometry_.GetAxialGeometry();
+    
     std::auto_ptr<ParallelSlices> result(new ParallelSlices);
 
     switch (projection)
     {
-    case VolumeProjection_Axial:
-      for (unsigned int z = 0; z < depth_; z++)
-      {
-        Vector origin = axialGeometry_.GetOrigin();
-        origin += static_cast<double>(z) * voxelDimensions_[2] * axialGeometry_.GetNormal();
+      case VolumeProjection_Axial:
+        for (unsigned int z = 0; z < depth_; z++)
+        {
+          Vector origin = axial.GetOrigin();
+          origin += static_cast<double>(z) * dimensions[2] * axial.GetNormal();
 
-        result->AddSlice(origin,
-                         axialGeometry_.GetAxisX(),
-                         axialGeometry_.GetAxisY());
-      }
-      break;
+          result->AddSlice(origin,
+                           axial.GetAxisX(),
+                           axial.GetAxisY());
+        }
+        break;
 
-    case VolumeProjection_Coronal:
-      for (unsigned int y = 0; y < height_; y++)
-      {
-        Vector origin = axialGeometry_.GetOrigin();
-        origin += static_cast<double>(y) * voxelDimensions_[1] * axialGeometry_.GetAxisY();
-        origin += static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal();
+      case VolumeProjection_Coronal:
+        for (unsigned int y = 0; y < height_; y++)
+        {
+          Vector origin = axial.GetOrigin();
+          origin += static_cast<double>(y) * dimensions[1] * axial.GetAxisY();
+          origin += static_cast<double>(depth_ - 1) * dimensions[2] * axial.GetNormal();
 
-        result->AddSlice(origin,
-                         axialGeometry_.GetAxisX(),
-                         -axialGeometry_.GetNormal());
-      }
-      break;
+          result->AddSlice(origin,
+                           axial.GetAxisX(),
+                           -axial.GetNormal());
+        }
+        break;
 
-    case VolumeProjection_Sagittal:
-      for (unsigned int x = 0; x < width_; x++)
-      {
-        Vector origin = axialGeometry_.GetOrigin();
-        origin += static_cast<double>(x) * voxelDimensions_[0] * axialGeometry_.GetAxisX();
-        origin += static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal();
+      case VolumeProjection_Sagittal:
+        for (unsigned int x = 0; x < width_; x++)
+        {
+          Vector origin = axial.GetOrigin();
+          origin += static_cast<double>(x) * dimensions[0] * axial.GetAxisX();
+          origin += static_cast<double>(depth_ - 1) * dimensions[2] * axial.GetNormal();
 
-        result->AddSlice(origin,
-                         axialGeometry_.GetAxisY(),
-                         -axialGeometry_.GetNormal());
-      }
-      break;
+          result->AddSlice(origin,
+                           axial.GetAxisY(),
+                           -axial.GetNormal());
+        }
+        break;
 
-    default:
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
     }
 
     return result.release();
@@ -276,24 +208,24 @@
 
     switch (slice.GetFormat())
     {
-    case Orthanc::PixelFormat_Grayscale8:
-    case Orthanc::PixelFormat_Grayscale16:
-    case Orthanc::PixelFormat_Grayscale32:
-    case Orthanc::PixelFormat_SignedGrayscale16:
-    {
-      int64_t a, b;
-      Orthanc::ImageProcessing::GetMinMaxIntegerValue(a, b, slice);
-      sliceMin = static_cast<float>(a);
-      sliceMax = static_cast<float>(b);
-      break;
-    }
+      case Orthanc::PixelFormat_Grayscale8:
+      case Orthanc::PixelFormat_Grayscale16:
+      case Orthanc::PixelFormat_Grayscale32:
+      case Orthanc::PixelFormat_SignedGrayscale16:
+      {
+        int64_t a, b;
+        Orthanc::ImageProcessing::GetMinMaxIntegerValue(a, b, slice);
+        sliceMin = static_cast<float>(a);
+        sliceMax = static_cast<float>(b);
+        break;
+      }
 
-    case Orthanc::PixelFormat_Float32:
-      Orthanc::ImageProcessing::GetMinMaxFloatValue(sliceMin, sliceMax, slice);
-      break;
+      case Orthanc::PixelFormat_Float32:
+        Orthanc::ImageProcessing::GetMinMaxFloatValue(sliceMin, sliceMax, slice);
+        break;
 
-    default:
-      return;
+      default:
+        return;
     }
 
     if (hasRange_)
@@ -367,8 +299,8 @@
         sagittal_->GetReadOnlyAccessor(accessor_);
         break;
 
-    default:
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
     }
   }
 
@@ -411,8 +343,8 @@
         sagittal_->GetWriteableAccessor(accessor_);
         break;
 
-    default:
-      throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
     }
   }
 
@@ -457,80 +389,4 @@
     const void* p = image_.GetConstRow(y + height_ * (depth_ - 1 - z));
     return reinterpret_cast<const uint16_t*>(p) [x];
   }
-
-
-  void ImageBuffer3D::UpdateGeometry()
-  {
-    Vector p = (axialGeometry_.GetOrigin() +
-                static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal());
-        
-    coronalGeometry_ = CoordinateSystem3D(p,
-                                          axialGeometry_.GetAxisX(),
-                                          -axialGeometry_.GetNormal());
-    
-    sagittalGeometry_ = CoordinateSystem3D(p,
-                                           axialGeometry_.GetAxisY(),
-                                           -axialGeometry_.GetNormal());
-
-    Vector origin = (
-      axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],
-                                                -0.5 * voxelDimensions_[1]) -
-      0.5 * voxelDimensions_[2] * axialGeometry_.GetNormal());
-
-    Vector scaling = (
-      axialGeometry_.GetAxisX() * voxelDimensions_[0] * static_cast<double>(width_) +
-      axialGeometry_.GetAxisY() * voxelDimensions_[1] * static_cast<double>(height_) +
-      axialGeometry_.GetNormal() * voxelDimensions_[2] * static_cast<double>(depth_));
-
-    transform_ = LinearAlgebra::Product(
-      GeometryToolbox::CreateTranslationMatrix(origin[0], origin[1], origin[2]),
-      GeometryToolbox::CreateScalingMatrix(scaling[0], scaling[1], scaling[2]));
-
-    LinearAlgebra::InvertMatrix(transformInverse_, transform_);
-  }
-
-  
-  Vector ImageBuffer3D::GetCoordinates(float x,
-                                       float y,
-                                       float z) const
-  {
-    Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
-
-    assert(LinearAlgebra::IsNear(p[3], 1));  // Affine transform, no perspective effect
-
-    // Back to non-homogeneous coordinates
-    return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
-  }
-
-
-  bool ImageBuffer3D::DetectProjection(VolumeProjection& projection,
-                                       const CoordinateSystem3D& plane) const
-  {
-    if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                    axialGeometry_.GetNormal()))
-    {
-      projection = VolumeProjection_Axial;
-      return true;
-    }
-
-    {
-      if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                      coronalGeometry_.GetNormal()))
-      {
-        projection = VolumeProjection_Coronal;
-        return true;
-      }
-    }
-
-    {
-      if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                      sagittalGeometry_.GetNormal()))
-      {
-        projection = VolumeProjection_Sagittal;
-        return true;
-      }
-    }
-
-    return false;
-  }
 }
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Volumes/ImageBuffer3D.h
--- a/Framework/Volumes/ImageBuffer3D.h	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Volumes/ImageBuffer3D.h	Thu May 16 16:47:46 2019 +0200
@@ -23,7 +23,7 @@
 
 #include "../StoneEnumerations.h"
 #include "../Layers/RenderStyle.h"
-#include "../Toolbox/CoordinateSystem3D.h"
+#include "../Toolbox/VolumeImageGeometry.h"
 #include "../Toolbox/DicomFrameConverter.h"
 #include "../Toolbox/ParallelSlices.h"
 
@@ -34,10 +34,7 @@
   class ImageBuffer3D : public boost::noncopyable
   {
   private:
-    CoordinateSystem3D     axialGeometry_;
-    CoordinateSystem3D     coronalGeometry_;
-    CoordinateSystem3D     sagittalGeometry_;
-    Vector                 voxelDimensions_;
+    VolumeImageGeometry    geometry_;  // TODO => Move this out of this class
     Orthanc::Image         image_;
     Orthanc::PixelFormat   format_;
     unsigned int           width_;
@@ -50,8 +47,6 @@
     Matrix                 transform_;
     Matrix                 transformInverse_;
 
-    void UpdateGeometry();
-
     void ExtendImageRange(const Orthanc::ImageAccessor& slice);
 
     void GetAxialSliceAccessor(Orthanc::ImageAccessor& target,
@@ -83,35 +78,16 @@
 
     void Clear();
 
-    // Set the geometry of the first axial slice (i.e. the one whose
-    // depth == 0)
-    void SetAxialGeometry(const CoordinateSystem3D& geometry);
-
-    const CoordinateSystem3D& GetAxialGeometry() const
+    VolumeImageGeometry& GetGeometry()
     {
-      return axialGeometry_;
-    }
-
-    const CoordinateSystem3D& GetCoronalGeometry() const
-    {
-      return coronalGeometry_;
+      return geometry_;
     }
 
-    const CoordinateSystem3D& GetSagittalGeometry() const
+    const VolumeImageGeometry& GetGeometry() const
     {
-      return sagittalGeometry_;
+      return geometry_;
     }
 
-    void SetVoxelDimensions(double x,
-                            double y,
-                            double z);
-
-    Vector GetVoxelDimensions(VolumeProjection projection) const;
-
-    void GetSliceSize(unsigned int& width,
-                      unsigned int& height,
-                      VolumeProjection projection);
-
     const Orthanc::ImageAccessor& GetInternalImage() const
     {
       return image_;
@@ -177,16 +153,7 @@
                                  unsigned int y,
                                  unsigned int z) const;
 
-    // Get the 3D position of a point in the volume, where x, y and z
-    // lie in the [0;1] range
-    Vector GetCoordinates(float x,
-                          float y,
-                          float z) const;
-
-    bool DetectProjection(VolumeProjection& projection,
-                          const CoordinateSystem3D& plane) const;
-
-
+    
     class SliceReader : public boost::noncopyable
     {
     private:
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/Volumes/VolumeReslicer.cpp
--- a/Framework/Volumes/VolumeReslicer.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/Volumes/VolumeReslicer.cpp	Thu May 16 16:47:46 2019 +0200
@@ -749,7 +749,8 @@
   {
     // Choose the default voxel size as the finest voxel dimension
     // of the source volumetric image
-    const OrthancStone::Vector dim = source.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
+    const OrthancStone::Vector dim =
+      source.GetGeometry().GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
     double voxelSize = dim[0];
     
     if (dim[1] < voxelSize)
diff -r dbc1d8bfc68a -r 7719eb852dd5 Framework/dev.h
--- a/Framework/dev.h	Thu May 16 16:01:36 2019 +0200
+++ b/Framework/dev.h	Thu May 16 16:47:46 2019 +0200
@@ -157,9 +157,9 @@
                 << "x" << loader_.GetSlicesCount() << " in " << Orthanc::EnumerationToString(format);
 
       image_.reset(new ImageBuffer3D(format, width, height, static_cast<unsigned int>(loader_.GetSlicesCount()), computeRange_));
-      image_->SetAxialGeometry(loader_.GetSlice(0).GetGeometry());
-      image_->SetVoxelDimensions(loader_.GetSlice(0).GetPixelSpacingX(),
-                                 loader_.GetSlice(0).GetPixelSpacingY(), spacingZ);
+      image_->GetGeometry().SetAxialGeometry(loader_.GetSlice(0).GetGeometry());
+      image_->GetGeometry().SetVoxelDimensions(loader_.GetSlice(0).GetPixelSpacingX(),
+                                               loader_.GetSlice(0).GetPixelSpacingY(), spacingZ);
       image_->Clear();
 
       downloadStack_.reset(new DownloadStack(static_cast<unsigned int>(loader_.GetSlicesCount())));
@@ -301,7 +301,7 @@
   };
 
 
-  class VolumeImageGeometry
+  class OLD_VolumeImageGeometry
   {
   private:
     unsigned int         width_;
@@ -403,8 +403,8 @@
     }
 
   public:
-    VolumeImageGeometry(const OrthancVolumeImage& volume,
-                        VolumeProjection projection)
+    OLD_VolumeImageGeometry(const OrthancVolumeImage& volume,
+                            VolumeProjection projection)
     {
       if (volume.GetSlicesCount() == 0)
       {
@@ -521,9 +521,9 @@
 
 
     OrthancVolumeImage&                 volume_;
-    std::auto_ptr<VolumeImageGeometry>  axialGeometry_;
-    std::auto_ptr<VolumeImageGeometry>  coronalGeometry_;
-    std::auto_ptr<VolumeImageGeometry>  sagittalGeometry_;
+    std::auto_ptr<OLD_VolumeImageGeometry>  axialGeometry_;
+    std::auto_ptr<OLD_VolumeImageGeometry>  coronalGeometry_;
+    std::auto_ptr<OLD_VolumeImageGeometry>  sagittalGeometry_;
 
 
     bool IsGeometryReady() const
@@ -536,9 +536,9 @@
       assert(&message.GetOrigin() == &volume_);
 
       // These 3 values are only used to speed up the IVolumeSlicer
-      axialGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Axial));
-      coronalGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Coronal));
-      sagittalGeometry_.reset(new VolumeImageGeometry(volume_, VolumeProjection_Sagittal));
+      axialGeometry_.reset(new OLD_VolumeImageGeometry(volume_, VolumeProjection_Axial));
+      coronalGeometry_.reset(new OLD_VolumeImageGeometry(volume_, VolumeProjection_Coronal));
+      sagittalGeometry_.reset(new OLD_VolumeImageGeometry(volume_, VolumeProjection_Sagittal));
 
       BroadcastMessage(IVolumeSlicer::GeometryReadyMessage(*this));
     }
@@ -567,7 +567,7 @@
       BroadcastMessage(IVolumeSlicer::ContentChangedMessage(*this));
     }
 
-    const VolumeImageGeometry& GetProjectionGeometry(VolumeProjection projection)
+    const OLD_VolumeImageGeometry& GetProjectionGeometry(VolumeProjection projection)
     {
       if (!IsGeometryReady())
       {
@@ -676,7 +676,7 @@
       if (IsGeometryReady() &&
           DetectProjection(projection, viewportSlice))
       {
-        const VolumeImageGeometry& geometry = GetProjectionGeometry(projection);
+        const OLD_VolumeImageGeometry& geometry = GetProjectionGeometry(projection);
 
         size_t closest;
 
@@ -716,7 +716,7 @@
   private:
     SliceViewerWidget&                  widget_;
     VolumeProjection                    projection_;
-    std::auto_ptr<VolumeImageGeometry>  slices_;
+    std::auto_ptr<OLD_VolumeImageGeometry>  slices_;
     size_t                              slice_;
 
   protected:
@@ -727,7 +727,7 @@
         const OrthancVolumeImage& image =
           dynamic_cast<const OrthancVolumeImage&>(message.GetOrigin());
 
-        slices_.reset(new VolumeImageGeometry(image, projection_));
+        slices_.reset(new OLD_VolumeImageGeometry(image, projection_));
         SetSlice(slices_->GetSlicesCount() / 2);
 
         widget_.FitContent();
diff -r dbc1d8bfc68a -r 7719eb852dd5 Resources/CMake/OrthancStoneConfiguration.cmake
--- a/Resources/CMake/OrthancStoneConfiguration.cmake	Thu May 16 16:01:36 2019 +0200
+++ b/Resources/CMake/OrthancStoneConfiguration.cmake	Thu May 16 16:47:46 2019 +0200
@@ -391,6 +391,7 @@
   ${ORTHANC_STONE_ROOT}/Framework/Toolbox/SlicesSorter.cpp
   ${ORTHANC_STONE_ROOT}/Framework/Toolbox/UndoRedoStack.cpp
   ${ORTHANC_STONE_ROOT}/Framework/Toolbox/ViewportGeometry.cpp
+  ${ORTHANC_STONE_ROOT}/Framework/Toolbox/VolumeImageGeometry.cpp
   ${ORTHANC_STONE_ROOT}/Framework/Viewport/CairoContext.cpp
   ${ORTHANC_STONE_ROOT}/Framework/Viewport/CairoSurface.cpp
   ${ORTHANC_STONE_ROOT}/Framework/Viewport/IMouseTracker.h
diff -r dbc1d8bfc68a -r 7719eb852dd5 Samples/Sdl/Loader.cpp
--- a/Samples/Sdl/Loader.cpp	Thu May 16 16:01:36 2019 +0200
+++ b/Samples/Sdl/Loader.cpp	Thu May 16 16:47:46 2019 +0200
@@ -1292,8 +1292,9 @@
                                                      parameters.GetImageInformation().GetHeight(),
                                                      slices.GetSlicesCount(), false /* don't compute range */));      
 
-        image_->SetAxialGeometry(slices.GetSliceGeometry(0));
-        image_->SetVoxelDimensions(parameters.GetPixelSpacingX(), parameters.GetPixelSpacingY(), spacingZ);
+        image_->GetGeometry().SetAxialGeometry(slices.GetSliceGeometry(0));
+        image_->GetGeometry().SetVoxelDimensions(parameters.GetPixelSpacingX(),
+                                                 parameters.GetPixelSpacingY(), spacingZ);
       }
       
       image_->Clear();
@@ -1458,6 +1459,12 @@
 
         that_.volume_.SetGeometry(slices);
 
+        {
+          OrthancStone::LinearAlgebra::Print(that_.volume_.GetImage().GetGeometry().GetCoordinates(0, 0, 0));
+          OrthancStone::LinearAlgebra::Print(that_.volume_.GetImage().GetGeometry().GetCoordinates(1, 1, 1));
+          return;
+        }
+
         for (size_t i = 0; i < that_.volume_.GetSlicesCount(); i++)
         {
           const DicomInstanceParameters& slice = that_.volume_.GetSliceParameters(i);