# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1642540135 -3600
# Node ID 6ce81914f7e47458737d5cf55be9823d6e9d7f67
# Parent  fe4befc9c2b0f9434791c0a0a0985be40e2dcdb7
added classes BucketAccumulator1D/2D

diff -r fe4befc9c2b0 -r 6ce81914f7e4 OrthancStone/Sources/Toolbox/DicomStructureSet.cpp
--- a/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp	Tue Jan 18 17:52:43 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/DicomStructureSet.cpp	Tue Jan 18 22:08:55 2022 +0100
@@ -1109,4 +1109,41 @@
       }
     }
   }
+
+
+  void DicomStructureSet::Test()
+  {
+    printf("OK\n");
+
+    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)
+      {
+        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);
+
+          // 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));
+        }
+      }
+    }
+  }
 }
diff -r fe4befc9c2b0 -r 6ce81914f7e4 OrthancStone/Sources/Toolbox/DicomStructureSet.h
--- a/OrthancStone/Sources/Toolbox/DicomStructureSet.h	Tue Jan 18 17:52:43 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/DicomStructureSet.h	Tue Jan 18 22:08:55 2022 +0100
@@ -241,5 +241,7 @@
     void GetStructurePoints(std::list< std::vector<Vector> >& target,
                             size_t structureIndex,
                             const std::string& sopInstanceUid) const;
+
+    void Test();
   };
 }
diff -r fe4befc9c2b0 -r 6ce81914f7e4 OrthancStone/Sources/Toolbox/LinearAlgebra.cpp
--- a/OrthancStone/Sources/Toolbox/LinearAlgebra.cpp	Tue Jan 18 17:52:43 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/LinearAlgebra.cpp	Tue Jan 18 22:08:55 2022 +0100
@@ -717,6 +717,49 @@
 
       return m;
     }
+
+
+    template <typename T>
+    static T ComputeMedianInternal(std::vector<T>& v)
+    {
+      if (v.size() == 0)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange, "Empty vector");
+      }
+      
+      const size_t middle = v.size() / 2;
+      nth_element(v.begin(), v.begin() + middle, v.end());
+
+      T median = v[middle];
+      
+      if (v.size() % 2 == 1)
+      {
+        return median;
+      }
+      else
+      {
+        /**
+         * Side-effect of "nth_element()": "All of the elements before
+         * this new nth element are less than or equal to the elements
+         * after the new nth element."
+         * https://en.cppreference.com/w/cpp/algorithm/nth_element
+         **/
+        
+        typename std::vector<T>::const_iterator m = std::max_element(v.begin(), v.begin() + middle);
+
+        return (*m + median) / static_cast<T>(2);
+      }
+    }
+
+    double ComputeMedian(std::vector<double>& v)
+    {
+      return ComputeMedianInternal<double>(v);
+    }
+
+    float ComputeMedian(std::vector<float>& v)
+    {
+      return ComputeMedianInternal<float>(v);
+    }
   }
 
   std::ostream& operator<<(std::ostream& s, const Vector& vec)
diff -r fe4befc9c2b0 -r 6ce81914f7e4 OrthancStone/Sources/Toolbox/LinearAlgebra.h
--- a/OrthancStone/Sources/Toolbox/LinearAlgebra.h	Tue Jan 18 17:52:43 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/LinearAlgebra.h	Tue Jan 18 22:08:55 2022 +0100
@@ -297,5 +297,9 @@
     bool IsShearMatrix(const Matrix& shear);  
 
     Matrix InvertShearMatrix(const Matrix& shear);
+
+    double ComputeMedian(std::vector<double>& v);
+
+    float ComputeMedian(std::vector<float>& v);
   };
 }
diff -r fe4befc9c2b0 -r 6ce81914f7e4 OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp
--- a/OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp	Tue Jan 18 17:52:43 2022 +0100
+++ b/OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp	Tue Jan 18 22:08:55 2022 +0100
@@ -25,6 +25,7 @@
 
 #include "../Sources/Toolbox/Internals/OrientedIntegerLine2D.h"
 #include "../Sources/Toolbox/Internals/RectanglesIntegerProjection.h"
+#include "../Sources/Toolbox/LinearAlgebra.h"
 #include "../Sources/Toolbox/SegmentTree.h"
 #include "../Sources/Toolbox/UnionOfRectangles.h"
 
@@ -973,3 +974,416 @@
     ASSERT_TRUE(contours.front()[12].IsEqual(OrthancStone::ScenePoint2D(2, 5)));
   }
 }
