一种测算放疗射线剂量分布和剂量目标函数的方法和系统与流程

本发明涉一种测算放疗射线剂量分布和剂量目标函数的方法和系统，属于医疗器械技术领域，具体地说是一种利用三维医学影像和分布式并行计算测算放疗射线剂量分布和剂量目标函数的方法和系统。

背景技术：

放射治疗计划系统(treatmentplanningsystem，tps)是放疗科物理师或者剂量师用来设计放疗计划的技术平台和重要软件类医疗器械。物理师或者剂量师将病人的三维影像和医生剂量处方输入tps，然后根据经验勾画靶区、重要器官和辅助器官、放置最佳射野、设置各种优化目标函数、最后进行计划的优化计算，计算出最合适的射线束分布及强度，从而得到放疗计划。如果计算结果不能满足处方要求，操作人员将进一步调整各种参数继续优化，直到得到满意结果为止。其中，测算放疗射线剂量分布从而计算剂量目标函数是评价放射治疗计划的重要步骤。此种评价被广泛应用于放射治疗计划的正向设计，以及放射治疗计划的逆向优化和质量保证上，是tps的支撑性技术和模块。

目前世界上通用的tps的优化技术采用的是基于多目标的非约束型优化技术。该技术使用现代优化算法，如cg(conjugategradient)算法或者遗传算法，通过最小化一个根据不同权重综合了多种相互竞争的临床目标和约束条件的成本函数来搜寻最优的调强放射治疗计划。这个基于多目标的非约束型优化问题可由以下公式来表达：

其中，p＝{pj，j＝1，...，ntarget，pk，k＝1，...，noar}为上述优化问题的参数集合，它们代表了各种剂量分布限制和各优化目标的权重；ntarget和noar分别为计划目标数目和重要器官数目，i为射线束强度分布，d(i)为根据射线束的强度分布得到的三维剂量分布，f为靶区(target)或者危机器官(oar)的剂量目标函数，i^*为最优的射线束强度分布。由于只有一个成本函数，该算法与传统的约束型优化技术相比具有计算速度较为快速的优点。如果与肿瘤放射医师的临床经验结合，该算法可以产生可行的治疗计划。可以看到，在上述计划优化过程中，每调整一次射线束强度分布，均要计算一次射线剂量分布和剂量目标函数。由此，射线剂量分布d(i)和剂量目标函数的计算效率对于优化的质量和效率至关重要。目前的商用tps一般先把直线加速器治疗机头的出射平面做网格化细分，然后采用实验数据插值法或者基于加速器机头模型的蒙特卡洛算法或者卷积叠加算法计算每一个网格出束在人体或者模体内的三维剂量沉积分布(即笔形束剂量沉积矩阵bi)，然后再根据出射平面上的光子通量调制分布对笔形束剂量沉积矩阵进行加权求和得到射线剂量分布d(i)。具体的计算公式如下：

其中，ii是光子通量平面上每个网格上的权重(即笔形束的权重)，k为笔形束的数量。由于在由笔型束叠加形成三维剂量分布的计算中每个笔形束剂量沉积矩阵bi均是三维矩阵(由于照射人体或者模体一般均被网格化为三维离散分布)，假设该矩阵的大小为l×m×n，则存储和计算的复杂度均为o(l×m×n×k)。对于一般的离散化剂量分布，l、m、n在100量级，而k在1000量级，由此得到的总存储和计算的复杂度均达到了10⁹量级。随着imrt和vmat等复杂放疗调强技术和基于人工智能的计划方法的应用，笔形束的数量、剂量分布的空间的精度不断提高和优化迭代次数随着精度和要求的提高而快速增加，射线剂量分布以及相应的剂量目标函数的计算所需的内存空间和时间还要迅速地指数级增长。此种增长给放射治疗计划系统性能的进一步提升带来困扰。如何提高射线剂量分布和剂量目标函数的计算效率和存储效率成为了放疗技术的新课题。

技术实现要素：

本发明所要解决的技术问题是针对放射治疗计划系统(tps)的上述技术现状，而提供一种测算放疗射线剂量分布和剂量目标函数的方法和系统。

本发明解决上述技术问题所采用的技术方案作为一种测算放疗射线剂量分布和剂量目标函数的方法，包括以下步骤：

