基于加權(quán)追蹤的醫(yī)學(xué)影像體繪制技術(shù)

閆 強(qiáng)，陳春曉，王 亮

(南京航空航天大學(xué)生物醫(yī)學(xué)工程系，江蘇 南京211106)

1 引言

計(jì)算機(jī)斷層掃描(Computed Tomography，CT)和磁共振成像(Magnetic Resonance Image，MRI)等生成的三維數(shù)據(jù)場(chǎng)可以通過體繪制技術(shù)在屏幕上顯示為二維圖像。體繪制技術(shù)能夠提供豐富的細(xì)節(jié)特征和三維的整體結(jié)構(gòu)，在臨床診斷和科學(xué)研究中具有廣泛的應(yīng)用。體繪制算法主要有光線投射算法、錯(cuò)切-變形算法、頻域繪制算法、拋雪球算法、光線追蹤算法等。目前，光線投射算法因需要的計(jì)算量較小且能夠提供較好的重建效果，得到普遍應(yīng)用。隨著計(jì)算機(jī)硬件技術(shù)的發(fā)展，實(shí)現(xiàn)可交互的光線追蹤體繪制成為了可能。光線追蹤算法可以模擬光在介質(zhì)中的真實(shí)傳播路徑，生成的圖像能反映陰影、景深等信息，視覺效果更加真實(shí)，有助于醫(yī)生更快速、更準(zhǔn)確的診斷[1]。

光線追蹤體繪制是一種基于物理的渲染方法，利用輻射傳輸方程求解光在介質(zhì)中傳播的能量變化，可以模擬現(xiàn)實(shí)中的光照效果。該方程由Chandrasekhar等人[2]提出，將輻射能變化與碰撞系數(shù)聯(lián)系起來。求解輻射傳輸方程的關(guān)鍵在于對(duì)自由路徑組成的光粒子軌跡的隨機(jī)構(gòu)造。對(duì)于簡(jiǎn)單的體數(shù)據(jù)，如體數(shù)據(jù)的消光系數(shù)設(shè)定為恒定值或?yàn)槎囗?xiàng)式形式變化時(shí)，可使用Ulam等人[3]提出的封閉式追蹤算法，此方法核心是對(duì)透射率累積分布函數(shù)進(jìn)行逆變換，進(jìn)而對(duì)光子自由路徑進(jìn)行采樣。雖然封閉式追蹤算法是無偏的，但其只能應(yīng)用于簡(jiǎn)單體數(shù)據(jù)，對(duì)于稍復(fù)雜的體數(shù)據(jù)可使用由Sutton等人[4]提出的規(guī)則追蹤算法，它將體數(shù)據(jù)看作為均勻的分塊，然后通過沿射線方向累積每一個(gè)塊內(nèi)的光學(xué)厚度來構(gòu)造自由路徑。規(guī)則追蹤算法中，當(dāng)塊數(shù)量變多時(shí)需要的計(jì)算量會(huì)迅速增加，追蹤速度變慢。降低追蹤所需計(jì)算量的一個(gè)直接方法是忽略邊界，Perlin等人[5]提出的光線步進(jìn)算法采用固定步長(zhǎng)沿光線方向行進(jìn)，可加快計(jì)算速度，但計(jì)算結(jié)果存在偏差。另一種方法是由Rabb等人[6]提出的，他們將核物理科學(xué)領(lǐng)域的Woodcock追蹤方法引入圖形學(xué)領(lǐng)域，該方法嚴(yán)格遵守物理過程，可以對(duì)非均勻體數(shù)據(jù)自由路徑進(jìn)行無偏采樣。Thomas等人[7]在2012年基于Woodcock追蹤構(gòu)建了一個(gè)可實(shí)時(shí)交互的單次散射蒙特卡羅光線追蹤體繪制系統(tǒng)，該系統(tǒng)模擬了真實(shí)鏡頭和光圈，實(shí)現(xiàn)了基于物理的體陰影，但它對(duì)于模型的大小有限制且無法進(jìn)行多重散射。Woodcock追蹤算法的局限是在稀疏的體數(shù)據(jù)中會(huì)產(chǎn)生較多的零碰撞，從而造成計(jì)算量的增加，針對(duì)這一局限性，Yue等人[8]使用kd樹，Szirmay-Kalos等人[8]使用網(wǎng)格，對(duì)體數(shù)據(jù)進(jìn)行分割，降低了 Woodcock追蹤算法中的零碰撞系數(shù)，減少了零碰撞次數(shù)，相應(yīng)減少了計(jì)算量。另一種減少計(jì)算量的方法是加權(quán)追蹤，可以應(yīng)用此方法來降低組合系數(shù)，減少零碰撞，Galtier等人[10]推導(dǎo)出加權(quán)追蹤的一般形式，并證明了加權(quán)追蹤的無偏性。Kutz等人[10]在Galtier的工作基礎(chǔ)上提出了多波長(zhǎng)追蹤算法和分解追蹤算法，該算法可以合并頻譜數(shù)據(jù)，消除對(duì)組合系數(shù)的限制，尤其針對(duì)稀疏體數(shù)據(jù)可大量減少對(duì)系數(shù)的查找次數(shù)，從而加快渲染速度。減少需要渲染的數(shù)據(jù)量同樣可以加快渲染速度，Rabbani[12]通過渲染像素的子集及使用基于GPU的共軛梯度求解器重建完整圖像，成功優(yōu)化了體渲染器的性能。針對(duì)大規(guī)模體數(shù)據(jù)，Markus[13]等人提出了稀疏跳躍算法，可以跳過體數(shù)據(jù)中的空白部分，對(duì)于大型、復(fù)雜體繪制有顯著的改進(jìn)。