+
+
+TEST(LinearAlgebra, ComputeMedian)
+{
+  {
+    std::vector<double> v;
+    ASSERT_THROW(OrthancStone::LinearAlgebra::ComputeMedian(v), Orthanc::OrthancException);
+    
+    v.push_back(1);
+    v.push_back(3);
+    v.push_back(3);
+    v.push_back(6);
+    v.push_back(7);
+    v.push_back(8);
+    v.push_back(9);
+    ASSERT_DOUBLE_EQ(6.0, OrthancStone::LinearAlgebra::ComputeMedian(v));
+  }
+
+  {
+    std::vector<double> v;
+    v.push_back(1);
+    v.push_back(2);
+    v.push_back(3);
+    v.push_back(4);
+    v.push_back(5);
+    v.push_back(6);
+    v.push_back(8);
+    v.push_back(9);
+    ASSERT_DOUBLE_EQ(4.5, OrthancStone::LinearAlgebra::ComputeMedian(v));
+  }  
+}
+
+
+
+
+
+
+namespace OrthancStone
+{
+  namespace Internals
+  {
+    class BucketMapper : public boost::noncopyable
+    {
+    private:
+      double  minValue_;
+      double  maxValue_;
+      size_t  bucketsCount_;
+
+    public:
+      BucketMapper(double minValue,
+                   double maxValue,
+                   size_t bucketsCount) :
+        minValue_(minValue),
+        maxValue_(maxValue),
+        bucketsCount_(bucketsCount)
+      {
+        if (minValue >= maxValue ||
+            bucketsCount <= 0)
+        {
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+        }    
+      }                   
+
+      size_t GetSize() const
+      {
+        return bucketsCount_;
+      }
+
+      void CheckIndex(size_t i) const
+      {
+        if (i >= bucketsCount_)
+        {
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+        }
+      }
+
+      double GetBucketLow(size_t i) const
+      {
+        CheckIndex(i);
+        double alpha = static_cast<double>(i) / static_cast<double>(bucketsCount_);
+        return (1.0 - alpha) * minValue_ + alpha * maxValue_;
+      }
+
+      double GetBucketHigh(size_t i) const
+      {    
+        CheckIndex(i);
+        double alpha = static_cast<double>(i + 1) / static_cast<double>(bucketsCount_);
+        return (1.0 - alpha) * minValue_ + alpha * maxValue_;
+      }
+
+      double GetBucketCenter(size_t i) const
+      {    
+        CheckIndex(i);
+        double alpha = (static_cast<double>(i) + 0.5) / static_cast<double>(bucketsCount_);
+        return (1.0 - alpha) * minValue_ + alpha * maxValue_;
+      }
+    
+      size_t GetBucketIndex(double value) const
+      {
+        if (value < minValue_ ||
+            value > maxValue_)
+        {
+          throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
+        }
+        else
+        {
+          double tmp = (value - minValue_) / (maxValue_ - minValue_) * static_cast<double>(bucketsCount_);
+          assert(tmp >= 0 && tmp <= static_cast<double>(bucketsCount_));
+
+          size_t bucket = static_cast<unsigned int>(std::floor(tmp));
+          if (bucket == bucketsCount_)  // This is the case if "value == maxValue_"
+          {
+            return bucketsCount_ - 1;
+          }
+          else
+          {
+            return bucket;
+          }
+        }
+      }    
+    };
+  }
+
+  
+  class BucketAccumulator1D : public boost::noncopyable
+  {
+  private:
+    struct Bucket
+    {
+      size_t             count_;
+      std::list<double>  values_;
+    };
+    
+    Internals::BucketMapper  mapper_;
+    std::vector<Bucket>      buckets_;
+    bool                     storeValues_;
+
+  public:
+    BucketAccumulator1D(double minValue,
+                        double maxValue,
+                        size_t countBuckets,
+                        bool storeValues) :
+      mapper_(minValue, maxValue, countBuckets),
+      buckets_(countBuckets),
+      storeValues_(storeValues)
+    {
+    }
+
+    size_t GetSize() const
+    {
+      return mapper_.GetSize();
+    }
+
+    double GetBucketLow(size_t i) const
+    {
+      return mapper_.GetBucketLow(i);
+    }
+
+    double GetBucketHigh(size_t i) const
+    {    
+      return mapper_.GetBucketHigh(i);
+    }
+
+    double GetBucketCenter(size_t i) const
+    {    
+      return mapper_.GetBucketCenter(i);
+    }
+    
+    size_t GetBucketContentSize(size_t i) const
+    {
+      mapper_.CheckIndex(i);
+      return buckets_[i].count_;
+    }
+
+    void AddValue(double value)
+    {
+      Bucket& bucket = buckets_[mapper_.GetBucketIndex(value)];
+
+      bucket.count_++;
+
+      if (storeValues_)
+      {
+        bucket.values_.push_back(value);
+      }
+    }
+
+    size_t FindBestBucket() const
+    {
+      size_t best = 0;
+
+      for (size_t i = 0; i < buckets_.size(); i++)
+      {
+        if (buckets_[i].count_ > buckets_[best].count_)
+        {
+          best = i;
+        }
+      }
+    
+      return best;
+    }
+
+    double ComputeBestCenter() const
+    {
+      return GetBucketCenter(FindBestBucket());
+    }
+
+    double ComputeBestMedian() const
+    {
+      if (!storeValues_)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+      }
+      
+      const std::list<double>& values = buckets_[FindBestBucket()].values_;
+
+      std::vector<double> v;
+      v.reserve(values.size());
+      for (std::list<double>::const_iterator it = values.begin(); it != values.end(); ++it)
+      {
+        v.push_back(*it);
+      }
+
+      return LinearAlgebra::ComputeMedian(v);
+    }
+  };
+
+
+  class BucketAccumulator2D : public boost::noncopyable
+  {
+  private:
+    struct Bucket
+    {
+      size_t             count_;
+      std::list<double>  valuesX_;
+      std::list<double>  valuesY_;
+    };
+    
+    Internals::BucketMapper  mapperX_;
+    Internals::BucketMapper  mapperY_;
+    std::vector<Bucket>      buckets_;
+    bool                     storeValues_;
+
+    size_t FindBestInternal() const
+    {
+      size_t best = 0;
+
+      for (size_t i = 0; i < buckets_.size(); i++)
+      {
+        if (buckets_[i].count_ > buckets_[best].count_)
+        {
+          best = i;
+        }
+      }
+
+      return best;
+    }
+
+  public:
+    BucketAccumulator2D(double minValueX,
+                        double maxValueX,
+                        size_t countBucketsX,
+                        double minValueY,
+                        double maxValueY,
+                        size_t countBucketsY,
+                        bool storeValues) :
+      mapperX_(minValueX, maxValueX, countBucketsX),
+      mapperY_(minValueY, maxValueY, countBucketsY),
+      buckets_(countBucketsX * countBucketsY),
+      storeValues_(storeValues)
+    {
+    }
+
+    void AddValue(double valueX,
+                  double valueY)
+    {
+      size_t x = mapperX_.GetBucketIndex(valueX);
+      size_t y = mapperY_.GetBucketIndex(valueY);
+      
+      Bucket& bucket = buckets_[x + y * mapperX_.GetSize()];
+
+      bucket.count_++;
+
+      if (storeValues_)
+      {
+        bucket.valuesX_.push_back(valueX);
+        bucket.valuesY_.push_back(valueY);
+      }
+    }
+
+    void FindBestBucket(size_t& bucketX,
+                        size_t& bucketY) const
+    {
+      size_t best = FindBestInternal();
+      bucketX = best % mapperX_.GetSize();
+      bucketY = best / mapperX_.GetSize();
+    }
+
+    void ComputeBestCenter(double& x,
+                           double& y) const
+    {
+      size_t bucketX, bucketY;
+      FindBestBucket(bucketX, bucketY);
+      x = mapperX_.GetBucketCenter(bucketX);
+      y = mapperY_.GetBucketCenter(bucketY);
+    }
+
+    void ComputeBestMedian(double& x,
+                           double& y) const
+    {
+      if (!storeValues_)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+      }
+      
+      const std::list<double>& valuesX = buckets_[FindBestInternal()].valuesX_;
+      const std::list<double>& valuesY = buckets_[FindBestInternal()].valuesY_;
+
+      std::vector<double> v;
+      v.reserve(valuesX.size());
+      for (std::list<double>::const_iterator it = valuesX.begin(); it != valuesX.end(); ++it)
+      {
+        v.push_back(*it);
+      }
+
+      x = LinearAlgebra::ComputeMedian(v);
+
+      v.clear();
+      v.reserve(valuesY.size());
+      for (std::list<double>::const_iterator it = valuesY.begin(); it != valuesY.end(); ++it)
+      {
+        v.push_back(*it);
+      }
+
+      y = LinearAlgebra::ComputeMedian(v);
+    }
+  };
+}
+
+
+TEST(BucketAccumulator1D, Basic)
+{
+  for (int storeValues = 0; storeValues <= 1; storeValues++)
+  {
+    OrthancStone::BucketAccumulator1D b(-10, 30, 4, storeValues != 0);
+    ASSERT_EQ(4u, b.GetSize());
+  
+    ASSERT_DOUBLE_EQ(-10.0, b.GetBucketLow(0));
+    ASSERT_DOUBLE_EQ(0.0, b.GetBucketLow(1));
+    ASSERT_DOUBLE_EQ(10.0, b.GetBucketLow(2));
+    ASSERT_DOUBLE_EQ(20.0, b.GetBucketLow(3));
+  
+    ASSERT_DOUBLE_EQ(0.0, b.GetBucketHigh(0));
+    ASSERT_DOUBLE_EQ(10.0, b.GetBucketHigh(1));
+    ASSERT_DOUBLE_EQ(20.0, b.GetBucketHigh(2));
+    ASSERT_DOUBLE_EQ(30.0, b.GetBucketHigh(3));
+  
+    ASSERT_DOUBLE_EQ(-5.0, b.GetBucketCenter(0));
+    ASSERT_DOUBLE_EQ(5.0, b.GetBucketCenter(1));
+    ASSERT_DOUBLE_EQ(15.0, b.GetBucketCenter(2));
+    ASSERT_DOUBLE_EQ(25.0, b.GetBucketCenter(3));
+
+    ASSERT_EQ(0u, b.GetBucketContentSize(0));
+    ASSERT_EQ(0u, b.GetBucketContentSize(1));
+    ASSERT_EQ(0u, b.GetBucketContentSize(2));
+    ASSERT_EQ(0u, b.GetBucketContentSize(3));
+
+    ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException);  // No data point
+
+    b.AddValue(-10.0);
+    b.AddValue(0.0);
+    b.AddValue(9.9999);
+    b.AddValue(10.0);
+    b.AddValue(20.0);
+    b.AddValue(29.9999);
+    b.AddValue(30.0);
+    ASSERT_THROW(b.AddValue(30.00001), Orthanc::OrthancException);
+
+    ASSERT_EQ(3u, b.FindBestBucket());
+    ASSERT_EQ(25.0, b.ComputeBestCenter());
+  
+    ASSERT_EQ(1u, b.GetBucketContentSize(0));
+    ASSERT_EQ(2u, b.GetBucketContentSize(1));
+    ASSERT_EQ(1u, b.GetBucketContentSize(2));
+    ASSERT_EQ(3u, b.GetBucketContentSize(3));
+
+    if (storeValues == 0)
+    {
+      ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException);
+    }
+    else
+    {
+      ASSERT_EQ(29.9999, b.ComputeBestMedian());
+    }
+  }
+}
+
+TEST(BucketAccumulator2D, Basic)
+{
+  for (int storeValues = 0; storeValues <= 1; storeValues++)
+  {
+    OrthancStone::BucketAccumulator2D b(-10, 30, 4,
+                                        0, 3, 3, storeValues != 0);
+
+    if (storeValues == 0)
+    {
+      //ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException);
+    }
+    else
+    {
+      //ASSERT_EQ(29.9999, b.ComputeBestMedian());
+    }
+  }
+}
diff -r fe4befc9c2b0 -r 6ce81914f7e4 RenderingPlugin/Sources/Plugin.cpp
--- a/RenderingPlugin/Sources/Plugin.cpp	Tue Jan 18 17:52:43 2022 +0100
+++ b/RenderingPlugin/Sources/Plugin.cpp	Tue Jan 18 22:08:55 2022 +0100
@@ -71,9 +71,9 @@
   class Accessor : public boost::noncopyable
   {
   private:
-    boost::mutex::scoped_lock               lock_;
-    std::string                             instanceId_;
-    const OrthancStone::DicomStructureSet*  rtstruct_;
+    boost::mutex::scoped_lock         lock_;
+    std::string                       instanceId_;
+    OrthancStone::DicomStructureSet*  rtstruct_;
 
   public:
     Accessor(DicomStructureCache& that,
@@ -113,7 +113,7 @@
       return rtstruct_ != NULL;
     }
 
-    const OrthancStone::DicomStructureSet& GetRtStruct() const
+    OrthancStone::DicomStructureSet& GetRtStruct() const
     {
       if (IsValid())
       {
@@ -765,6 +765,7 @@
 }
 
 
+
 OrthancPluginErrorCode OnChangeCallback(OrthancPluginChangeType changeType,
                                         OrthancPluginResourceType resourceType,
                                         const char* resourceId)
@@ -779,6 +780,13 @@
       
       break;
 
+    case OrthancPluginChangeType_OrthancStarted:
+    {
+      DicomStructureCache::Accessor accessor(DicomStructureCache::GetSingleton(), "54460695-ba3885ee-ddf61ac0-f028e31d-a6e474d9");
+      accessor.GetRtStruct().Test();
+      break;
+    }
+
     default:
       break;
   }