OTSU算法是由日本學者OTSU於1979年提出的一種對圖像進行二值化的高效算法。
基本介紹
- 外文名：Otsu
- 實質：大津法或最大類間方差法
- 所屬：算法
- 用於：計算機語言
OTSU算法（大津法或最大類間方差法）
一、Otsu最大類間方差法原理
利用閾值將原圖像分成前景，背景兩個圖象。
前景：用n1,csum,m1來表示在當前閾值下的前景的點數，質量矩，平均灰度
背景：用n2, sum-csum,m2來表示在當前閾值下的背景的點數，質量矩，平均灰度
當取最佳閾值時，背景應該與前景差別最大，關鍵在於如何選擇衡量差別的標準，而在otsu算法中這個衡量差別的標準就是最大類間方差，在本程式中類間方差用sb表示，最大類間方差用fmax
關於最大類間方差法（otsu）的性能:
類間方差法對噪音和目標大小十分敏感，它僅對類間方差為單峰的圖像產生較好的分割效果。
當目標與背景的大小比例懸殊時，類間方差準則函式可能呈現雙峰或多峰，此時效果不好，但是類間方差法是用時最少的。
最大類間方差法（otsu）的公式推導：
記t為前景與背景的分割閾值，前景點數占圖像比例為w0，平均灰度為u0；背景點數占圖像比例為w1，平均灰度為u1。
則圖像的總平均灰度為：u=w0*u0+w1*u1。
前景和背景圖象的方差：g=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)=w0*w1*(u0-u1)*(u0-u1),此公式為方差公式。
可參照機率論課本上面的g的公式也就是下面程式中的sb的表達式。當方差g最大時，可以認為此時前景和背景差異最大，此時的灰度t是最佳閾值sb = w0*w1*(u1-u0)*(u0-u1)
算法實現1：
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 int[256]; //圖象直方圖，共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表示 pixelNum[color]++; //相應的直方圖加1 } } //直方圖平滑化 for (k = 0; k <= 255; k++) { total = 0; for (t = -2; t <= 2; t++) //與附近2個灰度做平滑化，t值應取較小的值 { q = k + t; if (q < 0) //越界處理 q = 0; if (q > 255) q = 255; total = total + pixelNum[q]; //total為總和，累計值 } //平滑化，左邊2個+中間1個+右邊2個灰度，共5個，所以總和除以5，後面加0.5是用修正值 pixelNum[k] = (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)k * (double)pixelNum[k]; n += pixelNum[k]; //n為圖象總的點數，歸一化後就是累積機率 } fmax = -1.0; //類間方差sb不可能為負，所以fmax初始值為-1不影響計算的進行 n1 = 0; for (k = 0; k < 255; k++) //對每個灰度（從0到255）計算一次分割後的類間方差sb { n1 += pixelNum[k]; //n1為在當前閾值遍前景圖象的點數 if (n1 == 0) { continue; } //沒有分出前景後景 n2 = n - n1; //n2為背景圖象的點數 //n2為0表示全部都是後景圖象，與n1=0情況類似，之後的遍歷不可能使前景點數增加，所以此時可以退出循環 if (n2 == 0) { break; } csum += (double)k * pixelNum[k]; //前景的“灰度的值*其點數”的總和 m1 = csum / n1; //m1為前景的平均灰度 m2 = (sum - csum) / n2; //m2為背景的平均灰度 sb = (double)n1 * (double)n2 * (m1 - m2) * (m1 - m2); //sb為類間方差 if (sb > fmax) //如果算出的類間方差大於前一次算出的類間方差 { fmax = sb; //fmax始終為最大類間方差（otsu） threshValue = k; //取最大類間方差時對應的灰度的k就是最佳閾值 } } image.UnlockBits(bd); image.Dispose(); return threshValue;}算法實現2：
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 (i = y0 + 1; i < y0 + dy - 1; i++) { np = &image[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++; /* next pixel } }*/ for(j=y0;j<dy;j++) { for(i=0;i<dx;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]; // x*f(x) 質量矩 n += ihist[k]; //f(x) 質量 } if (!n) { // if n has no value, there is problems fprintf (stderr, "NOT NORMAL thresholdValue = 160\n"); 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); /**//* bbg: note: can be optimized. */ 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=%d\n",thresholdValue, gmin, gmax); return(thresholdValue);}