minimum error threshold 原理浅谈及c++实现

参考论文来自1986kittler，原理论证可以自行参考其论文，这里只关心其具体实现步骤。
看看如何从一副图像得到最小误差阈值法下的参数threshold。
首先对一幅图像，通过calcHist函数可以得到其直方图，图像的直方图则可以看做是不同像素值的点的概率密度分布。
总点数为image.rows * image.cols，对直方图中每一个柱条除以总点数得到其出现的先验概率。
显然，这些先验概率的概率之和为1。不妨将这些先验概率记作h(g) ,g为图像像素值取值为 [min_pixel,max_pixel]。
下一步则进行我们的遍历，取T在g[min_pixel,max_pixel]范围内的每一个像素值。那么T就把我们图像的像素点分为了两部分，一部分像素点像素值≤T记作A类，而另一部分像素点像素值大于T记作B类。
则A类像素出现的概率之和以及B类像素出现的概率之和分别为：
P1(T)=∑ h(g) g从min_pixel到T
P2(T)=∑ h(g) g从T+1到max_pixel
类似的我们要求出A、B类的平均值UA，UB以及各自的方差sigmaA，sigmaB。计算公式如下：
UA (T)= (1/ P1(T) )* (∑ g* h (g) ) g从min_pixel到T
UB (T)= (1/ P2(T) ) * (∑ g * h(g) ) g从T+1到max_pixel
sigmaA= (1/ P1(T) ) * {∑ [ (g-UA (T) )^2 ] * h(g)} g从min_pixel到T
sigmaB= (1/ P2(T) ) * {∑ [(g-UB(T))^2] * h(g)} g从T+1到max_pixel
恭喜，到这里我们已经得到了目标函数的所有参数，所谓阈值误差就是用一个目标函数去衡量，得分高误差就大，得分低误差就小。
通过遍历图像本身的所有阈值可能取值，找到对应最小目标函数取到的阈值，即为我们最小阈值误差法下得到的结果。目标函数如下：
J[T] = 1 + 2 * { P1 (T) * [log (sigmaA) - log ( P1 (T) )] + P2 (T) * [log (sigmaB) - log ( P2 (T) )]}
完成，现在只需要遍历J[T]找到T∈[min_pixel,max_pixel]的值返回输出即可。
好了，介绍完毕，接下来是具体的代码展示，
不需要的小伙伴可以自己动手实现一下。

 void minimum_error_threshold(cv::Mat &matInput,int* pixel); 
 //函数输入待处理图像mat矩阵，和指向需要返回整数像素pixel的指针。
 //最终返回值为得到的阈值

void  minimum_error_threshold(cv::Mat &matInput,int* pixel)
{
	Mat gray = matInput;
	int bins =256;
	int hist_size[]={bins};
	float range[]={0,256};
	const float * ranges[]={range};
	Mat   hist;
	int channels[]={0};
	calcHist(&gray,1,channels,Mat(),hist,1,hist_size,ranges,true,false);
	int min_pixel[1];
	int max_pixel[1];
	int hist_minpixel;
	int hist_maxpixel;
	for (hist_minpixel=0;hist_minpixel<256;hist_minpixel++)
	{
		if (hist.at<float>(hist_minpixel,0)!=0)
		{
			break;
		}
	}
	for (hist_maxpixel=255;hist_maxpixel>-1;hist_maxpixel--)
	{
		if (hist.at<float>(hist_maxpixel,0)!=0)
		{
			break;
		}
	}
	float a=1.0/512.0/512.0;
	hist=a*hist;
	min_pixel[0]=hist_minpixel;
	max_pixel[0]=hist_maxpixel;
	float J_t[256];
	int i;
	for (i=0;i<256;i++)
	{
		J_t[i]=0;
	}
	float p1_t=0,p2_t=0;
	float u1_t=0,u2_t=0;
	float sigma1=0,sigma2=0;
	int g=0;
	int j=0;
	for (i=min_pixel[0];i<max_pixel[0]+1;i++)
	{
		g=i;
		for (j=min_pixel[0];j<g+1;j++)
		{
			p1_t=p1_t+hist.at<float>(j,0);

		}
		for (j=min_pixel[0];j<g+1;j++)
		{
			u1_t=u1_t+((float)1.0/p1_t)*j*hist.at<float>(j,0);
		}
		for (j=min_pixel[0];j<g+1;j++)
		{
			sigma1=sigma1+((float)1.0/p1_t)*(j-u1_t)*(j-u1_t)*hist.at<float>(j,0);
		}
		for (j=g+1;j<max_pixel[0]+1;j++)
		{
			p2_t=p2_t+hist.at<float>(j,0);
		}
		for (j=g+1;j<max_pixel[0]+1;j++)
		{
			u2_t=u2_t+((float)1.0/p2_t)*j*hist.at<float>(j,0);
		}
		for (j=g+1;j<max_pixel[0]+1;j++)
		{
			sigma2=sigma2+((float)1.0/p2_t)*(j-u2_t)*(j-u2_t)*hist.at<float>(j,0);
		}
		J_t[i]=1+2*(p1_t*(log(sigma1)-log(p1_t))+p2_t*(log(sigma2)-log(p2_t)));
		if (J_t[i]<0)
		{
			J_t[i]=1000;
		}
		if (J_t[i]!=J_t[i])
		{
			J_t[i]=1001;
		}
		p1_t=0;
		p2_t=0;
		u1_t=0;
		u2_t=0;
		sigma1=0;
		sigma2=0;
	}
	float x=J_t[min_pixel[0]];
	int threshold=min_pixel[0];
	for (i=min_pixel[0];i<max_pixel[0]+1;i++)
	{
		if (x>J_t[i])
		{
			x=J_t[i];
			threshold=i;

		}
		cout<<J_t[i]<<endl;
	}
	*pixel=threshold;
}

得到阈值后直接调用OpenCV库里的threshold函数看看劳动后的成果如何吧！

收工：
力量似乎渐渐回来了。不多，但有用。