# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1642596628 -3600
# Node ID cdf91ad891a5e141ade56b64229ecf380d65a204
# Parent  3716d72161d21cb2ed501d916385372ed6870c7e
estimated geometry of rt-struct

diff -r 3716d72161d2 -r cdf91ad891a5 OrthancStone/Sources/Toolbox/BucketAccumulator1D.h
--- a/OrthancStone/Sources/Toolbox/BucketAccumulator1D.h	Wed Jan 19 12:32:15 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/BucketAccumulator1D.h	Wed Jan 19 13:50:28 2022 +0100
@@ -38,6 +38,11 @@
     {
       size_t             count_;
       std::list<double>  values_;
+
+      Bucket() :
+        count_(0)
+      {
+      }
     };
     
     Internals::BucketMapper  mapper_;
diff -r 3716d72161d2 -r cdf91ad891a5 OrthancStone/Sources/Toolbox/BucketAccumulator2D.cpp
--- a/OrthancStone/Sources/Toolbox/BucketAccumulator2D.cpp	Wed Jan 19 12:32:15 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/BucketAccumulator2D.cpp	Wed Jan 19 13:50:28 2022 +0100
@@ -176,4 +176,29 @@
 
     y = LinearAlgebra::ComputeMedian(v);
   }
+
+
+  void BucketAccumulator2D::Print(FILE* fp) const
+  {
+    fprintf(fp, "         ");
+    
+    for (size_t x = 0; x < mapperX_.GetSize(); x++)
+    {
+      fprintf(fp, "%7.2f ", mapperX_.GetBucketCenter(x));
+    }
+
+    fprintf(fp, "\n");
+    
+    for (size_t y = 0; y < mapperY_.GetSize(); y++)
+    {
+      fprintf(fp, "%7.2f: ", mapperY_.GetBucketCenter(y));
+
+      for (size_t x = 0; x < mapperX_.GetSize(); x++)
+      {
+        fprintf(fp, "%7lu ", GetBucketContentSize(x, y));
+      }
+
+      fprintf(fp, "\n");
+    }
+  }
 }
diff -r 3716d72161d2 -r cdf91ad891a5 OrthancStone/Sources/Toolbox/BucketAccumulator2D.h
--- a/OrthancStone/Sources/Toolbox/BucketAccumulator2D.h	Wed Jan 19 12:32:15 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/BucketAccumulator2D.h	Wed Jan 19 13:50:28 2022 +0100
@@ -26,6 +26,7 @@
 #include "Internals/BucketMapper.h"
 
 #include <list>
+#include <stdio.h>
 #include <vector>
 
 
@@ -39,6 +40,11 @@
       size_t             count_;
       std::list<double>  valuesX_;
       std::list<double>  valuesY_;
+
+      Bucket() :
+        count_(0)
+      {
+      }
     };
     
     Internals::BucketMapper  mapperX_;
@@ -119,5 +125,7 @@
 
     void ComputeBestMedian(double& x,
                            double& y) const;
+
+    void Print(FILE* fp) const;
   };
 }
diff -r 3716d72161d2 -r cdf91ad891a5 OrthancStone/Sources/Toolbox/DicomStructureSet.cpp
--- a/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp	Wed Jan 19 12:32:15 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp	Wed Jan 19 13:50:28 2022 +0100
@@ -29,6 +29,8 @@
 
 #include "OrthancDatasets/DicomDatasetReader.h"
 
+#include "BucketAccumulator2D.h"
+
 #include <Logging.h>
 #include <OrthancException.h>
 #include <Toolbox.h>
@@ -44,6 +46,7 @@
 
 #include <limits>
 #include <stdio.h>
+#include <boost/math/constants/constants.hpp>
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/point_xy.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
@@ -58,7 +61,7 @@
 #endif
 
 
-typedef boost::geometry::model::d2::point_xy<double> BoostPoint;
+  typedef boost::geometry::model::d2::point_xy<double> BoostPoint;
 typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon;
 typedef boost::geometry::model::multi_polygon<BoostPolygon>  BoostMultiPolygon;
 
@@ -625,6 +628,9 @@
         structures_[i].polygons_.push_back(polygon);
       }
     }
+
+    EstimateGeometry();
+    
 #if STONE_TIME_BLOCKING_OPS
     boost::posix_time::ptime timerEnd = boost::posix_time::microsec_clock::universal_time();
     boost::posix_time::time_duration duration = timerEnd - timerStart;
@@ -1111,10 +1117,98 @@
   }
 
 
