Canny边缘检测算法原理及其VC实现详解(三)

转载自：http://blog.csdn.net/likezhaobin/article/details/6892629

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

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

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

3.4 非极大值抑制

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

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

 

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所述的方案进行插值，然后判优，实现代码如下：

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