1 、原理
图象边缘就是图像颜色快速变化的位置,对于灰度图像来说,也就是灰度值有明显变化的位置。图像边缘信息主要集中在高频段,图像锐化或检测边缘实质就是高通滤波。数值微分可以求变化率,在图像上离散值求梯度,图像处理中有多种边缘检测(梯度)算子,常用的包括普通一阶差分,Robert算子(交叉差分),Sobel算子,二阶拉普拉斯算子等等,是基于寻找梯度强度。
Canny 边缘检测算法是John F. Canny 于1986年开发出来的一个多级边缘检测算法,也被很多人认为是边缘检测的 最优算法, 最优边缘检测的三个主要评价标准是:
低错误率: 标识出尽可能多的实际边缘,同时尽可能的减少噪声产生的误报。
高定位性: 标识出的边缘要与图像中的实际边缘尽可能接近。
最小响应: 图像中的边缘只能标识一次。
Canny算子求边缘点具体算法步骤如下:
1. 用高斯滤波器平滑图像.
2. 用一阶偏导有限差分计算梯度幅值和方向.
3. 对梯度幅值进行非极大值抑制.
4. 用双阈值算法检测和连接边缘.
2 、实现步骤
2.1、消除噪声
使用高斯平滑滤波器卷积降噪。下面显示了一个 size = 5 的高斯内核示例:
![44d433c2c7fdc61da9c935eddb13a15d.png](https://i-blog.csdnimg.cn/blog_migrate/534ae696ad2dd79ba683c3baef5e23b7.jpeg)
2.2、计算梯度幅值和方向
按照Sobel滤波器的步骤,计算水平和垂直方向的差分Gx和Gy:
![64f0821e94e3e7fb8a9c49a3943e183c.png](https://i-blog.csdnimg.cn/blog_migrate/75b454c2fecff17996a75bbd4db1f648.jpeg)
在vs中可以看到sobel像素值和形状:
![70c89b0fccf717614ec9091962910107.png](https://i-blog.csdnimg.cn/blog_migrate/990304fa9b4ee37f0f2d8ff4b4530e79.jpeg)
梯度幅值和方向为:
![d7212ce76ccda85614cf0127abe32a7c.png](https://i-blog.csdnimg.cn/blog_migrate/dc232bac7bd73929e4d379a959a89f6f.jpeg)
梯度方向近似到四个可能角度之一(一般 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的左右,我们可以用下面来说明这两种情况(方向相同和方向不同):
![2d57ca85f3188eb0fba2a14302845c1c.png](https://i-blog.csdnimg.cn/blog_migrate/1296cb1c09d2bc7126cf7d12e2b368b5.jpeg)
可以使用插值计算出真实梯度值:
![e5c24c6f09cc94122971966f71b3a6ea.png](https://i-blog.csdnimg.cn/blog_migrate/e8b76c9bc3758c31e1c588430b43ba7b.jpeg)
其中,插值计算方式为:dTemp1 = weight*g1 + (1-weight)*g2; dTemp2 = weight*g3 + (1-weight)*g4;
![7a9b3bba38f63b4609a16cf3396879b4.png](https://i-blog.csdnimg.cn/blog_migrate/e1d3894368b356193a53cf6d31c4ec7b.jpeg)
Matlab使用非常有技巧的方式来计算方向,如下不仅做了dx、dy的大小判断还做了方向的判定。
witch direction case 1 idx = find((iy<=0 & ix>-iy) | (iy>=0 & ix0 & -iy>=ix) | (ix<0 & -iy<=ix)); case 3 idx = find((ix<=0 & ix>iy) | (ix>=0 & ix0 & 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_(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(i + border, j + border); } } dx.at(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(0); uchar * dst_ptr = pThreadImg.ptr(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; } } } }}
4 、测试结论
测试1:左侧是原图,右侧是进行了sobel梯度计算和非极大值抑制后的图。
![7c2ac0882e1785bf3de9e4b59ce0110a.png](https://i-blog.csdnimg.cn/blog_migrate/41bbbaa6c98428b4bdbc9df4d238b242.jpeg)
可见右图,在企鹅轮廓内部还有孤立的点,放大后如下图。
![a87bb18bc193943ee3df57b796f8fe2d.png](https://i-blog.csdnimg.cn/blog_migrate/ed968fd280ae658faf9ca655eea6f39b.jpeg)
使用双阈值限定后如下图,内部点消失了。
![fc5c1a6d06e3fa5e0b3fafb6d52d4e20.png](https://i-blog.csdnimg.cn/blog_migrate/e9589c62980cbebcc192f5227583a970.jpeg)
测试2:
![250c721568333dd2ead4eb2174ac199a.png](https://i-blog.csdnimg.cn/blog_migrate/00e40f67df9d1210b121878db6f14182.jpeg)
测试3:
如下图,图2的双阈值计算梯度后最大梯度360,图3使用0.5倍高阈值,轮廓不连贯,可见阈值过高。改为0.2倍高阈值,结果如图4,改善了轮廓缺失问题。