最大类间方差法(大津法，OTSU)

最大类间方差法是由日本学者大津(Nobuyuki Otsu)于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。对于图像I(x,y),前景(即目标)和背景的分割阈值记作,属于前景的像素点数占整幅图像的比例记为,其平均灰度;背景像素点数占整幅图像的比例为,其平均灰度为。图像的总平均灰度记为,类间方差记为。

假设图像的背景较暗,并且图像的大小为,图像中像素的灰度值小于阈值的像素个数记作,像素灰度大于阈值的像素个数记作,则有:

       

采用遍历的方法得到使类间方差最大的阈值,即为所求。

代码区

[cpp] view plain copy

print?

1. OpenCV代码:  
2. int myOtsu(const IplImage *frame) //大津法求阈值  
3. {  
4. #define GrayScale 256   //frame灰度级  
5.     int width = frame->width;  
6.     int height = frame->height;  
7.     int pixelCount[GrayScale]={0};  
8.     float pixelPro[GrayScale]={0};  
9.     int i, j, pixelSum = width * height, threshold = 0;  
10.     uchar* data = (uchar*)frame->imageData;  
11.   
12.     //统计每个灰度级中像素的个数  
13.     for(i = 0; i < height; i++)  
14.     {  
15.         for(j = 0;j < width;j++)  
16.         {  
17.             pixelCount[(int)data[i * width + j]]++;  
18.         }  
19.     }  
20.   
21.     //计算每个灰度级的像素数目占整幅图像的比例  
22.     for(i = 0; i < GrayScale; i++)  
23.     {  
24.         pixelPro[i] = (float)pixelCount[i] / pixelSum;  
25.     }  
26.   
27.     //遍历灰度级[0,255],寻找合适的threshold  
28.     float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0;  
29.     for(i = 0; i < GrayScale; i++)  
30.     {  
31.         w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0;  
32.         for(j = 0; j < GrayScale; j++)  
33.         {  
34.             if(j <= i)   //背景部分  
35.             {  
36.                 w0 += pixelPro[j];  
37.                 u0tmp += j * pixelPro[j];  
38.             }  
39.             else   //前景部分  
40.             {  
41.                 w1 += pixelPro[j];  
42.                 u1tmp += j * pixelPro[j];  
43.             }  
44.         }  
45.         u0 = u0tmp / w0;  
46.         u1 = u1tmp / w1;  
47.         deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)) ;  
48.         if(deltaTmp > deltaMax)  
49.         {  
50.             deltaMax = deltaTmp;  
51.             threshold = i;  
52.         }  
53.     }  
54.     return threshold;  
55. }  

OpenCV代码:int 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* data = (uchar*)frame->imageData;//统计每个灰度级中像素的个数for(i = 0; i < height; i++){for(j = 0;j < width;j++){pixelCount[(int)data[i * width + j]]++;}}//计算每个灰度级的像素数目占整幅图像的比例for(i = 0; i < GrayScale; i++){pixelPro[i] = (float)pixelCount[i] / pixelSum;}//遍历灰度级[0,255],寻找合适的thresholdfloat 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;}

[plain] view plain copy

print?

1. MATLAB 代码1  
2.   
3. function binariedImage=OSTU(scoreImage)  
4.   
5. scoreImage=double(scoreImage);  
6. [height,length]=size(scoreImage);  
7. totalPixel=height*length;  
8. % maxPixNumber=max(max(scoreImage));  
9. % 这个地方为了以后计算方便 就不这样计算了 而是默认最大的为255  
10.   
11. pixelCount=zeros(1,256);%统计各个像素值的个数  
12. % 0-256  
13. for i=1:length  
14.     for j=1:height  
15.         number=scoreImage(j,i)+1;  
16.         pixelCount(number)=pixelCount(number)+1;  
17.     end  
18. end  
19.   
20. %概率  
21. pf=pixelCount/totalPixel;  
22. %前向累计概率密度函数  
23. cpf=zeros(1,256);  
24. cpf(1)=pf(1);  
25. for i=2:256  
26.     cpf(i)=cpf(i-1)+pf(i);  
27. end  
28. %后向累计概率密度函数  
29. bpf=zeros(1,256);  
30. bpf(256)=pf(256);  
31. for j=256:-1:2  
32.     bpf(j-1)=bpf(j)+pf(j-1);  
33. end   
34.   
35.   
36. %前向累计期望  
37. meanForward=zeros(1,256);  
38. meanForward(1)=1;  
39. for i=2:256  
40.     meanForward(i)=meanForward(i-1)*(cpf(i-1)/cpf(i))+(pf(i)/cpf(i))*i;  
41. end  
42.   
43. %后向累计期望  
44. meanBack=zeros(1,256);  
45. meanBack(max(max(scoreImage)))=max(max(scoreImage));  
46. for i=max(max(scoreImage)):-1:2  
47.     meanBack(i-1)=meanBack(i)*(bpf(i)/bpf(i-1))+(pf(i-1)/bpf(i-1))*(i-1);  
48. end  
49.   
50.   
51. delta=cpf.*bpf.*(meanForward-meanBack).^2;  
52. [value,index]=max(delta);  
53. % 返回的是逻辑的图像  
54. binariedImage=scoreImage>index;  

