canny边缘检测_图像算法在数值计算中的应用（1）：Canny边缘检测算法

引言

有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍在广泛应用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。在直角坐标系下，求解域差分网格通常为均匀的矩形，在表达非矩形物体边界时通常需要采用坐标映射或者网格点逼近，而网格点逼近较为简单快捷。

我们对有解析式的几何边界，构造有形状信息的布尔矩阵，可以实现边界形状的识别Python数值优化：使用Euler法求解二维热传导方程2。对于一般的图形通常没有解析边界，使有限差分方法适用于一般（平面）几何的关键是要有一个“形状”矩阵，其值表示域的外部，内部和边界。而图像就是用矩阵存储的，包含许多成熟的特征识别算法，因此考虑采用图像的方法来处理，主要利用的就是边缘检测算法。

Canny边缘检测算法 由计算机科学家 John F. Canny 于 1986 年提出，主要可以分为以下几个步骤

• 图像灰度化（降维处理）
• 高斯滤波（平滑和降噪）
• 计算图像梯度值和方向
• 应用非极大值抑制NMS
• 双阈值检测确定边界

灰度化

灰度化实际上是一种降维的操作，将三个通道的像素值转换为单通道数据，减少计算量（如不进行灰度化，也可以直接进行后面步骤），后面针对单通道进行计算。

def 

高斯滤波

高斯滤波主要使图像变得平滑，减少噪声，但同时也有可能增大了边缘的宽度。其作用原理和均值滤波器类似，都是取滤波器窗口内的像素的均值作为输出，而其系数按照高斯函数离散化，具体如下：

滤波核大小为（2k+1, 2k+1）， 这里

和

从0开始计数，与Python编程时保持下标一致，而MATLAB从1开始计数

def 

构造一个3x3的kernel窗口g_kernel = gaussian_kernel(1,3)

[[0.07511361 0.1238414  0.07511361]
 [0.1238414  0.20417996 0.1238414 ]
 [0.07511361 0.1238414  0.07511361]]

接下来用离散卷积的方法进行平滑滤波，即对原图扫描计算像素平均值

卷积扫描过程 https://mlnotebook.github.io/post/CNN1/

这里采用默认步长stride=1、padding=same，表示扫描每个像素点都做平滑，输出的矩阵大小和原图像大小相同。注意卷积运算是矩阵元素逆序相乘并求和，而python中np.multiply(A,B)函数（或点乘A*B）是数组对应元素顺序相乘，要利用这种乘法需要先将卷积核kernel旋转180°再进行运算，旋转函数np.rot90(m,k=1,axes=(0,1))，k>0代表逆时针旋转90°的次数，旋转2次即可

def 

平滑后峰值减小、‘半径’变大

计算梯度

边缘就是灰度值变化较大的的像素点的集合，灰度值变化程度和方向可以用梯度来表示。图片矩阵是离散的，这就可以按照向前差分格式计算梯度（一阶导数）：

均代表一个像素单位长度，转为矩阵运算形式如下

而梯度的幅值和方向则按照下列方式计算：

梯度代码实现如下，delta取值范围只有

，即只表达边界的单方向，不区分内外侧。

def 

非极大值抑制NMS

在高斯滤波后，边缘有可能被放大了。这个步骤使用一个规则来过滤不是边缘的点，使边缘的宽度尽可能为1个像素点：遍历梯度矩阵上的所有点，并保留边缘方向上具有极大值的像素，其余点将灰度值设为0。

如果输入的是一个二值图像（黑白图像），即像素值只有0和1或者0和255，其得到的灰度矩阵极大值很均匀，则可以统一采用简单的单阈值判断

。

如果是普通灰度图像，边缘的梯度是局部极大值而非全局极大值，需要采用NMS算法（Non-Maximum Suppression）。

NMS 需要将当前的像素的梯度，与其他方向进行比较。g1，g2，g3，g4 分别是 8个相邻节点中的 4 个点，如果 g是局部最大值，g点的梯度幅值要大于梯度方向直线与 g1g2，g4g3 两个交点的梯度幅值，即大于点 dTemp1 和 dTemp2 的梯度幅值，dTemp1和dTemp可以通过线性插值求得：

weight 

当y 方向梯度值比较大时，|dy|>|dx|，weight = |dx| / |dy|，g1，g2，g3，g4如下图所示：

当前位置 g[i, j] ，g2 = g[i-1, j]；g4 = g[i+1, j];

时（左），g1 = d[i-1, j-1]，g3 = d[i+1, j+1]；

时（右），g1 = d[i-1, j-1]，g3 = d[i+1, j+1]；

|dy|>|dx|时g1~g4的位置

当x方向梯度值比较大时，|dx|>|dy|，weight = |dx| / |dy|，g1，g2，g3，g4如下图所示：

当前位置 g[i, j] ，g2 = g[i, j-1]；g4 = g[i, j+1];

时（左），g1 = d[i-1, j-1]，g3 = d[i+1, j+1]；

时（右），g1 = d[i-1, j-1]，g3 = d[i+1, j+1]；

|dx|>|dy|时g1~g4的位置

根据以上分类，代码实现如下：

def 

提取梯度并细化边缘

由于计算梯度采用了向前差分格式计算，所以局部可能存在不连续点（边缘向左或上方移动一个像素单位）。

双阈值检测

确定上下两个阀值，位于minVal之上的都可以作为候选边缘，梯度大于maxVal的任何边缘肯定是真边缘，介于minVal和maxVal之间的像素点，如果它们连接到“真边缘”像素，则它们被视为边缘的一部分，否则也会被丢弃，这样就可能提高准确度。

G节点邻域的8个节点中的最大值大于maxVal，则该节点与“真边缘”像素连通，代码实现如下：

DT 

考虑到Python语言的特性，使用双重for-loop遍历整个图片会增加耗时和不必要的节点操作，而边缘节点只占整个数组的少量区域，因此将其转换为numpy数组切片和索引的操作以加快运算效率，采用np.where得到候选边缘的索引数组，使用单层for循环遍历较小的一维数组即可

def 

对于200x200大小的单通道图片，采用for循环遍历和数组索引方式的运行时间分别为0.018s和0.002s，效率提升接近10倍。

测试结果对比

把前面的函数封装一下，跟opencv做个对比

def 

OpenCV-Python中Canny函数的原型如下，image参数即为需要处理的单通道的灰度图，threshold1、threshold2 即对应下阈值 TL 上阈值 TH 。

import 

处理同一张图片看看效果：

if 

可以看到我们实现的canny边缘检测和opencv的canny边缘检测在采用相同高斯核平滑的情况下效果十分接近，但是执行速度上有比较明显的差别。

主要体现在滤波的卷积操作部分，Python中使用双重for循环在密集计算任务时相当‘缓慢’，而opencv在底层采用C语言优化和一些矩阵式运算，在CPU线程上已经十分高效。而python利用numpy的向量化运算也能大大提高速度，参见双阈值检测中的代码，并且numpy也有GPU版本Cupy和Numba，可以使用GPU加速大型向量运算。

参考

canny 算子python实现

图像处理基础(4)：高斯滤波器详解

OpenCV-Python教程（8、Canny边缘检测）