Ots算法(大律法或最大类间方差法)

1、Otsu 算法 ( 大律法或最大类间方差法 )一、 Otsu 最大类间方差法原理利用阈值将原图像分成前景，背景两个图象。前景：用 n1,csum,m1 来表示在当前阈值下的前景的点数，质量矩，平均灰度后景：用 n2, sum-csum,m2 来表示在当前阈值下的背景的点数，质量矩，平均灰度当取最佳阈值时，背景应该与前景差别最大，关键在于如何选择衡量差别的标准，而在 otsu 算法中这个衡量差别的标准就是最大类间 方差(英文简称 otsu ，这也就是这个算法名字的来源)，在本程序中类间方差用 sb 表示，最大类间方差用 fmax关于最大类间方差法( otsu )的性能 :类间方差法对噪音和目标大

2、小十分敏感，它仅对类间方差为单峰的图像产生较好的分割效果。当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好，但是类间方差法是用时最少的。最大类间方差法( otsu )的公式推导：记t为前景与背景的分割阈值，前景点数占图像比例为W0平均灰度为uO;背景点数占图像比例为w1，平均灰度为u1。则图像的总平均灰度为： u=w0*u0+w1*u1。前景和背景图象的方差： g=wO*(uO-u)*(uO-u)+w1*(u1-u)*(u1-u)=wO*w1*(uO-u1)*(uO-u1),此公式为方差公式。可参照概率论课本上面的g的公式也就是下面程序中的sb的表达式。当方差g最

3、大时，可以认为此时前景和背景差异最大，此时的灰度 t 是最佳阈值 sb = w1*w2*(u1-uO)*(uO-u1)算法实现 1 ：unsafe public int GetThreshValue(Bitmap image)BitmapData bd = image.LockBits(new Rectangle(O, O, image.Width, image.Height), ImageLockMode.WriteOnly, image.PixelFormat);byte* pt = (byte*)bd.ScanO; int pixelNum = new int256; byte colo

4、r;byte* pline;int n, n1, n2;int total;double m1, m2, sum, csum, fmax, sb; 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;/ 生成直方图

5、/图象直方图，共 256 个点/total 为总和，累计值/sb 为类间方差， fmax 存储最大方差值/ 阈值for (int i = O; i image.Height; i+)pline = pt + i * bd.Stride;for (int j = O; 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个灰度做平滑

6、化，t值应取较小的值q = k + t;if (q 255)q = 255;total = total + pixelNumq; /total 为总和，累计值/平滑化，左边2个+中间1个+右边2个灰度，共5个，所以总和除以 5,后面加0.5是用修正值pixelNumk = (int)(float)total / 5.0 + 0.5);/求阈值sum = csum = 0.0;n 二 0;/计算总的图象的点数和质量矩，为后面的计算做准备for (k = 0; k = 255; k+)/x*f(x)质量矩，也就是每个灰度的值乘以其点数(归一化后为概率)，sum为其总和sum += (double)

7、k * (double)pixelNumk;n += pixelNumk;fmax = -1.0;n1 = 0;for (k = 0; k fmax)fmax = sb; threshValue = k;如果算出的类间方差大于前一次算出的类间方差/fmax始终为最大类间方差(otsu)/取最大类间方差时对应的灰度的k就是最佳阈值 image.UnlockBits(bd);image.DisposeO; return threshValue;算法实现2：Otsu算法步骤如下：设图象包含L个灰度级(0,1丄-1)，灰度值为i的的象素点数为Ni，图象总的象素点数为N=N0+N1+.+N(L-1)。灰

8、度值为 i 的点的概率为：P(i) = N(i)/N.门限t将整幅图象分为暗区c1和亮区c2两类，则类间方差b是t的函数：b =a1*a2(u1-u2)A2式中，aj为类cj的面积与图象总面积之比，a仁sum(P(i) i-t, a2 = 1-a1;该法选择最佳门限tAuj 为类 cj 的均值，u1 = sum(i*P(i)/a1 0-t, u2 = sum(i*P(i)/a2, t+1-L-1, 使类间方差最大，即：令 u=u1-u2 , b b = maxa1(t)*a2(t) uA2代码实现：int otsu (Ipllmage *image, int rows, int cols, i

9、nt 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 ml, m2, sum, csum, fmax, sb;/对直方图置零memset(ihist, 0, sizeof(ihist);gmin=255; gmax=0;/生成直方图/*for (i = y0 + 1; i y0 + dy

10、- 1; i+)np = & imagei*cols+x0+1;for (j = x0 + 1; j gmax) gmax=* np;if(*np gmi n) gmi n=*np;np+; /* next pixel*/for(j=y0;jvdy;j+)for(i=0;idx;i+)unsigned char temp=CV_IMAGE_ELEM(image,uchar,j,i); ihisttemp+;/ set up everythingsum = csum = 0.0;n = 0;for (k = 0; k = 255; k+)sum += (double) k * (double) ihistk; / x*f(x) 质量矩n += ihistk;/f(x)质量if (!n)/ if n has no value, there is problemsfprintf (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 va