為了實(shí)現(xiàn)對(duì)CT、MRI等醫(yī)學(xué)體數(shù)據(jù)的光線追蹤多重散射體繪制，本文針對(duì)這類數(shù)據(jù)，設(shè)計(jì)了加權(quán)追蹤算法的加權(quán)系數(shù)，并通過數(shù)據(jù)預(yù)處理、細(xì)化包圍盒等方法，加快了追蹤速度。

2 加權(quán)追蹤方法

2.1 加權(quán)追蹤概述

將加權(quán)追蹤應(yīng)用到體繪制中需要構(gòu)造光傳輸?shù)穆窂讲⑶蠼廨椛鋫鬏敺匠獭?/p>

μs(xt)Ls(xt，ω)+Ld(xd，ω)]dt

(1)

其中L(x，ω)是x點(diǎn)沿ω方向d長(zhǎng)度點(diǎn)的輻射，d是從x到碰撞點(diǎn)的距離，ua是吸收系數(shù)，Le(xt，ω)是xt點(diǎn)沿ω方向的放射輻射，us是散射系數(shù)，Ls(xt，ω)是y點(diǎn)沿ω方向的散射輻射，T(t)指從點(diǎn)x沿ω方向t長(zhǎng)度點(diǎn)的透射率函數(shù)，它由Beer-Lambert定律推導(dǎo)得出

(2)

其中μt指消光系數(shù)，是μa與μs的和。透射率T(t)是衡量x到xt之間因介質(zhì)吸收和向外散射的凈減少因子。式(1)中積分符號(hào)內(nèi)的第三項(xiàng)為邊界項(xiàng)，它不依賴于t，因而可以對(duì)單獨(dú)積分，從而得到簡(jiǎn)化后的輻射傳輸方程

μs(xt)Ls(xt，ω)]dt+T(d)Ld(xd，ω)

(3)

對(duì)輻射傳輸方程求解時(shí)需要對(duì)光粒子的軌跡進(jìn)行構(gòu)造。式(3)中，光子在介質(zhì)中傳輸?shù)竭_(dá)t之前發(fā)生碰撞的概率滿足累積分布函數(shù)F(t)的定義

(4)

令F(t)=ξ，可得

(5)

即

(6)

對(duì)于均勻介質(zhì)，μt為常量，因而

即采樣距離為

(7)

然而，體數(shù)據(jù)一般都是非均勻的。在處理非均勻介質(zhì)時(shí)，可以在介質(zhì)中插入虛擬粒子，使得插入虛擬粒子之后的非均勻介質(zhì)變成均勻的，這樣就可以利用公式(7)所求的距離進(jìn)行采樣。插入虛擬粒子后會(huì)發(fā)生零碰撞，但光線行進(jìn)的路徑不變，如圖1所示。于是，輻射傳輸方程可表示為

μn(xt)L(xt，ω)]

(8)

圖1 加入虛擬粒子后光線行進(jìn)路徑

由于直接求解式(8)的計(jì)算量很大，實(shí)際求解時(shí)會(huì)結(jié)合蒙特卡羅采樣算法，得到加權(quán)追蹤體繪制方程的一般形式：

