otsu自适应阈值分割的算法描述和opencv实现,及其在肤色检测中的应用

原文地址如下：http://blog.csdn.net/onezeros/article/details/6136770

otsu算法选择使类间方差最大的灰度值为阈值，具有很好的效果
算法具体描述见otsu论文，或冈萨雷斯著名的数字图像处理那本书
这里给出程序流程：
1、计算直方图并归一化histogram
2、计算图像灰度均值avgValue.
3、计算直方图的零阶w[i]和一级矩u[i]
4、计算并找到最大的类间方差（between-class variance）
variance[i]=(avgValue*w[i]-u[i])*(avgValue*w[i]-u[i])/(w[i]*(1-w[i]))
对应此最大方差的灰度值即为要找的阈值
5、用找到的阈值二值化图像

我在代码中做了一些优化，所以算法描述的某些地方跟程序并不一致

 

otsu代码，先找阈值，继而二值化

[cpp]  view plain copy

1. // implementation of otsu algorithm  
2. // author: onezeros(@yahoo.cn)  
3. // reference: Rafael C. Gonzalez. Digital Image Processing Using MATLAB  
4. void cvThresholdOtsu(IplImage* src, IplImage* dst)  
5. {  
6.     int height=src->height;  
7.     int width=src->width;      
8.       
9.     //histogram  
10.     float histogram[256]={0};  
11.     for(int i=0;i<height;i++) {  
12.         unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;  
13.         for(int j=0;j<width;j++) {  
14.             histogram[*p++]++;  
15.         }  
16.     }  
17.     //normalize histogram  
18.     int size=height*width;  
19.     for(int i=0;i<256;i++) {  
20.         histogram[i]=histogram[i]/size;  
21.     }  
22.       
23.     //average pixel value  
24.     float avgValue=0;  
25.     for(int i=0;i<256;i++) {  
26.         avgValue+=i*histogram[i];  
27.     }  
28.   
29.     int threshold;    
30.     float maxVariance=0;  
31.     float w=0,u=0;  
32.     for(int i=0;i<256;i++) {  
33.         w+=histogram[i];  
34.         u+=i*histogram[i];  
35.   
36.         float t=avgValue*w-u;  
37.         float variance=t*t/(w*(1-w));  
38.         if(variance>maxVariance) {  
39.             maxVariance=variance;  
40.             threshold=i;  
41.         }  
42.     }  
43.   
44.     cvThreshold(src,dst,threshold,255,CV_THRESH_BINARY);  
45. }  

 

更多情况下我们并不需要对每一帧都是用otsu寻找阈值，于是可以先找到阈值，然后用找到的阈值处理后面的图像。下面这个函数重载了上面的，返回值就是阈值。只做了一点改变

[cpp]  view plain copy

1. // implementation of otsu algorithm  
2. // author: onezeros(@yahoo.cn)  
3. // reference: Rafael C. Gonzalez. Digital Image Processing Using MATLAB  
4. int cvThresholdOtsu(IplImage* src)  
5. {  
6.     int height=src->height;  
7.     int width=src->width;      
8.   
9.     //histogram  
10.     float histogram[256]={0};  
11.     for(int i=0;i<height;i++) {  
12.         unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;  
13.         for(int j=0;j<width;j++) {  
14.             histogram[*p++]++;  
15.         }  
16.     }  
17.     //normalize histogram  
18.     int size=height*width;  
19.     for(int i=0;i<256;i++) {  
20.         histogram[i]=histogram[i]/size;  
21.     }  
22.   
23.     //average pixel value  
24.     float avgValue=0;  
25.     for(int i=0;i<256;i++) {  
26.         avgValue+=i*histogram[i];  
27.     }  
28.   
29.     int threshold;    
30.     float maxVariance=0;  
31.     float w=0,u=0;  
32.     for(int i=0;i<256;i++) {  
33.         w+=histogram[i];  
34.         u+=i*histogram[i];  
35.   
36.         float t=avgValue*w-u;  
37.         float variance=t*t/(w*(1-w));  
38.         if(variance>maxVariance) {  
39.             maxVariance=variance;  
40.             threshold=i;  
41.         }  
42.     }  
43.   
44.     return threshold;  
45. }  

 

我在手的自动检测中使用这个方法，效果很好。

下面是使用上述两个函数的简单的主程序，可以试运行一下，如果处理视频，要保证第一帧时，手要在图像中。

 

[cpp]  view plain copy

1. #include <cv.h>  
2. #include <cxcore.h>  
3. #include <highgui.h>  
4. #pragma comment(lib,"cv210d.lib")  
5. #pragma comment(lib,"cxcore210d.lib")  
6. #pragma comment(lib,"highgui210d.lib")  
7.   
8. #include <iostream>  
9. using namespace std;  
10.   
11. int main(int argc, char** argv)  
12. {     
13. #ifdef VIDEO //video process  
14.     CvCapture* capture=cvCreateCameraCapture(-1);  
15.     if (!capture){  
16.         cout<<"failed to open camera"<<endl;  
17.         exit(0);  
18.     }  
19.   
20.     int threshold=-1;  
21.     IplImage* img;    
22.     while (img=cvQueryFrame(capture)){  
23.         cvShowImage("video",img);  
24.         cvCvtColor(img,img,CV_RGB2YCrCb);  
25.   
26.         IplImage* imgCb=cvCreateImage(cvGetSize(img),8,1);  
27.         cvSplit(img,NULL,NULL,imgCb,NULL);  
28.         if (threshold<0){  
29.             threshold=cvThresholdOtsu(imgCb);  
30.         }  
31.         //cvThresholdOtsu(imgCb,imgCb);  
32.         cvThreshold(imgCb,imgCb,threshold,255,CV_THRESH_BINARY);  
33.         cvErode(imgCb,imgCb);  
34.         cvDilate(imgCb,imgCb);  
35.           
36.         cvShowImage("object",imgCb);  
37.         cvReleaseImage(&imgCb);  
38.   
39.         if (cvWaitKey(3)==27){//esc  
40.             break;  
41.         }  
42.     }     
43.   
44.     cvReleaseCapture(&capture);  
45.       
46. #else //single image process  
47.     const char* filename=(argc>=2?argv[1]:"cr.jpg");  
48.     IplImage* img=cvLoadImage(filename,CV_LOAD_IMAGE_GRAYSCALE);  
49.   
50.     cvThresholdOtsu(img,img);  
51.     cvShowImage( "src", img );  
52.     char buf[256];  
53.     sprintf_s(buf,256,"%s.otsu.jpg",filename);  
54.     cvSaveImage(buf,img);  
55.   
56.     cvErode(img,img);  
57.     cvDilate(img,img);  
58.     cvShowImage( "dst", img );  
59.     sprintf_s(buf,256,"%s.otsu.processed.jpg",filename);  
60.     cvSaveImage(buf,img);  
61.   
62.     cvWaitKey(0);  
63. #endif  
64.       
65.     return 0;  
66. }  

 

 

效果图：

1、肤色cb分量

 

2、otsu自适应阈值分割效果

 

3、开运算后效果