计算机视觉入门之图像处理＜七＞：图像边缘检测——Canny边缘检测

往期文章回顾：
计算机视觉入门之图像处理<一>：图像处理基础概念
计算机视觉入门之图像处理<二>：图像处理基础概念
计算机视觉入门之图像处理<三>：图像插值方法
计算机视觉入门之图像处理<四>：图像直方图均衡化
计算机视觉入门之图像处理<五>：图像平滑处理
计算机视觉入门之图像处理<六>：图像锐化处理

---

前言

边缘检测是基于灰度突变的原理，在此基础上实现图像分割等操作。本篇文章主要讨论边缘检测方法——Canny边缘检测的基本实现原理和基本代码实现。

边缘模型

边缘模型主要分为台阶边缘模型、斜坡边缘模型和屋顶边缘模型。台阶边缘模型指的是在1个像素的距离上发生两个灰度级间理想的过渡（后续的Canny边缘检测算法中就采用该边缘模型进行推导），即这两个灰度级过渡时存在一条较为明显的分界线，如下图所示：

但是在实际情况中，数字图像可能被模糊或者是带有噪声，此时的边缘模型更接近斜坡边缘模型，灰度级的过渡存在一个较宽的区域，也就不存在一个较为明显的分界线，如下图所示：

而屋顶边缘模型是通过一个区域的线的模型，屋顶边缘的基底有该线的宽度和尖锐度决定，如下图所示：

基本边缘检测

为了检测图像中的边缘，灰度的变化可采用一阶或二阶导数完成，本文中仅讨论一阶一阶导数，详细原理见计算机视觉入门之图像处理＜六＞：图像锐化处理一文。
梯度指出了灰度变化最大的方向，作为一个向量，梯度$\nabla f$的大小表示为：$$M(x,y)=\sqrt{g_x^2+g_y^2}\text{\hspace{1em}(1)}$$式(1)中$g_x$和$g_y$分别表示水平方向和竖直方向的梯度，由$M(x,y)$构成的图像称为梯度图像，其大小与原图像大小相同。
梯度$\nabla f$的方向表示为：$$\alpha (x,y)=\arctan \left[\frac{g_y}{g_x}\right]\text{\hspace{1em}(2)}$$式(2)中的$\alpha (x,y)$同样是与原图像相同的图像，因为梯度表示灰度变化最大数的方向，因此在任意点$(x,y)$处的边缘方向总是与该点处的梯度方向正交，梯度的角度信息在接下来讨论的Canny边缘检测算法中具有重要作用。

Canny边缘检测

Canny边缘检测的优点：

• 低错误率。因为一般的边缘检测算子可能存在检测到伪边缘的情况，因此Canny算法检测到的边缘尽可能地是真实的边缘。
• 较好地定位边缘点。由检测器标记的边缘点与真实边缘点中心尽可能地接近。
• 单一的边缘响应。图像中的边缘只标记出一次。

接下来将分布讨论Canny边缘检测算法的原理及实现：

1.图像灰度化

对图像灰度化处理是图像处理的基本步骤，有利于简化计算，降低模型的复杂度，主要的灰度化方法在计算机视觉入门之图像处理＜二＞：图像处理基础概念一文中有较为详细的讨论和代码实现。
灰度化效果：

2.图像滤波处理

首先是为什么要进行滤波处理？因为边缘检测中的梯度是基于微分实现的，微分对图像中的灰度突变(一般为噪声)较为敏感，对图像进行滤波有利于减少伪边缘，提高结果的准确性，因此必须对图像图像进行滤波处理。在此处的滤波处理特指高斯滤波，为了实现单一的边缘响应，需要通过数值优化的方法来得到最优解从而达到单一边缘响应的目的，为了简化计算，实践证明，该种数值近似方法可以使用高斯平滑滤波的一阶微分近似得到，将边缘检测算子推广到二维的情况，高斯滤波公式可以表示为：$$G(x,y)=e^{-\frac{x^2+y^2}{2{\sigma }^2}}\text{\hspace{1em}(3)}$$式(3)中$x$和$y$分别表示图像坐标，$\sigma$表示标准差，式(3)有时也可带上规范化因子$\frac{1}{\sqrt{2\pi }\sigma }$

详细代码请参考计算机视觉入门之图像处理＜五＞：图像平滑处理中高斯滤波代码实现

• 高斯卷积核的大小将影响Canny边缘检测的性能，卷积核尺寸越大，检测器对噪声的敏感度就越低(敏感度低就说明就可以减少伪边缘的产生)，但在一定程度上也会增大边缘的定位误差，大量实验证明，使用5*5的高斯卷积核可以达到一个较好的trade-off.

滤波效果：

3.计算梯度幅值图像和角度图像

