Mercurial > hg > orthanc-stone

changeset 689:93a8949a1ef7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
VolumeImageGeometry::DetectSlice()
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 16 May 2019 18:33:57 +0200
parents ea1a5b963798
children 032a94cca5c4
files Framework/Toolbox/VolumeImageGeometry.cpp Framework/Toolbox/VolumeImageGeometry.h UnitTestsSources/UnitTestsMain.cpp
diffstat 3 files changed, 200 insertions(+), 28 deletions(-) [+]
line wrap: on
line diff
--- a/Framework/Toolbox/VolumeImageGeometry.cpp	Thu May 16 17:00:42 2019 +0200
+++ b/Framework/Toolbox/VolumeImageGeometry.cpp	Thu May 16 18:33:57 2019 +0200
@@ -39,7 +39,7 @@
     
     sagittalGeometry_ = CoordinateSystem3D(p,
                                            axialGeometry_.GetAxisY(),
-                                           -axialGeometry_.GetNormal());
+                                           axialGeometry_.GetNormal());
 
     Vector origin = (
       axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],
@@ -116,6 +116,25 @@
   }
 
 
+  const CoordinateSystem3D& VolumeImageGeometry::GetProjectionGeometry(VolumeProjection projection) const
+  {
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        return axialGeometry_;
+
+      case VolumeProjection_Coronal:
+        return coronalGeometry_;
+
+      case VolumeProjection_Sagittal:
+        return sagittalGeometry_;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+  }
+
+  
   Vector VolumeImageGeometry::GetVoxelDimensions(VolumeProjection projection) const
   {
     switch (projection)
@@ -135,26 +154,18 @@
   }
 
 
-  void VolumeImageGeometry::GetSliceSize(unsigned int& width,
-                                         unsigned int& height,
-                                         VolumeProjection projection)
+  unsigned int VolumeImageGeometry::GetProjectionWidth(VolumeProjection projection) const
   {
     switch (projection)
     {
       case VolumeProjection_Axial:
-        width = width_;
-        height = height_;
-        break;
+        return width_;
 
       case VolumeProjection_Coronal:
-        width = width_;
-        height = depth_;
-        break;
+        return width_;
 
       case VolumeProjection_Sagittal:
-        width = height_;
-        height = depth_;
-        break;
+        return height_;
 
       default:
         throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
@@ -162,6 +173,43 @@
   }
 
 
+  unsigned int VolumeImageGeometry::GetProjectionHeight(VolumeProjection projection) const
+  {
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        return height_;
+
+      case VolumeProjection_Coronal:
+        return depth_;
+
+      case VolumeProjection_Sagittal:
+        return depth_;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+  }
+
+
+  unsigned int VolumeImageGeometry::GetProjectionDepth(VolumeProjection projection) const
+  {
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        return depth_;
+
+      case VolumeProjection_Coronal:
+        return height_;
+
+      case VolumeProjection_Sagittal:
+        return width_;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }    
+  }
+
 
   Vector VolumeImageGeometry::GetCoordinates(float x,
                                              float y,
@@ -177,22 +225,19 @@
 
 
   bool VolumeImageGeometry::DetectProjection(VolumeProjection& projection,
-                                             const CoordinateSystem3D& plane) const
+                                             const Vector& planeNormal) const
   {
-    if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                    axialGeometry_.GetNormal()))
+    if (GeometryToolbox::IsParallel(planeNormal, axialGeometry_.GetNormal()))
     {
       projection = VolumeProjection_Axial;
       return true;
     }
-    else if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                         coronalGeometry_.GetNormal()))
+    else if (GeometryToolbox::IsParallel(planeNormal, coronalGeometry_.GetNormal()))
     {
       projection = VolumeProjection_Coronal;
       return true;
     }
-    else if (GeometryToolbox::IsParallel(plane.GetNormal(),
-                                         sagittalGeometry_.GetNormal()))
+    else if (GeometryToolbox::IsParallel(planeNormal, sagittalGeometry_.GetNormal()))
     {
       projection = VolumeProjection_Sagittal;
       return true;
@@ -202,4 +247,63 @@
       return false;
     }
   }
