灰度图像阈值化分割常见方法总结及VC实现

转载地址：http://blog.csdn.net/likezhaobin/article/details/6915755

在图像处理领域，二值图像运算量小，并且能够体现图像的关键特征，因此被广泛使用。将灰度图像变为二值图像的常用方法是选定阈值，然后将待处理图像的每个像素点进行单点处理，即将其灰度值与所设置的门限进行比对，从而得到二值化的黑白图。这样一种方式因为其直观性以及易于实现，已经在图像分割领域处于中心地位。本文主要对最近一段时间作者所学习的阈值化图像分割算法进行总结，全文描述了作者对每种算法的理解，并基于OpenCV和VC6.0对这些算法进行了实现。最终将源代码公开，希望大家一起进步。（本文的代码暂时没有考虑执行效率问题）

      首先给出待分割的图像如下：

1、Otsu法（最大类间方差法）

       该算法是日本人Otsu提出的一种动态阈值分割算法。它的主要思想是按照灰度特性将图像划分为背景和目标2部分，划分依据为选取门限值，使得背景和目标之间的方差最大。（背景和目标之间的类间方差越大，说明这两部分的差别越大，当部分目标被错划分为背景或部分背景错划分为目标都会导致这两部分差别变小。因此，使用类间方差最大的分割意味着错分概率最小。）这是该方法的主要思路。其主要的实现原理为如下：

      1）建立图像灰度直方图（共有L个灰度级，每个出现概率为p）

      2）计算背景和目标的出现概率，计算方法如下：

      上式中假设t为所选定的阈值，A代表背景（灰度级为0~N）,根据直方图中的元素可知，Pa为背景出现的概率，同理B为目标，Pb为目标出现的概率。

      3）计算A和B两个区域的类间方差如下：

        第一个表达式分别计算A和B区域的平均灰度值；

      第二个表达式计算灰度图像全局的灰度平均值；

      第三个表达式计算A、B两个区域的类间方差。

      4）以上几个步骤计算出了单个灰度值上的类间方差，因此最佳分割门限值应该是图像中能够使得A与B的类间灰度方差最大的灰度值。在程序中需要对每个出现的灰度值据此进行寻优。

      本人的VC实现代码如下。

/*****************************************************************************
*
* 函数名称：
*   OneDimentionOtsu()
*
* 输入参数:
*   pGrayMat:      二值图像数据
*            图形尺寸宽度
*   height:        图形尺寸高度
*   nTlreshold:    经过算法处理得到的二值化分割阈值
* 返回值:
*   无
* 函数说明:实现灰度图的二值化分割——最大类间方差法（Otsu算法，俗称大津算法）
*
****************************************************************************/

void CBinarizationDlg::OneDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)
{
    double nHistogram[256];         //灰度直方图
    double dVariance[256];          //类间方差
    int N = height*width;           //总像素数
    for(int i=0; i<256; i++)
    {
        nHistogram[i] = 0.0;
        dVariance[i] = 0.0;
    }
    for(i=0; i<height; i++)
    {
        for(int j=0; j<width; j++)
        {
            unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);
            nHistogram[nData]++;     //建立直方图
        }
    }
    double Pa=0.0;      //背景出现概率
    double Pb=0.0;      //目标出现概率
    double Wa=0.0;      //背景平均灰度值
    double Wb=0.0;      //目标平均灰度值
    double W0=0.0;      //全局平均灰度值
    double dData1=0.0,  dData2=0.0;
    for(i=0; i<256; i++)     //计算全局平均灰度
    {
        nHistogram[i] /= N;
        W0 += i*nHistogram[i];
    }
    for(i=0; i<256; i++)     //对每个灰度值计算类间方差
    {
        Pa += nHistogram[i];
        Pb = 1-Pa;
        dData1 += i*nHistogram[i];
        dData2 = W0-dData1;
        Wa = dData1/Pa;
        Wb = dData2/Pb;
        dVariance[i] = (Pa*Pb* pow((Wb-Wa), 2));
    }
    //遍历每个方差，求取类间最大方差所对应的灰度值
    double temp=0.0;
    for(i=0; i<256; i++)
    {
        if(dVariance[i]>temp)
        {
            temp = dVariance[i];
            nThreshold = i;
        }
    }
}

      阈值分割结果如下图，求解所得的阈值为116.

