# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1574251755 -3600
# Node ID 1644de437a7b16a9c962c9971ae88ac05c644a57
# Parent  ad4e21df4e40f14a601c93c1469932b197670c42
fixes related to swapped normal in sagittal geometry

diff -r ad4e21df4e40 -r 1644de437a7b Framework/Toolbox/DicomStructureSet.cpp
--- a/Framework/Toolbox/DicomStructureSet.cpp	Wed Nov 20 10:55:44 2019 +0100
+++ b/Framework/Toolbox/DicomStructureSet.cpp	Wed Nov 20 13:09:15 2019 +0100
@@ -373,8 +373,10 @@
     else if (GeometryToolbox::IsParallelOrOpposite
              (isOpposite, slice.GetNormal(), geometry_.GetAxisX()))
     {
-      // plane is constant X
+      // plane is constant X => Sagittal view (remember that in the
+      // sagittal projection, the normal must be swapped)
 
+      
       /*
       Please read the comments in the section above, by taking into account
       the fact that, in this case, the plane has a constant X, not Y (in 
@@ -427,10 +429,6 @@
         slice.ProjectPoint2(x1, y1, p1);
         slice.ProjectPoint2(x2, y2, p2);
 
-        // TODO WHY THIS???
-        y1 = -y1;
-        y2 = -y2;
-
         return true;
       }
     }
diff -r ad4e21df4e40 -r 1644de437a7b Framework/Toolbox/DicomStructureSet.h
--- a/Framework/Toolbox/DicomStructureSet.h	Wed Nov 20 10:55:44 2019 +0100
+++ b/Framework/Toolbox/DicomStructureSet.h	Wed Nov 20 13:09:15 2019 +0100
@@ -26,7 +26,7 @@
 #include "Extent2D.h"
 #include "../Scene2D/Color.h"
 
-#define USE_BOOST_UNION_FOR_POLYGONS 1
+//#define USE_BOOST_UNION_FOR_POLYGONS 1
 
 #include <Plugins/Samples/Common/FullOrthancDataset.h>
 
diff -r ad4e21df4e40 -r 1644de437a7b Framework/Toolbox/LinearAlgebra.cpp
--- a/Framework/Toolbox/LinearAlgebra.cpp	Wed Nov 20 10:55:44 2019 +0100
+++ b/Framework/Toolbox/LinearAlgebra.cpp	Wed Nov 20 13:09:15 2019 +0100
@@ -22,6 +22,7 @@
 #include "LinearAlgebra.h"
 
 #include "../StoneException.h"
+#include "GenericToolbox.h"
 
 #include <Core/Logging.h>
 #include <Core/OrthancException.h>
@@ -123,6 +124,12 @@
           return false;
         }
 
+#elif 1
+        if (!GenericToolbox::StringToDouble(target[i], items[i].c_str()))
+        {
+          return false;
+        }
+        
 #else
         /**
          * Fallback implementation using Boost (slower, but somewhat
diff -r ad4e21df4e40 -r 1644de437a7b Framework/Volumes/VolumeImageGeometry.cpp
--- a/Framework/Volumes/VolumeImageGeometry.cpp	Wed Nov 20 10:55:44 2019 +0100
+++ b/Framework/Volumes/VolumeImageGeometry.cpp	Wed Nov 20 13:09:15 2019 +0100
@@ -296,16 +296,18 @@
     {
       return false;
     }
-        
-    unsigned int d = static_cast<unsigned int>(std::floor(z));
-    if (d >= projectionDepth)
-    {
-      return false;
-    }
     else
     {
-      slice = d;
-      return true;
+      unsigned int d = static_cast<unsigned int>(std::floor(z));
+      if (d >= projectionDepth)
+      {
+        return false;
+      }
+      else
+      {
+        slice = d;
+        return true;
+      }
     }
   }
 
@@ -321,7 +323,18 @@
     Vector dim = GetVoxelDimensions(projection);
     CoordinateSystem3D plane = GetProjectionGeometry(projection);
 
-    plane.SetOrigin(plane.GetOrigin() + static_cast<double>(z) * plane.GetNormal() * dim[2]);
+    Vector normal = plane.GetNormal();
+    if (projection == VolumeProjection_Sagittal)
+    {
+      /**
+       * WARNING: In sagittal geometry, the normal points to REDUCING
+       * X-axis in the 3D world. This is necessary to keep the
+       * right-hand coordinate system. Hence the negation.
+       **/
+      normal = -normal;
+    }
+    
+    plane.SetOrigin(plane.GetOrigin() + static_cast<double>(z) * dim[2] * normal);
 
     return plane;
   }
