diff --git a/modules/img/alg/src/stat.cc b/modules/img/alg/src/stat.cc
index 8f43ca399e47fb9e6dac913624bd2ab020c4ec9e..2a9427abbaea05fdc290d0642ed6918ea6b455dd 100644
--- a/modules/img/alg/src/stat.cc
+++ b/modules/img/alg/src/stat.cc
@@ -41,11 +41,27 @@ template <typename T, class D>
void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
{
sum_=0.0;
- Real sum2=0.0,sum3=0.0,sum4=0.0;
+ mean_=0.0;
+
Real n = static_cast<Real>(isi.GetSize().GetVol());
- Vec3 sumcenter(0.0,0.0,0.0);
+ if(n==0.0){
+ var_=0.0;
+ std_dev_=0.0;
+ min_=0.0;
+ max_=0.0;
+ maxpos_=Point(0,0,0),
+ minpos_=Point(0,0,0),
+ rms_=0.0;
+ skewness_=0.0;
+ kurtosis_=0.0;
+ center_of_mass_=Vec3(0.0,0.0,0.0);
+ return;
+ }
- if(n==0.0) 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();
@@ -58,18 +74,23 @@ 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 oval=val;
+ 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;
- val*=oval;
- sum2+=val;
- val*=oval;
- sum3+=val;
- val*=oval;
- sum4+=val;
+ sum2+=val*val;
+ current_n+=1;
}
}
}
@@ -77,16 +98,15 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
max_=maxindex.first;
minpos_=minindex.second+isi.GetExtent().GetStart();
maxpos_=maxindex.second+isi.GetExtent().GetStart();
- mean_ = sum_/n;
- var_=sum2/n-mean_*mean_;
+ var_=m2/n;
std_dev_=sqrt(var_);
rms_=sqrt(sum2/n);
- if(std_dev_>0.0){
- skewness_=sqrt(n)*(sum3-3.0*mean_*sum2+3.0*mean_*mean_*sum_ -n*mean_*mean_*mean_)/(std_dev_*std_dev_*std_dev_);
- kurtosis_= n*(sum4-4.0*sum3*mean_+6.0*sum2*mean_*mean_-4.0*sum_*mean_*mean_*mean_+n*mean_*mean_*mean_*mean_)/(var_*var_)-3.0;
- }else{
+ if(m2<1e20){
skewness_=0.0;
- kurtosis_=0.0;
+ kurtosis_= 0.0;
+ }else{
+ skewness_=m3/sqrt(m2*m2*m2);
+ kurtosis_= (n*m4) / (m2*m2) - 3.0;
}
if(sum_!=0.0){
center_of_mass_=sumcenter/sum_+isi.GetExtent().GetStart().ToVec3();
@@ -98,46 +118,65 @@ void StatBase::VisitState(const ImageStateImpl<T,D>& isi)
void StatBase::VisitFunction(const Function& fnc)
{
sum_=0.0;
- Real sum2=0.0,sum3=0.0,sum4=0.0;
- ValIndex minindex(std::numeric_limits<Real>::max(),Point(0,0,0));
- ValIndex maxindex(-std::numeric_limits<Real>::max(),Point(0,0,0));
- Real n = (Real)(fnc.GetSize().GetVol());
- Vec3 sumcenter(0.0,0.0,0.0);
+ mean_=0.0;
- if(n==0.0) return;
+ Real n = (Real)(fnc.GetSize().GetVol());
+ if(n==0.0){
+ var_=0.0;
+ std_dev_=0.0;
+ min_=0.0;
+ max_=0.0;
+ maxpos_=Point(0,0,0),
+ minpos_=Point(0,0,0),
+ rms_=0.0;
+ skewness_=0.0;
+ kurtosis_=0.0;
+ center_of_mass_=Vec3(0.0,0.0,0.0);
+ 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();
for(ExtentIterator it(fnc.GetExtent());!it.AtEnd(); ++it) {
Real val=fnc.GetReal(it);
- Real oval=val;
+ 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;
- val*=oval;
- sum2+=val;
- val*=oval;
- sum3+=val;
- val*=oval;
- sum4+=val;
+ sum2+=val*val;
+ current_n+=1;
}
min_=minindex.first;
max_=maxindex.first;
minpos_=minindex.second;
maxpos_=maxindex.second;
- mean_ = sum_/n;
- var_=sum2/n-mean_*mean_;
+ var_=m2/n;
std_dev_=sqrt(var_);
rms_=sqrt(sum2/n);
- if(std_dev_>0.0){
- skewness_=sqrt(n)*(sum3-3*mean_*sum2+3*mean_*mean_*sum_-n*mean_*mean_*mean_)/(std_dev_*std_dev_*std_dev_);
- kurtosis_= n*(sum4-4*sum3*mean_+6*sum2*mean_*mean_-4*sum_*mean_*mean_*mean_+n*mean_*mean_*mean_*mean_)/(var_*var_)-3;
- }else{
+ if(m2<1e20){
skewness_=0.0;
- kurtosis_=0.0;
+ kurtosis_= 0.0;
+ }else{
+ skewness_=m3/sqrt(m2*m2*m2);
+ kurtosis_= (n*m4) / (m2*m2) - 3.0;
}
if(sum_!=0.0){
center_of_mass_=sumcenter/sum_;
diff --git a/modules/img/alg/src/stat.hh b/modules/img/alg/src/stat.hh
index 82757dc728e0fb8225055284eed646eeb7ac983f..450ffc3528de44dd2eb0df9d1d9cb7fbd2dc3de7 100644
--- a/modules/img/alg/src/stat.hh
+++ b/modules/img/alg/src/stat.hh
@@ -1,7 +1,7 @@
//------------------------------------------------------------------------------
// This file is part of the OpenStructure project <www.openstructure.org>
//
-// Copyright (C) 2008-2011 by the OpenStructure authors
+// Copyright (C) 2008-2010 by the OpenStructure authors
// Copyright (C) 2003-2010 by the IPLT authors
//
// This library is free software; you can redistribute it and/or modify it under
@@ -39,6 +39,11 @@ namespace ost { namespace img { namespace alg {
Since this algorithm is implemented as a combined image stage visitor
and algorithm, the main workhorse is this class StatBase, which will
act as the parent class of the actual algorithm class, Stat
+
+ Mean value, variance and standard deviation are calculated based on the one pass algorithm by Welford et al.:
+ B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420
+ The calculation of the higher order central moments is implemented according to Terriberry:
+ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online (http://people.xiph.org/~tterribe/notes/homs.html)
*/
class DLLEXPORT_IMG_ALG StatBase