OTSU算法提取图像阈值的C语言实现

OTSU算法也称最大类间差法，有时也称之为大津算法，被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。它是按图像的灰度特性,将图像分成背景和前景两部分。背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。

设灰度图像灰度级是L，则灰度范围为[0,L-1]，利用OTSU算法计算图像的最佳阈值为：
t = Max[w0(t) * (u0(t) - u)^2 + w1(t) * (u1(t) - u)^2)]
其中的变量说明：当分割的阈值为t时，w0为背景比例，u0为背景均值，w1为前景比例，u1为前景均值，u为整幅图像的均值。
使以上表达式值最大的t，即为分割图像的最佳阈值。

以下是一段在OpenCV中实现的C语言程序，即一个使用OTSU算法提取图像阈值的函数，输入参数为一个图像指针，返回分割该图像的最佳阈值。

[cpp]  view plain copy

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