canny算法的MATLAB和C++实现（转载）

图象的边缘是指图象局部区域亮度变化显著的部分，该区域的灰度剖面一般可以看作是一个阶跃，既从一个灰度值在很小的缓冲区域内急剧变化到另一个灰度相差较大的灰度值。图象的边缘部分集中了图象的大部分信息，图象边缘的确定与提取对于整个图象场景的识别与理解是非常重要的，同时也是图象分割所依赖的重要特征，边缘检测主要是图象的灰度变化的度量、检测和定位，自从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所示。

  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加权平均，权值参考OpenCV

  12. {

  13. for(intx=0; x<ColorImage->width; x++)

  14. {

  15. data1 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + x*3];     //B分量

  16. data2 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + x*3 + 1]; //G分量

  17. data3 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + x*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");

图3 灰度化后的图像

3.2 图像的高斯滤波

根据上面所讲的边缘检测过程，下一个步骤就是对图像进行高斯滤波。可根据之前博文描述的方法获取一维或者二维的高斯滤波核。因此进行图像高斯滤波可有两种实现方式，以下具体进行介绍。

首先定义该部分的通用变量：

  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] = (unsigned char)OpenCvGrayImage->imageData[j*nWidth+i];

  10. }

3.2.1 根据一维高斯核进行两次滤波

1）生成一维高斯滤波系数

  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. }


      

2）分别进行x向和y向的一维加权滤波，滤波后的数据保存在矩阵pCanny中


  1. for(i=0; i<nHeight; i++)                               //进行x向的高斯滤波(加权平均)

    2. 
  
  2. {

    3. 
  
  3. for(j=0; j<nWidth; j++)

    4. 
  
  4. {

    5. 
  
  5. double dSum = 0;

    6. 
  
  6. double dFilter=0;                                       //滤波中间值

    7. 
  
  7. for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)

    8. 
  
  8. {

    9. 
  
  9. if((j+nLimit)>=0 && (j+nLimit) < nWidth )       //图像不能超出边界

    10. 
  
  10. {

    11. 
  
  11. dFilter += (double)nImageData[i*nWidth+j+nLimit] * pdKernal_1[nCenter+nLimit];

    12. 
  
  12. dSum += pdKernal_1[nCenter+nLimit];

    13. 
  
  13. }

    14. 
  
  14. }

    15. 
  
  15. nData[i*nWidth+j] = dFilter/dSum;

    16. 
  
  16. }

    17. 
  
  17. }
    18. 

  18. for(i=0; i<nWidth; i++)                                //进行y向的高斯滤波(加权平均)

    19. 
  
  19. {

    20. 
  
  20. for(j=0; j<nHeight; j++)

    21. 
  
  21. {

    22. 
  
  22. double dSum = 0.0;

    23. 
  
  23. double dFilter=0;

    24. 
  
  24. for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)

    25. 
  
  25. {

    26. 
  
  26. if((j+nLimit)>=0 && (j+nLimit) < nHeight)       //图像不能超出边界

    27. 
  
  27. {

    28. 
  
  28. dFilter += (double)nData[(j+nLimit)*nWidth+i] * pdKernal_1[nCenter+nLimit];

    29. 
  
  29. dSum += pdKernal_1[nCenter+nLimit];

    30. 
  
  30. }

    31. 
  
  31. }

    32. 
  
  32. pCanny[j*nWidth+i] = (unsigned char)(int)dFilter/dSum;

    33. 
  
  33. }

    34. 
  
  34. }




3.2.2 根据二维高斯核进行滤波