2、一维交叉熵值法

      这种方法与类间最大方差很相似，是由Li和Lee应用了信息论中熵理论发展而来。首先简要介绍交叉熵的概念。

      对于两个分布P和Q，定义其信息交叉熵D如下：

      这代表的物理意义是两个分布之间信息理论距离，另外一种理解是，将分布P变为Q后所带来的信息变化。那么对于图像分割来说，如果要用分割图像来替换原来的图像，最优的分割依据应该就是使得两幅图像之间的交叉熵最小。以下对最小交叉熵法的过程进行简要总结。

      可以假设上文的P为源图像的灰度分布，Q为所得到的分割图像的灰度分布，其中：

      

      上式中H为统计直方图；

      N为图像总的像素点数；

      L为源图像总的灰度级数；

      P代表源图像，其每个元素代表每个灰度级上的灰度分布（平均灰度值）；

      Q为分割后的二值图像，两个u分别代表两个分割后的区域的平均灰度值，其中t为分割图像所采用的阈值。

      根据以上定义，以每个灰度级上的灰度和为计算量，可以很容易根据交叉熵的公式，推导出P和Q之间的交叉熵定量表达式：

      根据上文所述思路，使得D最小的t即为最小交叉熵意义下的最优阈值。

      作者VC实现代码如下。

/*****************************************************************************
*
* 函数名称：
*   MiniCross()
*
* 输入参数:
*   pGrayMat:      二值图像数据
*            图形尺寸宽度
*   height:        图形尺寸高度
*   nTlreshold:    经过算法处理得到的二值化分割阈值
* 返回值:
*   无
* 函数说明:实现灰度图的二值化分割——最小交叉熵算法
*
****************************************************************************/

void CBinarizationDlg::MiniCross(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)
{
    double dHistogram[256];         //灰度直方图
    double dEntropy[256];           //每个像素的交叉熵
    int N = height*width;           //总像素数
    for(int i=0; i<256; i++)
    {
        dHistogram[i] = 0.0;
        dEntropy[i] = 0.0;
    }
    for(i=0; i<height; i++)
    {
        for(int j=0; j<width; j++)
        {
            unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);
            dHistogram[nData]++;     //建立直方图
        }
    }
    double Pa=0.0;      //区域1平均灰度值
    double Pb=0.0;      //区域2平均灰度值
    double P0=0.0;      //全局平均灰度值
    double Wa=0.0;      //第一部分熵
    double Wb=0.0;      //第二部分的熵
    double dData1=0.0, dData2=0.0;  //中间值
    double dData3=0.0, dData4=0.0;  //中间值

    for(i=0; i<256; i++)     //计算全局平均灰度
    {
        dHistogram[i] /= N;
        P0 += i*dHistogram[i];
    }

    for(i=0; i<256; i++)
    {
        Wa=Wb=dData1=dData2=dData3=dData4=Pa=Pb=0.0;
        for(int j=0; j<256; j++)
        {
            if(j<=i)
            {
                dData1 += dHistogram[j];
                dData2 += j*dHistogram[j];
            }
            else
            {
                dData3 += dHistogram[j];
                dData4 += j*dHistogram[j];
            }
        }
        Pa = dData2/dData1;
        Pb = dData4/dData3;
        for(j=0; j<256; j++)
        {
            if(j<=i)
            {
                if((Pa!=0)&&(dHistogram[j]!=0))
                {
                    double d1 = log(dHistogram[j]/Pa);
                    Wa += j*dHistogram[j]*d1/log(2);
                }
            }
            else
            {
                if((Pb!=0)&&(dHistogram[j]!=0))
                {
                    double d2 = log(dHistogram[j]/Pb);
                    Wb += j*dHistogram[j]*d2/log(2);
                }
            }
        }
        dEntropy[i] = Wa+Wb;
    }
    //遍历熵值，求取最小交叉熵所对应的灰度值
    double temp=dEntropy[0];
    for(i=1; i<256; i++)
    {
        if(dEntropy[i]<temp)
        {
            temp = dEntropy[i];
            nThreshold = i;
        }
    }
}

      阈值分割结果如下图，求解所得的阈值为106.

