CN102036712A - 卷积/迭加剂量计算方法的、使用图形处理单元加速的用于放射治疗的实时剂量计算 - Google Patents

Abstract

一种用于放射疗法的系统，包括放射计划系统，其中所述放射计划系统包括并行处理器，所述并行处理器适于接收与具有预期放射治疗区域的身体有关的输入信息并输出用于为所述身体的所述预期放射治疗区域提供放射治疗的输出信息，其中所述并行处理器适于在确定用于提供放射治疗的所述输出信息时基于与所述身体有关的所述输入信息执行多个反向射线跟踪计算，所述多个反向射线跟踪计算中的每个包括：计算与所述身体的所述预期放射治疗区域的第一子区域对应的第一物理性质，该第一子区域被在所述源位置和所述预期放射治疗区域之间行进的射线穿过；以及在上述先提及的计算之后，计算与所述预期放射治疗区域的第二子区域对应的第二物理性质，该第二子区域与所述射线相交于比所述第一子区域更接近所述源位置的位置处。

Description

卷积/迭加剂量计算方法的、使用图形处理单元加速的用于放射治疗的实时剂量计算

相关申请的交叉引用

本申请要求于2008年5月8日提交的美国临时申请第61/126936号的优先权，因此通过引用将该申请的整体内容合并于此。

背景技术

1.发明领域

本发明涉及放射治疗系统，更具体而言，涉及体外射束放射治疗系统。

2.相关技术的讨论

放射治疗是医学使用放射来治疗恶性的细胞，例如癌细胞。这种放射可以具有电磁的形式，例如高能光子，或者具有微粒的形式，例如电子、质子、中子或者阿尔法粒子。

迄今为止，在现今实践中使用得最普遍的放射形式是高能光子。光子在人体组织中的吸收是由放射的能量及相关组织的原子结构所确定的。在放射肿瘤学中采用的能量的基本单位是电子伏特(eV)；10³eV＝1keV，10⁶eV＝1MeV。在组织中的光子吸收中可以包含三种相互作用：光电效应、康普顿效应(Compton effect)和电子偶的产生(pair production)。

在光电效应中，传入的光子将能量传递给被紧密束缚着的电子。光子几乎将它所有的能量传递给电子，然后不复存在。电子利用来自光子的大部分能量而脱离，并且开始使周围的分子离子化。这种相互作用取决于传入的光子的能量以及组织的原子数量；能量越低且原子数量越高，则光电效应越有可能发生。使光电效应在组织中占优势的能量范围约为10keV到25keV。

康普顿效应是用于癌症治疗的最重要的光子与组织的相互作用。在这种情况下，光子与“自由电子”(即没有被紧密束缚到原子的电子)碰撞。与光电效应不同，在康普顿相互作用中，光子与电子都被散射。虽然能量较低，但光子之后仍可以继续承受另外的相互作用。电子利用光子给它的能量开始离子化。康普顿相互作用的概率与传入的光子的能量成反比，并且独立于材料的原子数量。康普顿效应在25keV到25MeV的范围内占优势，并且因为大多数放射治疗在大约6MeV到20MeV的能量下执行，所以康普顿效应是临床上最经常发生的相互作用。

在电子偶的产生中，光子与原子核相互作用。光子把能量给原子核，并且在此过程中创造粒子的正电子-电子偶。正的电子(正电子)离子化直到其在正电子-电子湮没中与自由电子结合。这种正电子-电子湮没产生在相反方向上运动的两个光子。电子偶产生的概率与传入的光子的能量的对数成正比，并且取决于材料的原子数量。使电子偶占优势的能量范围为大于等于25MeV。这种相互作用在某种程度上在使用高能光子射束的常规放射治疗中发生。

随着高能线性加速器的出现，在治疗深度达到大约5cm的浅表肿瘤方面电子变为可行的选择。电子深度剂量特性的独特之处在于其产生高皮肤剂量，但是仅在几厘米后就表现出衰减。

电子在人体组织中的吸收极大的受到存在气孔(air cavities)和骨骼的影响。电子射束最通常的临床使用包括对皮肤病变(比如基底细胞癌，以及先前接收过光子放射的面积(比如乳癌病人在乳房肿瘤切除术或乳房切除术后的伤疤)的增大，以及在头部和颈部选择的结节区域)的治疗。

快速、精确的剂量计算算法对于放射疗法计划而言十分重要，因为该放射疗法计划是确保对特定病人给予所希望的剂量的唯一可用的方法。剂量计算包括两部分：源模型和传输模型。源模型提供入射流量。传输模型计算由入射流量产生的剂量，并且目前是性能瓶颈。三个主要的传输算法以增加的精确性/减少的性能的次序依次是笔状射束、迭加/卷积和蒙特卡洛(Monte Carlo)。迭加/卷积是计算外部射束放射疗法的放射剂量的当前临床标准方法。

近年来，通过使用强度调节提高了治疗质量。这种技术使用多叶准直器限定根据单一射束方向的多孔径，以提供改变穿过射束的放射强度的能力。这种技术允许使放射治疗符合目标的形状，并且在射束参数的数量急剧增长时避免临界的结构。为了确定多叶准直器最佳的一组设置，治疗计划系统必须通过剂量计算的多次迭代使射束参数的数量急剧增长的目标函数最优化。在实践中，为了对病人达到可能的最好结果，治疗计划者多次重复该最优化。因此，当对一组五个射束的单一最优化可能需要五分钟时，整个过程可能需要几小时来产生临床可接受的计划。这就限制了在临床工作流程中强度调节计划的数量和质量。