1）生成二维高斯滤波系数


  1. //////////////////////生成一维高斯滤波系数//////////////////////////////////

    2. 
  
  2. double* pdKernal_2 = new double[nWidowSize*nWidowSize]; //定义一维高斯核数组

    3. 
  
  3. double  dSum_2 = 0.0;                                   //求和，进行归一化

    4. 
  
  4. ///////////////////////二维高斯函数公式////////////////////////////////////

    5. 
  
  5. ////                         x*x+y*y                        ///////////////

    6. 
  
  6. ////                   -1*--------------                ///////////////

    7. 
  
  7. ////         1             2*Sigma*Sigma                ///////////////

    8. 
  
  8. ////   ---------------- e                                   ///////////////

    9. 
  
  9. ////   2*pi*Sigma*Sigma                                     ///////////////

    10. 
  
  10. ///////////////////////////////////////////////////////////////////////////

    11. 
  
  11. for(i=0; i<nWidowSize; i++)

    12. 
  
  12. {

    13. 
  
  13. for(int j=0; j<nWidowSize; j++)

    14. 
  
  14. {

    15. 
  
  15. int nDis_x = i-nCenter;

    16. 
  
  16. int nDis_y = j-nCenter;

    17. 
  
  17. pdKernal_2[i+j*nWidowSize]=exp(-(1/2)*(nDis_x*nDis_x+nDis_y*nDis_y)

    18. 
  
  18. /(nSigma*nSigma))/(2*3.1415926*nSigma*nSigma);

    19. 
  
  19. dSum_2 += pdKernal_2[i+j*nWidowSize];

    20. 
  
  20. }

    21. 
  
  21. }

    22. 
  
  22. for(i=0; i<nWidowSize; i++)

    23. 
  
  23. {

    24. 
  
  24. for(int j=0; j<nWidowSize; j++)                 //进行归一化

    25. 
  
  25. {

    26. 
  
  26. pdKernal_2[i+j*nWidowSize] /= dSum_2;

    27. 
  
  27. }

    28. 
  
  28. }


    

 2）采用高斯核进行高斯滤波，滤波后的数据保存在矩阵pCanny中


  1. int x;

    2. 
  
  2. int y;

    3. 
  
  3. for(i=0; i<nHeight; i++)

    4. 
  
  4. {

    5. 
  
  5. for(j=0; j<nWidth; j++)

    6. 
  
  6. {

    7. 
  
  7. double dFilter=0.0;

    8. 
  
  8. double dSum = 0.0;

    9. 
  
  9. for(x=(-nCenter); x<=nCenter; x++)                     //行

    10. 
  
  10. {

    11. 
  
  11. for(y=(-nCenter); y<=nCenter; y++)             //列

    12. 
  
  12. {

    13. 
  
  13. if( (j+x)>=0 && (j+x)<nWidth && (i+y)>=0 && (i+y)<nHeight) //判断边缘

    14. 
  
  14. {

    15. 
  
  15. dFilter += (double)nImageData [(i+y)*nWidth + (j+x)]

    16. 
  
  16. * pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];

    17. 
  
  17. dSum += pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];

    18. 
  
  18. }

    19. 
  
  19. }

    20. 
  
  20. }

    21. 
  
  21. pCanny[i*nWidth+j] = (unsigned char)dFilter/dSum;

    22. 
  
  22. }

    23. 
  
  23. }


3.3

图像增强——计算图像梯度及其方向


    

 根据上文分析可知，实现代码如下


  1. //////////////////同样可以用不同的检测器/////////////////////////

    2. 
  
  2. /////    P[i,j]=(S[i,j+1]-S[i,j]+S[i+1,j+1]-S[i+1,j])/2     /////

    3. 
  
  3. /////    Q[i,j]=(S[i,j]-S[i+1,j]+S[i,j+1]-S[i+1,j+1])/2     /////

    4. 
  
  4. /////////////////////////////////////////////////////////////////

    5. 
  
  5. double* P = new double[nWidth*nHeight];                 //x向偏导数

    6. 
  
  6. double* Q = new double[nWidth*nHeight];                 //y向偏导数

    7. 
  
  7. int* M = new int[nWidth*nHeight];                       //梯度幅值

    8. 
  
  8. double* Theta = new double[nWidth*nHeight];             //梯度方向

    9. 
  
  9. //计算x,y方向的偏导数

    10. 
  
  10. for(i=0; i<(nHeight-1); i++)

    11. 
  
  11. {

    12. 
  
  12. for(j=0; j<(nWidth-1); j++)

    13. 
  
  13. {

    14. 
  
  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. 
  
  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. 
  
  16. }

    17. 
  
  17. }

    18. 
  
  18. //计算梯度幅值和梯度的方向

    19. 
  
  19. for(i=0; i<nHeight; i++)

    20. 
  
  20. {

    21. 
  
  21. for(j=0; j<nWidth; j++)

    22. 
  
  22. {

    23. 
  
  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. 
  
  24. Theta[i*nWidth+j] = atan2(Q[i*nWidth+j], P[i*nWidth+j]) * 57.3;

    25. 
  
  25. if(Theta[i*nWidth+j] < 0)

    26. 
  
  26. Theta[i*nWidth+j] += 360;              //将这个角度转换到0~360范围

    27. 
  
  27. }

    28. 
  
  28. }


