CN112017166A - 一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 - Google Patents
一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 Download PDFInfo
- Publication number
- CN112017166A CN112017166A CN202010853607.XA CN202010853607A CN112017166A CN 112017166 A CN112017166 A CN 112017166A CN 202010853607 A CN202010853607 A CN 202010853607A CN 112017166 A CN112017166 A CN 112017166A
- Authority
- CN
- China
- Prior art keywords
- gpu
- dimensional
- cpu
- convolution
- dose distribution
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000009826 distribution Methods 0.000 title claims abstract description 67
- 238000001959 radiotherapy Methods 0.000 title claims abstract description 47
- 238000000034 method Methods 0.000 title claims abstract description 42
- 230000005855 radiation Effects 0.000 claims description 22
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 230000015654 memory Effects 0.000 claims description 16
- 239000007787 solid Substances 0.000 claims description 7
- 238000004891 communication Methods 0.000 claims description 6
- 239000003814 drug Substances 0.000 claims description 4
- 238000003384 imaging method Methods 0.000 claims description 4
- 238000004804 winding Methods 0.000 claims description 4
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 3
- 230000002093 peripheral effect Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 22
- 230000008569 process Effects 0.000 abstract description 7
- 238000012545 processing Methods 0.000 abstract description 4
- 238000010586 diagram Methods 0.000 description 11
- 238000005457 optimization Methods 0.000 description 10
- 230000006870 function Effects 0.000 description 8
- 238000002721 intensity-modulated radiation therapy Methods 0.000 description 6
- 230000008021 deposition Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 230000004907 flux Effects 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 210000000056 organ Anatomy 0.000 description 4
- 206010028980 Neoplasm Diseases 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000001678 irradiating effect Effects 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- XCWPUUGSGHNIDZ-UHFFFAOYSA-N Oxypertine Chemical compound C1=2C=C(OC)C(OC)=CC=2NC(C)=C1CCN(CC1)CCN1C1=CC=CC=C1 XCWPUUGSGHNIDZ-UHFFFAOYSA-N 0.000 description 1
- 231100000987 absorbed dose Toxicity 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 238000010894 electron beam technology Methods 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009290 primary effect Effects 0.000 description 1
- 230000002285 radioactive effect Effects 0.000 description 1
- 230000001225 therapeutic effect Effects 0.000 description 1
- 238000002560 therapeutic procedure Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/103—Treatment planning systems
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/1042—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy with spatial modulation of the radiation beam within the treatment head
- A61N5/1045—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy with spatial modulation of the radiation beam within the treatment head using a multi-leaf collimator, e.g. for intensity modulated radiation therapy or IMRT
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/1048—Monitoring, verifying, controlling systems and methods
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/1077—Beam delivery systems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
- G06F9/38—Concurrent instruction execution, e.g. pipeline or look ahead
- G06F9/3885—Concurrent instruction execution, e.g. pipeline or look ahead using a plurality of independent parallel functional units
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T1/00—General purpose image data processing
- G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H30/00—ICT specially adapted for the handling or processing of medical images
- G16H30/20—ICT specially adapted for the handling or processing of medical images for handling medical images, e.g. DICOM, HL7 or PACS
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20092—Interactive image processing based on input by user
- G06T2207/20104—Interactive definition of region of interest [ROI]
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- General Health & Medical Sciences (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Theoretical Computer Science (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Pathology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Animal Behavior & Ethology (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- Software Systems (AREA)
- Medical Informatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Primary Health Care (AREA)
- General Engineering & Computer Science (AREA)
- Epidemiology (AREA)
- Quality & Reliability (AREA)
- Radiation-Therapy Devices (AREA)
Abstract
本发明公开了一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，属于医疗器械技术领域。该方法通过卷积叠加方法计算笔形束剂量分布过程分解为多个一维卷积，并对每个一维卷积进行预采样，判断其计算量，并根据计算量对一维卷积进行归类，在分批并行处理多个一维卷积时，保证每个批次的一维卷积过程都具有相同的计算量，从而保证了每批并行处理线程的指令同步性，大大发挥了GPU的并行效率。同时，本方法将未能分配的具有不同计算量的剩余一维卷积过程分配至CPU进行多线程并行计算。由于CPU的线程没有对指令同步性设置要求，该部分计算可以高效的并行执行。而GPU与CPU计算的同时性进一步提高了混合集群的并行度和计算效率。
Description
技术领域
本发明涉一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，属于医疗器械技术领域，具体地说是一种利用CPU/GPU混合集群分布式并行计算放疗笔形束射线剂量分布的方法。
背景技术
放射治疗计划系统(Treatment Planning System，TPS)是放疗科物理师或者剂量师用来设计放疗计划的技术平台和重要软件类医疗器械。物理师或者剂量师将病人的三维影像和医生剂量处方输入TPS，然后根据经验勾画靶区、重要器官和辅助器官、放置最佳射野、设置各种优化目标函数、最后进行计划的优化计算，计算出最合适的射线束通量的强度分布，从而得到放疗计划。如果计算结果不能满足处方要求，操作人员将进一步调整各种参数继续优化，直到得到满意结果为止。目前世界上通用的TPS的优化技术采用的是基于多目标的非约束型优化技术。该技术使用现代优化算法，如准牛顿算法或者遗传算法，通过最小化一个根据不同权重综合了多种相互竞争的临床目标和约束条件的成本函数来搜寻最优的调强放射治疗计划。这个基于多目标的非约束型优化问题可由以下公式来表达：
其中，P＝{Pj，j＝1，...，Nt arg et，Pk，k＝1，...，NOAR}为上述优化问题的参数集合，它们代表了各种剂量分布限制和各优化目标的权重；Ntarget和NOAR分别为计划目标数目和重要器官数目，I为射线束强度分布，D(I)为根据射线束的强度分布得到的三维剂量分布，F为靶区(target)或者危机器官(OAR)的剂量目标函数，I*为最优的射线束强度分布。由于只有一个成本函数，该算法与传统的约束型优化技术相比具有计算速度较为快速的优点。如果与肿瘤放射医师的临床经验结合，该算法可以产生可行的治疗计划。目前的商用TPS一般先把直线加速器治疗机头的出射平面做网格化细分，然后采用实验数据插值法或者基于加速器机头模型的蒙特卡洛算法或者卷积叠加算法计算每一个网格出束在人体或者模体内的三维剂量沉积分布(即笔形束剂量沉积矩阵)，然后再根据出射平面上的光子通量调制分布对笔形束剂量沉积矩阵进行加权求和得到射线剂量分布D(I)。具体的计算公式如下：
其中，Ii是光子通量平面上每个网格上的权重(即笔形束的权重)，K为笔形束的数量。由于在由笔型束叠加形成三维剂量分布的计算中每个笔形束剂量沉积矩阵Bi均是三维矩阵(由于照射人体或者模体一般均被网格化为三维离散分布)，假设该矩阵的大小为L×M×N，则存储和计算的复杂度均为O(L×M×N×K)。对于一般的离散化剂量分布，L、M、N在100量级，而K在1000量级，由此得到的总存储和计算的复杂度均达到了109量级。随着IMRT和VMAT等复杂放疗调强技术和基于人工智能的计划方法的应用，笔形束的数量、剂量分布的空间的精度不断提高和优化迭代次数随着精度和要求的提高而快速增加，从而计算射线笔形束剂量分布所需的内存空间和时间还要迅速地指数级增长。此种增长给放射治疗计划系统性能的进一步提升带来困扰。
利用CPU集群进行并行计算是一种通用的利用并行计算提高射线笔形束剂量分布的计算效率和存储效率的方法。随着现代图形加速卡GPU的发展，基于通用图形加速卡(GPGPU)的并行计算逐渐发展成为一种低成本的拥有更高并行度的加速笔形束剂量分布的计算方法。然而，GPU的并行计算线程对于线程指令的同步性有很高的要求。如果并行计算的各线程不能做到指令同步，线程间互相等待的机制将大大降低并行效率。另外，现代CPU的运算速度也有了很大发展，由于CPU上的多线程计算没有GPU那样的对指令同步的要求，在GPU计算的同时闲置CPU是对计算资源的浪费。由此，如何同时将GPU/CPU在混合集群上利用起来进行并行计算，并发挥各自的优势弥补各自的劣势，成为计算放疗笔形束射线剂量分布的新问题。
发明内容
本发明所要解决的技术问题是针对放射治疗计划系统(TPS)的上述技术现状，而提供一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，包括以下步骤：
步骤1：从影像设备通过医学数字成像和通信标准(DICOM)协议接收人体或模体的三维医学影像数据；
步骤2：根据三维医学影像数据勾画和分割感兴趣区，并计算包含所有感兴趣区的边界盒(Bounding Box)。
步骤3：将放射治疗机头的每个射野的准直器围成的矩形分割成一系列网格，将通过每个网格的细小射线束记为笔形束Bn，n为笔形束编号；
步骤4：将所有笔形束均匀分配到混合集群计算节点上，并通过数据网络将笔形束的几何信息传递到该笔形束所对应的计算节点上，并为笔形束的剂量分布Dn(x，y，z)分配内存空间和显存空间，x，y，z为三维坐标点。
步骤5：利用集群节点的GPU通过光线追踪法计算每个笔形束Bn在模拟数字人体中的TERMA(单位质量总释放能量)Tn(x，y，z)，x，y，z为三维坐标点；
步骤6：将4π立体角分解为m个角度，用(n，m，i，j)元组表示的Bn的第m个方向的第(i，j)条卷积线，其中(i，j)为Bn在第m个方向上的卷积线的二维离散标识，并将该卷积线分配到集群的某一个GPU的某个流处理器线程上计算该卷积线上的第一个非零TERMA三维空间点的位置，以及从该位置沿该方向到与边界盒交点的线段所相交的三维空间离散腔胞数量Lnmij，并将它记为该卷积线的离散长度；
步骤7：计算Lnmij的的分布f(l)，根据该分布分批将具有相同离散长度的卷积线分配到集群的GPU的流处理器线程上并行计算沿该卷积线的一维卷积，
并将得到的三维剂量分布Egpu n(x，y，z)叠加到Dgpu n(x，y，z)；
步骤8：将未分配到GPU的所有剩余卷积线均匀分配到集群CPU的线程上并行计算该方向上的一维卷积，
并将得到的三维剂量分布Ecpu n(x，y，z)叠加到Dcpu n(x，y，z)；
步骤9：将每个GPU和CPU计算的笔形束剂量分布叠加起来得到最终的笔形束剂量分布：Dn(x，y，z)＝Dgpu n(x，y，z)+Dcpu n(x，y，z)。
上述的放疗射线包括光子线、电子线和质子线。
上述的集群计算节点至少包含一个或若干个CPU、一个或若干个GPU以及CPU对应的
内存和GPU对应的显存，而且集群计算节点之间通过数据网络连接；
上述的的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。
上述的的Bn在模拟数字人体中的单位质量总释放能量Tn(x，y，z)包括不同能级的TERMA的线性叠加，其中线性组合系数由放疗射线的能谱决定。
上述的的步骤7的分配算法至少包括如下子步骤：
步骤1：计算Lnmij的范围区间[lmin，lmax]，其中lmin为最小长度值，lmax为最大长度值；
步骤2：遍历区间[lmin，lmax]，对于该区间中的每一个长度l，计算f(l)/k的整数商s，其中f(l)为长度为l的卷积线的数量，k为每个GPU的流处理器数量；
步骤3：对于该区间中的每一个长度l，如果相应的s大于零，分配s×k条卷积线到s×k个GPU线程；
上述的的步骤7中的一维卷积算法的一维卷积核Km至少包括如下子步骤：
步骤1：计算蒙特卡罗方法模拟生成放疗射线点源在等效水中的三维剂量分布K3(x，y，z)作为三维卷积核；
步骤2：在第m个立体角方向对三维卷积核K3(x，y，z)进行插值采样得到Km(r)；
上述的的步骤7和步骤8可以在GPU和CPU上分别同时进行。
与现有技术相比，这种计算放疗笔形束射线剂量分布弥补了现有放射治疗计划系统的笔形束剂量分布的计算效率方面的不足。该方法通过卷积叠加方法计算笔形束剂量分布过程分解为多个一维卷积，并对每个一维卷积进行预采样，判断其计算量，并根据计算量对一维卷积进行归类，在分批并行处理多个一维卷积时，保证每个批次的一维卷积过程都具有相同的计算量，从而保证了每批并行处理线程的指令同步性，大大发挥了GPU的并行效率。同时，本方法将未能分配的具有不同计算量的剩余一维卷积过程分配至CPU进行多线程并行计算。由于CPU的线程没有对指令同步性设置要求，该部分计算可以高效的并行执行。而且，由于GPU与CPU计算的同时性进一步提高了混合集群的并行度和计算效率。
附图说明
图1是本发明实施例中的设备示意图；
图2是本发明实施例中的部署结构图；
图3是本发明实施例中步骤1-9的逻辑流程图；
图4是本发明实施例中步骤7中的的分配算法子步骤的逻辑流程图。
图5是本发明实施例中步骤7中的子步骤中的卷积线离散长度的分布统计图。
具体实施方式
下面结合附图，应用一组术语，来对本发明较佳的实施例进行详细阐述，以方便本领域技术人员理解本发明的优点和特征，从而对本发明的保护范围作出更清楚明确的界定。
本发明包含的创新概念具有多种实施例。因此不应该以这里关于较佳的实施例进行的详细阐述，作为本发明关于创新概念提出的权利要求的边界，而应当将这里的说明用于帮助本领域专业人员理解本发明所包含的创新概念。此外，示意图中物体的层面和部位的大小和相对大小会做适当变形以免重叠。
这里会以引用序号的形式来在示意图中指出一段说明的描述对象。但是，在本发明的实施例的示意图中，并非所有的组成部分都会编号。原因包括：1)相关领域中的公开的信息不会在这里做详细描述；2)与上下文重复的部分不会赘述。
这里在描述本发明时，会使用“包括”来引用本发明所含的对象。如果在上下文中没有明确说明，则这种引用是未尽的，表示描述中省略了部分信息，并且认为这种省略不会影响本领域专业人员对本发明的方法的理解。
这里的详细描述会使用一些带有先后关系的词汇来方便罗列需要描述的对象。这些词汇不应当被看作是对本发明的方法的组成和结构做出的限制，而只应当看作是为了方便区分所叙内容而暂时添加的标记，调换所叙内容的先后顺序，例如，将“第一项”与“第二项”这两个词组交换位置，不会影响这里对本发明的创新概念的描述。类似的，当这里使用“和/或”这样的词汇来连接一组陈述时，并不对这些陈述的对象的先后顺序的组合提出限制。任意改变这些对象的陈述顺序都是可以接受的。
除非这里作出明确的说明来赋予指定术语以独特的定义，否则在详细描述中使用的术语与本发明所属领域的专业人员的使用习惯保持一致。此外，这里的描述会使用一些日常使用的词汇来描述本发明。如果读者发现这些词汇在理想情况下的定义，或者非正式的场合下的使用习惯，与上下文不一致，此时除非给出了明确的定义，否则，应当认为这里在使用这些日常词汇时，已经基于通用的各类词典里的解释对这些词汇做了符合本领域使用习惯的调整。
以下结合附图对本发明的实施例作进一步详细描述。
本发明的一种测算放疗射线剂量分布和剂量目标函数的方法和系统的实施例的特征是：包括以下步骤：
步骤1：从影像设备通过医学数字成像和通信标准(DICOM)协议接收人体或模体的三维医学影像数据；
步骤2：根据三维医学影像数据勾画和分割感兴趣区，并计算包含所有感兴趣区的边界盒(Bounding Box)。
步骤3：将放射治疗机头的每个射野的准直器围成的矩形分割成一系列网格，将通过每个网格的细小射线束记为笔形束Bn，n为笔形束编号；
步骤4：将所有笔形束均匀分配到混合集群计算节点上，并通过数据网络将笔形束的几何信息传递到该笔形束所对应的计算节点上，并为笔形束的剂量分布Dn(x，y，z)分配内存空间和显存空间，x，y，z为三维坐标点。
步骤5：利用集群节点的GPU通过光线追踪法计算每个笔形束Bn在模拟数字人体中的TERMA(单位质量总释放能量)Tn(x，y，z)，x，y，z为三维坐标点；
步骤6：将4π立体角分解为m个角度，用(n，m，i，j)元组表示的Bn的第m个方向的第(i，j)条卷积线，其中(i，j)为Bn在第m个方向上的卷积线的二维离散标识，并将该卷积线分配到集群的某一个GPU的某个流处理器线程上计算该卷积线上的第一个非零TERMA三维空间点的位置，以及从该位置沿该方向到与边界盒交点的线段所相交的三维空间离散腔胞数量Lnmij，并将它记为该卷积线的离散长度；
步骤7：计算Lnmij的的分布f(l)，根据该分布分批将具有相同离散长度的卷积线分配到集群的GPU的流处理器线程上并行计算沿该卷积线的一维卷积，
并将得到的三维剂量分布Egpu n(x，y，z)叠加到Dgpu n(x，y，z)；
步骤8：将未分配到GPU的所有剩余卷积线均匀分配到集群CPU的线程上并行计算该方向上的一维卷积，
并将得到的三维剂量分布Ecpu n(x，y，z)叠加到Dcpu n(x，y，z)；
步骤9：将每个GPU和CPU计算的笔形束剂量分布叠加起来得到最终的笔形束剂量分布：Dn(x，y，z)＝Dgpu n(x，y，z)+Dcpu n(x，y，z)。
上述的放疗射线包括光子线、电子线和质子线。
上述的集群计算节点至少包含一个或若干个CPU、一个或若干个GPU以及CPU对应的
内存和GPU对应的显存，而且集群计算节点之间通过数据网络连接；
上述的的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。
上述的的Bn在模拟数字人体中的单位质量总释放能量Tn(x，y，z)包括不同能级的TERMA的线性叠加，其中线性组合系数由放疗射线的能谱决定。
上述的的步骤7的分配算法至少包括如下子步骤：
步骤1：计算Lnmij的范围区间[lmin，lmax]，其中lmin为最小长度值，lmax为最大长度值；
步骤2：遍历区间[lmin，lmax]，对于该区间中的每一个长度l，计算f(l)/k的整数商s，其中f(l)为长度为l的卷积线的数量，k为每个GPU的流处理器数量；
步骤3：对于该区间中的每一个长度l，如果相应的s大于零，分配s×k条卷积线到s×k个GPU线程；
上述的的步骤7中的一维卷积算法的一维卷积核Km至少包括如下子步骤：
步骤1：计算蒙特卡罗方法模拟生成放疗射线点源在等效水中的三维剂量分布K3(x，y，z)作为三维卷积核；
步骤2：在第m个立体角方向对三维卷积核K3(x，y，z)进行插值采样得到Km(r)；
上述的的步骤7和步骤8可以在GPU和CPU上分别同时进行。
图1是本发明实施例的设备示意图。在实施放疗时，放疗设备生成的射束沿着射束方向108照射位于治疗床106上的病人的体内的特定区域。治疗床106可用于帮助病人定位，包括：围绕治疗床转轴103转动，在水平面内沿着径向滑动等。在放疗计划中，治疗床106和病人的位置、加速器机头的位置等被看作位于一个三轴坐标系中。该坐标系的原点为等中心点101，竖直轴为治疗床转轴103，一个水平轴为机架转轴102，另一个水平轴为104。在实施放疗时，机架从竖直的治疗床转轴103转动一个角度
到射束方向108来实施放疗。由此，机架相对于靶区的相对位置由角度
决定，射出的射束经过准直器109的限制、沿着射束方向108对靶区进行照射并在病人体内形成剂量分布。当对这一分布做计算时，可以把由准直器围城的区域110分割成网格，由每一个网格通量照射靶区而形成的细小射束为笔形束107，即射束截面110由一组形如笔形束截面111的矩形组成。这样的射束可以看作是从加速器光源112发射的。
图2是本发明实施例中的部署结构图。装有CPU 203和GPU 204的计算节点201通过数据网络202组成集群，将所有需要计算剂量分布的笔形束206分配至每一个计算节点201。将每个笔形束剂量分布的三维卷积分解成多个方向的一维卷积。预估每个一维卷积线上的计算量，也就是非零离散长度。分批次将具有相同计算量的一维卷积线207传输给GPU 204进行计算一维卷积，以达到GPU 204所需的单指令多数据(SIMD)同步，从而提升GPU 204并行效率。剩余的一维卷积线208传递给CPU 203计算卷积，由于CPU 203的多线程并行计算对执行代码没有单指令多数据(SIMD)同步的要求，因此将取得较高的效率。同时，整个方法的效率也将从CPU 203与GPU 204之间的并行执行获得受益。
图3是本发明实施例中步骤1-11的逻辑流程图。从影像设备通过DICOM协议接收三维医学影像数据之步骤301可以通过任何一种高级计算机编程语言(如C++，Java，Python等)和主计算节点上的计算机操作系统数据通讯接口(如socket或者web service)实现。根据三维医学影像数据分割感兴趣区之步骤302可以运用已有的影像编辑软件工具或者软件算法互动式地或者自动地实现。将放射治疗机头的每个射野的准直器围成的矩形分割成一系列网格，将放射治疗机头的每个射野的准直器围成的矩形分割成一系列网格，将通过每个网格的细小射线束记为笔形束Bn之步骤303通过任何一种高级计算机编程语言实现。将所有笔形束均匀分配到混合集群计算节点上，并通过数据网络将笔形束的几何信息传递到该笔形束所对应的计算节点上，并为笔形束的剂量分布Dn(x，y，z)分配内存空间和显存空间之步骤304可以通过任何一种高级计算机编程语言和并行计算库(如MPI，MassivelyParallel Interface)里的数据广播通讯接口实现。利用集群节点的GPU通过光线追踪法计算每个笔形束Bn在模拟数字人体中的TERMA(单位质量总释放能量)Tn(x，y，z)之步骤305可以通过任何一种高级计算机编程语言实现。将4π立体角分解为m个角度，用(n，m，i，j)元组表示的Bn的第m个方向的第(i，j)条卷积线，其中(i，j)为Bn在第m个方向上的卷积线的二维离散标识，并将该卷积线分配到集群的某一个GPU的某个流处理器线程上计算该卷积线上的第一个非零TERMA三维空间点的位置，以及从该位置沿该方向到与边界盒交点的线段所相交的三维空间离散腔胞数量Lnmij，并将它记为该卷积线的离散长度之步骤306可以通过任何一种高级计算机编程语言实现。计算Lnmij的的分布f(l)，根据该分布将合适的卷积线分配到集群的GPU的流处理器线程上并行计算沿该卷积线的一维卷积，并将得到的三维剂量分布En(x，y，z)叠加到Dgpu n(x，y，z)之步骤307可以通过任何一种高级计算机编程语言实现。将未分配到GPU的所有剩余卷积线均匀分配到集群CPU的线程上并行计算该方向上的一维卷积，并将得到的三维剂量分布En(x，y，z)叠加到Dgpu n(x，y，z)之步骤308可以通过任何一种高级计算机编程语言实现。将每个GPU和CPU计算的笔形束剂量分布叠加起来得到最终的笔形束剂量分布之步骤309可以通过任何一种高级计算机编程语言结合通用并行计算库(如MPI)实现。在本实施流程的步骤中涉及从计算节点上的并行计算部分，可以采用CPU上的多进程或者多线程实现，也可以采用GPU上并行计算单元实现。例如若采用NVDIA的GPU，基础控制GPU的计算函数库可以采用CUDA实现。
图4是本发明实施例中步骤7中的的分配算法子流程的逻辑流程图。该子流程中的所有子步骤均可以由高级编程语言通过系统调用实现，例如若采用NVDIA的GPU，基础控制GPU的计算函数库可以采用CUDA实现。本流程对每批笔形束的所有卷积线计算离散长度的Lnmij的范围区间[lmin，lmax]，其中lmin为最小长度值，lmax为最大长度值(步骤401)；遍历区间[lmin，lmax]，对于该区间中的每一个长度l，计算f(l)/k的整数商s，其中f(l)为长度为l的卷积线的数量，k为每个GPU的线程数(步骤402)；对于该区间中的每一个长度l，如果相应的s大于零，分配s×k条卷积线到s×k个GPU线程(步骤403)。
图5是本发明实施例中步骤7中的子步骤中的卷积线离散长度的分布统计图。Lmin为最小长度值501，Lmax为最大长度值502，f(l)为长度为l的卷积线的数量503，k为每个GPU的线程数504，f(l)/k的整数商s 505，每批分配到GPU线程的卷积线条数506为s×k条，由于该数可以被k整除，而且每批卷积线的指令内容和数量均相同(取决于卷积线的离散长度)，这样分配可以最大化GPU的SIMD并行执行效率。
本发明所涉及的已有的专业名词解释定义如下：
调强放射治疗 调强放射治疗(Intensity-Modulated Radiation Therapy，IMRT)通常采用治疗计划预先确定的连续或离散的方式，调控光子或电子束的注量、相对于患者的射束方向和射野尺寸。IMRT的主要作用是提高剂量分布对计划靶区的适形度，同时使剂量尽量规避周围正常组织。
旋转容积调强治疗 调强放射治疗(Volumetric-Modulated Arc Therapy，VMAT)是一种通过同步调整三类参数来实现高精度放疗的方法。这些参数包括：机架旋转速度、用于形成射野轮廓的MLC叶片布局、放射剂量率。VMAT可以在一次机架旋转的过程里对目标区域形成精确刻画的三维剂量分布。VMAT与IMRT的区别在于，前者可以在机架沿着弧形轨迹运动的过程中随时照射目标区域，而不必像IMRT那样局限于预先给定的少数角度方向位置。
DICOM 医学数字成像和通信标准，医学数字成像和通信(Digital Imaging andCommunications in Medicine，缩写DICOM)标准是医学图像和相关信息的国际标准。这个标准定义了质量能满足临床需要的可用于数据交换的医学图像格式。当医学影像数据按照DICOM标准进行传输和储存时，其数据文件中包含的病人信息可按照DICOM标准文件解析算法解析为一系列二维影像，可供医生阅片。
感兴趣区 感兴趣区(Region of Interest，缩写ROI)是从三维影像中选择的一个特定三维区域，在肿瘤放射治疗中通常代表肿瘤靶区或者正常组织，通常由各层影像切片上的闭合轮廓线集合来表示。这些区域是影像分析关注的重点，在放疗中一般需要计算它们各自的剂量目标函数。
靶区 靶区(Target Volume)是指放射治疗中，准备向患者体内以治疗为目的的辐照一定吸收剂量的区域。
射野 由准直器确定的射线束的边界，并垂直于射线中心轴的射线束平面。有两种定义方法：一是几何学照射野，即放射源前表面经准直器在模体表面的投影；二是物理学照射野，即以射线束中心轴剂量为100％，照射野四边50％等剂量线之间的部分。
笔形束 一束狭窄的辐射束。一组界面为矩形笔形束可以组成一个具有较大或较复杂的射野；反之，一个矩形的或经过调制具有复杂截面的射束可分解为一组笔形束。
放射线路径距离 放射线路径距离(radiological path length)是指在放疗计算中，关于指定的线段从线段的一端到另一端对线段路径上的物质密度做积分得到的量。在实际计算中，积分路径的起点可以取为辐射源、被照射物体中发生散射的点或其他指定的点，积分路径的终点可以取为任意点。
本发明的最佳实施例已阐明，由本领域普通技术人员做出的各种变化或改型都不会脱离本发明的范围。
Claims (8)
1.一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：包括以下步骤：
步骤1：从影像设备通过医学数字成像和通信标准(DICOM)协议接收人体或模体的三维医学影像数据；
步骤2：根据三维医学影像数据勾画和分割感兴趣区，并计算包含所有感兴趣区的边界盒(Bounding Box)。
步骤3：将放射治疗机头的每个射野的准直器围成的矩形分割成一系列网格，将通过每个网格的细小射线束记为笔形束Bn，n为笔形束编号；
步骤4：将所有笔形束均匀分配到混合集群计算节点上，并通过数据网络将笔形束的几何信息传递到该笔形束所对应的计算节点上，并为笔形束的剂量分布Dn(x，y，z)分配内存空间和显存空间，x，y，z为三维坐标点。
步骤5：利用集群节点的GPU通过光线追踪法计算每个笔形束Bn在模拟数字人体中的TERMA(单位质量总释放能量)Tn(x，y，z)，x，y，z为三维坐标点；
步骤6：将4π立体角分解为m个角度，用(n，m，i，j)元组表示的Bn的第m个方向的第(i，j)条卷积线，其中(i，j)为Bn在第m个方向上的卷积线的二维离散标识，并将该卷积线分配到集群的某一个GPU的某个流处理器线程上计算该卷积线上的第一个非零TERMA三维空间点的位置，以及从该位置沿该方向到与边界盒交点的线段所相交的三维空间离散腔胞数量Lnmij，并将它记为该卷积线的离散长度；
步骤7：计算Lnmii的的分布f(l)，根据该分布分批将具有相同离散长度的卷积线分配到集群的GPU的流处理器线程上并行计算沿该卷积线的一维卷积，
并将得到的三维剂量分布Egpu n(x，y，z)叠加到Dgpu n(x，y，z)；
步骤8：将未分配到GPU的所有剩余卷积线均匀分配到集群CPU的线程上并行计算该方向上的一维卷积，
并将得到的三维剂量分布Ecpu n(x，y，z)叠加到Dcpu n(x，y，z)；
步骤9：将每个GPU和CPU计算的笔形束剂量分布叠加起来得到最终的笔形束剂量分布：Dn(x，y，z)＝Dgpu n(x，y，z)+Dcpu n(x，y，z)。
2.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的放疗射线包括光子线、电子线、质子线和重离子线。
3.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的集群计算节点至少包含一个或若干个CPU、一个或若干个GPU(若一个物理节点有多个GPU，可将该物理节点虚拟成多个节点，每个虚拟节点单独使用一个GPU)，以及CPU对应的内存和GPU对应的显存，而且集群计算节点之间通过数据网络连接。
4.根据权利要求1和3所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的数据网络包括并不限于：计算机外设互联总线、局域网、广域网、互联网、以及上述任何几种网络的组合。
5.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的Bn在模拟数字人体中的单位质量总释放能量Tn(x，y，z)包括不同能级的TERMA的线性叠加，其中线性组合系数由放疗射线的能谱决定。
6.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的步骤7的分配算法至少包括如下步骤：
步骤6.1：计算Lnmij的范围区间[lmin，lmax]，其中lmin为最小长度值，lmax为最大长度值；
步骤6.2：遍历区间[lmin，lmax]，对于该区间中的每一个长度l，计算f(l)/k的整数商s，其中f(l)为长度为l的卷积线的数量，k为每个GPU的流处理器数量。
步骤6.3：对于该区间中的每一个长度l，如果相应的s大于零，分配s×k条卷积线到s×k个GPU线程；
7.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：计算所述的步骤7中的一维卷积算法的一维卷积核Km至少包括如下步骤：
步骤7.1：计算蒙特卡罗方法模拟生成放疗射线点源在等效水中的三维剂量分布K3(x，y，z)作为三维卷积核；
步骤7.2：在第m个立体角方向对三维卷积核K3(x，y，z)进行插值采样得到Km(r)。
8.根据权利要求1所述的一种基于CPU/GPU混合集群的计算放疗笔形束射线剂量分布的方法，其特征是：所述的步骤7和步骤8可以在GPU和CPU上分别同时进行。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010853607.XA CN112017166B (zh) | 2020-08-17 | 2020-08-17 | 一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010853607.XA CN112017166B (zh) | 2020-08-17 | 2020-08-17 | 一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112017166A true CN112017166A (zh) | 2020-12-01 |
| CN112017166B CN112017166B (zh) | 2024-06-14 |
Family
ID=73505590
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010853607.XA Active CN112017166B (zh) | 2020-08-17 | 2020-08-17 | 一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112017166B (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101954148A (zh) * | 2010-09-15 | 2011-01-26 | 四川大学 | 放射治疗中基于gpu的剂量计算加速方法 |
| WO2011053802A2 (en) * | 2009-10-30 | 2011-05-05 | Tomotherapy Incorporated | Non-voxel-based broad-beam (nvbb) algorithm for intensity modulated radiation therapy dose calculation and plan optimization |
| CN102201036A (zh) * | 2011-05-16 | 2011-09-28 | 四川大学 | Gpu加速剂量计算中微分卷积积分算法的实现 |
| WO2012037472A2 (en) * | 2010-09-17 | 2012-03-22 | William Marsh Rice University | Gpu-based fast dose calculator for cancer therapy |
| CN109985316A (zh) * | 2017-12-29 | 2019-07-09 | 北京连心医疗科技有限公司 | 一种复杂射野的放疗剂量快速计算方法、设备和存储介质 |
| CN110404184A (zh) * | 2019-06-13 | 2019-11-05 | 苏州同调医学科技有限公司 | 一种测算放疗射线剂量分布和剂量目标函数的方法和系统 |
| CN110504016A (zh) * | 2018-05-18 | 2019-11-26 | 北京连心医疗科技有限公司 | 一种蒙特卡罗网格并行剂量计算方法、设备和存储介质 |
| CN111494815A (zh) * | 2020-05-14 | 2020-08-07 | 安徽慧软科技有限公司 | 基于混合变尺度模型的三维剂量计算方法、装置及介质 |
-
2020
- 2020-08-17 CN CN202010853607.XA patent/CN112017166B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011053802A2 (en) * | 2009-10-30 | 2011-05-05 | Tomotherapy Incorporated | Non-voxel-based broad-beam (nvbb) algorithm for intensity modulated radiation therapy dose calculation and plan optimization |
| CN101954148A (zh) * | 2010-09-15 | 2011-01-26 | 四川大学 | 放射治疗中基于gpu的剂量计算加速方法 |
| WO2012037472A2 (en) * | 2010-09-17 | 2012-03-22 | William Marsh Rice University | Gpu-based fast dose calculator for cancer therapy |
| CN102201036A (zh) * | 2011-05-16 | 2011-09-28 | 四川大学 | Gpu加速剂量计算中微分卷积积分算法的实现 |
| CN109985316A (zh) * | 2017-12-29 | 2019-07-09 | 北京连心医疗科技有限公司 | 一种复杂射野的放疗剂量快速计算方法、设备和存储介质 |
| CN110504016A (zh) * | 2018-05-18 | 2019-11-26 | 北京连心医疗科技有限公司 | 一种蒙特卡罗网格并行剂量计算方法、设备和存储介质 |
| CN110404184A (zh) * | 2019-06-13 | 2019-11-05 | 苏州同调医学科技有限公司 | 一种测算放疗射线剂量分布和剂量目标函数的方法和系统 |
| CN111494815A (zh) * | 2020-05-14 | 2020-08-07 | 安徽慧软科技有限公司 | 基于混合变尺度模型的三维剂量计算方法、装置及介质 |
Non-Patent Citations (2)
| Title |
|---|
| RINTARO FUJIMOTO等: "GPU-based fast pencil beam algorithm for proton therapy", 《PHYSICS IN MEDICINE & BIOLOGY》, vol. 56, no. 5, pages 1319, XP020191144, DOI: 10.1088/0031-9155/56/5/006 * |
| 王先良等: "GPU加速剂量计算中微分卷积/积分算法的实现", 《核技术》, vol. 36, no. 05, pages 58 - 63 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112017166B (zh) | 2024-06-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109069858B (zh) | 一种放射治疗系统及计算机可读存储装置 | |
| US11020615B2 (en) | Computing radiotherapy dose distribution | |
| US10507337B2 (en) | Radiotherapy treatment plan optimization workflow | |
| US20210339049A1 (en) | Radiotherapy treatment plans using differentiable dose functions | |
| EP3742702B1 (en) | Standardized cloud radiotherapy planning method and storage medium | |
| CN101120871A (zh) | 精确放射治疗计划系统 | |
| US9495513B2 (en) | GPU-based fast dose calculator for cancer therapy | |
| CN110404184A (zh) | 一种测算放疗射线剂量分布和剂量目标函数的方法和系统 | |
| Pastor-Serrano et al. | Millisecond speed deep learning based proton dose calculation with Monte Carlo accuracy | |
| CN114846476A (zh) | 训练用于放射疗法治疗计划的深度学习引擎 | |
| EP4230261A1 (en) | Dose volume histogram optimization based on quantile regression | |
| Pastor-Serrano et al. | Learning the physics of particle transport via transformers | |
| Jia et al. | Proton therapy dose calculations on GPU: advances and challenges | |
| CN112017166B (zh) | 一种基于cpu/gpu混合集群的计算放疗笔形束射线剂量分布的方法 | |
| EP4101502A1 (en) | Feature-space clustering for physiological cycle classification | |
| Endo | Creation, evolution, and future challenges of ion beam therapy from a medical physicist’s viewpoint (Part 2). Chapter 2. Biophysical model, treatment planning system and image guided radiotherapy | |
| Giordanengo et al. | RIDOS: A new system for online computation of the delivered dose distributions in scanning ion beam therapy | |
| Moreno et al. | Parallel gEUD Models for Accelerated IMRT Planning on Modern HPC Platforms | |
| CN110404186A (zh) | 一种设计非共面容积旋转调强放疗计划的方法 | |
| Matter | Technical implementation of daily adaptive proton therapy | |
| Neph | Accelerating Radiation Dose Calculation with High Performance Computing and Machine Learning for Large-scale Radiotherapy Treatment Planning | |
| Tseng | Dose Calculation Strategies Toward Efficient and Accurate Treatment Planning for Modern Radiotherapy | |
| Folkerts | Topics in Cancer Radiotherapy: Automated Treatment Planning and Quality Assurance | |
| Scholz | Development and evaluation of advanced dose calculations for modern radiation therapy techniques | |
| Quetin et al. | Deep learning for high-resolution dose prediction in high dose rate brachytherapy for breast cancer treatment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |