Canny多级边缘检测算法的C语言实现

1 Canny算法的实现流程
1.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加权平均，权值参考OpenCV
12.{
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.2 图像的高斯滤波
根据上面所讲的边缘检测过程，下一个步骤就是对图像进行高斯滤波。可根据之前博文描述的方法获取一维或者二维的高斯滤波核。因此进行图像高斯滤波可有两种实现方式，以下具体进行介绍。

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

[cpp] view plaincopyprint?

1.double nSigma = 0.4; //定义高斯函数的标准差
2.int nWidowSize = 1+2*ceil(3*nSigma); //定义滤波窗口的大小
3.int nCenter = (nWidowSize)/2; //定义滤波窗口中心的索引
double nSigma = 0.4; //定义高斯函数的标准差 int nWidowSize = 1+2*ceil(3*nSigma); //定义滤波窗口的大小 int nCenter = (nWidowSize)/2; //定义滤波窗口中心的索引

两种方法都需要用到的变量：

[cpp] view plaincopyprint?

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.}
int nWidth = OpenCvGrayImage->width; //获取图像的像素宽度 int nHeight = OpenCvGrayImage->height; //获取图像的像素高度 unsigned char* nImageData = new unsigned char[nWidth*nHeight]; //暂时保存图像中的数据 unsigned char*pCanny = new unsigned char[nWidth*nHeight]; //为平滑后的图像数据分配内存 double* nData = new double[nWidth*nHeight]; //两次平滑的中间数据 for(int j=0; j<nHeight; j++) //获取数据 { for(i=0; i<nWidth; i++) nImageData[j*nWidth+i] = (unsigned char)OpenCvGrayImage->imageData[j*nWidth+i]; }

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

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

[cpp] view plaincopyprint?

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.}
//////////////////////生成一维高斯滤波系数///////////////////////////// double* pdKernal_1 = new double[nWidowSize]; //定义一维高斯核数组 double dSum_1 = 0.0; //求和，用于进行归一化 ////////////////////////一维高斯函数公式////////////////////////////// //// x*x ///////////////// //// -1*---------------- ///////////////// //// 1 2*Sigma*Sigma ///////////////// //// ------------ e

///////////////// //// ///////////////// //// \/2*pi*Sigma ///////////////// ////////////////////////////////////////////////////////////////////// for(int i=0; i<nWidowSize; i++) { double n
Dis = (double)(i-nCenter); pdKernal_1[i] = exp(-(0.5)*nDis*nDis/(nSigma*nSigma))/(sqrt(2*3.14159)*nSigma); dSum_1 += pdKernal_1[i]; } for(i=0; i<nWidowSize; i++) { pdKernal_1[i] /= dSum_1; //进行归一化 }

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

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.}
for(i=0; i<nHeight; i++) //进行x向的高斯滤波(加权平均) { for(j=0; j<nWidth; j++) { double dSum = 0; double dFilter=0; //滤波中间值 for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++) { if((j+nLimit)>=0 && (j+nLimit) < nWidth ) //图像不能超出边界 { dFilter += (double)nImageData[i*nWidth+j+nLimit] * pdKernal_1[nCenter+nLimit]; dSum += pdKernal_1[nCenter+nLimit]; } } nData[i*nWidth+j] = dFilter/dSum; } } for(i=0; i<nWidth; i++) //进行y向的高斯滤波(加权平均) { for(j=0; j<nHeight; j++) { double dSum = 0.0; double dFilter=0; for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++) { if((j+nLimit)>=0 && (j+nLimit) < nHeight) //图像不能超出边界 { dFilter += (double)nData[(j+nLimit)*nWidth+i] * pdKernal_1[nCenter+nLimit]; dSum += pdKernal_1[nCenter+nLimit]; } } pCanny[j*nWidth+i] = (unsigned char)(int)

dFilter/dSum; } }

1.2.2 根据二维高斯核进行滤波
1）生成二维高斯滤波系数

[cpp] view plaincopyprint?

