templated StatAccumulator to allow selection of calculated moments

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

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);

modules/img/alg/src/stat_accumulator.cc

-//------------------------------------------------------------------------------
-// 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

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_;

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));