周 明， 嚴壯志， 黃 彬

上海大學通信與信息工程學院，上海200072

非線性圖像擴散的格子波爾茲曼(Lattice Boltzmann,LB)方法是近年發(fā)展起來的一種結合建模與數值計算的圖像處理技術，并成功應用于圖像去噪[1-6]、分割[7]以及修復[8]等領域.作為一種數值解法，LB方法物理意義清晰，邊界條件處理簡單，程序易于實現，且具有并行化和穩(wěn)定性等特點，顯然優(yōu)于傳統(tǒng)數值解法.盡管目前LB圖像處理方法較其他非線性偏微分方程的數值實現方法計算效率有所提高，但現有的處理速度遠遠不能滿足復雜圖像(例如醫(yī)學超聲圖像)處理的實時性要求[5,7-8].因此，發(fā)揮LB模型的邊界條件處理簡單、程序易于實現和并行化的特點，進一步提高計算效率是一個重要的研究課題[9].目前，由于多核處理器和GPU可以低成本使用，研究非線性擴散LB模型的并行算法及其在GPU上的實現具有重要的應用價值.

LB模型在GPU上的并行計算一般可分為通用GPU(general purpose GPU,GPGPU)和計算統(tǒng)一設備架構(compute unif ied device architecture,CUDA)兩種，且大多用于計算流體力學領域。前者包括混合流模擬[10]和三維實時流體模擬[11]等，后者包括圓形繞流大渦模擬[12]和二維、三維方腔流模擬[9,13-14]等.這些應用都獲得了較高的GPU-CPU加速比和每秒百萬格子更新(million Lattice update per second,MLUPS).例如，文獻[12]模擬圓形繞流大渦LB模型在nVIDIA GTX275上的計算，相對于CPU可得到18倍的加速比；文獻[13-14]則報道了在nVIDIA GTX280上模擬二維和三維方腔流的LB模型計算，相對于CPU的加速比可達50倍以上，每秒百萬格子更新可達300次以上.在圖像處理領域，僅有少量文獻報道了在GPGPU上實現LB模型的并行處理技術[6]，而通過CUDA實現圖像處理LB方法的研究則尚無文獻報道.

非線性圖像擴散LB模型具有清晰而明確的物理解釋，它以圖像特征和圖像先驗知識作為約束條件，能兼顧自下而上和自上而下的分析方法，故在圖像處理與分析領域具有很好的應用前景.而作為一種并行算法設計架構，CUDA具有易實現、可及性好和GPU-CPU加速比高等特點.為此，本文將聚焦于非線性圖像擴散LB模型，研究其CUDA并行算法的設計與實現.

1 非線性圖像擴散LB模型

文獻[1]提出了基于LB模型的非線性擴散方法并成功地應用于圖像去噪，文獻[2-5]進一步改進了非線性圖像擴散LB模型.本文針對文獻[5]提出的非線性圖像擴散LB模型，研究了CUDA算法.該模型的優(yōu)點在于將圖像邊緣特征嵌入松弛因子中，不但具有良好的邊緣保持能力，而且實現了大步長迭代，提高了計算效率.

假設圖像由離散化的網格組成，每個網格的值由粒子的分布函數Ii(i=0,1,2,···,b-1)和離散速度矢量ci(i=0,1,2,···,b-1)決定.根據ci方向數b的不同，一般把二維LB模型分為D2Q5和D2Q9兩類.LB模型在t時刻位于r處，沿ci方向的演化過程由遷移過程和碰撞過程組成[9]：

遷移過程

碰撞過程

式中，Ii(r,t)為粒子分布函數，(r,t)為平衡態(tài)分布函數，為發(fā)生遷移后的粒子分布函數，Δx為空間步長，Δt為時間步長，ω為松弛因子.

綜合以上兩個過程，可得到LB演化方程

松弛因子ω為

式中，C為迭代步長，g(▽Gσ?I)為在大步長迭代下的邊緣截止函數

式中，Gσ為方差為σ的高斯核；K為平滑閾值，用以控制擴散的范圍.