3、二维OTSU法

           这种方法是对类间最大方差法的扩展，将其从求两个一维分布最大类间方差扩充为求解类间离散度矩阵的迹的最大值，考虑像素点灰度级的基础上增加了对像素点邻域平均像素值的考虑。

      以下按照本人的理解对该方法的思路以及推倒过程进行分析：

      1）首先需要建立二维的灰度统计直方图P(f, g)；

      图像的灰度级为L级，那么其每个像素点的8邻域灰度平均值的灰度级也为L级，据此来构建直方图P。二维统计直方图的横轴为每个像素点的灰度值f(i, j)，纵坐标为同一个点对应的邻域平均值g(i, j) 其中(0≤i<height, 0≤j<width, 0≤f(i, j)<L)，而所对应的P(f,g)整幅图像中灰度值为f，邻域灰度均值为g的点的统计值占总像素点的比例（即为灰度值出现的联合概率密度）。其中P的每个元素满足如下公式：

      n为整幅图像中灰度值为f，邻域灰度均值为g的点的统计值；

      N为图像总的像素点个数；

      2）对于下图所示的二维统计直方图，t代表横坐标（灰度值），s代表纵坐标（像素点邻域的灰度均值）

      对已图像中的阈值点（t,s）来说，其灰度值t和其邻域内的灰度均值s不应该相差太多，如果t比s大很多（点位于上图中的II区域），说明像素的灰度值远远大于其临域的灰度均值，故而该点很可能是噪声点，反之如果t比s小很多（点位于途中的IV区域），即该点的像素值比其临域均值小很多，则说明是一个边缘点。据此我们在进行背景前景分割的时候忽略这些干扰因素，认为这两个区域内Pi,j=0。剩下的I区域和III区域则分别代表了前景和背景。以下据此来推导对于选定的阈值（t, s），进行离散度判据的最优推导。

      3）推导阈值（t, s）点处的离散度矩阵判据

      根据上文分析可知，由阈值（t, s）所分割的前景和背景出现的概率如下：

      定义两个中间变量，方便下面推导：

      据此，这两部分的灰度均值向量可以推导如下（两个分量分别根据灰度值以及每个点的灰度均值计算）：

      整幅图像的灰度均值向量为：

      与一维的大津法一样的思路，推导类间方差，这里是二维因此要用矩阵形式。参考一维法，可同样定义类间“方差”矩阵如下：

      为了在实现的时候，容易判断出这样一个矩阵的“最大值”，因此数学中采用矩阵的迹（对角线之和）来衡量矩阵的“大小”。因此以该矩阵的迹作为离散度测度，推导如下：

      这样就可以通过求解使得这个参数最大时的（t,s）即为所求得的最佳阈值组合。

      以下为具体算法实现过程：

      1）建立二维直方图

      2）对直方图进行遍历，计算每个（t,s）组合所得到的矩阵离散度，也就是一维大津法中所谓的最大类间方差。

      3）求得使“类间方差”最大的（t,s），由于t代表灰度值，s代表改点在其邻域内的灰度均值，因此本人认为在选择阈值时可以选择s为最佳，当然用t也可以，因为从求解结果可以看出，这这个数值往往很接近。

      具体的实现代码如下：

