Threshold.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 _Threshold_hpp
#define _Threshold_hpp

#include<ImLib3D/ImageProcessor.hpp>




namespace IP3D
{
    template <class ImageType>
    void SimpleThreshold(const ImageType& src,const typename ImageType::value_type &t,Mask3D& res)
    {
        res.Resize(src.Size());
        typename ImageType::const_iteratorFast sp;
        Mask3D::iteratorFast rp;
        for(sp=src.begin(),rp=res.begin();sp!=src.end();++sp,rp++){*rp= *sp>t ? 1 : 0;}
    }
};




namespace IP3D
{
    template <class ImageType>
    void SimpleThresholds2(const ImageType& src, 
                      const typename ImageType::value_type &t1,
                      const typename ImageType::value_type &t2,
                      Mask3D& res)
    {
        res.Resize(src.Size());
        typename ImageType::const_iteratorFast sp;
        Mask3D::iteratorFast rp;
        for(sp=src.begin(),rp=res.begin();sp!=src.end();++sp,rp++){*rp= (*sp>=t1 && *sp<=t2) ? 1 : 0;}
    }

};


namespace IP3D
{
    template <class ImageType>
    void LimitThreshold(const ImageType& src, 
                      const typename ImageType::value_type &t1,
                      const typename ImageType::value_type &t2,
                      ImageType& res)
    {
        res.Resize(src.Size());
        if (src.HasMask())
        {
            res.AddMask()=src.Mask();
            typename ImageType::const_iteratorFastMasked sp;
            typename ImageType::iteratorFastMasked rp;
            for(sp=src.begin(),rp=res.begin();sp!=src.end();++sp,rp++){*rp= min(max(*sp,t1),t2);}
        }
        else
        {
            typename ImageType::const_iteratorFast sp;
            typename ImageType::iteratorFast rp;
            for(sp=src.begin(),rp=res.begin();sp!=src.end();++sp,rp++){*rp= min(max(*sp,t1),t2);}
        }
    }

};



namespace IP3D
{
    template <class ImageType>
    void OtsuThresholds(const ImageType& src, Mask3D& res,int nbClasses=2,vector<double> *pthresholds=NULL,int nbBins=1000);
};


namespace IP3D
{
        void UniModalThreshold(const Image3Df& src, double& threshold, Mask3D* res=NULL);
};




#ifdef DEPRECATED

class Optimizer1D
{
    Function1Dd &f;
    double precision;
public:
    double Solve(double x0,double ygoal=0,double x1=0,double x2=0)
    {
        if(x2==0 && x1==0){x2=x0+precision;x1=x0-precision;}
        mprintf("starting interval %f - %f    : x0=%f\n",x1,x2,x0);

        double y1=f(x1);
        double y2=f(x2);
        mprintf("first pos: %f %f  - %f %f\n",x1,y1,x2,y2);
        for(int i=0;y2==y1 && i<100;i++)
        {
            x2=x1+(x2-x1)*2;
            x1=x2+(x1-x2)*2;
            y1=f(x1);
            y2=f(x2);
            mprintf("y1==y2!expanded to %f %f  - %f %f  \n",x1,y1,x2,y2);
        }
        if(y2==y1){ThrowError("Optimizer1D:Solve failed y2==y1");}
        while(fabs(x2-x1)>precision)
        {
            // new estimate of solution
            double x=x1+(ygoal-y1)*(x2-x1)/(y2-y1);
            double y=f(x);
            mprintf("on %f %f  - %f %f   new estimate:%f %f\n",x1,y1,x2,y2,x,y);
            // replace worst point by new estimate
            if(y==y1){x1=x;}
            else      
            if(y==y2){x2=x;}
            else
            if(fabs(x2-x)>fabs(x1-x) && y!=y1)
            {x2=x;y2=y;}
            else
            {x1=x;y1=y;}
            if(y1==y2){ThrowError("Optimizer1D:Solve failed y2==y1 in opt loop");}
        }       
        return (x2+x1)/2;
    }
    Optimizer1D(Function1Dd &_f,double _precision):
        f(_f),precision(_precision)
    {
    }
};

#endif// DEPRECATED



#include<ImLib3D/Threshold.hxx>

#endif // _Threshold_hpp

---