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 大小的像素點?？梢岳肗umPy中的切片功能完成圖像裁剪。

步驟 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]。

在 CUDA 程序架構中，主機端負責串行程序的運行及整體調度工作，而 GPU 端負責多個線程并行執行內核函數以完成計算任務。GPU 中的線程采用塊狀劃分方式，即所有線程以一維、二維及三維的方式組織。在本文的算法中，線程與圖像像素點一一對應。首先由 CPU 端將數據拷貝至 GPU 端，調用內核函數進行計算，計算完成后將結果拷貝回 CPU 端。基于 CUDA 的 Canny 邊緣檢測算法過程如圖 6 所示。

在圖像處理過程中，一個像素點可以通過與其對應的一個線程進行處理，多個線程可以視作對多個像素點并行處理。CUDA 架構下的線程是一種輕量級線程，適合處理小粒度的并行問題。基于這一特點，可以將圖像中的一個像素點與其對應的處理線程一一映射，從而對問題進行細粒度分解。Canny 邊緣檢測并行算法的詳細步驟如算法2所示。

算法2 ：基于CUDA的并行邊緣檢測算法

輸入：圖像數據，高斯卷積核大小Ksize=5，高斯模糊因子=0.1。

輸出：邊緣二值圖像。

步驟1：將圖像轉換為灰度圖，并進行零填充。

步驟2：生成高斯卷積核。根據輸入參數，生成大小為5×5，模糊因子為0.1的高斯卷積核G（x，y）。

步驟3：配置塊（block） 。設置每個塊的線程數TPB=16，即一個塊的大小為16×16，包含256 個線程（thread） 。調用GPU端計算每個線程的全局坐標，與其將要處理的對應像素點的關系，可表示為：

tdx， tdy = cuda.grid （2） （17）

式中：tdx， tdy為當前線程在全局空間中的位置，即其對應處理像素點在圖像中的位置。

步驟4：分配顯存空間，將灰度圖及高斯核數據傳輸至GPU端中的顯存。

步驟 5：在 GPU 端進行高斯濾波。在串行算法中，濾波前對圖像的上下左右邊界處進行了大小為pad的零填充，濾波完成后需要對圖像進行裁剪以恢復到原始大小。如果將濾波完成后的圖像傳回主機端進行裁剪再傳回 GPU 端，會造成傳輸延遲。因此，需要計算出濾波前后圖像像素點的對應關系，使得在濾波的同時便完成圖像裁剪，從而在 GPU 端完成整個過程，如圖 7 所示。

在圖中，箭頭所指向的即為濾波前后像素點的對應關系，其對應關系可表示為：

Ia [i，j ] = Ib [i + pad，j + pad ] （18）

式中：Ia表示進行高斯濾波后的圖像，Ib表示進行高斯濾波前的圖像，pad=Ksize/2。該步驟最終得到進行高斯濾波并且經過裁剪恢復大小后的圖像gauss_out。

步驟6：在GPU端進行梯度計算。首先在主機端生成Sobel梯度算子，即X方向差分算子與Y方向差分算子，并傳輸至GPU端，計算得到每一個線所對應處理的數據點，并行在GPU端進行X方向梯度算子與圖像的卷積運算及Y方向梯度算子與圖像的卷積運算，分9Y7qzhrCDFhFF6fov8XSYQ==別計算每一個像素點的X方向偏導數及Y方向偏導數。依據公式（5）和（6）計算出當前點梯度幅值及方向。最終得到圖像中所有點的X方向偏導數gradX、Y 方向偏導數gradY以及梯度幅值grad。

步驟 7：在 GPU 端進行非極大值抑制。根據梯度計算得到的結果，計算得到每一個線程所對應處理的數據點。所有線程并行根據梯度計算所得結果判斷當前線程所對應的數據點屬于圖 4 中的何種情況，然后依據該點Gx與Gy的方向是否相同，選取公式 （9） 至公式 （16） 計算沿當前點梯度方向兩側點的線性插值，并與當前點的梯度幅值進行比較。若當前點的梯度幅值較大，則保留該點，否則剔除。該步驟最終得到所有點經過非極大值抑制后的結果NMS。

步驟 8：在 GPU 端進行自適應雙閾值處理。在自適應雙閾值處理中，高閾值TH與低閾值TL的確定依據公式 （7） 和公式 （8）。首先需找出非極大值抑制后結果 NMS中梯度幅值的最大值。在串行算法中，可以直接調用 Numpy 類的max（）函數實現。計算得到每一個線程所對應處理的數據點。在 CUDA 核函數中，通過雙循環結構將 NMS中的值依次進行兩兩比較能夠找到最大值，但該算法的時間復雜度為O（n2） 。經測試發現，若采用此方案，并行算法的耗時不僅沒有降低，反而遠遠大于串行算法的耗時。如果將 NMS數據傳至主機端，通過 Numpy 類的max（）函數完成求取最大值，然后再將數據傳至 GPU 端，會產生傳輸延遲，從而降低算法整體效率。本文設計了一種稱為二維規約的方法來求取梯度幅值的最大值，如圖 8 所示。其中，stride 為規約步長，即每一輪比較中兩個對應比較線程之間的跨度。首先，對圖像中的垂直方向進行規約比較。對應的行中的線程進行兩兩比較，即第i 行中的所有線程分別獲取圖像中本行線程所對應處理的數據點以及圖像中第i+stride行的數據點進行梯度幅值比較。如果第i 行中點的梯度幅值小于第i+stride行中點的梯度幅值，則交換兩行的梯度幅值，否則不進行交換。由于線程是并發執行的，每一輪比較中所有對應線程的比較計算都是并發進行的。當垂直方向規約比較完成后，結果數據的首行中數據為所有行中的最大值。