这种临床工作流程的局限性延伸到更复杂的技术，例如容积调节的弧形疗法(Otto，K.，Med.Phys.35，310-317，2008)、强度调节的弧形疗法(Yu，C.X.，Phys.Med.Biol.40，1435-1449，1995)和自适应放射疗法(Yan，D.，Vicini，F.，Wong，J.，Martinez，A，Phys.Med.Biol.42，123-132，1997)。此外，这种临床工作流程的局限性禁止实时的放射疗法；每天扫描、重新计划和治疗每个病人的能力。从Ahnesjo等人(Ahnesjo，A.，Aspradakis，M，Phys.Med.Biol.44，R99-Rl55 1999)可得到对放射疗法中的剂量计算的全面的讨论。

因此，剂量计算的计算性能是放射疗法的治疗计划的质量的限制因素。传统地，通过更快的硬件已经实现了治疗质量方面的改进。但是，摩尔定律(Moore′s law)已改变。计算机不是每18个月速度加倍，而是处理内核(kernel)的数量加倍。并且随着处理器变为多核的，图形处理单元(GPU)的多核架构获得了运行一般用途的算法的灵活性。为了根据计算机硬件的近来趋势实现许诺的性能增益，在放射剂量计算中使用的传统的串行算法应该被并行算法取代。最近，Nucletron公司宣布了在其治疗计划系统中的GPU加速，不过尚未获得其公布的细节。然而，将现有的串行算法直接分隔来产生用于多处理核的多线程行不通。这是因为线程被外包(farm out)来基于相同输入数据计算相同放射剂量，从而容易出现读/写冲突。当覆写(write on write，WOW)冲突出现时，例如可能发生仅存储了最后的写入内容，从而导致剂量计算不精确。因此，在本领域中需要通过利用多处理核上的并行计算的改进的剂量计算。此外，仍然需要实时的剂量计算。

发明内容

本发明的一些实施方式提供了一种用于放射疗法的系统，包括放射计划系统，其中所述放射计划系统包括并行处理器，所述并行处理器适于接收与具有预期放射治疗区域的身体有关的输入信息并输出用于为所述身体的所述预期放射治疗区域提供放射治疗的输出信息，所述并行处理器适于在确定用于提供放射治疗的所述输出信息时基于与所述身体有关的所述输入信息执行多个反向射线跟踪计算，所述多个反向射线跟踪计算中的每个包括：计算与所述身体的所述预期放射治疗区域的第一子区域对应的第一物理性质，该第一子区域被在所述源位置和所述预期放射治疗区域之间行进的射线穿过；以及在上述先提及的计算之后，计算与所述预期放射治疗区域的第二子区域对应的第二物理性质，该第二子区域与所述射线相交于比所述第一子区域更接近所述源位置的位置处。

本发明的一些实施方式提供了一种确定放射疗法参数的方法，包括：获得关于身体的预期放射治疗区域的信息；基于所述信息，在从所述身体到源位置的方向上计算沿射线的第一子区域和第二子区域处的物理性质，所述射线在所述源位置和所述身体之间行进并且与所述身体相交于所述预期放射治疗区域之内，沿所述射线的所述第二子区域比沿所述射线的所述第一子区域更接近所述源位置；以及基于所述计算确定用于为所述预期放射治疗区域提供放射治疗的所述放射治疗参数。

本发明的一些实施方式提供了一种计算机可读介质，当由计算机执行时，使得所述计算机实行上述方法。

附图说明

从对说明书、附图和示例的讨论中，更多目标和优点将会变得显而易见。

图1示出本发明一种实施方式的示意图。

图2示出传统的正向射线跟踪。

图3示出根据本发明一种实施方式的反向射线跟踪。

图4示出使用市售的Pinnacle系统(Philips Radiation Oncology Systems Madison WI)和本发明一种实施方式以分析法计算出的沿中心轴的每单位质量释放的总能量(Total Energy Released per unit Mass，TERMA)，其中每个TERMA以10cm的深度归一化。

图5A示出使用正向跟踪计算出的有离散化伪影的TERMA的半影片段。

图5B示出使用本发明一种实施方式计算出的没有离散化伪影的同一TERMA的半影片段。

图6A和6B分别示出根据本发明一种实施方式使用固定步长和精确放射线距离计算出的TERMA的片段。

图7A和7B示出在迭加操作中倾斜的内核和不倾斜的内核。

图8A和8B分别示出根据传统的统一采样的存储器访问模式和所计算出的小野(5mm)剂量沉积片段的图。

图9A和9B分别示出使用多分辨率网格的根据本发明一种实施方式的存储器访问模式和所计算出的小野(5mm)剂量沉积片段的图。

图10示出由市售的Pinnacle系统(Philips Radiation Oncology Systems Madison WI)和本发明一种实施方式计算出的沿中心轴的吸收剂量图，其中每个吸收剂量图以10cm的深度归一化。

图11示出由市售的Pinnacle系统(Philips Radiation Oncology Systems Madison WI)和本发明一种实施方式计算出的在10cm的深度的吸收剂量分布图，其中每个吸收剂量分布图在中点处归一化。

图12示出根据本发明另一实施方式的一种方法。

具体实施方式

下面将详细讨论本发明的一些实施方式。在描述实施方式时，为了清楚起见而采用了特定的术语。然而，本发明并非意在受限于由此而选择的特定术语。本领域的技术人员可认识到在不偏离本发明的广义构思的情况下可以采用其它等同的部件和开发出其它方法。通过引用将所举出的所有参考文献合并于此，如同每个参考文献是单独合并的一样。