1.//////////////////////生成一维高斯滤波系数//////////////////////////////////
2.double* pdKernal_2 = new double[nWidowSize*nWidowSize]; //定义一维高斯核数组
3.double d
Sum_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.}
//////////////////////生成一维高斯滤波系数////////////////////////////////// double* pdKernal_2 = new double[nWidowSize*nWidowSize]; //定义一维高斯核数组 double dSum_2 = 0.0; //求和，进行归一化 ///////////////////////二维高斯函数公式//////////////////////////////////// //// x*x+y*y /////////////// //// -1*-------------- /////////////// //// 1 2*Sigma*Sigma /////////////// //// ---------------- e /////////////// //// 2*pi*Sigma*Sigma /////////////// /////////////////////////////////////////////////////////////////////////// for(i=0; i<nWidowSize; i++) { for(int j=0; j<nWidowSize; j++) { int nDis_x = i-nCenter; int nDis_y = j-nCenter; pdKernal_2[i+j*nWidowSize]=exp(-(1/2)*(nDis_x*nDis_x+nDis_y*nDis_y) /(nSigma*nSigma))/(2*3.1415926*nSigma*nSigma); dSum_2 += pdKernal_2[i+j*nWidowSize]; } } for(i=0; i<nWidowSize; i++) { for(int j=0; j<nWidowSize; j++) //进行归一化 { pdKernal_2[i+j*nWidowSize] /= dSum_2; } }
2）采用高斯核进行高斯滤波，滤波后的数据保存在矩阵pCanny中

[cpp] view plaincopyprint?

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)*n
WidowSize+(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.}
int x; int y; for(i=0; i<nHeight; i++) { for(j=0; j<nWidth; j++) { double dFilter=0.0; double dSum = 0.0; for(x=(-nCenter); x<=nCenter; x++) //行 { for(y=(-nCenter); y<=nCenter; y++) //列 { if( (j+x)>=0 && (j+x)<nWidth && (i+y)>=0 && (i+y)<nHeight) //判断边缘 { dFilter += (double)nImageData [(i+y)*nWidth + (j+x)] * pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)]; dSum += pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)]; } } } pCanny[i*nWidth+j] = (unsigned char)dFilter/dSum; } }

1.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(T

heta[i*nWidth+j] < 0)
26. Theta[i*nWidth+j] += 360; //将这个角度转换到0~360范围
27. }
28.}
//////////////////同样可以用不同的检测器///////////////////////// ///// P[i,j]=(S[i,j+1]-S[i,j]+S[i+1,j+1]-S[i+1,j])/2 ///// ///// Q[i,j]=(S[i,j]-S[i+1,j]+S[i,j+1]-S[i+1,j+1])/2 ///// ///////////////////////////////////////////////////////////////// double* P = new double[nWidth*nHeight]; //x向偏导数 double* Q = new double[nWidth*nHeight]; //y向偏导数 int* M = new int[nWidth*nHeight];
//梯度幅值 double* Theta = new double[nWidth*nHeight]; //梯度方向 //计算x,y方向的偏导数 for(i=0; i<(nHeight-1); i++) { for(j=0; j<(nWidth-1); j++) { 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; 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; } } //计算梯度幅值和梯度的方向 for(i=0; i<nHeight; i++) { for(j=0; j<nWidth; j++) { M[i*nWidth+j] = (int)(sqrt(P[i*nWidth+j]*P[i*nWidth+j] + Q[i*nWidth+j]*Q[i*nWidth+j])+0.5); Theta[i*nWidth+j] = atan2(Q[i*nWidth+j], P[i*nWidth+j]) * 57.3; if(Theta[i*nWidth+j] < 0) Theta[i*nWidth+j] += 360; //将这个角度转换到0~360范围 } }

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

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

[cpp] view plaincopyprint?

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; //插值的权重
unsigned char* N = new unsigned char[nWidth*nHeight]; //非极大值抑制结果 int g1=0, g2=0, g3=0, g4=0; //用于进行插值，得到亚像素点坐标值 double dTmp1=0.0, dTmp2=0.0; //保存两个亚像素点插值得到的灰度数据 double dWeight=0.0; //插值的权重 其次，对边界进行初始化：

[cpp] view plaincopyprint?

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.}
for(i=0; i<nWidth; i++) { N[i] = 0; N[(nHeight-1)*nWidth+i] = 0; } for(j=0; j<nHeight; j++) { N[j*nWidth] = 0; N[j*nWidth+(nWidth-1)] = 0; } 进行局部最大值寻找，根据

