基于GPU 的快速變指數(shù)PM 模型及BGA 射線圖像去噪

郭亞寧，白云蛟，張 權(quán)，劉 祎，桂志國

（中北大學(xué) 電子測試技術(shù)國家重點實驗室，山西 太原030051）

0 引 言

球柵陣列結(jié)構(gòu)的印刷電路板（Ball Grid Array，BGA）已經(jīng)廣泛地應(yīng)用于印刷電路板等制造領(lǐng)域，其可靠性直接影響著電路板產(chǎn)品的質(zhì)量．但目前檢測設(shè)備的管電流僅為微安級，而且采用微焦點X 射線源，這就使得射線穿過物體到達探測器的光子數(shù)很少，圖像中含有大量的噪聲，影響了缺陷識別和分析．所以，射線圖像去噪是缺陷檢測的關(guān)鍵步驟．

傳統(tǒng)的去噪方法如高斯濾波法，中值濾波法等，雖然能夠很好地去除噪聲，但模糊了圖像的邊緣和細節(jié)信息．1990年，P．Perona 和J．Malik（PM）首先提出了各項異性擴散模型［1］，有效地解決了保持邊緣和去除噪聲之間的矛盾．從此，研究 人 員 圍 繞PM 模 型 展 開 了 一 系 列 的 研 究［2－4］．2010年，S．M．Chao和D．M．Tsai［5］將局部方差引入到PM 模型中，較好地保持了圖像的細節(jié)信息．2012年，Z．C．Guo等提出了一種變指數(shù)各向異性擴散圖像去噪算法［6］，在有效地去除噪聲的同時，抑制了PM 模型的“階梯效應(yīng)”；但其算法復(fù)雜，尤其是在圖像尺寸較大時，去噪效率大大下降，所以在工程應(yīng)用中受到限制．

目前，隨著計算機硬件的發(fā)展，基于統(tǒng)一計算設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）［7］的硬件加速方法，越來越成為國內(nèi)外研究人員的研究熱點［8－10］．其中GPU 加速技術(shù)最典型的應(yīng)用就是CT 重建算法加速［11－13］，并且已經(jīng)取得很好的效果．本文充分利用文獻［6］中的變指數(shù)PM 模型優(yōu)越的去噪性能，結(jié)合GPU 加速技術(shù)，提出了一種基于GPU（Graphic Processing Unit）的快速變指數(shù)PM 去噪方法．該方法引入CUDA 并行架構(gòu)，使得運算過程在GPU 上并行化，縮短了算法的運行時間．最后采用BGA 射線圖像進行了實驗分析和驗證．實驗結(jié)果表明：所提方法在不影響去噪效果的前提下，提高了去噪效率．

1 基于變指數(shù)的PM 降噪模型

1．1 變指數(shù)PM 模型

原始PM 模型（擴散模型）雖然能夠去除噪聲，但會產(chǎn)生“階梯效應(yīng)”，而熱擴散方程在去噪的同時又會模糊邊緣細節(jié)等信息．因此，綜合考慮兩種模型的優(yōu)缺點，文獻［6］提出了新的擴散模型

其中，邊緣檢測算子α（x）取如下形式：

式中：f 為原始圖像；I 為演變圖像的梯度；▽為梯度算子；Gσ為標準差為σ 的高斯函數(shù)；K 為閾值常數(shù)；κ 為常數(shù)；λ 為調(diào)節(jié)系數(shù)，用來調(diào)節(jié)原噪聲圖像和平滑后圖像之間的保真度．這里，λ取

邊緣檢測算子α（x）有如下特性：α（s）＝2－2／（1＋κs）是一個遞增函數(shù)，α（0）＝0，而且在圖像平坦區(qū)域，此時α（x）→0，式（1）就變成了熱擴散方程，對圖像進行各向同性擴散；在圖像的邊界區(qū)域，｜▽GσI｜2→∞，此時α（x）→2，式（1）就變成了PM 模型，對圖像進行各向異性擴散．

