diff --git a/modules/img/alg/src/CMakeLists.txt b/modules/img/alg/src/CMakeLists.txt
index c181156dbfe99eb4ac8b6a4532fa359b8b161d45..5dc212da9f112a9eed2dafc53eb187254a550a95 100644
--- a/modules/img/alg/src/CMakeLists.txt
+++ b/modules/img/alg/src/CMakeLists.txt
@@ -25,7 +25,6 @@ power_spectrum.cc
randomize.cc
smooth_mask_image.cc
stat.cc
-stat_accumulator.cc
stat_min_max.cc
alg_transform.cc
alg_mirror.cc
diff --git a/modules/img/alg/src/stat.cc b/modules/img/alg/src/stat.cc
index f433cd6d69a7810a2e34ca043ad3257720cb9acc..01053ab7e28bbccdfaf0a7c4839633246d675995 100644
--- a/modules/img/alg/src/stat.cc
+++ b/modules/img/alg/src/stat.cc
@@ -64,7 +64,7 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
ValIndex maxindex(-std::numeric_limits<Real>::max(),Point(0,0,0));
min_ = std::numeric_limits<Real>::max();
max_ = -std::numeric_limits<Real>::max();
- StatAccumulator acc;
+ StatAccumulator<> acc;
int wi=isi.GetSize()[0];
int he=isi.GetSize()[1];
@@ -124,7 +124,7 @@ void StatBase::VisitFunction(const Function& fnc)
ValIndex maxindex(-std::numeric_limits<Real>::max(),Point(0,0,0));
min_ = std::numeric_limits<Real>::max();
max_ = -std::numeric_limits<Real>::max();
- StatAccumulator acc;
+ StatAccumulator<> acc;
for(ExtentIterator it(fnc.GetExtent());!it.AtEnd(); ++it) {
Real val=fnc.GetReal(it);
diff --git a/modules/img/alg/src/stat_accumulator.cc b/modules/img/alg/src/stat_accumulator.cc
deleted file mode 100644
index c4641d5504cafa83cfb788dec0ad21b59cead46e..0000000000000000000000000000000000000000
--- a/modules/img/alg/src/stat_accumulator.cc
+++ /dev/null
@@ -1,141 +0,0 @@
-//------------------------------------------------------------------------------
-// This file is part of the OpenStructure project <www.openstructure.org>
-//
-// Copyright (C) 2008-2011 by the OpenStructure authors
-//
-// This library is free software; you can redistribute it and/or modify it under
-// the terms of the GNU Lesser General Public License as published by the Free
-// Software Foundation; either version 3.0 of the License, or (at your option)
-// any later version.
-// This library is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
-// details.
-//
-// You should have received a copy of the GNU Lesser General Public License
-// along with this library; if not, write to the Free Software Foundation, Inc.,
-// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
-//------------------------------------------------------------------------------
-
-
-#include "stat_accumulator.hh"
-
-namespace ost { namespace img { namespace alg {
-
-StatAccumulator::StatAccumulator():
- mean_(0.0),
- sum_(0.0),
- sum2_(0.0),
- m2_(0.0),
- m3_(0.0),
- m4_(0.0),
- n_(1)
-
-{
-}
-
-void StatAccumulator::operator()(Real val)
-{
- Real delta = val - mean_;
- Real delta_n = delta / n_;
- Real delta_n2 = delta_n * delta_n;
- Real term = delta * delta_n * (n_-1);
-
- mean_ += delta_n;
- m4_ += term * delta_n2 * (n_*n_ - 3*n_ + 3) + 6 * delta_n2 * m2_ - 4 * delta_n * m3_;
- m3_ += term * delta_n * (n_ - 2) - 3 * delta_n * m2_;
- m2_ += term;
- sum_+=val;
- sum2_+=val*val;
- n_+=1;
-}
-
-StatAccumulator StatAccumulator::operator+(const StatAccumulator& acc2) const
-{
- StatAccumulator acc(acc2);
- acc+=*this;
- return acc;
-}
-
-StatAccumulator& StatAccumulator::operator+=(const StatAccumulator& acc)
-{
- if(acc.n_==1){
- return *this;
- }
- if(n_==1){
- mean_=acc.mean_;
- sum_=acc.sum_;
- sum2_=acc.sum2_;
- m2_=acc.m2_;
- m3_=acc.m3_;
- m4_=acc.m4_;
- n_=acc.n_;
- return *this;
- }
- Real na=n_-1;
- Real nb=acc.n_-1;
- Real nanb=na*nb;
- Real delta = acc.mean_ - mean_;
- Real delta_n = delta / (na+nb);
- Real delta_n2 = delta_n * delta_n;
-
- mean_ += nb*delta_n;
- m4_+=acc.m4_+delta*delta_n*delta_n2*nanb*(na*na-nanb+nb*nb)+6.0*delta_n2*(na*na*acc.m2_+nb*nb*m2_)+4.0*delta_n*(na*acc.m3_-nb*m3_);
- m3_+=acc.m3_+delta*delta_n2*nanb*(na-nb)+3.0*delta_n*(na*acc.m2_-nb*m2_);
- m2_ += acc.m2_+delta*delta_n*nanb;
- sum_+=acc.sum_;
- sum2_+=acc.sum2_;
- n_+=nb;
- return *this;
-}
-
-Real StatAccumulator::GetMean() const
-{
- return mean_;
-}
-
-Real StatAccumulator::GetSum() const
-{
- return sum_;
-}
-
-Real StatAccumulator::GetVariance() const
-{
- if(n_==1.0){
- return 0.0;
- }
- return m2_/(n_-1);
-}
-
-Real StatAccumulator::GetStandardDeviation() const
-{
- return sqrt(GetVariance());
-}
-
-Real StatAccumulator::GetRootMeanSquare() const
-{
- if(n_==1.0){
- return 0.0;
- }
- return sqrt(sum2_/(n_-1));
-}
-
-Real StatAccumulator::GetSkewness() const
-{
- if(m2_<1e-20){
- return 0.0;
- }else{
- return m3_/sqrt(m2_*m2_*m2_);
- }
-}
-
-Real StatAccumulator::GetKurtosis() const
-{
- if(m2_<1e-20){
- return 0.0;
- }else{
- return ((n_-1)*m4_) / (m2_*m2_);
- }
-}
-
-}}} //ns
diff --git a/modules/img/alg/src/stat_accumulator.hh b/modules/img/alg/src/stat_accumulator.hh
index 3cb5f06eb3c9882aa1481b637de626e11322d567..3483890b0ced8c3f5e71382f3273fcdcafc68ef0 100644
--- a/modules/img/alg/src/stat_accumulator.hh
+++ b/modules/img/alg/src/stat_accumulator.hh
@@ -25,21 +25,168 @@
namespace ost { namespace img { namespace alg {
+template<unsigned int MAX_MOMENT=4>
class DLLEXPORT_IMG_ALG StatAccumulator
{
public:
- StatAccumulator();
- void operator()(Real val);
- StatAccumulator operator+(const StatAccumulator& acc) const;
- StatAccumulator& operator+=(const StatAccumulator& acc);
-
- Real GetMean() const ;
- Real GetSum() const ;
- Real GetVariance() const;
- Real GetStandardDeviation() const ;
- Real GetRootMeanSquare() const;
- Real GetSkewness() const ;
- Real GetKurtosis() const;
+ StatAccumulator():
+ mean_(0.0),
+ sum_(0.0),
+ sum2_(0.0),
+ m2_(0.0),
+ m3_(0.0),
+ m4_(0.0),
+ n_(1)
+ {}
+
+ void operator()(Real val)
+ {
+ Real delta,delta_n,delta_n2,term;
+ if(MAX_MOMENT>0){
+ delta = val - mean_;
+ delta_n = delta / n_;
+ }
+ if(MAX_MOMENT>3){
+ delta_n2 = delta_n * delta_n;
+ }
+ if(MAX_MOMENT>1){
+ term = delta * delta_n * (n_-1);
+ }
+ if(MAX_MOMENT>3){
+ m4_ += term * delta_n2 * (n_*n_ - 3*n_ + 3) + 6 * delta_n2 * m2_ - 4 * delta_n * m3_;
+ }
+ if(MAX_MOMENT>2){
+ m3_ += term * delta_n * (n_ - 2) - 3 * delta_n * m2_;
+ }
+ if(MAX_MOMENT>1){
+ m2_ += term;
+ }
+ if(MAX_MOMENT>0){
+ mean_ += delta_n;
+ }
+ n_+=1;
+ sum_+=val;
+ sum2_+=val*val;
+ }
+
+ StatAccumulator operator+(const StatAccumulator& acc2) const
+ {
+ StatAccumulator acc(acc2);
+ acc+=*this;
+ return acc;
+ }
+
+ StatAccumulator& operator+=(const StatAccumulator& acc)
+ {
+ if(acc.n_==1){
+ return *this;
+ }
+ if(n_==1){
+ mean_=acc.mean_;
+ sum_=acc.sum_;
+ sum2_=acc.sum2_;
+ m2_=acc.m2_;
+ m3_=acc.m3_;
+ m4_=acc.m4_;
+ n_=acc.n_;
+ return *this;
+ }
+ Real delta,delta_n,delta_n2,na,nanb;
+ Real nb=acc.n_-1;
+ if(MAX_MOMENT>0){
+ na=n_-1;
+ delta = acc.mean_ - mean_;
+ delta_n = delta / (na+nb);
+ }
+ if(MAX_MOMENT>1){
+ nanb=na*nb;
+ }
+ if(MAX_MOMENT>2){
+ delta_n2 = delta_n * delta_n;
+ }
+ if(MAX_MOMENT>3){
+ m4_+=acc.m4_+delta*delta_n*delta_n2*nanb*(na*na-nanb+nb*nb)+6.0*delta_n2*(na*na*acc.m2_+nb*nb*m2_)+4.0*delta_n*(na*acc.m3_-nb*m3_);
+ }
+ if(MAX_MOMENT>2){
+ m3_+=acc.m3_+delta*delta_n2*nanb*(na-nb)+3.0*delta_n*(na*acc.m2_-nb*m2_);
+ }
+ if(MAX_MOMENT>1){
+ m2_ += acc.m2_+delta*delta_n*nanb;
+ }
+ if(MAX_MOMENT>0){
+ mean_ += nb*delta_n;
+ }
+ sum_+=acc.sum_;
+ sum2_+=acc.sum2_;
+ n_+=nb;
+ return *this;
+ }
+
+ Real GetNumSamples() const
+ {
+ return n_-1.0;
+ }
+
+ Real GetMean() const
+ {
+ if(MAX_MOMENT<1){
+ throw Error("Mean was not calculated.");
+ }
+ return mean_;
+ }
+
+ Real GetSum() const
+ {
+ return sum_;
+ }
+
+ Real GetVariance() const
+ {
+ if(MAX_MOMENT<2){
+ throw Error("Variance was not calculated.");
+ }
+ if(n_==1.0){
+ return 0.0;
+ }
+ return m2_/(n_-1);
+ }
+
+ Real GetStandardDeviation() const
+ {
+ return sqrt(GetVariance());
+ }
+
+ Real GetRootMeanSquare() const
+ {
+ if(n_==1.0){
+ return 0.0;
+ }
+ return sqrt(sum2_/(n_-1));
+ }
+
+ Real GetSkewness() const
+ {
+ if(MAX_MOMENT<3){
+ throw Error("Skewness was not calculated.");
+ }
+ if(m2_<1e-20){
+ return 0.0;
+ }else{
+ return m3_/sqrt(m2_*m2_*m2_);
+ }
+ }
+
+ Real GetKurtosis() const
+ {
+ if(MAX_MOMENT<4){
+ throw Error("Kurtosis was not calculated.");
+ }
+ if(m2_<1e-20){
+ return 0.0;
+ }else{
+ return ((n_-1)*m4_) / (m2_*m2_);
+ }
+ }
private:
Real mean_;
diff --git a/modules/img/alg/tests/test_stat.cc b/modules/img/alg/tests/test_stat.cc
index 592d00c9a907d591be0999790490471af88ce4da..9c9694a4caace6c7b79c14fd9d2fc2c9a1b57778 100644
--- a/modules/img/alg/tests/test_stat.cc
+++ b/modules/img/alg/tests/test_stat.cc
@@ -64,7 +64,7 @@ void test() {
BOOST_CHECK_CLOSE(stat.GetKurtosis(),Real(1.77),Real(0.01));
- StatAccumulator acc;
+ StatAccumulator<> acc;
for(int u=0;u<3;++u) {
for(int v=0;v<3;++v) {
acc(val[u][v]);
@@ -75,13 +75,13 @@ void test() {
BOOST_CHECK_CLOSE(acc.GetSkewness()+0.5,Real(0.5),Real(0.0001));
BOOST_CHECK_CLOSE(acc.GetKurtosis(),Real(1.77),Real(0.0001));
- StatAccumulator acc1,acc2,acc3;
+ StatAccumulator<> acc1,acc2,acc3;
for(int u=0;u<3;++u) {
acc1(val[u][0]);
acc2(val[u][1]);
acc3(val[u][2]);
}
- StatAccumulator acc_c=acc1+acc2+acc3;
+ StatAccumulator<> acc_c=acc1+acc2+acc3;
BOOST_CHECK_CLOSE(acc_c.GetMean(),Real(5.0),Real(0.0001));
BOOST_CHECK_CLOSE(acc_c.GetStandardDeviation(),Real(2.58198889747),Real(0.0001));
BOOST_CHECK_CLOSE(acc_c.GetSkewness()+0.5,Real(0.5),Real(0.0001));