Fermi架構(gòu)下各向異性擴散超聲斑點噪聲抑制

何興無

(1.成都師范學(xué)院網(wǎng)絡(luò)與信息管理中心，成都 611130;2.成都農(nóng)業(yè)科技職業(yè)學(xué)院電子信息分院，成都 611130)

1 概述

在醫(yī)學(xué)影像技術(shù)中，醫(yī)學(xué)超聲診斷技術(shù)因為其易用性，無侵入性，實時性，廉價性等特點和優(yōu)勢在現(xiàn)代的醫(yī)學(xué)診斷中，有著難以取代的作用［1］。然而在超聲回波信號成像處理過程中，由于需要對位于成像介質(zhì)內(nèi)的超聲隨機散射點產(chǎn)生的信號進(jìn)行疊加產(chǎn)生圖像信號，導(dǎo)致在超聲圖像上往往會出現(xiàn)被稱為斑點噪聲的不規(guī)則斑點，通常會讓人體組織難以識別，降低圖像的分辨率，從而影響圖像質(zhì)量。因此通過去除斑點噪聲來提高超聲圖像的質(zhì)量具有非常重要意義。

在斑點噪聲抑制算法中，傳統(tǒng)的顯式非線性擴散濾波方式可以達(dá)到比較好的期望效果，但僅限于有限的很小的時間步長范圍內(nèi)，對于全部時長則不穩(wěn)定。半隱式加性算子分裂的方式來離散化擴散方程［2］，使迭代過程對于任意步長都是穩(wěn)定的，并且在邊緣保持和噪聲抑制方面都有非常好的處理效果，但在傳統(tǒng)基于CPU處理時很難滿足實時系統(tǒng)的要求，這阻礙了其在實際臨床中的應(yīng)用。隨著多核時代的降臨，并行處理技術(shù)日益成為解決計算瓶頸的最重要手段，一些復(fù)雜計算也可以通過并行處理平臺進(jìn)行顯著加速。

隨著GPU技術(shù)的發(fā)展，其浮點數(shù)處理能力大大超過了CPU，特別是Fermi架構(gòu)GPU的出現(xiàn)，使這一新的處理平臺的計算能力顯著提高。相比較于早期的通用計算GPU設(shè)備，在保持圖形性能的前提下，F(xiàn)ermi架構(gòu)將通用計算技術(shù)提升到前所未有的高度，具體表現(xiàn)在［3］:計算資源更為豐富，每個多處理器擁有32個CUDA(Compute Unified Device Architecture，統(tǒng)一計算設(shè)備架構(gòu))核心(最基本的運算單元);其次，F(xiàn)ermi架構(gòu)引入了緩存機制，使GPU擁有真正意義的可讀寫L1緩存和L2緩存;浮點計算精度和速度也有大幅提高。

與傳統(tǒng)CPU方式不同，對于CUDA并行處理程序設(shè)計，要獲得比較高的性能，需要考慮幾個方面:首先是存儲器操作，CPU的內(nèi)存存取延遲主要是通過多級緩存來消除，而CUDA主要通過高度并行化和合并訪問的方式來達(dá)到高的帶寬利用率，即使在新一代Fermi架構(gòu)下讓CUDA程序盡可能滿足合并訪問，也會明顯地提高性能;其次，并行處理方式的設(shè)計，即如何提高運算并行度，盡可能的使GPU占用率高。

第一部分已經(jīng)介紹所作研究的相關(guān)背景和技術(shù)介紹，在第二部分中將詳細(xì)闡述設(shè)計的并行算法。實驗數(shù)據(jù)的結(jié)果及相關(guān)分析討論放在的第三部分;最后在第四部分給出了工作總結(jié)和對未來工作的展望。

2 Fermi架構(gòu)平臺下基于各向異性擴散斑點噪聲抑制并行處理算法

2.1 各向異性擴散斑點噪聲抑制算法概述

