Threshold.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/Image3D.hpp>
#include<ImLib3D/Threshold.hpp>
//  #include<ImLib3D/Display.hpp>
//  #include<ImLib3D/TestPatterns.hpp>
//  #include<ImLib3D/ConvenienceProcessors.hpp>

void
IP3D::UniModalThreshold(const Image3Df& src, double& threshold, Mask3D* res)
{
    int nbBins=500;

    vector<float> values;
    values.reserve(src.GetNVoxels());
    Image3Df::const_iteratorFastMasked p;
    for(p=src.begin();p!=src.end();++p){values.push_back(*p);}
    
    ProbabilityDistribution dist;
    dist.SetWithSamples(values,nbBins);

//  cout << "pos of max = " << max_element(bins.begin(), bins.end())-bins.begin() << endl;
    float max=dist.GetMax();
    float a0=dist.GetMaxPos();
    float a1, dummy;

    dist.GetRange(dummy, a1);
//      cout << "Histogram range ["<<dummy<<":"<<a1<<"]"<<endl;
//      cout <<"max=" << max << " at pos=" << a0 << endl;

    // y=b0+b1x is the slope for max histo value to max image value
    float b0, b1;
    b1=max/(a0-a1);
    b0=max-b1*a0;

//  iterator it  = max_element(bins.begin(), bins.end());
//      vector<float> datas;
//      copy(bins.begin(), bins.end(), back_inserter(datas));
//      Plot plot;
//      plot.SetStyle("lines");
//      plot.PlotData(datas,dummy, a1);
//      plot.PlotSlope(b1, b0);

    float maxdist=0;
    float maxdistindex=0;
    int startindex=(int)dist.BinNum(a0);
//  cout << "Start search index = " << startindex << endl;

    // find biggest distance of histogram perpendicular on the slope (b0,b1)
    for (int i=startindex; i<dist.GetNBins();i++)
    {
        float distance = (b1*dist.BinValue(i)+b0-dist[i]);
        if (distance>maxdist)
        {
            maxdist=distance;
            maxdistindex=i;
        }
    }

//  cout << "Found Threshold : " << dist.BinValue(maxdistindex) << endl;
    threshold=dist.BinValue((int)maxdistindex);
    if (res)
        IP3D::SimpleThreshold(src, (float)threshold, *res);
}

---