3.4 非极大值抑制



    

 根据上文所述的工作原理，这部分首先需要求解每个像素点在其邻域内的梯度方向的两个灰度值，然后判断是否为潜在的边缘，如果不是则将该点灰度值设置为0.


    

 首先定义相关的参数如下：


  1. unsigned char* N = new unsigned char[nWidth*nHeight];  //非极大值抑制结果

    2. 
  
  2. int g1=0, g2=0, g3=0, g4=0;                            //用于进行插值，得到亚像素点坐标值

    3. 
  
  3. double dTmp1=0.0, dTmp2=0.0;                           //保存两个亚像素点插值得到的灰度数据

    4. 
  
  4. double dWeight=0.0;                                    //插值的权重

  

   其次，对边界进行初始化：


  1. for(i=0; i<nWidth; i++)

    2. 
  
  2. {

    3. 
  
  3. N[i] = 0;

    4. 
  
  4. N[(nHeight-1)*nWidth+i] = 0;

    5. 
  
  5. }

    6. 
  
  6. for(j=0; j<nHeight; j++)

    7. 
  
  7. {

    8. 
  
  8. N[j*nWidth] = 0;

    9. 
  
  9. N[j*nWidth+(nWidth-1)] = 0;

    10. 
  
  10. }

    

 进行局部最大值寻找，根据上文图1所述的方案进行插值，然后判优，实现代码如下：


  1. for(i=1; i<(nWidth-1); i++)

    2. 
  
  2. {

    3. 
  
  3. for(j=1; j<(nHeight-1); j++)

    4. 
  
  4. {

    5. 
  
  5. int nPointIdx = i+j*nWidth;       //当前点在图像数组中的索引值

    6. 
  
  6. if(M[nPointIdx] == 0)

    7. 
  
  7. N[nPointIdx] = 0;         //如果当前梯度幅值为0，则不是局部最大对该点赋为0

    8. 
  
  8. else

    9. 
  
  9. {

    10. 
  
  10. ////////首先判断属于那种情况，然后根据情况插值///////

    11. 
  
  11. ////////////////////第一种情况///////////////////////

    12. 
  
  12. /////////       g1  g2                  /////////////

    13. 
  
  13. /////////           C                   /////////////

    14. 
  
  14. /////////           g3  g4              /////////////

    15. 
  
  15. /////////////////////////////////////////////////////

    16. 
  
  16. if( ((Theta[nPointIdx]>=90)&&(Theta[nPointIdx]<135)) ||

    17. 
  
  17. ((Theta[nPointIdx]>=270)&&(Theta[nPointIdx]<315)))

    18. 
  
  18. {

    19. 
  
  19. //////根据斜率和四个中间值进行插值求解

    20. 
  
  20. g1 = M[nPointIdx-nWidth-1];

    21. 
  
  21. g2 = M[nPointIdx-nWidth];

    22. 
  
  22. g3 = M[nPointIdx+nWidth];

    23. 
  
  23. g4 = M[nPointIdx+nWidth+1];

    24. 
  
  24. dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);   //反正切

    25. 
  
  25. dTmp1 = g1*dWeight+g2*(1-dWeight);

    26. 
  
  26. dTmp2 = g4*dWeight+g3*(1-dWeight);

    27. 
  
  27. }

    28. 
  
  28. ////////////////////第二种情况///////////////////////

    29. 
  
  29. /////////       g1                      /////////////

    30. 
  
  30. /////////       g2  C   g3              /////////////

    31. 
  
  31. /////////               g4              /////////////

    32. 
  
  32. /////////////////////////////////////////////////////

    33. 
  
  33. else if( ((Theta[nPointIdx]>=135)&&(Theta[nPointIdx]<180)) ||

    34. 
  
  34. ((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360)))

    35. 
  
  35. {

    36. 
  
  36. g1 = M[nPointIdx-nWidth-1];

    37. 
  
  37. g2 = M[nPointIdx-1];

    38. 
  
  38. g3 = M[nPointIdx+1];

    39. 
  
  39. g4 = M[nPointIdx+nWidth+1];

    40. 
  
  40. dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);   //正切

    41. 
  
  41. dTmp1 = g2*dWeight+g1*(1-dWeight);

    42. 
  
  42. dTmp2 = g4*dWeight+g3*(1-dWeight);

    43. 
  
  43. }

    44. 
  
  44. ////////////////////第三种情况///////////////////////

    45. 
  
  45. /////////           g1  g2              /////////////

    46. 
  
  46. /////////           C                   /////////////

    47. 
  
  47. /////////       g4  g3                  /////////////

    48. 
  
  48. /////////////////////////////////////////////////////

    49. 
  
  49. else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) ||

    50. 
  
  50. ((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270)))

    51. 
  
  51. {

    52. 
  
  52. g1 = M[nPointIdx-nWidth];

    53. 
  
  53. g2 = M[nPointIdx-nWidth+1];

    54. 
  
  54. g3 = M[nPointIdx+nWidth];

    55. 
  
  55. g4 = M[nPointIdx+nWidth-1];

    56. 
  
  56. dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);   //反正切

    57. 
  
  57. dTmp1 = g2*dWeight+g1*(1-dWeight);

    58. 
  
  58. dTmp2 = g3*dWeight+g4*(1-dWeight);

    59. 
  
  59. }

    60. 
  
  60. ////////////////////第四种情况///////////////////////

    61. 
  
  61. /////////               g1              /////////////

    62. 
  
  62. /////////       g4  C   g2              /////////////

    63. 
  
  63. /////////       g3                      /////////////

    64. 
  
  64. /////////////////////////////////////////////////////

    65. 
  
  65. else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) ||

    66. 
  
  66. ((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225)))

    67. 
  
  67. {

    68. 
  
  68. g1 = M[nPointIdx-nWidth+1];

    69. 
  
  69. g2 = M[nPointIdx+1];

    70. 
  
  70. g3 = M[nPointIdx+nWidth-1];

    71. 
  
  71. g4 = M[nPointIdx-1];

    72. 
  
  72. dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);   //正切

    73. 
  
  73. dTmp1 = g1*dWeight+g2*(1-dWeight);

    74. 
  
  74. dTmp2 = g3*dWeight+g4*(1-dWeight);

    75. 
  
  75. }

    76. 
  
  76. }

    77. 
  
  77. //////////进行局部最大值判断，并写入检测结果////////////////

    78. 
  
  78. if((M[nPointIdx]>=dTmp1) && (M[nPointIdx]>=dTmp2))

    79. 
  
  79. N[nPointIdx] = 128;

    80. 
  
  80. else

    81. 
  
  81. N[nPointIdx] = 0;

    82. 
  
  82. }

    83. 
  
  83. }


