OTSU算法学习 OTSU公式证明

转自http://www.cnblogs.com/guopengfei/p/4759569.html

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; }

 

参考了这篇文章 http://baike.baidu.com/link?url=ZgwdvFvH-oZhpE3panvp9kv0p7dn7KSnpc87v-AIBg5najnR4cVmIXwP_A_4nry7USDUZUuEYa-c5P09XOYoIa

 

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效果好一些,但不是很明显.这个算法的亮点在于考虑的点的附近区域的均值.