--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/OrthancStone/Resources/Documentation/Coordinates.txt	Wed May 12 19:57:50 2021 +0200
@@ -0,0 +1,24 @@
+
+Some notes about the coordinate systems
+=======================================
+
+* The 3D coordinates are expressed in the PATIENT system, *not* the
+  "gantry" coordinates. As a consequence, the tags "Image Position
+  Patient" (0020,0032) and "Image Orientation Patient" (0020,0037)
+  describe the 3D geometry of a slice.
+
+* The tag "Patient Position" (0018,5100) could be used to convert from
+  patient coordinates to gantry coordinates, but this would be useful
+  to device manufacturers, whereas Stone primarly deals with users of
+  modalities.
+
+* The "Image Position Patient" gives the 3D coordinates of the CENTER
+  of the voxel of slice.
+
+* In 2D compositors, the origin of a texture corresponds to the CORNER
+  of the texture (*not* to the center of the first pixel). Roughly
+  speaking, the operation "floor()" must be applied to move from
+  floating coordinates to pixel coordinates.
+
+* The classes deriving from "IVolumeSlicer" must pay to attention to
+  convert from centers of 3D voxels to bounding boxes of 2D pixels.

--- a/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.cpp	Wed May 12 17:43:51 2021 +0200
+++ b/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.cpp	Wed May 12 19:57:50 2021 +0200
@@ -165,10 +165,17 @@
                                                        const Vector& pixelOffsetX,
                                                        const Vector& pixelOffsetY)
   {
+    /**
+     * Shift from the center of the voxel (DICOM convention for 3D
+     * slices) to the corner of the voxel, because 2D textures are
+     * expressed relatively to their borders. (*)
+     **/
+    Vector p = origin + cuttingPlane.GetOrigin() - 0.5 * pixelOffsetX - 0.5 * pixelOffsetY;
+    
     double x0, y0, x1, y1, x2, y2;
-    cuttingPlane.ProjectPoint(x0, y0, origin + cuttingPlane.GetOrigin());
-    cuttingPlane.ProjectPoint(x1, y1, origin + cuttingPlane.GetOrigin() + pixelOffsetX);
-    cuttingPlane.ProjectPoint(x2, y2, origin + cuttingPlane.GetOrigin() + pixelOffsetY);
+    cuttingPlane.ProjectPoint(x0, y0, p);
+    cuttingPlane.ProjectPoint(x1, y1, p + pixelOffsetX);
+    cuttingPlane.ProjectPoint(x2, y2, p + pixelOffsetY);
 
     /**
 
@@ -215,7 +222,7 @@
         AffineTransform2D::CreateOffset(originX_, originY_),
         AffineTransform2D::CreateRotation(angle_),
         AffineTransform2D::CreateScaling(pixelSpacingX_, pixelSpacingY_),
-        AffineTransform2D::CreateOffset(-0.5, -0.5),
+        AffineTransform2D::CreateOffset(-0.5, -0.5),  // (*)
         AffineTransform2D::CreateFlip(flipX_, flipY_, width, height));
     }
     else

--- a/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.h	Wed May 12 17:43:51 2021 +0200
+++ b/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.h	Wed May 12 19:57:50 2021 +0200
@@ -59,6 +59,12 @@
   public:
     TextureBaseSceneLayer();
 
+  private:
+    /**
+     * TODO - The methods below could be removed, as well as the
+     * corresponding members
+     **/
+    
     // Center of the top-left pixel
     void SetOrigin(double x,
                    double y);
@@ -69,8 +75,6 @@
     // In radians
     void SetAngle(double angle);
 
-    void SetLinearInterpolation(bool isLinearInterpolation);
-
     void SetFlipX(bool flip);
     
     void SetFlipY(bool flip);
@@ -100,16 +104,6 @@
       return angle_;
     }
 
-    bool IsLinearInterpolation() const
-    {
-      return isLinearInterpolation_;
-    }
-
-    bool HasTexture() const
-    {
-      return (texture_.get() != NULL);
-    }
-
     bool IsFlipX() const
     {
       return flipX_;
@@ -120,15 +114,32 @@
       return flipY_;
     }
 
+  public:
+    bool IsLinearInterpolation() const
+    {
+      return isLinearInterpolation_;
+    }
+
+    void SetLinearInterpolation(bool isLinearInterpolation);
+
+    bool HasTexture() const
+    {
+      return (texture_.get() != NULL);
+    }
+
     const Orthanc::ImageAccessor& GetTexture() const;
 
     void SetTransform(const AffineTransform2D& transform);
 
     void ClearTransform();
 
