# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1667223766 -3600
# Node ID 187a261d7ae28b5f730aa917d372be0a4c644f0f
# Parent  20fa913272b78a688b66b776e8d25c842ae95eb3
computation of mean and stddev in rectangle probes

diff -r 20fa913272b7 -r 187a261d7ae2 Applications/Resources/RunCppCheck.sh
--- a/Applications/Resources/RunCppCheck.sh	Mon Oct 31 12:32:47 2022 +0100
+++ b/Applications/Resources/RunCppCheck.sh	Mon Oct 31 14:42:46 2022 +0100
@@ -30,8 +30,8 @@
 
 cat <<EOF > /tmp/cppcheck-suppressions.txt
 stlFindInsert:../../Applications/Samples/WebAssembly/SingleFrameViewer/SingleFrameViewerApplication.h
-stlFindInsert:../../Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp:508
-stlFindInsert:../../Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp:1151
+stlFindInsert:../../Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp:516
+stlFindInsert:../../Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp:1159
 unpreciseMathCall:../../OrthancStone/Sources/Scene2D/Internals/CairoFloatTextureRenderer.cpp
 unpreciseMathCall:../../OrthancStone/Sources/Scene2D/LookupTableTextureSceneLayer.cpp
 unreadVariable:../../OrthancStone/Sources/Viewport/SdlViewport.cpp:143
diff -r 20fa913272b7 -r 187a261d7ae2 OrthancStone/Sources/Scene2D/AnnotationsSceneLayer.cpp
--- a/OrthancStone/Sources/Scene2D/AnnotationsSceneLayer.cpp	Mon Oct 31 12:32:47 2022 +0100
+++ b/OrthancStone/Sources/Scene2D/AnnotationsSceneLayer.cpp	Mon Oct 31 14:42:46 2022 +0100
@@ -1524,60 +1524,91 @@
     Segment&  segment4_;
     Text&     label_;
 
-    void UpdateLabel()
+  protected:
+    virtual void UpdateProbeForLayer(const ISceneLayer& layer) ORTHANC_OVERRIDE
     {
-      TextSceneLayer content;
-
-      const double x1 = handle1_.GetCenter().GetX();
-      const double y1 = handle1_.GetCenter().GetY();
-      const double x2 = handle2_.GetCenter().GetX();
-      const double y2 = handle2_.GetCenter().GetY();
+      double x1 = handle1_.GetCenter().GetX();
+      double y1 = handle1_.GetCenter().GetY();
+      double x2 = handle2_.GetCenter().GetX();
+      double y2 = handle2_.GetCenter().GetY();
         
       // Put the label to the right of the right-most handle
       //const double y = std::min(y1, y2);
       const double y = (y1 + y2) / 2.0;
       if (x1 < x2)
       {
-        content.SetPosition(x2, y);
+        label_.SetPosition(x2, y);
       }
       else
       {
-        content.SetPosition(x1, y);
+        label_.SetPosition(x1, y);
       }
 
-      content.SetAnchor(BitmapAnchor_CenterLeft);
-      content.SetBorder(10);
-
-      const double area = std::abs(x1 - x2) * std::abs(y1 - y2);
-        
+      std::string text;
+      
       char buf[32];
 
-      switch (GetUnits())
+      if (GetUnits() == Units_Millimeters)
       {
-        case Units_Millimeters:
-          sprintf(buf, "%0.2f cm%c%c",
-                  area / 100.0,
-                  0xc2, 0xb2 /* two bytes corresponding to two power in UTF-8 */);
-          break;
+        const double area = std::abs(x1 - x2) * std::abs(y1 - y2);
 
-        case Units_Pixels:
-          // Don't report area (pixel-times-pixel is a strange unit)
-          sprintf(buf, "hello");
-          break;
-          
-        default:
-          throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
+        sprintf(buf, "%0.2f cm%c%c",
+                area / 100.0,
+                0xc2, 0xb2 /* two bytes corresponding to two power in UTF-8 */);
+        text = buf;
       }
 
-      content.SetText(buf);
+      if (layer.GetType() == ISceneLayer::Type_FloatTexture)
+      {
+        const TextureBaseSceneLayer& texture = dynamic_cast<const TextureBaseSceneLayer&>(layer);
+        const AffineTransform2D sceneToTexture = AffineTransform2D::Invert(texture.GetTransform());
+
+        const Orthanc::ImageAccessor& image = texture.GetTexture();
+        assert(image.GetFormat() == Orthanc::PixelFormat_Float32);
+
+        sceneToTexture.Apply(x1, y1);
+        sceneToTexture.Apply(x2, y2);
+        int ix1 = static_cast<int>(std::floor(x1));
+        int iy1 = static_cast<int>(std::floor(y1));
+        int ix2 = static_cast<int>(std::floor(x2));
+        int iy2 = static_cast<int>(std::floor(y2));
+
+        if (ix1 > ix2)
+        {
+          std::swap(ix1, ix2);
+        }
+
+        if (iy1 > iy2)
+        {
+          std::swap(iy1, iy2);
+        }
 
-      label_.SetContent(content);
-    }
-    
-  protected:
-    virtual void UpdateProbeForLayer(const ISceneLayer& layer) ORTHANC_OVERRIDE
-    {
-      // TODO
+        LinearAlgebra::OnlineVarianceEstimator estimator;
+
+        for (int y = std::max(0, iy1); y <= std::min(static_cast<int>(image.GetHeight()) - 1, iy2); y++)
+        {
+          int x = std::max(0, ix1);
+          
+          const float* p = reinterpret_cast<const float*>(image.GetConstRow(y)) + x;
+
+          for (; x <= std::min(static_cast<int>(image.GetWidth()) - 1, ix2); x++, p++)
+          {
+            estimator.AddSample(*p);
+          }
+        }
+
+        if (estimator.GetCount() > 0)
+        {
+          if (!text.empty())
+          {
+            text += "\n";
+          }
+          sprintf(buf, "Mean: %0.1f\nStdDev: %0.1f", estimator.GetMean(), estimator.GetStandardDeviation());
+          text += buf;
+        }
+      }
+      
+      label_.SetText(text);
     }
     
   public:
