ImageStats.hpp

Go to the documentation of this file.

/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program 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
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


#ifndef _ImageStats_hpp
#define _ImageStats_hpp

//#include<ImLib3D/Image3D.hpp>
#include<ImLib3D/ImageProcessor.hpp>
#include<map>
#include<algorithm>


namespace IP3D
{
    template <class ImageType>
    void Average(const ImageType& src,typename ImageType::value_type &average)
    {
        typename ImageType::const_iteratorFastMasked p;
        average= Zero<typename ImageType::value_type>();
        int ct=0;
        for(p=src.begin();p!=src.end();++p){ct++;average+=*p;}
        average/=ct;
    }
};
namespace IP3D
{
    template <class ImageType>
    void AverageAndVariance(const ImageType& src,double &average,double &variance)
    {
        typename ImageType::const_iteratorFastMasked p;
        average=0;variance=0;
        int ct=0;
        for(p=src.begin();p!=src.end();++p){ct++;average+=*p;variance+=(*p)*(*p);}
        average/=ct;
        variance=variance/ct-average*average;
    }
};

namespace IP3D
{
        void RobustAverageAndVariance(const Image3Df &src,double &average,double &variance);
};



namespace IP3D
{
    template <class ImageType>
    void Median(const ImageType& src,double &res,double part=.5,double precision=-.1)
    {
        
        vector<typename ImageType::value_type> values;
        values.reserve(src.GetNVoxels());
        for(typename ImageType::const_iteratorFastMasked p=src.begin();p!=src.end();++p)
        {
            values.push_back(*p);
        }
        int n=values.size();
        sort(values.begin(),values.end());
//          mprintf("CMedian result:%f %f %f %f %f\n",values[(int)(n*0)],values[(int)(n*.3)],
//                 values[(int)(n*.5)],values[(int)(n*.66)],values[(int)(n-1)]);
        res=(double)values[(int)(n*part)];
    }
};

class ProbabilityDistribution 
{
    typedef  float Value;
    vector<double> bins;
    float vmin,vmax;    

    long int accumulatedSegmentCount; // keep count of accumulated segments in order to accumulate further signals

    // cumulative prob distribution
    mutable bool cumulIsComputed;
    mutable vector<double> cumul;
    mutable bool cumulInverseIsComputed;
    mutable map<double,int> cumulInverse;

    void SetRange(const vector<Value> &samples);

public:
    static int ScottsRule(const vector<float>& vals); // for finding number of bins automatically

    vector<double> GetBins() const { return bins;}  
    int GetNBins() const {return bins.size();}
    void GetRange(float& rmin, float& rmax) const {rmin=vmin;rmax=vmax;}
    double GetMax() const {return *max_element(bins.begin(), bins.end());}
    double GetMaxPos() const {return BinValue(max_element(bins.begin(), bins.end())-bins.begin());}
    double operator[](int bin) const {return bins[bin];}
    double Density(int bin) const {return bins[bin]*bins.size()/(vmax-vmin);}
    void SetBins(float _vmin, float _vmax, int nbBins); 
    double BinValue(int binNum) const 
    {
        return vmin+binNum*(vmax-vmin)/bins.size();
    }

    double BinPos(Value v) const 
    {
        return (bins.size()*(v-vmin)/(vmax-vmin));
    }

    int BinNum(Value v) const 
    {
        // should use rint!
        int res=(int)(BinPos(v));
        if(res>=(int)bins.size() || res<0){ThrowError("ProbabilityDistribution::BinNum out of range:%f-> %d (min=%f max=%f)\n",v,res,vmin,vmax);}
        return res;
    }
protected:
    void ComputeCumul() const 
    {
        cumul.clear();
        cumul.assign(bins.size(),0);
        double v=0;
        for(uint i=0;i<bins.size();i++)
        {
            v+=bins[i];
            cumul[i]=v;
        }
        cumulIsComputed=true;
    }
    void ComputeInverseCumul() const
    {       
        cumulInverse.clear();
        if(!cumulIsComputed){ComputeCumul();}
        for(uint i=0;i<cumul.size();i++)
        {
            cumulInverse[cumul[i]]=i;
        }
        cumulInverseIsComputed=true;
    }
public:
    double CumulInverse(double v) const 
    {
        if(!cumulInverseIsComputed){ComputeInverseCumul();}
        map<double,int>::iterator res=cumulInverse.upper_bound(v);
        if(res==cumulInverse.end()){return vmax;}
        int cbin=(*res).second-1;
        return cbin*(vmax-vmin)/bins.size()+vmin;
    }
    double Cumul(double v) const 
    {
        if (!cumulIsComputed){ComputeCumul();}
        if (v <= vmin) return 0.0;
        if (v >= vmax) return 1.0;
        return cumul[BinNum(v)] - (BinValue(BinNum(v)+1) - BinValue(BinNum(v)))/(vmax-vmin)*bins.size()*bins[BinNum(v)];
    }

    bool NextIndex(vector<int> &indexes) 
    {
        for(uint i=0;i<indexes.size();i++)
        {
            indexes[i]++;
            if(indexes[i]<(int)bins.size()){return true;}
            indexes[i]=0;
        }
        return false;
    }

    class OtsuClass
    {
    public:
        int id;
        double moment0;
        double moment1;
        double moment2 ;
        int j0;
        int j1;
        OtsuClass():
            id(-1),moment0(0),moment1(0),moment2(0),j0(0),j1(0)
        {;}
        
    };

    OtsuClass ComputeClassStats(vector<OtsuClass> &classes,int classNum,vector<int> &indexes);
    void OtsuThreshholds(int nbClasses,vector<double> &res);

    void SetWithSamples(const vector<Value> &samples,int nbBins=1000);
    void AccumulateWithSamples(const vector<Value> &samples);

    void SetWithSignal(const vector<Value> &samples,int nbBins=1000);
    void AccumulateWithSignal(const vector<Value> &samples);

    void AccumulateImage3Dlinear3(const Image3Df &src);
    void SetWithImage3Dlinear3(const Image3Df &src,int nbBins);

    enum SetWithImage3DType {LINEAR3,LINEAR1,SAMPLES,FAST=LINEAR1,FASTEST=SAMPLES};
    void SetWithImage3D(const Image3Df &src,SetWithImage3DType type=LINEAR3,int nbBins=1000);

    void AccumulateWithImage3D(const Image3Df &src,SetWithImage3DType type=LINEAR3);

    void ShowWithGnuPlot(const string &label="") const ;
    
    ProbabilityDistribution():accumulatedSegmentCount(0),cumulIsComputed(false),cumulInverseIsComputed(false)
    {}
};

#endif // _ImageStats_hpp

---