上文图1所述的方案进行插值，然后判优，实现代码如下：

[cpp] view plaincopyprint?

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，则不是局部最大对该点赋为0
8. else
9. {
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. else
81. N[nPointIdx] = 0;
82. }
83.}
for(i=1; i<(nWidth-1); i++) { for(j=1; j<(nHeight-1); j++) { int nPointIdx = i+j*nWidth; //当前点在图像数组中的索引值 if(M[nPointIdx] == 0) N[nPointIdx] = 0; //如果当前梯度幅值为0，则不是局部最大对该点赋为0 else { ////////首先判断属于那种情况，然后根据情况插值/////// ////////////////////第一种情况/////////////////////// ///////// g1 g2 ///////////// ///////// C ///////////// ///////// g3
g4 ///////////// ///////////////////////////////////////////////////// if( ((Theta[nPointIdx]>=90)&&(Theta[nPointIdx]<135)) || ((Theta[nPointIdx]>=270)&&(Theta[nPointIdx]<315))) { //////根据斜率和四个中间值进行插值求解 g1 = M[nPointIdx-nWidth-1]; g2 = M[nPointIdx-nWidth]; g3 = M[nPointIdx+nWidth]; g4 = M[nPointIdx+nWidth+1]; dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]); //反正切 dTmp1 = g1*dWeight+g2*(1-dWeight); dTmp2 = g4*dWeight+g3*(1-dWeight); } ////////////////////第二种情况///////////

//////////// ///////// g1 ///////////// ///////// g2 C g3 ///////////// ///////// g4 ///////////// ///////////////////////////////////////////////////// else if( ((Theta[nPointIdx]>=135)&&(Theta[nPointIdx]<180)) || ((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360))) { g1 = M[nPointIdx-nWidth-1]; g2 = M[nPointIdx-1]; g3 = M[nPointIdx+1]; g4 = M[nPointIdx+nWidth+1]; dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]); //正切 dTmp1 = g2*dWeight+g1*(1-dWeight); dTmp2 = g4*dWeight+g3*(1-dWeight); } ////////////////////第三种情况/////////////////////// ///////// g1 g2 ///////////// ///////// C ///////////// ///////// g4 g3 ///////////// ///////////////////////////////////////////////////// else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) || ((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270))) { g1 = M[nPointIdx-nWidth]; g2 = M[nPointIdx-nWidth+1]; g3 = M[nPointIdx+nWidth]; g4 = M[nPointIdx+
nWidth-1]; dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]); //反正切 dTmp1 = g2*dWeight+g1*(1-dWeight); dTmp2 = g3*dWeight+g4*(1-dWeight); } ////////////////////第四种情况/////////////////////// ///////// g1 ///////////// ///////// g4 C g2 ///////////// ///////// g3 ///////////// ///////////////////////////////////////////////////// else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) || ((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225))) { g1 = M[nPointIdx-nWidth+1]; g2 = M[nPointIdx+1]; g3 = M[nPointIdx+nWidth-1]; g4 = M[nPointIdx-1]; dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]); //正切 dTmp1 = g1*dWeight+g2*(1-dWeight); dTmp2 = g3*dWeight+g4*(1-dWeight); } } //////////进行局部最大值判断，并写入检测结果//////////////// if((M[nPointIdx]>=dTmp1) && (M[nPointIdx]>=dTmp2)) N[nPointIdx] = 128; else N[nPointIdx] = 0; } }

1.5双阈值检测实现
1）定义相应参数如下
[cpp] view plaincopyprint?

1.int nHist[1024];
2.int nEdgeNum; //可能边界数
3.int nMaxMag = 0; //最大梯度数
4.int nHighCount;
int nHist[1024]; int nEdgeNum; //可能边界数 int nMaxMag = 0; //最大梯度数 int nHighCount;
2）构造灰度图的统计直方图，根据上文梯度幅值的计算公式可知，最大的梯度幅值为：



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

