灰度阈值分割算法.docx

分享到一键分享QQ空间新浪微博百度云收藏人人网腾讯微博百度相册开心网腾讯朋友百度贴吧豆瓣网搜狐微博百度新首页QQ好友和讯微博更多.百度分享登录|注册窗体顶端窗体底端收藏成功 确定收藏失败，请重新收藏 确定窗体顶端标题 标题不能为空网址 标签 摘要 公开 取消收藏 窗体底端分享资讯传PPT/文档提问题写博客传资源创建项目创建代码片设置昵称编辑自我介绍，让更多人了解你帐号设置退出社区博客论坛下载Share极客头条服务CODE英雄会活动CSTO俱乐部CTO俱乐部高校俱乐部奋斗斌斌的专栏解决项目中的琐碎细节问题b.zhao_目录视图摘要视图订阅有奖征资源，博文分享有内涵 人气博主的资源共享：老罗的Android之旅 微软Azure英雄会编程大赛题 关注CSDN社区微信，福利多多 社区问答：叶劲峰 游戏引擎架构 灰度图像阈值化分割常见方法总结及VC实现 分类： 图像处理 OpenCV 2011-11-11 23:20 7427人阅读 评论(14) 收藏 举报 算法图形byte图像处理扩展目录(?)+Otsu法最大类间方差法一维交叉熵值法二维OTSU法参考文献 在图像处理领域，二值图像运算量小，并且能够体现图像的关键特征，因此被广泛使用。将灰度图像变为二值图像的常用方法是选定阈值，然后将待处理图像的每个像素点进行单点处理，即将其灰度值与所设置的门限进行比对，从而得到二值化的黑白图。这样一种方式因为其直观性以及易于实现，已经在图像分割领域处于中心地位。本文主要对最近一段时间作者所学习的阈值化图像分割算法进行总结，全文描述了作者对每种算法的理解，并基于OpenCV和VC6.0对这些算法进行了实现。最终将源代码公开，希望大家一起进步。（本文的代码暂时没有考虑执行效率问题） 首先给出待分割的图像如下：1、Otsu法（最大类间方差法） 该算法是日本人Otsu提出的一种动态阈值分割算法。它的主要思想是按照灰度特性将图像划分为背景和目标2部分，划分依据为选取门限值，使得背景和目标之间的方差最大。（背景和目标之间的类间方差越大，说明这两部分的差别越大，当部分目标被错划分为背景或部分背景错划分为目标都会导致这两部分差别变小。因此，使用类间方差最大的分割意味着错分概率最小。）这是该方法的主要思路。其主要的实现原理为如下： 1）建立图像灰度直方图（共有L个灰度级，每个出现概率为p） 2）计算背景和目标的出现概率，计算方法如下： 上式中假设t为所选定的阈值，A代表背景（灰度级为0N）,根据直方图中的元素可知，Pa为背景出现的概率，同理B为目标，Pb为目标出现的概率。 3）计算A和B两个区域的类间方差如下： 第一个表达式分别计算A和B区域的平均灰度值； 第二个表达式计算灰度图像全局的灰度平均值； 第三个表达式计算A、B两个区域的类间方差。 4）以上几个步骤计算出了单个灰度值上的类间方差，因此最佳分割门限值应该是图像中能够使得A与B的类间灰度方差最大的灰度值。在程序中需要对每个出现的灰度值据此进行寻优。 本人的VC实现代码如下。cpp view plaincopyprint?/*函数名称：*OneDimentionOtsu()*输入参数:*pGrayMat:二值图像数据*width:图形尺寸宽度*height:图形尺寸高度*nTlreshold:经过算法处理得到的二值化分割阈值*返回值:*无*函数说明:实现灰度图的二值化分割最大类间方差法（Otsu算法，俗称大津算法）*/voidCBinarizationDlg:OneDimentionOtsu(CvMat*pGrayMat,intwidth,intheight,BYTE&nThreshold)doublenHistogram256;/灰度直方图 doubledVariance256;/类间方差 intN=height*width;/总像素数 for(inti=0;i256;i+)nHistogrami=0.0;dVariancei=0.0;for(i=0;iheight;i+)for(intj=0;jwidth;j+)unsignedcharnData=(unsignedchar)cvmGet(pGrayMat,i,j);nHistogramnData+;/建立直方图 doublePa=0.0;/背景出现概率 doublePb=0.0;/目标出现概率 doubleWa=0.0;/背景平均灰度值 doubleWb=0.0;/目标平均灰度值 doubleW0=0.0;/全局平均灰度值 doubledData1=0.0,dData2=0.0;for(i=0;i256;i+)/计算全局平均灰度 nHistogrami/=N;W0+=i*nHistogrami;for(i=0;i256;i+)/对每个灰度值计算类间方差 Pa+=nHistogrami;Pb=1-Pa;dData1+=i*nHistogrami;dData2=W0-dData1;Wa=dData1/Pa;Wb=dData2/Pb;dVariancei=(Pa*Pb*pow(Wb-Wa),2);/遍历每个方差，求取类间最大方差所对应的灰度值 doubletemp=0.0;for(i=0;itemp)temp=dVariancei;nThreshold=i;/* 函数名称：* OneDimentionOtsu()* 输入参数:* pGrayMat: 二值图像数据*width: 图形尺寸宽度*height: 图形尺寸高度*nTlreshold: 经过算法处理得到的二值化分割阈值* 返回值:* 无* 函数说明:实现灰度图的二值化分割最大类间方差法（Otsu算法，俗称大津算法）*/void CBinarizationDlg:OneDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)double nHistogram256; /灰度直方图double dVariance256; /类间方差int N = height*width; /总像素数for(int i=0; i256; i+)nHistogrami = 0.0;dVariancei = 0.0;for(i=0; iheight; i+)for(int j=0; jwidth; j+)unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);nHistogramnData+; /建立直方图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; i256; i+) /计算全局平均灰度nHistogrami /= N;W0 += i*nHistogrami;for(i=0; i256; i+) /对每个灰度值计算类间方差Pa += nHistogrami;Pb = 1-Pa;dData1 += i*nHistogrami;dData2 = W0-dData1;Wa = dData1/Pa;Wb = dData2/Pb;dVariancei = (Pa*Pb* pow(Wb-Wa), 2);/遍历每个方差，求取类间最大方差所对应的灰度值double temp=0.0;for(i=0; itemp)temp = dVariancei;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实现代码如下。 cpp view plaincopyprint?/*函数名称：*MiniCross()*输入参数:*pGrayMat:二值图像数据*width:图形尺寸宽度*height:图形尺寸高度*nTlreshold:经过算法处理得到的二值化分割阈值*返回值:*无*函数说明:实现灰度图的二值化分割最小交叉熵算法*/voidCBinarizationDlg:MiniCross(CvMat*pGrayMat,intwidth,intheight,BYTE&nThreshold)doubledHistogram256;/灰度直方图 doubledEntropy256;/每个像素的交叉熵 intN=height*width;/总像素数 for(inti=0;i256;i+)dHistogrami=0.0;dEntropyi=0.0;for(i=0;iheight;i+)for(intj=0;jwidth;j+)unsignedcharnData=(unsignedchar)cvmGet(pGrayMat,i,j);dHistogramnData+;/建立直方图 doublePa=0.0;/区域1平均灰度值 doublePb=0.0;/区域2平均灰度值 doubleP0=0.0;/全局平均灰度值 doubleWa=0.0;/第一部分熵 doubleWb=0.0;/第二部分的熵 doubledData1=0.0,dData2=0.0;/中间值 doubledData3=0.0,dData4=0.0;/中间值 for(i=0;i256;i+)/计算全局平均灰度 dHistogrami/=N;P0+=i*dHistogrami;for(i=0;i256;i+)Wa=Wb=dData1=dData2=dData3=dData4=Pa=Pb=0.0;for(intj=0;j256;j+)if(j=i)dData1+=dHistogramj;dData2+=j*dHistogramj;elsedData3+=dHistogramj;dData4+=j*dHistogramj;Pa=dData2/dData1;Pb=dData4/dData3;for(j=0;j256;j+)if(j=i)if(Pa!=0)&(dHistogramj!=0)doubled1=log(dHistogramj/Pa);Wa+=j*dHistogramj*d1/log(2);elseif(Pb!=0)&(dHistogramj!=0)doubled2=log(dHistogramj/Pb);Wb+=j*dHistogramj*d2/log(2);dEntropyi=Wa+Wb;/遍历熵值，求取最小交叉熵所对应的灰度值 doubletemp=dEntropy0;for(i=1;i256;i+)if(dEntropyitemp)temp=dEntropyi;nThreshold=i;/* 函数名称：* MiniCross()* 输入参数:* pGrayMat: 二值图像数据*width: 图形尺寸宽度*height: 图形尺寸高度*nTlreshold: 经过算法处理得到的二值化分割阈值* 返回值:* 无* 函数说明:实现灰度图的二值化分割最小交叉熵算法*/void CBinarizationDlg:MiniCross(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)double dHistogram256; /灰度直方图double dEntropy256; /每个像素的交叉熵int N = height*width; /总像素数for(int i=0; i256; i+)dHistogrami = 0.0;dEntropyi = 0.0;for(i=0; iheight; i+)for(int j=0; jwidth; j+)unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);dHistogramnData+; /建立直方图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; i256; i+) /计算全局平均灰度dHistogrami /= N;P0 += i*dHistogrami;for(i=0; i256; i+)Wa=Wb=dData1=dData2=dData3=dData4=Pa=Pb=0.0;for(int j=0; j256; j+)if(j=i)dData1 += dHistogramj;dData2 += j*dHistogramj;elsedData3 += dHistogramj;dData4 += j*dHistogramj;Pa = dData2/dData1;Pb = dData4/dData3;for(j=0; j256; j+)if(j=i)if(Pa!=0)&(dHistogramj!=0)double d1 = log(dHistogramj/Pa);Wa += j*dHistogramj*d1/log(2);elseif(Pb!=0)&(dHistogramj!=0)double d2 = log(dHistogramj/Pb);Wb += j*dHistogramj*d2/log(2);dEntropyi = Wa+Wb;/遍历熵值，求取最小交叉熵所对应的灰度值double temp=dEntropy0;for(i=1; i256; i+)if(dEntropyitemp)temp = dEntropyi;nThreshold = i; 阈值分割结果如下图，求解所得的阈值为106.3、二维OTSU法 这种方法是对类间最大方差法的扩展，将其从求两个一维分布最大类间方差扩充为求解类间离散度矩阵的迹的最大值，考虑像素点灰度级的基础上增加了对像素点邻域平均像素值的考虑。 以下按照本人的理解对该方法的思路以及推倒过程进行分析： 1）首先需要建立二维的灰度统计直方图P(f, g)； 图像的灰度级为L级，那么其每个像素点的8邻域灰度平均值的灰度级也为L级，据此来构建直方图P。二维统计直方图的横轴为每个像素点的灰度值f(i, j)，纵坐标为同一个点对应的邻域平均值g(i, j) 其中(0iheight, 0jwidth, 0f(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也可以，因为从求解结果可以看出，这这个数值往往很接近。 具体的实现代码如下：cpp view plaincopyprint?/*函数名称：*TwoDimentionOtsu()*输入参数:*pGrayMat:二值图像数据*width:图形尺寸宽度*height:图形尺寸高度*nTlreshold:经过算法处理得到的二值化分割阈值*返回值:*无*函数说明:实现灰度图的二值化分割最大类间方差法（二维Otsu算法）*备注：在构建二维直方图的时候，采用灰度点的3*3邻域均值*/voidCBinarizationDlg:TwoDimentionOtsu(CvMat*pGrayMat,intwidth,intheight,BYTE&nThreshold)doubledHistogram256256;/建立二维灰度直方图 doubledTrMatrix=0.0;/离散矩阵的迹 intN=height*width;/总像素数 for(inti=0;i256;i+)for(intj=0;j256;j+)dHistogramij=0.0;/初始化变量 for(i=0;iheight;i+)for(intj=0;jwidth;j+)unsignedcharnData1=(unsignedchar)cvmGet(pGrayMat,i,j);/当前的灰度值 unsignedcharnData2=0;intnData3=0;/注意9个值相加可能超过一个字节 for(intm=i-1;m=i+1;m+)for(intn=j-1;n=0)&(m=0)&(nwidth)nData3+=(unsignedchar)cvmGet(pGrayMat,m,n);/当前的灰度值 nData2=(unsignedchar)(nData3/9);/对于越界的索引值进行补零,邻域均值 dHistogramnData1nData2+;for(i=0;i256;i+)for(intj=0;j256;j+)dHistogramij/=N;/得到归一化的概率分布 doublePai=0.0;/目标区均值矢量i分量 doublePaj=0.0;/目标区均值矢量j分量 doublePbi=0.0;/背景区均值矢量i分量 doublePbj=0.0;/背景区均值矢量j分量 doublePti=0.0;/全局均值矢量i分量 doublePtj=0.0;/全局均值矢量j分量 doubleW0=0.0;/目标区的联合概率密度 doubleW1=0.0;/背景区的联合概率密度 doubledData1=0.0;doubledData2=0.0;doubledData3=0.0;doubledData4=0.0;/中间变量 intnThreshold_s=0;intnThreshold_t=0;doubletemp=0.0;/寻求最大值 for(i=0;i256;i+)for(intj=0;j256;j+)Pti+=i*dHistogramij;Ptj+=j*dHistogramij;for(i=0;i256;i+)for(intj=0;j256;j+)W0+=dHistogramij;dData1+=i*dHistogramij;dData2+=j*dHistogramij;W1=1-W0;dData3=Pti-dData1;dData4=Ptj-dData2;/*W1=dData3=dData4=0.0;/对内循环的数据进行初始化for(ints=i+1;s256;s+)for(intt=j+1;ttemp)temp=dTrMatrix;nThreshold_s=i;nThreshold_t=j;nThreshold=nThreshold_t;/返回结果中的灰度值 /nThreshold=100; /* 函数名称：* TwoDimentionOtsu()* 输入参数:* pGrayMat: 二值图像数据*width: 图形尺寸宽度*height: 图形尺寸高度*nTlreshold: 经过算法处理得到的二值化分割阈值* 返回值:* 无* 函数说明:实现灰度图的二值化分割最大类间方差法（二维Otsu算法）* 备注：在构建二维直方图的时候，采用灰度点的3*3邻域均值*/void CBinarizationDlg:TwoDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)double dHistogram256256; /建立二维灰度直方图double dTrMatrix = 0.0; /离散矩阵的迹int N = height*width; /总像素数for(int i=0; i256; i+)for(int j=0; j256; j+)dHistogramij = 0.0; /初始化变量for(i=0; iheight; i+)for(int j=0; jwidth; 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=0)&(m=0)&(nwidth)nData3 += (unsigned char)cvmGet(pGrayMat, m, n); /当前的灰度值nData2 = (unsigned char)(nData3/9); /对于越界的索引值进行补零,邻域均值dHistogramnData1nData2+;for(i=0; i256; i+)for(int j=0; j256; j+)dHistogramij /= 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; i256; i+)for(int j=0; j256; j+)Pti += i*dHistogramij;Ptj += j*dHistogramij;for(i=0; i256; i+)for(int j=0; j256; j+)W0 += dHistogramij;dData1 += i*dHistogramij;dData2 += j*dHistogramij;W1 = 1-W0;dData3 = Pti-dData1;dData4 = Ptj-dData2;/*W1=dData3=dData4=0.0; /对内循环的数据进行初始化for(int s=i+1; s256; s+)for(int t=j+1; t temp)temp = dTrMatrix;nThreshold_s = i;nThreshold_t = j;nThreshold = nThreshold_t; /返回结果中的灰度值/nThreshold = 100; 阈值分割结果如下图，求解所得的阈值为114.（s=114，t=117）参考文献 1 Nobuyuki Otsu. A Threshold SelectionMethod from Gray-Level His