图1是根据本发明一种实施方式的放射疗法系统100的示意图。放射疗法系统100包括：放射计划系统101，放射计划系统101进一步包括并行处理器102。并行处理器102适于接收与具有预期放射治疗区域的身体105有关的输入信息。并行处理器102还适于产生用于为身体105的预期放射治疗区域提供放射治疗的输出信息。并行处理器102适于在确定用于提供放射治疗的输出信息时基于所接收到的输入信息执行多个反向射线跟踪计算。每个反向射线跟踪包括：计算与身体的预期放射治疗区域的第一子区域对应的第一物理性质，该第一子区域被从源位置穿过预期放射治疗区域行进的射线穿过；以及在上述首先提及的计算之后，计算与预期放射治疗区域的第二子区域对应的第二物理性质，该第二子区域与上述射线相交于比第一子区域更接近源位置的位置处。放射计划系统101还可包括：存储器107、显示器108和输入/输出装置109。存储器107可以例如是硬盘驱动器、CD-ROM驱动器、DVD驱动器、闪存驱动器等。显示器108可以例如是液晶显示器(LCD)、阴极射线管(CRT)监视器、等离子显示器等。输入/输出装置109可包括例如鼠标、键盘、用于在网络或者数据总线上传输数据的接口等。

放射疗法系统100还可包括：放射治疗系统103，其与放射计划系统101进行通信。放射治疗系统103进一步包括放射源106。放射源106是直接发出放射束到待治疗的身体105上的源。放射源的实例可以包括：X射线源、伽马射线源、电子束源等。放射源106还可包括对射束进行校准的多叶准直器(MLC)。通过调节MLC的叶片的位置，剂量测定员可以使放射野与身体105的治疗区域的形状相匹配。在一些实施方式中可以包括其它的射束形状和/或轮廓。放射源106可以具有对应的源模型。放射系统103可以由放射治疗计划系统101控制，例如，以传递强度加以调节的放射能量，从而使得放射治疗符合身体105的预定放射治疗区域的形状。

放射治疗系统100还可包括诊断系统，该诊断系统与放射计划系统101进行通信，产生身体105的经验数据。该经验数据可被用作放射计划系统101和并行处理器102的输入信息，并且可被用于反向射线跟踪计算。诊断系统104包括传感器，以获得身体105的经验数据。诊断系统的实例可以是电子计算机X射线断层扫描技术(CT)扫描仪、核磁共振成像(MRI)扫描仪、正电子发射型断层显像(PET)扫描仪等。

身体105例如可以是人类或者动物。

示出的迭加/卷积算法用来产生对剂量分配的精确计算(Mackie，T.R.，Scrimger，J.W.，Battista，J.J.，Med.Phys.12，188-196，1985；Mackie，T.R.，Ahnesjo，A.，Dickof，P.，Snider，A，Use of Comp.In Rad.Ther.，107-1101987；Mackie，T.R.，Reckwerdt，P.J.，McNutt，T.R.，Gehring，M.，Sanders，C.，Proceedings of the 1996 AAPM Summer School，1996)。其包括两个阶段。第一，通过病人的密度表示传输入射流量来计算每个位置处的每单位质量释放的总能量(TERMA)。点r′处具有特定能量E的TERMA(即T_E(r′))被定义为将能量E的流量ψ_E(r′)用相对于水的密度ρ(r′)以及在点r′处的线性衰减μ_E(r′)加权，如下述方程1所示。

$T_E\left(r^{\′}\right)=\frac{{\mathrm{\μ}}_E\left(r^{\′}\right)}{\mathrm{\ρ}\left(r^{\′}\right)}{\mathrm{\ψ}}_E\left(r^{\′}\right)---\left(1\right)$

线性衰减系数μ_E(r′)也取决于原子材料。尽管康普顿散射在与放射疗法相关的百万伏能量范围内占优势，并且康普顿散射取决于电子密度而不是材料，但是，在临床上，对于正常人体组织，在标准CT数量和密度之间存在分段线性关系。通常而言，在此能量范围内非康普顿的相互作用被视为可以忽略。根据下述方程2，能量E在点r′处的流量ψ_E(r′)是由源焦点s和在方向r′上的能量E的入射流量ψ_E，0(r′)所决定的。

${\mathrm{\Ψ}}_E\left(r^{\′}\right)=\frac{{\mathrm{\ψ}}_{E,0}\left(r^{\′}\right)}{{\left|\right|r^{\′}-s\left|\right|}^2}{e}^{{\∫}_s^{r^{\′}}-{\mathrm{\μ}}_E\left(t\right)\mathrm{dt}}---\left(2\right)$

然后，迭加/卷积算法通过剂量沉积内核传播上述能量，以确定每个位置处的最终剂量。为了允许剂量沉积内核根据组织不均匀性来实际地进行调整，基于将迭加与卷积区别开来的下述方程3使用点之间的放射线距离d_ρ(r，r′)，。

$d_{\mathrm{\ρ}}(r,r^{\′})={\∫}_r^{r^{\′}}\mathrm{\ρ}\left(t\right)\mathrm{dt}---\left(3\right)$

根据下述方程4，点r处的剂量D(r)是由在TERMA体积上的积分并用取决于能量的剂量沉积内核K_E加权而计算出来的。标准压缩(collapsed)椎体内核是由在点和内核轴之间的放射线距离和相对角ω来标示的，并且缺少几何距离的平方的作用。

通常而言，不计算如方程4中的每个体素(voxel)的单能贡献。而是如基于方程5和6的方程7中一样，选择射线角度ω和方向v的离散组合，并且沿着使用的单个多能内核对其进行积分。由于距离平方效应，这由近似的指数内核衰减和任意单个远距离的体素急剧减少的贡献得到证明。射线表现出的立体角量上的增加抵消了距离平方效应。射线方向被选择来平衡几何形状和内核的能量因子。