步骤1：从影像设备通过dicom协议接收三维医学影像数据；

步骤2：根据三维医学影像数据分割感兴趣区；

步骤3：根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵；

步骤4：将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵，子区域和子矩阵的数量为可调用的从计算节点的数量；

步骤5：编制区域节点表，为每一个子区域对应唯一一个从计算节点；

步骤6：根据区域节点表，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵和感兴趣区几何数据通过数据网络传递到该子区域所对应的从计算节点上；

步骤7：根据区域节点表，将每一个子区域所对应的笔形束权重通过数据网络传递到该子区域所对应的从计算节点上；

步骤8：在每一个从计算节点上将所有笔形束的三维剂量沉积子矩阵进行矩阵元对矩阵元加权叠加，每个子矩阵的权重为该笔形束的射线通量，此步计算在所有从计算节点上的并行展开；

步骤9：在每一个从计算节点上计算每个危机器官的剂量直方曲线，此步计算在所有从计算节点上的并行展开；

步骤10：每一个从计算节点将所计算的剂量体积直方曲线传递到主计算节点；

步骤11：在主计算节点上分别将每一个感兴趣区的子区域剂量直方曲线合并计算成为该感兴趣区的全局剂量直方曲线；

步骤12：在主计算节点上根据每一个感兴趣区的剂量体积直方曲线计算该感兴趣区所关联的剂量目标函数。

上述的放疗射线包括光子线、电子线和质子线。

上述的从计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的主计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。

上述的剂量目标函数包括并不限于：最大剂量、最小剂量、平均剂量、最大剂量体积、最小剂量体积、最大等效均一剂量和最小等效均一剂量等。

本发明解决上述技术问题所采用的技术方案作为一种测算放疗射线剂量分布和剂量目标函数的系统，包括以下部分：

一个主计算节点，负责从影像设备通过dicom协议接收三维医学影像数据；根据三维医学影像数据分割感兴趣区；根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵；将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵，子区域和子矩阵的数量为可调用的从计算节点的数量；编制区域节点表，为每一个子区域对应唯一一个从计算节点；根据区域节点表，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵和感兴趣区几何数据通过数据网络传递到该子区域所对应的从计算节点上；根据区域节点表，将每一个子区域所对应的笔形束权重通过数据网络传递到该子区域所对应的从计算节点上；在从计算节点完成计算后，该主计算节点负责接收从每一个从计算节点所计算的局部剂量体积直方曲线，分别将每一个感兴趣区的局部剂量直方曲线合并计算成为该感兴趣区的全局剂量直方曲线，并根据每一个感兴趣区的全局剂量体积直方曲线计算该感兴趣区所关联的剂量目标函数；

若干个从计算节点，其中在每一个从计算节点上将所有笔形束的三维剂量沉积子矩阵进行矩阵元对矩阵元加权叠加得到该子区域内的剂量分布，每个子矩阵的权重为该笔形束的射线通量，此步计算在所有从计算节点上的并行展开；在每一个从计算节点上计算每个危机器官的局部剂量直方曲线，此步计算在所有从计算节点上的并行展开；

其中，主计算节点和所有从计算节点之间通过数据网络进行数据通讯。

上述的放疗射线包括光子线、电子线和质子线。

上述的从计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的主计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。

上述的剂量目标函数包括并不限于：最大剂量、最小剂量、平均剂量、最大剂量体积、最小剂量体积、最大等效均一剂量和最小等效均一剂量等。

与现有技术相比，这种测算放疗射线剂量分布和剂量目标函数的方法和系统弥补了现有放射治疗计划系统的在剂量分布和剂量目标函数的计算效率方面的不足。该方法通过将放射剂量计算区域分割成不重叠的若干子区域，并将每个区域分别放置在不同cpu或者gpu上进行平行计算，有效地提高了计算效率和存储效率。同时，这种计算方法和系统中的值得注意的一个优势为：为了计算tps所需的最终目标，即与优化目标相关的感兴趣区的目标函数，各个分布式计算节点并不向主计算节点传送每次优化迭代过程中产生的中间三维剂量分布数据，而是通过在每个cpu或gpu上并行计算该分布式计算节点所对应剂量区域的剂量分布所对应的感兴趣区的剂量体积直方曲线，然后将由一维数据组成的剂量体积直方曲线传送到主计算节点进行剂量体积直方曲线融合和归一化，然后通过融合后的剂量体积直方曲线计算总体的目标函数。这种仅向主计算节点传递少量综合性一维数据的方法有效减少了主从节点之间的数据通讯量，消除了并行计算中可能存在的潜在数据带宽瓶颈。

