【CV】Canny Edge Detection

已实现的功能简述及运行简要说明

功能简述

对输入的彩色图像实现Canny Edge检测算法，自己写函数实现图像梯度计算、非极大抑制、双阈值等步骤，显示中间处理结果及最终检测结果，并将结果输出成图像文件。最终边缘结果图覆盖在原彩色图像后，提取彩色边缘图。

运行简要说明

原始图像在img_path变量中存储，将在origin窗口中进行显示；高斯平滑滤波的结果在gaussian窗口中显示，在相同目录下存储为gaussian.jpg；图像梯度幅值结果在gradient窗口中进行显示，在相同目录下存储为gradient.jpg；非极大值抑制结果在nms窗口中显示，在相同目录下存储为nms.jpg；双阈值处理结果在dt窗口中显示，存储为dt.jpg；最终提取得到的彩色边缘图在edge窗口中显示，并存储为edge.jpg。

开发与运行环境

操作系统：Windows 10，64位
开发环境：Python 3.7.3
库环境：numpy 1.16.4, cv2 4.1.1

算法基本思路、原理

高斯滤波器平滑图像

Canny Edge 检测算法的第一个步骤是使用高斯滤波器进行滤波平滑操作，此步骤的目的是对原始图片进行模糊处理，减少原始图片的噪声，使得边缘信息更为明确。高斯滤波器是对连续高斯函数的离散近似，对高斯曲面进行离散采样和归一化得出。通过将原始图像与高斯滤波器进行卷积可以对图像实现高斯滤波平滑。

二维高斯函数如下所示：
$$
H(x,y)=e^{-\frac{x^2+y^2}{2\sigma^2}}
$$
对高斯函数进行离散化操作之后，得到 $(2k+1)\times(2k+1)$ 滤波器的计算公式如下：
$$
H[i,j] = \frac{1}{2\pi\sigma^2}e^{-\frac{(i-k-1)^2+(j-k-1)^2}{2\sigma^2}}
$$

边缘检测

边缘检测的一种方式是，用一阶偏导有限差分计算梯度幅值和方向检测边缘。在图像的边缘处会出现强度快速变化，因此可以通过计算图像的梯度，找出其中幅值来识别图像的边缘。由于图片信息是离散的，可以通过在单位像素点上图片强度的变化率对梯度进行离散近似。

对梯度进行离散近似之后的计算公式如下，分别计算x方向和y方向上的梯度：
$$
\frac{\partial f}{\partial x} =（\frac{f(x_{n+1},y_n)-f(x_n,y_n)}{\Delta x}+\frac{f(x_{n+1},y_{n+1})-f(x_n,y_{n+1})}{\Delta x}）/2
$$

$$
\frac{\partial f}{\partial y} =（\frac{f(x_n,y_{n+1})-f(x_n,y_n)}{\Delta y}+\frac{f(x_{n+1},y_{n+1})-f(x_{n+1},y_n)}{\Delta y}）/2
$$

计算梯度幅值的公式为：
$$
M=\sqrt{\frac{\partial f}{\partial x }^2 + \frac{\partial f}{\partial y} ^2}
$$
梯度方向，即函数最大变化率方向：
$$
\theta = arctan^{-1}(\frac{\partial f}{\partial y}/\frac{\partial f}{\partial x})
$$
另一种方式是使用sobel、prewitt等算子进行卷积计算。

非极大抑制（NMS）

去除幅值局部变化中非极大的点，使边缘变细。

在查阅资料以及学习的过程中发现有两种方式可以对梯度幅值进行非极大抑制。

第一种方式是实现方向角的离散化，即将当前坐标点的梯度方向（$\theta$）根据所在的区域离散化为0, 45, 90, 135四个值。如图所示：

比较梯度方向离散化后所在方向对应相邻的像素点与当前像素点的梯度幅值，若当前像素点较大则保留，否则丢弃。

