基于GPU加速的改進(jìn)的光線投射算法研究

摘 要： 傳統(tǒng)的醫(yī)學(xué)圖像體繪制算法大多通過CPU端實(shí)現(xiàn)，為解決傳統(tǒng)算法存在繪制時(shí)間較長(zhǎng)，交互不夠流暢，且使用平臺(tái)單一等問題，提出了一種基于圖形處理器（GPU）的醫(yī)學(xué)序列影像的實(shí)時(shí)體繪制技術(shù)。通過對(duì)不透明度值設(shè)置閾值，提前終止光線計(jì)算，并在光線遍歷體數(shù)據(jù)時(shí)，調(diào)整采樣間距，改進(jìn)了光線投射算法。實(shí)驗(yàn)結(jié)果表明，對(duì)其中一組最大數(shù)據(jù)（931張切片）的繪制時(shí)間為1.3 s左右，交互時(shí)的繪制幀數(shù)在20～40之間。在繪制時(shí)間方面，不僅比傳統(tǒng)的基于CPU的算法有大幅度提高，相比于前人的基于GPU的算法也有明顯改進(jìn)，加速比可以達(dá)到1.5左右。

關(guān)鍵詞： GPU； 光線投射算法； 體繪制； 采樣間距； 不透明度值

中圖分類號(hào)： TN911.73?34； TP391.9 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2016）21?0056?05

Study on improved ray casting algorithm based on GPU acceleration

TANG Zhenyu， Lü Xiaoqi， REN Guoyin

（School of Information Engineering， Inner Mongolia University of Science and Technology， Baotou 014010， China）

Abstract： The traditional volume rendering algorithms of medical image are mostly realized with CPU. In order to solve the problems of long rendering time， poor interaction and single platform existing in traditional algorithms， a real?time volume rendering technology of medical sequence image based on graphics process unit （GPU） is put forward. With the technology， the threshold of opacity value is set to terminate the light ray early， and the sampling interval is adjusted while the light ray traversing the volume data to improve the ray casting algorithm. The experimental results show that the volume rendering time of a group biggest data including 931 slices is about 1.3 s， and the rending frame is among 20～40 while interacting. The rendering time is faster than that of the traditional algorithm based on CPU， and the speed?up ratio can reach to about 1.5.

Keywords： GPU； ray casting algorithm； volume rendering； sampling interval； opacity value

0 引 言

近年來，計(jì)算機(jī)斷層技術(shù)（CT）、磁共振（MRI）和超聲（US）[1]等醫(yī)學(xué)影像技術(shù)快速發(fā)展。通過這些成像技術(shù)，醫(yī)生能夠獲得病人病變組織的二維切片圖像，輔助醫(yī)生進(jìn)行診斷和治療。但是，這些成像技術(shù)所提供的人體二維切片圖像只能表達(dá)人體病變組織的二維信息，病灶的大小和形狀只能憑醫(yī)生的經(jīng)驗(yàn)去估算，此方法缺乏準(zhǔn)確性和直觀性。醫(yī)學(xué)圖像三維可視化技術(shù)的發(fā)展極大地改善了這種狀況，醫(yī)學(xué)圖像三維可視化就是利用計(jì)算機(jī)將人體組織和器官的二維數(shù)字?jǐn)鄬忧衅ㄟ^插值技術(shù)將其轉(zhuǎn)化為具有直觀立體效果的圖像。

醫(yī)學(xué)圖像的三維可視化技術(shù)主要包括面繪制和體繪制[2]。與面繪制相比，體繪制重建出的圖像效果渲染質(zhì)量更高，更能清晰地展示體數(shù)據(jù)的內(nèi)部細(xì)節(jié)。體繪制中最常用的算法為光線投射算法，其繪制質(zhì)量高，繪制效果好，但計(jì)算量大，耗費(fèi)時(shí)間長(zhǎng)。在圖像質(zhì)量和繪制速度之間尋找一個(gè)最佳的解決方案，是目前三維可視化技術(shù)的研究熱點(diǎn)。隨著現(xiàn)代計(jì)算機(jī)圖形處理器（GPU）的快速發(fā)展，其浮點(diǎn)運(yùn)算速度已大大超過CPU，很多基于CPU實(shí)現(xiàn)的光線投射算法被加以改進(jìn)，以便能在GPU上運(yùn)行。文獻(xiàn)[3]提出了一種基于Pixel Shader 2.0 API的多步光線投射算法。文獻(xiàn)[4]實(shí)現(xiàn)了基于GPU的最大密度投影（MIP）繪制。盡管基于GPU的光線投射算法越來越成熟，繪制的效果也越來越好，但在繪制速度方面仍然有可提高的空間。本文基于GPU加速技術(shù)，實(shí)現(xiàn)了一種改進(jìn)的光線投射算法，該算法不僅繪制效果好，而且在繪制速度上有明顯提高。