各向異性擴散斑點噪聲抑制算法的核心是基于局部相干性，并且采用加性算子分裂方法進(jìn)行離散化完成快速斑點噪聲抑制?；诰植肯喔尚孕畔⒌陌腚[式各向異性擴散的超聲圖像斑點噪聲抑制算法［3］主要包括低通濾波、擴散系數(shù)表和三對角矩陣的生成及解三對角線性方程組三個環(huán)節(jié)。首先，各向異性擴散方程定義如下:

其中，擴散系數(shù)l(·)定義如下:

μ1和μ2是每點處結(jié)構(gòu)矩陣的特征值，結(jié)構(gòu)矩陣的定義為:J(▽Iσ)=(▽Iσ·▽ITσ)。▽Iσ是通過與一個二維高斯核進(jìn)行卷積(方差為σ)得到的模糊的原圖像I的梯度。把(1)式轉(zhuǎn)化為一個矩陣和向量的形式，并且采用加性算子分裂方法離散化，即得如下方程:

其中U是一個KL×KL的單位矩陣(K，L分別為圖像的高和寬)，λ是迭代步長，Tl(It)是一個三對角陣，Tl(It)=［tij(It)］。

此處，l是擴散系數(shù)。N(i)是一個像素在某一個維度上i的兩個相鄰點的集合。h為濾波窗口的大小。

2.2 數(shù)據(jù)準(zhǔn)備階段

利用GPU平臺進(jìn)行通用計算，相比較一般CPU的處理方式，算法的參數(shù)和數(shù)據(jù)都必須首先從主機端傳送到設(shè)備端才可以進(jìn)行處理，主機端和設(shè)備端都有多種具有不同性能特性的存儲器。因此在并行算法設(shè)計與實現(xiàn)時需要綜合考慮以減少數(shù)據(jù)傳輸?shù)拈_銷。本算法中主要的處理數(shù)據(jù)是圖像數(shù)據(jù)，而參數(shù)包括單個的系統(tǒng)參數(shù)如算法閾值，圖像大小等。這些數(shù)據(jù)的存儲器選擇策略如下:

在主機端，主要進(jìn)行的處理數(shù)據(jù)是每一幀圖像。而在CUDA平臺上為內(nèi)存開放了兩種類型的存儲器，分頁內(nèi)存和頁鎖定內(nèi)存。頁鎖定內(nèi)存的帶寬會比分頁內(nèi)存高出2倍多。頁鎖定內(nèi)存一共提供了4種工作模式。其中，寫結(jié)合模式可以使數(shù)據(jù)在通過PCI－e總線傳輸時不會被監(jiān)視，這能夠獲得高達(dá)40%的傳輸加速，是最適合CPU是只寫的情況。由于輸入數(shù)據(jù)對于CPU來講是只寫的，因此可以使用頁鎖定內(nèi)存中的這種寫結(jié)合模式［6］。

在設(shè)備端，由于Fermi架構(gòu)對全局存儲器引入了緩存機制，因此可以使用二維對齊的線性顯存空間存放來自主機端的輸入數(shù)據(jù)。另外在實現(xiàn)中，可以將一些預(yù)定義好的單個參數(shù)合并入一個數(shù)組放入常量存儲器供所有線程訪問，這樣不僅可以充分利用常量存儲器帶寬，也可以減少寄存器的占用率以提高整體運行的并行度。

2.3 低通濾波

在各向異性擴散方法中每次迭代都要先對原始圖像進(jìn)行一個低通濾波處理，常用的低通濾波器是高斯濾波器。它通常使用一個高斯核卷積實現(xiàn)。這里使用的是5×5的高斯核。對于一個2維的高斯卷積，為了減少冗余計算，可以將其分解為兩個1維高斯核分別沿X，Y軸作卷積的結(jié)果。同時，可以用一個均值濾波代替高斯濾波。通過測試，兩者的處理效果基本一致。但在CUDA程序設(shè)計時，采用均值濾波的處理能力更快，具體的并行實現(xiàn)方式是啟動兩個核函數(shù)分別進(jìn)行沿X方向和沿Y方向的處理。讓一個線程負(fù)責(zé)處理一行或一列，X方向的示意如下:

這里需要注意的是在沿X方向處理時，相鄰線程訪問的圖像像素點是隔開的，并不滿足合并訪問規(guī)則，全局存儲器在這種情況下效率極低。因此在進(jìn)行X方向濾波時，圖像數(shù)據(jù)應(yīng)綁定到紋理存儲器，然后將X方向一維濾波輸出寫入到全局存儲器中供Y方向濾波使用。在進(jìn)行沿Y方向濾波時，由于相鄰線程同時讀取圖像一行的像素值(物理上是相鄰并且對齊存儲的)，所以滿足了全局存儲器的合并訪問條件，實現(xiàn)最優(yōu)化訪存。

2.4 擴散系數(shù)表和三對角矩陣的生成計算

