ImageStats.cpp

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.
 */
#include<ImLib3D/ImageStats.hpp>
#include<ImLib3D/TestPatterns.hpp>
#include<ImLib3D/Arithmetic.hpp>
#include<ImLib3D/ConvenienceProcessors.hpp>

namespace IP3D
{
void 
RobustAverageAndVariance(const Image3Df &src,double &average,double &variance)
{
    // compute median as robust estimator of average
    IP3D::Median(src,average);
    // compute median absolute deviation (MAD) as estimator of std dev
    // MAD =1.4826 * median |xi- avg|
    Image3Df diff=src;
    diff-=average;
    IP3D::Abs(diff,diff);
    IP3D::Median(diff,variance);
    variance=1.4826*variance;
    // variance = stddev^2
    variance=variance*variance;
}
}

void 
ProbabilityDistribution::SetRange(const vector<Value> &samples)
{
    vmin=*min_element(samples.begin(),samples.end());
    vmax=*max_element(samples.begin(),samples.end());
    double grow=(vmax-vmin)*.0001;
    if(vmax==vmin){mprintf("ProbabilityDistribution::SetRange: WARNING vmax==vmin==%f\n",vmax);grow=.0001;}
    vmin-=grow;
    vmax+=grow;
}


// compute statistics for this class
// This just computes the first three moments of the distribution for pixels in class classNum
// Hiwever, it keeps memory of previously computed class stats so it doen't need to
// recompute everything each time
ProbabilityDistribution::OtsuClass
ProbabilityDistribution::ComputeClassStats(vector<OtsuClass> &classes,int classNum,vector<int> &indexes)
{
    double moment1=0,moment2=0,moment0=0;
    int j0=(classNum>0                   ? indexes[classNum-1] : 0               );
    int j1=(classNum<(int)indexes.size() ? indexes[classNum  ] : (int)bins.size());

//      if(dotalk)mprintf("class %d: from %f to %f\n",
//                       classNum,vmin+j0*((vmax-vmin)/bins.size()),vmin+j1*((vmax-vmin)/bins.size())   );

    OtsuClass &old=classes[classNum];
    // try using last statistics to avoid recomputing 
    if(old.id==classNum && old.j0==j0   && old.j1==j1)
    {
        moment0=old.moment0;
        moment1=old.moment1;
        moment2=old.moment2;                    
    }
    else
    if(old.id==classNum && old.j0==j0-1 && old.j1==j1)
    {
        int j=j0-1;
        double v=vmin+j*((vmax-vmin)/bins.size());
        moment0=old.moment0 -     bins[j];
        moment1=old.moment1 - v*  bins[j];
        moment2=old.moment2 - v*v*bins[j];                  
    }
    else
    if(old.id==classNum && old.j0==j0   && old.j1==j1-1)
    {
        int j=j1-1;
        double v=vmin+j*((vmax-vmin)/bins.size());
        moment0=old.moment0 +     bins[j];
        moment1=old.moment1 + v*  bins[j];
        moment2=old.moment2 + v*v*bins[j];                  
    }
    else
    {
        // *** we must recompute!****

//                  mprintf("recomp for (class %d):%d %d (old %d %d)  : indexes:",i,j0,j1,old.j0,old.j1);
//                  for(int zz=0;zz<indexes.size();zz++){mprintf("%d:%d:%f  ",zz,indexes[zz],thresholds[zz]);}
//                  mprintf("\n");

        for(int j=j0;j<j1;j++)
        {
            double v=vmin+j*((vmax-vmin)/bins.size());
            moment0+=    bins[j];
            moment1+=v*  bins[j];
            moment2+=v*v*bins[j];
        }
    }
    

    // save stats for avoiding recomputing 
    old.id=classNum;
    old.j0=j0;
    old.j1=j1;
    old.moment0=moment0;
    old.moment1=moment1;
    old.moment2=moment2;
    OtsuClass res;
    res.id=classNum;
    
    res.id=classNum;
    res.moment0=moment0;
    res.moment1=moment1;
    res.moment2=moment2;
    res.j0=j0;
    res.j1=j1;
    return res;
}


