Canny边缘检测原理和从零实现逐行代码详解

Canny边缘检测原理和逐行代码详解

理论基础

我们希望最后的图像是一个二值图，就是说只有两种情况，要么是边缘要么不是边缘。并且最后的边缘是一个很细的线，甚至只有一个像素点。

卷积算子

我们在进行边缘检测前需要执行一个去噪的步骤，去噪我们一般使用高斯核。

梯度大小和方向

观察图片我们可以发现，边一般的特征就是在边的两边（垂直方向）亮度变化比较大，这就类似与连续曲线的求导，最后我们再选取导数比较大的像素（即亮度变化比较大）即可。但是图像是一个一个像素点，是离散的数据，所以我们没法求导数，我们回顾一下导数的定义：
$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x-h)}{2h}$$
上面是中心差分，就是连续函数求导的意思。由于图像是离散数据，我们可以想象一下，如果我们令$h=1$，就是一个像素点，那么这个时候在图像中就是类似于$h\to 0$，因为一个像素点的差距就是最小的了，那么我们让某个像素点旁边两个像素相减就类似于对这个像素点求导了。由于一个像素点有两个方向，即水平方向和竖直方向，那我们就有一下两个求导算子，分别对水平方向和数值方向求导。

根据上面两个sobel算子$Gx$, $Gy$，算出所有像素点两个方向的导数后，然后再根据公式求出所有像素点的梯度，公式为：
$$gradient=\sqrt{{Gx}^2+{Gy}^2}$$
求梯度方向公式为：
$$angle(\theta )=ta{n}^{-1}(\frac{Gy}{Gx})$$
梯度大小很好理解，就是求出一个像素点的变化程度大小，因为对于变化程度，我们不仅要考虑水平方向，还要考虑垂直方向，所以我们要把它俩加起来，因为导数可能有正有负，最后我们想要一个正数，所以先平方再加再开根号。

为什么还要求方向呢？上面我们说了我们希望边缘就是确定边缘的位置在哪，我们把它标记出来就行，所以对于最后的图像的边缘点要是非常细的，只有一个像素点。为了完成这个目标，我们通过找到某个像素点变化最快的方向，在这个方向如果没有比这个点更大的梯度像素点，那么这个像素点就是边缘。如果有比这个点梯度更大的点，那么我们就认为这个点不是边缘，把它的大小置为0。

非极大值抑制

由于图像数据是离散的，我们求出梯度某个像素点的梯度方向后，可能出现改梯度方向没有像素点的情况，那这样就没有像素点和该像素点比较了，我们不希望这样的情况发生。

如上图，点C的梯度方向落在了g1和g2两点之间，g3和g4两点之间，那么我们就求g1和g2的加权和，权重是多少呢？离那个点(dTmp1,dTmp2)越近权重就越大，而且重要程度随着距离线性变化，知道了$\theta$，就很容易确定权重了。所以我们有公式：
$$dTmp1=tan\theta g_2+(1-tan\theta )g_1$$

边缘的连接

首先我们使用双阈值进行处理，设置一个最大值，一个最小值。大于最大值的像素点，我们认为该像素点一定是边缘。小于最小值的像素点我们认为它一定不是边缘，我们把它置为0. 一系列操作后，边缘点并不是连续的，可能断断续续的，造成不太美观，所以我们要把连续的点连再一起。怎么连呢？由于我们已经确定了一定为边缘的点。那么我们就从这些点开始进行连接，我们利用栈的做法，把这些点压栈，然后判断栈首的八领域内有没有可能为边缘的点，如果有就把这两个点连起来（赋值为1），然后将栈首元素弹出，直至栈空为止。

代码逐行详解

首先使用python实现我们的思路，后续会考虑使用c++实现。

普通卷积

第一步实现卷积，先使用三重循环实现卷积，下面再使用im2col的方法实现更加高效的卷积。

import numpy as np
import matplotlib.pyplot as plt

