基于曲光線跟蹤算法的超聲成像實(shí)時(shí)模擬研究

倪 東，陳思平，汪天富

醫(yī)學(xué)超聲關(guān)鍵技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，廣東省生物醫(yī)學(xué)信息檢測(cè)與超聲成像重點(diǎn)實(shí)驗(yàn)室，深圳大學(xué)醫(yī)學(xué)院，深圳518060

醫(yī)學(xué)超聲因具有無(wú)損、無(wú)輻射、實(shí)時(shí)、便捷等特點(diǎn)，被廣泛應(yīng)用于臨床診斷、治療和介入手術(shù)引導(dǎo)中. 在醫(yī)學(xué)超聲成像研究中，為優(yōu)化超聲探頭設(shè)計(jì)，Jensen［1-3］提出一種用于任意形、不等指長(zhǎng)換能器的脈沖壓力場(chǎng)的模擬方法，在工業(yè)界和學(xué)術(shù)界均取得巨大成功. 該方法產(chǎn)生的超聲模擬圖像雖然較為逼真，但是往往需要幾個(gè)小時(shí)才能產(chǎn)生一幅圖像.

近年來(lái)，越來(lái)越多的研究者開(kāi)始關(guān)注實(shí)時(shí)超聲成像模擬，Maggie 等［4-5］提出從真實(shí)超聲體數(shù)據(jù)中提取超聲圖像紋理來(lái)構(gòu)建超聲模擬圖像. Henry［6］和倪東［7］等提出利用虛擬的超聲探頭，從已獲取的寬視野三維超聲數(shù)據(jù)中獲取二維超聲圖像. 但是，構(gòu)建寬視野三維超聲數(shù)據(jù)需要魯棒的圖像配準(zhǔn)和拼接算法，或者需要利用定位器來(lái)確定不同三維超聲數(shù)據(jù)間的位置關(guān)系. Shams 等［8］提出通過(guò)合成離線獲取的散射圖和在線獲取的反射圖來(lái)獲取超聲圖像，其中散射圖通過(guò)Field II 軟件獲取. 這種方法的一個(gè)局限性是即使在一個(gè)有20 個(gè)CPU 的計(jì)算機(jī)服務(wù)器上，產(chǎn)生散射圖也需約30 h. Song［9］提出用光線跟蹤算法模擬聲束，并用于心臟超聲圖像的模擬. 這種方法用理想化的光線模擬聲束，光線不具有面積，降低了模擬圖像的真實(shí)性. 文獻(xiàn)［4-9］的另一個(gè)局限在于不能實(shí)時(shí)調(diào)整超聲成像的重要參數(shù)，如聲波頻率和焦距. 但在實(shí)際超聲掃查過(guò)程中，選擇最優(yōu)成像參數(shù)是一個(gè)重要的技巧.

斑點(diǎn)是由于聲波的隨機(jī)散射造成的，是超聲成像的固有特性，在病理診斷中有一定作用. 但由于聲波在人體內(nèi)傳輸?shù)膹?fù)雜性，完全模擬超聲斑點(diǎn)非常耗時(shí). 近年已有研究利用不同的概率分布函數(shù)來(lái)描述超聲斑點(diǎn)模式［10-14］. Johan［10］指出斑點(diǎn)可以點(diǎn)擴(kuò)散函數(shù)來(lái)顯示.

本研究提出一種基于曲光線跟蹤算法的超聲成像模擬新方法. 利用曲光線跟蹤算法和基本聲學(xué)模型，產(chǎn)生能量反射圖和衰減圖;用多條曲光線模擬聚焦聲束，通過(guò)調(diào)整光線參數(shù)模擬聲束頻率和焦距的變化;運(yùn)用瑞利概率分布函數(shù)產(chǎn)生超聲斑點(diǎn)模式，并利用高斯型點(diǎn)擴(kuò)散函數(shù)進(jìn)行顯示. 實(shí)驗(yàn)結(jié)果表明，該算法能滿足實(shí)時(shí)模擬超聲成像的要求.

1 聲束聚焦模型

