GPU加速的高精度位移估計方法及超聲彈性成像應(yīng)用

彭 博，黃 麗

?

GPU加速的高精度位移估計方法及超聲彈性成像應(yīng)用

彭 博，黃 麗

( 西南石油大學(xué) 計算機(jī)科學(xué)學(xué)院，成都 610500 )

高精度的位移估計方法對提高超聲彈性成像的質(zhì)量非常重要。在本研究中，通過英偉達(dá)公司的CUDA架構(gòu)實現(xiàn)了一種新穎的可以同時提高軸向和橫向運動估計精確度的位移估計方法在GPU上的高效并行計算。對比于原始方法在C Mex編譯條件下的實現(xiàn)，GPU實現(xiàn)的方法顯示最多可實現(xiàn)76X的加速，同時保持了較高的位移估計精度。

超聲彈性成像；運動追蹤；子采樣位移估計；圖形處理器；CUDA

0 引 言

基于超聲的彈性成像技術(shù)是一種非常有前途的成像模式，它可替代觸診來評估人體組織彈性模量的變化，是傳統(tǒng)超聲成像技術(shù)的有益補充。在超聲彈性成像過程中，首先追蹤變形前后的兩幀RF信號之間的位移，然后通過對位移進(jìn)行差分操作獲得應(yīng)變信息來表征組織的彈性信息。通常在許多人體器官系統(tǒng)中，正常組織與非正常組織的彈性模量具有非常大的差異，這使超聲彈性成像技術(shù)成功的應(yīng)用在鑒別良性或惡性乳腺腫瘤[1-2]，監(jiān)視肝腫瘤的熱消融治療[3-4]等病癥的臨床治療和診斷中。

高精度的軸向與橫向位移可以提高超聲彈性成像技術(shù)的準(zhǔn)確度。然而，在普通的超聲彈性成像系統(tǒng)中，通常不具有橫向位移估計功能或獲得的橫向位移估計精度遠(yuǎn)低于軸向位移估計精度。近來Jiang提出了一種可直接采用臨床成像系統(tǒng)獲得的超聲回波數(shù)據(jù)進(jìn)行散斑追蹤并能同時獲取高精度軸向和橫向散斑追蹤結(jié)果的替代算法[5](以下簡稱Coupled子采樣位移估計算法)。該算法基于醫(yī)學(xué)超聲成像系統(tǒng)的線性系統(tǒng)理論，通過推導(dǎo)變形前后的回波信號之間的互相關(guān)函數(shù)，其演示了log壓縮的互相關(guān)函數(shù)的等高線是一系列橢圓并共享同一個中心。這些橢圓的中心正對應(yīng)著互相關(guān)函數(shù)峰值。通過確定這些橢圓的中心就能同時獲取軸向和橫向位移。對比傳統(tǒng)的運動追蹤方法，該算法在計算機(jī)模擬體模、仿真組織體模和活體數(shù)據(jù)下，都獲得了非常高的軸向與橫向子采樣位移估計精度[5]。雖然該算法具有非常突出的位移估計性能，但是在CPU的計算環(huán)境下，其計算效率不高[5]，這使得該算法的臨床應(yīng)用受到極大的限制。

當(dāng)前GPU軟硬件的發(fā)展給超聲彈性成像技術(shù)帶來了更多的機(jī)會與動力，本文的研究目的主要是利用GPU多核處理器架構(gòu)實現(xiàn)“Coupled子采樣位移估計算法”的高效并行計算。通過分析該算法易于并行處理的特性，利用CUDA架構(gòu)實現(xiàn)該算法，從而達(dá)到通過GPU并行計算滿足該算法的實時計算的需求。

1 方 法

1.1 “Coupled子采樣位移估計算法”框架

二維塊匹配方法[6]是一種常用于超聲散斑運動追蹤的方法，它的基本思想是將變形前RF信號幀分成許多宏塊，然后對每個宏塊在變形后RF信號幀中在給定特定搜索范圍內(nèi)根據(jù)一定的匹配準(zhǔn)則(如互相關(guān)函數(shù))找出與當(dāng)前塊最相似的塊，即匹配塊，匹配塊與當(dāng)前塊的相對位移即為運動矢量。當(dāng)采用互相關(guān)函數(shù)作為匹配準(zhǔn)則時，由于互相關(guān)函數(shù)的最大值位置代表了變形前RF信號幀中某一宏塊在變形后RF信號幀中對應(yīng)的位置，因此傳統(tǒng)的處理方法只關(guān)心最大值及其位置，實際上并沒有充分利用其它在搜索范圍內(nèi)已經(jīng)計算出來的互相關(guān)函數(shù)值。