附图说明

图1是本发明实施例中的部署结构图；

图2是本发明实施例中步骤1-11的逻辑流程图；

图3是本发明实施例中步骤10中的子步骤的逻辑流程图。

具体实施方式

以下结合附图对本发明的实施例作进一步详细描述。

本发明的一种测算放疗射线剂量分布和剂量目标函数的方法和系统的实施例的特征是：包括以下步骤：

步骤1：从影像设备通过dicom协议接收三维医学影像数据；

步骤2：根据三维医学影像数据分割感兴趣区；

步骤3：根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵；

步骤4：将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵，子区域和子矩阵的数量为可调用的从计算节点的数量；

步骤5：编制区域节点表，为每一个子区域对应唯一一个从计算节点；

步骤6：根据区域节点表，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵和感兴趣区几何数据通过数据网络传递到该子区域所对应的从计算节点上；

步骤7：根据区域节点表，将笔形束权重通过数据网络传递到该子区域所对应的从计算节点上；

步骤8：在每一个从计算节点上将所有笔形束的三维剂量沉积子矩阵进行矩阵元对矩阵元加权叠加，每个子矩阵的权重为该笔形束的射线通量，此步计算在所有从计算节点上的并行展开；

步骤9：在每一个从计算节点上计算每个危机器官的剂量直方曲线，此步计算在所有从计算节点上的并行展升；

步骤10：每一个从计算节点将所计算的剂量体积直方曲线传递到主计算节点；

步骤11：在主计算节点上分别将每一个感兴趣区的子区域剂量直方曲线合并计算成为该roi的全局剂量直方曲线；

步骤12：在主计算节点上根据每一个感兴趣区的剂量体积直方曲线计算该感兴趣区所关联的剂量目标函数。

上述的放疗射线包括光子线、电子线和质子线。

上述的从计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的主计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。

上述的剂量目标函数包括并不限于：最大剂量、最小剂量、平均剂量、最大剂量体积、最小剂量体积、最大等效均一剂量和最小等效均一剂量等。

本发明的一种计算放疗射线剂量分布和剂量目标函数的系统，包括以下部分：

一个主计算节点，负责从影像设备通过dicom协议接收三维医学影像数据；根据三维医学影像数据分割感兴趣区；根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵；将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵，子区域和子矩阵的数量为可调用的从计算节点的数量；编制区域节点表，为每一个子区域对应唯一一个从计算节点；根据区域节点表，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵和感兴趣区几何数据通过数据网络传递到该子区域所对应的从计算节点上；根据区域节点表，将笔形束权重通过数据网络传递到该子区域所对应的从计算节点上；在从计算节点完成计算后，该主计算节点负责接收从每一个从计算节点所计算的局部剂量体积直方曲线，分别将每一个感兴趣区的局部剂量直方曲线合并计算成为该感兴趣区的全局剂量直方曲线，并根据每一个感兴趣区的全局剂量体积直方曲线计算该感兴趣区所关联的剂量目标函数；

若干个从计算节点，其中在每一个从计算节点上将所有笔形束的三维剂量沉积子矩阵进行矩阵元对矩阵元加权叠加得到该子区域内的剂量分布，每个子矩阵的权重为该笔形束的射线通量，此步计算在所有从计算节点上的并行展开；在每一个从计算节点上计算每个危机器官的局部剂量直方曲线，此步计算在所有从计算节点上的并行展开；

其中，主计算节点和所有从计算节点之间通过数据网络进行数据通讯。

上述的放疗射线包括光子线、电子线和质子线。

上述的从计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的主计算节点可以为cpu、gpu或者cpu与gpu的组合。

上述的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。

上述的剂量目标函数包括并不限于：最大剂量、最小剂量、平均剂量、最大剂量体积、最小剂量体积、最大等效均一剂量和最小等效均一剂量等。

