C++实现otsu算法

        类间方差最大化阈值分割算法 (Otsu)可以自动计算图像二值化时的阈值，otsu的中心思想是阈值T应使目标与背景两类的类间方差最大。具体原理如下：

        otsu是按图像的灰度特性，将图像分成背景和目标两部分。背景和目标之间的类间方差越大，说明构成图像的两部分的差别越大，当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此，使类间方差最大的分割意味着错分概率最小。对于图像I(x,y)，前景(即目标)和背景的分割阈值记作T，属于前景的像素点数占整幅图像的比例记为ω0，其平均灰度μ0；背景像素点数占整幅图像的比例为ω1，其平均灰度为μ1。图像的总平均灰度记为μ，类间方差记为g。假设图像的背景较暗，并且图像的大小为M×N，图像中像素的灰度值小于阈值T的像素个数记作N0，像素灰度大于阈值T的像素个数记作N1，则有：
　　　　　　ω0=N0/M×N          (1)
　　　　　　ω1=N1/M×N          (2)
　　　　　　N0+N1=M×N        (3)
　　　　　　ω0+ω1=1              (4)
　　　　　　μ=ω0*μ0+ω1*μ1   (5)
　　　　　　g=ω0(μ0-μ)^2+ω1(μ1-μ)^2    (6)
         将式(5)代入式(6),得到等价公式：
                        g=ω0ω1(μ0-μ1)^2     (7)
        采用遍历的方法得到使类间方差最大的阈值T，即为所求。

 

         Otsu算法的C++实现代码如下，传入参数为图像数据，返回值为二值化阈值。

int thresHoldOfOtsu(unsigned char img[HEIGHT][WIDTH])

{

       int i = 0;

       int j = 0;

       intthresHold = 0;//最佳二值化阈值

       unsignedint hist[256];//用来统计各个灰度值的像素个数

       intsumHist = WIDTH*HEIGHT;//总的像素个数

       float g =0.0;

       float maxg= 0.0;

       float u0 =0.0;

       float u1 =0.0;

       float w0 =0.0;

       float w1 =0.0;

       unsignedint N0 = 0;

       unsignedint N1 = 0;

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

       {

              hist[i]= 0;

       }

       //统计各个灰度值的像素个数

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

       {

              for( j = 0; j < WIDTH; j++)

              {

                     hist[img[i][j]]++;

              }

       }

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

       {

              N0= 0;

              //计算N0、u0、w0

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

              {

                     N0+= hist[j];

              }

              w0= N0*1.0/sumHist;

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

              {

                     if(N0 == 0)

                     {

                            u0= 0;

                     }

                     else

                     {

                            u0+= hist[j]*j*1.0/N0;

                     }

              }

              //计算N1、u1、w1

              N1= sumHist - N0;

              w1= 1 - w0;

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

              {

                     if(N1 == 0)

                     {

                            u1= 0;

                     }

                     else

                     {

                            u1+= hist[j]*j*1.0/N1;

                     }

              }

              g =w0*w1*(u0-u1)*(u0-u1);

              if(g>maxg)

              {

                     maxg= g;

                     thresHold= i;

              }

       }

       returnthresHold;

}