diff -r ad4e21df4e40 -r 1644de437a7b UnitTestsSources/UnitTestsMain.cpp
--- a/UnitTestsSources/UnitTestsMain.cpp	Wed Nov 20 10:55:44 2019 +0100
+++ b/UnitTestsSources/UnitTestsMain.cpp	Wed Nov 20 13:09:15 2019 +0100
@@ -570,17 +570,22 @@
 static bool IsEqualVector(OrthancStone::Vector a,
                           OrthancStone::Vector b)
 {
-  if (a.size() == 3 &&
-      b.size() == 3)
+  if (a.size() != b.size())
   {
-    OrthancStone::LinearAlgebra::NormalizeVector(a);
-    OrthancStone::LinearAlgebra::NormalizeVector(b);
-    return OrthancStone::LinearAlgebra::IsCloseToZero(boost::numeric::ublas::norm_2(a - b));
+    return false;
   }
   else
   {
-    return false;
-  } 
+    for (size_t i = 0; i < a.size(); i++)
+    {
+      if (!OrthancStone::LinearAlgebra::IsNear(a[i], b[i], 0.0001))
+      {
+        return false;
+      }
+    }
+
+    return true;
+  }
 }
 
 
@@ -642,11 +647,13 @@
 
 TEST(VolumeImageGeometry, Basic)
 {
-  OrthancStone::VolumeImageGeometry g;
+  using namespace OrthancStone;
+  
+  VolumeImageGeometry g;
   g.SetSizeInVoxels(10, 20, 30);
   g.SetVoxelDimensions(1, 2, 3);
 
-  OrthancStone::Vector p = g.GetCoordinates(0, 0, 0);
+  Vector p = g.GetCoordinates(0, 0, 0);
   ASSERT_EQ(3u, p.size());
   ASSERT_DOUBLE_EQ(-1.0 / 2.0, p[0]);
   ASSERT_DOUBLE_EQ(-2.0 / 2.0, p[1]);
@@ -657,66 +664,116 @@
   ASSERT_DOUBLE_EQ(-2.0 / 2.0 + 20.0 * 2.0, p[1]);
   ASSERT_DOUBLE_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]);
 
-  OrthancStone::VolumeProjection proj;
+  VolumeProjection proj;
   ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry().GetNormal()));
-  ASSERT_EQ(OrthancStone::VolumeProjection_Axial, proj);
+  ASSERT_EQ(VolumeProjection_Axial, proj);
   ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry().GetNormal()));
-  ASSERT_EQ(OrthancStone::VolumeProjection_Coronal, proj);
+  ASSERT_EQ(VolumeProjection_Coronal, proj);
   ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry().GetNormal()));
-  ASSERT_EQ(OrthancStone::VolumeProjection_Sagittal, proj);
+  ASSERT_EQ(VolumeProjection_Sagittal, proj);
 
-  ASSERT_EQ(10u, g.GetProjectionWidth(OrthancStone::VolumeProjection_Axial));
-  ASSERT_EQ(20u, g.GetProjectionHeight(OrthancStone::VolumeProjection_Axial));
-  ASSERT_EQ(30u, g.GetProjectionDepth(OrthancStone::VolumeProjection_Axial));
-  ASSERT_EQ(10u, g.GetProjectionWidth(OrthancStone::VolumeProjection_Coronal));
-  ASSERT_EQ(30u, g.GetProjectionHeight(OrthancStone::VolumeProjection_Coronal));
-  ASSERT_EQ(20u, g.GetProjectionDepth(OrthancStone::VolumeProjection_Coronal));
-  ASSERT_EQ(20u, g.GetProjectionWidth(OrthancStone::VolumeProjection_Sagittal));
-  ASSERT_EQ(30u, g.GetProjectionHeight(OrthancStone::VolumeProjection_Sagittal));
-  ASSERT_EQ(10u, g.GetProjectionDepth(OrthancStone::VolumeProjection_Sagittal));
+  ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Axial));
+  ASSERT_EQ(20u, g.GetProjectionHeight(VolumeProjection_Axial));
+  ASSERT_EQ(30u, g.GetProjectionDepth(VolumeProjection_Axial));
+  ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Coronal));
+  ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Coronal));
+  ASSERT_EQ(20u, g.GetProjectionDepth(VolumeProjection_Coronal));
+  ASSERT_EQ(20u, g.GetProjectionWidth(VolumeProjection_Sagittal));
+  ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Sagittal));
+  ASSERT_EQ(10u, g.GetProjectionDepth(VolumeProjection_Sagittal));
 
-  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
+  p = g.GetVoxelDimensions(VolumeProjection_Axial);
   ASSERT_EQ(3u, p.size());
   ASSERT_DOUBLE_EQ(1, p[0]);
   ASSERT_DOUBLE_EQ(2, p[1]);
   ASSERT_DOUBLE_EQ(3, p[2]);
-  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Coronal);
+  p = g.GetVoxelDimensions(VolumeProjection_Coronal);
   ASSERT_EQ(3u, p.size());
   ASSERT_DOUBLE_EQ(1, p[0]);
   ASSERT_DOUBLE_EQ(3, p[1]);
   ASSERT_DOUBLE_EQ(2, p[2]);
-  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Sagittal);
+  p = g.GetVoxelDimensions(VolumeProjection_Sagittal);
   ASSERT_EQ(3u, p.size());
   ASSERT_DOUBLE_EQ(2, p[0]);
   ASSERT_DOUBLE_EQ(3, p[1]);
   ASSERT_DOUBLE_EQ(1, p[2]);
 