在現(xiàn)代超聲成像系統(tǒng)中，探頭孔徑由一組陣元組成. 未聚焦的換能器形成的波束較寬，導(dǎo)致橫向分辨率較低. 因此，聲束聚焦技術(shù)被廣泛應(yīng)用于現(xiàn)代超聲成像系統(tǒng)中. Christensen［15］指出聲束在近場(chǎng)具有不規(guī)則形狀，難以定義近場(chǎng)聲束模型. 但在遠(yuǎn)場(chǎng)，聲束輪廓隨著距離的變化進(jìn)行線性傳播，并偏移固定的角度，如圖1(a). 偏移角φd反比于探頭直徑D，正比于聲波波長(zhǎng)λ，定義為

文獻(xiàn)［15］還指出聚焦區(qū)域橫向分辨率與聲束寬度相關(guān)，可近似為1.22λ/D. 提高聲波頻率可增加橫向分辨率，超聲圖像質(zhì)量與聲波頻率緊密相關(guān)，而在遠(yuǎn)場(chǎng)，超聲圖像因橫向分辨率的降低而顯得模糊. 基于以上聲波傳播的基本模型，本研究提出曲光線跟蹤模型.

圖1 聲束聚焦及曲光線跟蹤模型Fig.1 Illustration of the ray tracing model

2 曲光線模型

在傳統(tǒng)光線跟蹤模型中，理想光線用單個(gè)直線表示. 由于聲束具有一定寬度及聚焦特性，我們提出用多條曲線模擬聲束. 對(duì)于探頭發(fā)出的每根掃描線如圖1(b)均用多條曲線來(lái)模擬. 定義標(biāo)量n 為曲線的標(biāo)號(hào)，T 為探頭掃描線的方向矢量，單位矢量N 為方向矢量的正交矢量，標(biāo)量s 為聲波傳播距離，曲光線的方程定義為

其中，n = 0，±1，±2，…，±L;δn(s)是用于調(diào)整曲線形狀的函數(shù)，可由式(3)計(jì)得，

這里，A 是孔徑;F 是焦距;L 是曲線在矢量T 兩側(cè)的最大個(gè)數(shù);f 是聲波頻率;β 是用于控制曲線形狀的常量.

以下分析式(2)和式(3)的物理特性. 設(shè)n= L 對(duì)應(yīng)于圖1(c)中最外面的一條曲線，當(dāng)s ＜F時(shí)，隨著傳播距離s 的增加，通過(guò)控制β 值，δn(s)趨于最小值F，曲光線的曲率可由參數(shù)β 進(jìn)行調(diào)整，與β 值成反比;當(dāng)s ≥F 時(shí)，聲束處于遠(yuǎn)場(chǎng)，將式(2)對(duì)s 求導(dǎo)數(shù)，求得曲線的斜率，即

因此，本研究提出的曲光線模型在近場(chǎng)是曲線，其寬度隨著傳播距離的增加而減小，而在遠(yuǎn)場(chǎng)是直線，直線的斜率與聲束在遠(yuǎn)場(chǎng)的偏移角的正弦值相同. 在實(shí)際應(yīng)用中，聲束偏移角φd通常較小(＜5°)［16］，則

可見(jiàn)曲光線模型與聲束聚焦模型一致. 在式(3)中，當(dāng)s = F 時(shí)，δn(s)值最小，此時(shí)它就是聲束聚焦處的寬度.

圖2 曲光線跟蹤模型示意圖Fig.2 Illustration of different effects in ray tracing process

3 反射圖模擬

圖3 為為超聲圖像模擬示意圖. 其中，圖3(a)是合成的CT 圖像，圖3(b)是對(duì)目標(biāo)光線跟蹤得到的灰度圖M，在此基礎(chǔ)上，引入基本聲學(xué)模型，計(jì)算能量反射圖R 如圖3(c)和衰減圖S 如圖3(d).

