OpenCV实现验证otsu算法

转自：http://blog.csdn.net/augusdi/article/details/9012043

 

otsu法（最大类间方差法，有时也称之为大津算法）使用的是聚类的思想，把图像的灰度数按灰度级分成2个部分，使得两个部分之间的灰度值差异最大，每个部分之间的灰度差异最小，通过方差的计算来寻找一个合适的灰度级别 来划分。 所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响。因此,使类间方差最大的分割意味着错分概率最小。

设t为设定的阈值。

wo： 分开后 前景像素点数占图像的比例

uo： 分开后 前景像素点的平均灰度

w1：分开后 被景像素点数占图像的比例

u1： 分开后 被景像素点的平均灰度

u=w0*u0 + w1*u1 :图像总平均灰度

从L个灰度级遍历t，使得t为某个值的时候，前景和背景的方差最大， 则 这个 t 值便是我们要求得的阈值。

其中，方差的计算公式如下：

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

[ 此公式计算量较大，可以采用： g = wo * w1 * (uo - u1) * (uo - u1) ]

由于otsu算法是对图像的灰度级进行聚类，so 在执行otsu算法之前，需要计算该图像的灰度直方图。

按照上面的解释参考代码如下：

[cpp] view plain copy print ?

1. #include <stdio.h> 
2. #include <cv.h> 
3. #include <highgui.h> 
4. #include <math.h> 
5.  
6. #pragma comment(lib, "cv.lib") 
7. #pragma comment(lib, "cxcore.lib") 
8. #pragma comment(lib, "highgui.lib") 
9. int Otsu(IplImage* src); 
10.  
11. int _tmain(int argc, _TCHAR* argv[]) 
12. { 
13.     IplImage* img = cvLoadImage("c:\\aSa.jpg",0); 
14.     IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1); 
15.     int threshold = Otsu(img); 
16.  
17.     cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY); 
18.  
19.  
20.     cvNamedWindow( "img", 1 ); 
21.     cvShowImage("img", dst); 
22.  
23.  
24.     cvWaitKey(-1); 
25.  
26.     cvReleaseImage(&img); 
27.     cvReleaseImage(&dst); 
28.  
29.     cvDestroyWindow( "dst" ); 
30.     return 0; 
31. } 
32.  
33. int Otsu(IplImage* src)   
34. {   
35.     int height=src->height;   
36.     int width=src->width;       
37.     long size = height * width;  
38.  
39.     //histogram   
40.     float histogram[256] = {0};   
41.     for(int m=0; m < height; m++) 
42.     {   
43.         unsigned char* p=(unsignedchar*)src->imageData + src->widthStep * m;   
44.         for(int n = 0; n < width; n++)  
45.         {   
46.             histogram[int(*p++)]++;   
47.         }   
48.     }   
49.  
50.     int threshold;     
51.     long sum0 = 0, sum1 = 0;//存储前景的灰度总和和背景灰度总和 
52.     long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数 
53.     double w0 = 0, w1 = 0;//前景和背景所占整幅图像的比例 
54.     double u0 = 0, u1 = 0;  //前景和背景的平均灰度 
55.     double variance = 0;//最大类间方差 
56.     int i, j; 
57.     double u = 0; 
58.     double maxVariance = 0; 
59.     for(i = 1; i < 256; i++)//一次遍历每个像素 
60.     {   
61.         sum0 = 0; 
62.         sum1 = 0;  
63.         cnt0 = 0; 
64.         cnt1 = 0; 
65.         w0 = 0; 
66.         w1 = 0; 
67.         for(j = 0; j < i; j++) 
68.         { 
69.             cnt0 += histogram[j]; 
70.             sum0 += j * histogram[j]; 
71.         } 
72.  
73.         u0 = (double)sum0 /  cnt0;  
74.         w0 = (double)cnt0 / size; 
75.  
76.         for(j = i ; j <= 255; j++) 
77.         { 
78.             cnt1 += histogram[j]; 
79.             sum1 += j * histogram[j]; 
80.         } 
81.  
82.         u1 = (double)sum1 / cnt1; 
83.         w1 = 1 - w0; // (double)cnt1 / size; 
84.  
85.         u = u0 * w0 + u1 * w1; //图像的平均灰度 
86.         printf("u = %f\n", u); 
87.         //variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2); 
88.         variance =  w0 * w1 *  (u0 - u1) * (u0 - u1); 
89.         if(variance > maxVariance)  
90.         {   
91.             maxVariance = variance;   
92.             threshold = i;   
93.         }  
94.     }   
95.  
96.     printf("threshold = %d\n", threshold); 
97.     return threshold;   
98. }   

#include <stdio.h>
#include <cv.h>
#include <highgui.h>
#include <math.h>

#pragma comment(lib, "cv.lib")
#pragma comment(lib, "cxcore.lib")
#pragma comment(lib, "highgui.lib")
int Otsu(IplImage* src);

