CSDN亲密携手阿里云重磅推出云邮箱服务HTML5群组诚募管理员，“活跃之星”活动火热进行中Canny边缘检测算法原理及其VC实现详解(一)2011-10-20 21:39560人阅读评论(0)收藏举报图象的边缘是指图象局部区域亮度变化显著的部分，该区域的灰度剖面一般可以看作是一个阶跃，既从一个灰度值在很小的缓冲区域内急剧变化到另一个灰度相差较大的灰度值。

图象的边缘部分集中了图象的大部分信息，图象边缘的确定与提取对于整个图象场景的识别与理解是非常重要的，同时也是图象分割所依赖的重要特征，边缘检测主要是图象的灰度变化的度量、检测和定位，自从1959提出边缘检测以来，经过五十多年的发展，已有许多中不同的边缘检测方法。

根据作者的理解和实践，本文对边缘检测的原理进行了描述，在此基础上着重对Canny检测算法的实现进行详述。

本文所述内容均由编程验证而来，在实现过程中，有任何错误或者不足之处大家共同讨论（本文不讲述枯燥的理论证明和数学推导，仅仅从算法的实现以及改进上进行原理性和工程化的描述）。

1、边缘检测原理及步骤在之前的博文中，作者从一维函数的跃变检测开始，循序渐进的对二维图像边缘检测的基本原理进行了通俗化的描述。

结论是：实现图像的边缘检测，就是要用离散化梯度逼近函数根据二维灰度矩阵梯度向量来寻找图像灰度矩阵的灰度跃变位置，然后在图像中将这些位置的点连起来就构成了所谓的图像边缘（图像边缘在这里是一个统称，包括了二维图像上的边缘、角点、纹理等基元图）。

在实际情况中理想的灰度阶跃及其线条边缘图像是很少见到的，同时大多数的传感器件具有低频滤波特性，这样会使得阶跃边缘变为斜坡性边缘，看起来其中的强度变化不是瞬间的，而是跨越了一定的距离。

这就使得在边缘检测中首先要进行的工作是滤波。

1）滤波：边缘检测的算法主要是基于图像强度的一阶和二阶导数，但导数通常对噪声很敏感，因此必须采用滤波器来改善与噪声有关的边缘检测器的性能。

常见的滤波方法主要有高斯滤波，即采用离散化的高斯函数产生一组归一化的高斯核（具体见“高斯滤波原理及其编程离散化实现方法”一文），然后基于高斯核函数对图像灰度矩阵的每一点进行加权求和（具体程序实现见下文）。

2）增强：增强边缘的基础是确定图像各点邻域强度的变化值。

增强算法可以将图像灰度点邻域强度值有显著变化的点凸显出来。

在具体编程实现时，可通过计算梯度幅值来确定。

3）检测：经过增强的图像，往往邻域中有很多点的梯度值比较大，而在特定的应用中，这些点并不是我们要找的边缘点，所以应该采用某种方法来对这些点进行取舍。

实际工程中，常用的方法是通过阈值化方法来检测。

2、Canny边缘检测算法原理JohnCanny于1986年提出Canny算子，它与Marr（LoG）边缘检测方法类似，也属于是先平滑后求导数的方法。

本节对根据上述的边缘检测过程对Canny检测算法的原理进行介绍。

2.1 对原始图像进行灰度化Canny算法通常处理的图像为灰度图，因此如果摄像机获取的是彩色图像，那首先就得进行灰度化。

对一幅彩色图进行灰度化，就是根据图像各个通道的采样值进行加权平均。

以RGB格式的彩图为例，通常灰度化采用的方法主要有：方法1：Gray=(R+G+B)/3;方法2：Gray=0.299R+0.587G+0.114B;（这种参数考虑到了人眼的生理特点）注意1：至于其他格式的彩色图像，可以根据相应的转换关系转为RGB 然后再进行灰度化；注意2：在编程时要注意图像格式中RGB的顺序通常为BGR。