3.5双阈值检测实现



1）定义相应参数如下


  1. int nHist[1024];

    2. 
  
  2. int nEdgeNum;             //可能边界数

    3. 
  
  3. int nMaxMag = 0;          //最大梯度数

    4. 
  
  4. int nHighCount;


  

   2）构造灰度图的统计直方图，根据上文梯度幅值的计算公式可知，最大的梯度幅值为：



  

   因此设置nHist为1024足够。以下实现统计直方图：


  1. for(i=0;i<1024;i++)

    2. 
  
  2. nHist[i] = 0;

    3. 
  
  3. for(i=0; i<nHeight; i++)

    4. 
  
  4. {

    5. 
  
  5. for(j=0; j<nWidth; j++)

    6. 
  
  6. {

    7. 
  
  7. if(N[i*nWidth+j]==128)

    8. 
  
  8. nHist[M[i*nWidth+j]]++;

    9. 
  
  9. }

    10. 
  
  10. }


    

 3）获取最大梯度幅值及潜在边缘点个数


  1. nEdgeNum = nHist[0];

    2. 
  
  2. nMaxMag = 0;                    //获取最大的梯度值

    3. 
  
  3. for(i=1; i<1024; i++)           //统计经过“非最大值抑制”后有多少像素

    4. 
  
  4. {

    5. 
  
  5. if(nHist[i] != 0)       //梯度为0的点是不可能为边界点的

    6. 
  
  6. {

    7. 
  
  7. nMaxMag = i;

    8. 
  
  8. }

    9. 
  
  9. nEdgeNum += nHist[i];   //经过non-maximum suppression后有多少像素

    10. 
  
  10. }


  

   4）计算两个阈值


  1. double  dRatHigh = 0.79;

    2. 
  
  2. double  dThrHigh;

    3. 
  
  3. double  dThrLow;

    4. 
  
  4. double  dRatLow = 0.5;

    5. 
  
  5. nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);

    6. 
  
  6. j=1;

    7. 
  
  7. nEdgeNum = nHist[1];

    8. 
  
  8. while((j<(nMaxMag-1)) && (nEdgeNum < nHighCount))

    9. 
  
  9. {

    10. 
  
  10. j++;

    11. 
  
  11. nEdgeNum += nHist[j];

    12. 
  
  12. }

    13. 
  
  13. dThrHigh = j;                                   //高阈值

    14. 
  
  14. dThrLow = (int)((dThrHigh) * dRatLow + 0.5);    //低阈值


  

 

 这段代码的意思是，按照灰度值从低到高的顺序，选取前79%个灰度值中的最大的灰度值为高阈值，低阈值大约为高阈值的一半。这是根据经验数据的来的，至于更好地参数选取方法，作者后面会另文研究。
  

   5）进行边缘检测


  1. SIZE sz;

    2. 
  
  2. sz.cx = nWidth;

    3. 
  
  3. sz.cy = nHeight;

    4. 
  
  4. for(i=0; i<nHeight; i++)

    5. 
  
  5. {

    6. 
  
  6. for(j=0; j<nWidth; j++)

    7. 
  
  7. {

    8. 
  
  8. if((N[i*nWidth+j]==128) && (M[i*nWidth+j] >= dThrHigh))

    9. 
  
  9. {

    10. 
  
  10. N[i*nWidth+j] = 255;

    11. 
  
  11. TraceEdge(i, j, dThrLow, N, M, sz);

    12. 
  
  12. }

    13. 
  
  13. }

    14. 
  
  14. }


      

 以上代码在非极大值抑制产生的二值灰度矩阵的潜在点中按照高阈值寻找边缘，并以所找到的点为中心寻找邻域内满足低阈值的点，从而形成一个闭合的轮廓。然后对于不满足条件的点，可用如下代码直接删除掉。


  1. //将还没有设置为边界的点设置为非边界点

    2. 
  
  2. for(i=0; i<nHeight; i++)

    3. 
  
  3. {

    4. 
  
  4. for(j=0; j<nWidth; j++)

    5. 
  
  5. {

    6. 
  
  6. if(N[i*nWidth+j] != 255)

    7. 
  
  7. {

    8. 
  
  8. N[i*nWidth+j]  = 0 ;   // 设置为非边界点

    9. 
  
  9. }

    10. 
  
  10. }

    11. 
  
  11. }