接著進行水平方向規約比較。第j 個線程分別獲取圖像中本線程所對應處理的數據點以及第j+stride 個數據點進行梯度幅值比較。如果第j 個點的梯度幅值小于第j+stride個點的梯度幅值，則交換數據，否則不進行交換。當水平方向規約比較完成后，結果數據中的首元素為所有數據中的最大值。

在此方案中，垂直方向和水平方向進行規約比較的輪數各為log（n），算法的時間復雜度為O（log（n））。相較于前面的O（n2），時間復雜度大大降低。當得到梯度最大值后，每個線程并行執行后續的雙閾值處理，最終得到經過雙閾值處理的結果DT。

步驟 9：將結果傳送至主機端，輸出邊緣二值圖像。

3 實驗測試與結果分析

3.1 實驗環境

硬件平臺：處理器：Intel（R） Core（TM） i7-10710UCPU @ 1.10GHz 1.61 GHz；顯卡：NVIDIA GeForceMX350。

軟件平臺：操作系統：Windows 10；程序設計語言：Python 3.8.9；開發環境：PyCharm、CUDAToolkit10.2、Numba。

3.2 測試結果與分析

本文選取了圖像處理標準測試圖像集中的 air?plane、peppers 及 lena 三張圖像，并分別生成尺寸為256×256、512×512、1 024×1 024 和 2 048×2 048 共五張圖像。實例如圖 9 所示。

針對不同尺寸的圖像，首先進行 Canny 串行算法的運行耗時統計，然后進行并行算法的運行耗時統計，并計算并行算法相對于串行算法的加速比。為保證測試耗時的準確性，實驗結果取 10 次測試耗時的平均值。實驗結果如表 1 所示。

為直觀展示不同尺寸圖像下串行算法與在CUDA 架構 GPU 環境下并行算法的運行時間，繪制了柱狀圖 如圖10所示。從圖10中可以看出，Canny 并行邊緣檢測算法的運行時間相比于串行算法有了大幅度的降低。

為直觀展現并行算法在不同尺寸圖像下的加速性能，計算出算法在不同圖像尺寸上的檢測加速比，并繪制加速比曲線圖，如圖 11 所示。

從圖 11 中可以看出，在一定的圖像尺寸范圍內，加速比隨著圖像尺寸的增加呈上升趨勢，這表明本文的并行算法表現出了較好的數據可擴展性。在速度獲得大幅度提升的同時，Canny 邊緣檢測并行算法必須確保檢測精度。本文中 Canny 串行算法和并行算法的檢測效果如圖 12 所示。

從圖 12 中可以看出，本文的 Canny 并行算法相對于 Canny 串行算法而言，能夠完整地檢測到串行算法所能檢測到的邊緣，并且保留了更多的細節邊緣信息。這表明本文的并行算法能夠滿足檢測精度的要求。

4 結論

本文首先在 Python 環境中實現了自適應雙閾值的 Canny 邊緣檢測算法的串行版本，然后設計了一種基于 CUDA 架構的 Canny 邊緣檢測并行算法，并對其進行了優化，通過 Numba 編譯器實現。在尋找梯度幅值極大值的過程中，設計了一種基于 CUDA 架構的并行化二維規約的梯度幅值極大值尋找方法，使整個操作都在 GPU 中完成，全程只有一次主機端與 GPU 端的數據傳輸。實驗結果表明，在保證邊緣檢測精度的同時，基于 CUDA 的 Canny 邊緣檢測并行算法相比于串行算法能夠大幅度提升檢測效率，并且隨著圖像尺寸的增加，算法加速比也呈現出上升趨勢，表現出良好的數據可擴展性。在 2 048×2 048 圖像尺寸下，基于 CUDA 的 Canny 邊緣檢測并行算法的最高加速比可達 95.27，具備良好的數據可擴展性。

參考文獻：

[1] CANNY J.A computational approach to edge detection[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（6）：679-698.

[2] 唐斌，龍文.基于GPU+CPU的CANNY算子快速實現[J].液晶與顯示，2016，31（7）：714-720.

[3] 郭恒亮，柴曉楠，韓林，等.Canny邊緣檢測算法在飛騰平臺上的實現與優化[J].計算機工程，2021，47（7）：37-43.

[4] 劉振濤，燕必希，董明利，等.并行計算在動態攝影測量邊緣提取算法中應用[J].計算機工程與設計，2019，40（1）：97-102.

[5] 吳進，趙雋，李聰，等.機器視覺中邊緣檢測算法的SDSoC加速實現[J].計算機工程與應用，2019，55（12）：208-214.

[6] 藍貴文，吳昊錚，張強，等.一種基于CUDA的改進Canny邊緣檢測算法[J].桂林理工大學學報，2018，38（3）：494-500.

[7] 張帆，韓樹奎，張立國，等.Canny算法的GPU并行加速[J].中國光學，2017，10（6）：737-743.

[8] 查志勛，葉兵.結合圖像增強的實時邊緣檢測系統[J].現代電子技術，2022，45（14）：169-174.

【通聯編輯：唐一東】