基于图形处理器的射线追踪数字重建影像方法

刘 操，汪 俊，吴章文，勾成俊，侯 氢

（四川大学原子核科学技术研究所辐射物理及技术教育部重点实验室，四川 成都 610064）

0 引 言

数字重建影像（digitally reconstructed radiography，DRR）是肿瘤放射治疗计划系统（treatment planning system，TPS）里的重要组成部分。DRR图像可用于指导放射肿瘤医生和物理师确定肿瘤靶区位置和设定射线的方向以及射野形状[1]，还用于校正摆位误差和验证射野辐射注量或剂量分布。这些都是交互性的试错过程，然而通过计算机断层扫描图像（computed tomography，CT）重建DRR图像又是一个计算密集过程。因此，发展能快速获得DRR图像的方法具有实际应用价值。

加速DRR重建速度，可通过简化DRR重建的计算模型或运用计算机领域的新技术。基于可编程图形处理器（graphics processing unit，GPU）的并行计算，在图形图像及科学计算等领域得到广泛应用[2-5]。与传统的并行计算构架如集群计算相比，基于GPU的并行运算可在个人电脑上实现，且在相同计算效率情况下花费更少。

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

1 算法和GPU实现

就算法本身而言，基于中央处理器（central processing unit，CPU）和GPU的射线追踪数字影像重建过程具有相似的处理流程，如图1所示。然而，基于CPU的实现不需要CPU主内存与GPU显存之间的数据拷贝过程。总的来说，处理流程分为4个部分。第1部分，载入体模的CT数据，然后转化为体元的X线吸收系数u，转换公式[7-8]为

式中：uw——水的线性吸收系数；

HUi——第i个体元的CT值。

为便于讨论，将ui表示的体模称为初始体模。由于ui是体模在人体坐标系下体元的吸收系数，而射线追踪将在射束坐标系下进行；因此，第2部分就是在射束坐标系下创建虚拟体模，并计算出体元的吸收系数；获得虚拟体模后，第3部分则计算射线沿投影平面上的像素穿过虚拟体模到达X线源过程中的射线衰减；第4部分在设备上显示DRR图像。

图1 基于GPU实现射线追踪算法重建DRR图像的流程图

第2部分和第3部分是DRR重建中计算量最大的部分，接下来将详细说明这两部分在GPU上的实现过程。

1.1 计算吸收系数˜的GPU实现

当射束以一定角度倾斜照射体模时，创建一个表面与射束方向垂直的虚拟体模，有利于进行射线追踪，如图2所示。射束坐标系和人体坐标系在图2中也做了标示。虚拟体模中处于位置的体元α，在人体坐标系下的位置ra可以表示为

式中：θG——机头的转角；

θT——治疗床的转角；

θC——准直器的转角；

T（θG，θT，θC）——射束坐标系与人体坐标系间的3×3转换矩阵。

如果ra在初始体模中处在以其相邻8个体元中心作为顶点的立方体中，吸收系数可通过对初始体模中周围体元的吸收系数u进行线性插值获得。

图2 射线追踪示意图

通常，DRR图像重建过程中处理的体元数目很多，虚拟体模具有更多的体元。然而在CPU上实现上述插值运算，要对虚拟体模的每个体元进行循环。相比而言，GPU实现中的循环被并行的线程所替代。在GPU实现过程中，具有与体元相同数目的线程通过GPU核函数进行分配和运行[5]，每个线程完成对1个体元的插值运算。初始体模中体元的吸收系数u拷贝到GPU后，与纹理内存绑定。由于纹理单元具有硬件插值功能，文中对直接运用纹理单元做线性插值，称为硬件插值；对提取吸收系数u再进行线性插值，称为软件插值。两种插值方式除了声明不同外，具有相同的执行流程。在每个线程中，虚拟体模的体元在人体坐标系中的位置，通过式（2）计算获得。为缩短读取数据的时间，转换矩阵 T（θG，θT，θC）存放在GPU的常量内存中。

1.2 射线追踪过程的GPU实现

DRR图像是通过计算射线穿过虚拟体模后，透射到投影平面上的射线强度获得的。投影平面被离散化为像素单元，每个像素与1条射线相关联。像素i处的强度Fi计算式为

式中：I0——X线源的强度；

li——像素i到射线源的积分路径；

liα——第i条射线穿过体元α的几何长度。

像素值即灰度值的大小，定义为I0-Fi。对于DRR图像重建中的每个像素点，相关射线都从像素的中心位置向射线源进行追踪，追踪的步长与虚拟体模中体元的大小相同。与计算吸收系数u˜的插值运算类似，在GPU计算中，线程与像素点对应，完成式（3）表示的数值积分。

2 结果及分析