计算梯度幅值图像采用式(1)进行计算，有时也可以将式(1)简化为$M(x,y)=|g_x|+|g_y|$，角度图像采用式(2)进行计算，在计算梯度幅值图像和角度图像时为保证图像尺度与原图像相同，需要对高斯滤波后的图像进行边缘填充。
实现代码：

    def grad_img(self, gaussed_img):
        grad_img = np.zeros_like(gaussed_img)#建立一个空矩阵用来存储梯度图像
        padding_img = cv.copyMakeBorder(gaussed_img, 1, 1, 1, 1, borderType=cv.BORDER_CONSTANT, value=0)  # zero padding
        alpha_matrix = np.zeros(gaussed_img.shape)#用于存储梯度角
        sobel_x_img = np.zeros(gaussed_img.shape)#水平方向梯度图像
        sobel_y_img = np.zeros(gaussed_img.shape)#竖直方向梯度图像
        sobel_kernel_x = np.array([[-1, -2, -1],#水平Sobel算子
                                   [0, 0, 0],
                                   [1, 2, 1]])
        sobel_kernel_y = np.array([[-1, 0, 1],#竖直Sobel算子
                                   [-2, 0, 2],
                                   [-1, 0, 1]])
        for i in range(gaussed_img.shape[0]):
            for j in range(gaussed_img.shape[1]):
                m = i + 1
                n = j + 1
                arr = padding_img[m - 1: m + 2, n - 1: n + 2]#在零填充之后的图像中取3*3的邻域
                sobel_x_img[i, j]   = np.where(np.sum(arr * sobel_kernel_x)==0, 0.00001, np.sum(arr * sobel_kernel_x))#因为计算梯度角时该项位于分母不能为零
                sobel_y_img[i, j]   = np.sum(arr * sobel_kernel_y)
                grad_img[i, j]      = np.sqrt(np.square(sobel_x_img[i, j]) + np.square(sobel_y_img[i, j]))#计算梯度图像
                alpha_matrix[i, j]  = np.arctan(sobel_y_img[i, j] / sobel_x_img[i, j])*(180 / np.pi)#计算梯度角
                alpha_matrix[i, j]  = np.where(alpha_matrix[i, j] < 0, alpha_matrix[i, j] + 360, alpha_matrix[i, j])#使得梯度角为正
        return grad_img, alpha_matrix

梯度图像：

4.对梯度图像进行非极大值抑制

在计算机视觉入门之图像处理＜六＞：图像锐化处理一文中提到，一阶微分在灰度值斜坡过渡时不为零且存在较粗的边缘，较粗的边缘会增大边缘检测的误差，因此需要细化边缘，一种较为常用的方法是非极大值抑制(Non-Maximal Suppression)，即在梯度图像中寻找梯度方向上的最大值作为边缘，不是梯度方向上的最大值则抑制为零。因为梯度方向是灰度变化最大的方向，所以在该点处梯度值(该点处的一阶微分)最大且非零，因此可以利用该特性进行非极大值抑制的求解。

一阶微分最大。比较梯度图像中每一点的灰度值与梯度方向上至少两个梯度图像像素点灰度值的大小，根据上述大小关系来确定是否保留该点的灰度值。

如上图所示，左图表示的是梯度图像的$3\ast 3$邻域，假设c点为当前梯度图像中心点，蓝色指向方向表示的是c点的梯度方向，因为在梯度方向上并不存在实际的像素点，因此需要采用插值法求出梯度方向上两个虚拟点$dTmp1$和$dTmp2$的灰度值(详细求解方法请参考计算机视觉入门之图像处理＜三＞：图像插值方法)，然后比较这三点的灰度值大小。若梯度图像中c点灰度值均大于梯度方向上两点灰度值，则保持当前灰度值不变，否则，令该点灰度值为零，从而达到抑制的目的。(上图中的$\theta$并非真正的梯度方向角，梯度角是以x轴为基准，上图中水平方向是x轴，竖直方向是y轴，而实际OpenCV操作中竖直方向是x轴，水平方向是y轴，即将水平坐标系顺时针旋转90°就是图像坐标系)
上述方法中因为涉及到虚拟像素点灰度值的计算，会使得模型变复杂，为简化计算，可以$3\ast 3$邻域分为以水平、竖直、+45°和-45°线为边界的8部分区域，上图右图所示，判断梯度角在哪个角度范围内，就将上述讨论的虚拟点灰度值设为对应点的灰度值，比如当前梯度角为10°，因为-22.5°<10°<22.5°，则两虚拟点对应得灰度值取中心点对应竖直方向上两像素点的灰度值(这里以实际图像坐标系为基准)，确定虚拟像素点灰度值后，接下来方法同上。
实现代码：

    def non_maximum_supression(self, grad_img, alpha_matrix):
        nms_img = np.zeros_like(grad_img)#存储NMS之后的图像
        padding_img = cv.copyMakeBorder(grad_img, 1, 1, 1, 1, borderType=cv.BORDER_CONSTANT, value=0)#零填充
        for i in range(grad_img.shape[0]):
            for j in range(grad_img.shape[1]):
                m = i + 1
                n = j + 1
                arr = padding_img[m - 1: m + 2, n - 1: n + 2]#3*3邻域
                if (22.5 <= alpha_matrix[i, j] <= 67.5) or (202.5 < alpha_matrix[i, j] <= 247.5):
                        nms_img[i, j] = np.where(arr[1, 1] > max(arr[0, 0], arr[2, 2]), arr[1, 1], 0)

                elif (67.5 < alpha_matrix[i, j] <= 112.5) or (247.5 < alpha_matrix[i, j] <= 292.5):
                        nms_img[i, j] = np.where(arr[1, 1] > max(arr[1, 0], arr[1, 2]), arr[1, 1], 0)

                elif (112.5 < alpha_matrix[i, j] < 157.5) or (292.5 < alpha_matrix[i, j] <= 337.5):
                        nms_img[i, j] = np.where(arr[1, 1] > max(arr[2, 0], arr[0, 2]), arr[1, 1], 0)

                else:
                        nms_img[i, j] = np.where(arr[1, 1] > max(arr[0, 1], arr[2, 1]), arr[1, 1], 0)
        return  nms_img