綜上所述，在新的擴散模型中，在圖像的邊緣處，趨近于PM 模型，保持邊緣；在圖像的平坦區(qū)域，趨近于熱擴散方程．這樣新模型就能夠在去除噪聲的同時更好地保持邊緣和細節(jié)．

1．2 離散化表示變指數(shù)PM 降噪模型

式（1）離散化可表示為

式中：I0（x，y）和It（x，y）分別為原始圖像和第t次迭代時圖像坐標（x，y）處的灰度值．▽Iit（x，y）（i＝1，2，3，4）分別代表中心像素（x，y）的上、下、左、右4個鄰域的梯度．

Cit（x，y）為相應(yīng)的擴散系數(shù)，且是梯度▽Iit（x，y）的函數(shù)，

式中：K 為常數(shù)，α（x）取式（2），λ取值的離散形式為4

通過上述分析可知，新的擴散模型能很好地去除噪聲，同時保持邊緣細節(jié)信息．但是，該模型的計算量非常大，尤其是計算邊緣檢測算子α（x）和保真系數(shù)λ．因此，該模型隨著迭代次數(shù)的增加，計算量也在成倍地增加，限制了該模型在實際生產(chǎn)中的應(yīng)用．為此，本文引入CUDA 并行架構(gòu)，采用GPU 加速技術(shù)，對該算法進行了加速．高斯濾波；第二步，計算式（4）等號右邊第二項的值以及式（6）分子、分母的值；第三步，把第二步得到的結(jié)果值帶入式（4）中計算出最后圖像各像素的值．為了便于加速分析，首先給出變指數(shù)PM 降噪模型的CPU 實現(xiàn)偽代碼：

2 GPU 加速實現(xiàn)變指數(shù)PM 降噪模型

2．1 CUDA 編程模型

CUDA 程序在執(zhí)行過程中，主機代碼與設(shè)備代碼交替執(zhí)行，程序從主機端的串行代碼開始執(zhí)行，當運行至設(shè)備代碼時，調(diào)用內(nèi)核函數(shù)，并切換到設(shè)備端，啟動多個GPU 線程并行執(zhí)行設(shè)備代碼．設(shè)備完成計算后，返回主機線程，主機繼續(xù)執(zhí)行串行操作．重復(fù)以上過程，直到程序執(zhí)行完畢．編寫流程如圖1 所示．

圖1 CUDA 程序的編寫流程Fig．1 Writing process of CUDA procedure

2．2 GPU 加速實現(xiàn)過程

由上一節(jié)分析可知，變指數(shù)PM 降噪模型的計算主要分為三個過程：第一步，對原圖像進行

for（1：迭代次數(shù)）——迭代循環(huán)

｛

for（1：h－1）

for（1：w－1）｛do；｝——第一步 for（1：h－1）

for（1：w－1）｛do；｝——第二步

do；——求出保真項系數(shù)λ for（1：h－1）

for（1：w－1）｛do；｝——第三步 do；——邊界像素處理

｝

其中h和w 為圖像的高度和寬度．由上述代碼分析可知，總的計算過程可以看作是由一個大循環(huán)內(nèi)嵌套著三個嵌套的小循環(huán)組成．由于最外層的大循環(huán)是用來完成迭代計算的，而迭代計算的特點就是當前運算要用到上一步得到的值，實際上就是一個先后串行的執(zhí)行過程．因為CUDA并行架構(gòu)的主要出發(fā)點是為了處理相互之間并行計算且彼此之間沒有任何聯(lián)系的任務(wù)，所以可知，對迭代計算過程而言，并不適于并行運算，因此，這里的并行分解過程只能對該模型的內(nèi)部循環(huán)展開．