1 體繪制算法的主要內(nèi)容

1.1 可編程GPU的繪制管線

GPU渲染流程主要由圖形繪制管線描述，即“給定視點(diǎn)、照明模式、三維物體、紋理和光源等元素，繪制一幅二維圖像”。圖形繪制管線的三個(gè)階段為：應(yīng)用程序階段、幾何階段和光柵階段。

（1） 應(yīng)用程序階段。使用高級(jí)編程語言進(jìn)行開發(fā)，主要對(duì)內(nèi)存、CPU進(jìn)行操作，比如空間八叉樹、碰撞檢測(cè)、視錐裁剪、場(chǎng)景圖建立等經(jīng)典算法都是在這個(gè)階段執(zhí)行的。在該階段的末端，紋理坐標(biāo)、頂點(diǎn)坐標(biāo)、紋理、法向量等幾何體數(shù)據(jù)通過數(shù)據(jù)總線傳送到圖形硬件。

（2） 幾何階段。主要負(fù)責(zé)投影、裁剪、屏幕映射、光照以及頂點(diǎn)坐標(biāo)的變換，該階段主要基于GPU運(yùn)行。經(jīng)過投影和變換之后的紋理坐標(biāo)、顏色以及頂點(diǎn)坐標(biāo)都是通過該階段得到。

（3） 光柵階段。根據(jù)幾何階段輸出的數(shù)據(jù)，為每一個(gè)像素進(jìn)行配色，然后繪制得到完整圖像。光柵階段是對(duì)單個(gè)像素進(jìn)行操作，顏色緩沖器用來存儲(chǔ)每個(gè)像素的信息。

1.2 光線投射算法

光線投射（Ray Casting，RC）方法是一種以空間圖像為序列的DVR算法，RC算法的基本原理是：首先從成像平面的每一個(gè)像素，沿著視線方向發(fā)射一條光線穿過體數(shù)據(jù)，然后在此過程中對(duì)圖像序列等間距的進(jìn)行重采樣并獲得顏色值，對(duì)獲得的顏色值進(jìn)行累加，當(dāng)光線穿越完圖像序列，最后得到的累加值就是圖像最終的渲染顏色[5]。其主要流程如圖1所示。

RC算法的主要技術(shù)如下：

（1） 數(shù)據(jù)分類

三維數(shù)據(jù)場(chǎng)中的數(shù)據(jù)分類就是對(duì)多種不同的物質(zhì)用不同的數(shù)值表示。在醫(yī)學(xué)圖像中，比如腎臟、骨骼、肌肉、皮膚等不同的組織器官具有不同的密度，因此其圖像的CT值也不一樣。對(duì)數(shù)據(jù)進(jìn)行分類，是醫(yī)學(xué)圖像三維重建的一個(gè)十分重要的環(huán)節(jié)。只有分類準(zhǔn)確，才能在后續(xù)的圖像處理中得到較好的效果。一般來說，一個(gè)組織器官的圖像中往往會(huì)有不同的CT值，而具有相同CT值的像素也會(huì)存在于不同組織器官的圖像中，所以有時(shí)會(huì)需要人工對(duì)數(shù)據(jù)進(jìn)行分類。

本文采用閾值分界法對(duì)醫(yī)學(xué)圖像的三維體數(shù)據(jù)進(jìn)行分類。其原理就是根據(jù)人體器官和組織，為醫(yī)學(xué)三維數(shù)據(jù)場(chǎng)中的數(shù)據(jù)點(diǎn)設(shè)置一些閾值分界點(diǎn)[dii=1，2，…，n]。假設(shè)數(shù)據(jù)點(diǎn)的像素值為[f（xi，yj，zk），]那么可以把滿足式（1）的數(shù)據(jù)點(diǎn)歸入到同一類中。

（2） 顏色賦值

一般來說，體繪制的結(jié)果顯示的圖像要具有一定的透明效果，以便更直觀地展示內(nèi)部結(jié)構(gòu)，所以要給體素賦予一定的不透明度值[a，]當(dāng)[a=1]時(shí)表示圖像完全不透明，當(dāng)[a=0]時(shí)表示圖像完全透明。如果一個(gè)體素包含了多種物質(zhì)，可以利用各種物質(zhì)所占的百分比來表示該體素的顏色值。假設(shè)第[i]種物質(zhì)所占的比例為[pi，]那么第[i]種物質(zhì)的顏色值可以表示為[Ci=（aiRi，aiGi，][aiBi，ai），]該體素的最終顏色值就可以表示為：

