OTSU算法原理

OTSU算法原理及实现：

最大类间方差是由日本学者大津(Nobuyuki Otsu)于1979年提出，是一种自适应的阈值确定方法。算法假设图像像素能够根据阈值，被分成背景[background]和目标[objects]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大【用方差表达，具体公式见后】。OTSU的扩展算法，可进行多级阈值处理，称为“Multi Otsu method”【题外话】

设原始灰度级为M，灰度级为i的像素点个数为n_i，对灰度直方图进行归一化：

 

opencv实现代码：

int CLxrImage::Otsu(IplImage* src)  
{  
	int height=src->height;  
	int width=src->width;      

	//统计直方图  
	int   nHistogram[256] = {0};
	float fHistogram[256] = {0};  

	for(int i=0; i < height; i++)
	{  
		for(int j = 0; j < width; j++) 
		{  
			nHistogram[(unsigned char)src->imageData[i*src->widthStep+j]]++;  
		}  
	}

	//归一化直方图 
	int size = height * width;  
	for(int i = 0; i < 256; i++)
	{  
		fHistogram[i] = nHistogram[i] / (float)size;  
	}  

	//average pixel value  
	float avgValue=0;  
	for(int i=0; i < 256; i++)
	{  
		avgValue += i * fHistogram[i];  //整幅图像的平均灰度
	}   

	int threshold;    
	float maxVariance=0;  
	float w = 0, u = 0;  
	for(int i = 0; i < 256; i++) 
	{  
		w += fHistogram[i];			//假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例
		u += i * fHistogram[i];     // 灰度i 之前的像素(0~i)的平均灰度值： 前景像素的平均灰度值

		float t = avgValue * w - u;  
		float variance = t * t / (w * (1 - w) );  
		if(variance > maxVariance) 
		{  
			maxVariance = variance;  
			threshold = i;  
		}  
	}  
	return threshold;  
}
int CLxrImage::Otsu2(IplImage* src)
{
	int i,j,nThresh;
	int nHistogram[256] = {0};
	double fStdHistogram[256] = {0.0};
	double fGrayAccu[256] = {0.0};
	double fGrayAve[256] = {0.0};
	double fAverage = 0;
	double fTemp = 0;
	double fMax = 0;
	//统计直方图
	for(i = 0; i <src->height; i++)
	{
		for(j = 0; j <src->width; j++)
		{
			nHistogram[(unsigned char)src->imageData[i*src->width+j]] ++;
		}
	}
	//归一化直方图
	for(i = 0; i <= 255;i++)
	{
		fStdHistogram[i] = nHistogram[i]/(double)(src->width * src->height);
	}
	for(i=0;i<=255;i++)
	{
		for(j=0;j<=i;j++)
		{
			fGrayAccu[i] += fStdHistogram[j];//所有灰度级，关于w0的数组						
			fGrayAve[i] += j*fStdHistogram[j];//所有灰度级，关于u(t)的数组
		}
		fAverage += i*fStdHistogram[i];	
	}

	//计算OSTU
	for(i=0;i<=255;i++)
	{
		fTemp=(fAverage*fGrayAccu[i]-fGrayAve[i])*(fAverage*fGrayAccu[i]-fGrayAve[i])/(fGrayAccu[i]*(1-fGrayAccu[i]));
		if(fTemp>fMax)
		{
			fMax=fTemp;
			nThresh=i;
		}
	}
	return nThresh;
}