進一步分析上面的代碼，可知內(nèi)部三個嵌套的小循環(huán)之間存在先后順序：只有當?shù)谝徊街醒h(huán)執(zhí)行完后才能執(zhí)行第二步循環(huán)，當?shù)诙降难h(huán)執(zhí)行完后才能執(zhí)行第三步的循環(huán)．由于這些性質(zhì)，這里的并行分解只能逐步地展開：先對第一步進行循環(huán)分解，然后對第二步進行循環(huán)分解，最后循環(huán)分解第三步．

2．2．1 對第一步高斯濾波算法的CUDA 并行分解

高斯濾波器是用待處理像素相鄰范圍內(nèi)的幾個像素的加權(quán)平均值來代替該點的像素值，而每個相鄰的像素點的權(quán)值是根據(jù)該點與待處理點的距離決定的，通常是單調(diào)增減的．對一塊數(shù)據(jù)進行高斯濾波處理時，算法對每個像素的處理無任何數(shù)據(jù)相關(guān)，所以本文中把對每個像素的處理映射到每個線程中．

本文的濾波過程使用二維零均值離散高斯函數(shù)作平滑濾波器．由于高斯函數(shù)的可分離性，二維高斯函數(shù)卷積可以分兩步進行：首先將圖像與一維高斯函數(shù)進行卷積，然后將卷積結(jié)果與方向垂直的相同一維高斯函數(shù)卷積．在GPU 上實現(xiàn)高斯濾波時，不能將這兩步合并為一個內(nèi)核處理過程，因為對圖像進行Y 方向濾波時需要用到經(jīng)X方向的濾波結(jié)果，所以在GPU 上處理時需要寫兩個內(nèi)核函數(shù)，分別對兩個垂直方向進行處理．先調(diào)用對圖像X 方向濾波的內(nèi)核函數(shù)，然后進行CPU 與GPU 的同步，確保X 方向的內(nèi)核濾波函數(shù)完成后，再調(diào)用對圖像Y 方向濾波的內(nèi)核函數(shù)．

由于在高斯濾波的處理過程中，濾波器的系數(shù)是不會改變的，而且每個線程的處理都要用到該濾波器．為了提高運行效率，要充分利用硬件里的更快的存儲器，本文選擇把一維濾波器的內(nèi)容傳遞到常量存儲器里面，這里需要用到函數(shù)cudaMemcpyToSymbol（），常量存儲器帶有緩存機制，這樣能加速每個線程訪問高斯濾波器，節(jié)約了算法的運行時間．

2．2．2 對第二步過程的CUDA 并行分解

上述所示變指數(shù)PM 降噪模型的CPU 實現(xiàn)偽代碼中第二部分的計算，主要是為了求出式（4）等號右邊第二項的值，以及式（6）分子、分母的值，以方便下一步計算得到保真項系數(shù)λ．觀察式（4）和式（6）可以發(fā)現(xiàn)，兩個公式的計算過程都是以像素為單位進行的，而且在式（4）中等號右邊第二項算出一個值的同時，式（6）中分子和分母對應(yīng)像素的某一項值也能計算出來．這樣，就可以把計算式（4）中等號右面第二項的值和計算式（6）中分子和分母對應(yīng)項的值放在同一個循環(huán)體中進行，等循環(huán)完成后，各個像素對應(yīng)的值都能夠算出來．再觀察式（6）可知，分子和分母都是上一步求的值的總和，為了求得保真項系數(shù)λ，需要先分別做求和運算，再做除法運算．對于求和運算，在CPU代碼實現(xiàn)時，可以每求得一項值后，緊接著累加到一個變量中，由于CPU 循環(huán)中實際上執(zhí)行的是串行運算，所以當循環(huán)結(jié)束后，最后的累加和就是所有各項之和．但是，如果想在GPU 端得到相同的結(jié)果，就必須并行計算得到所有各項的值后，再進行累加．因為GPU 端代碼執(zhí)行時，一個線程負責計算一個像素的值，各線程之間并行執(zhí)行，所以就不能像在CPU 端執(zhí)行一樣，算出結(jié)果后直接累加到一個變量上，否則得到的結(jié)果很可能不對．因此，在對第二部分做并行分解時，必須寫成兩個內(nèi)核函數(shù)，一個用來求各項的值，另一個用來求各項值的累加和．其中，在求各項值的累加和時使用并行歸約求和算法，充分利用共享存儲器可以被同一個線程塊中所有線程共享的特點，并且進行了深層次的優(yōu)化，因為在CUDA 存儲器模型中，GPU 能夠在4個時鐘周期內(nèi)訪問共享存儲器，這遠低于訪問全局存儲器的400～600個時鐘周期，所以這樣可以極大地提高速度．

