大津法(otsu),中间像素统计法，kittler全局阈值图像二值化

图像的二值化有很多方法，这里介绍的三种是全局阈值的二值化方法。

这里给个原图：

1.大津法

最大类间方差法是由日本学者大津展之于1979年提出的，是一种自适应的阈值确定的方法,又叫大津法，简称OTSU。它是按图像的灰度特性,将图像分成背景和目标两部分。背景和目标之间的类间方差越大，说明构成图像的两部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。（详细算法）

大津算法的源代码非本人所写，见下面参考文献，这个算法实现的很好。大津算法的所用的时间是2.0ms。阈值为117.处理结果如下：

2.中间像素统计法

这个算法的具体名称我也不清楚，姑且这么叫吧。取像素数中值的灰度作为阈值，以此阈值进行二值化。

这个算法所用的时间1.3ms。阈值为130.处理结果如下：

3.kittler算法

是一种快速的全局阈值法。它的中心思想是，计算整幅图像的梯度灰度的平均值，以此平均值做为阈值。百科上说它比大津法速度快很多倍，而且效果也差不多，但是实际情况却并不理想，难道是我程序写的有问题？

这个算法运行时间有3.5ms，所确定的阈值是122.  处理结果如下：

源代码如下：

// 20130624_Binaryzation.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include  "highgui.h"
#include<windows.h>

//为了测试程序运行时间所定义的变量
LARGE_INTEGER freq, start, end;
char info[256];
double dTime = 0;

void midValueBinary( IplImage *imge)
{
	int *histogram = new int [256 ];
	int height = imge->height;
	int width = imge->width;
	int widthStep = imge->widthStep;
	char *pData = imge->imageData;

	int t_nTotalPiexl = imge -> height * imge -> width;

	int count = 0;

	for ( int i = 0; i < 256; ++i )
		histogram[i] = 0;

	//统计直方图中的相应的灰度的数量
	for ( int i = 0; i < height; ++i )
		for( int j = 0; j < width; ++j )
		{
			int pos = i * ( widthStep) + j;
			++histogram[ ( uchar ) pData[ pos ] ];
		}

		int threshold = 0;  //阈值大小
		for ( ; threshold < 256; ++threshold )
		{
			if ( count > t_nTotalPiexl / 2 )  //比总像素数一半大的灰度即为阈值
				break;
			count += histogram[threshold];
		}

		printf( " threshold的大小： %d ", threshold );

		for ( int i = 0; i < height ; ++i )
			for( int j = 0; j <width; ++j )
			{
				int curPos = i *widthStep + j;
				pData[curPos] = (uchar)pData[curPos] > threshold ? 255 : 0;
			}
			
			delete histogram;
}

void Kittler( IplImage *imge )
{
	int sumGrads = 0;
	int sumGrayGrads = 0;
	char *pData = imge->imageData;
	int height = imge ->height;
	int width = imge->width;
	int widthStep = imge->widthStep;
 
	for ( int i = 1; i < height - 1; i++ )
		for( int j = 1; j < width - 1; j++ )
		{
			int curPos = i * widthStep + j;
			int leftPos = curPos - 1; 
			int rightPos = curPos + 1;
            int upPos = curPos -widthStep;
			int downPos =curPos + widthStep;
			//求水平或垂直方向的最大梯度
			int Grads=MAX(abs(pData[leftPos]-pData[rightPos]),abs(pData[upPos] -pData[downPos]));
			sumGrads += Grads;
			//梯度与当前点灰度的积
			sumGrayGrads += Grads*((uchar)pData[curPos]);
		
	}
   int threshold = sumGrayGrads / sumGrads;
	printf("%d\n",threshold);

	for ( int i = 0; i < height; ++i )
		for( int j = 0; j <width; ++j )
		{
			int curPos = i *widthStep + j;
			pData[curPos] = (uchar)pData[curPos] >threshold ? 255 : 0;
		}
}
void myOtsu(const IplImage *frame) //大津法求阈值
{
    #define GrayScale 256	//frame灰度级
	int width = frame->width;
	int height = frame->height;
	int pixelCount[GrayScale]={0};
	float pixelPro[GrayScale]={0};
	int i, j, pixelSum = width * height, threshold = 0;
	uchar* pData = (uchar*)frame->imageData;

	//统计每个灰度级中像素的个数
	for(i = 0; i < height; i++)
	{
		for(j = 0;j < width;j++)
		{
			pixelCount[(int)pData[i * width + j]]++;
		}
	}

	//计算每个灰度级的像素数目占整幅图像的比例
	for(i = 0; i < GrayScale; i++)
	{
		pixelPro[i] = (float)pixelCount[i] / pixelSum;
	}

	//遍历灰度级[0,255],寻找合适的threshold
	float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0;
	for(i = 0; i < GrayScale; i++)
	{
		w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0;
		for(j = 0; j < GrayScale; j++)
		{
			if(j <= i)   //背景部分
			{
				w0 += pixelPro[j];
				u0tmp += j * pixelPro[j];
			}
			else   //前景部分
			{
				w1 += pixelPro[j];
				u1tmp += j * pixelPro[j];
			}
		}
		u0 = u0tmp / w0;
		u1 = u1tmp / w1;
		deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)) ;
		if(deltaTmp > deltaMax)
		{
			deltaMax = deltaTmp;
			threshold = i;
		}
	}
	//return threshold;
	printf("阈值为：%d \n ", threshold);
	for ( int i = 0; i < height ; ++i )
		for( int j = 0; j <width; ++j )
		{
			int curPos = i *frame ->widthStep + j;
			pData[curPos] = (uchar)pData[curPos] >threshold ? 255 : 0;
		}
}
int main(int argc, char* argv[])
{
	IplImage *imge  = cvLoadImage("lena.bmp", -1);

	//测量程序运行时间开始
	QueryPerformanceFrequency(&freq);
	QueryPerformanceCounter(&start);

	//要运行的函数
//midValueBinary(imge);
//Kittler(imge);
myOtsu(imge);

	//测量程序运行时间结束
	QueryPerformanceCounter(&end);
	dTime = (double)(end.QuadPart - start.QuadPart);
	dTime = 1000*dTime / freq.QuadPart;
	printf( " time: %.2f ms.\n", dTime);


	cvShowImage("lena.bmp", imge);

			 // wait for a key
	cvWaitKey(0);

			 // release the image
	cvReleaseImage(&imge);
	
	return 0;
}

参考文献：

http://blog.csdn.net/daxiamit/article/details/6619075

http://baike.baidu.cn/view/1532602.htm

http://zh.wikipedia.org/wiki/%E5%A4%A7%E6%B4%A5%E7%AE%97%E6%B3%95