Numba 下自适应双阈值的Canny 边缘检测并行算法

摘要：针对自适应双阈值的Canny边缘检测算法中运算量大、耗时长的特点，提出了一种基于CUDA架构的GPU环境下的自适应双阈值的Canny边缘检测并行算法。针对算法过程中的寻找梯度幅值极大值操作，设计了一种基于CUDA的并行化梯度幅值极大值统计方案，并使用Anaconda推出的Python编译器Numba实现基于CUDA架构的Canny并行算法。在不同尺寸图像上进行串行算法与并行算法的对比实验，实验结果表明，并行算法在能够满足检测精度的同时可获得较高加速比，最高加速比可达95.27，加速效果显著，具有良好的数据可扩展性。

关键词： 边缘检测；图形处理器；统一计算设备架构；并行算法；自适应双阈值

中图分类号：TP391 文献标识码：A

文章编号：1009-3044（2024）31-0034-06

开放科学（资源服务）标识码（OSID） ：

0 引言

边缘检测被广泛应用于目标物体识别、目标追踪、农作物检测、医学成像等领域，在计算机视觉领域中具有非常重要的地位。在边缘检测算法中，由JohnF. Canny于1986年提出的Canny边缘检测算法[1]，由于其良好的检测效果，一直作为一种标准的边缘检测算法，此后也出现了多种基于Canny算法的改进算法。Canny算法及其改进算法近年来被广泛应用于图像的边缘检测中。然而，Canny边缘检测算法存在计算量大、耗时较多及计算效率低的问题。

针对Canny边缘检测中的效率问题，唐斌等人[2]提出了一种基于GPU+CPU的快速实现Canny算子的方法，获得的最高加速比为5.39；郭恒亮等人[3]设计了一种面向FT-M7002 平台的Canny 梯度计算并行算法，将运行速度提升了1.490至2.112倍；刘振涛等人[4]设计了一种将Hyper-Q技术应用于双站位相机图像Canny边缘提取算法的方案，实现了11.6倍的速度提升；吴进等人[5]基于SDSoC开发环境，选用ZC706作为开发平台对Canny边缘检测进行加速，显著降低了检测用时；蓝贵文等人[6]设计了一种基于CUDA的改进Canny边缘检测算法，获得的最高加速比为24；张帆等人[7]利用GPU流处理器数量众多的优势以及强大的多线程并发执行能力对Canny算法进行并行加速，使执行速度从80毫秒降至6毫秒以内；查志勋等人[8]利用FPGA在并行和高速处理数据方面的特点，采用结合灰度变换的改进型Canny算法，设计出了一种实时边缘检测系统。

非自适应双阈值的Canny边缘检测算法在双阈值过滤操作中需手动指定高低阈值，因此一般情况下该边缘检测算法不具备很强的自适应性。相反，在自适应双阈值Canny算法中，高低阈值是根据梯度幅值大小自动计算的，这样不仅能检测到更多边缘细节，还具有较强的自适应性。然而，在高低阈值的计算过程中需要涉及寻找梯度幅值极大值的操作，然而根据GPU架构的特点，GPU并不擅长执行极值等统计操作。针对这一问题，一种可行的解决方案是在程序运行到计算高低阈值时将数据传回主机，由主机完成计算后再将结果传回GPU端。这样会增加一次数据在主机端与设备端之间的传输，从而延长整个检测算法的耗时。因此，本文提出了一种基于CUDA 架构的GPU环境中并行化二维规约梯度幅值极大值统计方案，使整个操作在GPU中完成，避免了GPU与主机间的数据来回传输，从而提高了算法的整体运行效率。

Python语言是一种广泛应用的脚本语言，涵盖数据分析、自然语言处理、机器学习、科学计算以及推荐系统构建等多个领域[9]。近年来，随着人工智能与大数据的快速发展，Python的应用逐渐增加，已成为最受欢迎的编程语言之一。本文使用Python语言，首先在CPU主机端实现Canny串行算法，然后对算法进行并行化设计，利用由Anaconda 推出的Python编译器Numba实现基于CUDA架构的Canny并行算法，并将所提出的并行算法封装为函数，提供调用接口，以便用户能够直接使用Python语言进行调用，从而提升边缘检测效率。

