canny算子及边缘提取原理

整理于知乎：canny中的非极大抑制

和csdn：canny算子原理

算法原理：

(1)去噪

第一步是对原始数据与高斯 mask 作卷积，得到的图像与原始图像相比有些轻微的模糊（blurred）。

(2) 用 一阶偏导的有限差分来计算梯度的幅值和方向(用一对正交的微分滤波器做卷积(如prewitt滤波)，得到包含水平和竖直方向上的导数图像H V，然后针对两图上的每个像素，对应的去求梯度方向和幅度。方向用v/h并取反正切，得到角度值，幅度则是v h的平方和开根号。) 。

(3) 对梯度幅值进行非极大值抑制。

仅仅得到全局的梯度并不足以确定边缘，因此为确定边缘，必须保留局部梯度最大的点，而抑制非极大值。（non-maxima suppression,NMS）

解决方法：利用梯度的方向。

(4) 用双阈值算法检测并连接边缘。

解决方法：双阈值算法。双阈值算法对非极大值抑制后的图象作用两个阈值τ1和τ2，且 2τ1≈τ2。从而可以 得到两个阈值边缘图象N1［i,j］和N2［i，j］(操作是：分别 <τ1，pixel=0，得到N1[i，j]； <τ2，pixel=0，得到N2[i，j])。 由于N2［i，j］使用高阈值得到，因而含有很少的假边缘，但有间断(不闭合，因为阈值可能相对太大了...)。 双阈值法要在N2［i，j］中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在N1［i,j］的8邻点位置寻找可以连接到轮廓上的边缘，这样， 算法不断地在N1［i,j］中收集边缘，直到将N2［i,j］连接起来为止( 以N2为基础，利用N1辅助，完成N2中轮廓线的连接...对着N1，把N2里面需要补全的补上.)。

现在解释一下(3)的 非最大拟制原理：

进行非最大抑制就是：寻找局部的像素最大值点，可以剔除很多非边缘点。

非最大值抑制只在0、90、45、135四个梯度方向上进行，根据gx、gy(x y 方向上的梯度，gx更大的话，那么梯度线(下图蓝色的线)就更贴近x轴..反之gy更大则贴近y轴 )  的大小比较，可以初步的确定用来插值得到dTmp1、dTmp2值的g1g2g3g4到底是哪几个点。

进行非极大值抑制，就 首先要确定像素点C的灰度值在其8值邻域内是最大的，然后再在它的梯度方向上再找是否有更大的 。图1中 蓝色的线条方向为C点的梯度方向，这样就可以确定其局部的最大值肯定分布在这条线上，也即除了C点外，梯度方向的交点dTmp1和dTmp2这两个点的值也可能会是局部最大值。 因此，判断C点灰度与这两个点的灰度大小即可判断C点是否为其邻域内的局部最大灰度点。 如果 经过判断， C点灰度值小于这两个点中的任一个，那就说明C点不是局部极大值，那么则可以排除C点为边缘。这就是非极大值抑制的工作原理。

图1 图2

在理解的过程中需要注意以下两点：

1.中非最大抑制是回答这样一个问题：“当前的梯度值在梯度方向上是一个局部最大值吗？”所以,要把当前位置的梯度值与梯度方向上两侧的梯度值进行比较。

2.梯度方向垂直于边缘方向。但实际上，我们只能得到C点邻域的8个点的值，而dTmp1和dTmp2

并不在其中，要得到这两个值就需要对dTmp1和dTmp2两点进行线性插值，也即根据图1中的g1和g2对dTmp1进行插值，根据g3和g4对dTmp2进行插值(以得到dTmp1、dTmp2两个位置处的像素值)，这要用到其梯度方向，这也是Canny算法中要求解梯度方向矩阵Thita的原因(算法的第二步)。

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

代码+注释：

