大津法(OTSU法)
大津法由大津于1979年提出,对图像Image,记t为前景与背景的分割阈值,前景点数占图像比例为w0,平均灰度为u0;背景点数占图像比例为w1,平均灰度为u1。图像的总平均灰度为:u=w0*u0+w1*u1。从最小灰度值到最大灰度值遍历t,当t使得值g=w0*(u0-u)2+w1*(u1-u)2 最大时t即为分割的最佳阈值。对大津法可作如下理解:该式实际上就是类间方差值,阈值t分割出的前景和背景两部分构成了整幅图像,而前景取值u0,概率为 w0,背景取值u1,概率为w1,总均值为u,根据方差的定义即得该式。因方差是灰度分布均匀性的一种度量,方差值越大,说明构成图像的两部分差别越大, 当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小,因此使类间方差最大的分割意味着错分概率最小。
直接应用大津法计算量较大,因此我们在实现时采用了等价的公式g=w0*w1*(u0-u1)2。部分计算过程如下:
-
- for intCurrentLevel:=0 to intArrLen do
- begin
- if intSclGrayLevel[intCurrentLevel]=0 then
- continue
- else
- begin
-
- intCount:=0;
- intSumPels:=0;
-
- for intLoop:=0 to intCurrentLevel do
- begin
- intCount:=intCount+intSclGrayLevel[intLoop];
- intSumPels:=intSumPels+intSumPelsArr[intLoop];
- end;
- w0:=intCount/intSize;
- u0:=intSumPels/intCount;
- w1:=1-w0;
- if intSize-intCount<>0 then
- u1:=(intTotalPels-intSumPels)/(intSize-intCount)
- else
- u1:=0;
-
- RlTempO:=w0*w1*(u0-u1)*(u0-u1);
- if RlTempO>RlMaxO then
- begin
- RlMaxO:=RlTempO;
- Result:=intCurrentLevel;
- end;
- end;
我们在测试中发现:大津法选取出来的阈值非常理想,对各种情况的表现都较为良好。虽然它在很多情况下都不是最佳的分割,但分割质量通常都有一定的保障,可以说是最稳定的分割。由上可知,大津算法是一种较为通用的分割算法。在它的思想的启迪下,人们进一步提出了多种类似的评估阈值的算法,具体可参加【5】、【6】等。
OTSU的算法,很好用,好不容易才找到的。
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
- {
- unsigned char *np;
- int thresholdValue=1;
- int ihist[256];
-
- int i, j, k;
- int n, n1, n2, gmin, gmax;
- double m1, m2, sum, csum, fmax, sb;
-
-
- memset(ihist, 0, sizeof(ihist));
-
- gmin=255; gmax=0;
-
- for (i = y0 + 1; i < y0 + dy - 1; i++) {
- np = [i*cols+x0+1];
- for (j = x0 + 1; j < x0 + dx - 1; j++) {
- ihist[*np]++;
- if(*np > gmax) gmax=*np;
- if(*np < gmin) gmin=*np;
- np++;
- }
- }
-
-
- sum = csum = 0.0;
- n = 0;
-
- for (k = 0; k <= 255; k++) {
- sum += (double) k * (double) ihist[k];
- n += ihist[k];
- }
-
- if (!n) {
-
- fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
- return (160);
- }
-
-
- fmax = -1.0;
- n1 = 0;
- for (k = 0; k < 255; k++) {
- n1 += ihist[k];
- if (!n1) { continue; }
- n2 = n - n1;
- if (n2 == 0) { break; }
- csum += (double) k *ihist[k];
- m1 = csum / n1;
- m2 = (sum - csum) / n2;
- sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
-
- if (sb > fmax) {
- fmax = sb;
- thresholdValue = k;
- }
- }
-
-
-
-
- if ( vvv & 1 )
- fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
- thresholdValue, gmin, gmax);
-
- return(thresholdValue);
- }
(Ryan Dibble) |