[cpp] view plaincopyprint?

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.}
for(i=0;i<1024;i++) nHist[i] = 0; for(i=0; i<nHeight; i++) { for(j=0; j<nWidth; j++) { if(N[i*nWidth+j]==128) nHist[M[i*nWidth+j]]++; } }

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

[cpp] view plaincopyprint?

1.nEdgeNum = n

Hist[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.}
nEdgeNum = nHist[0]; nMaxMag = 0; //获取最大的梯度值 for(i=1; i<1024; i++) //统计经过“非最大值抑制”后有多少像素 { if(nHist[i] != 0) //梯度为0的点是不可能为边界点的 { nMaxMag = i; } nEdgeNum += nHist[i]; //经过non-maximum suppression后有多少像素 }
4）计算两个阈值

[cpp] view plaincopyprint?

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); //低阈值
double dRatHigh = 0.79; double dThrHigh; double dThrLow; double dRatLow = 0.5; nHighCount = (int)(dRatHigh * nEdgeNum + 0.5); j=1; nEdgeNum = nHist[1]; while((j<(nMaxMag-1)) && (nEdgeNum < nHighCount)) { j++; nEdgeNum += nHist[j]; } dThrHigh = j; //高阈值 dThrLow = (int)((dThrHigh) * dRatLow + 0.5); //低阈值
这段代码的意思是，按照灰度值从低到高的顺序，选取前79%个灰度值中的最大的灰度值为高阈值，低阈值大约为高阈值的一半。这是根据经验数据的来的，至于更好地参数选取方法，作者后面会另文研究。
5）进行边缘检测

[cpp] view plaincopyprint?

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.}
SIZE sz; sz.cx = nWidth; sz.cy = nHeight; for(i=0; i<nHeight; i++) { for(j=0; j<nWidth; j++) { if((N[i*nWidth+j]==128) && (M[i*nWidth+j] >= dThrHigh)) { N[i*nWidth+j] = 255; TraceEdge(i, j, dThrLow, N, M, sz); } } }



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

[cpp] view plaincopyprint?

1.//将还没有设置为边界的点设置为非边界点
2.for(i=0; i<nHeight; i++)

3.{
4. for(j=0; j<nWidth; j++)
5. {
6. if(N[i*nWidth+j] != 255)
7. {
8. N[i*nWidth+j] = 0 ; // 设置为非边界点
9. }
10. }
11.}
//将还没有设置为边界的点设置为非边界点 for(i=0; i<nHeight; i++) { for(j=0; j<nWidth; j++) { if(N[i*nWidth+j] != 255) { N[i*nWidth+j] = 0 ; // 设置为非边界点 } } }

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

[cpp] view plaincopyprint?

1.void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz)
2.{
3. //对8邻域像素进行查询
4. int xNum[8] = {1,1,0,-1,-1,-1,0,1};
5. int yNum[8] = {0,1,1,1,0,-1,-1,-1};
6. LONG yy,xx,k;
7. for(k=0;k<8;k++)
8. {
9. yy = y+yNum[k];
10. xx = x+xNum[k];
11. if(pResult[yy*sz.cx+xx]==128 && pMag[yy*sz.cx+xx]>=nThrLow )
12. {
13. //该点设为边界点
14. pResult[yy*sz.cx+xx] = 255;
15. //以该点为中
心再进行跟踪
16. TraceEdge(yy,xx,nThrLow,pResult,pMag,sz);
17. }
18. }
19.}
void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz) { //对8邻域像素进行查询 int xNum[8] = {1,1,0,-1,-1,-1,0,1}; int yNum[8] = {0,1,1,1,0,-1,-1,-1}; LONG yy,xx,k; for(k=0;k<8;k++) { yy = y+yNum[k]; xx = x+xNum[k]; if(pResult[yy*sz.cx+xx]==128 && pMag[yy*sz.cx+xx]>=nThrLow ) { //该点设为边界点 pResult[yy*sz.cx+xx] = 255; //以该点为中心再进行跟踪 TraceEdge(yy,xx,nThrLow,pResult,pMag,sz); } } }
以上就从原理上实现了整个Canny算法。其检测效果如图4所示。注意：以上代码仅为作者理解所为，目的是验证本人对算法的理解，暂时没有考虑到代码的执行效率的问题。