OTSU算法

OTSU算法

OTSU算法以最佳门限将图像灰度直方图分割成两部分，使两部分类间方差取最大值，即分离性最大。

设图像灰度级 ，第 级象素 个,总象素，则第级灰度出现的概率为 。

设灰度门限值为 ，则图像像素按灰度级被分为两类:

，

图像总平均灰度级：

类的平均灰度级为: ，像素数为:

类的平均灰度级为:   , 像素数为:

两部分图像所占比例分别为:

      

对，均值作处理:

     

图像总均值可化为:  

类间方差: 

化为:     

从变化，使 最大的 即为所求之最佳门限。

称为目标选择函数。

 

 

 

 

 

My code：

       int i;          //图像灰度级 from 0－－255

       double n[256];     //第i级像素n[i]个

       double Pro[256];   //第i级灰度出现的概率

       double Tavegray = 0;   //图像总平均灰度级

       double avegray0[256];   //C0类的平均灰度级

       double avegray1[256];   //C1类的平均灰度级

       int N = 0;          //总像素

       int N0 = 0;         //C0类的像素数

       int N1 = 0;         //C1类的像素数

       double w0 = 0, w1 =0;     //分别为两部分图像所占比例

       double mean0 = 0, mean1 =0, mean = 0;  //对C0，C1 均值作处理，及图像总均值

       double var[256];        //类间方差

       int k = 0;      //灰度门限值

       double w[256];     //w[k]当门限值为k时，C0类图像所占的比值

       int Threshold;  //门限值

      

       //数组初始化

       for( i = 0; i <256; i++ )

       {

              n[i] = 0;

              Pro[i] = 0;

              avegray0[i] = 0;

              avegray1[i] = 0;

              var[i] = 0;

              w[i] = 0;

       }

       int I;

      

       LPSTRlpDibHdr = (LPSTR)::GlobalLock(hDib);

       LPBITMAPINFOHEADERlpbmi = (LPBITMAPINFOHEADER)lpDibHdr; //得到指向位图的信息头结构的指针

       LPSTRlpDibBits = lpDibHdr - lpbmi->biSize; //得到指向位图数据的指针

 

       for( int y = 0; y <lpbmi->biHeight; y++ )

       {

              for( int x = 0; x < lpbmi->biWidth; x++ )

              {

                     COLORREFc;

                     BYTE*pDate = (BYTE*)lpDibBits + (lpbmi->biHeight-y) * getpiont( hDib )+ ( x ) *3; //((lpbmi->biWidth*24+31)/32*4)

                     c= RGB( pDate[2], pDate[1], pDate[0]);

              //     COLORREF cNew;

              //     BYTE* pDateNew = (BYTE*)lpNewDibBits +(lpNewbmi->biHeight -y) * getpiont(hDib) + ( x ) * 3; //((lpbmi->biWidth*24+31)/32*4)

       //   cNew = RGB( pDateNew[2], pDateNew[1], pDateNew[0]);

                     I= ( pDate[2] + pDate[1] + pDate[0] ) / 3;

                     n[I]++;

              }

       }

       //总象素

       for( i = 0; i <256; i++ )

       {

              N += n[i];

       }

       //第i级灰度出现的概率

       for( i = 0; i <256; i++ )

       {

              Pro[i] = ( n[i] / N );

       }

       //图像总平均灰度级

       for( i = 1; i <=256; i++ )

       {

              Tavegray += ( i * Pro[i - 1] );

       }

      

  for ( k = 0; k< 256; k++ )

  {

       for( i = 1; i < k+1; i++ )

       {

              avegray0[k] += i*Pro[i];//C0类的平均灰度级

       //     N0 +=n[i];            //C0类像素数

              w0 += Pro[i];

       }

       w[k]= w0;//C0部分图像所占比例

      

//     avegray1[k]= Tavegray - avegray0[k];//C1类的平均灰度级和像素数

//     N1= N - N0;                   //C1类像素数

       w1= 1 - w[k];    //C0部分图像所占比例

       if ( w[k] != 0 )  //对C0均值作处理

       {

       mean0 = avegray0[k]/ w[k]; 

       }

       else mean0 = 0;

           

       mean1= ( Tavegray - avegray0[k] ) / ( 1 - w[k] );//对C1均值作处理

       mean= ( w0 * mean0 ) + ( w1 * mean1 );

       var[k]= (( ( mean * w[k] ) - avegray0[k]) * ( ( mean * w[k] ) - avegray0[k])) / (w[k] * ( 1 - w[k]));

      

}    

 

  int temp = 0;

  for( k = 0; k< 256; k++ )

  {

       if ( var[k] >temp)

       {

              temp = var[k];

              Threshold = k;

       }

}

      

   

 

 

 

OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。下面的代码最早由Ryan Dibble提供，此后经过多人Joerg.Schulenburg,R.Z.Liu 等修改，补正。
 
算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。划分点就是求得的阈值。
  
 parameter:  *image             --- buffer for image
                     rows, cols        --- size of image
                     x0, y0, dx, dy   --- region of vector used for computing threshold
                     vvv                 --- debug option, is 0, no debug information outputed
 */
/*======================================================================*/
/*   OTSU global thresholdingroutine                                                */
/*   takes a 2D unsigned char array pointer, number of rows,and        */
/*   number of cols in the array. returns the value of thethreshold       */
/*======================================================================*/
int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, intdy, int vvv)
{

  unsigned char*np;      // 图像指针
  int thresholdValue=1; // 阈值
  intihist[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 */
    }
  }

  // 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 thresholdingmethod
  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 ourthresholding value

  // debug code to displaythresholding values
  if ( vvv & 1 )
  fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%dgmax=%d\n",
     thresholdValue, gmin, gmax);

  return(thresholdValue);
}

matlab的OTSU算法；

function threshold=ostu(filename);

 

x=imread('filename');

%figure;

%imshow(x);

[m,n]=size(x);

N=m*n;

num=zeros(1,256);

p=zeros(1,256);

 

for i=1:m

for j=1:n

num(x(i,j)+1)=num(x(i,j)+1)+1;

end

end

 

for i=0:255;

p(i+1)=num(i+1)/N;

end

 

totalmean=0;

for i=0:255;

totalmean=totalmean+i*p(i+1);

end

 

maxvar=0;

 

for k=0:255

kk=k+1;

zerosth=sum(p(1:kk));

 

firsth=0;

for h=0:k

firsth=firsth+h*p(h+1);

end

 

var=totalmean*zerosth-firsth;

var=var*var;

var=var/(zerosth*(1-zerosth)+0.01);

var=sqrt(var);

if(var>maxvar)

maxvar=var;

point=k;

end

 

end

 

threshold=point;

 

网址：http://wenku.baidu.com/link?url=tOgi4aWIUlfxOSYWzXFRQnicPsTGDa048Yj9g_jySVmInlwnyhdCJp6qFVCQtYAqscDnnu7G4PCu-NDrw8uJ08Kt34gEptNX3GmcCWE-gE7