@@ -1594,8 +1625,13 @@
       segment4_(AddTypedPrimitive<Segment>(new Segment(*this, p1.GetX(), p2.GetY(), p1.GetX(), p1.GetY()))),
       label_(AddTypedPrimitive<Text>(new Text(that, *this)))
     {
+      TextSceneLayer content;
+      content.SetAnchor(BitmapAnchor_CenterLeft);
+      content.SetBorder(10);
+      content.SetText("?");
+
+      label_.SetContent(content);
       label_.SetColor(COLOR_TEXT);
-      UpdateLabel();
     }
 
     virtual unsigned int GetHandlesCount() const ORTHANC_OVERRIDE
@@ -1686,11 +1722,7 @@
         throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
       }
       
-      UpdateLabel();
-    }
-
-    virtual void UpdateProbe(const Scene2D& scene) ORTHANC_OVERRIDE
-    {
+      TagProbeAsChanged();
     }
 
     virtual void Serialize(Json::Value& target) ORTHANC_OVERRIDE
diff -r 20fa913272b7 -r 187a261d7ae2 OrthancStone/Sources/Toolbox/LinearAlgebra.cpp
--- a/OrthancStone/Sources/Toolbox/LinearAlgebra.cpp	Mon Oct 31 12:32:47 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/LinearAlgebra.cpp	Mon Oct 31 14:42:46 2022 +0100
@@ -41,6 +41,59 @@
 {
   namespace LinearAlgebra
   {
+    void OnlineVarianceEstimator::AddSample(double value)
+    {
+      count_++;
+      sum_ += value;
+      sumOfSquares_ += value * value;
+    }
+
+
+    void OnlineVarianceEstimator::Clear()
+    {
+      count_ = 0;
+      sum_ = 0;
+      sumOfSquares_ = 0;
+    }
+      
+
+    double OnlineVarianceEstimator::GetMean() const
+    {
+      if (count_ == 0)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+      }
+      else
+      {
+        return sum_ / static_cast<double>(count_);
+      }
+    }
+      
+
+    double OnlineVarianceEstimator::GetVariance() const
+    {
+      if (count_ == 0)
+      {
+        throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
+      }        
+      else if (count_ == 1)
+      {
+        return 0;
+      }
+      else
+      {
+        const double n = static_cast<double>(count_);
+        return (sumOfSquares_ * n - sum_ * sum_) / (n * (n - 1.0));
+      }
+    }
+
+      
+    double OnlineVarianceEstimator::GetStandardDeviation() const
+    {
+      return sqrt(GetVariance());
+    }
+    
+
     void Print(const Vector& v)
     {
       for (size_t i = 0; i < v.size(); i++)
diff -r 20fa913272b7 -r 187a261d7ae2 OrthancStone/Sources/Toolbox/LinearAlgebra.h
--- a/OrthancStone/Sources/Toolbox/LinearAlgebra.h	Mon Oct 31 12:32:47 2022 +0100
+++ b/OrthancStone/Sources/Toolbox/LinearAlgebra.h	Mon Oct 31 14:42:46 2022 +0100
@@ -46,6 +46,35 @@
 
   namespace LinearAlgebra
   {
+    class OnlineVarianceEstimator
+    {
+    private:
+      unsigned int  count_;
+      double        sum_;
+      double        sumOfSquares_;
+
+    public:
+      OnlineVarianceEstimator()
+      {
+        Clear();
+      }
+
+      unsigned int GetCount() const
+      {
+        return count_;
+      }
+
+      void AddSample(double value);
+
+      void Clear();
+
+      double GetMean() const;  // Same as "mean()" in Matlab/Octave
+
+      double GetVariance() const;  // Same as "var()" in Matlab/Octave
+
+      double GetStandardDeviation() const;  // Same as "std()" in Matlab/Octave
+    };
+    
     void Print(const Vector& v);
 
     void Print(const Matrix& m);
@@ -301,5 +330,5 @@
     double ComputeMedian(std::vector<double>& v);
 
     float ComputeMedian(std::vector<float>& v);
-  };
+  }
 }