$T\left(r^{\′}\right)=\underset{E}{\mathrm{\Σ}}{T}_E\left(r^{\′}\right)---\left(5\right)$

$K(d_r,\mathrm{\ω})=\underset{E}{\mathrm{\Σ}}{K}_E(d_r,\mathrm{\ω})---\left(6\right)$

$D\left(r\right)\≈\underset{\mathrm{\ω},v}{\mathrm{\Σ}}\∫T(r+\mathrm{tv})K(d_{\mathrm{\ρ}}(r,r+\mathrm{tv}),\mathrm{\ω})\mathrm{dt}---\left(7\right)$

传统地，根据方程1和2，TERMA是由投射一组将入射流量沉积至体积的射线来计算的。为了数字的精确性，大约四条射线应该穿过每个TERMA体素，应该在病人的表面处开始衰减，并且应该通过贡献给该体素的射线的总长度来对每个体素的流量进行归一化。这种归一化去除了由于发散源造成的正常反向平方流量的下降。因此必须再次应用该归一化，且通常将其应用于TERMA网格(grid)。临床上可接受的速度增强以不使剂量沉积内核倾斜来与每个体素处的射线轴对准。在这种情况下，将发散校正应用于剂量网格会更精确(Papanikolaou，N.，Mackie，T.R.，Meger-Wells，C.，Gehring，M.，Reckwerdt，P.，Med.Phys.20，1327-1336，1993)。

TERMA强烈取决于射束光谱。该光谱关于射束轴旋转对称，并且随着材料的深度而变得稳定。稳定的光谱指的是：由于较低的能量成分被优先散射出射束外并被吸收，因此光谱中较高的能量成分占优势。传统地，使用查找表来将衰减模型化，当该查找表与线性衰减查找表结合时，该查找表具有深度轴、密度轴和偏轴角。该表还需要使用固定步骤射线投射算法，该算法避免了在每一步骤都要进行传统的、成本昂贵的对指数的评估。尽管出于性能原因在临床上接受了该查找表，但该查找表假定介质是均匀的。由于不均匀的组织优先使不同的光谱衰减，因此这是不正确的。此外，固定步长和离散射线导致数字的伪象。

剂量沉积内核还依赖于每个体素处的能量谱。然而，这种效果是可以忽略的，并且单个多能内核已被展示对临床使用而言足够精确且该内核是当前医疗的标准。多能内核是通过对单能内核谱进行组合来创建的，单能内核是使用蒙特卡洛模拟产生的，该模拟强制单能光子在水球体中心处互相作用并且对沉积在该球体每处的剂量进行计数(Ahnesjo，A.，Andreo，P.，Brahme，A.，Acta.Oncol.，26，49-56，1987；Mackie，T.R.，Bielajew，A.F.，Rogers，D.W.O.，Battista，J.J.，Phys.Med.Biol.33，1-20，1988)。

图形处理算法明显地适于并行或者图形处理单元(GPU)的架构。存在适应于传统图形系统的多个可视化算法，例如射线投射、射线跟踪和体积的可视化算法。然而，当可视化算法关注线路性质(比如其整数或最大化)时，用于放射疗法的剂量计算根本上与线路和体积的相互作用密切相关。这种区别致使多个为可视化开发的并行射线跟踪算法不适用于放射疗法。尽管如在(Kriiger，J.，Westermann，R.，2003)中所做的轻易改变先前GPU的执行是不切实际的，并且将不会得到精确的结果，但传统的TERMA算法的一部分可以与体积射线投射粗略地相似。因为剂量沉积主要处理电子的相互作用，因此其与可视化算法根本不同。

根据本发明的一种实施方式，将NVIDA的统一计算设备架构(CUDA)软件开发环境与Digital Mars的D编程语言结合起来用于执行本发明，可使迭加/卷积算法适应于GPU。以下描述该适应的细节。

如同所讨论的，迭加/卷积是两阶段的算法：首先计算TERMA，然后迭加剂量沉积内核。

在TERMA计算中，可以执行标准的正向TERMA算法。该方法需要将分散的积累进行重叠，这在并行运行时将导致读写冲突。当线程同时尝试读、处理然后写入相同的存储器时，发生读写冲突，导致仅有最后的更新被记录。这极大地减少了对每个体素有作用的射线的有效数量，使其低于数字精确性所需的极限。然而，可以使用源的扩散性来创造射线组，保证一直有一个体素是分离的。当单个的组大到足够在GPU上有效运行时，这种串行化引起过多的GPU呼叫开销。传统的3D谱衰减查找表超过结构缓存值，从而降低了性能。通过单个地使分开的谱箱衰减可以提升性能，使用共享存储器来减少寄存器的使用。这减少了对小的2D结构的线性衰减系数的查找，这些线性衰减系数由能量和密度使其参数化。这可以由在GPU上的硬件加速指数来实现。当在不均匀的介质中谱正确地衰减时，也可以提高精确度。然而，数字精确性所需的射线的数量与分辨率的关系并不清楚，离散化的效果会很明显，如同将要讨论的那样，并且计算会包括总的剂量计算时间的很大一部分。

根据本发明的一些实施方式，通过反向TERMA算法可以避免这些问题。通过对方程(1)、(2)和(5)进行整理，将入射流量ψ_E，0(r′)分为谱权w_E和净流量因子ψ_E(r′)，并且定义衰减因子，获得下述方程8。

