Fermi架構下各向異性擴散超聲斑點噪聲抑制

何興無

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

1 概述

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 低通濾波

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

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

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

由公式(3)可知，本算法需要進行多次迭代才能得到比較好的去噪效果，在每一次迭代中，擴散系數(shù)是不同的。如果在每次迭代過程中都要計算擴散系數(shù)，就會大大降低算法的效率［7］。于是在設計中采用生成擴散系數(shù)表，每次迭代用查表的方式得到擴散系數(shù)，這樣就減少了大量的冗余計算。因為擴散系數(shù)與圖像上某點的梯度有關，而梯度分量的取值范圍為－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ù)來決定紋理拾取的工作條件，并利用cudaBindTexture函數(shù)將擴散系數(shù)表綁定到紋理存儲器。

2.5 解三對角線性方程組

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

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

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

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

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

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

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

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

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

2.6 圖像顯示

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

3 實驗結果與分析

實驗平臺為 2.20GHz的 AMD Althlon(tm)644200+，操作系統(tǒng)為Windows XP。GPU為 NVIDIA Geforce GTX 560 Ti，顯存為 1GB，核心頻率 1.645GHz，14個多處理器，使用4.0版本的CUDA toolkit及對應的開發(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處理結果。(e)和(f)是斑點噪聲抑制后GPU的處理結果。通過對比知道CPU和GPU的處理結果完全保持一致，圖像像素值誤差為0，其中算法迭代次數(shù)為4，步長設置的是1.5。

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

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

表1 性能比較

4 結 束 語

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

［1］李治安.臨床超聲影像學［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］夏春蘭，石丹，劉東權.基于CUDA的超聲B模式成像［J］.計算機應用研究，2011，28(6):2011 －2015.

［5］范正娟，劉東權.基于CUDA的超聲彩色血流成像［J］.計算機應用，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］孫世新.并行算法及其應用［M］.北京:機械工業(yè)出版社，2005.