采用CUDA C编程，实现了基于GPU的DRR射线追踪重建算法。采用头颈部体模（模型1）和骨盆体模（模型2），对GPU实现进行了性能测试。测试模型的机头转角分别选取0°，90°和-45°，准直器和治疗床的转角都为0°。测试平台是拥有IntelCore2 Q8400的CPU和NVIDIA GeForce G210显卡的个人电脑，程序的编译和运行环境为Windows操作系统。

2.1 计算吸收系数的性能分析

表1数据明显看出，GPU实现的插值运算速度比CPU的快很多。对于2个测试对象，在GPU上通过软件插值，加速倍数TCPU/T*GPU达到14～20倍，不同射束方向加速倍数有所不同。如果采用GPU的硬件插值，加速倍数还会进一步提高25%～40%。另外，花费在传输吸收系数u上的时间TGPU-T*GPU，占程序在GPU上执行总时间的大半部分，有时甚至比计算吸收系数u˜的核函数执行时间T*GPU还长。然而，在TPS中，吸收系数u只需要传输一次，而后无论如何改变射束方向都不需要数据传输，这使得传输时间的比重也降低了。

2.2 射线追踪过程的性能分析

在获得虚拟体模后，对射线追踪在GPU上的实现过程进行了单独测试。表2列出了射线追踪的GPU和CPU实现的用时比较。从表中可以看出，加速倍数TCPU/T*GPU达到了13～16倍。这里TGPU包含了把GPU上计算得到的DRR图像（即Fi）拷贝到CPU上的用时。由于DRR图像数据是二维数据，数据传输时间不到10%。

表1 插值运算在GPU（GeForce G210）和CPU上实现的用时比较

表2 射线追踪在GPU（GeForce G210）和CPU上实现的用时比较

2.3 DRR图像重建质量分析

目前，GPU和CPU上的浮点运算并不具有相同精度[9]，因而有必要对GPU和CPU重建出的DRR图像质量进行比较。而且，GPU实现中的软件插值和硬件插值过程也可能具有不同精度，同样有必要比较GPU实现中采用不同插值过程计算虚拟体模吸收系数获得的DRR图像质量。

图3显示了射束方向为0°时，模型2的DRR图像。图 3（a）由 CPU 实现获得；图 3（b）和图 3（c）分别是采用软件插值和硬件插值的GPU实现获得。3幅图像很难被区分开，表明GPU上重建获得的DRR图像，满足治疗计划阶段确定肿瘤靶区位置以及设定射线方向和射野形状的要求。

图3 模型2的DRR重建图像

图4 当机头转角θG=90°时，模型1的DRR重建图像

考虑到DRR图像被用于自动校正摆位或设置误差，以及验证剂量分布，对CPU和GPU重建的DRR图像质量进行了定量分析。当以CPU重建的DRR图像作为参考标准，GPU重建的DRR图像仅有数量不到0.1%的像元的相对误差超过1%。

图4显示了机头转角θG=90°情况下，模型1由CPU重建的DRR图像。其中亮点标出了GPU重建DRR图像与CPU重建DRR图像相对误差≥0.5%的像元。可以看出，在GPU上软件插值重建的DRR图像只有很少的像元相对误差≥0.5%，而硬件插值重建的DRR图像具有更多的误差像元，而且，所有误差像元都是分布在人体的外轮廓。

图 5 显示了机头转角 θG=-45°，0°和 90°情况下模型1的DRR图像。图像中的亮点则表示了在GPU上通过硬件插值重建的DRR图像和CPU重建的DRR图像，相对误差不小于0.5%的像元，同样也发现大量的误差像元出现在人体的外轮廓上。在机头转角θG=-45°时的DRR图像上，有一些误差像元出现在人体的左肩，这是由于离散化的体模表面是有限大小的体元而不是平滑的。因此，当射线强度的衰减主要来自于表面体元时（例如射线与表面几近相切），该射线对应的DRR图像像素值对插值过程和射线追踪过程中的浮点运算精度就比较敏感。

图5 当机头角θG=-45°、0°和90°时，模型1的DRR重建图像

3 结束语

提出了基于GPU并行计算的DRR图像重建方法，在拥有16个处理核的低端GPU GeForce G210上，在保证数字重建图像质量的情况下，重建效率大大提高，重建的速度足以应用到实时治疗计划中。当然，要达到实时治疗计划，还需要对整个治疗计划系统中的其他计算模块（如剂量计算等）进行加速。在GeForce G210上测试完基于GPU的应用后，又在拥有72个处理核的GeForce GT320上做了相同的测试，对于衰减系数的计算和射线追踪这2个过程，加速倍数能分别提高到50倍和90倍。随着更新的具有更多处理核的GPU的问世，必将促进实时治疗计划的产生。

[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，Frenkel 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.