$A\left(r^{\′}\right)=\frac{1}{{\left|\right|r^{\′}-s\left|\right|}^2}\underset{E}{\mathrm{\Σ}}{w}_E\left(r^{\′}\right)\frac{{\mathrm{\μ}}_E}{\mathrm{\ρ}\left(r^{\′}\right)}{e}^{{\∫}_s^{r^{\′}}-{\mathrm{\μ}}_E\left(t\right)\mathrm{dt}}---\left(8\right)$

图2示出了传统的正向射线跟踪。图3示出了根据本发明一种实施方式的反向射线跟踪。通过从每个TERMA体素背向源投射射线，沿着路线聚集净衰减，当到达病人的边界时，可以获得早期的射线终结，从而提高性能和精确度。在计算上述定义的衰减因子时，预先编译的查找表可以用于数量μ_E/ρ。例如，基于放射治疗区域的以Hounsfield为单位的CT图像，可以从查找表查找这个数量。谱权w_E例如可以包含在用于放射治疗的放射源的模型中，并且每个谱箱可以针对给定的射线独立地衰减。

然后可以根据下述方程9通过将针对每个体素的净衰减乘以从源朝向体素的净入射流量来计算TERMA。

T(r′)＝ψ₀(r′)A(r′)            (9)

图4示出使用市售的Pinnacle系统(Philips Radiation Oncology Systems Madison WI)和本发明一种实施方式以分析法计算出的沿中心轴的TERMA能量，其中每个TERMA能量以10cm的深度归一化。

尽管方程9是O(n⁴)的算法，但由于与用于标准正向方法的O(n^-3)不同，每个线程仅向其自身的体素写入，因此避免了读写冲突。该算法允许对结构存储器联合进行读访问，从而大幅度地提高了存储性能。当仅改变流量野时，只要在相互使用或者强度调节最优化期间，就可以在多种临床分辨率下使用标准正向和反向射线跟踪算法，并且可以计算该计算方法的高额衰减量。根据本发明的一些实施方式，这种预先计算可以提供多个量级的性能改进。

大小	正向	反向	近似的衰减	强度调节
					323	23ms	2ms	1ms	1ms
643	44ms	18ms	2ms	1ms
					1283	150ms	220ms	26ms	3ms
2563	1869ms	3317ms	454ms	13ms

表1.标准正向、新的反向和近似的TERMA计算方法的性能。也包括在强度调节期间有用的仅更新入射流量野的性能。

表1中比较了多种分辨率下基于衰减量的TERMA算法相对于传统射线投射算法的性能。两种方法均物理结合了对多谱衰减的校正并且针对GPU进行了优化。在正向射线跟踪计算中使用的射线投射可能有读/写冲突。射线投射的可量测性表现出最佳位置和所需的参数调整，以维持穿过每个体素的射线的必要数量。这种调整费力且很不直观。通过反向射线跟踪计算的基于衰减量的TERMA计算表现出经验性能O(n^3.76)，该性能在其理论的O(n⁴)上有微小的改进。根据本发明一些实施方式的反向射线跟踪计算不会经受读/写冲突。基于快速的、近似的放射线深度的衰减算法使得对标准均匀的材料近似提供约为8倍的加速。在强度调节下介绍了如通常在治疗计划中的强度调节中使用的那样仅更新入射流量的性能。

反向射线跟踪计算的一个优点是消除由于射线离散化导致的伪像。这是因为，与将每个体素存储器访问的射线减少到一个从而消除相关的离散化伪像的光栅化方法(Arnanatides，J.，Woo.，A.，Eurographics′87，Conference Proceedings，1987)相似，在反向射线跟踪方法中可以使用确切的放射线路径。图5A和5B分别示出了根据本发明一种实施方式使用固定步长和精确的放射线距离所计算出的TERMA片段。

此外，物理校正的多谱衰减(Liu，H.H.，Mackie，T.R.，McCullough，E.C.，Med.Phys.24，1729-1741，1997)可被应用于反向射线跟踪计算。对于相互射束角度变化的情况，可以实现如通常用于当前临床系统那样的无需物理校正多谱衰减的更简单的变型，从而可以得到显著的性能改进。在TERMA计算中，物理校正的多谱衰减(Liu，H.H.，Mackie，T.R.，McCullough，E.C.，Med.Phys.24，1729-1741，1997)可被用于反向射线计算。当前的临床标准是固定步骤的、快速的、近似的正向方法，该方法没有物理校正的多谱衰减。正如根据本发明的一些实施方式所指出的，缓存衰减量的使用可被用来加速TERMA的计算。图6A示出了使用没有多谱衰减的标准正向射线跟踪计算出的有离散化伪像的TERMA的半影片段。图6B示出了有多谱衰减的根据本发明一种实施方式计算出的没有离散化伪像的同一TERMA的半影片段。

在计算出TERMA之后，可以应用用于计算被吸收的放射能量的量的剂量沉积内核的迭加。该迭加具有两个标准的方程式。正向方程将剂量从TERMA体素传播到周围的剂量体素。该正向方程需要计算给每个病人的体素的剂量并且会经受读写冲突，因为每个剂量体素都有多个TERMA体素做贡献。反向内核方程聚集从周围的TERMA体素给剂量体素的贡献。由于仅计算了感兴趣的体积的剂量，因此这种计算是高效的。因为剂量沉积内核是可逆的，所以这是可行的。

图7A和7B示出了在迭加操作期间倾斜的内核和非倾斜的内核。严格来说，在标准迭加中使用倾斜的内核破坏了内核可逆的假定。然而，在已知源的距离远和内核快速下降的情况下，可逆性仍是一个合理的假定从而在临床中得到使用。虽然如此，内核的快速下降也产生了在通常的临床分辨率下数字采样的问题。有两种标准的可选方案。根据下述方程10，累积内核(cumulative kernel，CK)(Ahnesjo，A.，Med.Phys.16，577-92，1989)表示从射线段到点的剂量沉积。

