图像二值化方法--OTSU（最大类间方差法）

来源：互联网发布：农村淘宝 app 下载编辑：程序博客网时间：2024/05/22 00:29

前面学习了直方图双峰法：图像二值化方法中的阈值法

最大类间方差法（OTSU）是找到自适应阈值的常用方法。原理参考了冈萨雷斯的《数字图像处理》。

以下是自己写的函数：

//----获取灰度图in的OTSU阈值--int Segment::otsuMat(Mat in){int i,j;int temp;//第一类均值，第二类均值，全局均值，mk=p1*m1, 第一类概率，第二类概率double m1,m2,mG,mk,p1,p2;int hist[256] = {0};double pro_hist[256] = {0.0};double cov;double maxcov=0.0;int maxthread=0;int row = in.rows;int col = in.cols;//统计每个灰度的数量for(i=0; i<row; ++i) {for(j=0; j<col; ++j) {temp = in.at<uchar>(i,j);hist[temp]++;}}//计算每个灰度级占图像的概率for(i=0; i<256; ++i)pro_hist[i] = (double)hist[i]/(double)(row*col);//计算平均灰度值mG = 0.0;for(i=0; i<256; ++i)mG += i * pro_hist[i];//统计前景和背景的平均灰度值，并计算类间方差for(i=0; i<256; ++i){m1=0.0; m2=0.0; mk=0.0; p1=0.0; p2=0.0;for(j=0; j<i; ++j) {p1 += pro_hist[j];mk += j*pro_hist[j];}m1 = mk/p1;  //mk=p1*m1,是一个中间值p2 = 1 - p1;  //p1+p2=1;m2 = (mG - mk)/p2;  //mG=p1*m1+p2*m2;//计算类间方差cov = p1*p2*(m1-m2)*(m1-m2);if(cov>maxcov) {maxcov = cov;maxthread = i;}}//cout<<maxthread<<endl;return maxthread;}

每次需要的时候就直接复制粘贴，效果常常很满意，就感觉要二值化OTSU法总可以用上，但是今天发现不好用了。我需要把一组图片分离出目标和背景，选取了Lab空间中的L通道，直接使用OTSU法。本来这一组图背景很相似，但是结果发现分割后有一些图片错误太明显。查看了阈值，发现相差蛮多的。例如下面的四张图。

灰度图：

对应的二值图：

图1,2,4的结果都可以接受，但是图3的结果突然就差了很远，即使我在OTSU法算得的阈值基础上加加减减也没办法兼顾到每一张图。

我列出了灰度图的均值和方差，以及OTSU算得的阈值

这四张图的均值和方差很接近，但是OTSU算得的阈值却变化较大，尤其是图3，导致二值化结果不正确。

觉得自己对“最大类间方差”理解有误，至少，它不能解决这一类问题。它的结果 out of my control. 琢磨琢磨，到底是怎么回事？

最后通过均值和标准差设置了阈值，效果可控。

另外，要使用OTSU，也可以直接采用opencv 库函数：

double threshold(InputArray src, OutputArray dst, double thresh, double maxval, int type)

当 type = CV_THRESH_OTSU 时, thresh 的值会被忽略，所以可以任取一个值。

在手册上，type 中没有列出OTSU方法，但是在函数源代码中是有的：

double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double maxval, int type ){    Mat src = _src.getMat();    bool use_otsu = (type & THRESH_OTSU) != 0;    type &= THRESH_MASK;    if( use_otsu )    {        CV_Assert( src.type() == CV_8UC1 );        thresh = getThreshVal_Otsu_8u(src);    }    _dst.create( src.size(), src.type() );    Mat dst = _dst.getMat();    if( src.depth() == CV_8U )    {        int ithresh = cvFloor(thresh);        thresh = ithresh;        int imaxval = cvRound(maxval);        if( type == THRESH_TRUNC )            imaxval = ithresh;        imaxval = saturate_cast<uchar>(imaxval);        if( ithresh < 0 || ithresh >= 255 )        {            if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||                ((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < 0) ||                (type == THRESH_TOZERO && ithresh >= 255) )            {                int v = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :                        type == THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :                        /*type == THRESH_TRUNC ? imaxval :*/ 0;                dst.setTo(v);            }            else                src.copyTo(dst);            return thresh;        }        thresh = ithresh;        maxval = imaxval;    }    else if( src.depth() == CV_16S )    {        int ithresh = cvFloor(thresh);        thresh = ithresh;        int imaxval = cvRound(maxval);        if( type == THRESH_TRUNC )            imaxval = ithresh;        imaxval = saturate_cast<short>(imaxval);        if( ithresh < SHRT_MIN || ithresh >= SHRT_MAX )        {            if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||               ((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < SHRT_MIN) ||               (type == THRESH_TOZERO && ithresh >= SHRT_MAX) )            {                int v = type == THRESH_BINARY ? (ithresh >= SHRT_MAX ? 0 : imaxval) :                type == THRESH_BINARY_INV ? (ithresh >= SHRT_MAX ? imaxval : 0) :                /*type == THRESH_TRUNC ? imaxval :*/ 0;                dst.setTo(v);            }            else                src.copyTo(dst);            return thresh;        }        thresh = ithresh;        maxval = imaxval;    }    else if( src.depth() == CV_32F )        ;    else        CV_Error( CV_StsUnsupportedFormat, "" );    parallel_for_(Range(0, dst.rows),                  ThresholdRunner(src, dst, thresh, maxval, type),                  dst.total()/(double)(1<<16));    return thresh;}

其中 type 取值：

/* Threshold types */enum{    CV_THRESH_BINARY      =0,  /* value = value > threshold ? max_value : 0       */    CV_THRESH_BINARY_INV  =1,  /* value = value > threshold ? 0 : max_value       */    CV_THRESH_TRUNC       =2,  /* value = value > threshold ? threshold : value   */    CV_THRESH_TOZERO      =3,  /* value = value > threshold ? value : 0           */    CV_THRESH_TOZERO_INV  =4,  /* value = value > threshold ? 0 : value           */    CV_THRESH_MASK        =7,    CV_THRESH_OTSU        =8  /* use Otsu algorithm to choose the optimal threshold value;                                 combine the flag with one of the above CV_THRESH_* values */};

其中 getThreshVal_Otsu_8u():

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

1 0