（3） 明暗計(jì)算

因?yàn)轶w繪制中沒有面的信息，因此體繪制中要求先算出法向量才能進(jìn)行明暗計(jì)算。在體繪制中數(shù)據(jù)點(diǎn)的法向量用梯度值代替。假設(shè)三維數(shù)據(jù)場(chǎng)中某個(gè)數(shù)據(jù)點(diǎn)的函數(shù)值為[f（xi，yj，zk），]那么就可采用如式（3）所示的中心差分法求出此數(shù)據(jù)點(diǎn)處的梯度值。

[Grad_x=f（xi+1，yj，zk）-f（xi-1，yj，zk）2ΔxGrad_y=f（xi，yj+1，zk）-f（x，yj-1，zk）2ΔyGrad_z=f（xi，yj，zk+1）-f（xi，yj，zk-1）2Δz] （3）

各數(shù)據(jù)點(diǎn)的梯度值求出后，就可以利用光照模型計(jì)算出各個(gè)數(shù)據(jù)點(diǎn)的漫反射分量，也就更能夠突顯體數(shù)據(jù)中的邊界面。

（4） 圖像合成

圖像合成是RC算法的最后一步。其主要原理為：根據(jù)像素點(diǎn)所發(fā)出的射線，此射線穿過三維體數(shù)據(jù)并進(jìn)行采樣，然后計(jì)算該采樣點(diǎn)的不透明度值和顏色值，最后通過相關(guān)公式合成該像素點(diǎn)的最終顏色值。要想形成最終的圖像，需要計(jì)算屏幕中所有的像素點(diǎn)的顏色值。體繪制中主要有兩種圖像的合成方法。

① 由前向后的圖像合成

該算法通過沿著射線把各采樣點(diǎn)的不透明度值和顏色值由前往后的合成在一起，得到最終的圖像。假設(shè)[aout]和[Cout]分別是經(jīng)過第1個(gè)體素后的不透明度值和顏色值，[ain]和[Cin]分別是進(jìn)入第[i]個(gè)體素的不透明度值和顏色值，[aout]和[Cnow]分別為第[i]個(gè)體素的不透明度值和顏色值，則有：

[Coutaout=Cinain+Cnowanow（1-ain）aout=ain+aout（1-ain）] （4）

② 由后向前的圖像合成

[Cout=Cin（1-anow）+Cnowanow] （5）

可以看出，這兩種圖像合成方法符合光學(xué)的吸收和反射模型，實(shí)質(zhì)上是一致的。其中，由前向后圖像合成過程中阻光度是逐漸增大的。如果阻光度的值接近于1，就表明該像素點(diǎn)幾乎為不透明，那么就可以把后面的像素點(diǎn)忽略掉。所以，由前向后的圖像合成方法效率更高，速度較快。本文采用由前向后的圖像合成方法。

2 改進(jìn)的光線投射算法

采用RC算法生成的圖像質(zhì)量較高，繪制效果好，但計(jì)算量大，耗費(fèi)時(shí)間長(zhǎng)。即使采用高性能的GPU進(jìn)行圖形硬件加速，有時(shí)依然達(dá)不到所要求繪制的時(shí)間。這時(shí)就需要在算法方面進(jìn)行加速優(yōu)化，以滿足實(shí)時(shí)性的要求。目前最常用的兩種方法為：無效體素剔除法（Empty Skipping and Occlusion Clipping）[3，6]和提前光線終止法（Early Ray Termination）[5?6]。本文實(shí)現(xiàn)的是提前光線終止法。

2.1 無效體素剔除法

一些體素值轉(zhuǎn)化為光學(xué)特性后，并不會(huì)顯示，也不會(huì)影響最終圖像的合成，這些體素稱為空體素。三維體數(shù)據(jù)中一般都存在一定數(shù)量的空體素。所以，通過剔除冗余體素，不讓其進(jìn)入到GPU的渲染管線中，從而達(dá)到加速繪制的效果。

2.2 提前光線終止法

