大津法理解

最近在看图像分割，大津法在阈值分割上表现还是很不错的，这里记录一下原理和算法实现
参考大津法理解http://blog.163.com/yuyang_tech/blog/static/216050083201302113341762/
自适应阈值算法（大津法）http://blog.csdn.net/guoyk1990/article/details/7606032
简介：

大津法（OTSU）是一种确定图像二值化分割阈值的算法，由日本学者大津于1979年提出。从大津法的原理上来讲，该方法又称作最大类间方差法，因为按照大津法求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大。

原理：

对于图像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)　这就是类间方差
采用遍历的方法得到使类间方差g最大的阈值T,即为所求。

opencv源码实现：

    sdasda

    static double    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;    }