第二种方式是在梯度方向上假设存在额外的像素点，如图所示：

假设a,b是额外生成的像素，可通过线性插值得到a,b的颜色信息。

例如对于b：
$$
tan \theta = \frac{Gy}{Gx}
$$

$$
d_1 = \frac{Gx}{Gy}, \;\; d_2 = 1-d_1
$$

$$
b=\frac{d_1}{d_1+d_2} \times g3 + \frac{d_2}{d_1+d_2} \times g4
$$

比较当前像素点c与新像素点a,b的梯度关系，若当前像素点比新的两个像素点都大则认为是极值，保留，否则丢弃。

双阈值检测与边缘连接

设定两个边缘阈值TL和TH，梯度大于TH的被认为是真边缘，保留；低于TL的被认为非边缘，丢弃；介于两者之间的则根据连通性进行判断，若其与真边缘联通，则认为是边缘，否则丢弃。

提取彩色边缘

读取生成的边缘图像信息，若为边缘，则显示原图像的像素点，否则不显示。最终得到边缘结果覆盖在原彩色图像上得到的彩色边缘图。

具体实现

高斯平滑滤波

使用参数为1.4，大小为5*5的高斯滤波器进行滤波，先使用离散化近似的高斯公示计算高斯核每一位置的值，然后进行归一化处理，将处理完成的高斯核与原图像进行卷积，得到模糊处理后的图像。

实现方式如下：

for i in range(5):
    for j in range(5):
        g_ker[i,j] = np.exp(-1*(np.square(i-3)+np.square(j-3))/(2*np.square(sigma)))/(2*np.pi*np.square(sigma))
        g_sum = g_sum + g_ker[i,j]
g_ker = g_ker / g_sum #归一化
g_img = cv2.filter2D(img, -1, g_ker) #卷积

边缘检测

在实验过程中，尝试使用了sobel算子卷积和梯度幅值计算两种方式进行边缘检测。

sobel算子卷积的实现方法如下：

#sobel算子
sobelker_x = np.array([[-1,0,1], [-2,0,2], [-1,0,1]])
sobelker_y = np.array([[-1,-2,-1], [0,0,0], [1,2,1]])
#卷积
grad_x = cv2.convertScaleAbs(cv2.filter2D(img, -1, sobelker_x))
grad_y = cv2.convertScaleAbs(cv2.filter2D(img, -1, sobelker_y))
#合成
grad_img = cv2.addWeighted(grad_x, 0.5, grad_y, 0.5, 0)

一阶偏导有限差分计算梯度幅值及梯度方向的实现方法如下，为了方便后续处理，将角度均转为0-180范围内的角度值：

for i in range(width-1):
    for j in range(height-1):
        gx[i,j] = (int(g_img[i,j]) - int(g_img[i+1,j]) + int(g_img[i,j+1]) - int(g_img[i+1,j+1])) / 2
        gy[i,j] = (int(g_img[i,j]) - int(g_img[i,j+1]) + int(g_img[i+1,j]) - int(g_img[i+1,j+1])) / 2
        m[i,j] = np.sqrt(np.square(gx[i,j])+np.square(gy[i,j]))
        if gx[i,j] != 0:
            theta[i,j] = np.arctan(gy[i,j] / gx[i,j]) * 180 / np.pi
            if(theta[i,j] < 0):
                theta[i,j] = theta[i,j] + 180
        else:
            theta[i,j] = 90

应用于后续操作之后，发现使用slbel算子的效果不如梯度幅值。并且由于下一步可能需要用到梯度方向，因此在最终的项目种保留梯度幅值计算检测边缘。

非极大抑制（NMS）

尝试了在算法原理部分所述的两种实现方式，实现方法分别如下：

角度离散化：