void NonMaxSuppress(int*pMag,int* pGradX,int*pGradY,SIZE sz,LPBYTE pNSRst)  
{  
    LONG x,y;  
    int nPos;  
    // the component of the gradient   
    int gx,gy;  
    // the temp varialbe   
    int g1,g2,g3,g4;  
    double weight;  
    double dTemp,dTemp1,dTemp2;  
    //设置图像边缘为不可能的分界点   
    for(x=0;x<sz.cx;x++)  
    {  
        pNSRst[x] = 0;  
        pNSRst[(sz.cy-1)*sz.cx+x] = 0;  
          
    }  
    for(y=0;y<sz.cy;y++)  
    {  
        pNSRst[y*sz.cx] = 0;  
        pNSRst[y*sz.cx + sz.cx-1] = 0;  
    }  
      
    for (y=1;y<sz.cy-1;y++)  
    {  
        for (x=1;x<sz.cx-1;x++)  
        {  
            nPos=y*sz.cx+x;  
            // if pMag[nPos]==0, then nPos is not the edge point   
            if (pMag[nPos]==0)  
            {  
                pNSRst[nPos]=0;  
            }  
            else  
            {  
                // the gradient of current point   
                dTemp=pMag[nPos];  
                // x,y 方向导数   
                gx=pGradX[nPos];  
                gy=pGradY[nPos];  
                //如果方向导数y分量比x分量大，说明导数方向趋向于y分量,即更贴近y轴   
                if (abs(gy)>abs(gx))  
                {  
                    // calculate the factor of interplation   
                    weight=fabs(gx)/fabs(gy);  
                    g2 = pMag[nPos-sz.cx];  // 确定g2出现在c(中心点)的上一行   
                    g4 = pMag[nPos+sz.cx];  // 确定g4出现在c(中心点)的下一行   
                     
                    //C 为当前像素，与g1-g4 的位置关系为：   
                    //g1 g2   
                    //   C   
                    //   g4 g3   
                    if(gx*gy>0) //如果x,y两个方向导数的符号相同  
                    {  
                        g1 = pMag[nPos-sz.cx-1];  
                        g3 = pMag[nPos+sz.cx+1];  
                    }   //对的，画个图可以很好的理解。 左上角为原点，x，y导数相同即导数线跨越2 4象限 如上图1                     
                    //如果x,y两个方向的方向导数方向相反   
                    //C是当前像素，与g1-g4的关系为：   
                    //    g2 g1   
                    //    C   
                    // g3 g4   
                    else  //左上角为原点，x，y导数相反即导数线跨越1 3象限  如上图2

                    {  
                        g1 = pMag[nPos-sz.cx+1];  
                        g3 = pMag[nPos+sz.cx-1];  
                    }  
                }  
                else  
                {  
                    //插值比例   
                    weight = fabs(gy)/fabs(gx);                   
                    g2 = pMag[nPos+1]; //后一列   
                    g4 = pMag[nPos-1];  // 前一列                 
                    //如果x,y两个方向的方向导数符号相同   
                    //当前像素C与 g1-g4的关系为   
                    // g3   
                    // g4 C g2   
                    //       g1   
                    if(gx * gy > 0)  
                    {  
                        g1 = pMag[nPos+sz.cx+1];  
                        g3 = pMag[nPos-sz.cx-1];  
                    }  
                      
                    //如果x,y两个方向导数的方向相反   
                    // C与g1-g4的关系为   
                    // g1   
                    // g4 C g2   
                    //      g3   
                    else  
                    {  
                        g1 = pMag[nPos-sz.cx+1];  
                        g3 = pMag[nPos+sz.cx-1];  
                    }  
                }  
                  
                dTemp1 = weight*g1 + (1-weight)*g2;  
                dTemp2 = weight*g3 + (1-weight)*g4;               
                //当前像素的梯度是局部的最大值   
                //该点可能是边界点   
                if(dTemp>=dTemp1 && dTemp>=dTemp2)  
                {  
                    pNSRst[nPos] = 128;  
                }  
                else  
                {  
                    //不可能是边界点   
                    pNSRst[nPos] = 0;  
                }             
            }  
        }  
    }  
}  

以上代码来自：canny详解   我加了点注释而已。