图1是本发明实施例中的部署结构图。装有cpu/gpu的主计算节点101负责从影像设备102通过数据网络112采用dicom协议接收三维医学影像数据103；根据三维医学影像数据分割感兴趣区；根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵104；主计算节点101将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵105，子区域和子矩阵的数量为可调用的从计算节点107的数量；然后，编制区域节点表106，为每一个子区域对应唯一一个从计算节点上107，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵105和感兴趣区几何数据106通过数据网络113传递到该子区域所对应的从计算节点109上；将每一个子区域所对应的所有笔形束的权重107通过数据网络113传递到该子区域所对应的从计算节点109上；在每一个从计算节点109上将所有笔形束的三维剂量沉积子矩阵105进行矩阵元对矩阵元加权叠加，每个子矩阵的权重为该笔形束的射线通量，此步计算在所有从计算节点109上的并行展开；在每一个从计算节点109上计算每个危机器官的局部剂量直方曲线110，此步计算在所有从计算节点上109的并行展开；每一个从计算节点将所计算的剂量体积直方曲线传递到主计算节点101；在主计算节点101上分别将每一个感兴趣区的局部剂量直方曲线110合并计算成为该roi的全局剂量直方曲线111；在主计算节点上根据每一个感兴趣区的全局剂量体积直方曲线111计算该感兴趣区所关联的剂量目标函数112。从图1中可以看出，由于剂量目标函数112仅仅依赖于全局剂量体积直方曲线111，因此为了获得最终优化计划需要的剂量目标函数112，无需将各子区域的剂量传递至主计算节点101，只需传递各子区域的局部剂量直方曲线110，通过在主计算节点101合成全局剂量直方曲线111即可。而局部剂量直方曲线110是经过对三维剂量分布进行综合统计得到的一维数据，所需传输的数据量极少，从而避免了分布式并行计算中的数据带宽瓶颈。在一般实施例中，笔形束的三维剂量沉积子矩阵105和感兴趣区几何数据106一般不变，只需向各从计算节点传送一次即可，因此虽然三维剂量沉积子矩阵有一定的数据量，但是不会给网络造成瓶颈。所有笔形束的权重107在计划评估实施例中，也只分发一次即可。而在计划优化的实施例中，笔形束的权重107就需要被分发多次，剂量目标函数也需要多次计算，在这种情况下，由于只传递一维局部剂量直方曲线数据110而不再传递三维剂量分布而带来的传输高效性就更加显著了。

图2是本发明实施例中步骤1-11的逻辑流程图。从影像设备通过dicom协议接收三维医学影像数据步骤201可以通过任何一种高级计算机编程语言(如c++，java，python等)和主计算节点上的计算机操作系统数据通讯接口(如socket或者webservice)实现。根据三维医学影像数据分割感兴趣区步骤202可以运用已有的影像编辑软件工具或者软件算法互动式地或者自动地实现。根据三维医学影像数据计算医用加速器治疗机头射野开口平面上的所有笔形束的三维剂量沉积矩阵步骤203可以通过已有的蒙特卡洛算法或者卷积叠加方法实现。将放射剂量计算区域分割成不重叠的若干子区域，将每一个笔形束的三维剂量沉积矩阵依据子区域分割成若干子矩阵步骤204可以通过任何一种高级计算机编程语言实现。编制区域节点表，为每一个子区域对应唯一一个从计算节点步骤205可以通过任何一种高级计算机编程语言实现。根据区域节点表，将每一个子区域所对应的所有笔形束的三维剂量沉积子矩阵和感兴趣区几何数据通过数据网络传递到该子区域所对应的从计算节点上步骤206，以及根据区域节点表，将笔形束权重通过数据网络传递到该子区域所对应的从计算节点上207可以通过任何一种高级计算机编程语言和并行计算库(如mpi，massivelyparallelinterface)里的数据广播通讯接口实现。在每一个从计算节点上将所有笔形束的三维剂量沉积子矩阵进行矩阵元对矩阵元加权叠加得到该子区域内的剂量分布，每个子矩阵的权重为该笔形束的射线通量步骤208；以及在每一个从计算节点上计算每个危机器官的剂量直方曲线步骤209可以通过任何一种高级计算机编程语言结合通用并行计算库(如mpi)实现。每一个从计算节点将所计算的剂量体积直方曲线传递到主计算节点步骤210可以通过任何一种高级计算机编程语言和通用并行计算库(如mpi)里的数据通讯接口实现。在主计算节点上分别将每一个感兴趣区的局部剂量直方曲线合并计算成为该roi的全局剂量直方曲线步骤211；以及在主计算节点上根据每一个感兴趣区的剂量体积直方曲线计算该感兴趣区所关联的剂量目标函数步骤212可以通过任何一种高级计算机编程语言实现。在本实施流程的步骤中涉及从计算节点上的并行计算部分，可以采用cpu上的多进程或者多线程实现，也可以采用gpu上并行计算单元实现。例如若采用nvdia的gpu，基础控制gpu的计算函数库可以采用cuda实现。

