使用 Python 从头开始编写 Canny 边缘检测算法。

原始图像（左）和检测到的边缘（右）| 图片由作者提供

一、说明

在本文中，我将解释有关 Canny 边缘检测的所有内容，以及在不使用一些预先编写的库的情况下对算法进行编码，以便您能够了解真正发生的情况。

二、关于Canny算子和边缘检测

Canny 是一种多级边缘检测算法，可以检测图像中的各种边缘。

但是，等等......为什么我们需要检测图像中的边缘？

作为人类（我假设你是人类），我们的大脑可以毫无问题地检测任何图像中的边缘，但为了在计算机上自动执行任务，我们必须使用可以完成该任务的程序。

一些必须检测给定数据中的边缘的实际应用程序的示例，

医学影像
指纹识别
在自动驾驶汽车中
卫星成像
ETC，…

在检测边缘时，canny 是唯一的选择吗？

不是的，方法有很多，分别是

索贝尔边缘检测
普威特边缘检测
拉普拉斯边缘检测
ETC，…

尽管有多种不同的方法，但Canny 边缘检测是一种广泛使用的图像边缘检测技术。

该算法由John F. Canny 于 1986 年开发，现已成为图像处理的标准技术。

正如我们之前所说的 Canny 是一种多阶段算法，这意味着它是一种包含许多其他算法的算法，我们讨论的 Sobel 边缘检测就是 Canny 中的一种这样的算法。

三、Canny边缘检测的步骤

灰度转换
降噪
梯度计算
非极大值抑制
双阈值和迟滞

3.1. 灰度转换

RGB 彩色图像（左）和灰度图像（右）| 图片由作者提供

要将 HSV、YUV 或 RGB 色阶转换为灰度，我们会这样：

“呵呵，这么简单，取每个通道的平均值就可以了。”

理论上，该公式是100%正确的，但平均法并不能按预期工作。

原因是人脑对 RGB 的反应不同。眼睛对绿光最敏感，对红光不太敏感，对蓝光最不敏感。因此，当我们生成灰度图像时，三种颜色应该具有不同的权重。

这给我们带来了另一种方法，称为加权方法，也称为光度方法。

实现起来相当简单，

权重是根据它们的波长计算的，因此改进后的公式如下所示，

Grayscale = 0.2989 * R + 0.5870 * G + 0.1140 * B

但我认为在我们的上下文（边缘检测）中，我们也可以使用平均方法，尽管我遵循惯例。

import numpy as np


def to_gray(img: np.ndarray, format: str):
    '''
    Algorithm:
    >>> 0.2989 * R + 0.5870 * G + 0.1140 * B 
    - Returns a gray image
    '''

    r_coef = 0.2989
    g_coef = 0.5870
    b_coef = 0.1140

    if format.lower() == 'bgr':
        b, g, r = img[..., 0], img[..., 1], img[..., 2]
        return r_coef * r + g_coef * g + b_coef * b
    elif format.lower() == 'rgb':
        r, g, b = img[..., 0], img[..., 1], img[..., 2]
        return r_coef * r + g_coef * g + b_coef * b
    else:
        raise Exception('Unsupported value in parameter \'format\'')

3.2 降噪

图片进行模糊化，为什么要这样做？最根本目的是，营造一个图片信息连续化处理，即将台阶函数变成逐步趋近，避免对微积分运算（尤其高阶运算）造成中断问题。此处理可以极大防止不可微点的存在。

灰度图像（左）和高斯模糊图像（右）| 图片由作者提供

模糊有助于平滑图像并减少像素强度中小的随机变化的影响。

有很多用于模糊的内核，

高斯滤波器
箱式过滤器
均值滤波器
中值滤波器

这里我们将使用高斯滤波器进行模糊，其实我已经写了一篇完整的文章了，所以我不会详细解释它。

编码：在 Python 中从头开始进行高斯模糊操作。

编码底层概念而不调用一些自定义函数，就这样，您将接触到什么......

该过程是将图像与高斯核矩阵进行卷积运算，得到对应给定图像的模糊图像。

3.3.梯度计算

因此，灰度转换和模糊是预处理阶段。

在这一步中，我们将使用 Sobel 滤波器，因此这一步实际上是 Sobel 边缘检测方法。

在解释细节之前，我将一次性告诉您完整的故事。

该过程使用Sobel 滤波器，它是导数的离散近似，我们将对上面生成的图像（模糊）执行卷积运算，

