OTSU(大津算法)

转载自https://blog.youkuaiyun.com/piaoxuezhong/article/details/78302893

在图像处理中Otsu方法，是以 Nobuyuki otsu 的名字命名的（日本人，大津展之），常用于基于图像分割的聚类。该算法的理论依据是：假定图像包含两类像素（前景像素和背景像素），直方图为双峰直方图，然后计算使得两类像素能分开的最佳阈值（类内方差），或等价的间类间方差最大。

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)
　　　　　　g=ω0ω1(μ0-μ1)^2 　　 (7)

采用遍历的方法使得类间方差g最大的阈值T，即为所求。Ostu方法可以形象地理解为：求取直方图有两个峰值的图像中那两个峰值之间的低谷值 T 。

Otsu算法实现：

    getThreshVal_Otsu_8u( const Mat& _src )  
    {  
        Size size = _src.size();  
        if( _src.isContinuous() )  
        {  
            size.width *= size.height;  
            size.height = 1;  
        }  
        const int N = 256;  
        int i, j, h[N] = {0};  
        for( i = 0; i < size.height; i++ )  
        {  
            const uchar* src = _src.data + _src.step*i;  
            j = 0;  
            #if CV_ENABLE_UNROLLED  
            for( ; j <= size.width - 4; j += 4 )  
            {  
                int v0 = src[j], v1 = src[j+1];  
                h[v0]++; h[v1]++;  
                v0 = src[j+2]; v1 = src[j+3];  
                h[v0]++; h[v1]++;  
            }  
            #endif  
            for( ; j < size.width; j++ )  
                h[src[j]]++;  
        }  

        double mu = 0, scale = 1./(size.width*size.height);  
        for( i = 0; i < N; i++ )  
            mu += i*(double)h[i];  

        mu *= scale;  
        double mu1 = 0, q1 = 0;  
        double max_sigma = 0, max_val = 0;  

        for( i = 0; i < N; i++ )  
        {  
            double p_i, q2, mu2, sigma;  

            p_i = h[i]*scale;  
            mu1 *= q1;  
            q1 += p_i;  
            q2 = 1. - q1;  

            if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )  
                continue;  

            mu1 = (mu1 + i*p_i)/q1;  
            mu2 = (mu - q1*mu1)/q2;  
            sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);  
            if( sigma > max_sigma )  
            {  
                max_sigma = sigma;  
                max_val = i;  
            }  
        }  

        return max_val;  
    }