通過式（5）可以看出，在像素混合迭代求解的過程中，不透明度[a]的值是逐漸增大的，當(dāng)[a]的值接近1時(shí)，表示該光線對(duì)像素值的計(jì)算結(jié)果趨近于飽和，即使繼續(xù)采樣，最終的像素值也幾乎不會(huì)受到影響。因此，可以采用提前終止光線計(jì)算的方法，這樣可以提高整體繪制的時(shí)間。可以設(shè)置一個(gè)閾值，只要[a]的值超過這個(gè)閾值，就終止光線，本文設(shè)置的閾值為0.95。同時(shí)，與自適應(yīng)加速采樣法結(jié)合，即當(dāng)光線遍歷體數(shù)據(jù)時(shí)，如果碰到無效的采樣點(diǎn)，就適當(dāng)增加采樣的間距，降低重采樣的頻率，如果碰到有效采樣點(diǎn)，就恢復(fù)重采樣頻率。通過將提前光線終止法與自適應(yīng)采樣加速法結(jié)合，可以進(jìn)一步提高整體繪制的時(shí)間。流程圖如圖2所示。

（1） 通過投射光線對(duì)體數(shù)據(jù)進(jìn)行重采樣，獲得體素的值并進(jìn)行數(shù)據(jù)分類。

（2） 判斷體素值是否有效。如果無效，則增大采樣間距，并轉(zhuǎn)到（5）；如果有效，則采用原采樣間距，并進(jìn)行下一步。

（3） 對(duì)采樣點(diǎn)進(jìn)行顏色賦值和明暗計(jì)算，并進(jìn)行圖像合成。

（4） 判斷不透明度值[a]是否大于閾值，如果沒有大于閾值，則繼續(xù)重采樣；如果大于閾值，則轉(zhuǎn)到（6）。

（5） 判斷光線是否穿出包圍體，如果沒有穿出，則繼續(xù)重采樣；如果穿出，則進(jìn)行下一步。

（6） 根據(jù)累加的顏色值，合成最終圖像。

3 GPU加速光線投射算法

RC算法在GPU上實(shí)現(xiàn)的關(guān)鍵是能夠利用GPU完成對(duì)體數(shù)據(jù)的遍歷、采樣和計(jì)算，并將得到的采樣值進(jìn)行混合，以產(chǎn)生最終的重建結(jié)果[7?9]。此部分通過Cg語言和OpenGL編寫的頂點(diǎn)程序和片段程序?qū)崿F(xiàn)。Cg語言是為GPU編程設(shè)計(jì)的高級(jí)著色語言，而OpenGL擴(kuò)展則提供了可操作GPU的編程接口。頂點(diǎn)程序主要負(fù)責(zé)射線發(fā)出的方向和設(shè)置頂點(diǎn)的位置以及變換頂點(diǎn)坐標(biāo)、頂點(diǎn)顏色和紋理坐標(biāo)。而片段程序主要負(fù)責(zé)核心算法，即改進(jìn)的光線投射算法。通過可編程GPU編寫著色器，可以非常靈活方便地實(shí)現(xiàn)各種繪制效果。本文通過結(jié)合GPU加速技術(shù)和改進(jìn)的光線投射算法實(shí)現(xiàn)了最大密度投影繪制（Maximum Intensity Projection，MIP）[10]以及合成體繪制（包括骨骼、皮膚、肌肉等）[11]。不同的繪制方法分別對(duì)應(yīng)不同的組織結(jié)構(gòu)。

3.1 最大密度投影繪制

最大密度投影就是通過計(jì)算每條射線所遇到的最大像素密度確定最終圖像的顏色值和不透明度值。圖像效果可以反映組織密度的差異，對(duì)比度高，常應(yīng)用于密度相對(duì)較高的組織結(jié)構(gòu)。用基于GPU的光線投射算法實(shí)現(xiàn)最大密度投影算法的原理為：計(jì)算發(fā)射光線與體數(shù)據(jù)相交時(shí)的最大密度值，并把值賦予該像素。如圖3所示。

MIP繪制的效果如圖4所示。

3.2 合成體繪制

合成體繪制是最常用的一種光線投射方法。其原理為：通過頂點(diǎn)程序把體數(shù)據(jù)值轉(zhuǎn)換成顏色值和不透明度值等光學(xué)屬性，并將這些屬性值合成到屏幕上的像素中形成最終圖像。如果改變顏色值和不透明度值的參數(shù)，就能實(shí)現(xiàn)不同的繪制效果，比如骨骼、肌肉等。

繪制顯示效果如圖5所示。

4 繪制速度對(duì)比

實(shí)驗(yàn)數(shù)據(jù)采用Philips公司Brilliance iCT設(shè)備產(chǎn)生的人體數(shù)據(jù)。計(jì)算機(jī)配置為Intel[?] CoreTM i7?3770 CPU @3.40 GHz；顯卡為AMD HD 7470，1 GB顯存；內(nèi)存為3.49 GB；Windows 7 32位操作系統(tǒng)；開發(fā)平臺(tái)為Eclipse Luna Service Release1（4.4.1）。繪制時(shí)間通過測(cè)量十次取平均值得到，繪制幀數(shù)采用第三方軟件Fraps測(cè)量得到。

