Kirsch算子地質(zhì)圖像邊緣檢測算法并行化研究?

田宸瑋，王雪純，楊嘉能，錢育蓉?

(1. 新疆大學(xué) 軟件學(xué)院，新疆 烏魯木齊830008；2. 新疆維吾爾自治區(qū)信號檢測與處理重點(diǎn)實(shí)驗(yàn)室，新疆 烏魯木齊830008；3. 新疆大學(xué) 軟件工程重點(diǎn)實(shí)驗(yàn)室，新疆 烏魯木齊830008；4. 中國石油吐哈油田公司勘探開發(fā)研究院，新疆 哈密839009)

0 引言

近年來，飛速發(fā)展的數(shù)字圖像處理技術(shù)普遍應(yīng)用于各個(gè)領(lǐng)域，在地質(zhì)勘探領(lǐng)域中也展現(xiàn)出巨大的潛力，用于檢測及分析各種地質(zhì)學(xué)信息[1,2]．在視頻圖像處理領(lǐng)域，邊緣檢測也是較為基本的方法，常見的邊緣檢測算子包括Sobel算子、Prewitt算子、Roberts算子、LoG算子、Kirsch算子、Canny算子[3,4]．Sobel算子和Prewitt算子在邊緣定位方面有較好的效果，但是像素邊緣容易產(chǎn)生過量[5?7]；Roberts算子同樣在定位邊緣時(shí)有較高的準(zhǔn)確度，但是在抗噪聲方面表現(xiàn)較差[8]；LoG算子在噪聲抑制時(shí)會平滑較為尖銳的邊緣，所以在弱邊緣檢測方面表現(xiàn)較差[9]；Canny算子檢測具有較高的檢測準(zhǔn)確率以及信噪比，且圖像邊緣的連續(xù)性和完整性較好，但是偽邊緣的存在、抗噪聲能力較差等問題依舊存在[10,11]；Kirsch算子的邊緣連續(xù)性和完整性較差，優(yōu)點(diǎn)在于能較好的抑制噪聲[12]．因此，本文提出CUDA（Compute Unified Device Architecture，計(jì)算統(tǒng)一設(shè)備架構(gòu)）下基于自動(dòng)閾值的Kirsch算子邊緣檢測算法．

Kirsch算子通過對原始圖像進(jìn)行8個(gè)方向的邊緣檢測，對這8個(gè)方向的噪聲具有一定平滑的作用．該算法需要對每一個(gè)像素進(jìn)行相應(yīng)的計(jì)算[13]，屬于典型的SIMT（Single Instruction Multiple Threads，單指令多線程）計(jì)算模式．考慮到邊緣檢測算子需要計(jì)算整個(gè)圖像、處理效率不理想的問題，使用CUDA用于設(shè)計(jì)和優(yōu)化并行算法．近年來，CUDA表現(xiàn)出很強(qiáng)的并行計(jì)算性能．GPU架構(gòu)是圍繞具有彈性的SM（Streaming Multiprocessor，流多處理器）陣列構(gòu)建的，在邊緣檢測方面應(yīng)用也較為廣泛，如Adhir Jain[14]使用CUDA進(jìn)行Sobel算子的邊緣檢測；Yam-Uicab R[15]使用CUDA加速霍夫變換；張晗[16]基于CUDA使用Prewitt算子檢測地質(zhì)圖像邊緣．本文通過驗(yàn)證算法的計(jì)算效率以及在并行傳輸開銷方面對所提出算法進(jìn)行測試．

1 改進(jìn)后的邊緣檢測算法

Kirsch算子的邊緣檢測算法是對圖像像素使用順時(shí)針循環(huán)求梯度的方法進(jìn)行邊緣檢測，將其中方向響應(yīng)最大者作為邊緣幅度圖像的邊緣．八個(gè)方向的卷積模板如下所示:

邊緣梯度選取大小為：

其中：i方向?qū)?yīng)的模板為Wi；坐標(biāo)（m,n）處的中心3?3像素區(qū)域?yàn)镻mn．

Kirsch算子的邊緣檢測算法要對原始圖像的每一個(gè)像素點(diǎn)進(jìn)行卷積運(yùn)算，故存在邊緣問題（最外圍一圈像素，周圍像素個(gè)數(shù)不為8），使用擴(kuò)展法來減少判斷邊緣的計(jì)算量，在圖像邊緣填充模板尺寸二分之一的像素邊緣，根據(jù)式（1）可知，計(jì)算單個(gè)像素的響應(yīng)度所需的總計(jì)算量為：加法56（7?8）次，乘法16（2?8）次，還有求比較最大響應(yīng)的7次比較運(yùn)算，計(jì)算一副圖像像素為x?y的圖像計(jì)算量為：加法56（x?y）次，乘法16（x?y）次，比較運(yùn)算7（x?y）次，故不符合實(shí)時(shí)圖像處理系統(tǒng)的要求．為減少計(jì)算量，通過模板Wi對每個(gè)像素區(qū)域進(jìn)行卷積求和的方法，得出計(jì)算各個(gè)方向的二范數(shù)公式：