在D2Q9模型中，平衡態(tài)分布函數為

松弛因子ω為

由式(3)～(8)可以看出，該LB模型的計算具有良好的局部性，非常適用于并行計算.以D2Q5為例，圖像去噪的非線性擴散LB算法流程如下：

確定迭代次數N

While n＜N do

設置初始平衡態(tài)函數：

根據式(5)計算松弛因子

根據式(1)計算遷移過程

根據式(2)計算碰撞過程

更新節(jié)點分布函數

更新平衡態(tài)分布函數

n=n+1

End While

輸出

2 CUDA編程模型及CUDA算法設計

2.1 CUDA編程模型

CUDA程序的執(zhí)行如圖1(a)所示，即讓主機(CPU)中每個核函數(kernel)按照線程網格(grid)的概念在顯卡硬件(GPU)上執(zhí)行；每個線程網格可以包含多個線程塊(block)；每個線程塊又可以包含多個線程(thread)，線程為最小的執(zhí)行單元.在硬件執(zhí)行的時候，如圖1(b)所示，讓每個線程對應一個流處理器(SP)；每個線程塊對應一個多流處理器(SM)；線程網格被以線程塊為單位分配到SM上執(zhí)行.

圖1(c)為GPU的存儲器資源.在Fermi架構下，芯片上有大量線程私有的32位寄存器和線程塊私有的64 kB共享內存，它們具有低延遲性.其中同一個線程塊內共享內存的數據是可共享的.芯片外有全局內存即顯存，雖然數據存儲空間大，但對其訪問有400～600個時鐘周期的延遲.在全局內存中有一塊區(qū)域大小為768 kB的紋理內存，紋理內存具有緩存功能.在CUDA程序中，數據是由CPU傳到GPU的全局內存上.全局內存的優(yōu)化思路是合并對齊訪問隱藏延遲.合并對齊訪問條件為：①線程訪問的數據長度為4、8或16字節(jié)；②被訪問地址構成一片連續(xù)內存空間；③第N個線程訪問第N個全局內存地址；④起始全局內存地址對齊到所存儲數據長度的16倍處.

圖1 CUDA軟硬件架構及GPU存儲器資源[15]Figure 1 CUDA architecture and GPU memory resources[15]

2.2 CUDA算法設計

LB去噪方法的CUDA算法流程如下：

①分別在CPU和GPU上分配內存，初始化CPU和GPU上的相關數據；

②如圖2所示，設置線程網格和線程塊的大小，并在GPU端進行LB去噪方法的并行計算；

③當計算滿足迭代次數n時結束計算，將數據拷貝回CPU做后續(xù)處理，釋放GPU和CPU的內存.

圖2 假設圖像大小為9×9，線程網格設置為(3×9)，線程塊設置為(3×1)，線程數據所對應的內存存儲Figure 2 Working process of global memory when image size is 9×9,grid is 3×9,and block is 3×1

由于LB模型的計算只涉及當前節(jié)點及周圍節(jié)點的信息，而每個格點的計算又相對簡單，因此CUDA程序性能的瓶頸主要是GPU的數據訪問.由2.1節(jié)可知，數據均存儲在全局內存中，是否滿足合并對齊訪問條件是影響程序性能的關鍵.碰撞過程僅涉及當前節(jié)點的計算，數據訪問容易滿足合并對齊訪問條件；而遷移過程需要將計算后的數據進行移動，在需要向左或向右移動的方向，數據訪問的起始地址會發(fā)生偏移(不滿足合并對齊訪問條件的第4條)，如圖3所示.

圖3 遷移過程的數據訪問地址發(fā)生偏移，如線程0的數據遷移到線程1上Figure 3 Diagram of shift during data access in propagation steps:the data of thread 0 shifts to thread 1

2.2.1 CUDA算法

針對上述遷移過程數據訪問的特點，本文設計了以下3種CUDA算法：

