Otsu算法步骤如下:
设图象包含L个灰度级(0,1…,L-1),灰度值为i的的象素点数为Ni
,图象总的象素点数为N=N0+N1+...+N(L-1)。灰度值为i的点的概率为:
P(i) = N(i)/N.
门限t将整幅图象分为暗区c1和亮区c2两类,则类间方差σ是t的函数:
σ=a1*a2(u1-u2)^2 (2)
式中,aj 为类cj的面积与图象总面积之比,a1 = sum(P(i)) i->t, a2 =
1-a1; uj为类cj的均值,u1 = sum(i*P(i))/a1 0->t,
u2 = sum(i*P(i))/a2, t+1->L-1
该法选择最佳门限t^ 使类间方差最大,即:
令Δu=u1-u2,σb = max{a1(t)*a2(t)Δu^2}
代码实现:
int otsu (IplImage *image, int rows, int cols, int x0, int y0, int
dx, int dy, int vvv)
{
unsigned char *np; // 图像指针
int thresholdValue=1; // 阈值
int ihist[256]; // 图像直方图,256个点
int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;
// 对直方图置零
memset(ihist, 0, sizeof(ihist));
gmin=255; gmax=0;
// 生成直方图
for(j=y0;j
{
for(i=0;i
{
unsigned char temp=CV_IMAGE_ELEM(image,uchar,j,i);
ihist[temp]++;
}
}
// set up everything
sum = csum = 0.0;
n = 0;
for (k = 0; k <= 255; k++) {
sum += (double) k * (double) ihist[k];
n += ihist[k];
}
if (!n) {
// if n has no value, there is problems
fprintf (stderr, "NOT NORMAL thresholdValue = 160n");
return (160);
}
// do the otsu global thresholding method
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;
}
}
// at this point we have our thresholding value
// debug code to display thresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d
gmax=%dn",thresholdValue, gmin, gmax);
return(thresholdValue);
}