閆 強，陳春曉，王 亮

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

1 引言

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

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

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

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點沿ω方向d長度點的輻射，d是從x到碰撞點的距離，ua是吸收系數(shù)，Le(xt，ω)是xt點沿ω方向的放射輻射，us是散射系數(shù)，Ls(xt，ω)是y點沿ω方向的散射輻射，T(t)指從點x沿ω方向t長度點的透射率函數(shù)，它由Beer-Lambert定律推導(dǎo)得出

(2)

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

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

(3)

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

(4)

令F(t)=ξ，可得

(5)

即

(6)

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

即采樣距離為

(7)

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

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

(8)

圖1 加入虛擬粒子后光線行進路徑

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

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

(9)

其中

(10)

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

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

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

(11)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實驗結(jié)果與分析

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

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

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

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

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

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

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

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

5 總結(jié)與展望

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

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