void 
ProbabilityDistribution::OtsuThreshholds(int nbClasses,vector<double> &res)
{
    uint i;
    res.assign(nbClasses-1,0);
    vector<double> thresholds;
    vector<int> indexes;
    indexes.assign(nbClasses-1,0);
    thresholds.assign(nbClasses-1,0);
    double minS2=-1;
    Periodic talk(.5);
    vector<OtsuClass> classes;
    OtsuClass dud;
    classes.assign(nbClasses,dud);

    // try out all possible threshold combinations
    while(NextIndex(indexes))
    {
        // check if thresholds are in order t1<t2<t3...
        for(i=1;i<indexes.size();i++){if(indexes[i]<=indexes[i-1]) break;}
        if(i<indexes.size()){continue;}// bad order, try again!

        int dotalk=talk();
        // compute thresholds from indexes
        for(i=0;i<indexes.size();i++){thresholds[i]=vmin+indexes[i]*((vmax-vmin)/bins.size());}
        if(dotalk)
        {
            mprintf("considering thresholds:");
            for(i=0;i<indexes.size();i++){mprintf("%d:%f  ",i,thresholds[i]);}
            mprintf("\n");
        }

        // now compute the within - class variance using theese thresholds
        double s2=0;
        for(int classNum=0;classNum<nbClasses;classNum++)
        {
            // compute statistics for this class
            OtsuClass c=ComputeClassStats(classes,classNum,indexes);
            // empty class... ignore this segmentation
            if(c.moment0==0){s2=2*minS2+1;break;}

            double classProb=c.moment0;
            double average  =c.moment1/classProb;
            // variance within this class
            double variance =c.moment2/classProb-average*average;

            // compute s2
            s2+=variance*classProb;
            if(dotalk)mprintf("class %d: from %f to %f: classProb:%f average:%f variance:%f contr:%f\n",
                             classNum,vmin+c.j0*((vmax-vmin)/bins.size()),vmin+c.j1*((vmax-vmin)/bins.size()),
                             classProb,average,variance,variance*classProb);
        }
        // check if this threshold combination is optimal
        if(s2<minS2 || minS2<0)
        {// it is!
            minS2=s2;
            res=thresholds;
        }
        if(dotalk)
        {
            mprintf("TOT S2:%f BEST:%f  : ",s2,minS2);
            for(uint i=0;i<res.size();i++){mprintf("%d:%f  ",i,res[i]);}
            mprintf("\n");
        }
    }
    mprintf("OtsuThresholds result:");
    for(i=0;i<res.size();i++){mprintf("%d:%f  ",i,res[i]);}
    mprintf("\n");
}
void 
ProbabilityDistribution::ShowWithGnuPlot(const string &label) const
{
    vector<double> values(bins.size());
    for(uint i=0;i<bins.size();i++)
    {
        values[i]=(i*(vmax-vmin)/bins.size())+vmin;
    }
    GnuPlot(values,bins,label,"histogram");
}

void 
ProbabilityDistribution::SetWithSamples(const vector<Value> &samples,int nbBins)
{
    bins.assign(nbBins,0);
    SetRange(samples);
    accumulatedSegmentCount = 0;

    AccumulateWithSamples(samples);
}

void 
ProbabilityDistribution::AccumulateWithSamples(const vector<Value> &samples)
{
    // We first make histogram, then accumulate before we again normalize to make probability
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]*=accumulatedSegmentCount;
    }

    for(uint i=0;i<samples.size();i++)
    {
        bins[BinNum(samples[i])]+=1.0;
        accumulatedSegmentCount++;
    }
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]/=accumulatedSegmentCount;
    }

//      double totproba=0;
//      for(int tt=0;tt<bins.size();tt++)
//      {
//          totproba+=bins[tt];
//      }
//  mprintf("totproba:%f\n",totproba);
}