$\mathrm{CK}=(x,\mathrm{\ω})={\∫}_0^{x}K(t,\mathrm{\ω})\mathrm{dt}---\left(10\right)$

累积-累积内核(CCK)(Lu，W.，Olivera，G.H.，Chen，M.，Reckwerdt，P.J.，Mackie，T.R.，Phys.Med.Biol.50，655-680，2005)表示从射线段到射线段的剂量沉积。

$\mathrm{CCK}(x,\mathrm{\ω})={\∫}_0^x\mathrm{CK}(t,\mathrm{\ω})\mathrm{dt}---\left(11\right)$

对于特定的放射线深度d_r和角度ω，这两者都源于标准的点到点内核K(d_r，ω)的积分：

${\∫}_x^{x+\mathrm{\Δx}}K(v,\mathrm{\ω})\mathrm{dv}=\mathrm{CK}(x+\mathrm{\Δx},\mathrm{\ω})-\mathrm{CK}(x,\mathrm{\ω})---\left(12\right)$

${\∫}_0^{\mathrm{\Δs}}\mathrm{CK}(x+u,\mathrm{\ω})\mathrm{du}=\mathrm{CCK}(x+\mathrm{\Δs},\mathrm{\ω})-\mathrm{CCK}(x,\mathrm{\ω})---\left(13\right)$

${\∫}_0^{\mathrm{\Δs}}{\∫}_{x+v}^{x+\mathrm{\Δx}+v}K(u,\mathrm{\ω})\mathrm{dudv}=(\mathrm{CCK}(x+\mathrm{\Δs}+\mathrm{\Δx},\mathrm{\ω})-\mathrm{CCK}(x+\mathrm{\Δx},\mathrm{\ω}))-$

$(\mathrm{CCK}(x+\mathrm{\Δs},\mathrm{\ω})-\mathrm{CCK}(x,\mathrm{\ω}))---\left(14\right)$

尽管更精确(尤其在粗糙的分辨率下)，但CCK方程传统上比CK方程慢50％。然而，缓存存储器访问和提供专用的线性内插硬件的GPU的结构单元使得能够使用CCK方程且其性能降低可忽略。因此，使用累积-累积内核(Lu，W.，Olivera，G.H.，Chen，M.，Reckwerdt，P.J.，Mackie，T.R.，Phys.Med.Biol.50，655-680，2005)能够增强迭加计算。

通过移动所有索引越过一个体素，串行CPU的执行允许重新使用射线投射索引计算。但是，这阻止了内核倾斜，从而导致大的偏轴角误差，如将在表2中讨论的那样。相反，在GPU上执行的本发明的一些实施方式允许待计算的内核倾斜和非倾斜。内核倾斜传统上导致300％的性能损失(Liu，H.H.，Mackie，T.R.，McCullough，E.C.，Med.Phys.24，1729-1741，1997)，但是在GPU上执行的本发明的一些实施方式的性能损失仅为19％。

由于迭加操作在粗糙的分辨率下维持精确性并且内核在边缘处表现出快速的下降，根据本发明的一些实施方式，可以采用将每个射线近似为真实的立体角的多分辨率迭加算法。与射线不同，立体角的宽度随着几何距离而增加。在离散量中，射线的宽度与体素的宽度是成正比的。因此，通过随着几何距离增加体素的宽度，射线可以近似为立体角。这种近似也以对数方式增加了步长，这通过将计算复杂度从O(ωDT^1/3)减小到O(ωDlog(T^1/3))而提高了性能；其中ω是角的数量，D是剂量体素数量，而T是TERMA体素的数量。

与标准方法相比，根据本发明一种实施方式的多分辨率迭加表现出令人感兴趣的精确性折衷。表2将根据本发明一些实施方式的多分辨率迭加算法的精确性在多个野尺寸和内核射线采样的情况下与标准方法相对比。由于很少有TERMA被射线几何地错过，因此根据本发明一些实施方式的多分辨率迭加算法通常在半影和用于小野尺寸的低剂量区域中能够执行得更好。然而，由于大步长可以引起射束界限被模糊，因此高剂量区域的精确性可以轻微地折衷。使用相同的步长、但不使用多分辨率数据结构的作为本发明的实施方式的一种变型的多分辨率迭加算法的一种变型可表现出缓存性能下降以及可平均增加60％的平均误差。

表2.平均沉积剂量相对于最大沉积剂量D_max的误差，其野尺寸从1cm到23cm。通过具有被限定为具有比0.3D_max更大的剂量梯度的半影区域和比0.2D_max更低的低剂量区域的区域来减小误差。使用4608射线的倾斜内核采样来计算野是方形的、大小被限定在深度为10cm的参考剂量沉积。笔状射束的精确性是通过将迭加内核在半径上截短3cm来近似的。2％到5％的绝对放射量测定误差在临床上是可接受的(Ahnesjo，A.，Aspradakis，M.，Phys.Med.Biol.44，R99-R155，1999)。

通过使用根据本发明一些实施方式的多分辨率网格可以减少当由于射线采样稀疏而导致漏掉整个射束时产生的小野伪像。图8A和8B分别示出根据传统的统一采样的存储器访问模式和所计算出的小野(5mm)剂量沉积片段的图。图9A和9B分别示出使用多分辨率网格的根据本发明一种实施方式的存储器访问模式和所计算出的小野(5mm)剂量沉积片段的图。

然而，更大的步长可能减小剂量沉积内核的精确性。当邻近的体素横穿不同的粗糙分辨率的体素时，产生伪像。应注意本实施具有固有的同向性，这将减少不均匀角度采样的益处。