算法1 借助紋理內存實現遷移過程的數據訪問.紋理內存具有高速緩存，支持二維尋址，一個數據的“上下左右”數據都能讀入緩存，數據讀取時不需要滿足合并對齊規(guī)則.因此紋理內存十分適合用于圖像處理和非對齊的數據訪問[15].

算法2 借助共享內存實現遷移過程的數據訪問.共享內存只有1個時鐘周期延遲且同一個線程塊內數據可共享，因此比較適合圖像處理以提高效率.

算法3 直接使用全局內存實現遷移過程的數據訪問.這種算法的優(yōu)點是只需一個核函數即可完成，減少了對全局內存的訪問次數，且全局內存資源豐富，可有效提高多流處理器的占用率.

此外，由于在傳統(tǒng)算法2和3中，當采用一塊數組內存時，容易發(fā)生數據覆蓋.如圖4所示，數據由0傳到1，再由1傳到2，使原來1中的數據被覆蓋而產生錯誤.為此，在本文的算法2和3中，分別采用兩個相同大小的數組來儲存格點的分布函數，分別設為A和B.對于奇數迭代步，數據從A中讀取，經LB去噪方法處理后存儲到B；對于偶數迭代步，則反過來執(zhí)行.

圖4 遷移過程的數據覆蓋Figure 4 Data sample covering problem in propagation steps

2.2.2 CUDA算法的實現

圖5具體介紹算法1的實現.首先定義紋理變量，然后使用cudaBind Texture2D()將全局內存中的圖像數據綁定到紋理對象，最后使用tex2D()函數訪問該紋理對象，就可利用紋理內存的特性對圖像數據進行操作.

圖5借助紋理內存實現遷移過程Figure 5 Propagation steps using texture memory

圖6 具體介紹算法2的實現，本文選取了粒子向右方向的遷移.首先將最右端的粒子分布值移到最左端，其次將其余的粒子依次向右移，這樣在遷移過程就滿足合并對齊訪問條件.經過上述處理，雖然滿足了合并對齊訪問，但線程塊兩端的數據沒有正確的遷移.因此，需要另一個核函數來處理線程塊兩端的數據，將每個線程塊最左端的數據移到右邊線程塊的最左端，從而正確地完成了粒子的向右遷移.同樣，算法2也可以通過粒子向左方向的遷移來實現.

圖6 借助共享內存實現粒子向右方向的遷移Figure 6 Propagation steps from left to right using shared memory

算法3的實現按照算法流程即可.數據傳輸到GPU后，根據非線性擴散LB算法流程實施：即先計算松弛因子ω，接著根據式(1)完成遷移過程以及式(2)完成碰撞過程，再更新Ii(r,t)和，滿足迭代次數后停止計算.

3 實驗及討論

3.1 實驗平臺

為了驗證CUDA算法在圖像處理中的可行性，本文搭建一個實驗平臺.平臺包括Xeon W3550的CPU和8 GB內存.為了分析不同GPU的性能，本文分別使用Quadro4000和GT550M兩款GPU在此平臺上做實驗.程序編譯運行環(huán)境為CUDA4.0及VS2010.兩款GPU均為Fermi架構.

在正式實驗之前，本文首先做一個前置實驗來選取合適的線程塊尺寸(線程塊將對GPU計算效率有影響，見圖2). 在前置實驗中選取D2Q9(一種常見的二維LB模型)為對象.其中針對算法1，為了滿足blockDim.y≥2(線程塊在縱坐標的尺寸，因為這個算法使用二維紋理)，選用(32×4)，(64×2)，(64×4)，(32×8)，(128×2)等5種線程塊尺寸對512×512大小的圖像做實驗.針對算法2和算法3，一般將blockDim.y設置為1[9,14]，選取(32×1)，(64×1)，(128×1)，(256×1)等4種線程塊尺寸對512×512大小的圖像做實驗.

圖7給出了本文CUDA算法每次迭代時間與不同線程塊尺寸的關系.從圖7中可以看出，不同的線程塊尺寸，程序性能有差異；算法1～3的最優(yōu)線程塊尺寸分別為(64×2)、(128×1)、(128×1).因此，這3個線程塊將被運用到正式的實驗中去.