-    // Initialize a transform that maps a texture slice in 3D, to a
-    // cutting plane (the cutting plane should be parallel to the 3D
-    // slice). The "pixelOffsetX/Y" must take pixel spacing into account.
+    /**
+     * Initialize a transform that maps a texture slice in 3D, to a
+     * cutting plane (the cutting plane should be parallel to the 3D
+     * slice). The "pixelOffsetX/Y" must take pixel spacing into
+     * account. This method automatically converts from voxel centers
+     * (3D) to pixel corners (2D).
+     **/
     void SetCuttingPlaneTransform(const CoordinateSystem3D& cuttingPlane,
                                   const Vector& origin,        // coordinates of the center of the voxel
                                   const Vector& pixelOffsetX,  // 3D offset from (0,0) voxel to (1,0) voxel

--- a/OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.cpp	Wed May 12 17:43:51 2021 +0200
+++ b/OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.cpp	Wed May 12 19:57:50 2021 +0200
@@ -95,11 +95,6 @@
     
     const CoordinateSystem3D& system = volume_.GetGeometry().GetProjectionGeometry(projection_);
 
-    /**
-     * TODO => There was a shift of (0.5, 0.5) introduced by
-     * TextureBaseSceneLayer::GetTransform(). Is it an error?
-     **/
-    
     Vector pixelSpacing = volume_.GetGeometry().GetVoxelDimensions(projection_);
 
     texture->SetCuttingPlaneTransform(cuttingPlane, system.GetOrigin(),

--- a/OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.cpp	Wed May 12 17:43:51 2021 +0200
+++ b/OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.cpp	Wed May 12 19:57:50 2021 +0200
@@ -87,11 +87,14 @@
         const Vector p1 = cuttingPlane.MapSliceToWorldCoordinates(x1, y1);        
         const Vector p2 = cuttingPlane.MapSliceToWorldCoordinates(x1, y2);
 
+        // The "0.5" shift is to move from the corner of voxel to the center of the voxel
+        
         if (1)
         {
-          texture->SetCuttingPlaneTransform(cuttingPlane, p1,
-                                            s * cuttingPlane.GetAxisX(),
-                                            s * cuttingPlane.GetAxisY());
+          texture->SetCuttingPlaneTransform(
+            cuttingPlane, p1 + 0.5 * s * cuttingPlane.GetAxisX() + 0.5 * s * cuttingPlane.GetAxisY(),
+            s * cuttingPlane.GetAxisX(),
+            s * cuttingPlane.GetAxisY());
         }
         else
         {
@@ -100,9 +103,10 @@
            * possible for the X axis?
            **/
           
-          texture->SetCuttingPlaneTransform(cuttingPlane, p2,
-                                            s * cuttingPlane.GetAxisX(),
-                                            -s * cuttingPlane.GetAxisY());
+          texture->SetCuttingPlaneTransform(
+            cuttingPlane, p2 + 0.5 * s * cuttingPlane.GetAxisX() + 0.5 * s * cuttingPlane.GetAxisY(),
+            s * cuttingPlane.GetAxisX(),
+            -s * cuttingPlane.GetAxisY());
         }
 
 #else

--- a/UnitTestsSources/VolumeRenderingTests.cpp	Wed May 12 17:43:51 2021 +0200
+++ b/UnitTestsSources/VolumeRenderingTests.cpp	Wed May 12 19:57:50 2021 +0200
@@ -21,6 +21,7 @@
 
 #include "../OrthancStone/Sources/Scene2D/CairoCompositor.h"
 #include "../OrthancStone/Sources/Scene2D/CopyStyleConfigurator.h"
+#include "../OrthancStone/Sources/Scene2D/ColorTextureSceneLayer.h"
 #include "../OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.h"
 #include "../OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.h"
 
@@ -29,7 +30,7 @@
 
 #include <gtest/gtest.h>
 
-TEST(VolumeRendering, Basic)
+TEST(VolumeRendering, Axial)
 {
   Orthanc::DicomMap dicom;
   dicom.SetValue(Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, "study", false);
@@ -43,7 +44,7 @@
   OrthancStone::VolumeImageGeometry geometry;
   geometry.SetSizeInVoxels(3, 3, 1);
   geometry.SetAxialGeometry(axial);
-  
+
   boost::shared_ptr<OrthancStone::DicomVolumeImage> volume(new OrthancStone::DicomVolumeImage);
   volume->Initialize(geometry, Orthanc::PixelFormat_Grayscale8, false);
   volume->SetDicomParameters(OrthancStone::DicomInstanceParameters(dicom));
@@ -62,6 +63,11 @@
     }
   }
 
+  OrthancStone::Vector v = volume->GetGeometry().GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
+  ASSERT_FLOAT_EQ(1, v[0]);
+  ASSERT_FLOAT_EQ(1, v[1]);
+  ASSERT_FLOAT_EQ(1, v[2]);
+  
   OrthancStone::CoordinateSystem3D viewpoint;
 
   for (unsigned int mode = 0; mode < 2; mode++)
@@ -83,6 +89,15 @@
     OrthancStone::CopyStyleConfigurator configurator;
     std::unique_ptr<OrthancStone::ISceneLayer> layer(slice->CreateSceneLayer(&configurator, viewpoint));
 
+    ASSERT_EQ(OrthancStone::ISceneLayer::Type_FloatTexture, layer->GetType());
+
+    OrthancStone::Extent2D box;
+    layer->GetBoundingBox(box);
+    ASSERT_FLOAT_EQ(-1.0f, box.GetX1());
+    ASSERT_FLOAT_EQ(-1.0f, box.GetY1());
+    ASSERT_FLOAT_EQ(2.0f, box.GetX2());
+    ASSERT_FLOAT_EQ(2.0f, box.GetY2());
+    
     {
       const Orthanc::ImageAccessor& a = dynamic_cast<OrthancStone::TextureBaseSceneLayer&>(*layer).GetTexture();
       Orthanc::Image i(Orthanc::PixelFormat_Grayscale8, a.GetWidth(), a.GetHeight(), false);
@@ -144,4 +159,3 @@
     Orthanc::ImageTraits<Orthanc::PixelFormat_RGB24>::GetPixel(pixel, j, 4, 4);  ASSERT_EQ(200, pixel.red_);
   }
 }
-