〈μs(xt)Ls(xt，ω)〉Ps+〈μn(xt)L(xt，ω)〉Pn]

(9)

其中

(10)

式(9)中Pa、Ps、Pn表示交互概率，只要不違反概率的定義，可以任意設(shè)計(jì)。

2.2 加權(quán)追蹤參數(shù)設(shè)計(jì)和圖像預(yù)處理

在醫(yī)學(xué)體數(shù)據(jù)中，一般不存在發(fā)射光，因而將發(fā)射項(xiàng)系數(shù)置為0，將式(9)中的發(fā)射項(xiàng)消除。散射項(xiàng)和零碰撞項(xiàng)各自發(fā)生的概率由消光系數(shù)所占的比率決定。于是，加權(quán)追蹤的系數(shù)可以設(shè)置為

(11)

其中μn加絕對(duì)值是因?yàn)棣蘮可以為負(fù)值，這樣可以減小μn與μs組成的總消光系數(shù)，從而減少零碰撞發(fā)生的次數(shù)。雖然物理上負(fù)的消光系數(shù)難以解釋，但在數(shù)學(xué)上可以證明，取負(fù)的消光系數(shù)后式(9)的計(jì)算結(jié)果仍然是無偏的[10]。

為了減少噪聲對(duì)透射率的影響，以CT圖像為例，針對(duì)醫(yī)學(xué)體數(shù)據(jù)的背景噪聲，論文對(duì)CT數(shù)據(jù)進(jìn)行了預(yù)處理：

1)用3×3的核對(duì)CT數(shù)據(jù)進(jìn)行中值濾波；

2)利用幀差法，對(duì)CT數(shù)據(jù)間隔5張的兩張圖像求差分，獲取差分圖像直方圖的第一個(gè)峰值點(diǎn)(x，y)，其灰度值為d(x，y)；

3)計(jì)算差分圖像中灰度值為d(x，y)的所有點(diǎn)在原始CT圖像對(duì)應(yīng)位置處灰度的平均值，并以此作為提取背景的參考值；

4)計(jì)算體數(shù)據(jù)的直方圖，在背景參考值附近尋找峰值點(diǎn)，峰值點(diǎn)后的第一個(gè)谷值即為提取背景的閾值；

5)將CT圖像中所有灰度小于背景閾值的置為0。

在體繪制中，遍歷整個(gè)數(shù)據(jù)進(jìn)行處理是導(dǎo)致繪制速度慢的主要原因，因此，需要對(duì)背景去噪處理后的CT數(shù)據(jù)創(chuàng)建包圍盒。包圍盒是指包含所有體數(shù)據(jù)的最小長(zhǎng)方體。在醫(yī)學(xué)體繪制中，很多區(qū)域?qū)ψ罱K繪制結(jié)果不會(huì)產(chǎn)生影響，通常稱為空域。針對(duì)這部分空域，論文借鑒層次包圍盒的思想，在體數(shù)據(jù)整體的大包圍盒上從四周開始進(jìn)行擴(kuò)展，首先在x，y方向同步擴(kuò)展，在接觸到非零體數(shù)據(jù)后分別在x，y方向擴(kuò)展，找出面積最大包圍盒，最后刪去包圍盒的重疊部分。將獲取的小包圍盒稱為空包圍盒，除去空包圍盒的部分稱為實(shí)包圍盒。實(shí)包圍盒便是實(shí)際需要繪制的區(qū)域。

3 加權(quán)追蹤實(shí)現(xiàn)流程

程序采用c++語言，利用QT 4.8.6、VTK 5.8.0函數(shù)庫(kù)實(shí)現(xiàn)，基于GPU的并行加權(quán)追蹤算法采用CUDASDK 8.0實(shí)現(xiàn)。

加權(quán)追蹤的流程圖如圖2所示，體數(shù)據(jù)載入和顯示窗口利用VTK實(shí)現(xiàn)，實(shí)現(xiàn)步驟如下：

1)生成光線：由攝像機(jī)位置與屏幕像素點(diǎn)的連線確定光線方向，光線沿此射線行進(jìn)；

2)判斷光線是否與主包圍盒相交，將不相交的光線拋棄；

3)判斷是否光線是否與實(shí)包圍盒相交，若相交則記錄下交點(diǎn)，若不相交，則需要求取光線穿過空包圍盒的出射點(diǎn)，再判斷它是否與實(shí)包圍盒相交，如果與實(shí)包圍盒相交，則記錄交點(diǎn)；

