1 、原理

图象边缘就是图像颜色快速变化的位置，对于灰度图像来说，也就是灰度值有明显变化的位置。图像边缘信息主要集中在高频段，图像锐化或检测边缘实质就是高通滤波。数值微分可以求变化率，在图像上离散值求梯度，图像处理中有多种边缘检测（梯度）算子，常用的包括普通一阶差分，Robert算子（交叉差分），Sobel算子，二阶拉普拉斯算子等等，是基于寻找梯度强度。

Canny 边缘检测算法是John F. Canny 于1986年开发出来的一个多级边缘检测算法，也被很多人认为是边缘检测的 最优算法, 最优边缘检测的三个主要评价标准是:

低错误率: 标识出尽可能多的实际边缘，同时尽可能的减少噪声产生的误报。

高定位性: 标识出的边缘要与图像中的实际边缘尽可能接近。

最小响应: 图像中的边缘只能标识一次。

Canny算子求边缘点具体算法步骤如下：

1. 用高斯滤波器平滑图像．

2. 用一阶偏导有限差分计算梯度幅值和方向.

3. 对梯度幅值进行非极大值抑制.

4. 用双阈值算法检测和连接边缘．

2 、实现步骤

2.1、消除噪声

使用高斯平滑滤波器卷积降噪。下面显示了一个 size = 5 的高斯内核示例:

2.2、计算梯度幅值和方向

按照Sobel滤波器的步骤，计算水平和垂直方向的差分Gx和Gy：

在vs中可以看到sobel像素值和形状：

梯度幅值和方向为：

梯度方向近似到四个可能角度之一(一般 0, 45, 90, 135)。

2.3、非极大值抑制

非极大值抑制是指寻找像素点局部最大值。sobel算子检测出来的边缘太粗了，我们需要抑制那些梯度不够大的像素点，只保留最大的梯度，从而达到瘦边的目的。沿着梯度方向，比较它前面和后面的梯度值，梯度不够大的像素点很可能是某一条边缘的过渡点，排除非边缘像素，最后保留了一些细线。

在John Canny提出的Canny算子的论文中，非最大值抑制就只是在0、90、45、135四个梯度方向上进行的，每个像素点梯度方向按照相近程度用这四个方向来代替。梯度向量的每个四分之一圆被45°线分成两种情况，一种情况是倾向于水平，另一种倾向于竖直，一共 8 个方向。这种情况下，非最大值抑制所比较的相邻两个像素就是：

1） 0：左边 和 右边

2） 45：右上 和 左下

3） 90：上边 和 下边

4）135：左上 和 右下

这样做的好处是简单，但是这种简化的方法无法达到最好的效果，因为自然图像中的边缘梯度方向不一定是沿着这四个方向的，即梯度方向的线并没有落在8邻域坐标点上。因此，就有很大的必要进行插值，找出在一个像素点上最能吻合其所在梯度方向的两侧的像素值。

如果|gx|>|gy|，这说明该点的梯度方向更靠近X轴方向，所以g2和g4则在C的左右，我们可以用下面来说明这两种情况（方向相同和方向不同）：

可以使用插值计算出真实梯度值：

其中，插值计算方式为：dTemp1 = weight*g1 + (1-weight)*g2; dTemp2 = weight*g3 + (1-weight)*g4;

Matlab使用非常有技巧的方式来计算方向，如下不仅做了dx、dy的大小判断还做了方向的判定。

witch direction    case 1        idx = find((iy<=0 & ix>-iy)  | (iy>=0 & ix<-iy));    case 2        idx = find((ix>0 & -iy>=ix)  | (ix<0 & -iy<=ix));    case 3        idx = find((ix<=0 & ix>iy) | (ix>=0 & ix<iy));    case 4        idx = find((iy<0 & ix<=iy) | (iy>0 & ix>=iy));end

2.4、双阈值检测和区域连通

最后一步，Canny 使用了滞后阈值，滞后阈值需要两个阈值(高阈值和低阈值)。如果边缘像素的梯度值高于高阈值，则将其标记为强边缘像素；如果边缘像素的梯度值小于高阈值并且大于低阈值，则将其标记为弱边缘像素；如果边缘像素的梯度值小于低阈值，则会被抑制。阈值的选择取决于给定输入图像的内容。Canny 推荐的 高:低 阈值比在 2:1 到3:1之间。

3、代码实现

3.1 计算梯度