其中TraceEdge函数为一个嵌套函数，用于在每个像素点的邻域内寻找满足条件的点。其实现代码如下：


  1. void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz)

    2. 
  
  2. {

    3. 
  
  3. //对8邻域像素进行查询

    4. 
  
  4. int xNum[8] = {1,1,0,-1,-1,-1,0,1};

    5. 
  
  5. int yNum[8] = {0,1,1,1,0,-1,-1,-1};

    6. 
  
  6. LONG yy,xx,k;

    7. 
  
  7. for(k=0;k<8;k++)

    8. 
  
  8. {

    9. 
  
  9. yy = y+yNum[k];

    10. 
  
  10. xx = x+xNum[k];

    11. 
  
  11. if(pResult[yy*sz.cx+xx]==128 && pMag[yy*sz.cx+xx]>=nThrLow )

    12. 
  
  12. {

    13. 
  
  13. //该点设为边界点

    14. 
  
  14. pResult[yy*sz.cx+xx] = 255;

    15. 
  
  15. //以该点为中心再进行跟踪

    16. 
  
  16. TraceEdge(yy,xx,nThrLow,pResult,pMag,sz);

    17. 
  
  17. }

    18. 
  
  18. }

    19. 
  
  19. }


以上就从原理上实现了整个Canny算法。其检测效果如图4所示。注意：以上代码仅为作者理解所为，目的是验证本人对算法的理解，暂时没有考虑到代码的执行效率的问题。




