# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1558015296 -7200
# Node ID dbc1d8bfc68a51652a5047d9ac7cc85538456e1a
# Parent  9723fceccb9f58bd6c28031bd19550e47f85df51
reorganizing ImageBuffer3D

diff -r 9723fceccb9f -r dbc1d8bfc68a Framework/Toolbox/SlicesSorter.h
--- a/Framework/Toolbox/SlicesSorter.h	Thu May 16 14:26:11 2019 +0200
+++ b/Framework/Toolbox/SlicesSorter.h	Thu May 16 16:01:36 2019 +0200
@@ -84,7 +84,8 @@
     // slices. This is notably the case if all the slices are not
     // parallel to the reference normal that will be selected.
     bool Sort();
-    
+
+    // TODO - Remove this
     bool LookupClosestSlice(size_t& index,
                             double& distance,
                             const CoordinateSystem3D& slice) const;
diff -r 9723fceccb9f -r dbc1d8bfc68a Framework/Volumes/ImageBuffer3D.cpp
--- a/Framework/Volumes/ImageBuffer3D.cpp	Thu May 16 14:26:11 2019 +0200
+++ b/Framework/Volumes/ImageBuffer3D.cpp	Thu May 16 16:01:36 2019 +0200
@@ -21,6 +21,8 @@
 
 #include "ImageBuffer3D.h"
 
+#include "../Toolbox/GeometryToolbox.h"
+
 #include <Core/Images/ImageProcessing.h>
 #include <Core/Logging.h>
 #include <Core/OrthancException.h>
@@ -121,6 +123,8 @@
     LOG(INFO) << "Created a 3D image of size " << width << "x" << height
               << "x" << depth << " in " << Orthanc::EnumerationToString(format)
               << " (" << (GetEstimatedMemorySize() / (1024ll * 1024ll)) << "MB)";
+
+    UpdateGeometry();
   }
 
 
@@ -133,6 +137,7 @@
   void ImageBuffer3D::SetAxialGeometry(const CoordinateSystem3D& geometry)
   {
     axialGeometry_ = geometry;
+    UpdateGeometry();
   }
 
 
@@ -146,35 +151,30 @@
     {
       throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
     }
-
+    else
     {
       LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
+      UpdateGeometry();
     }
   }
 
 
   Vector ImageBuffer3D::GetVoxelDimensions(VolumeProjection projection) const
   {
-    Vector result;
     switch (projection)
     {
     case VolumeProjection_Axial:
-      result = voxelDimensions_;
-      break;
+      return voxelDimensions_;
 
     case VolumeProjection_Coronal:
-      LinearAlgebra::AssignVector(result, voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
-      break;
+      return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]);
 
     case VolumeProjection_Sagittal:
-      LinearAlgebra::AssignVector(result, voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
-      break;
+      return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]);
 
     default:
       throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
     }
-
-    return result;
   }
 
 
@@ -459,20 +459,78 @@
   }
 
 
+  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 ps = GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
+    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;
+    }
 
-    const CoordinateSystem3D& axial = GetAxialGeometry();
-    
-    Vector origin = (axial.MapSliceToWorldCoordinates(-0.5 * ps[0], -0.5 * ps[1]) -
-        0.5 * ps[2] * axial.GetNormal());
+    {
+      if (GeometryToolbox::IsParallel(plane.GetNormal(),
+                                      coronalGeometry_.GetNormal()))
+      {
+        projection = VolumeProjection_Coronal;
+        return true;
+      }
+    }
 
-    return (origin +
-            axial.GetAxisX() * ps[0] * x * static_cast<double>(GetWidth()) +
-        axial.GetAxisY() * ps[1] * y * static_cast<double>(GetHeight()) +
-        axial.GetNormal() * ps[2] * z * static_cast<double>(GetDepth()));
+    {
+      if (GeometryToolbox::IsParallel(plane.GetNormal(),
+                                      sagittalGeometry_.GetNormal()))
+      {
+        projection = VolumeProjection_Sagittal;
+        return true;
+      }
+    }
+
+    return false;
   }
 }