MATLAB 代码1function binariedImage=OSTU(scoreImage)scoreImage=double(scoreImage);[height,length]=size(scoreImage);totalPixel=height*length;% maxPixNumber=max(max(scoreImage));% 这个地方为了以后计算方便 就不这样计算了 而是默认最大的为255pixelCount=zeros(1,256);%统计各个像素值的个数% 0-256for i=1:length    for j=1:height        number=scoreImage(j,i)+1;        pixelCount(number)=pixelCount(number)+1;    endend%概率pf=pixelCount/totalPixel;%前向累计概率密度函数cpf=zeros(1,256);cpf(1)=pf(1);for i=2:256    cpf(i)=cpf(i-1)+pf(i);end%后向累计概率密度函数bpf=zeros(1,256);bpf(256)=pf(256);for j=256:-1:2    bpf(j-1)=bpf(j)+pf(j-1);end %前向累计期望meanForward=zeros(1,256);meanForward(1)=1;for i=2:256    meanForward(i)=meanForward(i-1)*(cpf(i-1)/cpf(i))+(pf(i)/cpf(i))*i;end%后向累计期望meanBack=zeros(1,256);meanBack(max(max(scoreImage)))=max(max(scoreImage));for i=max(max(scoreImage)):-1:2    meanBack(i-1)=meanBack(i)*(bpf(i)/bpf(i-1))+(pf(i-1)/bpf(i-1))*(i-1);enddelta=cpf.*bpf.*(meanForward-meanBack).^2;[value,index]=max(delta);% 返回的是逻辑的图像binariedImage=scoreImage>index;

[plain] view plain copy

print?

1.   

[plain] view plain copy

print?

1. MATLAB 代码2  
2.   
3. function binariedImage=OSTU(scoreImage)  
4.   
5.   
6. [height,length]=size(scoreImage);  
7. totalNum=height*length;  
8. pixelCount=zeros(1,256);%统计各个像素值的个数  
9. for i=1:length  
10.     for j=1:height  
11.         number=scoreImage(j,i)+1;  
12.         pixelCount(number)=pixelCount(number)+1;  
13.     end  
14. end  
15. pf=pixelCount/totalNum;  
16. deltaMax=0;  
17. a=1:256;  
18. maxmax=max(max(scoreImage));  
19. for i=1:maxmax  
20.       
21.     w1=sum(pf(1:i));  
22.     w2=1-w1;  
23.     mean1=sum(pf(1:i).*a(1:i))/sum(pf(1:i));  
24.     mean2=sum(pf(i+1:maxmax).*a(i+1:maxmax))/sum(pf(i+1:maxmax));  
25.     delta=w1*w2*(mean1-mean2)^2;  
26.     if delta>deltaMax  
27.         deltaMax=delta;  
28.         threshold=i;  
29.     end  
30.       
31. end  
32.   
33. % 返回的是逻辑的图像  
34. binariedImage=scoreImage>threshold;  

MATLAB 代码2function binariedImage=OSTU(scoreImage)[height,length]=size(scoreImage);totalNum=height*length;pixelCount=zeros(1,256);%统计各个像素值的个数for i=1:length    for j=1:height        number=scoreImage(j,i)+1;        pixelCount(number)=pixelCount(number)+1;    endendpf=pixelCount/totalNum;deltaMax=0;a=1:256;maxmax=max(max(scoreImage));for i=1:maxmax        w1=sum(pf(1:i));    w2=1-w1;    mean1=sum(pf(1:i).*a(1:i))/sum(pf(1:i));    mean2=sum(pf(i+1:maxmax).*a(i+1:maxmax))/sum(pf(i+1:maxmax));    delta=w1*w2*(mean1-mean2)^2;    if delta>deltaMax        deltaMax=delta;        threshold=i;    end    end% 返回的是逻辑的图像binariedImage=scoreImage>threshold;