聲波在傳播過(guò)程中，由一種介質(zhì)到達(dá)另一種介質(zhì)，在兩種介質(zhì)的分界面上，由于聲阻抗的不同，聲波會(huì)發(fā)生方向和能量的變化:一部分聲波被反射回原介質(zhì)中，稱為反射波;另一部分聲波透過(guò)界面在另一種介質(zhì)中繼續(xù)傳播，稱為折射波. 這種現(xiàn)象將導(dǎo)致聲波能量的衰減. 由于超聲成像主要利用反射波檢測(cè)信號(hào)，并且折射波的模擬非常困難和耗時(shí)［4］，因此，本研究?jī)H考慮聲波的反射和由此引起的能量衰減.

聲波反射與聲阻抗密切相關(guān). 聲阻抗Z = ρc.其中，ρ 為介質(zhì)密度;c 是聲速. 聲速在人體軟組織中的傳播速度約為1 540 m/s. 設(shè)聲速為常數(shù)，聲阻抗正比于組織密度. 利用文獻(xiàn)［17］提出的CT圖像亨氏單位(Hounsfield unit)與組織密度之間存在映射關(guān)系，可由圖像灰度求取組織聲阻抗，并根據(jù)聲波能量反射系數(shù)公式(Z2-Z1)2/(Z2+Z1)2求得反射圖.

聲波衰減與聲波所經(jīng)過(guò)的組織和聲波頻率相關(guān)，文獻(xiàn)［15］指出，聲波能量I 隨傳輸距離的變化呈指數(shù)衰減，即

其中，γ 是與組織性質(zhì)相關(guān)的衰減系數(shù);f 是聲波頻率;τ 是經(jīng)驗(yàn)常數(shù)，本研究設(shè)τ = 1.

圖3 超聲圖像模擬中間結(jié)果示意圖Fig.3 Intermediate results

4 斑點(diǎn)模式模擬

斑點(diǎn)是超聲圖像由于聲波隨機(jī)散射造成的固有特性，在病理診斷中有一定作用. 但要完全從物理機(jī)制上模擬斑點(diǎn)模式是非常困難和耗時(shí)的，不適用于超聲成像的實(shí)時(shí)模擬. 瑞利分布是描述超聲斑點(diǎn)模式最常用的概率分布之一，我們采用瑞利概率分布函數(shù)產(chǎn)生超聲斑點(diǎn)模式，并利用高斯型點(diǎn)擴(kuò)散函數(shù)進(jìn)行顯示，起到平滑斑點(diǎn)噪聲的效果. 為達(dá)到實(shí)時(shí)模擬的效果，斑點(diǎn)模擬和顯示均用圖形處理器(graphics processing unit，GPU)進(jìn)行加速.

最終的超聲圖像模擬定義為

其中，S、M 和R 分別為能量衰減圖、灰度圖和反射圖;P 為高斯卷積函數(shù);A 為圖像灰度值的瑞利概率分布值;δx和δy是高斯卷積函數(shù)的參數(shù)，分別對(duì)應(yīng)橫向分辨率和縱向分辨率. 超聲圖像模擬示意圖見(jiàn)圖3(e).

5 實(shí)驗(yàn)結(jié)果分析

為驗(yàn)證算法的可行性，我們對(duì)合成圖像、解剖結(jié)構(gòu)模型以及CT 圖像進(jìn)行了超聲成像模擬，并比較了用不同的聲波頻率和聚焦位置獲得的模擬結(jié)果. 圖4 中的白色箭頭表示焦距的位置.

首先，將自行合成的體模作為超聲成像模擬的輸入，圖4(a)中(i)是分辨率體模圖像，圖4(b)中(i)是另一種體模圖像. 圖4(a)其他圖像采用不同的聲波頻率對(duì)輸入體模進(jìn)行成像:(ii)為0.5 MHz，(iii)為1 MHz，(iv)為3 MHz. 結(jié)果顯示，隨著頻率的增加，圖像的橫向分辨率得到提高，但由于高頻率導(dǎo)致聲波傳輸能量的衰減增多，灰度逐步變小. 圖4(b)采用不同的聲波頻率和焦距進(jìn)行成像:(ii)為0.5 MHz，1.5 cm;(iii)為0.5 MHz，3 cm;(iv)為1 MHz，3 cm. 從結(jié)果可見(jiàn)，隨著頻率的升高圖像的橫向分辨率得到提高，且在焦點(diǎn)位置的圖像分辨率高于其他位置的圖像.