2．2．3 對第三步過程的CUDA 并行分解

在第三個循環(huán)體中，只需要把上兩步得到的結(jié)果帶入式（4）后，計算得到的值就為新模型下圖像的運算結(jié)果．由于每次計算都是以像素為單位進行的，各像素之間的運算相互獨立，互不影響，所以可以通過開啟與圖像像素數(shù)等量的線程來施行并行分解運算．本文中采用的具體方法是開啟與圖像高度（＜1024）相同的線程塊數(shù)block，即開啟h個block，每個線程塊中開啟與圖像寬度（＜1 024）相同的線程數(shù)，即每個block中開啟w 個線程．這樣，程序就可以根據(jù)圖像的實際大小開啟適量的線程，每個線程負責計算一個像素值，可以獲得很高的加速效果．

至此，整個加速過程分析完成，GPU 端代碼就是按照以上給出的各個核函數(shù)的順序依次執(zhí)行，最后把計算結(jié)果輸送到主機端．

3 實驗結(jié)果與分析

為了驗證所提加速方法的效果和正確性，本文使用表1 所示測試環(huán)境對大小為508508 和10161016的BGA 射線圖像進行了測試．表2 給出了變指數(shù)PM 模型在CPU 和GPU 下不同迭代次數(shù)時的運行時間．

表1 測試環(huán)境參數(shù)Tab．1 Parameters of test environment

按行分析表2 可知，在相同迭代次數(shù)的情況下，當圖像大小由508508 變成10161016 時，CPU 的運行時間變?yōu)樵瓉淼? 倍，而應(yīng)用GPU的運行時間變?yōu)樵瓉淼?倍，這是因為基于CPU的運算是串行運算，而基于GPU 的運算是并行運算，同時執(zhí)行二分之一的運算量，而且加速比也成為原來的2倍．可見，圖像越大，加速比越大，本文算法的加速優(yōu)勢越明顯．當圖像的大小相同時，CPU 和GPU 的運行時間隨著迭代次數(shù)的增加而增加，但加速比基本穩(wěn)定．可見，本文算法的加速優(yōu)勢沒有受到迭代次數(shù)的影響．總之，本文所提方法能夠有效地提高去噪效率．

表2 CPU 與GPU 處理變指數(shù)PM 模型的時間Tab．2 Time of variable exponent PM model processed by CPU and GPU

為了進一步說明算法的去噪性能，圖2 給出了大小為1 016×1 016的BGA 射線圖像，圖3給出了圖2 中圖像黑框區(qū)域的放大圖像，以及該區(qū)域經(jīng)過PM 模型，文獻［5］中模型，變指數(shù)PM 模型以及本文所提模型處理后的放大圖像．其中，本文的快速變指數(shù)模型迭代次數(shù)為5 次．通過圖3（b），（c）和（d）可以看出，變指數(shù)模型有效地抑制了PM 模型的“階梯效應(yīng)”和文獻［5］中模型存在的“斑點”，在保持邊緣的同時很好地去除了噪聲．比較圖3（d）和（e）可知，兩幅圖像的去噪效果基本相同，其平均絕對誤差小于10－5（由于GPU不支持雙精度數(shù)據(jù)計算，所以計算結(jié)果與CPU 相比存在微小差距）．可見，基于GPU 的快速變指數(shù)PM 模型能夠在保持有效去除噪聲的同時加快算法的運行速度，是一種實用的快速去噪方法．

