基于圖形處理器的射線追蹤數(shù)字重建影像方法

劉 操，汪 俊，吳章文，勾成俊，侯 氫

（四川大學(xué)原子核科學(xué)技術(shù)研究所輻射物理及技術(shù)教育部重點實驗室，四川 成都 610064）

0 引 言

數(shù)字重建影像（digitally reconstructed radiography，DRR）是腫瘤放射治療計劃系統(tǒng)（treatment planning system，TPS）里的重要組成部分。DRR圖像可用于指導(dǎo)放射腫瘤醫(yī)生和物理師確定腫瘤靶區(qū)位置和設(shè)定射線的方向以及射野形狀[1]，還用于校正擺位誤差和驗證射野輻射注量或劑量分布。這些都是交互性的試錯過程，然而通過計算機斷層掃描圖像（computed tomography，CT）重建DRR圖像又是一個計算密集過程。因此，發(fā)展能快速獲得DRR圖像的方法具有實際應(yīng)用價值。

加速DRR重建速度，可通過簡化DRR重建的計算模型或運用計算機領(lǐng)域的新技術(shù)。基于可編程圖形處理器（graphics processing unit，GPU）的并行計算，在圖形圖像及科學(xué)計算等領(lǐng)域得到廣泛應(yīng)用[2-5]。與傳統(tǒng)的并行計算構(gòu)架如集群計算相比，基于GPU的并行運算可在個人電腦上實現(xiàn)，且在相同計算效率情況下花費更少。

在各種DRR重建算法中，射線追蹤算法有更好的計算精度[6-7]，然而也是計算量最大的算法。本文提出了基于GPU并行的射線追蹤算法用于DRR影像重建。

1 算法和GPU實現(xiàn)

就算法本身而言，基于中央處理器（central processing unit，CPU）和GPU的射線追蹤數(shù)字影像重建過程具有相似的處理流程，如圖1所示。然而，基于CPU的實現(xiàn)不需要CPU主內(nèi)存與GPU顯存之間的數(shù)據(jù)拷貝過程。總的來說，處理流程分為4個部分。第1部分，載入體模的CT數(shù)據(jù)，然后轉(zhuǎn)化為體元的X線吸收系數(shù)u，轉(zhuǎn)換公式[7-8]為

式中：uw——水的線性吸收系數(shù)；

HUi——第i個體元的CT值。

為便于討論，將ui表示的體模稱為初始體模。由于ui是體模在人體坐標系下體元的吸收系數(shù)，而射線追蹤將在射束坐標系下進行；因此，第2部分就是在射束坐標系下創(chuàng)建虛擬體模，并計算出體元的吸收系數(shù)；獲得虛擬體模后，第3部分則計算射線沿投影平面上的像素穿過虛擬體模到達X線源過程中的射線衰減；第4部分在設(shè)備上顯示DRR圖像。

圖1 基于GPU實現(xiàn)射線追蹤算法重建DRR圖像的流程圖

第2部分和第3部分是DRR重建中計算量最大的部分，接下來將詳細說明這兩部分在GPU上的實現(xiàn)過程。

1.1 計算吸收系數(shù)?的GPU實現(xiàn)

當射束以一定角度傾斜照射體模時，創(chuàng)建一個表面與射束方向垂直的虛擬體模，有利于進行射線追蹤，如圖2所示。射束坐標系和人體坐標系在圖2中也做了標示。虛擬體模中處于位置的體元α，在人體坐標系下的位置ra可以表示為

式中：θG——機頭的轉(zhuǎn)角；

θT——治療床的轉(zhuǎn)角；

θC——準直器的轉(zhuǎn)角；

T（θG，θT，θC）——射束坐標系與人體坐標系間的3×3轉(zhuǎn)換矩陣。

如果ra在初始體模中處在以其相鄰8個體元中心作為頂點的立方體中，吸收系數(shù)可通過對初始體模中周圍體元的吸收系數(shù)u進行線性插值獲得。

圖2 射線追蹤示意圖