int _tmain(int argc, _TCHAR* argv[])
{
	IplImage* img = cvLoadImage("c:\\aSa.jpg",0);
	IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
	int threshold = Otsu(img);

	cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);


	cvNamedWindow( "img", 1 );
	cvShowImage("img", dst);


	cvWaitKey(-1);

	cvReleaseImage(&img);
	cvReleaseImage(&dst);

	cvDestroyWindow( "dst" );
	return 0;
}

int Otsu(IplImage* src)  
{  
	int height=src->height;  
	int width=src->width;      
	long size = height * width; 

	//histogram  
	float histogram[256] = {0};  
	for(int m=0; m < height; m++)
	{  
		unsigned char* p=(unsigned char*)src->imageData + src->widthStep * m;  
		for(int n = 0; n < width; n++) 
		{  
			histogram[int(*p++)]++;  
		}  
	}  

	int threshold;    
	long sum0 = 0, sum1 = 0; //存储前景的灰度总和和背景灰度总和
	long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数
	double w0 = 0, w1 = 0; //前景和背景所占整幅图像的比例
	double u0 = 0, u1 = 0;  //前景和背景的平均灰度
	double variance = 0; //最大类间方差
	int i, j;
	double u = 0;
	double maxVariance = 0;
	for(i = 1; i < 256; i++) //一次遍历每个像素
	{  
		sum0 = 0;
		sum1 = 0; 
		cnt0 = 0;
		cnt1 = 0;
		w0 = 0;
	    w1 = 0;
		for(j = 0; j < i; j++)
		{
			cnt0 += histogram[j];
			sum0 += j * histogram[j];
		}

		u0 = (double)sum0 /  cnt0; 
		w0 = (double)cnt0 / size;

		for(j = i ; j <= 255; j++)
		{
			cnt1 += histogram[j];
			sum1 += j * histogram[j];
		}

		u1 = (double)sum1 / cnt1;
		w1 = 1 - w0; // (double)cnt1 / size;

		u = u0 * w0 + u1 * w1; //图像的平均灰度
		printf("u = %f\n", u);
		//variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);
		variance =  w0 * w1 *  (u0 - u1) * (u0 - u1);
		if(variance > maxVariance) 
		{  
			maxVariance = variance;  
			threshold = i;  
		} 
	}  

	printf("threshold = %d\n", threshold);
	return threshold;  
}  

把w1写成w0 ··害我debug 了好久~~总是不认真，脑袋浑浑噩噩的···这都看不出来。。。。

==================

对了，之前搜集的一个otsu的算法，代码如下（由于时间太久了，不知道出处了。。。膜拜大牛哈）

[cpp] view plain copy print ?

1. <p>#include <stdio.h> 
2. #include <cv.h> 
3. #include <highgui.h> 
4. #include <math.h></p><p>#pragma comment(lib, "cv.lib") 
5. #pragma comment(lib, "cxcore.lib") 
6. #pragma comment(lib, "highgui.lib") 
7. </p>int Otsu(IplImage* src); 
8.  
9. int _tmain(int argc, _TCHAR* argv[]) 
10. { 
11.     IplImage* img = cvLoadImage("c:\\aSa.jpg",0); 
12.     IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1); 
13.     int threshold = Otsu(img); 
14.     printf("threshold = %d\n", threshold); 
15.     cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY); 
16.  
17.     cvNamedWindow( "img", 1 ); 
18.     cvShowImage("img", dst); 
19.  
20.  
21.     cvWaitKey(-1); 
22.  
23.     cvReleaseImage(&img); 
24.     cvReleaseImage(&dst); 
25.      
26.     cvDestroyWindow( "dst" ); 
27.     return 0; 
28. } 
29.  
30. int Otsu(IplImage* src)   
31. {   
32.     int height=src->height;   
33.     int width=src->width;       
34.  
35.     //histogram   
36.     float histogram[256] = {0};   
37.     for(int i=0; i < height; i++) 
38.     {   
39.         unsigned char* p=(unsignedchar*)src->imageData + src->widthStep * i;   
40.         for(int j = 0; j < width; j++)  
41.         {   
42.             histogram[*p++]++;   
43.         }   
44.     }   
45.     //normalize histogram   
46.     int size = height * width;   
47.     for(int i = 0; i < 256; i++) 
48.     {   
49.         histogram[i] = histogram[i] / size;   
50.     }   
51.  
52.     //average pixel value   
53.     float avgValue=0;   
54.     for(int i=0; i < 256; i++) 
55.     {   
56.         avgValue += i * histogram[i];  //整幅图像的平均灰度 
57.     }    
58.  
59.     int threshold;     
60.     float maxVariance=0;   
61.     float w = 0, u = 0;   
62.     for(int i = 0; i < 256; i++)  
63.     {   
64.         w += histogram[i];  //假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例 
65.         u += i * histogram[i];  // 灰度i 之前的像素(0~i)的平均灰度值： 前景像素的平均灰度值 
66.  
67.         float t = avgValue * w - u;   
68.         float variance = t * t / (w * (1 - w) );   
69.         if(variance > maxVariance)  
70.         {   
71.             maxVariance = variance;   
72.             threshold = i;   
73.         }   
74.     }   
75.  
76.     return threshold;   
77. }