-  void DicomStructureSet::Test()
+  void DicomStructureSet::EstimateGeometry()
   {
-    printf("OK\n");
+    static const double PI = boost::math::constants::pi<double>();
+
+    BucketAccumulator2D accumulator(0, PI, 9,   /* for acos() */
+                                    -PI, PI, 9, /* for atan() */
+                                    true /* store values */);
+
+    unsigned int countPolygons = 0;
+    for (size_t i = 0; i < structures_.size(); i++)
+    {
+      const Polygons& polygons = structures_[i].polygons_;
+
+      for (Polygons::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
+      {
+        countPolygons++;
+        
+        const Points& points = it->GetPoints();
+        
+        if (points.size() >= 3)
+        {
+          // Compute the normal of the polygon using 3 successive points
+          Vector normal;
+          bool valid = false;
+
+          for (size_t i = 0; i + 2 < points.size(); i++)
+          {
+            const Vector& a = points[i];
+            const Vector& b = points[i + 1];
+            const Vector& c = points[i + 2];
+            LinearAlgebra::CrossProduct(normal, b - a, c - a);  // (*)
+            LinearAlgebra::NormalizeVector(normal);
+
+            if (LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal), 1.0))  // (**)
+            {
+              valid = true;
+              break;
+            }
+          }
 
+          if (valid)
+          {
+            // Check that all the points of the polygon lie in the plane defined by the normal
+            double d1 = GeometryToolbox::ProjectAlongNormal(points[0], normal);
+          
+            for (size_t i = 1; i < points.size(); i++)
+            {
+              double d2 = GeometryToolbox::ProjectAlongNormal(points[i], normal);
+            
+              if (!LinearAlgebra::IsNear(d1, d2))
+              {
+                valid = false;
+                break;
+              }
+            }
+          }
+
+          if (valid)
+          {
+            if (normal[2] < 0)
+            {
+              normal = -normal;
+            }
+
+            // The normal is a non-zero unit vector because of (*) and (**), so "r == 1"
+            // https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
+            
+            const double theta = acos(normal[2]);
+            const double phi = atan(normal[1]);
+            accumulator.AddValue(theta, phi);
+          }
+        }
+      }
+    }
+
+    size_t bestX, bestY;
+    accumulator.FindBestBucket(bestX, bestY);
+
+    if (accumulator.GetBucketContentSize(bestX, bestY) > 0)
+    {
+      double normalTheta, normalPhi;
+      accumulator.ComputeBestMedian(normalTheta, normalPhi);
+      
+      // Back to (x,y,z) normal, taking "r == 1"
+      // https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
+      double sinTheta = sin(normalTheta);
+      estimatedNormal_ = LinearAlgebra::CreateVector(sinTheta * cos(normalPhi), sinTheta * sin(normalPhi), cos(normalTheta));
+    }
+      
+    std::vector<double> polygonsProjection;
+    polygonsProjection.reserve(countPolygons);
+    
     for (size_t i = 0; i < structures_.size(); i++)
     {
       const Polygons& polygons = structures_[i].polygons_;
@@ -1122,28 +1216,31 @@
       for (Polygons::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
       {
         const Points& points = it->GetPoints();
-        
-        if (points.size() >= 3)
-        {
-          const Vector& a = points[0];
-          const Vector& b = points[1];
-          const Vector& c = points[2];
-          Vector n;
-          LinearAlgebra::CrossProduct(n, b - a, c - a);
-          LinearAlgebra::NormalizeVector(n);
-          if (n[2] < 0)
-            n = -n;
-          //LinearAlgebra::Print(n);
+        polygonsProjection.push_back(GeometryToolbox::ProjectAlongNormal(points[0], estimatedNormal_));
+      }
+    }
+
+    std::sort(polygonsProjection.begin(), polygonsProjection.end());
+
+    std::vector<double> deltas;
+    deltas.reserve(polygonsProjection.size());
 
-          // https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates#Vector_fields_2
-          double r = 1.0;
-          double theta = acos(n[2]);
-          double phi = atan(n[1]);
-          printf("%.02f %.02f %.02f => ", n[0], n[1], n[2]);
-          printf("%.02f %.02f =>", theta, phi);
-          printf("%.02f %.02f %.02f\n", r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta));
-        }
+    for (size_t i = 0; i + 1 < polygonsProjection.size(); i++)
+    {
+      if (!LinearAlgebra::IsNear(polygonsProjection[i], polygonsProjection[i + 1]))
+      {
+        assert(polygonsProjection[i + 1] > polygonsProjection[i]);
+        deltas.push_back(polygonsProjection[i + 1] - polygonsProjection[i]);
       }
     }
+
+    if (deltas.empty())
+    {
+      estimatedSliceThickness_ = 1;
+    }
+    else
+    {
+      estimatedSliceThickness_ = LinearAlgebra::ComputeMedian(deltas);
+    }
   }
 }
diff -r 3716d72161d2 -r cdf91ad891a5 OrthancStone/Sources/Toolbox/DicomStructureSet.h
--- a/OrthancStone/Sources/Toolbox/DicomStructureSet.h	Wed Jan 19 12:32:15 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/DicomStructureSet.h	Wed Jan 19 13:50:28 2022 +0100
@@ -152,6 +152,8 @@
 
     Structures        structures_;
     ReferencedSlices  referencedSlices_;
+    Vector            estimatedNormal_;
+    double            estimatedSliceThickness_;
 
     void Setup(const IDicomDataset& dataset);
     
@@ -168,6 +170,8 @@
       const Structure& structure,
       const CoordinateSystem3D& slice) const;
 
+    void EstimateGeometry();
+    
   public:
     explicit DicomStructureSet(const FullOrthancDataset& instance)
     {
@@ -242,6 +246,14 @@
                             size_t structureIndex,
                             const std::string& sopInstanceUid) const;
 
-    void Test();
+    const Vector& GetEstimatedNormal() const
+    {
+      return estimatedNormal_;
+    }
+
+    const double GetEstimatedSliceThickness() const
+    {
+      return estimatedSliceThickness_;
+    }
   };
 }
diff -r 3716d72161d2 -r cdf91ad891a5 RenderingPlugin/Sources/Plugin.cpp
--- a/RenderingPlugin/Sources/Plugin.cpp	Wed Jan 19 12:32:15 2022 +0100
+++ b/RenderingPlugin/Sources/Plugin.cpp	Wed Jan 19 13:50:28 2022 +0100
@@ -783,7 +783,8 @@
     case OrthancPluginChangeType_OrthancStarted:
     {
       DicomStructureCache::Accessor accessor(DicomStructureCache::GetSingleton(), "54460695-ba3885ee-ddf61ac0-f028e31d-a6e474d9");
-      accessor.GetRtStruct().Test();
+      OrthancStone::LinearAlgebra::Print(accessor.GetRtStruct().GetEstimatedNormal());
+      printf("Slice thickness: %f\n", accessor.GetRtStruct().GetEstimatedSliceThickness());
       break;
     }