Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

ImageStats.hpp

Go to the documentation of this file.
00001 /* ImLib3D
00002  * Copyright (c) 2001, ULP-IPB Strasbourg.
00003  *
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2 of the License, or (at
00007  * your option) any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017  */
00019 
00020 
00021 #ifndef _ImageStats_hpp
00022 #define _ImageStats_hpp
00023 
00024 //#include<ImLib3D/Image3D.hpp>
00025 #include<ImLib3D/ImageProcessor.hpp>
00026 #include<map>
00027 #include<algorithm>
00028 
00029 
00030 namespace IP3D
00031 {
00039     template <class ImageType>
00040     void Average(const ImageType& src,typename ImageType::value_type &average)
00041     {
00042         typename ImageType::const_iteratorFastMasked p;
00043         average= Zero<typename ImageType::value_type>();
00044         int ct=0;
00045         for(p=src.begin();p!=src.end();++p){ct++;average+=*p;}
00046         average/=ct;
00047     }
00048 };
00049 namespace IP3D
00050 {
00059     template <class ImageType>
00060     void AverageAndVariance(const ImageType& src,double &average,double &variance)
00061     {
00062         typename ImageType::const_iteratorFastMasked p;
00063         average=0;variance=0;
00064         int ct=0;
00065         for(p=src.begin();p!=src.end();++p){ct++;average+=*p;variance+=(*p)*(*p);}
00066         average/=ct;
00067         variance=variance/ct-average*average;
00068     }
00069 };
00070 
00071 namespace IP3D
00072 {
00082         void RobustAverageAndVariance(const Image3Df &src,double &average,double &variance);
00083 };
00084 
00085 
00087 
00089 namespace IP3D
00090 {
00101     template <class ImageType>
00102     void Median(const ImageType& src,double &res,double part=.5,double precision=-.1)
00103     {
00104         
00105         vector<typename ImageType::value_type> values;
00106         values.reserve(src.GetNVoxels());
00107         for(typename ImageType::const_iteratorFastMasked p=src.begin();p!=src.end();++p)
00108         {
00109             values.push_back(*p);
00110         }
00111         int n=values.size();
00112         sort(values.begin(),values.end());
00113 //          mprintf("CMedian result:%f %f %f %f %f\n",values[(int)(n*0)],values[(int)(n*.3)],
00114 //                 values[(int)(n*.5)],values[(int)(n*.66)],values[(int)(n-1)]);
00115         res=(double)values[(int)(n*part)];
00116     }
00117 };
00118 
00120 class ProbabilityDistribution 
00121 {
00122     typedef  float Value;
00123     vector<double> bins;
00124     float vmin,vmax;    
00125 
00126     long int accumulatedSegmentCount; // keep count of accumulated segments in order to accumulate further signals
00127 
00128     // cumulative prob distribution
00129     mutable bool cumulIsComputed;
00130     mutable vector<double> cumul;
00131     mutable bool cumulInverseIsComputed;
00132     mutable map<double,int> cumulInverse;
00133 
00134     void SetRange(const vector<Value> &samples);
00135 
00136 public:
00137     static int ScottsRule(const vector<float>& vals); // for finding number of bins automatically
00138 
00139     vector<double> GetBins() const { return bins;}  
00140     int GetNBins() const {return bins.size();}
00141     void GetRange(float& rmin, float& rmax) const {rmin=vmin;rmax=vmax;}
00142     double GetMax() const {return *max_element(bins.begin(), bins.end());}
00143     double GetMaxPos() const {return BinValue(max_element(bins.begin(), bins.end())-bins.begin());}
00144     double operator[](int bin) const {return bins[bin];}
00145     double Density(int bin) const {return bins[bin]*bins.size()/(vmax-vmin);}
00147     void SetBins(float _vmin, float _vmax, int nbBins); 
00149     double BinValue(int binNum) const 
00150     {
00151         return vmin+binNum*(vmax-vmin)/bins.size();
00152     }
00153 
00155     double BinPos(Value v) const 
00156     {
00157         return (bins.size()*(v-vmin)/(vmax-vmin));
00158     }
00159 
00161     int BinNum(Value v) const 
00162     {
00163         // should use rint!
00164         int res=(int)(BinPos(v));
00165         if(res>=(int)bins.size() || res<0){ThrowError("ProbabilityDistribution::BinNum out of range:%f-> %d (min=%f max=%f)\n",v,res,vmin,vmax);}
00166         return res;
00167     }
00168 protected:
00169     void ComputeCumul() const 
00170     {
00171         cumul.clear();
00172         cumul.assign(bins.size(),0);
00173         double v=0;
00174         for(uint i=0;i<bins.size();i++)
00175         {
00176             v+=bins[i];
00177             cumul[i]=v;
00178         }
00179         cumulIsComputed=true;
00180     }
00181     void ComputeInverseCumul() const
00182     {       
00183         cumulInverse.clear();
00184         if(!cumulIsComputed){ComputeCumul();}
00185         for(uint i=0;i<cumul.size();i++)
00186         {
00187             cumulInverse[cumul[i]]=i;
00188         }
00189         cumulInverseIsComputed=true;
00190     }
00191 public:
00193     double CumulInverse(double v) const 
00194     {
00195         if(!cumulInverseIsComputed){ComputeInverseCumul();}
00196         map<double,int>::iterator res=cumulInverse.upper_bound(v);
00197         if(res==cumulInverse.end()){return vmax;}
00198         int cbin=(*res).second-1;
00199         return cbin*(vmax-vmin)/bins.size()+vmin;
00200     }
00202     double Cumul(double v) const 
00203     {
00204         if (!cumulIsComputed){ComputeCumul();}
00205         if (v <= vmin) return 0.0;
00206         if (v >= vmax) return 1.0;
00207         return cumul[BinNum(v)] - (BinValue(BinNum(v)+1) - BinValue(BinNum(v)))/(vmax-vmin)*bins.size()*bins[BinNum(v)];
00208     }
00209 
00210     bool NextIndex(vector<int> &indexes) 
00211     {
00212         for(uint i=0;i<indexes.size();i++)
00213         {
00214             indexes[i]++;
00215             if(indexes[i]<(int)bins.size()){return true;}
00216             indexes[i]=0;
00217         }
00218         return false;
00219     }
00220 
00221     class OtsuClass
00222     {
00223     public:
00224         int id;
00225         double moment0;
00226         double moment1;
00227         double moment2 ;
00228         int j0;
00229         int j1;
00230         OtsuClass():
00231             id(-1),moment0(0),moment1(0),moment2(0),j0(0),j1(0)
00232         {;}
00233         
00234     };
00235 
00236     OtsuClass ComputeClassStats(vector<OtsuClass> &classes,int classNum,vector<int> &indexes);
00237     void OtsuThreshholds(int nbClasses,vector<double> &res);
00238 
00240     void SetWithSamples(const vector<Value> &samples,int nbBins=1000);
00242     void AccumulateWithSamples(const vector<Value> &samples);
00243 
00245     void SetWithSignal(const vector<Value> &samples,int nbBins=1000);
00247     void AccumulateWithSignal(const vector<Value> &samples);
00248 
00250     void AccumulateImage3Dlinear3(const Image3Df &src);
00251     void SetWithImage3Dlinear3(const Image3Df &src,int nbBins);
00252 
00254     enum SetWithImage3DType {LINEAR3,LINEAR1,SAMPLES,FAST=LINEAR1,FASTEST=SAMPLES};
00256     void SetWithImage3D(const Image3Df &src,SetWithImage3DType type=LINEAR3,int nbBins=1000);
00257 
00259     void AccumulateWithImage3D(const Image3Df &src,SetWithImage3DType type=LINEAR3);
00260 
00262     void ShowWithGnuPlot(const string &label="") const ;
00263     
00264     ProbabilityDistribution():accumulatedSegmentCount(0),cumulIsComputed(false),cumulInverseIsComputed(false)
00265     {}
00266 };
00267 
00268 #endif // _ImageStats_hpp
00269 

Generated on Fri Jun 17 13:36:04 2005 for ImLib3D by  doxygen 1.4.2