多源γ 射束立体定向放射治疗计划系统的剂量计算加速算法研究

【作 者】曹国刚

上海世鹏实验室科技发展有限公司，上海，200072

1951年，瑞典著名神经外科专家Leksell教授提出立体定向放射外科(Stereotactic Radiosurgery)的概念。1968年，Leksell教授与Larsson、Rexed等核物理专家以及放射治疗学专家、解剖学家通力合作，在瑞典的Uppsala大学首先研制成功世界上第一台伽玛刀，翻开了立体定向放射神经外科的新纪元[1-2]。立体定向放射治疗是指采用立体定向手段进行的高精度放射治疗的技术。与传统的放射治疗相比，它在增加肿瘤照射剂量的同时，分散了对健康组织的照射剂量，大大提高了肿瘤放射治疗的临床疗效。经过各国专家几十年的努力，立体定向放射治疗现已形成了静态式、旋转式和弧形等中心三大类方式，出现了头部伽玛刀、旋转式头部和体部伽玛刀和头部和体部X刀等多种装置。

放射治疗的质量取决于治疗时的吸收剂量分布，而剂量计算是目前能全面给出吸收剂量分布的唯一方法，因此高精度地计算人体的三维吸收剂量分布就成为评价治疗计划优劣的先决条件。放射治疗计划系统中，剂量计算总是在计算速度和计算精度间寻找平衡点。基于Monte Carlo粒子输运理论的剂量计算是目前最精确的算法，但计算所需的时间是惊人的，极少能用于临床[3-4]；基于Collapsed Cone卷积叠加模型的剂量计算方法，常用于X射线计划系统[5-6]；而查表插值法等近似模型则适合 γ - 射线的性质及照射方法[7]。

本文以陀螺刀（一种陀螺旋转式多源 γ 射束立体定向放射治疗系统）治疗计划系统为例，以查表插值法近似模型为基础，介绍一种基于锥形坐标系和多核并行化的剂量计算加速算法。该治疗计划系统由上海世鹏实验室科技发展有限公司自主研发，已通过国家相关机构检验。

1 多源 γ 射束立体定向放射治疗计划系统

多源 γ 射束立体定向放射治疗系统(Multi-Source Stereotactic Radiotherapy System with Gamma Beam)，俗称伽玛刀(Gamma Knife)，采用几何聚焦的原理，利用立体定向装置、CT模拟定位、计算机立体成像与三维治疗计划系统等技术手段，将由多个钴-60放射源经准直器形成的很细的γ-射线在三维空间聚焦于一点，要求该点与病灶完全重合，形成一个高剂量分布区，称之为靶点(shot)。由于靶点上的能量极高，可以在短时间内将靶点击毁；而靶点之外的健康组织所受的照射剂量却很小，从而达到了比手术切除更好的治疗效果[2,8]。

治疗计划系统（treatment planning system,TPS）是伽玛刀的重要组成部分，是一种融合了诸如图像处理和分析、数值计算和可视化、计算机图形学、软件工程、放射物理学以及肿瘤学等多种学科的复杂软件系统，医生藉此可以方便迅速地制定出合理的治疗计划方案。治疗计划系统一般由以下几个功能模块组成：剂量学数据管理、病例数据管理、病例数据导入和导出、框架定位、图像配准与融合、轮廓勾画、虚拟人体建模、计划规划与优化、剂量计算与计划评估、治疗计划输出与存档等[9]。

目前的伽玛刀治疗计划系统均采用正向设计、手工规划的方法，计划者负责选择靶点相关的参数，由计算机进行剂量计算。这是一个反复试误(try and error)的过程，在找到合适的方案之前要一直进行下去。这种方法中剂量计算速度决定计划设计时间[10]，因此提高剂量计算速度成为计划系统的关键。

2 治疗计划系统中的剂量计算

所谓剂量计算，就是计算粒子或光子射线照射到人体（或体模）中能量吸收的空间分布。由于粒子或光子线与人体组织的相互作用比较复杂，射线的种类、能量和照射方式不同，其形成的剂量分布会有很大的差别。在伽玛刀治疗计划系统中，根据 γ - 射线的性质及照射方法，采用剂量计算近似模型，用实际测量的剂量学数据进行查表插值的方法计算剂量分布的近似估计。该方法主要包括单源剂量计算、多源剂量分布叠加和多靶点剂量分布叠加。

2.1 单源剂量计算方法

单个静态 γ - 射线源对病人体内某一点贡献的剂量可以用下式表示：