-  ASSERT_EQ(0, (int) OrthancStone::VolumeProjection_Axial);
-  ASSERT_EQ(1, (int) OrthancStone::VolumeProjection_Coronal);
-  ASSERT_EQ(2, (int) OrthancStone::VolumeProjection_Sagittal);
+  // Loop over all the voxels of the volume
+  for (unsigned int z = 0; z < g.GetDepth(); z++)
+  {
+    const float zz = (0.5f + static_cast<float>(z)) / static_cast<float>(g.GetDepth());  // Z-center of the voxel
+    
+    for (unsigned int y = 0; y < g.GetHeight(); y++)
+    {
+      const float yy = (0.5f + static_cast<float>(y)) / static_cast<float>(g.GetHeight());  // Y-center of the voxel
+
+      for (unsigned int x = 0; x < g.GetWidth(); x++)
+      {
+        const float xx = (0.5f + static_cast<float>(x)) / static_cast<float>(g.GetWidth());  // X-center of the voxel
+
+        const float sx = 1.0f;
+        const float sy = 2.0f;
+        const float sz = 3.0f;
+        
+        Vector p = g.GetCoordinates(xx, yy, zz);
+
+        Vector q = (g.GetAxialGeometry().MapSliceToWorldCoordinates(
+                      static_cast<double>(x) * sx,
+                      static_cast<double>(y) * sy) +
+                    z * sz * g.GetAxialGeometry().GetNormal());
+        ASSERT_TRUE(IsEqualVector(p, q));
+        
+        q = (g.GetCoronalGeometry().MapSliceToWorldCoordinates(
+               static_cast<double>(x) * sx,
+               static_cast<double>(g.GetDepth() - 1 - z) * sz) +
+             y * sy * g.GetCoronalGeometry().GetNormal());
+        ASSERT_TRUE(IsEqualVector(p, q));
+
+        /**
+         * WARNING: In sagittal geometry, the normal points to
+         * REDUCING X-axis in the 3D world. This is necessary to keep
+         * the right-hand coordinate system. Hence the "-".
+         **/
+        q = (g.GetSagittalGeometry().MapSliceToWorldCoordinates(
+               static_cast<double>(y) * sy,
+               static_cast<double>(g.GetDepth() - 1 - z) * sz) +
+             x * sx * (-g.GetSagittalGeometry().GetNormal()));
+        ASSERT_TRUE(IsEqualVector(p, q));
+      }
+    }
+  }
+
+  ASSERT_EQ(0, (int) VolumeProjection_Axial);
+  ASSERT_EQ(1, (int) VolumeProjection_Coronal);
+  ASSERT_EQ(2, (int) VolumeProjection_Sagittal);
   
   for (int p = 0; p < 3; p++)
   {
-    printf("=========== %d\n", p);
-    
-    OrthancStone::VolumeProjection projection = (OrthancStone::VolumeProjection) p;
-    const OrthancStone::CoordinateSystem3D& s = g.GetProjectionGeometry(projection);
+    VolumeProjection projection = (VolumeProjection) p;
+    const CoordinateSystem3D& s = g.GetProjectionGeometry(projection);
     
     ASSERT_THROW(g.GetProjectionSlice(projection, g.GetProjectionDepth(projection)), Orthanc::OrthancException);
 
     for (unsigned int i = 0; i < g.GetProjectionDepth(projection); i++)
     {
-      printf("------------- %d\n", i);
-    
-      OrthancStone::CoordinateSystem3D plane = g.GetProjectionSlice(projection, i);
+      CoordinateSystem3D plane = g.GetProjectionSlice(projection, i);
 
-      ASSERT_TRUE(IsEqualVector(plane.GetOrigin(), s.GetOrigin() + static_cast<double>(i) * 
-                                s.GetNormal() * g.GetVoxelDimensions(projection)[2]));
+      if (projection == VolumeProjection_Sagittal)
+      {
+        ASSERT_TRUE(IsEqualVector(plane.GetOrigin(), s.GetOrigin() + static_cast<double>(i) * 
+                                  (-s.GetNormal()) * g.GetVoxelDimensions(projection)[2]));
+      }
+      else
+      {
+        ASSERT_TRUE(IsEqualVector(plane.GetOrigin(), s.GetOrigin() + static_cast<double>(i) * 
+                                  s.GetNormal() * g.GetVoxelDimensions(projection)[2]));
+      }
+      
       ASSERT_TRUE(IsEqualVector(plane.GetAxisX(), s.GetAxisX()));
       ASSERT_TRUE(IsEqualVector(plane.GetAxisY(), s.GetAxisY()));
 
       unsigned int slice;
-      OrthancStone::VolumeProjection q;
+      VolumeProjection q;
       ASSERT_TRUE(g.DetectSlice(q, slice, plane));
       ASSERT_EQ(projection, q);
       ASSERT_EQ(i, slice);