/* 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.
 */
/*! \file   Threshold.hpp */

#ifndef _Threshold_hpp
#define _Threshold_hpp

#include<ImLib3D/ImageProcessor.hpp>



//! Simple Threshold 
/*! Result=1 if value>threshold  
 */
namespace IP3D
{
/*! #&SimpleThreshold&#     Simple Threshold.
        Result=1 if value>threshold. Result is a Mask3D.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Threshold     type:Im3DValue
      <LI>Result mask   type:Mask3D
     </UL>
~*/
	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;}
	}
};



//! Simple Threshold 
/*! Result=1 if t1<=value<=t2
 */
namespace IP3D
{
/*! #&SimpleThresholds2&#     Simple Threshold between 2 values 
        Result is 1 for values between t1 and t2, 0 otherwise.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Threshold 1   type:Im3DValue
      <LI>Threshold 2   type:Im3DValue
      <LI>Result mask   type:Mask3D
     </UL>
~*/
	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;}
	}

};

//! Limit 
/*! Result=t1 if value<=t1    result=t2 t2<=value
 */
namespace IP3D
{
/*! #&LimitThreshold&#     Limits range of image values to t1 t2 
        If image value less than t1 or image value greater 
        than t2 then image value equals t1 or t2 respectively, otherwise 
        conserves original value. Keeps mask
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Lower limit   type:Im3DValue
      <LI>Upper limit   type:Im3DValue
      <LI>Result image  type:ImageType
     </UL>
~*/
	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);}
		}
	}

};


//! Otsu Thresholding
/*! Result=labled image (at most 255 classes)
 */
namespace IP3D
{
/*! #&OtsuThresholds&#     Use otsu histogram partitioning to automatically determine thresholds and apply 
    them on an image.  
        The algorithm brute-force searches all possible threshhold combinations to 
        find the one that minimizes interclass variance. This can be long if the 
        number of classes is high.
     <UL>Args:
      <LI>Source image                 type:ImageType
      <LI>Result image                 type:Mask3D
      <LI>Number of classes            type:int
      <LI>Result thresholds            type:vector<double>
      <LI>Number of bins in histogram  type:int
     </UL>
~*/
	template <class ImageType>
	void OtsuThresholds(const ImageType& src, Mask3D& res,int nbClasses=2,vector<double> *pthresholds=NULL,int nbBins=1000);
};


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




#ifdef DEPRECATED

//! Simple optimization,Given a function f, finds x / f(x)=y
class Optimizer1D
{
	Function1Dd &f;
	double precision;
public:
	//! simple works for monotonous functions
	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
