最大类间方差法.doc

最大类间方差法是由日本学者大津于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。对于图像I(x,y),前景(即目标)和背景的分割阈值记作T,属于前景的像素点数占整幅图像的比例记为0,其平均灰度0;背景像素点数占整幅图像的比例为1,其平均灰度为1。图像的总平均灰度记为,类间方差记为g。假设图像的背景较暗,并且图像的大小为MN,图像中像素的灰度值小于阈值T的像素个数记作N0,像素灰度大于阈值T的像素个数记作N1,则有:0=N0/ MN (1)1=N1/ MN (2)N0+N1=MN (3)0+1=1 (4)=0*0+1*1 (5)g=0(0-)2+1(1-)2 (6)将式(5)代入式(6),得到等价公式: g=01(0-1)2 (7)采用遍历的方法得到使类间方差最大的阈值T,即为所求。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 = maxa1(t)*a2(t)u2OTSU方法计算图像二值化的自适应阈值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 ihist256; / 图像直方图，256个点int i, j, k; / various countersint 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 = &imagei*cols+x0+1;for (j = x0 + 1; j gmax) gmax=*np;if(*np gmin) gmin=*np;np+;/ set up everythingsum = csum = 0.0;n = 0;for (k = 0; k = 255; k+) sum += (double) k * (double) ihistk;n += ihistk;if (!n) / if n has no value, there is problems.fprintf (stderr, NOT NORMAL thresholdValue = 160n);return (160);/ do the otsu global thresholding methodfmax = -1.0;n1 = 0;for (k = 0; k fmax) fmax = sb;thresholdValue = k;/ at this point we have our thresholding value/ debug code to display thresholding valuesif ( vvv & 1 )fprintf(stderr,# OTSU: thresholdValue = %d gmin=%d gmax=%dn,thresholdValue, gmin, gmax);return(thresholdValue);/另外一个/unsafe public int GetThreshValue(Bitmap image)BitmapData bd = image.LockBits(new Rectangle(0, 0, image.Width, image.Height), ImageLockMode.WriteOnly, image.PixelFormat);byte* pt = (byte*)bd.Scan0;int pixelNum = new int256;/图象直方图，共256个点byte color;byte* pline;int n, n1, n2;int total;/total为总和，累计值double m1, m2, sum, csum, fmax, sb;/sb为类间方差，fmax存储最大方差值int k, t, q;int threshValue = 1;/ 阈值int step = 1;switch (image.PixelFormat)case PixelFormat.Format24bppRgb:step = 3;break;case PixelFormat.Format32bppArgb:step = 4;break;case PixelFormat.Format8bppIndexed:step = 1;break;/生成直方图for (int i = 0; i image.Height; i+)pline = pt + i * bd.Stride;for (int j = 0; j image.Width; j+)color = *(pline + j * step);/返回各个点的颜色，以RGB表示pixelNumcolor+;/相应的直方图加1/直方图平滑化for (k = 0; k = 255; k+)total = 0;for (t = -2; t = 2; t+)/与附近2个灰度做平滑化，t值应取较小的值q = k + t;if (q 255)q = 255;total = total + pixelNumq;/total为总和，累计值pixelNumk = (int)(float)total / 5.0 + 0.5);/平滑化，左边2个+中间1个+右边2个灰度，共5个，所以总和除以5，后面加0.5是用修正值/求阈值sum = csum = 0.0;n = 0;/计算总的图象的点数和质量矩，为后面的计算做准备for (k = 0; k = 255; k+)sum += (double)k * (double)pixelNumk;/x*f(x)质量矩，也就是每个灰度的值乘以其点数（归一化后为概率），sum为其总和n += pixelNumk;/n为图象总的点数，归一化后就是累积概率fmax = -1.0;/类间方差sb不可能为负，所以fmax初始值为-1不影响计算的进行n1 = 0;for (k