1 Canny 算法

Canny边缘检测的基本原理有三条[1]：

1） 良好的信噪比，即将非边缘点判定为边缘点的概率应低，同时将边缘点判定为非边缘点的概率也应低。

2） 高的定位性能，即检测出的边缘点应尽可能位于实际边缘的中心。

3） 对单一边缘仅产生唯一响应，即单个边缘产生多个响应的概率应低，且虚假响应的边缘应得到最大程度的抑制。

Canny边缘检测算法主要由四个步骤组成，依次为：高斯滤波、梯度幅值和方向计算、非极大值抑制以及自适应双阈值过滤处理。这四个步骤均基于灰度图像完成，当读入的图像为彩色图像时，首先应将其进行灰度化处理。

1.1 高斯滤波

由于图像中存在的噪声点会对边缘检测产生很大影响，因此在进行边缘检测前，需要先去除灰度图中的噪声点。在Canny算法中，采用高斯滤波对图像进行平滑，以去除图像中的噪声点。首先，利用二维高斯函数生成高斯卷积核，其二维高斯函数如下所示：

式中：σ为高斯模糊因子，使用所生成的二维高斯卷积核与图像进行卷积操作，完成高斯滤波，公式如下所示：

Is （x，y ） = I （x，y ） × G （x，y ） （2）

式中： I （x，y ）为原始图像；Is （x，y ）为经过卷积后的图像。

1.2 计算梯度

在对图像完成高斯滤波后，接下来进行梯度计算。梯度计算包括梯度幅值的计算和梯度方向的计算。在图像处理中，通常使用一阶有限差分近似来求取灰度值的梯度值。具体的公式如下所示：

式中：I 为图像灰度值，Gx为X方向的偏导数，Gy为Y方向的偏导数，M 为当前点梯度的幅值，θ 是当前点梯度的方向，即角度。图像中梯度方向与边缘方向互相垂直，其关系如图1所示。

1.3 非极大值抑制

计算出图像的梯度后，还须进一步选择局部真实边缘点，非极大值抑制能够保留局部最大梯度值点而抑制所有非局部最大梯度值的点。即只保留梯度变化中最剧烈的点。方法是将梯度方向依次分为0°、45°、90°、135°４个方向。将当前像素点M（i，j）的梯度幅值与其梯度方向上两个像素点的梯度幅值进行比较，若当前像素点的梯度幅值最大，则保留该点梯度幅值，否则将该点梯度幅值置零。

1.4 自适应双阈值

经过非极大抑制后，图像中仍存在许多噪声点。在 Canny 算法中，使用了一种双阈值处理技术，即分别设置高阈值TH和低阈值TL来区分边缘像素点。如果边缘点的梯度值大于高阈值TH，则判定该点为强边缘点并予以保留。若边缘点的梯度值小于高阈值TH且大于低阈值TL，则标记为弱边缘点。小于低阈值的点则被抑制，即将该点的梯度幅值置为零。

在自适应双阈值处理中，高低阈值的选择按照3∶1的比例来选取。即首先找到图像中所有像素点的梯度幅值的最大值，然后高阈值TH和低阈值TL的计算公式如下所示：

TH = MAX （M [i，j ]） × 0.3 （7）

TL = MAX （M [i，j ]） × 0.1 （8）

对于被标记为弱边缘的像素点，需要进行进一步处理。具体方法是检查该弱边缘点的 8 邻域中的所有像素点。如果该邻域中存在被标记为强边缘的点，则该弱边缘点被判定为真边缘并予以保留；否则，该点被判定为假边缘点并被剔除，即将该点的梯度幅值置为零。

