基于加權追蹤的醫學影像體繪制技術

閆 強，陳春曉，王 亮

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

1 引言

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

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

為了實現對CT、MRI等醫學體數據的光線追蹤多重散射體繪制，本文針對這類數據，設計了加權追蹤算法的加權系數，并通過數據預處理、細化包圍盒等方法，加快了追蹤速度。

2 加權追蹤方法

2.1 加權追蹤概述

將加權追蹤應用到體繪制中需要構造光傳輸的路徑并求解輻射傳輸方程。

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

(1)

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

(2)

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

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

(3)

對輻射傳輸方程求解時需要對光粒子的軌跡進行構造。式(3)中，光子在介質中傳輸到達t之前發生碰撞的概率滿足累積分布函數F(t)的定義

(4)

令F(t)=ξ，可得

(5)

即

(6)

對于均勻介質，μt為常量，因而

即采樣距離為

(7)

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

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

(8)

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

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

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

(9)

其中

(10)

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

2.2 加權追蹤參數設計和圖像預處理

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

(11)

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

為了減少噪聲對透射率的影響，以CT圖像為例，針對醫學體數據的背景噪聲，論文對CT數據進行了預處理：

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

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

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

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

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

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

3 加權追蹤實現流程

程序采用c++語言，利用QT 4.8.6、VTK 5.8.0函數庫實現，基于GPU的并行加權追蹤算法采用CUDASDK 8.0實現。

加權追蹤的流程圖如圖2所示，體數據載入和顯示窗口利用VTK實現，實現步驟如下：

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

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

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

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

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

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

圖2 加權追蹤流程圖

4 實驗結果與分析

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

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

表1 CT數據

圖3 加權追蹤單次和多次散射對比

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

表2 加權追蹤與Woodcock追蹤的時間對比

經數據預處理和細化包圍盒優化前后加權追蹤體繪制結果如表3所示。因為頭顱模型中存在較大的空包圍盒，因而速度的提升比例達到了15.07%。

表3 采用優化方法前后加權追蹤體繪制時間對比

5 總結與展望

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

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