因此，我们将得到两个X 和 Y 梯度，梯度是指向图像强度变化率最大方向的向量，

使用这两个梯度，我们生成一个总梯度，如果我们绘制该梯度，我们会得到一个由图像边缘组成的图像，

此外，我们还将把X 和 Y 梯度之间的角度 (theta)存储在变量中以供进一步使用。

Sobel X 和 Y 滤波器 | 图片由作者提供

在开始讨论之前，我想问你什么是边缘，我们如何将图像的特定部分称为边缘？

我们可以说边缘是颜色突然变化的部分。

因此，在阅读其余部分时请记住这一点。

如果我们使用 Sobel 滤波器 X 对图像进行卷积运算，我们将得到另一个矩阵，该矩阵显示 x 方向上的颜色变化。

此外，在图像上对 Sobel 滤波器 Y 进行卷积会产生另一个矩阵，该矩阵显示 y 方向上的颜色变化。

但是，这到底是怎么发生的呢？
其实很简单，看下面的图就知道了。

索贝尔滤波器如何生成 X 和 Y 梯度的图示 | 图片由作者提供

现在，如果我们在真实图像上执行此操作，我们会得到如下所示的结果，

Y&X 的变化 | 图片由作者提供

为了找到总变化，我们利用毕达哥拉斯定理。

以Base为 X 梯度，Altitude为 Y 梯度，因此斜边将导致总变化。

斜边的 eqn | 图片由作者提供

θ 的计算公式如下：

eqn 计算 theta | 图片由作者提供

可以使用 numpy 在 python 中计算，如下所示。

theta = np.arctan(Gradient_Y / (Gradient_X + np.finfo(float).eps))
# np.finfo(float).eps is added to tackle the division by zero error

但是，如果您看下图，您会发现，即使我们以这种方式获得了正确的角度，方向性也可能会丢失。

反正切方向性问题的示例 | 图片由作者提供

所以你找到的角度不会在整个圆的范围内。

这样，如果我们计算 theta 就会像下面这样：

使用 atan 检测到的角度矩阵的图片表示 | 图片由作者提供

代替，

使用 atan2 检测到的角度矩阵的图片表示 | 图片由作者提供

为了解决这个问题，数学家在计算时添加了一些条件，因此新方法称为atan2。

根据维基百科，条件如下：

从atan 转换为atan2 的条件| 图片由作者提供

因此，在 atan2 中，我们没有给出 y 和 x 之间的比率，而是给出两者，以便算法可以对其进行调整。

在使用 numpy 的 python 中，您可以按如下方式使用 atan2，

theta = np.arctan2(Gradient_Y, Gradient_X)

最后，您了解了所有细节，并找到了边缘。

这是计算梯度和 theta 的代码（theta 将在接下来的阶段中使用）。

G = np.sqrt((Gradient_X ** 2.0)+(Gradient_Y ** 2.0))
''' or simply do the following '''
G = np.hypot(Gradient_X, Gradient_Y)

G = G / G.max() * 255 # the total gradient; ie, the edge detection result
theta = np.arctan2(Gradient_Y, Gradient_X)

sobel边缘检测结果| 图片由作者提供

现在我要问你一个问题，我们结束了吗？

实际上我们不能这么说，因为根据您计划如何处理结果，结果可能会有所不同。

我们目前遇到的问题是，

Sobel 边缘检测算法还可以找到边缘的厚度。
它不是二值图像，而是灰度图像。
而且它还有很多我们可以减少的噪音。

为了解决这些问题，我们将使用非极大值算法抑制厚度，并进行阈值处理以使其更好

3.4.非极大值抑制

非极大值抑制结果| 图片由作者提供

在这个阶段，我们将利用theta，

其过程是，循环所有像素点，取当前像素点的两个相邻像素点进行比较，判断当前像素点的强度是否大于这两个相邻像素点的强度，如果是则继续，如果不是则设置当前像素强度为 0。

这里我们要把重点放在取相邻像素上，我们不能只取当前像素周围的一些随机像素，而是必须根据角度（theta），

这个想法是取几乎垂直于主像素角度方向的像素。

例如，如果 theta 介于 22.5° 和 67.5° 之间，则采用如下所示的像素，其中(i, j) 是当前像素。

用作 (i,j) 像素的相邻像素的像素的图示 | 图片由作者提供

下面是该过程的另一个说明，箭头表示为角度，我们的目标是仅检测灰色阴影像素。