圖7 CUDA算法每次迭代時間與不同線程塊尺寸的關系Figure 7 Relationships between caculation time of CUDA algorithms and different blocks

3.2 實驗方法

為了驗證基于LB模型的CUDA算法在圖像處理中的有效性和效率，本文分別針對圖像去噪質量、計算效率和真實醫(yī)學圖像處理效果進行3個實驗.

3.2.1 評價圖像去噪質量的實驗方法

本實驗的目的是定量評價分別利用CPU和GPU實現非線性圖像擴散去噪的質量.實驗采用Lena圖像作為標準圖像來測試.在該標準圖像中加入均值為0和方差為0.01、0.03、0.05、0.07、0.09的高斯噪聲.采用圖像的峰值信噪比(PSNR)作為評價標準.PSNR定義為

式中，I(r)為未加入噪聲的原始圖像，ρ(r,N)為去噪后的圖像，m和n分別為圖像的行數和列數.實驗中分別用CPU和GPU對上述加噪圖像進行處理，以處理后圖像的PSNR為指標對圖像去噪質量進行評價，實驗結果見圖8.

3.2.2 評價計算效率的實驗方法

本實驗用來比較分析上述CUDA算法的性能及兩款GPU實現的程序性能.實驗采用256×256、512×512、1 024×1 024大小的圖片做處理.實驗中針對上述3種大小的圖片，分別用CPU和兩款GPU對它們進行處理，以每次迭代運行時間和MLUPS值為指標，對CUDA算法的計算效率進行評價，實驗結果見表1.

3.2.3 真實圖像去噪驗證的實驗方法

本實驗進一步利用上海腫瘤醫(yī)院超聲科提供的腫瘤圖像(256×256大小)，來驗證基于LB模型的CUDA算法在真實醫(yī)學圖像中的處理效果.為便于比較，將CUDA算法與非線性擴散去噪差分方法中應用廣泛的加性算子分裂算法(AOS)進行實驗對比.實驗中取步長C=2，迭代次數為8，邊緣截止函數閾值K=4.實驗以兩種算法處理圖像后是否存在偽紋理和對圖像邊緣的保持能力為評級指標，處理后圖像的細節(jié)見圖9.

3.3 實驗結果及討論

3.3.1 圖像去噪質量的結果

圖8為CUDA算法處理圖像后獲得的PSNR隨噪聲方差增加的變化曲線.從圖8中可以看出，3種CUDA算法的PSNR隨噪聲方差增加有相同的下降趨勢.由于CPU在硬件執(zhí)行、浮點計算精度方面與GPU存在差異，CPU和GPU上實現的PSNR略有偏差，但總體保持一致，具有相同的去噪效果[15].

圖8 CUDA算法處理圖像后的PSNR與噪聲方差的關系Figure 8 Relationships between PSNR of CUDA algorithms and noise variance

3.3.2 計算效率的結果

表1為D2Q9在CPU和兩款GPU上實現的每次迭代運行時間和MLUPS值，每行為圖像大小，每列為CPU和GPU上實現的每次迭代運行時間和MLUPS值，其中GPU分為Quadro4000和GT550M，兩款GPU下又分為3種CUDA算法(算法1～3).本文中所提到的加速比為CPU和GPU上的運行時間比，如GPU-CPU加速比為CPU和GPU的運行時間比.

表1 D2Q9在CPU和兩款GPU上實現的每次迭代運行時間和MLUPS值Table 1 Calculation time per iteration and MLUPS value of D2Q9 amongst two GPUs and CPU

