自适应阈值算法（大津阈值法）

                最大类间方差法是由日本学者大津于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。对于图像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,即为所求。

由于，当图像在254或255灰度值上没有像素点时，求平均灰度时会出现0/0的情况，为避免抛出异常，可在当出现前景像素数为零时，跳出循环。

以下是c++代码（用到OpenCV）：

  
  

    
    

     
     

    
    
    
    

     
     

      
      int otsuThreshold(IplImage* img)
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
      {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      int T = 
      
      0;
      
      //阈值
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      int height = img->height;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      int width  = img->width;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      int step      = img->widthStep;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      int channels  = img->nChannels;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       uchar* data  = (uchar*)img->imageData;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double gSum0;
      
      //第一类灰度总值
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double gSum1;
      
      //第二类灰度总值
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double N0 = 
      
      0;
      
      //前景像素数
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double N1 = 
      
      0;
      
      //背景像素数
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double u0 = 
      
      0;
      
      //前景像素平均灰度
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double u1 = 
      
      0;
      
      //背景像素平均灰度
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double w0 = 
      
      0;
      
      //前景像素点数占整幅图像的比例为ω0
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double w1 = 
      
      0;
      
      //背景像素点数占整幅图像的比例为ω1
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double u = 
      
      0;
      
      //总平均灰度
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double tempg = 
      
      -1;
      
      //临时类间方差
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double g = 
      
      -1;
      
      //类间方差
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double Histogram[
      
      256]={
      
      0};
      
      // = new double[256];//灰度直方图
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      double N = width*height;
      
      //总像素数
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      for(
      
      int i=
      
      0;i<height;i++)
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       {
      
      //计算直方图
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        
      
      for(
      
      int j=
      
      0;j<width;j++)
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         
      
      double temp =data[i*step + j * 
      
      3] * 
      
      0.114 + data[i*step + j * 
      
      3+
      
      1] * 
      
      0.587 + data[i*step + j * 
      
      3+
      
      2] * 
      
      0.299;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         temp = temp<
      
      0? 
      
      0:temp;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         temp = temp>
      
      255? 
      
      255:temp;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         Histogram[(
      
      int)temp]++;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        } 
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       }
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      //计算阈值
     
     
    
    

    
    

     
     

    
    
    
    

     
     

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

    
    

     
     

    
    
    
    

     
     

      
       {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        gSum0 = 
      
      0;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        gSum1 = 
      
      0;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        N0 += Histogram[i];   
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        N1 = N-N0;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        
      
      if(
      
      0==N1)
      
      break;
      
      //当出现前景无像素点时，跳出循环
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        w0 = N0/N;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        w1 = 
      
      1-w0;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

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

    
    

     
     

    
    
    
    

     
     

      
        {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         gSum0 += j*Histogram[j];
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        }
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        u0 = gSum0/N0;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

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

    
    

     
     

    
    
    
    

     
     

      
        {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         gSum1 += k*Histogram[k];
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        }
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        u1 = gSum1/N1;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        
      
      //u = w0*u0 + w1*u1;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

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

    
    

     
     

    
    
    
    

     
     

      
        
      
      if (tempg<g)
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        {
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         tempg = g;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
         T = i;
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
        }
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       }
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
       
      
      return T; 
     
     
    
    

    
    

     
     

    
    
    
    

     
     

      
      }