通過醫(yī)學(xué)超聲系統(tǒng)的線性系統(tǒng)理論[7]，非常容易演示在上述運動追蹤過程中，變形前后的兩幀RF信號在通過互相關(guān)函數(shù)計算時，延時其log壓縮的互相關(guān)函數(shù)峰值附近可表示為

這里，1和2是相關(guān)于超聲系統(tǒng)的參數(shù)，通常為常數(shù)。D和D表示一個均勻的網(wǎng)格，可以得到一個對應(yīng)的二維互相關(guān)系數(shù)矩陣。式(1)指出，在互相關(guān)的峰值附近，基于log壓縮的互相關(guān)能被近似為一個二階的多項式表面。因此子采樣位移追蹤的目的就是發(fā)現(xiàn)對應(yīng)于互相關(guān)函數(shù)的真正最大值的位置。擬合互相關(guān)系數(shù)峰值附近的互相關(guān)系數(shù)值來估計位移有一階和二階曲線兩種方法。但“Coupled子采樣位移估計算法”不同于擬合曲線方法，它將需要擬合的互相關(guān)系數(shù)看成是一個橢圓。通過定位橢圓中心等價于尋找互相關(guān)函數(shù)的理論峰值，也就是當(dāng)和，互相關(guān)函數(shù)獲得最大值。簡而言之，這個方法首先選擇互相關(guān)峰值附近的一個互相關(guān)系數(shù)值，然后選擇作為等高線的值并擬合為一個橢圓，擬合得到的橢圓中心點在數(shù)學(xué)角度上等于未知子采樣位移。

1.2 “Coupled子采樣位移估計算法”的實現(xiàn)

“Coupled子采樣位移估計算法”包含四個步驟：1) 采用標(biāo)準(zhǔn)的二維塊匹配算法進(jìn)行初始的整數(shù)位移估計，實現(xiàn)細(xì)節(jié)可參考[6]；2) 利用步驟1得到的整數(shù)位移對變形后的RF信號幀進(jìn)行運動補償，提高計算互相關(guān)系數(shù)的相關(guān)性。在一個新的搜索區(qū)域上如[-3:0.5:3,垂直方向]′[-3:0.5:3,水平方向]計算新的互相關(guān)函數(shù)矩陣；3) 選擇合適的等高線值計算該等高線上位置的坐標(biāo)；4) 擬合所選擇的等高線的坐標(biāo)成一個橢圓[8]，橢圓方程中的中心點坐標(biāo)就是子采樣位移矢量。在第三步中，等高值的計算式[5]為

1.3 GPU編程環(huán)境

GPU編程環(huán)境下的并行計算編程可簡單分為兩個階段：1) 計算任務(wù)并行化，即是同時啟動多個GPU線程進(jìn)行計算以提高計算效率，這是算法可并行計算的必要條件。2) 訪存優(yōu)化，與CPU編程環(huán)境下不同，GPU編程環(huán)境下對于訪存優(yōu)化的要求更高。當(dāng)前在英偉達(dá)的CUDA編程架構(gòu)中，合理使用高速存儲器(寄存器，共享內(nèi)存)是提高GPU程序執(zhí)行效率的一種有效手段[10-11]。除了寄存器和共享內(nèi)存，常量存儲器和紋理映射也是一種提高訪存效率的有效內(nèi)存訪問手段。

1.4 “Coupled子采樣位移估計算法”的GPU實現(xiàn)

“Coupled子采樣位移估計算法”在塊匹配方法的框架下實現(xiàn)并行計算。在實現(xiàn)過程中，通過以下幾個策略來提高該算法的計算效率。第一，在互相關(guān)函數(shù)的計算中，軸向的互相關(guān)窗口長度和互相關(guān)窗口中的回波線數(shù)量及軸向的重疊率，這些參數(shù)決定了位移估計點數(shù)量。假定最終得到位移圖像包含rows行和cols列，rows和cols表示搜索區(qū)域的尺寸，由于這些估計點的計算過程不需要互相通信。相應(yīng)的kernel函數(shù)可以啟動(rows′cols)′(rows′cols)個GPU線程滿足算法并行執(zhí)行的需求。為保證內(nèi)存訪問的一致性，即相鄰的線程訪問相近的內(nèi)存區(qū)域，盡量滿足合并內(nèi)存事務(wù)而采用一維的線程結(jié)構(gòu)。第二，RF數(shù)據(jù)的存儲到紋理內(nèi)存，英偉達(dá)GPU的紋理內(nèi)存技術(shù)通過硬件加速數(shù)據(jù)的插值過程，避免了計算非整數(shù)位置時的插值過程。最后，一些重要的變量(如軸向與橫向搜索范圍，互相關(guān)計算核尺寸)存儲到GPU常量內(nèi)存中，以實現(xiàn)快速訪問。