图4 边缘检测结果


4、扩展


首先看一下OpenCV中cvCanny函数对该图像的处理结果，如图5所示。




图5 OpenCV中的Canny边缘检测结果


对比图4和图5可以发现，作者自己实现的边缘检测效果没有OpenCV的好，具体体现在：1）丢失了一些真的边缘；2）增加了一些假的边缘。


经过对整个算法的来回检查，初步推断主要的问题可能在于在进行灰度矩阵梯度幅值计算式所采用的模板算子性能不是太好，还有就是关于两个阈值的选取方法。关于这两个方面的改进研究，后文阐述。


5、总结


本文是过去一段时间，对图像边缘检测方法学习的总结。主要阐述了Canny算法的工作原理，实现过程，在此基础上基于VC6.0实现了该算法，并给出了效果图。最后，通过对比发现本文的实现方法虽然能够实现边缘检测，但效果还不是很理想，今后将在阈值选取原则和梯度幅值算子两个方面进行改进。


(MATLAB版)


我可没直接调用系统函数，要是那样就太水了。其实我的matlab代码很容易就能翻译成c/c++的。


　　canny边缘检测一共四个部分：


　　1.对原图像高斯平滑


　　2.对高斯平滑后的图像进行sobel边缘检测。这里需要求横的和竖的还有联合的，所以一共三个需要sobel边缘检测图像。


　　3.对联合的sobel检测图像进行非极大抑制


　　4.连接边缘点并进行滞后阈值处理。


下面是代码：


main.m





clear all;

close all;

clc;


img=imread(‘lena.jpg‘);

imshow(img);

[m n]=size(img);

img=double(img);


%%canny边缘检测的前两步相对不复杂，所以我就直接调用系统函数了

%%高斯滤波

w=fspecial(‘gaussian‘,[5 5]);

img=imfilter(img,w,‘replicate‘);

figure;

imshow(uint8(img))


%%sobel边缘检测

w=fspecial(‘sobel‘);

img_w=imfilter(img,w,‘replicate‘);      %求横边缘

w=w‘;

img_h=imfilter(img,w,‘replicate‘);      %求竖边缘

img=sqrt(img_w.^2+img_h.^2);        %注意这里不是简单的求平均，而是平方和在开方。我曾经好长一段时间都搞错了

figure;

imshow(uint8(img))


%%下面是非极大抑制

new_edge=zeros(m,n);

for i=2:m-1

    for j=2:n-1

        Mx=img_w(i,j);

        My=img_h(i,j);


if My~=0

            o=atan(Mx/My);      %边缘的法线弧度

        elseif My==0 && Mx>0

            o=pi/2;

        else

            o=-pi/2;

        end