for i in range(1,width - 2):
    for j in range(1, height - 2):
        if ( theta[i,j] < 22.5 or (theta[i,j] >= 157.5 and theta[i,j] < 180) ):
            theta[i,j] = 0
        elif ( theta[i,j] >= 22.5 and theta[i,j] < 67.5 ):
            theta[i,j] = 45
        elif ( theta[i,j] >= 67.5 and theta[i,j] < 112.5 ) :
            theta[i,j] = 90
        elif ( theta[i,j] >= 112.5 and theta[i,j] < 157.5):
            theta[i,j] = 135

for i in range(1, width-2):
    for j in range(1, height-2):
        if theta[i,j] == 0 and m[i,j] == np.max([m[i,j],m[i+1,j],m[i-1,j]]):
            nms_img[i,j] = m[i,j]
        if theta[i,j] == 45 and m[i,j] == np.max([m[i,j],m[i-1,j-1],m[i+1,j+1]]):
            nms_img[i,j] = m[i,j]
        if theta[i,j] == 90 and  m[i,j] == np.max([m[i,j],m[i,j+1],m[i,j-1]]):
            nms_img[i,j] = m[i,j]
        if theta[i,j] == 45 and m[i,j] == np.max([m[i,j],m[i-1,j+1],m[i+1,j-1]]):
            nms_img[i,j] = m[i,j]

插值生成额外像素点：

for i in range(1, width-2):
    for j in range(1, height-2):
        if theta[i,j] <= 45:
            if gx[i,j] != 0:
                w = np.abs(gy[i,j]) / np.abs(gx[i,j])
            g1 = (1-w) * m[i-1, j-1] + w * m[i, j-1]
            g2 = (1-w) * m[i+1, j+1] + w * m[i, j+1]
        elif theta[i,j] > 45 and theta[i,j] <= 90:
            if gy[i,j] != 0:
                w = np.abs(gx[i,j]) / np.abs(gy[i,j])
            g1 = (1-w) * m[i-1, j-1] + w * m[i-1, j]
            g2 = (1-w) * m[i+1, j+1] + w * m[i+1, j]
        elif theta[i,j] > 90 and theta[i,j] <= 135:
            if gy[i,j] != 0:
                w = np.abs(gx[i,j]) / np.abs(gy[i,j])
            g1 = (1-w) * m[i-1, j+1] + w * m[i-1, j]
            g2 = (1-w) * m[i+1, j-1] + w * m[i+1, j]
        else:
            if gx[i,j] != 0:
                w = np.abs(gy[i,j]) / np.abs(gx[i,j])
            g1 = (1-w) * m[i+1, j-1] + w * m[i, j-1]
            g2 = (1-w) * m[i-1, j+1] + w * m[i, j+1]
        if m[i,j] > g1 and m[i,j] > g2:
            nms_img[i,j] = m[i,j]

试验发现两种方法效果差不多，在最终项目中保留了插值的结果。

双阈值检测与边缘连接

进行尝试之后，设定的两个阈值分别是0.1*np.max(nms_img)和0.19*np.max(nms_img)，阈值调试的过程及效果比较将在实验结果部分展示。

双阈值检测部分的实现方式如下：

for i in range(1, width-1):
    for j in range(1, height-1):
        if nms_img[i,j] < TL:
            dt_img[i,j] = 0
        elif nms_img[i,j] > TH:
            dt_img[i,j] = 255
        elif nms_img[i+1,j] > TH or nms_img[i-1,j] > TH or nms_img[i,j+1] > TH or nms_img[i,j-1] > TH or nms_img[i-1,j-1] > TH or nms_img[i-1,j+1] > TH or nms_img[i+1,j-1] > TH or nms_img[i+1,j+1] > TH:
            dt_img[i,j] = 255

提取彩色边缘

实现方法如下：

img_edge = np.zeros(img.shape, dtype="uint8")
img_edge = cv2.cvtColor(img_edge, cv2.COLOR_GRAY2BGR)
for i in range(width):
    for j in range(height):
        if dt_img[i,j] == 255:
            img_edge[i,j] = img_color[i,j]