4)路徑生成：利用加權(quán)的消光系數(shù)生成路徑長(zhǎng)度s，從交點(diǎn)處沿射線方向行進(jìn)s，然后采用俄羅斯輪盤賭方法判斷光線是否與該點(diǎn)發(fā)生碰撞，如果發(fā)生碰撞，利用式(9)和(11)計(jì)算出輻射量，并加入該光線的累積輻射當(dāng)中，之后開始下一個(gè)碰撞檢測(cè)，否則光線繼續(xù)沿此方向行進(jìn)路徑長(zhǎng)度s；

5)設(shè)計(jì)累計(jì)5次碰撞后停止光線行進(jìn)。碰撞次數(shù)的增多，光線對(duì)輻射度的貢獻(xiàn)越來越小，而且會(huì)導(dǎo)致計(jì)算時(shí)間大量增加；

6)在一輪光線追蹤終止之后，將生成的圖像進(jìn)行濾波、顏色映射及顯示。

圖2 加權(quán)追蹤流程圖

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

本文實(shí)現(xiàn)了基于加權(quán)追蹤體繪制技術(shù)，利用預(yù)處理算法和包圍盒細(xì)分算法對(duì)加權(quán)追蹤算法進(jìn)行優(yōu)化。論文對(duì)比了多次散射與單次散射的繪制結(jié)果，加權(quán)追蹤與Woodcock追蹤算法多次散射的繪制結(jié)果，以及加權(quán)追蹤優(yōu)化前后加的繪制結(jié)果。實(shí)驗(yàn)平臺(tái)：I7 9700 CPU，Nvidia 1660Ti GPU，內(nèi)存16GB。實(shí)驗(yàn)所使用的CT數(shù)據(jù)如表1所示，分別采用了頭顱模型、軀干模型、小鼠模型，體數(shù)據(jù)規(guī)模從6Mb到133Mb不等。

設(shè)定同樣的相機(jī)位置、傳遞函數(shù)、光源，迭代20次后。多次散射和單次散射加權(quán)追蹤結(jié)果如圖3所示。結(jié)果顯示，在同樣迭代次數(shù)下，多重散射得到的圖像更加平滑，相比之下單次散射的結(jié)果噪聲更多。

表1 CT數(shù)據(jù)

圖3 加權(quán)追蹤單次和多次散射對(duì)比

加權(quán)追蹤與Woodcock追蹤多次散射體繪制的結(jié)果如表2所示。兩種方法生成圖像的平均結(jié)構(gòu)相似性達(dá)到了0.9783，在視覺感知上幾乎不存在差別，但本文實(shí)現(xiàn)的加權(quán)追蹤平均耗費(fèi)時(shí)間比Woodcock追蹤減少了12.063%。

表2 加權(quán)追蹤與Woodcock追蹤的時(shí)間對(duì)比

經(jīng)數(shù)據(jù)預(yù)處理和細(xì)化包圍盒優(yōu)化前后加權(quán)追蹤體繪制結(jié)果如表3所示。因?yàn)轭^顱模型中存在較大的空包圍盒，因而速度的提升比例達(dá)到了15.07%。

表3 采用優(yōu)化方法前后加權(quán)追蹤體繪制時(shí)間對(duì)比

5 總結(jié)與展望

本文根據(jù)醫(yī)學(xué)體數(shù)據(jù)的特性設(shè)計(jì)了加權(quán)追蹤系數(shù)，且提出了基于幀差法進(jìn)行數(shù)據(jù)預(yù)處理和細(xì)化包圍盒的多重散射加權(quán)追蹤體繪制方法。實(shí)驗(yàn)結(jié)果表明，多重散射加權(quán)追蹤體繪制對(duì)體數(shù)據(jù)的內(nèi)部細(xì)節(jié)和外部特點(diǎn)均有很好的展示，而且多重散射比單次散射更能得到光滑、細(xì)致逼真的圖像。對(duì)比Woodcock追蹤方法，加權(quán)追蹤在繪制速度上有很大提升。

針對(duì)目前廣泛使用的光線投射方法而言，光線追蹤方法的繪制結(jié)果更加逼真。本文基于加權(quán)追蹤實(shí)現(xiàn)了光線追蹤的體繪制，但繪制速度仍不能滿足實(shí)時(shí)要求。因而，如何加快光線追蹤體繪制算法的速度，仍需要進(jìn)一步探索。