其中：Fi(m,n)為P(m+j,n+k)與Wij,k卷積運(yùn)算得到的二范數(shù)；Wij,k為i方向?qū)?yīng)模板坐標(biāo)(j,k)處的元素；i=(0,π/4,π/2,3π/4,π,5π/4,3π/2,7π/4)．

最后按照垂直組合的規(guī)則將8個(gè)方向結(jié)合成四組，分別計(jì)算梯度的二范數(shù)：

其中：i=(1,2,...8),ji=π/4?(i?1), ki= mod((3π/2+π/4?(i?1)),2π)，取‖Gi‖2i=(1,2,...8)的最大值為該點(diǎn)像素的梯度，即：

傳統(tǒng)的Kirsch算子在保持邊緣細(xì)節(jié)和抑制噪聲方面有很好的效果，但邊緣連續(xù)性和完整性較差．針對Kirsch算子的缺點(diǎn)，使用整幅圖像的平均灰度作為閾值對其進(jìn)行二值化，閾值(TH)如式(5) 所示：

其中：M為圖像寬度，N為圖像高度．

遍歷整個(gè)圖像，將每個(gè)像素點(diǎn)與得到的閾值進(jìn)行對比后更新，更新規(guī)則如式(6)所示：

分析(1)～(5)式可知，算法執(zhí)行過程中對每一個(gè)像素點(diǎn)進(jìn)行獨(dú)立操作，故存在并行的可能性．設(shè)計(jì)較好的并行算法可達(dá)到減少計(jì)算時(shí)間、增大算法實(shí)時(shí)可能性的目的．

2 基于CUDA對改進(jìn)后的算法并行分析

2.1 算法總流程

每一個(gè)像素使用單個(gè)線程進(jìn)行計(jì)算，多線程可以實(shí)現(xiàn)大量像素的并行計(jì)算，極大提高了計(jì)算效率．閾值自適應(yīng)的Kirsch算子邊緣檢測并行計(jì)算流程如圖1所示．

圖1 算法總流程Fig 1 Algorithm process

Step1Host讀入圖像到內(nèi)存；

Step2Host將數(shù)據(jù)發(fā)送給Device全局內(nèi)存；

Step3Host啟動(dòng)Device核函數(shù)，此時(shí)Device多線程執(zhí)行；

Step4 Device處理結(jié)果回傳到Host內(nèi)存；

Step5Host生成最終圖像．

在Step3中，處理的數(shù)據(jù)為圖像，圖像是由二維點(diǎn)集構(gòu)成的，根據(jù)圖2所示映射結(jié)構(gòu)創(chuàng)建二級線程模式，Grid為二維線程格子，Block為二維線程塊．

圖2 二級線程模式Fig 2 Secondary thread mode

2.2 累計(jì)灰度值以及計(jì)算閾值的并行化

該步驟主要設(shè)計(jì)思想是累加各個(gè)線程塊計(jì)算所得到的局部閾值，最后匯總得到總的閾值．具體實(shí)現(xiàn)思路是先把Block內(nèi)所有的線程從全局內(nèi)存中將圖像像素讀取到共享內(nèi)存，此時(shí)塊內(nèi)的線程需要進(jìn)行同步操作來確保全部線程已經(jīng)讀取像素，最后通過累加各個(gè)線程塊得到的局部閾值來確定最終的全局灰度閾值．

2.3 算法優(yōu)化部分

二維卷積操作同樣會遇到與一維卷積操作一樣的邊界問題．本文使用拓展法（如圖3所示，選取到的二維像素點(diǎn)在周圍添加卷積核大小二分之一的像素填充并賦值為0）來解決邊界判斷問題．使用此方法卷積操作無需考慮邊界問題，但是需要調(diào)整計(jì)算索引．算法1展示了主要的處理流程．

圖3 填充邊緣Fig 3 Fill the edge

算法1：

（1）設(shè)置共享內(nèi)存空間

（3）獲取線程塊中線程坐標(biāo)

（4）int tx = threadIdx.x;

（5）int ty = threadIdx.y;

（6）每個(gè)線程相對于所有線程的坐標(biāo)

（9）使用拓展法，調(diào)整正確索引

（12）int size = width * height;

（13）for (循環(huán)圖像的通道數(shù)) {

（14）讀取每個(gè)像素點(diǎn)并將拓展的像素進(jìn)行賦值

（17）}

（18）else{

（19）Data[ty][tx] = 0.f;

（20）}

（21）同步等待所有線程完成操作