2.2 对图像进行高斯滤波图像高斯滤波的实现可以用两个一维高斯核分别两次加权实现，也可以通过一个二维高斯核一次卷积实现。

1）高斯核实现上式为离散化的一维高斯函数，确定参数就可以得到一维核向量。

上式为离散化的二维高斯函数，确定参数就可以得到二维核向量。

注意1：关于参数Sigma的取值详见上篇博文。

注意2：在求的高斯核后，要对整个核进行归一化处理。

2）图像高斯滤波对图像进行高斯滤波，听起来很玄乎，其实就是根据待滤波的像素点及其邻域点的灰度值按照一定的参数规则进行加权平均。

这样可以有效滤去理想图像中叠加的高频噪声。

通常滤波和边缘检测是矛盾的概念，抑制了噪声会使得图像边缘模糊，这回增加边缘定位的不确定性；而如果要提高边缘检测的灵敏度，同时对噪声也提高了灵敏度。

实际工程经验表明，高斯函数确定的核可以在抗噪声干扰和边缘检测精确定位之间提供较好的折衷方案。

这就是所谓的高斯图像滤波，具体实现代码见下文。

2.3 用一阶偏导的有限差分来计算梯度的幅值和方向关于图像灰度值得梯度可使用一阶有限差分来进行近似，这样就可以得图像在x和y方向上偏导数的两个矩阵。

常用的梯度算子有如下几种：1）Roberts算子上式为其x和y方向偏导数计算模板，可用数学公式表达其每个点的梯度幅值为：2）Sobel算子上式三个矩阵分别为该算子的x向卷积模板、y向卷积模板以及待处理点的邻域点标记矩阵，据此可用数学公式表达其每个点的梯度幅值为：3）Prewitt算子和Sobel算子原理一样，在此仅给出其卷积模板。

4）Canny算法所采用的方法在本文实现的Canny算法中所采用的卷积算子比较简单，表达如下：其x向、y向的一阶偏导数矩阵，梯度幅值以及梯度方向的数学表达式为：求出这几个矩阵后，就可以进行下一步的检测过程。

2.4 对梯度幅值进行非极大值抑制图像梯度幅值矩阵中的元素值越大，说明图像中该点的梯度值越大，但这不不能说明该点就是边缘（这仅仅是属于图像增强的过程）。

在Canny 算法中，非极大值抑制是进行边缘检测的重要步骤，通俗意义上是指寻找像素点局部最大值，将非极大值点所对应的灰度值置为0，这样可以剔除掉一大部分非边缘的点（这是本人的理解）。

图1 非极大值抑制原理根据图1 可知，要进行非极大值抑制，就首先要确定像素点C的灰度值在其8值邻域内是否为最大。

图1中蓝色的线条方向为C点的梯度方向，这样就可以确定其局部的最大值肯定分布在这条线上，也即出了C点外，梯度方向的交点dTmp1和dTmp2这两个点的值也可能会是局部最大值。

因此，判断C点灰度与这两个点灰度大小即可判断C点是否为其邻域内的局部最大灰度点。

如果经过判断，C点灰度值小于这两个点中的任一个，那就说明C点不是局部极大值，那么则可以排除C 点为边缘。

这就是非极大值抑制的工作原理。

作者认为，在理解的过程中需要注意以下两点：1）中非最大抑制是回答这样一个问题：“当前的梯度值在梯度方向上是一个局部最大值吗？” 所以,要把当前位置的梯度值与梯度方向上两侧的梯度值进行比较；2）梯度方向垂直于边缘方向。

但实际上，我们只能得到C点邻域的8个点的值，而dTmp1和dTmp2并不在其中，要得到这两个值就需要对该两个点两端的已知灰度进行线性插值，也即根据图1中的g1和g2对dTmp1进行插值，根据g3和g4对dTmp2进行插值，这要用到其梯度方向，这是上文Canny算法中要求解梯度方向矩阵Thita的原因。

完成非极大值抑制后，会得到一个二值图像，非边缘的点灰度值均为0，可能为边缘的局部灰度极大值点可设置其灰度为128。