void 
ProbabilityDistribution::SetWithSignal(const vector<Value> &samples,int nbBins)
{
    bins.assign(nbBins,0);
    SetRange(samples);
    accumulatedSegmentCount = 0;

    AccumulateWithSignal(samples);
}

void 
ProbabilityDistribution::AccumulateWithSignal(const vector<Value> &samples)
{
    // We first make histogram, then accumulate before we again normalize to make probability
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]*=accumulatedSegmentCount;
    }

    // for each sample pair fill all the bins whose values span both samples
    for(uint i=0;i<samples.size()-1;i++)
    {
        Value v0=samples[i];
        Value v1=samples[i+1];
        Value t=max(v0,v1);
        v0=min(v0,v1);
        v1=t;

        if( BinNum(v0)==BinNum(v1))
        {// just one bin:
            bins[BinNum(v0)]+=1.0;
        }
        else
        {
            // more than one bin:
            double j0=BinPos(v0);
            double j1=BinPos(v1);
            //          double dprob=1.0/(samples.size()*(j1-j0));
            double dweight=1.0/(j1-j0);
            int cj0=(int)ceil(j0);
            int fj1=(int)floor(j1);
                
            double dy;
            // first bin
            dy=cj0-j0;
            bins[(int)j0]+=dy*dweight;
            // middle bins
            for(int j=cj0;j<fj1;j++)
            {
                bins[j]+=dweight;
            }
            // last bin
            dy=j1-fj1;
            bins[(int)j1]+=dy*dweight;
        }
        accumulatedSegmentCount++;
    }
    // renormalize histogram to obtain proba
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]/=accumulatedSegmentCount;
    }
    // display totproba
    double totproba=0;
    for(uint tt=0;tt<bins.size();tt++)
    {
        totproba+=bins[tt];
    }
    mprintf("totproba:%f bins:%d\n",totproba,bins.size());
}

void 
ProbabilityDistribution::SetWithImage3D(const Image3Df &image,SetWithImage3DType type,int nbBins)
{
    if(type==LINEAR3){SetWithImage3Dlinear3(image,nbBins);return;}
    vector<float> samples;
    samples.reserve(image.GetNVoxels());
    for(Image3Df::const_iteratorFastMasked p=image.begin();p!=image.end();++p)
    {
        samples.push_back(*p);
    }
    if(type==LINEAR1){SetWithSignal(samples,nbBins);return;}
    if(type==SAMPLES){SetWithSamples(samples,nbBins);return;}
    ThrowError("ProbabilityDistribution::SetWithImage3D unknown type:%d\n",type);
}

void 
ProbabilityDistribution::AccumulateWithImage3D(const Image3Df &image,SetWithImage3DType type)
{
    if(type==LINEAR3){AccumulateImage3Dlinear3(image);return;}
    vector<float> samples;
    samples.reserve(image.GetNVoxels());
    for(Image3Df::const_iteratorFastMasked p=image.begin();p!=image.end();++p)
    {
        samples.push_back(*p);
    }
    if(type==LINEAR1){AccumulateWithSignal(samples);return;}
    if(type==SAMPLES){AccumulateWithSamples(samples);return;}
    ThrowError("ProbabilityDistribution::SetWithImage3D unknown type:%d\n",type);
}


void 
ProbabilityDistribution::SetWithImage3Dlinear3(const Image3Df &src,int nbBins)
{
    Vect3Di srcSize=src.SizeV();
    IP3D::FindMinMax(src,vmin,vmax);
    double grow=(vmax-vmin)*.0001;
    if(vmax==vmin){mprintf("ProbabilityDistribution::SetWithImage3D: WARNING vmax==vmin==%f\n",vmax);grow=.0001;}
    vmin-=grow;
    vmax+=grow;
    bins.assign(nbBins,0);
    accumulatedSegmentCount = 0;

    AccumulateImage3Dlinear3(src);
}