/*****************************************************************************
*
* 函数名称：
*   TwoDimentionOtsu()
*
* 输入参数:
*   pGrayMat:      二值图像数据
*            图形尺寸宽度
*   height:        图形尺寸高度
*   nTlreshold:    经过算法处理得到的二值化分割阈值
* 返回值:
*   无
* 函数说明:实现灰度图的二值化分割——最大类间方差法（二维Otsu算法）
* 备注：在构建二维直方图的时候，采用灰度点的3*3邻域均值
******************************************************************************/
void CBinarizationDlg::TwoDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)
{
    double dHistogram[256][256];        //建立二维灰度直方图
    double dTrMatrix = 0.0;             //离散矩阵的迹
    int N = height*width;               //总像素数
    for(int i=0; i<256; i++)
    {
        for(int j=0; j<256; j++)
            dHistogram[i][j] = 0.0;      //初始化变量
    }
    for(i=0; i<height; i++)
    {
        for(int j=0; j<width; j++)
        {
            unsigned char nData1 = (unsigned char)cvmGet(pGrayMat, i, j);  //当前的灰度值
            unsigned char nData2 = 0;
            int nData3 = 0;         //注意9个值相加可能超过一个字节
            for(int m=i-1; m<=i+1; m++)
            {
                for(int n=j-1; n<=j+1; n++)
                {
                    if((m>=0)&&(m<height)&&(n>=0)&&(n<width))
                        nData3 += (unsigned char)cvmGet(pGrayMat, m, n); //当前的灰度值
                }
            }
            nData2 = (unsigned char)(nData3/9);    //对于越界的索引值进行补零,邻域均值
            dHistogram[nData1][nData2]++;
        }
    }
    for(i=0; i<256; i++)
        for(int j=0; j<256; j++)
            dHistogram[i][j] /= N;  //得到归一化的概率分布

    double Pai = 0.0;      //目标区均值矢量i分量
    double Paj = 0.0;      //目标区均值矢量j分量
    double Pbi = 0.0;      //背景区均值矢量i分量
    double Pbj = 0.0;      //背景区均值矢量j分量
    double Pti = 0.0;      //全局均值矢量i分量
    double Ptj = 0.0;      //全局均值矢量j分量
    double W0 = 0.0;       //目标区的联合概率密度
    double W1 = 0.0;       //背景区的联合概率密度
    double dData1 = 0.0;
    double dData2 = 0.0;
    double dData3 = 0.0;
    double dData4 = 0.0;   //中间变量
    int nThreshold_s = 0;
    int nThreshold_t = 0;
    double temp = 0.0;     //寻求最大值
    for(i=0; i<256; i++)
    {
        for(int j=0; j<256; j++)
        {
            Pti += i*dHistogram[i][j];
            Ptj += j*dHistogram[i][j];
        }
    }
    for(i=0; i<256; i++)
    {
        for(int j=0; j<256; j++)
        {
            W0 += dHistogram[i][j];
            dData1 += i*dHistogram[i][j];
            dData2 += j*dHistogram[i][j];

            W1 = 1-W0;
            dData3 = Pti-dData1;
            dData4 = Ptj-dData2;
/*          W1=dData3=dData4=0.0;   //对内循环的数据进行初始化
            for(int s=i+1; s<256; s++)
            {
                for(int t=j+1; t<256; t++)
                {
                    W1 += dHistogram[s][t];
                    dData3 += s*dHistogram[s][t];  //方法2
                    dData4 += t*dHistogram[s][t];  //也可以增加循环进行计算
                }
            }*/

            Pai = dData1/W0;
            Paj = dData2/W0;
            Pbi = dData3/W1;
            Pbj = dData4/W1;   // 得到两个均值向量，用4个分量表示
            dTrMatrix = ((W0*Pti-dData1)*(W0*Pti-dData1)+(W0*Ptj-dData1)*(W0*Ptj-dData2))/(W0*W1);
            if(dTrMatrix > temp)
            {
                temp = dTrMatrix;
                nThreshold_s = i;
                nThreshold_t = j;
            }
        }
    }
    nThreshold = nThreshold_t;   //返回结果中的灰度值
    //nThreshold = 100;
}

      阈值分割结果如下图，求解所得的阈值为114.（s=114，t=117）