opencv最大类间方差法（otsu)实现（三种算法比较）

OTSU算法：就是计算出灰度图最佳阈值的算法

1.先对灰度图进行直方图计算并归一化处理，得到0-255之间每个像素在灰度图中出现的概率，即表示为某个像素在灰度图中出现了n个，灰度图总的像素点为N个，则这个像素的出现概率为Pi=n/N

2.每个灰度图可以由阈值k将灰度图分为A，B两大类，很容易得到A，B类在灰度图中的出现概率以及灰度均值

3.计算灰度图A,B类得类间方差，在最佳阈值K处，求得的类间方差最大，也就是类间方差最大的那个时刻的阈值就为灰度图的最佳阈值

代码实现：

以下三种算法时间复杂度依次降低，前两种算法自己所写，最后一种算法从网上看到大神写的，膜拜！当然如果你想偷懒，可是更直接点用cvThreshold一步实现cvThreshold(srcCopy,thresholdImage,0,255,CV_THRESH_OTSU)，花费时间与算法三相当。

算法一：

处理一张640*480的图片花费时间大概在8000毫秒左右，非常低效，但是是一种非常直接的，显而易见的实现方法。

int otsu(IplImage* src)
{
	CvScalar cs = cvSum(src);                             //图像总灰度值
	int Width = src->width;
	int Height = src->height;
	double U_t = 1.0*cs.val[0]/(Width*Height);                  //图像的平均灰度
	int threshold = 0;
	double delta = 0;
	for(int k=0; k<256; k++)
	{
		double sumOfGray = 0;
		int num = 0;
		for(int i=0; i < Width; i++)
		{
			for(int j=0; j < Height;j++)
			{
				int t = cvGet2D(src, j,i).val[0];
				if(t < k)
				{
					num++;             //统计灰度小于k的像素点的个数
					sumOfGray += t; //灰度值小于k的像素点的灰度和
				}
			}
		}
		double w0 = 1.0*num/(Width*Height);                     //灰度值小于阈值k的像素的概率
		double w1 = 1.0 - w0;                               //灰度值大于k像素的概率
		double u0 = 1.0*sumOfGray/num;                     //灰度值小于阈值k的像素的平均灰度值
		double u1 = 1.0*(cs.val[0]-sumOfGray)/(Width*Height-num); //灰度值大于阈值k的像素的平均灰度值

		double delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //类间方差
		
		if(delta_tmp > delta)
		{
			delta = delta_tmp;
			threshold = k;
		}
	}
	cout<<cs.val[0]<<endl;
	cout<<U_t<<endl;
	return threshold;
}

算法二：

借用opencv中的直方图函数，大大减少循环数，同一张图片花费时间大概在16毫秒左右。

int otsu(IplImage* src)  
{   
    int hist_size = 256;    //直方图尺寸    
    int hist_height = 256;    
    float range[] = {0,255};  //灰度级的范围    
    float* ranges[]={range};    
    //创建一维直方图，统计图像在[0 255]像素的均匀分布    
    CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);    
    //计算灰度图像的一维直方图    
    cvCalcHist(&src,gray_hist,0,0);    
    //归一化直方图    
    cvNormalizeHist(gray_hist,1.0);    
      
    int Width = src->width;  
    int Height = src->height;   
    int threshold = 0;  
    double delta = 0;
	CvScalar sumOfGray = cvSum(src);
	double U_t = sumOfGray.val[0]/(Width*Height);
	
    for(int k=0; k<256; k++)  
    {  
        double u0 = 0,  u1 = 0, w0 = 0, w1 = 0, delta_tmp = 0, u00 = 0, u11 = 0;  
         
        for(int i=0; i<k; i++)  
        {  
            u00 += cvQueryHistValue_1D(gray_hist,i)*i;   //灰度小于阈值k的像素的平均灰度值  
            w0 += cvQueryHistValue_1D(gray_hist,i);     //灰度小于阈值k的像素的概率  
        }
		
		u0 = u00/w0;

        for(int j=k; j<256; j++)  
        {  
            u11 += cvQueryHistValue_1D(gray_hist,j)*j;   //灰度大于阈值k的像素的平均灰度值  
            w1 += cvQueryHistValue_1D(gray_hist,j);     //灰度小于阈值k的像素的概率  
        }

		u1 = u11/w1;

        delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //类间方差  
        if(delta_tmp > delta)  
        {  
            delta = delta_tmp;  
            threshold = k;  
        }  
    } 

	cvReleaseHist(&gray_hist);
    return threshold;   
}

算法三：

一次循环实现，其中的公式需要自己用笔推导一下，非常精妙！完成一次运算的时间毫秒已经不能计量了，显示为0！

int otsu(IplImage* src)  
{   
    int hist_size = 256;    //直方图尺寸    
    int hist_height = 256;    
    float range[] = {0,255};  //灰度级的范围    
    float* ranges[]={range};    
    //创建一维直方图，统计图像在[0 255]像素的均匀分布    
    CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);    
    //计算灰度图像的一维直方图    
    cvCalcHist(&src,gray_hist,0,0);    
    //归一化直方图    
    cvNormalizeHist(gray_hist,1.0);    
      
    int Width = src->width;  
    int Height = src->height;   
    int threshold = 0;  
    double delta = 0;
	double U_t = 0;
	for(int m=0; m<256; m++)
	{
		U_t += cvQueryHistValue_1D(gray_hist,m)*m;
	}

	
	double u = 0, w = 0; 
    for(int k=0; k<256; k++)  
    {  
            
        u += cvQueryHistValue_1D(gray_hist,k)*k;   //  
        w += cvQueryHistValue_1D(gray_hist,k);      //灰度大于阈值k的像素的概率  

		double t = U_t * w - u;
		double delta_tmp = t * t / (w * (1 - w) ); 
		
        if(delta_tmp > delta)  
        {  
            delta = delta_tmp;  
            threshold = k;  
        }  
    } 

	cvReleaseHist(&gray_hist);
    return threshold;   
}