void 
ProbabilityDistribution::AccumulateImage3Dlinear3(const Image3Df &src)
{
    if (bins.size() == 0)
        ThrowError("Histogram must be set before accumulating");
    
    Vect3Di srcSize=src.SizeV();
    int nbBins = bins.size();

    // We first make histogram, then accumulate before we again normalize to make probability
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]*=accumulatedSegmentCount;
    }

    int segmentCount=0;
    for(int dir=0;dir<3;dir++)
    {
        int d0=(dir+1)%3;
        int d1=(dir+2)%3;
        int d2=dir;
        float minDSrc=.001*(vmax-vmin)/(double)nbBins;
        Vect3Di dpos(0,0,0);dpos(d2)=1;
        Vect3Di pos(0,0,0);
        for(int x0=0;x0<srcSize(d0);x0++)
        {
            pos(d0)=x0;
            for(int x1=0;x1<srcSize(d1);x1++)
            {
                pos(d1)=x1;
                for(int x2=0;x2<srcSize(d2)-1;x2++)
                {
                    pos(d2)=x2;
                    Vect3Di pos0=pos;
                    Vect3Di pos1=pos+dpos;
                    // skip if we're on unmasked
                    if(src.properties.mask && !src.Mask()(pos0)){continue;}
                    if(src.properties.mask && !src.Mask()(pos1)){continue;}
                    segmentCount++;
                    float src0=src(pos0);
                    float src1=src(pos1);
//                  pos0.Show("pos0:","   ");
//                  pos1.Show("pos1:","\n");
//                  printf("src0:%f     src1:%f\n",src0,src1);
                                                    
                    // swap so that src1 > src0
                    if(src1<src0){swap(src0,src1);}
                    if(src1-src0<minDSrc){src1=src0+minDSrc;}// make sure dsrc!=0
                    int srcBin0=BinNum(src0);
                    int srcBin1=BinNum(src1);
//                  printf("src0:%f     src1:%f srcBin0:%d srcBin1:%d\n",src0,src1,srcBin0,srcBin1);
                    for(int s=srcBin0;s<=srcBin1;s++)
                    {
//                      float srcvb0=max(src0,(float)BinValue(s));
//                      float srcvb1=min(src1,(float)BinValue(s+1));
//                      float srcBinProba=(srcvb1-srcvb0)/(src1-src0);
//                      float srcBinProba=srcBin1-srcBin0;
                        bins[s]+=1.0/(1+srcBin1-srcBin0);
                    }
                }
            }
        }
    }
    accumulatedSegmentCount += segmentCount;
    for(uint tt=0;tt<bins.size();tt++)
    {
        bins[tt]/=accumulatedSegmentCount;
    }
    // display totproba
    double totproba=0;
    for(uint tt=0;tt<bins.size();tt++)
    {
        totproba+=bins[tt];
    }
    mprintf("totproba2:%f bins:%d\n",totproba,bins.size());

}


int 
ProbabilityDistribution::ScottsRule(const vector<float>& vals)
{
    float N=vals.size();
    double av=0;
    double var=0;
    for (int i=0;i<N;i++)
    {
        av+=vals[i]; 
        var+=vals[i]*vals[i];
    }
    av/=N;
    var/=(N-1);
    var-=N/(N-1)*av*av;

    double s=sqrt(var);
    double w=3.5*s/pow((double)N, 1/3.0);
    float maximum = *max_element(vals.begin(), vals.end());
    float minimum = *min_element(vals.begin(), vals.end());
    int nBins=(int)((maximum-minimum)/w);
    return nBins;
}

void 
ProbabilityDistribution::SetBins(float _vmin, float _vmax, int nbBins)
{
    if (bins.size() > 0) 
        bins.clear();
    vmin=_vmin; 
    vmax=_vmax;
    bins.assign(nbBins,0);
    accumulatedSegmentCount = 0; 
}

---