由表1可以看出，在Quadro4000上實現，算法1的GPU-CPU加速比為70+，MLUPS值為300+；算法2的GPU-CPU加速比最低也有40+，MLUPS值最低也有200+；算法3的GPU-CPU加速比可達90+，MLUPS值可達400+.在GT550M上實現，算法3的GPU-CPU加速比為30+，MLUPS值為150+；算法1和算法2的MLUPS值也有110+.本文的MLUPS值低于流體LB方法[9,14]，原因是LB去噪方法在碰撞過程需進行梯度和濾波運算，數據訪問相對復雜.即便如此，D2Q9在兩款GPU上MLUPS值最低為110+，最高可達400+；而文獻[6]在GPGPU上實現D2Q9的LB去噪方法，其MLUPS值僅為7.相比之下，本文CUDA上實現的計算效率得到了大幅度的提高.

從表1中還可以看出，3種算法的效率由高到低依次是算法3、算法1、算法2.原因如下：算法1和2雖然解決了全局內存的非合并對齊訪問，但算法1中需要紋理綁定，且多個核函數增加對全局內存的訪問次數，影響程序性能.算法2中需要使用if語句來判斷線程塊兩端數據的遷移和邊界處理，這樣就使線程束內的線程走向不同的分支，所需要的時間將是不同分支之和，分支嚴重影響效率.隨著計算規(guī)模的增大，線程塊對共享存儲器使用率增大使得加速比逐漸提高.算法3更高效，原因是使用全局內存可使本去噪方法在同一個核函數內執(zhí)行，避免了多個核函數對全局內存數據的多次訪問；同時線程塊不受資源限制，可使SM中SP的占用率更高.

3.3.3 真實圖像去噪驗證結果

圖9為CUDA算法與AOS處理超聲圖像后的細節(jié)對比，圖9(a)和9(d)為原超聲圖像，圖9(b)和9(e)為AOS處理后的圖像，圖9(c)和9(f)為CUDA算法(LB去噪方法)處理后的圖像.由圖9可以發(fā)現，AOS處理后的圖像會產生偽紋理(圖9(b)和9(e)箭頭所示處)，而CUDA算法不會出現偽紋理(圖9(c)和9(f)箭頭所示處)，具有良好的邊緣保持能力.

3.3.4 關于兩款GPU性能對比的結果和討論

從理論上來說，GPU計算性能主要由單精度計算能力和帶寬決定[16-17].在不考慮帶寬等因素的影響下，GPU計算性能為

式中，μ為換算系數，是一個常數；P為計算性能；f為峰值單精度浮點值.此外，當顯存帶寬提升一倍，約可提升30%的計算性.表2為實驗中選取的兩款GPU的主要技術指標.由表2可以得出Quadro4000和GT550M的峰值單精度浮點值和帶寬之比分別為1.71和3，于是可推出Quadro4000計算性能約為GT550M的2.74倍，與兩款GPU的流處理器數之比相當.

表2 Quadro4000和GT550M主要技術指標Table 2 Main specif ications of Quadro4000 and GT550M

圖9 CUDA算法與AOS處理超聲圖像后的細節(jié)對比Figur e 9 Comparisons of CUDA algorithm with AOS in detail of ultrasound images

由表1可以看出，算法1和3的GPU-GPU加速比在2.7～3.1之間，和上述推論吻合；算法2偏低的原因是if語句使束內線程走向不同分支而串行執(zhí)行，對程序的并行性產生了不確定性因素.因此，在保證程序并行性的前提下，綜合考慮GPU的單精度浮點計算能力和帶寬，不同GPU的計算性能和流處理器數量成正比.

4 結語

本文針對非線性圖像擴散LB模型研究其并行實現的3種CUDA算法，并通過實驗驗證了三種算法在圖像處理中的可行性.驗證指標包括處理后圖像的PSNR，加速比和MLUPS值.實驗結果表明，在保證去噪質量的前提下，CUDA并行算法比已有的CPU算法計算效率大幅提高.其中直接使用全局內存實現遷移過程的CUDA算法(算法3)在效率上明顯優(yōu)于借助紋理內存(算法1)和共享內存(算法2)的CUDA算法.本文還通過醫(yī)學超聲圖像去噪，驗證CUDA算法(LB去噪方法)在大步長迭代下具有良好的保持邊緣能力.此外，本文進一步使用GT550M和Quadro4000做實驗對比，通過以單精度浮點峰值性能和帶寬為指標，驗證了在保證程序并行性的前提下，GPU的計算性能與其流處理器的數目成正比.