图3是本发明实施例中步骤10中实施剂量体积直方图融合计算的逻辑流程图。该子流程中的所有子步骤均可以由高级编程语言通过系统调用实现。本流程利用了各计算区域计算的局部剂量体积直方曲线的线性可加性。所有局部微分剂量体积直方曲线301依据公式302叠加起来得到全局微分剂量体积直方曲线303，再以及公式304作归一化得到归一化的全局微分剂量体积直方曲线305，最后用公司306计算出归一化的全局积分剂量体积直方曲线307。其中，求和限s为剂量离散化后的离散区间的最高标号。

本发明所涉及的已有的专业名词解释定义如下：

调强放射治疗调强放射治疗(intensity-modulatedradiationtherapy，imrt)通常采用治疗计划预先确定的连续或离散的方式，调控光子或电子束的注量、相对于患者的射束方向和射野尺寸。imrt的主要作用是提高剂量分布对计划靶区的适形度，同时使剂量尽量规避周围正常组织。

dicom医学数字成像和通信标准，医学数字成像和通信(digitalimagingandcommunicationsinmedicine，缩写dicom)标准是医学图像和相关信息的国际标准。这个标准定义了质量能满足临床需要的可用于数据交换的医学图像格式。当医学影像数据按照dicom标准进行传输和储存时，其数据文件中包含的病人信息可按照dicom标准文件解析算法解析为一系列二维影像，可供医生阅片。

感兴趣区感兴趣区(regionofinterest，缩写roi)是从三维影像中选择的一个特定三维区域，在肿瘤放射治疗中通常代表肿瘤靶区或者正常组织，通常由各层影像切片上的闭合轮廓线集合来表示。这些区域是影像分析关注的重点，在放疗中一般需要计算它们各自的剂量目标函数。

靶区靶区(targetvolume)是指放射治疗中，准备向患者体内以治疗为目的的辐照一定吸收剂量的区域。

剂量目标函数剂量目标函数是指某个射线剂量的函数，该函数反映某一个感兴趣区中的剂量分布之特征，例如：最大剂量(maxdose)、最小剂量(mindose)、平均剂量(meandose)、最大剂量体积(maxdvh)、最小剂量体积(mindvh)、最大等效均一剂量(maxeud)和最小等效均一剂量(mineud)等。

等效均一剂量等效均一剂量(equivalentuniformdose，缩写eud)，指一个不均匀剂量分布的生物等效均匀剂量。计算公式如下：

其中，n是感兴趣区解剖结构三维格点的数量；di是第i个三维格点的剂量；a是反映剂量体积生物效应的参数，不同感兴趣区(肿瘤或正常组织)有不同的取值。

剂量体积直方曲线剂量体积直方曲线(dosevolumehistogram，缩写dvh)的横轴为剂量，纵轴为体积百分比；有微分曲线和积分曲线两种。微分剂量体积直方曲线反映在某个微小剂量区间内的剂量在某感兴趣区的百分体积占比；积分剂量体积直方曲线反映在大于某个值的剂量在某感兴趣区的百分体积占比。它们均为综合反映剂量在某个感兴趣区内分布特征的剂量学量，由它们可以计算得到各种剂量目标函数。

本发明的最佳实施例已阐明，由本领域普通技术人员做出的各种变化或改型都不会脱离本发明的范围。