templated data type for StatAccumulator

modules/img/alg/src/stat_accumulator.hh

@@ -28,16 +28,28 @@
 
 namespace ost { namespace img { namespace alg {
 
+namespace {
 
-template<unsigned int MAX_MOMENT=4>
+template<class D>
+struct MinMax{
+  static bool less_cmp_(const D& v1, const D& v2){return v1<v2;}
+};
+template<>
+struct MinMax<Complex>{
+  static bool less_cmp_(const Complex& v1, const Complex& v2){return abs(v1)<abs(v2);}
+};
+
+} //anon ns
+
+template<unsigned int MAX_MOMENT=4,typename DATATYPE=Real>
 class  StatAccumulator
 {
 public:
   StatAccumulator():
     sum_(0.0),
     sum2_(0.0),
-    max_(-std::numeric_limits<Real>::max()),
-    min_(std::numeric_limits<Real>::max()),
+    max_(-std::numeric_limits<DATATYPE>::max()),
+    min_(std::numeric_limits<DATATYPE>::max()),
     m_(),
     w_(0.0),
     n_(0)
@@ -48,7 +60,7 @@ public:
 
   }
 
-  StatAccumulator(Real val, Real w=1.0):
+  StatAccumulator(DATATYPE val, Real w=1.0):
     sum_(val),
     sum2_(val*val),
     max_(val),
@@ -65,19 +77,19 @@ public:
     }
   }
 
-  void operator()(Real val, Real w=1.0)
+  void operator()(DATATYPE val, Real w=1.0)
   {
     *this+=StatAccumulator(val,w);
   }
 
-  StatAccumulator<MAX_MOMENT> operator+(const StatAccumulator<MAX_MOMENT>& acc2) const
+  StatAccumulator<MAX_MOMENT,DATATYPE> operator+(const StatAccumulator<MAX_MOMENT,DATATYPE>& acc2) const
   {
-    StatAccumulator<MAX_MOMENT> acc(acc2);
+    StatAccumulator<MAX_MOMENT,DATATYPE> acc(acc2);
     acc+=*this;
     return acc;
   }
 
-  StatAccumulator<MAX_MOMENT>& operator+=(const StatAccumulator<MAX_MOMENT>& acc)
+  StatAccumulator<MAX_MOMENT,DATATYPE>& operator+=(const StatAccumulator<MAX_MOMENT,DATATYPE>& acc)
   {
     if(0.0>=w_){
       *this=acc;
@@ -91,25 +103,25 @@ public:
 
     sum_+=acc.sum_;
     sum2_+=acc.sum2_;
-    max_=std::max<Real>(max_,acc.max_);
-    min_=std::min<Real>(min_,acc.min_);
+    max_=std::max<DATATYPE>(max_,acc.max_,MinMax<DATATYPE>::less_cmp_);
+    min_=std::min<DATATYPE>(min_,acc.min_,MinMax<DATATYPE>::less_cmp_);
     if(MAX_MOMENT>0){
-      Real delta=acc.m_[0]-m_[0];
-      Real delta_w=delta/wn;
+      DATATYPE delta=acc.m_[0]-m_[0];
+      DATATYPE delta_w=delta/wn;
       if(MAX_MOMENT>2){
-        Real w1w2_delta_w=w_*acc.w_*delta_w;
-        Real w1w2_delta_wp=w1w2_delta_w*w1w2_delta_w;
+        DATATYPE w1w2_delta_w=w_*acc.w_*delta_w;
+        DATATYPE w1w2_delta_wp=w1w2_delta_w*w1w2_delta_w;
         Real iw2=1.0/acc.w_;
         Real iw2pm1=iw2;
         Real miw1=-1.0/w_;
         Real miw1pm1=miw1;
-        Real mn[MAX_MOMENT]; // only MAX_MOMENT-2 values needed but overdimensioned to allow compilation for MAX_MOMENT==0 (even though it gets kicked out in the end by the dead code elimination)
+        DATATYPE mn[MAX_MOMENT]; // only MAX_MOMENT-2 values needed but overdimensioned to allow compilation for MAX_MOMENT==0 (even though it gets kicked out in the end by the dead code elimination)
         for(unsigned int p=3;p<=MAX_MOMENT;++p){
           w1w2_delta_wp*=w1w2_delta_w;
           iw2pm1*=iw2;
           miw1pm1*=miw1;
-          Real delta_wk=1.0;
-          Real s=0.0;
+          DATATYPE delta_wk=1.0;
+          DATATYPE s=0.0;
           Real mw2k=1.0;
           Real w1k=1.0;
           for(unsigned int k=1;k<=p-2;++k){
@@ -138,12 +150,12 @@ public:
     return n_;
   }
 
-  Real GetMaximum() const
+  DATATYPE GetMaximum() const
   {
     return max_;
   }
 
-  Real GetMinimum() const
+  DATATYPE GetMinimum() const
   {
     return min_;
   }
@@ -153,7 +165,7 @@ public:
     return w_;
   }
 
-  Real GetMoment(unsigned int i) const
+  DATATYPE GetMoment(unsigned int i) const
   {
     if(1>i){
       throw Error("Invalid moment.");
@@ -164,7 +176,7 @@ public:
     return m_[i-1];
   }
 
-  Real GetMean() const
+  DATATYPE GetMean() const
   {
     if(MAX_MOMENT<1){
       throw Error("Mean was not calculated.");
@@ -172,12 +184,12 @@ public:
     return m_[0];
   }
 
-  Real GetSum() const
+  DATATYPE GetSum() const
   {
     return sum_;
   }
 
-  Real GetVariance() const
+  DATATYPE GetVariance() const
   {
     if(MAX_MOMENT<2){
       throw Error("Variance was not calculated.");
@@ -188,12 +200,12 @@ public:
     return m_[1]/(w_);
   }
 
-  Real GetStandardDeviation() const
+  DATATYPE GetStandardDeviation() const
   {
     return sqrt(GetVariance());
   }
 
-  Real GetRootMeanSquare() const
+  DATATYPE GetRootMeanSquare() const
   {
     if(n_==0.0){
       return 0.0;
@@ -201,7 +213,7 @@ public:
     return sqrt(sum2_/(w_));
   }
 
-  Real GetSkewness() const
+  DATATYPE GetSkewness() const
   {
     if(MAX_MOMENT<3){
       throw Error("Skewness was not calculated.");
@@ -213,7 +225,7 @@ public:
     }
   }
 
-  Real GetKurtosis() const
+  DATATYPE GetKurtosis() const
   {
     if(MAX_MOMENT<4){
       throw Error("Kurtosis was not calculated.");
@@ -226,11 +238,11 @@ public:
   }
 
 private:
-  Real sum_;
-  Real sum2_;
-  Real max_;
-  Real min_;
-  Real m_[MAX_MOMENT];
+  DATATYPE sum_;
+  DATATYPE sum2_;
+  DATATYPE max_;
+  DATATYPE min_;
+  DATATYPE m_[MAX_MOMENT];
   Real w_;
   unsigned int n_;
 };