diff -r 9723fceccb9f -r dbc1d8bfc68a Framework/Volumes/ImageBuffer3D.h
--- a/Framework/Volumes/ImageBuffer3D.h	Thu May 16 14:26:11 2019 +0200
+++ b/Framework/Volumes/ImageBuffer3D.h	Thu May 16 16:01:36 2019 +0200
@@ -35,6 +35,8 @@
   {
   private:
     CoordinateSystem3D     axialGeometry_;
+    CoordinateSystem3D     coronalGeometry_;
+    CoordinateSystem3D     sagittalGeometry_;
     Vector                 voxelDimensions_;
     Orthanc::Image         image_;
     Orthanc::PixelFormat   format_;
@@ -45,6 +47,10 @@
     bool                   hasRange_;
     float                  minValue_;
     float                  maxValue_;
+    Matrix                 transform_;
+    Matrix                 transformInverse_;
+
+    void UpdateGeometry();
 
     void ExtendImageRange(const Orthanc::ImageAccessor& slice);
 
@@ -86,6 +92,16 @@
       return axialGeometry_;
     }
 
+    const CoordinateSystem3D& GetCoronalGeometry() const
+    {
+      return coronalGeometry_;
+    }
+
+    const CoordinateSystem3D& GetSagittalGeometry() const
+    {
+      return sagittalGeometry_;
+    }
+
     void SetVoxelDimensions(double x,
                             double y,
                             double z);
@@ -121,6 +137,7 @@
       return format_;
     }
 
+    // TODO - Remove
     ParallelSlices* GetGeometry(VolumeProjection projection) const;
     
     uint64_t GetEstimatedMemorySize() const;
@@ -166,6 +183,9 @@
                           float y,
                           float z) const;
 
+    bool DetectProjection(VolumeProjection& projection,
+                          const CoordinateSystem3D& plane) const;
+
 
     class SliceReader : public boost::noncopyable
     {
diff -r 9723fceccb9f -r dbc1d8bfc68a Samples/Sdl/Loader.cpp
--- a/Samples/Sdl/Loader.cpp	Thu May 16 14:26:11 2019 +0200
+++ b/Samples/Sdl/Loader.cpp	Thu May 16 16:01:36 2019 +0200
@@ -25,8 +25,9 @@
 #include "../../Framework/Messages/MessageBroker.h"
 #include "../../Framework/StoneInitialization.h"
 #include "../../Framework/Toolbox/GeometryToolbox.h"
+#include "../../Framework/Toolbox/SlicesSorter.h"
 #include "../../Framework/Volumes/ImageBuffer3D.h"
-#include "../../Framework/Toolbox/SlicesSorter.h"
+#include "../../Framework/Scene2D/ScenePoint2D.h"
 
 // From Orthanc framework
 #include <Core/Compression/GzipCompressor.h>
@@ -44,8 +45,8 @@
 #include <Core/Logging.h>
 #include <Core/MultiThreading/SharedMessageQueue.h>
 #include <Core/OrthancException.h>
+#include <Core/SystemToolbox.h>
 #include <Core/Toolbox.h>
-#include <Core/SystemToolbox.h>
 
 #include <json/reader.h>
 #include <json/value.h>
@@ -101,6 +102,17 @@
 
 
 
+  class IVolumeSlicer : public boost::noncopyable
+  {
+  public:
+    virtual ~IVolumeSlicer()
+    {
+    }
+
+    virtual void SetViewportPlane(const OrthancStone::CoordinateSystem3D& plane) = 0;
+  };
+
+
 
   class OracleCommandWithPayload : public IOracleCommand
   {
@@ -1444,11 +1456,11 @@
           slices.AddSlice(geometry, instance.release());
         }
 
-        that_.image_.SetGeometry(slices);
+        that_.volume_.SetGeometry(slices);
 
-        for (size_t i = 0; i < that_.image_.GetSlicesCount(); i++)
+        for (size_t i = 0; i < that_.volume_.GetSlicesCount(); i++)
         {
-          const DicomInstanceParameters& slice = that_.image_.GetSliceParameters(i);
+          const DicomInstanceParameters& slice = that_.volume_.GetSliceParameters(i);
           
           const std::string& instance = slice.GetOrthancInstanceIdentifier();
           if (instance.empty())
@@ -1490,7 +1502,7 @@
 #endif
 
           command->SetExpectedPixelFormat(slice.GetExpectedPixelFormat());
-          command->SetPayload(new LoadSliceImage(that_.image_, i));
+          command->SetPayload(new LoadSliceImage(that_.volume_, i));
 
           that_.oracle_.Schedule(that_, command.release());
         }
@@ -1521,7 +1533,7 @@
 
     IOracle&          oracle_;
     bool              active_;
-    DicomVolumeImage  image_;
+    DicomVolumeImage  volume_;
 
   public:
     VolumeSeriesOrthancLoader(IOracle& oracle,