通常，DRR圖像重建過程中處理的體元數(shù)目很多，虛擬體模具有更多的體元。然而在CPU上實現(xiàn)上述插值運算，要對虛擬體模的每個體元進行循環(huán)。相比而言，GPU實現(xiàn)中的循環(huán)被并行的線程所替代。在GPU實現(xiàn)過程中，具有與體元相同數(shù)目的線程通過GPU核函數(shù)進行分配和運行[5]，每個線程完成對1個體元的插值運算。初始體模中體元的吸收系數(shù)u拷貝到GPU后，與紋理內(nèi)存綁定。由于紋理單元具有硬件插值功能，文中對直接運用紋理單元做線性插值，稱為硬件插值；對提取吸收系數(shù)u再進行線性插值，稱為軟件插值。兩種插值方式除了聲明不同外，具有相同的執(zhí)行流程。在每個線程中，虛擬體模的體元在人體坐標系中的位置，通過式（2）計算獲得。為縮短讀取數(shù)據(jù)的時間，轉(zhuǎn)換矩陣 T（θG，θT，θC）存放在GPU的常量內(nèi)存中。

1.2 射線追蹤過程的GPU實現(xiàn)

DRR圖像是通過計算射線穿過虛擬體模后，透射到投影平面上的射線強度獲得的。投影平面被離散化為像素單元，每個像素與1條射線相關(guān)聯(lián)。像素i處的強度Fi計算式為

式中：I0——X線源的強度；

li——像素i到射線源的積分路徑；

liα——第i條射線穿過體元α的幾何長度。

像素值即灰度值的大小，定義為I0-Fi。對于DRR圖像重建中的每個像素點，相關(guān)射線都從像素的中心位置向射線源進行追蹤，追蹤的步長與虛擬體模中體元的大小相同。與計算吸收系數(shù)u?的插值運算類似，在GPU計算中，線程與像素點對應(yīng)，完成式（3）表示的數(shù)值積分。

2 結(jié)果及分析

采用CUDA C編程，實現(xiàn)了基于GPU的DRR射線追蹤重建算法。采用頭頸部體模（模型1）和骨盆體模（模型2），對GPU實現(xiàn)進行了性能測試。測試模型的機頭轉(zhuǎn)角分別選取0°，90°和-45°，準直器和治療床的轉(zhuǎn)角都為0°。測試平臺是擁有IntelCore2 Q8400的CPU和NVIDIA GeForce G210顯卡的個人電腦，程序的編譯和運行環(huán)境為Windows操作系統(tǒng)。

2.1 計算吸收系數(shù)的性能分析

表1數(shù)據(jù)明顯看出，GPU實現(xiàn)的插值運算速度比CPU的快很多。對于2個測試對象，在GPU上通過軟件插值，加速倍數(shù)TCPU/T*GPU達到14～20倍，不同射束方向加速倍數(shù)有所不同。如果采用GPU的硬件插值，加速倍數(shù)還會進一步提高25%～40%。另外，花費在傳輸吸收系數(shù)u上的時間TGPU-T*GPU，占程序在GPU上執(zhí)行總時間的大半部分，有時甚至比計算吸收系數(shù)u?的核函數(shù)執(zhí)行時間T*GPU還長。然而，在TPS中，吸收系數(shù)u只需要傳輸一次，而后無論如何改變射束方向都不需要數(shù)據(jù)傳輸，這使得傳輸時間的比重也降低了。

2.2 射線追蹤過程的性能分析

在獲得虛擬體模后，對射線追蹤在GPU上的實現(xiàn)過程進行了單獨測試。表2列出了射線追蹤的GPU和CPU實現(xiàn)的用時比較。從表中可以看出，加速倍數(shù)TCPU/T*GPU達到了13～16倍。這里TGPU包含了把GPU上計算得到的DRR圖像（即Fi）拷貝到CPU上的用時。由于DRR圖像數(shù)據(jù)是二維數(shù)據(jù)，數(shù)據(jù)傳輸時間不到10%。

表1 插值運算在GPU（GeForce G210）和CPU上實現(xiàn)的用時比較

表2 射線追蹤在GPU（GeForce G210）和CPU上實現(xiàn)的用時比較

2.3 DRR圖像重建質(zhì)量分析

目前，GPU和CPU上的浮點運算并不具有相同精度[9]，因而有必要對GPU和CPU重建出的DRR圖像質(zhì)量進行比較。而且，GPU實現(xiàn)中的軟件插值和硬件插值過程也可能具有不同精度，同樣有必要比較GPU實現(xiàn)中采用不同插值過程計算虛擬體模吸收系數(shù)獲得的DRR圖像質(zhì)量。

圖3顯示了射束方向為0°時，模型2的DRR圖像。圖 3（a）由 CPU 實現(xiàn)獲得；圖 3（b）和圖 3（c）分別是采用軟件插值和硬件插值的GPU實現(xiàn)獲得。3幅圖像很難被區(qū)分開，表明GPU上重建獲得的DRR圖像，滿足治療計劃階段確定腫瘤靶區(qū)位置以及設(shè)定射線方向和射野形狀的要求。

