基于Python的Canny边缘检测算法实战

这篇文章将解释有关Canny边缘检测，以及不使用预先编写的库编写该算法，以便我们了解Canny边缘检测的原理。

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

作为人类（我假设你是人类），我们的大脑在任何图像中都可以轻松检测到边缘，但是为了在计算机上自动执行此任务，我们必须使用可以执行该任务的程序。

以下是必须在给定数据中检测边缘的一些实际应用示例：

• 医学成像
• 指纹识别
• 在自动驾驶汽车中
• 卫星成像
• 等等……

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

不，还有相当多的方法，例如：

• Sobel边缘检测
• Prewitt边缘检测
• Laplacian边缘检测
• ……

尽管存在不同的方法，但在图像中检测边缘时，Canny边缘检测是一种广泛使用的技术。该算法由John F. Canny于1986年开发，并且自那时以来已成为图像处理中的标准技术。正如我们之前所说，Canny是一个多阶段算法，这意味着它是由许多其他算法组成的算法，我们讨论过的Sobel边缘检测就是Canny中的一个这样的算法。

Canny边缘检测的步骤：

1. 灰度转换
2. 降噪
3. 梯度计算
4. 非最大抑制
5. 双阈值和滞后阈值

1.灰度转换

要将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\'')

2. 降噪

模糊有助于平滑图像并减小像素强度中的小随机变化的影响。有许多用于模糊的卷积核：

• 高斯滤波器
• 方框滤波器
• 均值滤波器
• 中值滤波器

这里我们将使用高斯滤波器进行模糊处理，该过程是对图像使用高斯核矩阵执行卷积操作，从而得到与给定图像相对应的模糊图像。

3. 梯度计算

因此，灰度转换和模糊处理是预处理阶段。在这一步中，我们将使用Sobel滤波器，因此此步骤实际上是Sobel边缘检测方法。该过程使用Sobel滤波器，它是导数的离散近似，我们将对上述生成的图像（模糊）执行卷积操作，因此，我们将得到两个X和Y梯度，这些梯度是指向图像强度变化最大的方向的向量，使用这两个梯度，我们生成一个单一的总梯度，如果绘制该梯度，我们将得到一个由图像边缘组成的图像，我们还将将X和Y梯度之间的角度（theta）存储在一个变量中以供进一步使用。

在深入研究之前，我想问一下边缘是什么，我们如何称呼图像的特定部分是边缘？
我们可以说边缘是颜色突变的地方。

如果我们对图像使用Sobel Filter X执行卷积操作，我们将获得另一个矩阵，该矩阵显示x方向上的颜色变化。此外，在图像上卷积Sobel Filter Y会导致另一个矩阵，该矩阵显示y方向上的颜色变化。

但是，这究竟是如何发生的？
非常简单，只需查看下面的插图

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

为了找到总变化，我们使用毕达哥拉斯定理。以X梯度为底，Y梯度为高，这样斜边就是总变化。

角度可以按照下面的方式计算：

可以使用Python中的NumPy按照以下方式计算：

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，结果将如下：

而不是：

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

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

因此，在atan2中，我们不是给出y和x之间的比率，而是将它们都给出，以便算法可以进行调整。在Python中，使用NumPy，您可以如下使用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)

现在我必须问一个问题，我们结束了吗？

实际上我们不能这么说，因为根据您打算对结果做什么，它可能会有所不同。我们目前面临的问题有：

1.Sobel边缘检测算法还找到了边缘的厚度。

2.它不是二进制图像，而是灰度图像。

3.它也有许多噪音，我们可以减少。

为了解决这些问题，我们将使用非极大值算法来抑制厚度，并使用阈值进行进一步改进。

4. 非极大值抑制

在这个阶段，我们将利用theta。该过程是，我们遍历所有像素，并取当前像素的两个相邻像素并将其与它们比较，以找出当前像素的强度是否大于两个相邻像素，如果是，则继续，否则将当前像素强度设置为0。

在这里，我们必须专注于采用相邻像素，我们不能只是在当前像素周围采用一些随机像素，而必须根据角度theta进行。这个想法是取与主像素的角度几乎垂直方向的像素。例如，如果theta在22.5°和67.5°之间，则采用如下所示的像素，其中（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

仍然只是减小了边缘的宽度，我们必须使边框颜色一致并消除一些噪声。

5. 双阈值和滞后阈值

双阈值用于识别图像中的强边缘和弱边缘。梯度幅度高于高阈值的像素被视为强边缘，因此我们为其分配像素值255。梯度幅度低于低阈值的像素被视为非边缘，因此它们被分配像素值0。

梯度幅度在低和高阈值之间的像素被视为弱边缘。仅当它们连接到强边缘时，这些像素才被分配像素值255，否则将其分配为0，这个过程被称为滞后阈值。以下是计算双阈值处理过程的代码：

def double_threshold(image, low_threshold_ratio=0.05, high_threshold_ratio=0.09):
    high_threshold = image.max() * high_threshold_ratio
    low_threshold = high_threshold * low_threshold_ratio


    M, N = image.shape
    result = np.zeros((M, N), dtype=np.int32)


    strong_i, strong_j = np.where(image >= high_threshold)
    zeros_i, zeros_j = np.where(image < low_threshold)


    weak_i, weak_j = np.where((image <= high_threshold) & (image >= low_threshold))


    result[strong_i, strong_j] = 255
    result[zeros_i, zeros_j] = 0
    result[weak_i, weak_j] = 50  # 用于滞后阈值


    return result

以下是滞后阈值处理过程的代码：

def hysteresis(image, weak=50):
    M, N = image.shape
    result = np.zeros((M, N), dtype=np.int32)


    strong_i, strong_j = np.where(image == 255)
    weak_i, weak_j = np.where(image == weak)


    for i, j in zip(strong_i, strong_j):
        result[i, j] = 255


    for i, j in zip(weak_i, weak_j):
        if any(image[x, y] == 255 for x in range(i-1, i+2) for y in range(j-1, j+2)):
            result[i, j] = 255
        else:
            result[i, j] = 0


    return result

结论

希望现在您对Canny边缘检测算法有了清晰的理解。从头开始编写代码并不是一种不好的做法，它实际上比仅仅解释更好地帮助您理解事物。但是，当您在实际项目中工作时，您不必从头开始编写代码，然后您可以使用诸如OpenCV之类的库来帮助您。在Python中使用OpenCV，您可以生成如下的高斯模糊图像

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)