encapsulated stat functionality into separate accumulator to make it easier to re-use

modules/img/alg/src/CMakeLists.txt

@@ -25,6 +25,7 @@ power_spectrum.cc
 randomize.cc
 smooth_mask_image.cc
 stat.cc
+stat_accumulator.cc
 stat_min_max.cc
 alg_transform.cc
 alg_mirror.cc
@@ -75,6 +76,7 @@ power_spectrum.hh
 randomize.hh
 smooth_mask_image.hh
 stat.hh
+stat_accumulator.hh
 stat_min_max.hh
 threshold.hh
 transcendentals.hh

modules/img/alg/src/stat.cc

@@ -30,6 +30,7 @@
 #include <ost/img/value_util.hh>
 
 #include "stat.hh"
+#include "stat_accumulator.hh"
 
 namespace ost { namespace img { namespace alg {
 
@@ -58,14 +59,12 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
     return;
   }
 
-  Real sum2=0.0;
   Vec3 sumcenter(0.0,0.0,0.0);
-  Real current_n=1.0;
-  Real m2=0.0,m3=0.0,m4=0.0;
   ValIndex minindex(std::numeric_limits<Real>::max(),Point(0,0,0));
   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;
 
   int wi=isi.GetSize()[0];
   int he=isi.GetSize()[1];
@@ -74,23 +73,11 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
     for(int j=0;j<he;++j) {
       for(int k=0;k<de;++k) {
         Real val=Val2Val<T,Real>(isi.Value(Index(i,j,k)));
-        Real delta = val - mean_;
-        Real delta_n = delta / current_n;
-        Real delta_n2 = delta_n * delta_n;
-        Real term = delta * delta_n * (current_n-1);
-
-        mean_ += delta_n;
-        m4 += term * delta_n2 * (current_n*current_n - 3*current_n + 3) + 6 * delta_n2 * m2 - 4 * delta_n * m3;
-        m3 += term * delta_n * (current_n - 2) - 3 * delta_n * m2;
-        m2 += term;
-
         ValIndex vi(val,Point(i,j,k));
         minindex=std::min<ValIndex>(vi,minindex);
         maxindex=std::max<ValIndex>(vi,maxindex);
         sumcenter+=Vec3(i,j,k)*val;
-        sum_+=val;
-        sum2+=val*val;
-        current_n+=1;
+        acc(val);
       }
     }
   }
@@ -98,16 +85,13 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
   max_=maxindex.first;
   minpos_=minindex.second+isi.GetExtent().GetStart();
   maxpos_=maxindex.second+isi.GetExtent().GetStart();  
-  var_=m2/n;
-  std_dev_=sqrt(var_);
-  rms_=sqrt(sum2/n);
-  if(m2<1e20){
-    skewness_=0.0;
-    kurtosis_= 0.0;
-  }else{
-    skewness_=m3/sqrt(m2*m2*m2);
-    kurtosis_= (n*m4) / (m2*m2) - 3.0;
-  }
+  var_=acc.GetVariance();
+  std_dev_=acc.GetStandardDeviation();
+  rms_=acc.GetRootMeanSquare();
+  skewness_=acc.GetSkewness();
+  kurtosis_= acc.GetKurtosis();
+  sum_=acc.GetSum();
+  mean_=acc.GetMean();
   if(sum_!=0.0){
     center_of_mass_=sumcenter/sum_+isi.GetExtent().GetStart().ToVec3();
   }else{
@@ -135,49 +119,32 @@ void StatBase::VisitFunction(const Function& fnc)
     return;
   }
 
-  Real sum2=0.0;
   Vec3 sumcenter(0.0,0.0,0.0);
-  Real current_n=1.0;
-  Real m2=0.0,m3=0.0,m4=0.0;
   ValIndex minindex(std::numeric_limits<Real>::max(),Point(0,0,0));
   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;
 
   for(ExtentIterator it(fnc.GetExtent());!it.AtEnd(); ++it) {
     Real val=fnc.GetReal(it);
-    Real delta = val - mean_;
-    Real delta_n = delta / current_n;
-    Real delta_n2 = delta_n * delta_n;
-    Real term = delta * delta_n * (current_n-1);
-
-    mean_ += delta_n;
-    m4 += term * delta_n2 * (current_n*current_n - 3*current_n + 3) + 6 * delta_n2 * m2 - 4 * delta_n * m3;
-    m3 += term * delta_n * (current_n - 2) - 3 * delta_n * m2;
-    m2 += term;
-
     ValIndex vi(val,it);
     minindex=std::min<ValIndex>(vi,minindex);
     maxindex=std::max<ValIndex>(vi,maxindex);
     sumcenter+=Point(it).ToVec3()*val;
-    sum_+=val;
-    sum2+=val*val;
-    current_n+=1;
+    acc(val);
   }
   min_=minindex.first;
   max_=maxindex.first;
   minpos_=minindex.second;
-  maxpos_=maxindex.second;  
-  var_=m2/n;
-  std_dev_=sqrt(var_);
-  rms_=sqrt(sum2/n);
-  if(m2<1e20){
-    skewness_=0.0;
-    kurtosis_= 0.0;
-  }else{
-    skewness_=m3/sqrt(m2*m2*m2);
-    kurtosis_= (n*m4) / (m2*m2) - 3.0;
-  }
+  maxpos_=maxindex.second;
+  var_=acc.GetVariance();
+  std_dev_=acc.GetStandardDeviation();
+  rms_=acc.GetRootMeanSquare();
+  skewness_=acc.GetSkewness();
+  kurtosis_= acc.GetKurtosis();
+  sum_=acc.GetSum();
+  mean_=acc.GetMean();
   if(sum_!=0.0){
     center_of_mass_=sumcenter/sum_;
   }else{

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

+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+
+#ifndef OST_STAT_ACCUMULATOR_HH
+#define OST_STAT_ACCUMULATOR_HH
+
+#include <ost/base.hh>
+#include <ost/img/alg/module_config.hh>
+
+namespace ost { namespace img { namespace alg {
+
+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;
+
+private:
+  Real mean_;
+  Real sum_;
+  Real sum2_;
+  Real m2_;
+  Real m3_;
+  Real m4_;
+  Real n_;
+};
+
+}}} //ns
+#endif // OST_STAT_ACCUMULATOR_HH

modules/img/alg/tests/test_stat.cc

@@ -29,6 +29,7 @@
 #include <ost/img/image.hh>
 
 #include <ost/img/alg/stat.hh>
+#include <ost/img/alg/stat_accumulator.hh>
 
 namespace test_stat {
 
@@ -51,7 +52,40 @@ void test() {
   im.Apply(stat);
   BOOST_CHECK_CLOSE(stat.GetMean(),Real(5.0),Real(0.0001));
   BOOST_CHECK_CLOSE(stat.GetStandardDeviation(),Real(2.58198889747),Real(0.0001));
-  
+  BOOST_CHECK_CLOSE(stat.GetSkewness()+0.5,Real(0.5),Real(0.0001));
+  BOOST_CHECK_CLOSE(stat.GetKurtosis(),Real(1.77),Real(0.0001));
+
+  // check for rounding errors
+  im+=10000.0;
+  im.Apply(stat);
+  BOOST_CHECK_CLOSE(stat.GetMean(),Real(10005.0),Real(0.0001));
+  BOOST_CHECK_CLOSE(stat.GetStandardDeviation(),Real(2.58198889747),Real(0.01));
+  BOOST_CHECK_CLOSE(stat.GetSkewness()+0.5,Real(0.5),Real(0.01));
+  BOOST_CHECK_CLOSE(stat.GetKurtosis(),Real(1.77),Real(0.01));
+
+
+  StatAccumulator acc;
+  for(int u=0;u<3;++u) {
+    for(int v=0;v<3;++v) {
+      acc(val[u][v]);
+    }
+  }
+  BOOST_CHECK_CLOSE(acc.GetMean(),Real(5.0),Real(0.0001));
+  BOOST_CHECK_CLOSE(acc.GetStandardDeviation(),Real(2.58198889747),Real(0.0001));
+  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;
+  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;
+  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));
+  BOOST_CHECK_CLOSE(acc_c.GetKurtosis(),Real(1.77),Real(0.0001));
 }
 
 } // namespace