彭 博，黃 麗

?

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

彭 博，黃 麗

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

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

超聲彈性成像；運(yùn)動(dòng)追蹤；子采樣位移估計(jì)；圖形處理器；CUDA

0 引 言

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

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

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

1 方 法

1.1 “Coupled子采樣位移估計(jì)算法”框架

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

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

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

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

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

1.3 GPU編程環(huán)境

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

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

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

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

1.5 應(yīng)變計(jì)算

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

1.6 驗(yàn)證

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

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

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

2 結(jié) 果

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

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

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

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