使用测定体积的最大强度投影(MIP)映射(Williams，L.，SIGGRAPH Comput.Graph.17，3，1-11，1983)，可以实施多分辨率算法，因为其计算高效且具有良好的缓存性能。分辨率的改变可被限于之前每个步骤的最大值并且可仅允许其在分辨率较粗糙的体素边界处发生。这样可以防止TERMA体素多次给同一剂量体素做贡献。由于增加倾斜将会增加所需的寄存器以致超过性能阈值，从而将GPU占用从25％降为17％，所以内核倾斜不会被合并。

可以使用多个策略来优化统一计算设备架构(CUDA)的性能。CUDA的执行模型是线程的3D块中的2D网格，该线程执行一项函数(称为内核)。每个块代表一个并行工作单元，并且因此在尺寸上有限制。使用NVIDIA的CUDA占用计算器优化了块线程计算。为了体积处理，可以使用在x和y方向上的线程到体素的1∶1的映射。通过以增加的z索引使体素函数循环可以处理z方向。步幅可以是z块大小，其维持了线程空间的凝聚力。通过立方体状块尺寸增强的提高了的线程空间的凝聚力可以减小缓存遗漏从而提升性能。通常在结构中可以缓存所有的输入阵列数据，这在迭加的情况下可以提供约2倍的性能提升。可以使用共享的存储器来缓存阵列多分辨率体积结构。共享的存储器还可以提供用于多谱TERMA计算的寄存器方面的性能改进。最多可以选择21个能量箱，因为它们不仅足够用于高能量射束而且没有存储冲突。

在根据本发明实施方式的TERMA计算和迭加这两个阶段中，与线性光栅化相似(Amanatides，J.，Woo.，A.，Eurographics′87，Conference Proceedings，1987)，可以用精确地放射线路径方法代替标准的固定步长射线投射算法。

图10示出由市售的Pinnacle系统(Philips Radiation Oncology Systems Madison WI)和本发明一种实施方式计算出的沿中心轴的吸收剂量图，其中每个吸收剂量图以10cm的深度归一化。图11示出由市售的Pinnacle系统和本发明一种实施方式计算出的在10cm的深度的吸收剂量分布图，其中每个吸收剂量分布图在中点处归一化。

对传输算法(比如迭加/卷积)的定量分析可以由对来自用于放射治疗的放射源的源模型处的入射流量的强烈依赖而复杂化(Mohan，R.，Chui，C.，Lidofsky，L.，Med.Phys.12，592-597，1985；McNutt，T.R.：Dose Calculations，Philips white paper)。源模型转而需要被优化，以使测量剂量和计算剂量之间的误差最小。因此，在传输算法中，源模型经常可以掩盖误差。此外，市售系统可以包括单独的电子污染模型。在使用了本发明一种实施方式引导进行的初步实验中，仅考虑简单的源模型。初步实验产生了与那些使用市售治疗计划系统的模型参数相似的结果，同时提供了如表3中所示的数量级的速度上的改进。

表3.在立方水假体上的剂量引擎性能的比较。速率是每秒钟的体积(VPS)。倾斜的32射线内核比不倾斜的80内核更精确(参见表2)。自适应的多网格方法(未实施)在临床数据上具有2倍的加速因子。

在一些实验中，使用基于来自Varian 6EX线性加速器的临床数据的Pinnacle(Philips-Madison，WI)治疗计划系统产生了临床参考数据，提供剂量沉积、TERMA、传输阵列、谱、质量衰减表和单一能量剂量沉积内核。所有实验在AMD处理器254(双核，2.8GHz)上运行。使用标准迭加/卷积引擎，在没有其它程序运行的情况下至少重复十次计时实验。

使用CPU的高性能硬件计数器，将针对本发明实施方式的计时结果重复多次。实验在有单个NVIDIA GeForce GTX 280的3GHz的奔腾4上执行。80-射线内核使用了10顶点、用8方位角。32-射线内核使用了4个顶点乘以8个方位角方向。倾斜的72-射线内核使用了6个顶点乘以12个方位角方向，这提供了更高的精确性(参见表3)。总的剂量引擎性能是由源模型、TERMA和迭加性能的总和所确定的，不包括可再次使用的计算，比如修平滤波器结构和TERMA衰减量。

所有试验都在边长为25.6cm的立方水假体上执行。试验在代表标准临床工作量的具有64³个体素的体积上、和在具有128³个体素的额外的高分辨率体积上执行。作为本发明一种实施方式的多分辨率迭加方法利用测定体积的mip-map(一种纹理映射技术)并且可以比传统的迭加执行起来快2倍到3倍以及如示出的那样的具有对高分辨率调节得更好的性能。

图12示出根据本发明另一实施方式的一种方法。本方法包括：获得关于身体的预期放射治疗区域的输入信息1201；在方框1202中，在从所述身体到源位置的方向上，基于输入信息1201计算沿射线的第一和第二子区域的物理性质，射线从源位置出发穿过身体并且与身体交叉于预期放射治疗区域，沿射线的第二子区域比沿射线的第一子区域更接近于源位置；以及在方框1203中，确定用于为预期放射治疗区域提供放射治疗的输出信息1204。

在方框1202中，基于在第一子区域处计算的物理性质计算第二子区域处的物理性质。物理性质是指表示从入射的放射束提取的能量的相对数量的衰减因子。如上讨论的衰减因子是一个实例。此外，输入信息1201可以包括来自身体的经验数据。测量数据可以包括身体的预期治疗区域的至少一个电子计算机X射线断层扫描技术(CT)图像。此外，在方框1202中可以使用至少一个预编译的查找表。