（23）計(jì)算四對二范數(shù)，取最大值作為梯度值，完成后將結(jié)果寫回全局內(nèi)存

（25）}

3 實(shí)驗(yàn)結(jié)果與分析

計(jì)算機(jī)平臺包含支持CUDA加速的顯卡．主機(jī)系統(tǒng)配置為IntelCorei7-7820 HK，時(shí)鐘頻率2.90 GHz，內(nèi)存大小16 GBRAM，GPU規(guī)格見表1．

表1 NVIDIA GeForceGTX 1070 規(guī)格Tab 1 NVIDIA GeForceGTX 1070 parameter

3.1 線程塊索引與數(shù)據(jù)傳輸開銷分析

調(diào)整線程塊索引主要針對數(shù)據(jù)在共享內(nèi)存時(shí)使用，為了降低線程的發(fā)散程度，表2所示分別在全局內(nèi)存和共享內(nèi)存上做了對比，說明本算法使用共享內(nèi)存是更好的選擇．主要原因是：1.硬件方面，共享內(nèi)存位于設(shè)備上，而使用全局內(nèi)存計(jì)算需要頻繁地在主機(jī)和設(shè)備之間執(zhí)行拷貝操作，這無疑大大增加了算法執(zhí)行的時(shí)間；2.軟件層面，通過線程索引的調(diào)整降低了線程的發(fā)散程度．

表2 全局內(nèi)存和共享內(nèi)存讀取時(shí)間對比Tab 2 Global memory and share memory read time compare

3.2 線程規(guī)模對計(jì)算性能的影響

NVIDIAGeForce 1070中，由于每次都需要計(jì)算幅度值以及方向，故在單個(gè)線程塊中需要使用三個(gè)區(qū)域，其中兩個(gè)保存上個(gè)任務(wù)的處理結(jié)果和本次的處理結(jié)果，另一個(gè)用來保存梯度方向．所以每個(gè)拓展塊的容量不能超過16 KB．如表1所示，GPU硬件規(guī)定warp大小為32，故選擇線程塊大小為16×16，包含16×16個(gè)線程．

表3 文獻(xiàn)[17]并行與串行算法執(zhí)行時(shí)間(ms)與加速比(倍)Tab 3 Parallel and serial algorithm execution time(ms) and acceleration ratio (times) in literature[17]

表4 文獻(xiàn)[18]并行與串行算法執(zhí)行時(shí)間(ms)與加速比(倍)Tab 4 Parallel and serial algorithm execution time(ms) and acceleration ratio (times) in literature[18]

為驗(yàn)證算法并行處理的優(yōu)勢，選取不同尺寸圖片對比CPU串行執(zhí)行時(shí)間和GPU并行執(zhí)行時(shí)間，并與文獻(xiàn)[14]、文獻(xiàn)[15]的實(shí)驗(yàn)對比，表明本文提出的算法可以較大縮短執(zhí)行時(shí)間、提高檢測效率，見表3～表5．

表5 本文并行與串行算法執(zhí)行時(shí)間(ms)與加速比(倍)Tab 5 Parallel and serial algorithm execution time(ms) and acceleration ratio (times) in this paper

所用實(shí)驗(yàn)數(shù)據(jù)選取不同分辨率的遙感圖像測試并行算法的檢測效果，如圖4所示．本文提出的算法能夠清晰的分辨出地質(zhì)的紋路．

圖4 遙感地質(zhì)圖像Fig 4 Remote sensing geological image

圖5 算法加速效果與傳輸時(shí)間Fig 5 Algorithm execution effect and transmission time

4 結(jié)束語

本文主要完成了CUDA下自適應(yīng)性閾值Kirsch算子邊緣檢測算法并行優(yōu)化研究，對優(yōu)化傳統(tǒng)Kirsch算子的邊緣檢測算法分析，使用拓展法節(jié)省邊緣判斷計(jì)算所消耗的時(shí)間，將原有8個(gè)方向的計(jì)算減少為計(jì)算4個(gè)方向的二范數(shù)；其次，使用平均灰度作為二值化閾值避免人工設(shè)定閾值隨機(jī)性較大的問題．并行方面從線程塊大小的設(shè)置以及傳輸開銷方面分析算法瓶頸．實(shí)驗(yàn)驗(yàn)證如圖5所示，針對圖像尺寸大于2 048*2 048的地質(zhì)圖像，在考慮傳輸開銷的情況下算法的加速效果保持在80倍以上，不考慮傳輸開銷可保持在300倍以上．可見使用CUDA加速的部分瓶頸在于主機(jī)與設(shè)備的傳輸開銷所占用時(shí)間比例過大，之后會對CUDA的傳輸開銷做進(jìn)一步研究．圖像的處理獲得了較為滿意的結(jié)果，提高了算法的執(zhí)行效率，適用于大規(guī)模地質(zhì)圖像的邊緣檢測．