耳部數(shù)據(jù)集的大小為828 MB，頭部數(shù)據(jù)集的大小為477 MB，腳踝數(shù)據(jù)集的大小為75 MB。通過表1可以看出，繪制時(shí)間上，基于GPU的算法明顯快于基于CPU的算法，而本文實(shí)現(xiàn)的基于GPU改進(jìn)的光線投射算法要優(yōu)于傳統(tǒng)的GPU算法，最大加速比達(dá)到1.79。同時(shí)可以看出，數(shù)據(jù)集越大，基于GPU的加速效果越明顯。

5 結(jié) 語

隨著現(xiàn)代醫(yī)學(xué)成像設(shè)備的發(fā)展，其所生成的醫(yī)學(xué)圖像數(shù)據(jù)越來越龐大，傳統(tǒng)的基于CPU的體繪制算法已不能滿足醫(yī)學(xué)臨床需求，基于GPU的醫(yī)學(xué)圖像處理算法逐漸成為目前研究的熱門領(lǐng)域。本文在Eclipse平臺(tái)下，基于GPU加速，結(jié)合改進(jìn)的光線投射算法，實(shí)現(xiàn)了最大密度投影繪制以及合成體繪制。Eclipse是一個(gè)兼容Windows系統(tǒng)和Linux系統(tǒng)的平臺(tái)，基于Eclipse平臺(tái)開發(fā)的程序移植也比較方便。通過結(jié)論可以看出，繪制效果已能滿足醫(yī)學(xué)臨床實(shí)踐的需求，且繪制速度上比傳統(tǒng)的GPU算法有明顯的提升。下一步的研究工作主要為完善程序的功能，實(shí)現(xiàn)更多的優(yōu)化算法，并對(duì)各種算法進(jìn)行優(yōu)缺點(diǎn)的比較。

參考文獻(xiàn)

[1] KUTTER O， SHAMS R， NAVAB N. Visualization and GPU?accelerated simulation of medical ultrasound from CT images [J]. Computer methods and programs in biomedicine， 2009， 94（3）： 250?266.

[2] 常旖旎，魯雯，聶生東.醫(yī)學(xué)圖像三維可視化技術(shù)及其應(yīng)用[J].中國(guó)醫(yī)學(xué)物理學(xué)雜志，2012，29（2）：3254?3258.

[3] LI W， MUELLER K， KAUFMAN A. Empty space skipping and occlusion clipping for texture?based volume rendering [C]// Proceedings of 2003 the 14th IEEE Visualization Conference. Washington， D. C.： IEEE， 2003： 1?8.

[4] KRUGER J， WESTERMANN R. Acceleration techniques for GPU?based volume rendering [C]// Proceedings of 2003 the 14th IEEE Visualization Conference. Washington， D. C.： IEEE， 2003： 287?292.

[5] LEVOY M. Efficient ray tracing of volume data [J]. ACM tran?sactions on graphics， 1990， 9（3）： 245?261.

[6] HADWIGER M， LJUNG P， SALAMA C R， et al. Advanced illumination techniques for GPU volume raycasting [C]// Procee?dings of 2008 ACM Siggraph Asia Courses. [S.l.]： ACM， 2008： 163?166.

[7] 胡劍，王明泉，侯慧琳，等.使用 GPU 編程的工業(yè)CT斷層圖像三維可視化技術(shù)[J].應(yīng)用科技，2011，38（5）：61?65.

[8] 康健超，康寶生，馮筠，等.一種改進(jìn)的基于 GPU 編程的光線投射算法[J].計(jì)算機(jī)工程與應(yīng)用，2012，48（1）：199?201.

[9] 林碧英，陳達(dá)峰.基于GPU的三維數(shù)據(jù)場(chǎng)等值面繪制[J].中國(guó)電力教育，2008（z3）：402?404.

[10] 諸葛斌，馮煥清，周荷琴.基于有序體數(shù)據(jù)的最大密度投影算法[J].中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2003，33（4）：473?479.

[11] 呂曉琪，張傳亭，侯賀，等.基于圖形處理器加速光線投射算法的多功能體繪制技術(shù)[J].計(jì)算機(jī)應(yīng)用，2014，34（1）：135?138.

[12] WANG X H， GOOD W F. Real?time stereographic rendering and display of medical images with programmable GPUs [J]. Computerized medical imaging and graphics， 2008， 32（2）： 118?123.