圖4 合成體模超聲成像模擬圖像Fig.4 Effect of acoustic frequency and beam focus on phantom data

為進(jìn)一步說(shuō)明本算法的可行性，將本研究結(jié)果與超聲成像模擬的權(quán)威軟件——Field II 生成的結(jié)果進(jìn)行對(duì)比. 在文獻(xiàn)［1］ 中介紹了兩種超聲體模:腎臟體模和胎兒體模. 由于文獻(xiàn)［1］將這兩種體模描述為散射子的分布，我們首先基于解剖結(jié)構(gòu)信息將散射子分布圖轉(zhuǎn)為灰度圖，并將其作為超聲成像模擬的輸入. 圖5(a)和圖6(a)是Field II 軟件的生成結(jié)果，在圖5(b)中，由于聲束聚焦位置在上端，下面的圖像略顯模糊，而圖5(c)中，由于聲束聚焦位置在下端的圖像區(qū)域，則顯示了較好的細(xì)節(jié)，圖像也較清晰. 圖(6)顯示的是胎兒體模的超聲模擬圖，由于焦距的不同，圖6(b)和圖6(c)中的胎兒手指清晰度有較大差異.

圖5 腎臟體模超聲成像模擬圖Fig.5 Comparison of simulated kidney image

圖6 胎兒體模超聲成像模擬圖Fig.6 Comparison of simulated fetus image

最后，我們將模擬結(jié)果與真實(shí)的超聲圖像進(jìn)行了對(duì)比. 把一個(gè)實(shí)際的肝臟體模(CIRS Model 057)分別進(jìn)行CT 成像和超聲3 維成像，并將超聲三維圖像與CT 圖像進(jìn)行配準(zhǔn). 選取CT 圖像的1 個(gè)切面作為超聲成像模擬的輸入，并將模擬結(jié)果與配準(zhǔn)后相同位置的真實(shí)超聲圖像進(jìn)行對(duì)比. 圖7(a)是1幅CT 圖像，圖7(b)是用本算法產(chǎn)生的模擬結(jié)果，圖7(c)是用超聲掃描儀獲取的真實(shí)超聲圖像，兩者之間的相似性說(shuō)明本算法可行.

另外，我們還對(duì)模擬成像的速度進(jìn)行了評(píng)估.模擬軟件在臺(tái)式計(jì)算機(jī)上進(jìn)行測(cè)試(CUP 為Penium 4 dual core 2，2.8 GHz;內(nèi)存為2 Gbit;GPU 為Nvidia 8800). 通過(guò)采用GPU 對(duì)斑點(diǎn)的生成和顯示進(jìn)行加速，對(duì)于128 ×300 (128 根掃描線，每線300 個(gè)采樣點(diǎn))的超聲圖像，模擬成像速度可達(dá)約20 幀/s. 因此，本研究在超聲引導(dǎo)下的虛擬手術(shù)和計(jì)劃，超聲診斷培訓(xùn)等方面有較好應(yīng)用前景.

圖7 肝臟體模超聲成像模擬圖Fig.7 Simulated ultrasound image from CT volume

結(jié) 語(yǔ)

本研究提出一種實(shí)時(shí)模擬超聲成像的新方法，使用曲光線跟蹤模型模擬聲束，可實(shí)時(shí)調(diào)整聲波的頻率和焦點(diǎn)位置，同時(shí)利用概率分布函數(shù)產(chǎn)生超聲斑點(diǎn)，并用GPU 進(jìn)行加速顯示. 將本研究算法產(chǎn)生的結(jié)果與Field II 軟件和實(shí)際超聲掃描儀獲得的結(jié)果進(jìn)行了對(duì)比，結(jié)果的相似性說(shuō)明了本研究算法的可行性，本算法的速度也達(dá)到了實(shí)時(shí)性的要求.下一步，我們將基于聲學(xué)理論，研究斑點(diǎn)噪聲的概率分布模式、聲波散射和折射現(xiàn)象的模擬，并研究臨床病例CT 數(shù)據(jù)的超聲成像模擬，最終將研究成果用于超聲診斷培訓(xùn)、多模態(tài)圖像配準(zhǔn)以及超聲引導(dǎo)下的虛擬手術(shù)等領(lǐng)域.