圖3 模型2的DRR重建圖像

圖4 當機頭轉(zhuǎn)角θG=90°時，模型1的DRR重建圖像

考慮到DRR圖像被用于自動校正擺位或設(shè)置誤差，以及驗證劑量分布，對CPU和GPU重建的DRR圖像質(zhì)量進行了定量分析。當以CPU重建的DRR圖像作為參考標準，GPU重建的DRR圖像僅有數(shù)量不到0.1%的像元的相對誤差超過1%。

圖4顯示了機頭轉(zhuǎn)角θG=90°情況下，模型1由CPU重建的DRR圖像。其中亮點標出了GPU重建DRR圖像與CPU重建DRR圖像相對誤差≥0.5%的像元。可以看出，在GPU上軟件插值重建的DRR圖像只有很少的像元相對誤差≥0.5%，而硬件插值重建的DRR圖像具有更多的誤差像元，而且，所有誤差像元都是分布在人體的外輪廓。

圖 5 顯示了機頭轉(zhuǎn)角 θG=-45°，0°和 90°情況下模型1的DRR圖像。圖像中的亮點則表示了在GPU上通過硬件插值重建的DRR圖像和CPU重建的DRR圖像，相對誤差不小于0.5%的像元，同樣也發(fā)現(xiàn)大量的誤差像元出現(xiàn)在人體的外輪廓上。在機頭轉(zhuǎn)角θG=-45°時的DRR圖像上，有一些誤差像元出現(xiàn)在人體的左肩，這是由于離散化的體模表面是有限大小的體元而不是平滑的。因此，當射線強度的衰減主要來自于表面體元時（例如射線與表面幾近相切），該射線對應(yīng)的DRR圖像像素值對插值過程和射線追蹤過程中的浮點運算精度就比較敏感。

圖5 當機頭角θG=-45°、0°和90°時，模型1的DRR重建圖像

3 結(jié)束語

提出了基于GPU并行計算的DRR圖像重建方法，在擁有16個處理核的低端GPU GeForce G210上，在保證數(shù)字重建圖像質(zhì)量的情況下，重建效率大大提高，重建的速度足以應(yīng)用到實時治療計劃中。當然，要達到實時治療計劃，還需要對整個治療計劃系統(tǒng)中的其他計算模塊（如劑量計算等）進行加速。在GeForce G210上測試完基于GPU的應(yīng)用后，又在擁有72個處理核的GeForce GT320上做了相同的測試，對于衰減系數(shù)的計算和射線追蹤這2個過程，加速倍數(shù)能分別提高到50倍和90倍。隨著更新的具有更多處理核的GPU的問世，必將促進實時治療計劃的產(chǎn)生。

[1]Khan F M，Potish R A.Treatment planning in radiation oncology[M].Philadelphia：Lippincott Williams&Wilkins，2000：11-29.

[2]Garland M， Grand S L，Nickolls J，et al.Parallel computing experiences with CUDA[J].IEEE Micro，2008（28）：13-27.

[3]Frishman Y，Tal A.Online dynamic graph drawing[J].IEEE Trans Visualization Comput Graphics，2008，14（4）：727-740.

[4]Meel J A V，Arnold A，F(xiàn)renkel D，et al.Harvesting graphics power for MD simulations[J].Mol Simul，2008（34）：259-266.

[5]NVIDIA Corporation.NVIDIA CUDA compute unified device architecture[EB/OL].[2011-08-11].http：//developer.download.nvidia.com/compute/cuda/3_1/toolkit/docs/VIDIA_CUDA_C_Programing Guide_3.1.pdf 2010.

[6]Greef M，Crezee J，Eijk J，et al.Accelerated ray tracing for radiotherapy dose calculations on a GPU[J].Med Phys，2009，36（9）：4095-4102.

[7]Wang J，Wu Z，Zhang C，et al.Enhanced digitally reconstructed radiographs generated by ray tracing[J].J Biomed Eng，2005，22（1）：125-128.

[8]Killoran J H，Baldini E H，Beard C J，et al.A technique for optimization of digitally reconstructed radiographs of the chest in virtual simulation[J].Int J Radiat Oncol Biol Phys，2001，49（1）：231-239.

[9]Anderson J A，Lorenz C D，Travesset A.General purpose moleculardynamicssimulationsfully implemented on graphics processing units[J].J Compt Phys，2008（227）：5342-5359.