卷积 / 迭加剂量计算方法的、 使用图形处理单元加速的用 于放射治疗的实时剂量计算 相关申请的交叉引用
本申请要求于 2008 年 5 月 8 日提交的美国临时申请第 61/126936 号的优先权， 因 此通过引用将该申请的整体内容合并于此。
背景技术 1. 发明领域
本发明涉及放射治疗系统， 更具体而言， 涉及体外射束放射治疗系统。
2. 相关技术的讨论
放射治疗是医学使用放射来治疗恶性的细胞， 例如癌细胞。这种放射可以具有电 磁的形式， 例如高能光子， 或者具有微粒的形式， 例如电子、 质子、 中子或者阿尔法粒子。 迄今为止， 在现今实践中使用得最普遍的放射形式是高能光子。光子在人体组 织中的吸收是由放射的能量及相关组织的原子结构所确定的。在放射肿瘤学中采用的能 量的基本单位是电子伏特 (eV) ； 103eV ＝ 1keV， 106eV ＝ 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( 即 TE(r′ )) 被定义为将能量 E 的流量 ψE(r′ ) 用相对于水的密度 ρ(r′ ) 以及在点 r′处 的线性衰减 μE(r′ ) 加权， 如下述方程 1 所示。
线性衰减系数 μE(r ′ ) 也取决于原子材料。尽管康普顿散射在与放射疗法相 关的百万伏能量范围内占优势， 并且康普顿散射取决于电子密度而不是材料， 但是， 在临床 上， 对于正常人体组织， 在标准 CT 数量和密度之间存在分段线性关系。通常而言， 在此能量 范围内非康普顿的相互作用被视为可以忽略。根据下述方程 2， 能量 E 在点 r′处的流量 ψE(r′ ) 是由源焦点 s 和在方向 r′上的能量 E 的入射流量 ψE， 0(r′ ) 所决定的。
然后， 迭加 / 卷积算法通过剂量沉积内核传播上述能量， 以确定每个位置处的最 终剂量。为了允许剂量沉积内核根据组织不均匀性来实际地进行调整， 基于将迭加与卷积 区别开来的下述方程 3 使用点之间的放射线距离 dρ(r， r′ )， 。
根据下述方程 4， 点 r 处的剂量 D(r) 是由在 TERMA 体积上的积分并用取决于能量 的剂量沉积内核 KE 加权而计算出来的。标准压缩 (collapsed) 椎体内核是由在点和内核 轴之间的放射线距离和相对角 ω 来标示的， 并且缺少几何距离的平方的作用。
通常而言， 不计算如方程 4 中的每个体素 (voxel) 的单能贡献。 而是如基于方程 5 和 6 的方程 7 中一样， 选择射线角度 ω 和方向 v 的离散组合， 并且沿着使用的单个多能内 核对其进行积分。由于距离平方效应， 这由近似的指数内核衰减和任意单个远距离的体素 急剧减少的贡献得到证明。射线表现出的立体角量上的增加抵消了距离平方效应。射线方 向被选择来平衡几何形状和内核的能量因子。
传统地， 根据方程 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′ ) 分为谱权 wE 和净流量因子 ψE(r′ )， 并且定义衰减因子， 获得下述方程 8。
图 2 示出了传统的正向射线跟踪。图 3 示出了根据本发明一种实施方式的反向射 线跟踪。通过从每个 TERMA 体素背向源投射射线， 沿着路线聚集净衰减， 当到达病人的边界 时， 可以获得早期的射线终结， 从而提高性能和精确度。在计算上述定义的衰减因子时， 预 先编译的查找表可以用于数量 μE/ρ。例如， 基于放射治疗区域的以 Hounsfield 为单位的 CT 图像， 可以从查找表查找这个数量。谱权 wE 例如可以包含在用于放射治疗的放射源的模 型中， 并且每个谱箱可以针对给定的射线独立地衰减。
然后可以根据下述方程 9 通过将针对每个体素的净衰减乘以从源朝向体素的净 入射流量来计算 TERMA。
T(r′ ) ＝ ψ0(r′ )A(r′ ) (9)
图 4 示 出 使 用 市 售 的 Pinnacle 系 统 (Philips Radiation Oncology Systems Madison WI) 和本发明一种实施方式以分析法计算出的沿中心轴的 TERMA 能量， 其中每个 TERMA 能量以 10cm 的深度归一化。
尽管方程 9 是 O(n4) 的算法， 但由于与用于标准正向方法的 O(n-3) 不同， 每个线程 仅向其自身的体素写入， 因此避免了读写冲突。 该算法允许对结构存储器联合进行读访问， 从而大幅度地提高了存储性能。当仅改变流量野时， 只要在相互使用或者强度调节最优化 期间， 就可以在多种临床分辨率下使用标准正向和反向射线跟踪算法， 并且可以计算该计 算方法的高额衰减量。根据本发明的一些实施方式， 这种预先计算可以提供多个量级的性 能改进。10102036712 A CN 102036715
说大小 323 643 1283 2563 正向 23ms 44ms 150ms 1869ms 反向 2ms 18ms 220ms 3317ms明书近似的衰减 1ms 2ms 26ms 454ms 强度调节 1ms 1ms 3ms 13ms8/13 页表 1. 标准正向、 新的反向和近似的 TERMA 计算方法的性能。也包括在强度调节期 间有用的仅更新入射流量野的性能。
表 1 中比较了多种分辨率下基于衰减量的 TERMA 算法相对于传统射线投射算法的 性能。两种方法均物理结合了对多谱衰减的校正并且针对 GPU 进行了优化。在正向射线跟 踪计算中使用的射线投射可能有读 / 写冲突。射线投射的可量测性表现出最佳位置和所 需的参数调整， 以维持穿过每个体素的射线的必要数量。这种调整费力且很不直观。通过 该性能在其理论的 反向射线跟踪计算的基于衰减量的 TERMA 计算表现出经验性能 O(n3.76)， 4 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) 表示 从射线段到点的剂量沉积。
累积 - 累积内核 (CCK)(Lu， W.， Olivera， G.H.， Chen， M.， Reckwerdt， P.J.， Mackie， T.R.， Phys.Med.Biol.50， 655-680， 2005) 表示从射线段到射线段的剂量沉积。
对于特定的放射线深度 dr 和角度 ω， 这两者都源于标准的点到点内核 K(dr， ω) 的积分 ：
尽管更精确 ( 尤其在粗糙的分辨率下 )， 但 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(ωDT1/3) 减小到 O(ωDlog(T1/3)) 而提高了性能 ； 其中 ω 是角的数量， D 是剂量体素数量， 而 T 是 TERMA 体素 的数量。
与标准方法相比， 根据本发明一种实施方式的多分辨率迭加表现出令人感兴趣的 精确性折衷。表 2 将根据本发明一些实施方式的多分辨率迭加算法的精确性在多个野尺寸 和内核射线采样的情况下与标准方法相对比。由于很少有 TERMA 被射线几何地错过， 因此 根据本发明一些实施方式的多分辨率迭加算法通常在半影和用于小野尺寸的低剂量区域 中能够执行得更好。 然而， 由于大步长可以引起射束界限被模糊， 因此高剂量区域的精确性 可以轻微地折衷。使用相同的步长、 但不使用多分辨率数据结构的作为本发明的实施方式 的一种变型的多分辨率迭加算法的一种变型可表现出缓存性能下降以及可平均增加 60％ 的平均误差。
表 2. 平均沉积剂量相对于最大沉积剂量 Dmax 的误差， 其野尺寸从 1cm 到 23cm。通 过具有被限定为具有比 0.3Dmax 更大的剂量梯度的半影区域和比 0.2Dmax 更低的低剂量区域 的区域来减小误差。使用 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 的立方水假体上执行。试验在代表标准临床工作量 的具有 643 个体素的体积上、 和在具有 1283 个体素的额外的高分辨率体积上执行。作为本 发明一种实施方式的多分辨率迭加方法利用测定体积的 mip-map( 一种纹理映射技术 ) 并 且可以比传统的迭加执行起来快 2 倍到 3 倍以及如示出的那样的具有对高分辨率调节得更 好的性能。
图 12 示出根据本发明另一实施方式的一种方法。本方法包括 ： 获得关于身体的 预期放射治疗区域的输入信息 1201 ； 在方框 1202 中， 在从所述身体到源位置的方向上， 基 于输入信息 1201 计算沿射线的第一和第二子区域的物理性质， 射线从源位置出发穿过身 体并且与身体交叉于预期放射治疗区域， 沿射线的第二子区域比沿射线的第一子区域更接 近于源位置 ； 以及在方框 1203 中， 确定用于为预期放射治疗区域提供放射治疗的输出信息 1204。
在方框 1202 中， 基于在第一子区域处计算的物理性质计算第二子区域处的物理 性质。物理性质是指表示从入射的放射束提取的能量的相对数量的衰减因子。如上讨论的
衰减因子是一个实例。此外， 输入信息 1201 可以包括来自身体的经验数据。测量数据可以 包括身体的预期治疗区域的至少一个电子计算机 X 射线断层扫描技术 (CT) 图像。此外， 在 方框 1202 中可以使用至少一个预编译的查找表。
方框 1203 还包括计算被预期放射治疗区域吸收的放射能量的量的步骤 1204。使 用统一采样、 可变采样、 多分辨率网格或者其变型中的至少一个以覆盖预期放射治疗区域， 可以计算出被吸收的放射能量。此外， 多分辨率网格被指定使得来自源位置穿过身体行进 的射线可以近似于立体角。
然而， 本发明的另一实施方式可以包括计算机可读介质， 其包括软件， 当该软件被 计算机执行时使计算机执行上述方法。
在描述本发明的实施方式时， 为了清楚起见而采用了特定的术语。 然而， 本发明并 非意在受限于由此而选择的特定术语。 本领域的技术人员可想到在不偏离本发明的情况下 可根据上述构思对本发明的上述实施方式做出修改或变型。因此可以理解， 可以在所附权 利要求和与其等同方案的范围内以与所明确描述的方式不同的方式来实施本发明。