NMS处理后的效果：

5.双阈值处理和连接边缘

在完成非极大值抑制后，得到的图像中仍然包括由噪声或其他原因造成的伪边缘，因此需要进一步处理。所谓的双阈值处理就是根据实际情况需要设置一个灰度高阈值和一个灰度低阈值对NMS后的图像进行过滤，使得得到的边缘尽可能是真实的边缘。设经过NMS处理后的图像的当前像素点灰度值为$f$，灰度高阈值为$f_H$，灰度低阈值为$f_L$，则双阈值的基本处理步骤为：

1. $f$>$f_H$，将当前像素点标记为强边缘像素；
2. $f_L$<$f$<$f_H$，将当前像素点标记为弱边缘像素；
3. $f$<$f_L$，将当前像素点灰度置零。

经过上述操作后的强边缘像素被假设为有效边缘像素(可以将其灰度值设为255)，接下来是对弱边缘像素的处理，假如弱边缘像素点的邻域(该邻域可以是4-邻域也可以是8-邻域)内存在强边缘像素点，则也将其标记(同样可以将其灰度值设为255)，否则就将其灰度置零。

• 边缘检测的效果较大程度上受双阈值大小的影响，因此应该根据实际需要进行设置。

实现代码：

    def threshold_process(self, nms_img, low_threshold, high_threshold):
    #这里设置的低阈值为grad_img.mean() * 0.5， 高阈值为低阈值的3倍
        flag_matrix = np.zeros(nms_img.shape)#建立一个标记矩阵，与NMS处理之后的图像同大小
        threshold_img = np.zeros_like(nms_img)#用于存储阈值处理之后的图像
        for i in range(nms_img.shape[0]):
            for j in range(nms_img.shape[1]):
                if nms_img[i, j] >= high_threshold:
                    flag_matrix[i, j] = 1#标记为强边缘像素点
                elif (nms_img[i, j] < high_threshold) and (nms_img[i, j] >= low_threshold):
                    flag_matrix[i, j] = 0.5#标记为弱边缘像素点
                else:
                    flag_matrix[i, j] = 0#非边缘像素点置零
        padding_flag_matrix = cv.copyMakeBorder(flag_matrix, 1, 1, 1, 1, borderType=cv.BORDER_CONSTANT, value=0)

        for i in range(flag_matrix.shape[0]):#根据标记矩阵来进一步确定弱边缘像素是否有效
            for j in range(flag_matrix.shape[1]):
                m = i + 1
                n = j + 1
                arr = padding_flag_matrix[m - 1: m + 2, n - 1: n + 2]  # 3*3邻域
                if padding_flag_matrix[m, n] == 1:
                    threshold_img[i, j] = 255
                elif padding_flag_matrix[m, n] == 0.5:
                    if 1 in arr:#判断邻域像素矩阵内是否有强边缘像素点
                        threshold_img[i, j] = 255
                    else:
                        threshold_img[i, j] = 0
                else:
                    threshold_img[i, j] = 0
        return threshold_img

边缘检测效果：

Canny边缘检测步骤总结：

1. 图像灰度化；
2. 图像滤波处理；
3. 计算梯度幅值图像和角度图像；
4. 对梯度图像进行非极大值抑制；
5. 双阈值处理和连接边缘。

[Until：2021年3月3日，…………]