diff --git a/modules/img/alg/src/stat_accumulator.hh b/modules/img/alg/src/stat_accumulator.hh index 6dfc10aa7854efebcc3fab1359fe732fdb18c67d..cb4530eac56c0c133d58581359b19540d25c685d 100644 --- a/modules/img/alg/src/stat_accumulator.hh +++ b/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_; };