在CUDA編程結(jié)構(gòu)中，kernel函數(shù)實際就是SIMD(單指令多數(shù)據(jù))的實現(xiàn)過程。而device函數(shù)可以認(rèn)為就是一個單線程函數(shù)。因為每一個估計點的計算是獨立的，因此，每一步計算過程中kernel函數(shù)的線程配置都可以啟動(rows′cols)個GPU線程滿足算法并行執(zhí)行的需求。第一步與第二步中搜索區(qū)域(rows′cols)上的互相關(guān)系數(shù)值計算可實現(xiàn)完全的并行計算，這兩步的并行計算程度較高。第三步與第四步中計算等高線上位置的坐標(biāo)的過程與擬合橢圓并求解該橢圓方程的過程均設(shè)置為device函數(shù)，并在線程配置為(rows′cols)的kernel函數(shù)中調(diào)用，這兩步的并行計算程度較低。

1.5 應(yīng)變計算

超聲彈性成像的應(yīng)變來源于對位移計算差分，軸向與橫向應(yīng)變的估計公式為和。局部軸向應(yīng)變和橫向應(yīng)變的計算通過一個低通數(shù)字差分濾波器實現(xiàn)[12]。

1.6 驗證

本研究中，計算平臺的操作系統(tǒng)為Windows 7 64-bit，計算CPU為Xeon E3-1220 V2 CPU @3.1GHz，計算GPU為NVIDIA Telsa K20c，主機(jī)內(nèi)存為8 GB。

仿真模擬數(shù)據(jù)集使用文獻(xiàn)[13]中提出的模擬壓縮前后的生物組織變化數(shù)學(xué)模型模擬獲得。仿真的散射子模型內(nèi)部包含一個硬球，模型的寬度為4 cm，深度為4 cm，硬球的半徑為0.5 cm，位于模型內(nèi)部中間位置。這個模型由200 000散射子構(gòu)成，散射子的強(qiáng)度符合高斯分布。

體模實驗采用CIRS公司(諾福克，弗吉尼亞州，美國)專用于彈性成像研究的Model 049彈性體模。實驗使用的超聲系統(tǒng)是iMago C21超聲機(jī)(Saset醫(yī)療保健公司，成都，中國)，系統(tǒng)RF信號采樣頻率設(shè)為40 MHz。實驗使用的探頭型號為SA5L38B的128陣元的線陣探頭，中心頻率為5 MHz，75%分?jǐn)?shù)階帶寬。實驗中成像對象為彈性體模中直徑為10 mm、彈性模量是63 kPa的硬包容物(Type IV)。在數(shù)據(jù)采集過程中未使用任何額外控制設(shè)備，徒手保持一恒定的速度將探頭進(jìn)行軸向壓縮/釋放。

2 結(jié) 果

表1顯示了計算不同大小位移圖像時，CPU和GPU實現(xiàn)的計算時間。對于位移估計的計算來講，單精度浮點數(shù)已經(jīng)足夠，因此本文中所有的GPU實現(xiàn)都默認(rèn)為采用單精度浮點數(shù)。其中與計算效率相關(guān)的主要參數(shù)值分別為，互相關(guān)窗口為61′11，整數(shù)位移的搜索范圍為[-5:1:5]′[-5:1:5]，子采樣互相關(guān)函數(shù)的計算范圍為[-3:0.25:3]′[-3:0.25:3]。從表1可以得出，本文所實現(xiàn)的GPU方法計算效率得到了較大的提高，GPU方法最多獲得了76X的加速。

表1 GPU實現(xiàn)與CPU實現(xiàn)的計算效率比較

Table 1 Comparison of running time for two different implementations s

The displacement of the image sizeCPU (Matlab C Mex code)GPU (Cuda C)Speed up 50′5034.4490.52066X 50′10068.8360.93074X 50′150103.3191.35476X 50′200107.7371.91970X 50′250172.1452.47869X 50′300206.4783.03268X

圖1顯示了“Couple子采樣位移估計算法”與二階拋物線擬合算法[14]在仿真組織體模下產(chǎn)生的位移及對應(yīng)的應(yīng)變圖。雖然兩種方法的軸向位移無明顯差距，但生成的軸向應(yīng)變圖可以明顯看到兩者依然具有一些差異。

圖1 仿真數(shù)據(jù)上兩種方法計算得到的位移和應(yīng)變對比