2 Canny 算法的并行化设计

2.1 算法步骤

在进行 Canny 算法的并行化之前，首先需要在主机端设计并实现串行算法。本文中串行算法的步骤如算法 1 所示。

算法1： Canny串行算法。

输入：图像数据，高斯卷积核尺寸Ksize=5，高斯模糊因子=0.1。

输出：边缘二值图像。

步骤1：将图像转换为灰度图，并进行零填充。在进行高斯滤波操作之前，为了使图像边界处的点也能够参与卷积运算，需要对图像的上下左右边界各进行pad=Ksize/2的零填充。本文中设置Ksize=5，因此pad=2，如图2所示。

步骤 2：生成大小为5×5，模糊因子为 0.1的高斯卷积核。

步骤 3：进行高斯卷积运算。

步骤 4：图像裁剪。由于在步骤 1 中对图像进行了零填充，因此在卷积完成后，需要将图像恢复为原始大小。为此，需要将上下左右边界处各裁剪掉pad 大小的像素点。可以利用NumPy中的切片功能完成图像裁剪。

步骤 5：梯度计算。本文采用边缘差分算子 So⁃bel 来计算每个像素点在水平方向的差分Gx和垂直方向的差分Gy，进而计算出梯度的幅值和方向。Sobel 卷积因子如图 3 所示。

步骤 6：非极大值抑制。在 Canny 原算法中，非极大值抑制仅在0°、45°、90°、135°四个梯度方向上进行。这样做的优点是简单，但如果边缘点的梯度方向不沿这四个方向，结果可能不够准确。然而，由于在实际存储中，像素点是离散的二维矩阵，故沿当前点梯度方向两侧的点不一定存在，这类点称作亚像素点。为了处理这类情况，可以使用线性插值法来计算亚像素点的梯度值，如图 4 所示。

在图 4 中，（1） 和 （2） 两种情况表示|Gy|>|Gx|，即该点的梯度方向更靠近 Y 轴方向，因此权重w = |Gx|/|Gy|。其中，情况 （1） 表示Gx与Gy的方向相同，其插值可以表示为：

gu[i，j]= w × g[i - 1，j - 1]+ （1 - w） × g[i - 1，j] （9）

gd[i，j]= w × g[i + 1，j + 1]+ （1 - w） × g[i + 1，j] （10）

情况（2）表示Gx与Gy方向相反，插值表示为：

gu[i，j]= w × g[i - 1，j + 1]+ （1 - w） × g[i - 1，j] （11）

gd[i，j]= w × g[i + 1，j - 1]+ （1 - w） × g[i + 1，j] （12）

在图4中，（3）和（4）两种情况表示|Gx|>|Gy|，即该点梯度方向更靠近X轴方向， 则权重w = |Gy|/|Gx|。其中，情况（3）表示Gx与Gy的方向相同，其插值可表示为：

gl[i，j]= w × g[i + 1，j - 1]+ （1 - w） × g[i + 1，j] （13）

gr [i，j]= w × g[i - 1，j + 1]+ （1 - w） × g[i，j + 1] （14）

情况（4）表示Gx与Gy方向相反，插值表示为：

gl[i，j]= w × g[i - 1，j - 1]+ （1 - w） × g[i，j - 1] （15）

gr [i，j]= w × g[i + 1，j + 1]+ （1 - w） × g[i，j + 1] （16）

若Gx、Gy都为0，则说明该点不是边缘点。

步骤7：自适应双阈值处理。首先，对经过非极大值抑制后的所有点进行两两比较，找到梯度值最大值Gmax，根据公式（7）和（8）计算出高阈值TH，低阈值TL。然后，依次对每一个点进行双阈值判断，流程如图5 所示。

2.2 CUDA 并行算法及优化

计算机统一设备架构（CUDA） 是由 NVIDIA 公司推出的一种 GPU 并行计算架构，能够利用 GPU 的计算能力提供一套高效的密集型数据计算解决方案[10]。