diff -r 20fa913272b7 -r 187a261d7ae2 OrthancStone/UnitTestsSources/GenericToolboxTests.cpp
--- a/OrthancStone/UnitTestsSources/GenericToolboxTests.cpp	Mon Oct 31 12:32:47 2022 +0100
+++ b/OrthancStone/UnitTestsSources/GenericToolboxTests.cpp	Mon Oct 31 14:42:46 2022 +0100
@@ -4499,3 +4499,28 @@
   ASSERT_DOUBLE_EQ(9093     , v[10]);
   ASSERT_DOUBLE_EQ(0        , v[11]);
 }
+
+
+TEST(LinearAlgebra, OnlineVarianceEstimator)
+{
+  OrthancStone::LinearAlgebra::OnlineVarianceEstimator e;
+  ASSERT_EQ(0u, e.GetCount());
+  ASSERT_THROW(e.GetMean(), Orthanc::OrthancException);
+  ASSERT_THROW(e.GetVariance(), Orthanc::OrthancException);
+  ASSERT_THROW(e.GetStandardDeviation(), Orthanc::OrthancException);
+
+  e.AddSample(42);
+  ASSERT_EQ(1u, e.GetCount());
+  ASSERT_DOUBLE_EQ(42.0, e.GetMean());
+  ASSERT_DOUBLE_EQ(0.0, e.GetVariance());
+  ASSERT_DOUBLE_EQ(0.0, e.GetStandardDeviation());
+
+  e.Clear();
+  e.AddSample(87.9);
+  e.AddSample(-82.4);
+  e.AddSample(17.3);
+  ASSERT_EQ(3u, e.GetCount());
+  ASSERT_DOUBLE_EQ(7.6, e.GetMean());
+  ASSERT_DOUBLE_EQ(7321.09, e.GetVariance());
+  ASSERT_DOUBLE_EQ(85.5633683301447, e.GetStandardDeviation());  
+}