Otsu算法

Otsu算法：最大类间方差法（大津算法），是一种确定阈值的算法。

之所以称为最大类间方差法是因为，用该阈值进行的图像固定阈值二值化，类间方差最大，它是按图像的灰度特性,将图像分成背景和前景两部分，使类间方差最大的分割意味着错分概率最小。

原理：

对于图像I(x,y)，前景(即目标)和背景的分割阈值记作T，属于前景的像素点数占整幅图像的比例记为ω₀，其平均灰度μ₀；背景像素点数占整幅图像的比例为ω₁，其平均灰度为μ₁。图像的总平均灰度记为μ，类间方差记为g。
假设图像的背景较暗，并且图像的大小为M×N，图像中像素的灰度值小于阈值T的像素个数记作N₀，像素灰度大于阈值T的像素个数记作N₁，则有：
　　　　　　(1) ω₀=N₀/ (M×N)
　　　　　　(2) ω₁=N₁/ (M×N)
　　　　　　(3) N₀ + N₁ = M×N
　　　　　　(4) ω₀ + ω₁ = 1
　　　　　　(5) μ = ω₀ * μ₀ + ω₁ * μ₁
　　　　　　(6) g = ω₀ * (μ₀ - μ)² + ω₁ * (μ₁ - μ)²
将式(5)代入式(6)，得到等价公式：
　　　　　　(7) g = ω₀ *ω₁ * (μ₀ - μ₁)²
采用遍历的方法得到使类间方差g最大的阈值T。

算法评价：

优点：算法简单，当目标与背景的面积相差不大时，能够有效地对图像进行分割。
缺点：当图像中的目标与背景的面积相差很大时，表现为直方图没有明显的双峰，或者两个峰的大小相差很大，分割效果不佳，或者目标与背景的灰度有较大的重叠时也不能准确的将目标与背景分开。
原因：该方法忽略了图像的空间信息，同时将图像的灰度分布作为分割图像的依据，对噪声也相当敏感。

程序：

opencv可以直接调用这种算法，threshold(gray, dst, 0, 255, CV_THRESH_OTSU);

C++

int otsu(IplImage* src)  
{   
    int hist_size = 256;    //直方图尺寸    
    int hist_height = 256;    
    float range[] = {0,255};  //灰度级的范围    
    float* ranges[]={range};    
    //创建一维直方图，统计图像在[0 255]像素的均匀分布    
    CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);    
    //计算灰度图像的一维直方图    
    cvCalcHist(&src,gray_hist,0,0);    
    //归一化直方图    
    cvNormalizeHist(gray_hist,1.0);    
      
    int Width = src->width;  
    int Height = src->height;   
    int threshold = 0;  
    double delta = 0;
    double u = 0;
    for(int m=0; m<256; m++)
    {
        u += cvQueryHistValue_1D(gray_hist,m)*m;  // cvGetReal1D 均值
    }

    
    double v_0 = 0, w_0 = 0; 
    for(int k=0; k<256; k++)  
    {  
            
        v_0 += cvQueryHistValue_1D(gray_hist,k)*k;   // cvGetReal1D  v_0 = w_0 * u_0
        w_0 += cvQueryHistValue_1D(gray_hist,k);      //灰度大于阈值k的像素的概率  

        double t = u * w_0 - v_0;
        double delta_tmp = t * t / (w * (1 - w) ); 
        
        if(delta_tmp > delta)  
        {  
            delta = delta_tmp;  
            threshold = k;  
        }  
    } 

    cvReleaseHist(&gray_hist);
    return threshold;   
} 

Python

import numpy as np

def otsu_threshold(im):
    width, height = im.size
    img_size = width * height
    pixel_counts = np.zeros(256)
    for x in range(width):
        for y in range(height):
            pixel = im.getpixel((x, y))
            pixel_counts[pixel] = pixel_counts[pixel] + 1
    # 得到图片的以0-255索引的像素值个数列表
    s_max = (0, -10)
    for threshold in range(256):
        # 遍历所有阈值
        # 更新
        n_0 = sum(pixel_counts[:threshold])  # 得到阈值以下像素个数
        n_1 = sum(pixel_counts[threshold:])   # 得到阈值以上像素个数
        
        w_0 = n_0 / img_size
        w_1 = n_1 / img_size
        # 得到阈值下所有像素的平均灰度
        u_0 = sum([i * pixel_counts[i] for i in range(0, threshold)]) / n_0 if n_0 > 0 else 0
        
        # 得到阈值上所有像素的平均灰度
        u_1 = sum([i * pixel_counts[i] for i in range(threshold, 256)]) / n_1 if n_1 > 0 else 0
        
        # 总平均灰度
        u = w_0 * u_0 + w_1 * u_1

        # 类间方差
        g = w_0 * (u_0 - u) * (u_0 - u) + w_1 * (u_1 - u) * (u_1 - u)

        # 类间方差等价公式
        # g = w_0 * w_1 * (u_0 * u_1) * (u_0 * u_1)

        # 取最大的
        if g > s_max[1]:
            s_max = (threshold, g)
    return s_max[0]