根据下文的具体测试图像可以看出，这样一个检测结果还是包含了很多由噪声及其他原因造成的假边缘。

因此还需要进一步的处理。

2.5 用双阈值算法检测和连接边缘Canny算法中减少假边缘数量的方法是采用双阈值法。

选择两个阈值（关于阈值的选取方法在扩展中进行讨论），根据高阈值得到一个边缘图像，这样一个图像含有很少的假边缘，但是由于阈值较高，产生的图像边缘可能不闭合，未解决这样一个问题采用了另外一个低阈值。

在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，该算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合。

以上即为整个Canny边缘检测算法的原理分析，接下来我们进行VC 下的算法实现和效果分析。

3、Canny算法的实现流程由于本文主要目的在于学习和实现算法，而对于图像读取、视频获取等内容不进行阐述。

因此选用OpenCV算法库作为其他功能的实现途径（关于OpenCV的使用，作者将另文表述）。

首先展现本文将要处理的彩色图片。

图2 待处理的图像3.1 图像读取和灰度化编程时采用上文所描述的第二种方法来实现图像的灰度化。

其中ptr数组中保存的灰度化后的图像数据。

具体的灰度化后的效果如图3所示。

[cpp]view plaincopyprint?1.IplImage* ColorImage = cvLoadImage( "12.jpg", -1 ); //读入图像，获取彩图指针2.IplImage* OpenCvGrayImage; //定义变换后的灰度图指针3.unsigned char* ptr; //指向图像的数据首地址4.if (ColorImage == NULL)5.return;6.int i = ColorImage->width * ColorImage->height;7.BYTE data1; //中间过程变量8.BYTE data2;9.BYTE data3;10.ptr = new unsigned char[i];11.for(intj=0; j<ColorImage->height; j++) //对RGB加权平均，权值参考OpenCV12.{13.for(intx=0; x<ColorImage->width; x++)14.{15.data1 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3]; //B分量16.data2 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3 + 1]; //G分量17.data3 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3 + 2]; //R分量18.ptr[j*ColorImage->width+x]=(BYTE)(0.072169*data1 + 0.715160*data2 +0.212671*data3);19.}20.}21.OpenCvGrayImage=cvCreateImageHeader(cvGetSize(ColorImage), ColorImage->depth,1);22.cvSetData(GrayImage,ptr, GrayImage->widthStep); //根据数据生成灰度图23.cvNamedWindow("GrayImage",CV_WINDOW_AUTOSIZE);24.cvShowImage("GrayImage",OpenCvGrayImage); //显示灰度图25.cvWaitKey(0);26.cvDestroyWindow("GrayImage");1.double nSigma = 0.4; //定义高斯函数的标准差2.int nWidowSize = 1+2*ceil(3*nSigma); //定义滤波窗口的大小3.int nCenter = (nWidowSize)/2; //定义滤波窗口中心的索引1.int nWidth = OpenCvGrayImage->width; //获取图像的像素宽度2.int nHeight = OpenCvGrayImage->height; //获取图像的像素高度3.unsigned char* nImageData = new unsigned char[nWidth*nHeight]; //暂时保存图像中的数据4.unsigned char*pCanny = new unsigned char[nWidth*nHeight]; //为平滑后的图像数据分配内存5.double* nData = new double[nWidth*nHeight]; //两次平滑的中间数据6.for(int j=0; j<nHeight; j++) //获取数据7.{8.for(i=0; i<nWidth; i++)9.nImageData[j*nWidth+i] = (unsignedchar)OpenCvGrayImage->imageData[j*nWidth+i];10.}1.//////////////////////生成一维高斯滤波系数/////////////////////////////2.double* pdKernal_1 = new double[nWidowSize]; //定义一维高斯核数组3.double dSum_1 = 0.0; //求和，用于进行归一化4.////////////////////////一维高斯函数公式//////////////////////////////5.//// x*x /////////////////6.//// -1*---------------- /////////////////7.//// 1 2*Sigma*Sigma /////////////////8.//// ------------ e /////////////////9.//// /////////////////10.//// \/2*pi*Sigma /////////////////11.//////////////////////////////////////////////////////////////////////12.for(int i=0; i<nWidowSize; i++)13.{14.double nDis = (double)(i-nCenter);15.pdKernal_1[i] =exp(-(0.5)*nDis*nDis/(nSigma*nSigma))/(sqrt(2*3.14159)*nSigma);16.dSum_1 += pdKernal_1[i];17.}18.for(i=0; i<nWidowSize; i++)19.{20.pdKernal_1[i] /= dSum_1; //进行归一化21.}1.for(i=0; i<nHeight; i++) //进行x向的高斯滤波(加权平均)2.{3.for(j=0; j<nWidth; j++)4.{5.double dSum = 0;6.double dFilter=0; //滤波中间值7.for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)8.{9.if((j+nLimit)>=0 && (j+nLimit) < nWidth ) //图像不能超出边界10.{11.dFilter += (double)nImageData[i*nWidth+j+nLimit] * pdKernal_1[nCenter+nLimit];12.dSum += pdKernal_1[nCenter+nLimit];13.}14.}15.nData[i*nWidth+j] = dFilter/dSum;16.}17.}18.19.for(i=0; i<nWidth; i++) //进行y向的高斯滤波(加权平均)20.{21.for(j=0; j<nHeight; j++)22.{23.double dSum = 0.0;24.double dFilter=0;25.for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)26.{27.if((j+nLimit)>=0 && (j+nLimit) < nHeight) //图像不能超出边界28.{29.dFilter += (double)nData[(j+nLimit)*nWidth+i] * pdKernal_1[nCenter+nLimit];30.dSum += pdKernal_1[nCenter+nLimit];31.}32.}33.pCanny[j*nWidth+i] = (unsigned char)(int)dFilter/dSum;34.}35.}3.2.2 根据二维高斯核进行滤波1）生成二维高斯滤波系数[cpp]view plaincopyprint?1.//////////////////////生成一维高斯滤波系数//////////////////////////////////2.double* pdKernal_2 = new double[nWidowSize*nWidowSize]; //定义一维高斯核数组3.double dSum_2 = 0.0; //求和，进行归一化4.///////////////////////二维高斯函数公式////////////////////////////////////5.//// x*x+y*y ///////////////6.//// -1*-------------- ///////////////7.//// 1 2*Sigma*Sigma ///////////////8.//// ---------------- e ///////////////9.//// 2*pi*Sigma*Sigma ///////////////10.///////////////////////////////////////////////////////////////////////////11.for(i=0; i<nWidowSize; i++)12.{13.for(int j=0; j<nWidowSize; j++)14.{15.int nDis_x = i-nCenter;16.int nDis_y = j-nCenter;17.pdKernal_2[i+j*nWidowSize]=exp(-(1/2)*(nDis_x*nDis_x+nDis_y*nDis_y)18./(nSigma*nSigma))/(2*3.1415926*nSigma*nSigma);19.dSum_2 += pdKernal_2[i+j*nWidowSize];20.}21.}22.for(i=0; i<nWidowSize; i++)23.{24.for(int j=0; j<nWidowSize; j++) //进行归一化25.{26.pdKernal_2[i+j*nWidowSize] /= dSum_2;27.}28.}1.int x;2.int y;3.for(i=0; i<nHeight; i++)4.{5.for(j=0; j<nWidth; j++)6.{7.double dFilter=0.0;8.double dSum = 0.0;9.for(x=(-nCenter); x<=nCenter; x++) //行10.{11.for(y=(-nCenter); y<=nCenter; y++) //列12.{13.if( (j+x)>=0 && (j+x)<nWidth && (i+y)>=0 && (i+y)<nHeight) //判断边缘14.{15.dFilter += (double)nImageData [(i+y)*nWidth + (j+x)]16.* pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];17.dSum += pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];18.}19.}20.}21.pCanny[i*nWidth+j] = (unsigned char)dFilter/dSum;22.}23.}3.3 图像增强——计算图像梯度及其方向根据上文分析可知，实现代码如下[cpp]view plaincopyprint?1.//////////////////同样可以用不同的检测器/////////////////////////2.///// P[i,j]=(S[i,j+1]-S[i,j]+S[i+1,j+1]-S[i+1,j])/2 /////3.///// Q[i,j]=(S[i,j]-S[i+1,j]+S[i,j+1]-S[i+1,j+1])/2 /////4./////////////////////////////////////////////////////////////////5.double* P = new double[nWidth*nHeight]; //x向偏导数6.double* Q = new double[nWidth*nHeight]; //y向偏导数7.int* M = new int[nWidth*nHeight]; //梯度幅值8.double* Theta = new double[nWidth*nHeight]; //梯度方向9.//计算x,y方向的偏导数10.for(i=0; i<(nHeight-1); i++)11.{12.for(j=0; j<(nWidth-1); j++)13.{14.P[i*nWidth+j] = (double)(pCanny[i*nWidth + min(j+1, nWidth-1)] -pCanny[i*nWidth+j] + pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+j])/2;15.Q[i*nWidth+j] = (double)(pCanny[i*nWidth+j] - pCanny[min(i+1,nHeight-1)*nWidth+j] + pCanny[i*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)])/2;16.}17.}18.//计算梯度幅值和梯度的方向19.for(i=0; i<nHeight; i++)20.{21.for(j=0; j<nWidth; j++)22.{23.M[i*nWidth+j] = (int)(sqrt(P[i*nWidth+j]*P[i*nWidth+j] +Q[i*nWidth+j]*Q[i*nWidth+j])+0.5);24.Theta[i*nWidth+j] = atan2(Q[i*nWidth+j], P[i*nWidth+j]) * 57.3;25.if(Theta[i*nWidth+j] < 0)26.Theta[i*nWidth+j] += 360; //将这个角度转换到0~360范围27.}28.}1.unsigned char* N = new unsigned char[nWidth*nHeight]; //非极大值抑制结果2.int g1=0, g2=0, g3=0, g4=0; //用于进行插值，得到亚像素点坐标值3.double dTmp1=0.0, dTmp2=0.0; //保存两个亚像素点插值得到的灰度数据4.double dWeight=0.0; //插值的权重1.for(i=0; i<nWidth; i++)2.{3.N[i] = 0;4.N[(nHeight-1)*nWidth+i] = 0;5.}6.for(j=0; j<nHeight; j++)7.{8.N[j*nWidth] = 0;9.N[j*nWidth+(nWidth-1)] = 0;10.}1.for(i=1; i<(nWidth-1); i++)2.{3.for(j=1; j<(nHeight-1); j++)4.{5.int nPointIdx = i+j*nWidth; //当前点在图像数组中的索引值6.if(M[nPointIdx] == 0)7.N[nPointIdx] = 0; //如果当前梯度幅值为0，则不是局部最大对该点赋为08.else9.{10.////////首先判断属于那种情况，然后根据情况插值///////11.////////////////////第一种情况///////////////////////12.///////// g1 g2 /////////////13.///////// C /////////////14.///////// g3 g4 /////////////15./////////////////////////////////////////////////////16.if( ((Theta[nPointIdx]>=90)&&(Theta[nPointIdx]<135)) ||17.((Theta[nPointIdx]>=270)&&(Theta[nPointIdx]<315)))18.{19.//////根据斜率和四个中间值进行插值求解20.g1 = M[nPointIdx-nWidth-1];21.g2 = M[nPointIdx-nWidth];22.g3 = M[nPointIdx+nWidth];23.g4 = M[nPointIdx+nWidth+1];24.dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]); //反正切25.dTmp1 = g1*dWeight+g2*(1-dWeight);26.dTmp2 = g4*dWeight+g3*(1-dWeight);27.}28.////////////////////第二种情况///////////////////////29.///////// g1 /////////////30.///////// g2 C g3 /////////////31.///////// g4 /////////////32./////////////////////////////////////////////////////33.else if( ((Theta[nPointIdx]>=135)&&(Theta[nPointIdx]<180)) ||34.((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360)))35.{36.g1 = M[nPointIdx-nWidth-1];37.g2 = M[nPointIdx-1];38.g3 = M[nPointIdx+1];39.g4 = M[nPointIdx+nWidth+1];40.dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]); //正切41.dTmp1 = g2*dWeight+g1*(1-dWeight);42.dTmp2 = g4*dWeight+g3*(1-dWeight);43.}44.////////////////////第三种情况///////////////////////45.///////// g1 g2 /////////////46.///////// C /////////////47.///////// g4 g3 /////////////48./////////////////////////////////////////////////////49.else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) ||50.((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270)))51.{52.g1 = M[nPointIdx-nWidth];53.g2 = M[nPointIdx-nWidth+1];54.g3 = M[nPointIdx+nWidth];55.g4 = M[nPointIdx+nWidth-1];56.dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]); //反正切57.dTmp1 = g2*dWeight+g1*(1-dWeight);58.dTmp2 = g3*dWeight+g4*(1-dWeight);59.}60.////////////////////第四种情况///////////////////////61.///////// g1 /////////////62.///////// g4 C g2 /////////////63.///////// g3 /////////////64./////////////////////////////////////////////////////65.else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) ||66.((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225)))67.{68.g1 = M[nPointIdx-nWidth+1];69.g2 = M[nPointIdx+1];70.g3 = M[nPointIdx+nWidth-1];71.g4 = M[nPointIdx-1];72.dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]); //正切73.dTmp1 = g1*dWeight+g2*(1-dWeight);74.dTmp2 = g3*dWeight+g4*(1-dWeight);75.}76.}77.//////////进行局部最大值判断，并写入检测结果////////////////78.if((M[nPointIdx]>=dTmp1) && (M[nPointIdx]>=dTmp2))79.N[nPointIdx] = 128;80.else81.N[nPointIdx] = 0;82.}83.}1.int nHist[1024];2.int nEdgeNum; //可能边界数3.int nMaxMag = 0; //最大梯度数4.int nHighCount;1.for(i=0;i<1024;i++)2.nHist[i] = 0;3.for(i=0; i<nHeight; i++)4.{5.for(j=0; j<nWidth; j++)6.{7.if(N[i*nWidth+j]==128)8.nHist[M[i*nWidth+j]]++;9.}10.}[cpp]view plaincopyprint?1.nEdgeNum = nHist[0];2.nMaxMag = 0; //获取最大的梯度值3.for(i=1; i<1024; i++) //统计经过“非最大值抑制”后有多少像素4.{5.if(nHist[i] != 0) //梯度为0的点是不可能为边界点的6.{7.nMaxMag = i;8.}9.nEdgeNum += nHist[i]; //经过non-maximum suppression后有多少像素10.}1.double dRatHigh = 0.79;2.double dThrHigh;3.double dThrLow;4.double dRatLow = 0.5;5.nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);6.j=1;7.nEdgeNum = nHist[1];8.while((j<(nMaxMag-1)) && (nEdgeNum < nHighCount))9.{10.j++;11.nEdgeNum += nHist[j];12.}13.dThrHigh = j; //高阈值14.dThrLow = (int)((dThrHigh) * dRatLow + 0.5); //低阈值1.SIZE sz;2.sz.cx = nWidth;3.sz.cy = nHeight;4.for(i=0; i<nHeight; i++)5.{6.for(j=0; j<nWidth; j++)7.{8.if((N[i*nWidth+j]==128) && (M[i*nWidth+j] >= dThrHigh))9.{10.N[i*nWidth+j] = 255;11.TraceEdge(i, j, dThrLow, N, M, sz);12.}13.}14.}成一个闭合的轮廓。