由公式(3)可知，本算法需要進(jìn)行多次迭代才能得到比較好的去噪效果，在每一次迭代中，擴散系數(shù)是不同的。如果在每次迭代過程中都要計算擴散系數(shù)，就會大大降低算法的效率［7］。于是在設(shè)計中采用生成擴散系數(shù)表，每次迭代用查表的方式得到擴散系數(shù)，這樣就減少了大量的冗余計算。因為擴散系數(shù)與圖像上某點的梯度有關(guān)，而梯度分量的取值范圍為－255～255，所以生成的擴散系數(shù)表的大小為511×511。同時，利用式(4)查表可以得到三對角矩陣。

值得注意的是，為了減少數(shù)據(jù)傳輸時間，將擴散系數(shù)表的生成放在了GPU端，由多線程核函數(shù)完成。同時這里的擴散系數(shù)表在查表時是具有一定隨機性的，因此采用二維紋理存儲器來存儲它，這樣在查表的過程中就可以通過紋理緩存機制實現(xiàn)比較好的帶寬利用率。二維紋理綁定的具體實現(xiàn)如下:

使用cudaMallocArray函數(shù)分配專為紋理拾取而優(yōu)化的顯存存儲空間CUDA Array空間存放數(shù)據(jù)。然后設(shè)置紋理通道參數(shù)來決定紋理拾取的工作條件，并利用cudaBindTexture函數(shù)將擴散系數(shù)表綁定到紋理存儲器。

2.5 解三對角線性方程組

本算法需要要解出的三對角線性方程組有如下形式:

其中 ai，bi，ci和 xi，i=1，…，n，為已知。而解上式的三對角線性方程組是實現(xiàn)加性算子分裂方式離散化各向異性擴散方法的關(guān)鍵。在CPU上，一般使用Thomas算法求解三對角線性方程組［8］。Thomas算法是一種高效的串行算法，但是它并不能完全發(fā)揮GPU的并行計算能力，因此需要尋找一種能夠并行化的求解方法。使用循環(huán)約化的方法來解這個三對角線性系統(tǒng)［5］，可以在GPU上并行的實現(xiàn)。

采用GPU進(jìn)行處理的循環(huán)約化的方法主要由四個部分組成:

(1)對圖像的每一行建立一個三對角矩陣。

(2)并行地解三對角矩陣。

(3)對每一列重復(fù)以上兩個步驟，即對圖像的每一列建立三對角矩陣。

(4)循環(huán)約化法使用的是高斯消去法迭代地并行消去三對角矩陣中的奇數(shù)行。其偽代碼如下:

在GPU平臺上實現(xiàn)時，設(shè)計兩個CUDA核函數(shù)分別處理行向和列向。在線程結(jié)構(gòu)方面，讓一個線程對應(yīng)與一行或一列，由于是在Fermi平臺下，每個線程塊內(nèi)設(shè)計應(yīng)至少為2個warp塊，即64個線程，對于此處測試所用的平臺，為了盡可能利用多處理器，塊內(nèi)線程數(shù)設(shè)置較小，在其它應(yīng)用場合最好根據(jù)圖像數(shù)據(jù)規(guī)模合理設(shè)置線程結(jié)構(gòu)，以盡可能利用GPU計算資源。

在存儲器使用方面，在Fermi架構(gòu)下全局存儲器已提供了L1緩存機制，因此在核函數(shù)處理設(shè)計時沒有使用共享存儲器作為中間運算的載體。而是直接在全局存儲器完成，這里為了讓L1更好的工作，使用cudaFuncSetCacheConfig函數(shù)將這兩個核函數(shù)的執(zhí)行配置設(shè)為L1優(yōu)先模式，即讓L1擁有48KB的空間，參數(shù)項為cudaFuncCachePreferL1。

2.6 圖像顯示

經(jīng)過圖像后處理的超聲圖像顯示數(shù)據(jù)存放在顯存上，要進(jìn)行顯示，既可以寫回主機端然后進(jìn)入圖形學(xué)管線進(jìn)行繪制，也可以直接利用CUDA技術(shù)與圖形學(xué)管線的互操作接口如OpenGL進(jìn)行資源共享。后者可以將顯存中存放的數(shù)據(jù)直接寫入到圖形學(xué)紋理體素空間中進(jìn)行圖像渲染顯示，這就避免內(nèi)存和顯存進(jìn)行數(shù)據(jù)傳輸?shù)臅r間消耗。具體的實現(xiàn)方式是首先初始化OpenGL設(shè)備，然后進(jìn)行紋理資源創(chuàng)建與注冊，并設(shè)定資源的讀寫屬性。在把CUDA中處理完畢的圖像數(shù)據(jù)寫入到OpenGL紋理空間，就可以進(jìn)行圖像渲染，最后釋放占有的資源。

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