def convolve(image, kernel, padding, stride):
    """
    :param image: 要卷积的图像
    :param kernel: 卷积核
    :param padding: 输入要上下左右四个方向要padding的宽度，例如[1,1,1,1]代表四个方向padding宽度为1
    :param stride: 代表卷积的步长
    :return: 卷积后的图像
    """
    result = None  # 保存最后的结果，先在函数域声明一下
    img_size = image.shape  # 输入图像的形状
    filter_size = kernel.shape  # 滤波核的形状，传统方法不同于深度学习，滤波核是二维的，深度学习卷积核一般是三维的

    if len(img_size) == 3:  # 判断输入图像维度，如果是三维的，一般指彩色图像
        result = []  # 用来存储卷积之后的结果
        for i in range(0, img_size[-1]):  # 遍历通道
            channel = []  # 创建一个通道列表，用来存储通道数据
            padded_img = np.pad(image[:, :, i], ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')  # 对每个通道进行padding
            # 遍历行维度，添加一个列表到channel尾部，用来存放生成数据
            for j in range(0, img_size[0], stride):  # 遍历行
                channel.append([])  # 添加一个空列表到通道列表中，空列表用来存储卷积后的每一行的数据
                for k in range(0, img_size[1], stride):  # 遍历列
                    val = (kernel * padded_img[j:j + filter_size[0], k:k + filter_size[1]]).sum()  # 将滤波核与遍历图像相乘
                    channel[-1].append(val)  # 将得到的数据添加到刚刚的空列表中
            result.append(channel)  # 最后将通道列表添加到最后的结果中

    elif len(img_size) == 2:  # 如果输入图像是二维的，一般指灰度图或二值图
        result = []  # 存储卷积之后的结果
        padded_img = np.pad(image, ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')  # 进行padding
        for j in range(0, img_size[0], stride):  # 遍历行
            result.append([])  # 因为图像是二维的，就不像三维的需要一个通道维度，我们直接添加一个空列表，作为输出结果的行
            for k in range(0, img_size[-1], stride):  # 遍历列
                val = (kernel * padded_img[j:j + filter_size[0], k:k + filter_size[1]]).sum()  # 计算卷积结果
                result[-1].append(val)  # 将结果添加到行中

    return result

im2col方式卷积

下面介绍一下im2col的方式：

上面的输入图像是三维的，但是我们用的是二维数据，不过不影响我们进行理解。将图像按照一块一块的方式变成一行，然后放在一块，滤波器变成列，我们只有一个滤波器就只有一列，然后按照矩阵相乘的方法进行运算，由于矩阵相乘速度非常快，所以这种方法卷积速度就很快，下面会进行实验验证准确性以及相比普通卷积速度快了多少。先放出im2col的代码。

def im2col(image, filter_h, filter_w, padding, stride):
    """
    只支持二维图像，传统方法一般都是用灰度图或二值图
    :param image: 要转换的图像
    :param filter_h: 滤波的长
    :param filter_w: 滤波的宽
    :param padding: padding的宽度
    :param stride: 卷积的步长
    :return: 要转换的图像
    """

    matrix = []  # 保存最后的输出结果
    img_shape = image.shape  # 读取输入图像形状
    image = np.pad(image, ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')  # 对输入图像进行pad
    for j in range(0, img_shape[0], stride):  # 遍历行
        for k in range(0, img_shape[1], stride):  # 遍历列
            col = image[j:j + filter_w, k:k + filter_h].reshape(1, -1)  # 找到对应位置的数据（二维的），将其reshape成一行数据
            matrix.append(col)  # 将reshape的数据添加到结果列表中

    return np.array(matrix)

上面只是将图像im2col的代码，想完成最后的卷积运算，我们可以这样：

sobel_x = sobel_kernel_x.reshape(-1, 1)
outs = np.dot(col_img*255, sobel_x).reshape(img_shape[0], img_shape[1])

假设sobel_kernel_x是一个滤波核，我们先将它变成一列，然后和刚刚已经im2col的图像进行矩阵相乘，最后在reshape回原图像的尺寸进行了.

下面就验证im2col卷积的正确性以及速度了。

从上我们可以看出，使用普通卷积和使用im2col然后矩阵相乘的方法完成的卷积最后得到的结果是一模一样的。运算时间一副512*512像素的图像，普通卷积花费时间为：0.8410706520080566。im2col然后矩阵相乘运算时间为：0.37093663215637207。可以看到提升效果还是很明显的。另外我们关注一下sobel算子x方向和y方向图像的却别，x方向图像呈竖的条纹状，y方向是成横的条纹状的。

非极大值抑制（Non Maximum suppress）

def non_maximum_suppress(grad_l, grad_d):
    """
    :param grad_l: 图像的梯度大小
    :param grad_d: 图像的梯度方向
    :return: 抑制之后图像的梯度大小
    """
    h = grad_l.shape[1]  # 获取图像的行数
    w = grad_l.shape[0]  # 获取图像的列数
    result = np.zeros((h, w))  # 创建和图像同样形状的全零矩阵，用来存放结果数据

    for i in range(1, h - 2):  # 遍历行
        for j in range(1, w - 2):  # 遍历列

            eight_neibor = grad_l[i - 1:i + 2, j - 1:j + 2]  # 取出图像梯度大小的八领域
            if 0 <= grad_d[i][j] <= 1:  # 如果梯度大小为[0,1]之间，即角度在0-45和180-225
                dTmp1 = (1 - grad_d[i][j]) * eight_neibor[1][2] + grad_d[i][j] * eight_neibor[0][2]  # 按照公式计算dTmp1
                dTmp2 = (1 - grad_d[i][j]) * eight_neibor[2][0] + grad_d[i][j] * eight_neibor[1][0]  # 按照公式计算dTmp2
                if grad_l[i][j] > dTmp1 and grad_l[i][j] > dTmp2:  # 如果这个中心点比dTmp1和dTmp2都大，那么这个点就是极大值点，
                    # 我们就把它放到结果矩阵对应的【i，j】位置
                    result[i][j] = grad_l[i][j]  # 放到结果矩阵中
            if grad_d[i][j] > 1:  # 角度在45-90和225-270
                dTmp1 = (1 - 1 / grad_d[i][j]) * eight_neibor[0][1] + 1 / grad_d[i][j] * eight_neibor[0][2]
                dTmp2 = (1 - 1 / grad_d[i][j]) * eight_neibor[2][1] + 1 / grad_d[i][j] * eight_neibor[2][0]
                if grad_l[i][j] > dTmp1 and grad_l[i][j] > dTmp2:
                    result[i][j] = grad_l[i][j]
            if grad_d[i][j] < -1:  # 角度在90-135和270-315
                dTmp1 = (1 - 1 / grad_d[i][j]) * eight_neibor[0][1] + 1 / grad_d[i][j] * eight_neibor[0][0]
                dTmp2 = (1 - 1 / grad_d[i][j]) * eight_neibor[2][1] + 1 / grad_d[i][j] * eight_neibor[2][2]
                if grad_l[i][j] > dTmp1 and grad_l[i][j] > dTmp2:
                    result[i][j] = grad_l[i][j]
            if -1 <= grad_d[i][j] <= 0:  # 角度在135-180和315-360
                dTmp1 = (1 - grad_d[i][j]) * eight_neibor[1][0] + grad_d[i][j] * eight_neibor[0][0]
                dTmp2 = (1 - grad_d[i][j]) * eight_neibor[1][2] + grad_d[i][j] * eight_neibor[2][2]
                if grad_l[i][j] > dTmp1 and grad_l[i][j] > dTmp2:
                    result[i][j] = grad_l[i][j]

    return result  # 返回最后的结果

双阈值处理

def twothrehold(inp, low, high):
    """
    :param inp: 输入要双阈值抑制的图像
    :param high: 高阈值，
    :param low: 低阈值，
    :return:
    """
    inp[inp >= high] = 255  # 比这个阈值大的一定是边缘点，都记为255
    inp[inp <= low] = 0  # 比这个阈值低的一定不是边缘点，都记为0

    return input

边缘连接

def connect(inp, low, high):
    """
    对处于阈值之间的像素点进行进一步判断，如果该点旁边有一个肯定是边缘点的，那么我们就将这个点视为边缘点，这样就把断断续续的边缘点连接起来
    了。具体怎么做呢？利用栈，把通过双阈值的已经确定为边缘点的点全部压入栈中，然后弹栈，将弹出的元素的八邻域中在双阈值之间的点视为边缘点，
    再将这些点压入栈中。再弹栈，重复操作，直至栈空。
    :param inp:  图像梯度大小
    :param low:  高阈值
    :param high: 低阈值
    :return:
    """
    st = []  # stack 栈
    img_shape = inp.shape
    cood = []  # 存放像素点坐标

    for i in range(1, img_shape[0]):  # 遍历行
        cood.append([])  # 图像是二维的，所以该空列表用来存放行的坐标
        for j in range(1, img_shape[1]):  # 遍历列
            cood[-1].append([i, j])  # 将坐标添加到cood中
            if inp[i][j] == 255:  # 如果这个梯度值为255，那么它一定是边缘点，那么将其入栈
                st.append([i, j])  # 入栈
    cood = np.array(cood)

    while len(st) != 0:  # 判断栈是否为空
        d = st.pop()  # 弹出栈尾元素
        cood_eight_neibor = cood[d[0] - 1:d[0] + 2, d[1] - 1:d[1] + 2, :]  # 找到这个元素的八领域坐标

        for i in range(3):
            for j in range(3):
                if low < inp[cood_eight_neibor[i][j][0], cood_eight_neibor[i][j][1]] < high:  # 如果八领域内梯度值在阈值之间
                    inp[cood_eight_neibor[i][j][0], cood_eight_neibor[i][j][1]] = 255  # 一定是边缘点，记为255
                    st.append(list(cood_eight_neibor[i][j]))  # 再将其入栈

    inp[inp != 255] = 0  # 剩下的点都不是边缘点，记为0

    return inp

在边缘检测之前我们一般会使用高斯模糊核对图像进行模糊，消除噪声影响，高斯核生成函数为：

def gaussian_kernel(size, sigma):
    kernel = np.zeros((size, size))
    mid = size // 2
    for i in range(size):
        for j in range(size):
            kernel[i][j] = (1 / (2 * np.pi * sigma ** 2)) * np.exp(
                -((i - mid) ** 2 + (j - mid) ** 2) / (2 * sigma ** 2))

    return kernel / np.sum(kernel)

最后我们得到最终的结果：

opencv官方结果这个样子：

可以看到我们自己的代码也取得了较好的结果，但是和opencv相比还是有差距的，代码还是有优化空间的，最值得优化的地方应该是边缘连接这个步骤，应该有更好的算法实现边缘链接。