diff --git a/modules/img/alg/src/CMakeLists.txt b/modules/img/alg/src/CMakeLists.txt
index 11fba15d4e3ffc675692c8ad0d73bb3c1738a5b7..c181156dbfe99eb4ac8b6a4532fa359b8b161d45 100644
--- a/modules/img/alg/src/CMakeLists.txt
+++ b/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
diff --git a/modules/img/alg/src/stat.cc b/modules/img/alg/src/stat.cc
index 2a9427abbaea05fdc290d0642ed6918ea6b455dd..f433cd6d69a7810a2e34ca043ad3257720cb9acc 100644
--- a/modules/img/alg/src/stat.cc
+++ b/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{
diff --git a/modules/img/alg/src/stat_accumulator.cc b/modules/img/alg/src/stat_accumulator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c4641d5504cafa83cfb788dec0ad21b59cead46e
--- /dev/null
+++ b/modules/img/alg/src/stat_accumulator.cc
@@ -0,0 +1,141 @@
+//------------------------------------------------------------------------------
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..3cb5f06eb3c9882aa1481b637de626e11322d567
--- /dev/null
+++ b/modules/img/alg/src/stat_accumulator.hh
@@ -0,0 +1,55 @@
+//------------------------------------------------------------------------------
+// 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
diff --git a/modules/img/alg/tests/test_stat.cc b/modules/img/alg/tests/test_stat.cc
index 3f5e6d2ed6a3a1114f59458361021b7f9a0ccdad..592d00c9a907d591be0999790490471af88ce4da 100644
--- a/modules/img/alg/tests/test_stat.cc
+++ b/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