實驗平臺為 2.20GHz的 AMD Althlon(tm)644200+，操作系統(tǒng)為Windows XP。GPU為 NVIDIA Geforce GTX 560 Ti，顯存為 1GB，核心頻率 1.645GHz，14個多處理器，使用4.0版本的CUDA toolkit及對應(yīng)的開發(fā)包。編程環(huán)境為Visual Studio 2010。

為了測試提出的并行實現(xiàn)算法的處理效果和運行效率，使用由數(shù)字超聲掃描器在超聲系統(tǒng)iMago C21上采集得到的人體數(shù)據(jù)作為研究工作的實驗數(shù)據(jù)。圖1(a)和(b)顯示的是由3.5MHz凸陣掃描器采集得到的人體組織的超聲原圖像，(c)和(d)是斑點噪聲抑制后CPU處理結(jié)果。(e)和(f)是斑點噪聲抑制后GPU的處理結(jié)果。通過對比知道CPU和GPU的處理結(jié)果完全保持一致，圖像像素值誤差為0，其中算法迭代次數(shù)為4，步長設(shè)置的是1.5。

圖1 采用超聲凸陣掃描器采集得到的人體組織圖像分別測試CPU和GPU的去噪效果對比圖

表1給出了所研究的斑點噪聲算法CPU和GPU處理的性能比較。程序的運行時間是圖像數(shù)據(jù)進(jìn)入顯存后進(jìn)行噪聲抑制處理的全部運行時間。從表1可以明顯看出，基于Fermi架構(gòu)GPU并行處理算法相比較于傳統(tǒng)的CPU串行處理提高了大約120倍。加速比也說明了數(shù)據(jù)計算規(guī)模越密集，GPU并行加速效果越明顯。

表1 性能比較

4 結(jié) 束 語

實驗結(jié)果顯示了基于Fermi架構(gòu)GPU的超聲斑點噪聲抑制并行實現(xiàn)方法得到的圖像去噪質(zhì)量和通用CPU處理的結(jié)果保持基本一致，而在時間性能方面達(dá)到了實時系統(tǒng)的處理要求，得到大約120多倍的加速效果。在這樣的GPU處理效率下，使用各向異性的方式進(jìn)行超聲圖像斑點噪聲抑制不僅可以取得比較好的圖像噪聲抑制效果，也可以滿足臨床超聲檢測系統(tǒng)的實時處理要求。另一方面，從加速比的分析來看，也說明了GPU并行處理平臺對于計算密集的情況加速效果更好，這與Gustafson加速比模型分析一致［9］。

［1］李治安.臨床超聲影像學(xué)［M］.北京:人民衛(wèi)生出版社，2003.

［2］J.Weichert，B.Romeny，M.A.Viergever.Efficient and reliable schemes for nonlinear diffusion filtering［J］.Image Processing，1998，7(3):398 －410.

［3］NVIDIA Corporation.CUDA ProgrammingGuide4.0［S/OL］.［2011 －05 －06］.http://www.nvdia.com.

［4］夏春蘭，石丹，劉東權(quán).基于CUDA的超聲B模式成像［J］.計算機應(yīng)用研究，2011，28(6):2011 －2015.

［5］范正娟，劉東權(quán).基于CUDA的超聲彩色血流成像［J］.計算機應(yīng)用，2011，31(3):856 －859.

［6］NVIDIA Co..CUDA API Reference Manual 4.0［M］.Santa Clara，CA，2011.

［7］Wang，Bo，Tan，Chaowei，Liu Dong C.Local Coherence based Fast Speckle Reducing Anisotropic Diffusion［C］.Proceedings of 2008 International Pre－Olympic Congress on Computer Science，Nanjing:IACSS，2008:304 －309.

［8］Krüger J.，Westermann R..Linear Algebra Operators for GPU Implementation of Numerical Algorithms［J］.ACM Transactions on Graphics.2003，22(3):908 －916.

［9］孫世新.并行算法及其應(yīng)用［M］.北京:機械工業(yè)出版社，2005.