OTSU算法學習 OTSU公式證明
1 otsu的公式如下,如果當前閾值為t,
w0 前景點所占比例
w1 = 1- w0 背景點所占比例
u0 = 前景灰度均值
u1 = 背景灰度均值
u = w0*u0 + w1*u1 全局灰度均值
g = w0(u0-u)*(u0-u) + w1(u1-u)*(u1-u) = w0*(1 – w0)*(u0 - u1)* (u0 - u1)
目標函數為g, g越大,t就是越好的閾值.為什么采用這個函數作為判別依據,直觀是這個函數反映了前景和背景的差值.
差值越大,閾值越好.
下面是一段證明g的推導的matlab代碼
| syms w0 u0 u1 %w0 前景均值 u0 前景灰度均值 u1 背景灰度均值 %背景均值 w1 =1- w0; %全局灰度均值 u=w0*u0+w1*u1; %目標函數 g=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u); %化簡的形式 g1 =w0*w1*(u0-u1)*(u0-u1); %因式展開 a1 = expand(g)%結果是 - u0^2*w0^2 + u0^2*w0 + 2*u0*u1*w0^2 - 2*u0*u1*w0 - u1^2*w0^2 + u1^2*w0 a2 = expand(g)%結果是 - u0^2*w0^2 + u0^2*w0 + 2*u0*u1*w0^2 - 2*u0*u1*w0 - u1^2*w0^2 + u1^2*w0 %對g進行因式分解 a2 = factor(g)%結果 -w0*(u0 - u1)^2*(w0 - 1) |
這里是matlab初等代數運算的講解 http://wenku.baidu.com/link?url=SODqdtPjbNLhKPEvCjsHkOhMi9LMb34qIrnp9_QRBKUNqPLGLxRCuLJgL2sp1vhLk55b6hpp242-RTCVp6ma_7a7-0imT3WVyBcsTmQ-5HS
2 關於最大類間方差法（otsu）的性能:
類間方差法對噪音和目標大小十分敏感，它僅對類間方差為單峰的圖像產生較好的分割效果。
當目標與背景的大小比例懸殊時，類間方差准則函數可能呈現雙峰或多峰，此時效果不好，但是類間方差法是用時最少的。
3 代碼實現
| 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; }
|
4 二維otsu算法
下圖是二維otsu的立體圖,是對右邊的A進行二維直方圖統計得到的圖像, 遍歷區域為5*5.
這是對應的matlab代碼
| %統計二維直方圖 i 當前點的亮度 j n*n鄰域均值亮度 function hist2 = hist2Function(image, n)
%初始化255*2565矩陣 hist2 = zeros(255,255);
[height, width]= size(image);
for i =1:height for j =1:width data = image(i,j);
tempSum =0.0; for l =-n:1:n for m =-n:1:n x = i + l; y = j+m;
if x <1 x =1; elseif x > width x = width; end
if y <1 y =1; elseif y > height y = height; end tempSum = tempSum + double(image(y,x)); end end
tempSum = tempSum /((2*n+1)*(2*n+1));
hist2(data,floor(tempSum))= hist2(data,floor(tempSum))+1; end end
%加載圖像 imagea = imread('a.bmp'); %顯示圖像 %imshow(imagea); %顯示直方圖 %figure;imhist(imagea); %計算二維直方圖 hist2 = hist2Function(imagea,2); %顯示二維直方圖 [x,y]=meshgrid(1:1:255); mesh(x,y,hist2)
|
<灰度圖象的二維Otsu自動閾值分割法.pdf> 這篇文章講解的不錯.文章這里有下載http://download.csdn.net/detail/wisdomfriend/9046341
下面用數學語言表達一下
i :表示亮度的維度
j : 表示點區域均值的維度
w0: 表示在閾值(s,t)時 所占的比例
w1: 表示在閾值(s,t)時, 所占的比例
u0(u0i, u0j): 表示在閾值(s,t)時時 的均值.u0時2維的
u1(u1i, u1j): 表示在閾值(s,t)時的均值.u1時2維的
uT: 全局均值
和一維otsut函數類似的目標函數
sb = w0(u0-uT)*(u0-uT)’ + w1(u1-uT)*(u1-uT)’
= w0[(u0i-uTi)* (u0i-uTi) + (u0j-uTj)* (u0j-uTj)] + w1[(u1i-uTi)* (u1i-uTi) + (u1j-uTj)* (u1j-uTj)]
這里是代碼實現.出自這篇文章:http://blog.csdn.net/yao_wust/article/details/23531031
| int histogram[256][256]; double p_histogram[256][256]; double Pst_0[256][256];//Pst_0用來存儲概率分布情況 double Xst_0[256][256];//存儲x方向上的均值矢量 int OTSU2d(IplImage * src) { int height = src->height; int width = src->width; long pixel = height * width;
int i,j;
for(i =0;i <256;i++)//初始化直方圖 { for(j =0; j <256;j++) histogram[i][j]=0; }
IplImage * temp = cvCreateImage(cvGetSize(src),8,1); cvSmooth(src,temp,CV_BLUR,3,0);
for(i =0;i < height;i++)//計算直方圖 { for(j =0; j < width;j++) { int data1 = cvGetReal2D(src,i,j); int data2 = cvGetReal2D(temp,i,j); histogram[data1][data2]++; } }
for(i =0; i <256;i++)//直方圖歸一化 for(j =0; j <256;j++) p_histogram[i][j]=(histogram[i][j]*1.0)/(pixel*1.0);
Pst_0[0][0]= p_histogram[0][0]; for(i =0;i <256;i++)//計算概率分布情況 for(j =0;j <256;j++) { double temp =0.0; if(i-1>=0) temp = temp + Pst_0[i-1][j]; if(j-1>=0) temp = temp + Pst_0[i][j-1]; if(i-1>=0&& j-1>=0) temp = temp - Pst_0[i-1][j-1]; temp = temp + p_histogram[i][j]; Pst_0[i][j]= temp; }
Xst_0[0][0]=0* Pst_0[0][0]; for(i =0; i <256;i++)//計算x方向上的均值矢量 for(j =0; j <256;j++) { double temp =0.0; if(i-1>=0) temp = temp + Xst_0[i-1][j]; if(j-1>=0) temp = temp + Xst_0[i][j-1]; if(i-1>=0&& j-1>=0) temp = temp - Xst_0[i-1][j-1]; temp = temp + i * p_histogram[i][j]; Xst_0[i][j]= temp; }
double Yst_0[256][256];//存儲y方向上的均值矢量 Yst_0[0][0]=0* Pst_0[0][0]; for(i =0; i <256;i++)//計算y方向上的均值矢量 for(j =0; j <256;j++) { double temp =0.0; if(i-1>=0) temp = temp + Yst_0[i-1][j]; if(j-1>=0) temp = temp + Yst_0[i][j-1]; if(i-1>=0&& j-1>=0) temp = temp - Yst_0[i-1][j-1]; temp = temp + j * p_histogram[i][j]; Yst_0[i][j]= temp; }
int threshold1; int threshold2; double variance =0.0; double maxvariance =0.0; for(i =0;i <256;i++)//計算類間離散測度 for(j =0;j <256;j++) { longdouble p0 = Pst_0[i][j]; longdouble v0 = pow(((Xst_0[i][j]/p0)-Xst_0[255][255]),2)+ pow(((Yst_0[i][j]/p0)-Yst_0[255][255]),2); longdouble p1 = Pst_0[255][255]-Pst_0[255][j]-Pst_0[i][255]+Pst_0[i][j]; longdouble vi = Xst_0[255][255]-Xst_0[255][j]-Xst_0[i][255]+Xst_0[i][j]; longdouble vj = Yst_0[255][255]-Yst_0[255][j]-Yst_0[i][255]+Yst_0[i][j]; longdouble v1 = pow(((vi/p1)-Xst_0[255][255]),2)+pow(((vj/p1)-Yst_0[255][255]),2);
variance = p0*v0+p1*v1;
if(variance > maxvariance) { maxvariance = variance; threshold1 = i; threshold2 = j; } }
//printf("%d %d",threshold1,threshold2);
return(threshold1+threshold2)/2;
}
|
總結: 二維otsu算法得到一個閾值,然后對圖像做二值化.仍然不能解決光照不均勻二值化的問題.比一維otsu效果好一些,但不是很明顯.這個算法的亮點在於考慮的點的附近區域的均值.