/** Sobel 梯度计算*/Mat gradients(Mat &img, Mat &sobel){    int W = img.cols;    int H = img.rows;    Mat dx = Mat_<int>(img.size());    int border = (int)sobel.rows / 2;    for (int r = border; r < H - border; r++)    {        for (int c = border; c < W - border; c++)        {            float tmp = 0;            for (int i = -border; i <= border; i++) {                for (int j = -border; j <= border; j++) {                    tmp += (int)img.data[(r + i)*W + c + j] * sobel.at<int>(i + border, j + border);                 }            }            dx.at<int>(r, c) = tmp;        }    }    return dx;}

3.2计算非极大值抑制（详细推导过程见参考文献文章）

/*fucntion: non-maximum suppressioninput:pMag:   pointer to Magnitude,pGradX: gradient of x-directionpGradY: gradient of y-directionsz: size of pMag (width = size.cx, height = size.cy)limit: limitationoutput:pNSRst: result of non-maximum suppression*/void NonMaxSuppress(int *pMag, int * pGradX, int *pGradY, Size sz, int *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.width; x++)    {        pNSRst[x] = 0;        pNSRst[(sz.height - 1)*sz.width + x] = 0;    }    for (y = 0; y < sz.height; y++)    {        pNSRst[y*sz.width] = 0;        pNSRst[y*sz.width + sz.width - 1] = 0;    }    for (y = 1; y < sz.height - 1; y++)    {        for (x = 1; x < sz.width - 1; x++)        {            nPos = y * sz.width + 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分量                if (abs(gy) > abs(gx))                {                    // calculate the factor of interplation                    weight = fabs(gx) / fabs(gy);                    g2 = pMag[nPos - sz.width];  // 上一行                    g4 = pMag[nPos + sz.width];  // 下一行                    //如果x,y两个方向导数的符号相同                    //C 为当前像素，与g1-g4 的位置关系为：                    //g1 g2                    //   C                    //   g4 g3                    if (gx*gy > 0)                    {                        g1 = pMag[nPos - sz.width - 1];                        g3 = pMag[nPos + sz.width + 1];                    }                    //如果x,y两个方向的方向导数方向相反                    //C是当前像素，与g1-g4的关系为：                    //    g2 g1                    //    C                    // g3 g4                    else                    {                        g1 = pMag[nPos - sz.width + 1];                        g3 = pMag[nPos + sz.width - 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.width + 1];                        g3 = pMag[nPos - sz.width - 1];                    }                    //如果x,y两个方向导数的方向相反                    // C与g1-g4的关系为                    // g1                    // g4 C g2                    //      g3                    else                    {                        g1 = pMag[nPos - sz.width + 1];                        g3 = pMag[nPos + sz.width - 1];                    }                }                dTemp1 = weight * g1 + (1 - weight)*g2;                dTemp2 = weight * g3 + (1 - weight)*g4;                if(dTemp )                //当前像素的梯度是局部的最大值                //该点可能是边界点                if (dTemp >= dTemp1 && dTemp >= dTemp2)                {                    pNSRst[nPos] = dTemp;                }                else                {                    //不可能是边界点                    pNSRst[nPos] = 0;                }            }        }    }}

3.3 双阈值检测和边缘连接

void duble_threshold(Mat &pMag, Mat &pThreadImg, float threshold){    double maxv;    int * img_ptr = pMag.ptr<int>(0);    uchar * dst_ptr = pThreadImg.ptr<uchar>(0);    minMaxLoc(pMag, 0, &maxv, 0, 0);    cout << "max" << maxv << endl;    int TL = 0.333 * threshold *maxv; // 1/3 of TH    int TH = threshold *maxv;    int w = pMag.cols;    int h = pMag.rows;    for (int r = 1; r < pMag.rows; r++)    {        for (int c = 1; c < pMag.cols; c++)        {            int tmp = img_ptr[r*w + c];            if (tmp < TL) {                dst_ptr[r*w + c] = 0;            }            else if (tmp >= TH) {                dst_ptr[r*w + c] = 255;            }            else {                bool connect = false;                for(int i=-1; i<=1 && connect == false; i++)                    for (int j = -1; j <= 1 && connect == false; j++)                    {                        if (img_ptr[r + i, c + j] >= TH)                         {                            dst_ptr[r*w + c] = 255;                            connect = true;                            break;                        }                        else  dst_ptr[r*w + c] = 0;                    }            }        }    }}

测试1：左侧是原图，右侧是进行了sobel梯度计算和非极大值抑制后的图。

可见右图，在企鹅轮廓内部还有孤立的点，放大后如下图。

使用双阈值限定后如下图，内部点消失了。

测试2：

测试3：

如下图，图2的双阈值计算梯度后最大梯度360，图3使用0.5倍高阈值，轮廓不连贯，可见阈值过高。改为0.2倍高阈值，结果如图4，改善了轮廓缺失问题。