| 作者图片 |

下面是非极大值抑制的代码，

'''Non Max Suppression'''

M, N = G.shape
Z = np.zeros((M,N), dtype=np.int32) # resultant image
angle = theta * 180. / np.pi        # max -> 180, min -> -180
angle[angle < 0] += 180             # max -> 180, min -> 0

for i in range(1,M-1):
    for j in range(1,N-1):
        q = 255
        r = 255
        
        if (0 <= angle[i,j] < 22.5) or (157.5 <= angle[i,j] <= 180):
            r = G[i, j-1]
            q = G[i, j+1]

        elif (22.5 <= angle[i,j] < 67.5):
            r = G[i-1, j+1]
            q = G[i+1, j-1]

        elif (67.5 <= angle[i,j] < 112.5):
            r = G[i-1, j]
            q = G[i+1, j]

        elif (112.5 <= angle[i,j] < 157.5):
            r = G[i+1, j+1]
            q = G[i-1, j-1]

        if (G[i,j] >= q) and (G[i,j] >= r):
            Z[i,j] = G[i,j]
        else:
            Z[i,j] = 0

尽管如此，我们只是减少了边缘的宽度，我们必须使边框颜色一致并消除一些噪音。

3.5. 双阈值和迟滞

双阈值用于识别图像中的强边缘和弱边缘。

梯度幅度高于高阈值的像素被视为强边缘，因此我们为其分配像素值 255。梯度幅度低于低阈值的像素被视为非边缘，因此它们的像素值为 0。

梯度幅度在低阈值和高阈值之间的像素被视为弱边缘。这些像素只分配像素值255，如果它们连接到强边缘，否则分配像素值0给它，这个过程称为滞后。

下面是计算双阈值处理的代码。

def threshold(img, lowThresholdRatio=0.05, highThresholdRatio=0.09):
    '''
    Double threshold
    '''
    
    highThreshold = img.max() * highThresholdRatio;
    lowThreshold = highThreshold * lowThresholdRatio;
    
    M, N = img.shape
    res = np.zeros((M,N), dtype=np.int32)
    
    weak = np.int32(25)
    strong = np.int32(255)
    
    strong_i, strong_j = np.where(img >= highThreshold)
    zeros_i, zeros_j = np.where(img < lowThreshold)
    
    weak_i, weak_j = np.where((img <= highThreshold) & (img >= lowThreshold))
    
    res[strong_i, strong_j] = strong
    res[weak_i, weak_j] = weak
    
    return (res, weak, strong)

双阈值结果 | 图片由作者提供

下面是滞后过程的代码。

def hysteresis(img, weak, strong=255):
    M, N = img.shape

    for i in range(1, M-1):
        for j in range(1, N-1):
            if (img[i, j] == weak):
                if (
                    (img[i+1, j-1] == strong) or (img[i+1, j] == strong) or
                    (img[i+1, j+1] == strong) or (img[i, j-1] == strong) or
                    (img[i, j+1] == strong) or (img[i-1, j-1] == strong) or
                    (img[i-1, j] == strong) or (img[i-1, j+1] == strong)
                ):
                    img[i, j] = strong
                else:
                    img[i, j] = 0
    return img

滞后结果（最终图像）| 图片由作者提供

完整的代码可以在我的 GitHub 存储库中找到，

GitHub — rohit-krish/CVFS：从“从头开始”编码计算机视觉相关算法。

您目前无法执行该操作。您使用另一个选项卡或窗口登录。您在另一个选项卡中退出或...

四、结论

希望现在您对 Canny 边缘检测算法有了清晰的了解。

从头开始编码并不是一个坏习惯，

它实际上可以帮助你更好地理解事物，而不仅仅是解释。

但是，当您处理实际项目时，您不必从头开始编码，那么您可以使用 OpenCV 等库来提供帮助。

在使用 OpenCV 的 Python 中，您可以生成如下高斯模糊图像，

import cv2

img = cv2.imread(<img_path>)

# to grayscale
imgGray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# guassian blur
imgBlur = cv2.GaussianBlur(imgGray, (13, 13), 0)

# canny edge detection
imgCanny = cv2.Canny(imgBlur, 50, 50)

五、参考

Canny Edge Detector
https://towardsdatascience.com/canny-edge-detection-step-by-step-in-python-computer-vision-b49c3a2d8123
Image Processing 101 Chapter 1.3: Color Space Conversion | Dynamsoft Blog