From 0e7e75d3f55a470d45577fde970342d4342db106 Mon Sep 17 00:00:00 2001
From: ccolin <contact@ccolin.fr>
Date: Fri, 12 Nov 2021 21:23:50 +0100
Subject: [PATCH] better stats functions

---
 src/curvature.cpp | 21 ++++++--------
 src/curvature.h   | 73 -----------------------------------------------
 src/util.h        | 28 ++++++++++++++++++
 3 files changed, 37 insertions(+), 85 deletions(-)

diff --git a/src/curvature.cpp b/src/curvature.cpp
index 936ca81..18da97d 100644
--- a/src/curvature.cpp
+++ b/src/curvature.cpp
@@ -128,20 +128,17 @@ void Courbures::set_fixed_colors() {
 
 
 void Courbures::set_K_colors() {
-    MyStats<double> dist;
-    for (VertexHandle vh : _mesh.vertices()) {
-        dist.push_back(_mesh.property(vprop_K, vh));
-    }
-    double m = dist.mean();
-	double sigma = dist.stdev(m);
-    std::cout << "K min : " << dist.min()
-			  << " - max : " << dist.max() << std::endl ;
-    std::cout << "K mean : " << m << " - sigma : " << sigma << std::endl;
+	auto prop = _mesh.property(vprop_K);
+	auto prop_data = prop.data_vector();
+    double avg = average(prop_data.begin(), prop_data.end());
+	double stddev = standard_deviation
+		(prop_data.begin(), prop_data.end(), avg);
+    std::cout << "K mean : " << avg << " - sigma : " << stddev << std::endl;
 
     for (VertexHandle vh : _mesh.vertices()) {
-        double tmp  = _mesh.property(vprop_K, vh) - m;
-		tmp = clamp<double>(tmp, -sigma, +sigma);
-		tmp = (tmp + sigma) / (2 * sigma);
+        double tmp = _mesh.property(vprop_K, vh) - avg;
+		tmp = clamp<double>(tmp, -stddev, +stddev);
+		tmp = (tmp + stddev) / (2 * stddev);
 		_mesh.set_color(vh, MyMesh::Color(tmp, .23, .5));
     }
 }
diff --git a/src/curvature.h b/src/curvature.h
index db9b3d3..7f93aca 100644
--- a/src/curvature.h
+++ b/src/curvature.h
@@ -5,79 +5,6 @@
 #include "quad_patch.h"
 
 
-template <typename T>
-class MyStats {
-private:
-    std::vector<T> _distrib ;
-public:
-    MyStats () {} ;
-
-    void push_back (T data)
-		{
-			_distrib.push_back(data) ;
-		}
-
-    T min ()
-		{
-			T tmp = _distrib.at(0) ;
-			for (size_t i=1 ; i<_distrib.size(); i++)
-				{
-					if (_distrib.at(i) < tmp)
-						tmp = _distrib.at(i) ;
-				}
-			return tmp ;
-		}
-
-    T max ()
-		{
-			T tmp = _distrib.at(0) ;
-			for (size_t i=1 ; i<_distrib.size(); i++)
-				{
-					if (_distrib.at(i) > tmp)
-						tmp = _distrib.at(i) ;
-				}
-			return tmp ;
-		}
-
-    T mean ()
-		{
-			T acc(0) ;
-			std::cout << "acc : " << acc << std::endl;
-			if (_distrib.size() > 0)
-				{
-					for(size_t i=0; i<_distrib.size(); i++)
-						{
-							acc += _distrib.at(i) ;
-						}
-					return acc/(_distrib.size()) ;
-				}
-			else
-				return acc ;
-		}
-
-    T stdev ()
-		{
-			return stdev(mean());
-		}
-
-    T stdev (T m)
-		{
-			T acc(0), tmp ;
-			if (_distrib.size() > 0)
-				{
-					for(size_t i=0; i<_distrib.size(); i++)
-						{
-							tmp = _distrib.at(i) - m ;
-							acc += tmp * tmp ;
-						}
-					return sqrt(acc/(_distrib.size())) ;
-				}
-			else
-				return acc ;
-		}
-};
-
-
 class Courbures {
 private:
     MyMesh &_mesh ;
diff --git a/src/util.h b/src/util.h
index 068ebf7..10985a2 100644
--- a/src/util.h
+++ b/src/util.h
@@ -73,4 +73,32 @@ constexpr const T& clamp(const T& v, const T& lo, const T& hi) {
 }
 
 
+template <typename ForwardIt>
+constexpr typename ForwardIt::value_type average(const ForwardIt first,
+												 const ForwardIt last) {
+	typedef typename ForwardIt::value_type value_type;
+	if (first == last) return 0;
+	value_type average = 0;
+	for (ForwardIt it = first; it != last; ++it) {
+		average += *it;
+	}
+	return average / (last - first);
+}
+
+
+template <typename ForwardIt>
+constexpr  typename ForwardIt::value_type
+standard_deviation(const ForwardIt first,
+				   const ForwardIt last,
+				   typename ForwardIt::value_type average) {
+	typedef typename ForwardIt::value_type value_type;
+	if (first == last) return 0;
+	value_type accumulator = 0;
+	for (ForwardIt it = first; it != last; ++it) {
+		accumulator += (*it - average) * (*it - average);
+	}
+	return sqrt(accumulator / (last - first));
+}
+
+
 #endif