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 – 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.Width, image.Height), ImageLockMode.WriteOnly, image.PixelFormat); byte* pt =(byte*)bd.Scan0; 个点 byte color; byte* pline; int n, n1, n2; int total;//total为总和，累计值 double m1, m2, sum, csum, fmax, sb;//sb为类间方差，fmax存储最大方差值 int k, t, q; ;// 阈值 ; switch(image.PixelFormat) { case PixelFormat.Format24bppRgb: step ; break; case PixelFormat.Format32bppArgb: step ; break; case PixelFormat.Format8bppIndexed: step ; break; } //生成直方图 ; i < image.Height; i++) { pline = pt + i * bd.Stride; ; j < image.Width; j++) { color =*(pline + j * step);//返回各个点的颜色，以RGB表示 pixelNum } } //直方图平滑化 ; k++) { total ; 个灰度做平滑化，t值应取较小的值 { q = k + t; )//越界处理 q ; ) q ; total = total + pixelNum[q];//total为总和，累计值 } //个+中间1个+右边2个灰度，共5个，所以总和除以5，后面加0.5是用修正值 pixelNum[k]=(int)((float)total /5.0+0.5); } //求阈值 sum = csum =0.0; n ; //计算总的图象的点数和质量矩，为后面的计算做准备 ; k++) { //x*f(x)质量矩，也就是每个灰度的值乘以其点数（归一化后为概率），sum为其总和 sum +=(double)k *(double)pixelNum[k]; n += pixelNum[k];//n为图象总的点数，归一化后就是累积概率 } fmax =–1.0;//类间方差sb不可能为负，所以fmax初始值为-1不影响计算的进行 n1 ; 到255）计算一次分割后的类间方差sb { n1 += pixelNum[k];//n1为在当前阈值遍前景图象的点数 ){continue;}//没有分出前景后景 n2 = n – n1;//n2为背景图象的点数 //n2表示全部都是后景图象，与n1=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 ); [height, width]= size(image); :height :width data = image(i,j); tempSum =0.0; :n :n x = i + l; y = j+m;
x ; elseif x > width x = width; end
y ; elseif y > height y = height; end tempSum = tempSum + double(image(y,x)); end end tempSum )); hist2; end end %加载图像 imagea = imread(‘a.bmp’); %显示图像 %imshow(imagea); %显示直方图 %figure;imhist(imagea); %计算二维直方图 hist2 ); %显示二维直方图 ); 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
|
]; ]; ];//Pst_0用来存储概率分布情况 ];//存储x方向上的均值矢量 int OTSU2d(IplImage * src) { int height = src->height; int width = src->width; long pixel = height * width; int i,j; ;i++)//初始化直方图 { ;j++) histogram; } IplImage ); cvSmooth); ;i < height;i++)//计算直方图 { ; j < width;j++) { int data1 = cvGetReal2D(src,i,j); int data2 = cvGetReal2D(temp,i,j); histogram[data1][data2]++; } } ;i++)//直方图归一化 ;j++) p_histogram[i][j]=(histogram[i][j]*1.0)/(pixel*1.0); Pst_0]; ;i++)//计算概率分布情况 ;j++) { double temp =0.0; ) temp ][j]; ) temp ]; ) temp ]; temp = temp + p_histogram[i][j]; Pst_0[i][j]= temp; } Xst_0]; ;i++)//计算x方向上的均值矢量 ;j++) { double temp =0.0; ) temp ][j]; ) temp ]; ) temp ]; temp = temp + i * p_histogram[i][j]; Xst_0[i][j]= temp; } ];//存储y方向上的均值矢量 Yst_0]; ;i++)//计算y方向上的均值矢量 ;j++) { double temp =0.0; ) temp ][j]; ) temp ]; ) temp ]; temp = temp + j * p_histogram[i][j]; Yst_0[i][j]= temp; } int threshold1; int threshold2; double variance =0.0; double maxvariance =0.0; ;i++)//计算类间离散测度 ;j++) { longdouble p0 = Pst_0[i][j]; ); ]+Pst_0[i][j]; ]+Xst_0[i][j]; ]+Yst_0[i][j]; ); variance = p0*v0+p1*v1; if(variance > maxvariance) { maxvariance = variance; threshold1 = i; threshold2 = j; } } //printf(“%d %d”,threshold1,threshold2); ; } |
总结: 二维otsu算法得到一个阈值,然后对图像做二值化.仍然不能解决光照不均匀二值化的问题.比一维otsu效果好一些,但不是很明显.这个算法的亮点在于考虑的点的附近区域的均值.