式中：D(x,y,z,p,n)为当源在位置p，使用第n个准直器照射时，体内某点(x,y,z)的剂量值；x,y,z为病人体内某点的坐标；p为源位置；n为准直器编号；DF(n)为第n个准直器的标定剂量率；TMR(deff,n)为第n个准直器，深度为deff时的组织最大剂量比(tissue maximum ratio,TMR)；deff为计算点的有效深度(包含非均匀组织的电子密度校正）；OAR(d,r,n)为第n个准直器，离轴距离为r，深度为d的离轴比(off axis ratio,为深度校正几何因子；SAD为源轴距，即源至等中心的距离；△d为计算点在中心轴上投影距离等中心的距离，远离源方向时取正，反之取负。

2.2 多源剂量计算方法

对于多个放射源或放射源为运动的情况，采用离散累加的方式计算剂量，所以计算点的实际剂量是多个放射源位置的剂量贡献之和：

其中，D(x,y,z,n)为使用第n个准直器照射时，体内某点(x,y,z)的剂量值；Wp是源在位置p的权重。

2.3 多靶点剂量计算方法

如果有m个靶点，病人体内的剂量是多个靶点的加权和，即：

其中，D(x,y,z)为体内某点(x,y,z)的剂量值，wi是第i个靶点的权重，ni为第i个靶点的准直器编号。

3 剂量计算加速算法

传统剂量计算方法通常在直角坐标系中进行，不同计算点需重复计算深度、有效深度、离轴距离等中间结果多次，且这些中间结果无法重复使用，导致计算速度慢。因此，我们考虑建立一个“锥形”坐标系，重复利用中间结果，提高计算速度。此外，剂量计算过程中，需多次计算源在不同位置的剂量分布，是一个典型的计算密集型过程，十分耗时。随着多核处理器的可用性和性价比的不断提高，越来越多的应用将采用多核架构，多核并行化的剂量计算的研究更具有实用价值。

3.1 基于锥形坐标系的剂量计算

首先建立锥形坐标系。通过空间几何变换将机器坐标系变换到以源为原点的直角坐标系(x,y,z)，锥形坐标系(d,φ,θ)与上述直角坐标系关系如图1。两个坐标系的转换公式如下：

其中，当z=0时，θ=0，φ=0；当x=0，y>0时,θ=90；当x=0，y<0时，θ=270。

在锥形坐标系中进行剂量计算时，可以减少对deff、d、r和△d等中间结果的大量重复计算，加快计算速度。在该坐标系下进行剂量计算后，只需将其结果转换回直角坐标系下即可。

图1 锥形坐标系与直角坐标系Fig.1 Cone coordingate system and cartesian coordinate system

3.2 剂量计算的多核并行化

多核技术是近年来全球计算机技术发展的重要内容，自从英特尔在2006年底发布了全球第一款主流四核处理器至强5300后，计算机行业宣告正式进入多核时代。多核计算将成为一种普及的计算模式，并将深刻影响企业和消费者的使用模式[11-12]。为了提高剂量计算速度，本文提出一种多核并行化的剂量计算算法，以加快剂量计算速度。

多核并行化将每个单源剂量计算分配给一个虚拟的处理器(核)，所有处理器共享输入数据，形成数据并行的算法。每个虚拟的处理器处理对应的计算过程，无需关心其他的计算过程，但需在对共享数据进行读写时，采用同步处理保证对共享数据的互斥访问。

执行共享数据的多核并行计算时，线程不是孤立的，各线程间存在依赖关系，线程间对数据的共享可能会导致错误的行为和结果。因此，应用多线程实现并行计算必需引入同步机制协调各线程，以便安全共享数据。线程的同步机制包括临界区(Critical sections)、互斥量(Mutex)、信号量(Semaphores)等。其中，临界区和互斥量用来串行化多个线程共享数据的访问，其思想是每个独占性地访问一个资源的线程，可在访问那个资源之前锁定临界区，访问之后解除锁定。如果另一个线程也试图锁定该临界区，则该线程将被阻塞直到该临界区空闲。信号量保存有代表可用资源数量的资源数，锁定信号量会减少资源数，释放信号量则增加资源数。只有在线程试图锁定资源数为0的信号量时，线程才会被阻塞。信号量可以用来同步化同一进程中的线程，也可同步化不同进程中的线程。当系统中有多个同种类型的资源共享时，就可利用信号量在多个线程间分配这些资源。

基于多核处理器的并行化剂量计算是共享输入数据的并行化算法，其并行程序与串行程序十分相似，仅在同时计算累加剂量分布数据时用临界区保证互斥访问。

4 实验结果与分析

4.1 实验条件

实验采用直径160 mm的有机玻璃球模，并将根据该球模产生的仿真CT数据输入TPS；采用标称尺寸为12 mm，标称吸收剂量率为200cGy/min的聚焦野剂量学数据；将靶点置于球模中心，滚筒旋转120o，持续照射5分钟。采用不同方法进行剂量计算，比较不同方法的点剂量误差、面积重合率和计算速度等参数。

实验环境：HP xw6600图形工作站；2颗英特尔至强E5405四核处理器，4G内存；Windows XP Professional SP3。

4.2 评价方式

分别在直角坐标系和锥形坐标系下计算剂量分布，对比剂量分布中的点剂量误差和面积重合率。用点剂量误差评价某个固定点在两种计算方法中绝对剂量的差异，而面积重合率评价两种计算方法中相对剂量分布的差异。

点剂量绝对误差和相对误差定义分别如下：

其中，△D为某点的剂量绝对误差，△DR为该点的剂量相对误差，D1和D2分别为在直角坐标系和锥形坐标系下计算所得的该点剂量值。

面积重合率定义如下：

其中，△S为某指定等剂量线的面积重合率，S1和S2分别为在直角坐标系和锥形坐标系下计算所得的该等剂量线所包围区域的面积，SC为S1和S2两部分面积的重合部分。

4.3 不同坐标系下计算结果与分析

根据过球模中心CT层面的剂量分布结果，统计点剂量误差和面积重合率。两种坐标系下剂量计算结果的点剂量误差和面积重合率如表1和表2，在两种坐标系下靶点计算量与计算时间结果如表3，表中各符号含义同4.2节。

从上述数据可以看出，两种计算方法在点剂量误差和面积重合率两个方面的差异均非常小，可以忽略不计；且新方法可以极大提高计算速度，与原有方法相比，该方法可以节约65%的时间。

表1 点剂量误差Tab.1 Error of point dose

表2 面积重合率Tab.2 Overlap rate of area

表3 靶点剂量和计算时间Tab.3 Point dose of shot center and calculation time

表4 多核并行计算时间与速度Tab.4 Time and speed in multi-core based parallelized dose calculation

4.4 多核并行化剂量计算结果与分析

表4列出了两种计算方法使用不同数量核计算所需的时间和速度。其中“传统”表示不采用多核并行的传统方法，表中分别将两种方法的传统计算方法速度定义为1.0。

图2为计算速度与核数量之间的关系，图中横坐标为剂量计算使用核的数量，纵坐标为计算速度，虚线、点虚线、实线分别代表直角坐标系、锥形坐标系和理想情况。

图2 多核并行化剂量计算中速度与核数量的关系Fig.2 Relationship of speed and number of cores in multi-core based parallelized dose calculation

从上述图表可以看出，多核并行化的剂量计算只使用一个核进行计算时，与不采用多核并行的传统计算方法耗时相当，无明显差异；多核并行化的剂量计算速度与所使用核的数量接近正比，但比理想情况稍差，原因是多线程运算时，线程间切换和数据同步需要占用少量的时间。因此，多核并行化的剂量计算方法更进一步提高了剂量计算速度，其速度提升与使用的核接近正比，当采用8核进行并行计算时，其计算速度可提升7-8倍。

5 小结

本文提出一种基于锥形坐标系和多核并行化的剂量计算加速算法。基于锥形坐标系，该方法在保证计算精度的前提下可以极大提高剂量计算速度；同时利用多核并行化计算方法，还可进一步提高剂量计算的速度，其速度提升与使用的核接近正比。从实验结果和分析看出，该方法可以将计算时间从传统计算方法的584.3秒减少到30.3秒，节约95%的时间，即速度提升到传统方法的20倍左右。在实际应用中，还可以通过减小计算区域和降低计算网格精度等方法，得到剂量计算的快算结果，在可接受范围内降低计算精度，进一步提高计算速度。

[1]陈炳桓.立体定向放射[M].北京: 北京出版社,1994

[2]严玉龙.立体定向放射疗法的治疗规划——适形放疗与立体定向放射外科治疗规划的理论与技术研究[D].东南大学博士论文,1997

[3]Kawrakow I,Fippel M.Investigation of variance reduction techniques for Monte Carlo photon dose calculation using XVMC[J].Phys.Med.Biol,2000,45: 2163-2183.

[4]Kawrakow I,Fippel M.VMC++,a MC algorithm optimized for electron and photon beam dose calculations for RTP[C].Proceedings of the 22th annual EMBS international conference,2000,1490-1493.

[5]Miften M,Wiesmeyer M,Monthofer S,et al.Implement at ion of FFT convolution and multigrid superposition models in the FOCUS RTP system[J].Phys.Med.Biol.,2000,45(4): 817- 833.

[6]Zhengdong Zhou,Wei Song.A Study of Dose Calculation Based on Collapsed Cone Convolution Algorithm[C].BMEI2010,2010,1331-1335.

[7]罗浩,郭吉丰,刘涵.CT机器人伽玛刀中基于三维区域生长算法的剂量快速搜索计算研究[J].中国医疗器械杂志,2005,29(3):161-163.

[8]于金明,殷蔚伯,李宝生主编.肿瘤精确放射治疗学[M].山东科学技术出版社,2004

[9]於文雪,李松毅,罗立民.钴-60三维精确适形放射治疗计划系统的研制.中国医疗器械杂志[J],2007,31(6): 407-410.

[10]周正东,罗立民,舒华忠,等.放射治疗计划的优化方法[J].中国医疗器械杂志,2007,31(6): 391-394.

[11]CAO Guogang,LUO Limin,RONG Chengcheng.Multicore-based parallelized differential evolution for medical image registration[C].MIPPR2009,2009,7497: 74972L1-74972L7.

[12]曹国刚.计算机辅助手术导航关键技术研究[D].东南大学博士论文,2010.