[1]JAw ERTHB,LINP,SINZINGERE.Lattice Boltzmann models for anisotropic diffusion of images[J].Journal of Mathematical Imaging and Vision,1999,11(3):231-237.

[2]陳玉，嚴壯志，錢躍竑.基于格子波爾茲曼模型的圖像去噪[J].電子學報，2009,37(3):574-580.

CHENYu,YANZhuangzhi,QIANYuehong.The Lattice Boltzmann method based image denoising[J].Acta Electronica Sinica,2009,37(3):574-580.(in Chinese)

[3]CHANGQ S,YANG T.A Lattice Boltzmann method for image denoising[J].IEEE Transactions on Image Processing,2009,12(18):2797-2802.

[4]ZHANGW H,SHI B C.Application of Lattice Boltzmann method to image f iltering[J].Journal of Mathematical Imaging and Vision,2012,43:135-142.

[5]王志強，嚴壯志，錢躍竑.圖像非線性擴散去噪的格子波爾茲曼方法[J].應用科學學報，2010,28(4):367-373.

WANG Zhiqiang,YAN Zhuangzhi,QIAN Yuehong.Nonlinear diffusion for image denoising using Lattice Boltzmann method[J].Journal of Applied Sciences,2010,28(4):367-373.(in Chinese)

[6]ZHAOY.Lattice Boltzmann based PDE solver on the GPU[J].The Visual Computer,2008,24:323-333.

[7]WANG Z Q,YAN Z Z,CHEN G.Lattice Boltzmann method of active contour for image segmentation[C]//Sixth International Conference in Image and Graphics(ICIG),2011:338-343.

[8]張蕊，嚴壯志，劉瑋.圖像修復的格子波爾茲曼方法[J].電子測量技術，2011,34(3):46-65.

ZHANGRui,YANZhuangzhi,LIUWei.Lattice Boltzmann method based image inpainting[J].Electronic Measurement Technology,2011,34(3):46-65.(in Chinese)

[9]KUZNIKF,OBRECHTC,RUSAOUENG.LBMbased flowsimulationusingGPUcomputingprocessor[J].Computers&MathematicswithApplications,2010,59(7):2380-2392.

[10]朱紅斌，劉學慧，柳有權.基于LatticeBoltzmann模型的液混合流模擬[J].計算機學報，2006,29:2071-2079.

ZHUHongbin,LIUXuehui,LIUYouquan.Binary mixturessimulationbasedonLatticeBoltzmann method[J].ChineseJournalofComputers,2006,29:2071-2079.(inChinese)[

11]柳有權，劉學慧，吳恩華.基于GPU帶有復雜邊界的三維實時流體模擬[J].軟件學報，2006,17:568-576.

LIUYouquan,LIUXuehui,WUEnhua.Real-time3D fluidsimulationonGPUwithcomplexobstacles[J].JournalofSoftware,2006,17:568-576.(inChinese)

[12]ZHOUH,MOGY,WUF.GPUimplementation ofLatticeBoltzmannmethodforflowswithcurved boundaries[J].ComputerMethodsinAppliedMechanicsandEngineering,2012,225-228:65-73.

[13]OBRECHTC,KUZNIKF,TOURANCHEAUB.Anew approachtotheLatticeBoltzmannmethodfor graphicsprocessingunits[J].Computers&MathematicswithApplications,2011,61:3628-363.

[14]TOLKEJ.ImplementationofaLatticeBoltzmann kernelusingthecomputeunifieddevicearchitecture developedbynVidia[J].ComputingandVisualizationinScience,2010,13:29-39.

[15]陳曙暉，熊淑華.大規(guī)模并行處理編程實戰(zhàn)[M].北京：清華大學出版社，2010:106-116.

[16]http://www.realworldtech.com/amd-nvidia-gpu-per formance/.

[17]http://www.realworldtech.com/gpu-memory-bandw idth/.