00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OST_STAT_ACCUMULATOR_HH
00021 #define OST_STAT_ACCUMULATOR_HH
00022
00023 #include <limits>
00024 #include <boost/math/special_functions/binomial.hpp>
00025 #include <ost/base.hh>
00026 #include <ost/message.hh>
00027 #include <ost/img/alg/module_config.hh>
00028
00029 namespace ost { namespace img { namespace alg {
00030
00031 namespace {
00032
00033 template<class D>
00034 struct MinMax{
00035 static bool less_cmp_(const D& v1, const D& v2){return v1<v2;}
00036 };
00037 template<>
00038 struct MinMax<Complex>{
00039 static bool less_cmp_(const Complex& v1, const Complex& v2){return abs(v1)<abs(v2);}
00040 };
00041
00042 }
00043
00044 template<unsigned int MAX_MOMENT=4,typename DATATYPE=Real>
00045 class StatAccumulator
00046 {
00047 public:
00048 StatAccumulator():
00049 sum_(0.0),
00050 sum2_(0.0),
00051 max_(-std::numeric_limits<DATATYPE>::max()),
00052 min_(std::numeric_limits<DATATYPE>::max()),
00053 m_(),
00054 w_(0.0),
00055 n_(0)
00056 {
00057 for(unsigned int i=0;i<MAX_MOMENT;++i){
00058 m_[i]=0.0;
00059 }
00060
00061 }
00062
00063 StatAccumulator(DATATYPE val, Real w=1.0):
00064 sum_(val),
00065 sum2_(val*val),
00066 max_(val),
00067 min_(val),
00068 m_(),
00069 w_(w),
00070 n_(1)
00071 {
00072 if(MAX_MOMENT>0){
00073 m_[0]=val;
00074 for(unsigned int i=1;i<MAX_MOMENT;++i){
00075 m_[i]=0.0;
00076 }
00077 }
00078 }
00079
00080 void operator()(DATATYPE val, Real w=1.0)
00081 {
00082 *this+=StatAccumulator(val,w);
00083 }
00084
00085 StatAccumulator<MAX_MOMENT,DATATYPE> operator+(const StatAccumulator<MAX_MOMENT,DATATYPE>& acc2) const
00086 {
00087 StatAccumulator<MAX_MOMENT,DATATYPE> acc(acc2);
00088 acc+=*this;
00089 return acc;
00090 }
00091
00092 StatAccumulator<MAX_MOMENT,DATATYPE>& operator+=(const StatAccumulator<MAX_MOMENT,DATATYPE>& acc)
00093 {
00094 if(0.0>=w_){
00095 *this=acc;
00096 return *this;
00097 }
00098 if(0.0>=acc.w_){
00099 return *this;
00100 }
00101 n_+=acc.n_;
00102 Real wn=w_+acc.w_;
00103
00104 sum_+=acc.sum_;
00105 sum2_+=acc.sum2_;
00106 max_=std::max<DATATYPE>(max_,acc.max_,MinMax<DATATYPE>::less_cmp_);
00107 min_=std::min<DATATYPE>(min_,acc.min_,MinMax<DATATYPE>::less_cmp_);
00108 if(MAX_MOMENT>0){
00109 DATATYPE delta=acc.m_[0]-m_[0];
00110 DATATYPE delta_w=delta/wn;
00111 if(MAX_MOMENT>2){
00112 DATATYPE w1w2_delta_w=w_*acc.w_*delta_w;
00113 DATATYPE w1w2_delta_wp=w1w2_delta_w*w1w2_delta_w;
00114 Real iw2=1.0/acc.w_;
00115 Real iw2pm1=iw2;
00116 Real miw1=-1.0/w_;
00117 Real miw1pm1=miw1;
00118 DATATYPE mn[MAX_MOMENT];
00119 for(unsigned int p=3;p<=MAX_MOMENT;++p){
00120 w1w2_delta_wp*=w1w2_delta_w;
00121 iw2pm1*=iw2;
00122 miw1pm1*=miw1;
00123 DATATYPE delta_wk=1.0;
00124 DATATYPE s=0.0;
00125 Real mw2k=1.0;
00126 Real w1k=1.0;
00127 for(unsigned int k=1;k<=p-2;++k){
00128 w1k*=w_;
00129 mw2k*=-acc.w_;
00130 delta_wk*=delta_w;
00131 s+=boost::math::binomial_coefficient<Real>(p, k)*(mw2k*m_[p-k-1]+w1k*acc.m_[p-k-1])*delta_wk;
00132 }
00133 mn[p-3]=acc.m_[p-1]+s+w1w2_delta_wp*(iw2pm1-miw1pm1);
00134 }
00135 for(unsigned int p=3;p<=MAX_MOMENT;++p){
00136 m_[p-1]+=mn[p-3];
00137 }
00138 }
00139 if(MAX_MOMENT>1){
00140 m_[1]+=acc.m_[1]+delta_w*delta*acc.w_*w_;
00141 }
00142 m_[0]+=delta_w*acc.w_;
00143 w_=wn;
00144 }
00145 return *this;
00146 }
00147
00148 Real GetNumSamples() const
00149 {
00150 return n_;
00151 }
00152
00153 DATATYPE GetMaximum() const
00154 {
00155 return max_;
00156 }
00157
00158 DATATYPE GetMinimum() const
00159 {
00160 return min_;
00161 }
00162
00163 Real GetWeight() const
00164 {
00165 return w_;
00166 }
00167
00168 DATATYPE GetMoment(unsigned int i) const
00169 {
00170 if(1>i){
00171 throw Error("Invalid moment.");
00172 }
00173 if(MAX_MOMENT<i){
00174 throw Error("Moment was not calculated.");
00175 }
00176 return m_[i-1];
00177 }
00178
00179 DATATYPE GetMean() const
00180 {
00181 if(MAX_MOMENT<1){
00182 throw Error("Mean was not calculated.");
00183 }
00184 return m_[0];
00185 }
00186
00187 DATATYPE GetSum() const
00188 {
00189 return sum_;
00190 }
00191
00192 DATATYPE GetVariance() const
00193 {
00194 if(MAX_MOMENT<2){
00195 throw Error("Variance was not calculated.");
00196 }
00197 if(n_==0.0){
00198 return 0.0;
00199 }
00200 return m_[1]/(w_);
00201 }
00202
00203 DATATYPE GetStandardDeviation() const
00204 {
00205 return sqrt(GetVariance());
00206 }
00207
00208 DATATYPE GetRootMeanSquare() const
00209 {
00210 if(n_==0.0){
00211 return 0.0;
00212 }
00213 return sqrt(sum2_/(w_));
00214 }
00215
00216 DATATYPE GetSkewness() const
00217 {
00218 if(MAX_MOMENT<3){
00219 throw Error("Skewness was not calculated.");
00220 }
00221 if(m_[1]<1e-20){
00222 return 0.0;
00223 }else{
00224 return m_[2]/sqrt(m_[1]*m_[1]*m_[1]);
00225 }
00226 }
00227
00228 DATATYPE GetKurtosis() const
00229 {
00230 if(MAX_MOMENT<4){
00231 throw Error("Kurtosis was not calculated.");
00232 }
00233 if(m_[1]<1e-20){
00234 return 0.0;
00235 }else{
00236 return w_*m_[3] / (m_[1]*m_[1]);
00237 }
00238 }
00239
00240 private:
00241 DATATYPE sum_;
00242 DATATYPE sum2_;
00243 DATATYPE max_;
00244 DATATYPE min_;
00245 DATATYPE m_[MAX_MOMENT];
00246 Real w_;
00247 unsigned int n_;
00248 };
00249
00250 }}}
00251 #endif // OST_STAT_ACCUMULATOR_HH