+
+  
+  bool VolumeImageGeometry::DetectSlice(VolumeProjection& projection,
+                                        unsigned int& slice,
+                                        const CoordinateSystem3D& plane) const
+  {
+    if (!DetectProjection(projection, plane.GetNormal()))
+    {
+      return false;
+    }
+
+    // Transforms the coordinates of the origin of the plane, into the
+    // coordinates of the axial geometry
+    const Vector& origin = plane.GetOrigin();
+    Vector p = LinearAlgebra::Product(
+      transformInverse_,
+      LinearAlgebra::CreateVector(origin[0], origin[1], origin[2], 1));
+
+    assert(LinearAlgebra::IsNear(p[3], 1));
+
+    double z;
+
+    switch (projection)
+    {
+      case VolumeProjection_Axial:
+        z = p[2];
+        break;
+
+      case VolumeProjection_Coronal:
+        z = p[1];
+        break;
+
+      case VolumeProjection_Sagittal:
+        z = p[0];
+        break;
+
+      default:
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+    }
+
+    const unsigned int projectionDepth = GetProjectionDepth(projection);
+    
+    z *= static_cast<double>(projectionDepth);
+    if (z < 0)
+    {
+      return false;
+    }
+        
+    unsigned int d = static_cast<unsigned int>(std::floor(z));
+    if (d >= projectionDepth)
+    {
+      return false;
+    }
+    else
+    {
+      slice = d;
+      return true;
+    }
+  }
 }
--- a/Framework/Toolbox/VolumeImageGeometry.h	Thu May 16 17:00:42 2019 +0200
+++ b/Framework/Toolbox/VolumeImageGeometry.h	Thu May 16 18:33:57 2019 +0200
@@ -74,6 +74,18 @@
       return sagittalGeometry_;
     }
 
+    const CoordinateSystem3D& GetProjectionGeometry(VolumeProjection projection) const;
+    
+    const Matrix& GetTransform() const
+    {
+      return transform_;
+    }
+
+    const Matrix& GetTransformInverse() const
+    {
+      return transformInverse_;
+    }
+
     void SetSize(unsigned int width,
                  unsigned int height,
                  unsigned int depth);
@@ -88,10 +100,12 @@
 
     Vector GetVoxelDimensions(VolumeProjection projection) const;
 
-    void GetSliceSize(unsigned int& width,
-                      unsigned int& height,
-                      VolumeProjection projection);
-    
+    unsigned int GetProjectionWidth(VolumeProjection projection) const;
+
+    unsigned int GetProjectionHeight(VolumeProjection projection) const;
+
+    unsigned int GetProjectionDepth(VolumeProjection projection) 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,
@@ -99,6 +113,10 @@
                           float z) const;
 
     bool DetectProjection(VolumeProjection& projection,
-                          const CoordinateSystem3D& plane) const;
+                          const Vector& planeNormal) const;
+
+    bool DetectSlice(VolumeProjection& projection,
+                     unsigned int& slice,
+                     const CoordinateSystem3D& plane) const;
   };
 }
--- a/UnitTestsSources/UnitTestsMain.cpp	Thu May 16 17:00:42 2019 +0200
+++ b/UnitTestsSources/UnitTestsMain.cpp	Thu May 16 18:33:57 2019 +0200
@@ -748,12 +748,62 @@
   ASSERT_FLOAT_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]);
 
   OrthancStone::VolumeProjection proj;
-  ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry()));
+  ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry().GetNormal()));
   ASSERT_EQ(OrthancStone::VolumeProjection_Axial, proj);
-  ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry()));
+  ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry().GetNormal()));
   ASSERT_EQ(OrthancStone::VolumeProjection_Coronal, proj);
-  ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry()));
+  ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry().GetNormal()));
   ASSERT_EQ(OrthancStone::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));
+
+  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Axial);
+  ASSERT_EQ(3u, p.size());
+  ASSERT_FLOAT_EQ(1, p[0]);
+  ASSERT_FLOAT_EQ(2, p[1]);
+  ASSERT_FLOAT_EQ(3, p[2]);
+  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Coronal);
+  ASSERT_EQ(3u, p.size());
+  ASSERT_FLOAT_EQ(1, p[0]);
+  ASSERT_FLOAT_EQ(3, p[1]);
+  ASSERT_FLOAT_EQ(2, p[2]);
+  p = g.GetVoxelDimensions(OrthancStone::VolumeProjection_Sagittal);
+  ASSERT_EQ(3u, p.size());
+  ASSERT_FLOAT_EQ(2, p[0]);
+  ASSERT_FLOAT_EQ(3, p[1]);
+  ASSERT_FLOAT_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);
+  
+  for (int p = 0; p < 3; p++)
+  {
+    OrthancStone::VolumeProjection projection = (OrthancStone::VolumeProjection) p;
+    const OrthancStone::CoordinateSystem3D& s = g.GetProjectionGeometry(projection);
+    
+    for (unsigned int i = 0; i < g.GetProjectionDepth(projection); i++)
+    {
+      OrthancStone::CoordinateSystem3D plane(
+        s.GetOrigin() + static_cast<double>(i) * s.GetNormal() * g.GetVoxelDimensions(projection)[2],
+        s.GetAxisX(),
+        s.GetAxisY());
+
+      unsigned int slice;
+      OrthancStone::VolumeProjection q;
+      ASSERT_TRUE(g.DetectSlice(q, slice, plane));
+      ASSERT_EQ(projection, q);
+      ASSERT_EQ(i, slice);
+    }
+  }
 }
 
 int main(int argc, char **argv)