方框1203还包括计算被预期放射治疗区域吸收的放射能量的量的步骤1204。使用统一采样、可变采样、多分辨率网格或者其变型中的至少一个以覆盖预期放射治疗区域，可以计算出被吸收的放射能量。此外，多分辨率网格被指定使得来自源位置穿过身体行进的射线可以近似于立体角。

然而，本发明的另一实施方式可以包括计算机可读介质，其包括软件，当该软件被计算机执行时使计算机执行上述方法。

在描述本发明的实施方式时，为了清楚起见而采用了特定的术语。然而，本发明并非意在受限于由此而选择的特定术语。本领域的技术人员可想到在不偏离本发明的情况下可根据上述构思对本发明的上述实施方式做出修改或变型。因此可以理解，可以在所附权利要求和与其等同方案的范围内以与所明确描述的方式不同的方式来实施本发明。

Claims (27)

1.一种用于放射疗法的系统，包括放射计划系统，

其中所述放射计划系统包括并行处理器，所述并行处理器适于接收与具有预期放射治疗区域的身体有关的输入信息并输出用于为所述身体的所述预期放射治疗区域提供放射治疗的输出信息，

其中所述并行处理器适于在确定用于提供放射治疗的所述输出信息时基于与所述身体有关的所述输入信息执行多个反向射线跟踪计算，所述多个反向射线跟踪计算中的每个包括：

计算与所述身体的所述预期放射治疗区域的第一子区域对应的第一物理性质，该第一子区域与在所述源位置和所述预期放射治疗区域之间行进射线相交；以及

在所述先提及的计算之后，计算与所述预期放射治疗区域的第二子区域对应的第二物理性质，该第二子区域与所述射线相交于比所述第一子区域更接近所述源位置的位置处。

2.如权利要求1所述的系统，其中，计算所述第二子区域处的所述物理性质的所述计算取决于来自计算所述第一子区域处的所述物理性质的所述计算的结果。

3.如权利要求1所述的系统，其中所述多个射线跟踪计算是由所述并行处理器基本上同时执行的至少两个反向射线跟踪计算。

4.如权利要求1所述的系统，其中，所述第一物理性质和第二物理性质对应于表示从入射的放射束提取的能量的相对量的衰减因子。

5.如权利要求1所述的系统，其中，与具有所述预期放射治疗区域的所述身体相关的所述输入信息包括来自所述身体的经验数据。

6.如权利要求5所述的系统，其中，所述经验数据包括所述预期放射治疗区域的电子计算机X射线断层扫描技术(CT)图像。

7.如权利要求1所述的系统，其中，与具有所述预期放射治疗区域的所述身体相关的所述输入信息包括预编译的查找表。

8.如权利要求1所述的系统，其中，所述并行处理器还适于基于所述输入信息计算会由放射治疗的所述预期区域所吸收的放射能量的量。

9.如权利要求1所述的系统，其中，所述第一子区域和第二子区域在对应于数量基本相同的体素的尺寸和形状方面基本上相等。

10.如权利要求1所述的系统，其中，所述第一子区域比所述第二子区域更大。

11.如权利要求1所述的系统，进一步包括：放射治疗系统，其与所述放射计划系统进行通信，所述放射治疗系统包括放射源。

12.如权利要求1所述的系统，其中，所述放射源适用于调节强度的放射治疗，以使所述放射治疗符合所述身体的所述预期放射治疗区域的形状。

13.如权利要求1所述的系统，其中，所述放射源是具有至少1MeV的平均光子能量的高能光子源。

14.如权利要求1所述的系统，还包括诊断系统，所述诊断系统与所述放射计划系统进行通信，所述诊断系统包括传感器以从所述身体获得要输入到所述放射计划系统的经验数据。

15.如权利要求14所述的系统，其中，所述诊断系统是CT扫描仪。

16.如权利要求1所述的系统，其中，所述并行处理器是通用图形处理单元。

17.如权利要求1所述的系统，其中，所述身体是人类或者动物中的至少一个。

18.一种确定放射治疗参数的方法，包括：

获得关于身体的预期放射治疗区域的信息；

基于所述信息，在从所述身体到源位置的方向上计算沿射线的第一子区域和第二子区域处的物理性质，所述射线在所述源位置和所述身体之间行进并且与所述身体相交于所述预期放射治疗区域内，沿所述射线的所述第二子区域比沿所述射线的所述第一子区域更接近所述源位置；以及

基于所述计算确定用于为所述预期放射治疗区域提供放射治疗的所述放射治疗参数。

19.如权利要求18所述的方法，其中，所述放射治疗参数的所述确定还包括：计算所述预期放射治疗区域将吸收的放射能量的量。

20.如权利要求19所述的方法，其中，所述第一子区域和第二子区域在对应于数量基本相同的体素的尺寸和形状方面基本上相等。

21.如权利要求20所述的方法，其中，所述第一子区域比所述第二子区域大。

22.如权利要求18所述的方法，其中，所述第二子区域处的所述物理性质是基于所述第一子区域处计算出的所述物理性质来计算的。

23.如权利要求18所述的方法，其中，所述物理性质是表示从入射的放射束提取的能量的相对量的衰减因子。

24.如权利要求18所述的方法，其中，所述信息包括来自具有所述预期放射治疗区域的所述身体的经验数据。

25.如权利要求24所述的系统，其中，所述经验数据包括所述预期放射治疗区域的电子计算机X射线断层扫描技术(CT)图像数据。

26.如权利要求18所述的系统，其中，所述信息包括预编译的查找表。

27.一种计算机可读介质包含软件，当该软件由计算机执行时使所述计算机执行如权利要求18至26之一所述的方法。

Images