%Mx处用My和img进行插值

        adds=get_coords(o);      %边缘像素法线一侧求得的两点坐标，插值需要

        M1=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %插值后得到的像素，用此像素和当前像素比较

        adds=get_coords(o+pi);　　%边缘法线另一侧求得的两点坐标，插值需要

        M2=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %另一侧插值得到的像素，同样和当前像素比较


isbigger=(Mx*img(i,j)>M1)*(Mx*img(i,j)>=M2)+(Mx*img(i,j)<M1)*(Mx*img(i,j)<=M2); %如果当前点比两边点都大置1


if isbigger

           new_edge(i,j)=img(i,j);

        end

    end

end

figure;

imshow(uint8(new_edge))


%%下面是滞后阈值处理

up=120;     %上阈值

low=100;    %下阈值

set(0,‘RecursionLimit‘,10000);  %设置最大递归深度

for i=1:m

    for j=1:n

      if new_edge(i,j)>up &&new_edge(i,j)~=255  %判断上阈值

            new_edge(i,j)=255;

            new_edge=connect(new_edge,i,j,low);

      end

    end

end

figure;

imshow(new_edge==255)






get_coords.m





function re=get_coords(angle)       %angle是边缘法线角度，返回法线前后两点

    sigma=0.000000001;

    x1=ceil(cos(angle+pi/8)*sqrt(2)-0.5-sigma);

    y1=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);

    x2=ceil(cos(angle-pi/8)*sqrt(2)-0.5-sigma);

    y2=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);

    re=[x1 y1 x2 y2];


end






connect.m





function nedge=connect(nedge,y,x,low)       %种子定位后的连通分析

    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八连通搜寻

    [m n]=size(nedge);

    for k=1:8

        yy=y+neighbour(k,1);

        xx=x+neighbour(k,2);

        if yy>=1 &&yy<=m &&xx>=1 && xx<=n

            if nedge(yy,xx)>=low && nedge(yy,xx)~=255   %判断下阈值

                nedge(yy,xx)=255;

                nedge=connect(nedge,yy,xx,low);

            end

        end

    end


end






每步运行效果：


原图


高斯模糊后


sobel边缘检测后


非极大抑制后


上阈值120，下阈值100检测结果。


其实应该还有一个sigma变量，这个是控制高斯模板用的，如果自己做模板当然需要sigma了，这里就不需要了。至于如何做高斯模板，看这里。


