Ostu算法

        图像二值化算法是图像处理的基础。一般来说，二值化算法可以分为两个类别：全局二值化和局部二值化。全局二值化是指通过某种算法找到一个全局的阈值T，对图像中坐标为（x,y）的像素值做如下处理：f(x,y)>T则f(x,y)=255,else f(x,y)=0
        Ostu就是这样一种全局二值化算法，又叫最大类间方差。因为该算法会遍历图像中任意一个像素值i，计算当其为阈值时，图像的前景和背景图像（并不一定是真正的前景和背景,只是我们把当前小于i的记做背景，大于i的记做前景）的方差值。当方差值达到最大时，我们认为此时的i是该图像的全局阈值。
说明： 
前景点数占图像比例：w0；平均灰度：u0
背景点数占图像比例：w1，平局灰度：u1
图像平均灰度：u=w0*u0 + w1*u1
前景和背景图像方差：

g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)

  =wo * w1 * (uo - u1) * (uo - u1)  

int otsu(Mat dst){

	int i, j;
	int tmp;

	double u0, u1, w0, w1, u, uk;

	double cov;
	double maxcov = 0.0;
	int maxthread = 0;

	int hst[MAX_GRAY_VALUE] = { 0 };
	double pro_hst[MAX_GRAY_VALUE] = { 0.0 };

	int height = dst.cols;
	int width = dst.rows;

	//统计每个灰度的数量
	for (i = 0; i<width; i++){
		for (j = 0; j<height; j++){
			tmp = dst.at<uchar>(i, j);
			hst[tmp]++;
		}
	}

	//计算每个灰度级占图像中的概率
	for (i = MIN_GRAY_VALUE; i<MAX_GRAY_VALUE; i++)
		pro_hst[i] = (double)hst[i] / (double)(width*height);

	//计算全局平均灰度值
	u = 0.0;
	for (i = MIN_GRAY_VALUE; i<MAX_GRAY_VALUE; i++)
		u += i*pro_hst[i];

	//double det = 0.0;
	//for (i = MIN_GRAY_VALUE; i< MAX_GRAY_VALUE; i++)
		//det += (i - u)*(i - u)*pro_hst[i];

	//统计前景和背景的平均灰度值，并计算类间方差

	for (i = MIN_GRAY_VALUE; i<MAX_GRAY_VALUE; i++){

		w0 = 0.0; w1 = 0.0; u0 = 0.0; u1 = 0.0; uk = 0.0;

		for (j = MIN_GRAY_VALUE; j < i; j++){

			uk += j*pro_hst[j];
			w0 += pro_hst[j];//前景占图像画幅比例

		}
		u0 = uk / w0;//前景平均灰度
		w1 = 1 - w0;//背景占图像画幅比例
		u1 = (u - uk) / (1 - w0);//背景平均灰度

		//计算类间方差
		cov = w0*w1*(u1 - u0)*(u1 - u0);

		if (cov > maxcov)
		{
			maxcov = cov;
			maxthread = i;
		}
	}

	cout << maxthread << endl;
	return maxthread;
}


int main(){
	int width, height;
	int i, j;
	Mat obj = imread("E:/VS2013/face/xuelian/png/1.png");

	Mat dst;
	cvtColor(obj, dst, CV_RGB2GRAY);//灰度化
	
	height = dst.cols;//计算图像高度
	width = dst.rows;//计算图像宽度
	int thd = otsu(dst);

	imshow("origin img", dst);

	for (i = 0; i < width; i++)
	for (j = 0; j< height; j++)
	if (dst.at<uchar>(i, j) > thd)
		dst.at<uchar>(i, j) = MAX_GRAY_VALUE-1;
	else
		dst.at<uchar>(i, j) = MIN_GRAY_VALUE;
	imshow("ostued", dst);
	waitKey(0);
	return 0;
}