圖2 BGA 射線圖像Fig．2 BGA X－ray image

圖3 不同方法對BGA 射線圖像去噪的局部放大圖Fig．3 Local enlarged image of BGA X－ray image using different methods

4 結(jié) 論

在BGA 射線圖像去噪中，為了充分利用變指數(shù)PM 模型能夠有效去除并且保持邊緣的優(yōu)勢，提出了一種基于GPU 的快速變指數(shù)PM 去噪方法．本文采用GPU 硬件加速工具，利用多線程并行運算功能，對該算法進行了加速．實驗證明，在達到與基于CPU 算法相同降噪效果的同時，提高了算法的去噪效率．

［1］Perona P，Malik J．Scale－space and edge detection using anisotropic diffusion［J］．IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12（7）：629－639．

［2］夏石川，桂志國，張權(quán)，等．基于偏微分方程的BGA射線圖像去噪方法［J］．中北大學(xué)學(xué)報（自然科學(xué)版），2013，34（6）：667－672．Xia Shichuan，Gui Zhiguo，Zhang Quan，et al．Denoising methods based on partial differential equation for BGA X－Ray image［J］．Journal of North University of China（Natural Science Edition），2013，34（6）：667－672．（in Chinese）

［3］Yang M，Liang J，Zhang J，ea al．Non－local means theory based Perona－Malik model for image denosing［J］．Neurocomputing，2013，120：262－267．

［4］翟東海，魚江，段維夏，等．米字型各向異性擴散模型的圖像去噪算法［J］．計算機應(yīng)用，2014，34（5）：1494－1498．Zhai Donghai，Yu Jiang，Duan Weixia，et al．Improved image donoising algorithm using UK－flag shaped anisotropic diffusion model［J］．Journal of Computer Applications，2014，34（5）：1494－1498．（in Chi－nese）

［5］Chao S M，Tsai D M．An improved anisotropic diffusion model for detail and edge－preserving smoothing［J］．Pattern Recognition Letters，2010，31（13）：2012－2023．

［6］Guo Z C，Sun J B，Zhang D Z，et al．Adaptive Perona－Malik model based on the variable exponent for image denoising［J］．IEEE Transactions on Image Processing，2012，21（3）：958－967．

［7］劉金碩，鄧娟，周崢，等．基于CUDA 的并行程序設(shè)計［M］．北京：科學(xué)出版社，2014：46－47．

［8］Yang C T，Huang C L，Lin C F．Hybrid CUDA，OpenMP，and MPI parallel programming on multicore GPU clusters［J］．Computer Physics Communications，2011，182（1）：266－269．

［9］Thurley M J，Danell V．Fast morphological image processing open－source extensions for GPU processing with CUDA［J］．IEEE Journal of Selected Topics in Signal Processing，2012，6（7）：849－855．

［10］Jacobsen D A，Senocak I．Multi－level parallelism for incompressible flow computations on GPU clusters［J］．Parallel Computing，2013，39（1）：1－20．

［11］Ehlke M，Ramm H，Lamecker H，et al．Fast Generation of Virtual X－ray Images for Reconstruction of 3D Anatomy［J］．IEEE Transactions on Visualization and Computer Graphics，2013，19（12）：2673－2682．

［12］Blas J G，Abella M，Isaila F，et al．Surfing the optimization space of a multiple－GPU parallel implementation of a X－ray tomography reconstruction algorithm［J］．Journal of Systems and Software，2014，95：166－175．

［13］Beasley D G，Marques A C，Alves L C，et al．GPUaccelerated reconstruction methods for Proton Induced X－Ray Emission Tomography［J］．Radiation Physics and Chemistry，2014，95：251－253．