我主要参考了《特征提取与图像处理》这本书。
				


				
时间： 2024-11-10 19:55:04 

		


		


	

	
	

		


		
canny算法的MATLAB和C++实现（转载）的相关文章


								

		

			

                实现Canny算法原理 python
			

		

		

									

				创作很累,如果您觉得对您有帮助,请点赞支持,感谢! 一. 总的算法流程: ① 使用高斯滤波器滤波 ② 使用 Sobel 滤波器滤波获得在 x 和 y 方向上的输出,在此基础上求出边缘的强度和边缘的角度 edge 为边缘强度,tan 为边缘角度 ↑ ③ 对边缘角度进行量化处理 对边缘角度进行量化处理算法 ↑ ④ 根据边缘角度对边缘强度进行非极大值抑制(Non-maximum suppression),使图像边缘变得更细 非极大值抑制算法:0°时取(x,y).(x+1,y).(x-1,y) 中的最大			

		

	

				

		

			

                【最短路】求两点间最短路的Floyd算法及其matlab实现
			

		

		

									

				代码来源:<图论算法及其matlab实现>(北京航空航天出版社) P22 此代码返回第一个点和最后一个点之间最短路径,以及最短路径的长度. 代码如下: 1 function [P,u ] = Floyd(W) 2 %W表示权值矩阵 3 %P表示最短路径 4 %u表示最短路的权和 5 n=length(W); 6 U=W; 7 m=1; 8 9 while m<=n %判断是否满足停止条件 10 for i=1:n 11 for j=1:n 12 if U(i,j)>U(i,m)+U			

		

	

				

		

			

                【最短路】求两点间最短路径的改进的Dijkstra算法及其matlab实现
			

		

		

									

				代码来源:<图论算法及其matlab实现>(北京航空航天出版社) P18 书中提出了基于经典Dijkstra算法改进的两种算法. 其中算法Ⅱ的效率较高. 代码如下: 1 function a=Dijk(a) 2 %a(输入量)表示图的邻接矩阵 3 %a(输出量)表示所求最短路径的距离矩阵 4 5 %建立邻接矩阵,若不还是对称矩阵,则变为对称矩阵 6 n=length(a); 7 for i=2:n 8 for j=1:(i-1) 9 a(i,j)=a(j,i); 10 end 11 end 1			

		

	

				

		

			

                群智能优化算法-测试函数matlab源码
			

		

		

									

				群智能优化算法测试函数matlab源代码 global M; creatematrix(2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %画ackley图. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ackley x from[-5 5] % x=-5:0.01:5; % [x,y]=meshgrid(x); % temp1=x.^2+y.^2; % temp2=cos(2*pi*x)+cos(2*pi*y); % z=20+			

		

	

				

		

			

                de Casteljau算法的matlab实现
			

		

		

									

				一直在写c++程序,不过对于一些作图程序来说,还是MATLAB比较实在. de Casteljau算法是作贝塞尔曲线的一种高效的算法,其思想就是对[0,1]区间中所有的t,通过n个控制顶点不断递推得到一个顶点:下面是我的代码实现: function deCasteljau(P,Q) %P is 1*n matrix for X %Q is 1*n matrix for Y m=size(P); n=m(2); x=zeros(1,101); y=zeros(1,101); p=zeros(n);			

		

	

				

		

			

                【最短路】求最大可靠路的算法及其matlab实现
			

		

		

									

				内容来源:<图论算法及其matlab实现>(北京航空航天出版社) P34 [算法用途] 求图中两顶点间的最大可靠路. 代码如下: 1 %最大可靠路的算法 2 %调用Floyd文件 3 function [P p f]=p_pathf(A) 4 %f=0表示找到路,否则f=1 5 %初始化 6 7 [m n]=size(A); 8 f=0; 9 B=zeros(m,n); 10 %对原矩阵进行转换 11 for i=1:m 12 for j=1:n 13 if A(i,j)>0&A			

		

	

				

		

			

                【最短路】必须通过两个指定点的最短路径算法及其matlab实现
			

		

		

									

				代码来源:<图论算法及其matlab实现>(北京航空航天出版社) P25 首先,函数n2shorf用来计算任意两点之间最短路径长度及最短路经过的节点 需输入起点.终点 1 function [ P u] = n2shorf( W,k1,k2) 2 3 %W表示权值矩阵 4 %k1表示始点,k2表示终点 5 %P表示两顶点之间最短路,u表示最短路的距离 6 7 8 %初始化 9 n=length(W); 10 U=W; 11 m=1; 12 %Floyd算法思想 13 while m<=n			

		

	

				

		

			

                FCM算法的matlab程序（初步）
			

		

		

									

				FCM算法的matlab程序 在https://www.cnblogs.com/kailugaji/p/9648430.html文章中已经介绍了FCM算法,现在用matlab程序实现它. 作者:凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/ 1.采用iris数据库 iris_data.txt 5.1 3.5 1.4 0.2 4.9 3 1.4 0.2 4.7 3.2 1.3 0.2 4.6 3.1 1.5 0.2 5 3.6 1.4 0.2 5.4 3.			

		

	

				

		

			

                GMM算法的matlab程序（初步）
			

		

		

									

				GMM算法的matlab程序 在https://www.cnblogs.com/kailugaji/p/9648508.html文章中已经介绍了GMM算法,现在用matlab程序实现它. 作者:凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/ 1.采用iris数据库 iris_data.txt 5.1 3.5 1.4 0.2 4.9 3 1.4 0.2 4.7 3.2 1.3 0.2 4.6 3.1 1.5 0.2 5 3.6 1.4 0.2 5.4 3.