/References:

［1］Jensen J A. Field:a program for simulating ultrasound systems ［J］. Medical ＆ Biological Engineering ＆ Computing，1996，34(1):351-353.

［2］Jensen J A，Gandhi D，O'Brien W D Jr. Ultrasound fields in an attenuating medium ［C］// Proceedings of the IEEE Ultrasonics symposium. Vancouver (Canada):IEEE Press，1993，2:943-946.

［3］Jensen J A，Nikolov S. Fast simulation of ultrasound images ［C］// Proceedings of the IEEE Ultrasonics symposium.Vancouver (Canada):IEEE Press，2000:1721-1724.

［4］Magee D，Zhu Y，Ratnalingam R，et al. An augmented reality simulator for ultrasound guided needle placement training ［J］. Journal of Shenzhen University Science and Engineering，2007，45(10):957-967.

［5］Zhu Y，Magee D，Ratnalingam R，et al. A training system for ultrasound-guided needle insertion procedures［C］// Medical Image Computing and Computer-Assisted Intervention.Berlin:Springer，2007:566-574.

［6］Henry D，Troccaz J，Bosson J，et al. Ultrasound imaging simulation:application to the diagnosis of deep venous thromboses of lower limbs ［C］// Medical Image Computing and Computer-Assisted Intervention. Berlin:Springer，1998:1032-1040.

［7］NI Dong，CHAN Wing-ying，QIN Jing，et al. A virtual reality simulator for ultrasound guided organ biopsy training［J］. IEEE on Computer Graphics and Application，2011，31(2):36-48.

［8］Shams R，Hartley R，Navab N. Real-time simulation of medical ultrasound from CT images ［C］// Medical Image Computing and Computer-Assisted Intervention. Berlin:Springer，2008:734-741.

［9］Song M，Haralick R，Sheehan F. Ultrasound imaging simulation and echocardiographic image synthesis ［C］//Proceedings of International Conference on Image.Vancouver:IEEE Press，2000:420-423.

［10］Johan M T. Speckle formation，analysis and processing applied to ultrasound tissue characterization ［J］. Physics for Medical Imaging Applications，2007，240(3):151-176.

［11］Palhano Xavier de Fontes，Guillermo Andrade Barroso，Pierrick Coupé，et al. Real time ultrasound image denoising ［J］. Journal of Real-time Image Processing，2011，6(1):15-22.

［12］Goodman J W. Some fundamental properties of speckle［J］. Journal of the Optical Society of America B，1976，66(11):1145-1150.

［13］CHEN Si-ping，ZENG Si-ning，WANG Tian-fu，et al.Numerical simulation on acoustic radiation force induced by ultrasonic coded excitation ［J］. Journal of Shenzhen University Science and Engineering，2011，28(2):165-171.(in Chinese)陳思平，曾斯寧，汪天富，等. 超聲編碼激勵(lì)產(chǎn)生聲輻射力的數(shù)值模擬［J］. 深圳大學(xué)學(xué)報(bào)理工版，2011，28(2):165-171.

［14］Ogier A，Hellier P. A modified total variation denoisting method in the context of 3D ultrasound images ［C］//Medical Image Computing and Computer-Assisted Intervention.Berlin:Springer，200:70-77.

［15］Christensen D A. Ultrasound Bioinstrumentation ［M］.1st ed. New York:John Wiley ＆ Sons Inc，1988:89-112.

［16］Angelsen B A J，Waag R C. Ultrasound imaging:waves，signals，and signal processing ［J］. Journal of the Acoustical Society of America，2007，121(4):1820.

［17］Wein W，Brunke S，Khamene A，et al. Automatic CTultrasound registration for diagnostic imaging and imageguided intervention ［J］. Medical Image Analysis，2008，12(5):577-585.