基于GPU蒙特卡洛算法的磁场下光子和电子剂量计算方法与流程

本发明属于计算机信息技术在核技术领域的应用，尤其涉及到医学放射治疗技术领域，具体的说是一种基于gpu蒙特卡洛算法的磁场下光子和电子剂量计算方法，可用于磁共振实时引导放射治疗mrigrt的治疗计划系统中的剂量计算。

背景技术：

癌症是人类死亡的主要原因之一，大约有70％的癌症患者在治疗癌症的过程中需要使用放射治疗，约有40％的癌症可以用放射治疗根治。实际患者接受的放射治疗方案是通过治疗计划系统来设计的。在治疗计划系统中，医生会根据病人的ct勾画出靶区器官与危及器官的轮廓，并设置好靶区器官的处方剂量及危及器官的剂量限值，然后这些信息会传送到物理师工作站，物理师在上述基础上，为病人制定的具体的放疗方案。这个方案主要是指设置病人实际放疗时放射源的能量、照射方向、照射野数目和权重。通过治疗计划系统计算出物理师制定的放疗方案是否满足医生规定的靶区剂量和危及器官剂量要求。若满足则按照物理师方案中提到的放射源对患者进行照射。反之则重新制定治疗计划，修改放射源数目、照射方向和权重，直到剂量计算的靶区剂量和危及器官剂量达到医生设定的要求。剂量计算主要计算的是放射源与人体模型相互作用时，在人体模型中沉积的剂量。人体模型主要是通过患者ct数据获取而来，放射源主要是参照上述物理师为病人设置的放射源的能量、照射方向、照射野数目和权重。

目前的治疗计划系统中的剂量计算主要是采用解析法来计算靶区剂量和危及器官剂量，解析法适合处理均匀介质的剂量计算问题，但是患者的器官是非均匀介质，解析法在处理非均匀介质的剂量计算问题时存在较大误差。另一方面，治疗计划系统主要支持ct的数据，制定放射治疗计划方案。导致病人必须通过ct扫描才能获取人体模型数据，对于一些运动器官，为了获得它们的清晰图像，体内需植入标记物，这让不少的病人很痛苦。临床数据表明ct获得的软组织分辨率不高，运动器官即使植入标记物，图像仍然不清晰，且在ct扫描时患者受到了非治疗需要的辐射剂量。

随着放疗技术与成像技术的发展，国际上提出用mri引导放射治疗，即在治疗前、治疗中通过mri来引导整个放射治疗的过程。mri可实时跟踪器官运动变化，无辐射剂量，并能获得丰富的解剖结构、功能性信息。但mri中存在磁场，而目前的治疗计划系统剂量计算并未将磁场的因素考虑，这将使得目前的治疗计划系统无法准确计算出磁场存在时的辐射剂量。另外通过pet也可以获取患者的图像信息，但目前的治疗计划系统暂时不支持pet成像的数据，用于剂量计算。

蒙特卡洛方法一直是公认计算辐射剂量的“黄金准则”。蒙特卡洛方法能解决在治疗计划系统中的剂量计算遇到的非均匀介质问题。但蒙特卡洛方法是基于统计的思想，需进行大量模拟才能保证结果的准确性，所以用蒙特卡洛方法解决治疗计划系统中的剂量计算问题，将花费很长的时间才能获得剂量结果。正是因为蒙特卡洛方法花费时间长，即使它能准确计算出剂量结果，该方法也一直无法应用于临床。

随着计算机的发展，基于kepler构架和maxwell构架的gpu在当前的科学计算中使用十分广泛。将已有的cpu大型模拟程序如geant4、mcnp移植到gpu上，需手动改写代码，且gpu上十几gb内存分配到几千个线程，对于基于蒙特卡洛方法的程序而言，这是远远不足的。即使移植成功，也会因为gpu上的线程分歧、内存延时，c++多态功能导致编译器无法进行函数内联导致运行速度很慢。因而，目前若将基于蒙特卡洛方法的放射源粒子输运程序与gpu相结合，只能用于特定用途的蒙特卡洛程序模拟。

另一方面，基于kepler构架和maxwell构架的gpu在科学计算时存在缺点。它们缺乏固有双精度原子加法运算的硬件支持。传统的解决方法是使用nvidia的“比较-交换”算法从软件上效仿双精度原子加法运算。这个算法的特点是步骤不定长，gpu的线程竞争越激烈，所需步骤越多。当前gpu硬件上采用“单指令-多线程”技术，一个线程包内的32个线程在时间上以“锁步”(lockstep)方式执行指令。这样导致的结果是：如果32个线程有某些需要更新同一处全局内存，则会产生激烈的线程竞争，“比较-交换”算法所需步骤急剧增加，运算时间大幅增长。因此，将gpu与蒙特卡洛程序相结合，若不解决该问题，即使是用于特定用途的蒙特卡洛程序模拟，也无法获得快速且准确的剂量结果。

技术实现要素：

本发明是为了解决上述现有技术存在的不足之处，提出一种基于gpu蒙特卡洛算法的磁场下光子和电子剂量计算方法，以期能快速且准确的计算出磁场作用下光子和电子的辐射剂量，可用于磁共振实时引导放射治疗mrigrt的治疗计划系统中的剂量计算，进而提高mrigrt治疗计划系统中剂量计算的准确性与速度，改善放射治疗的效果。

本发明为解决技术问题采用如下技术方案：

本发明一种基于gpu蒙特卡洛算法的磁场下光子和电子剂量计算方法的特点是按如下步骤进行：

步骤1：采集数据

步骤1.1、获取放射治疗加速器的照射源数据并进行处理，得到放射源信息其中，e表示源能量，表示源位置，表示源发射方向；

步骤1.2、获取人体解剖结构的图像数据并重建人体模型；获取核磁共振仪的磁场强度数据

步骤1.3、获取光子和电子分别与物质发生反应的核数据并进行处理，得到所述核数据的宏观截面数据σ，对所述宏观截面数据按照能量的高低进行降序排序，得到排序后的光子宏观截面数据σ^p和电子宏观截面数据σ^e；

步骤2、确定gpu的最优线程数和输运任务的批次；

步骤2.1、利用runtimeattribute程序接口获得gpu中每个线程所需寄存器的数目r；则gpu中每个流多处理器工作在满载状态的最小线程个数为r表示每个流多处理器上的寄存器个数；从而得到gpu工作在满载状态所需线程总数为t＝mt，m表示gpu中流多处理器的个数；

步骤2.2、设置放射源粒子的数目为n，并将n个放射源粒子的输运任务划分为t个批次，使得每个批次上以串行地方式待计算放射源粒子的个数为

步骤3、利用蒙特卡洛算法计算每个批次在磁场作用下的光子和电子辐射剂量；

步骤3.1、定义每一批抽取放射源粒子的次数为w，并初始化w＝1；

步骤3.2、定义第w次抽取放射源粒子时的输运次数为u，并初始化u＝0；

步骤3.3、利用随机数生成器从放射源信息s中第w次抽取第w个放射源粒子sw；所述第w个放射源粒子sw第u次输运的状态为：并判断第w次抽取第w个放射源粒子的类型，若为光子记为则执行步骤3.4；若为电子记为则执行步骤3.5；

步骤3.4、基于排序后的光子宏观截面数据σ^p抽取第u+1次输运时的运动步长和运动方向再执行步骤3.7；

步骤3.5、基于排序后的电子宏观截面数据σ^e抽取第u+1次输运时的运动步长和运动方向

步骤3.6、判断所述电子是否处于人体模型的磁场区，若是，先将第w个电子沿运动方向移动距离再将第w个电子沿着式(1)修正的运动方向移动否则，仍然采用所述运动方向

式(1)中，为第w个电子第u+1次输运时修正后的新方向；norm{}为归一化算符；为第u+1次对第w个电子进行抽样所得到的步长；q为电子电荷数；c为真空中光速，m为电子静质量，为第w个电子在第u次输运时的能量；为电子第u+1次输运时未修正的运动方向；为第w个电子在人体模型中第u+1次输运时所处位置的磁场强度；

步骤3.7、基于排序后的光子宏观截面数据σ^p和电子宏观截面数据σ^e、所述第w个放射源粒子sw第u次输运的状态运动步长、运动方向或修正的运动方向，对第w个放射源粒子sw的反应类型进行抽样，得到第w个放射源粒子sw在人体模型中进行第u+1次输运的状态

步骤3.8、计算第w次抽取的第w个放射源粒子sw人体模型中进行第u+1次输运时沉积的剂量

步骤3.9、将u+1赋值给u，并返回步骤3.4执行，直到能量低于截止能量或第w个放射源粒子sw超出人体模型的边界为止；从而统计得到第w次抽取的第w个放射源粒子sw在人体模型中沉积的剂量dosew；

步骤3.10、将w+1赋值给w，并返回步骤3.1执行，直到w＞n为止；从而统计得到每个批次n次抽取的n个放射源粒子在人体模型中沉积的剂量dose；

步骤4、基于gpu快速原子加法统计剂量结果：

步骤4.1、将t个批次所获得沉积的剂量以三维矩阵的形式存入gpu的全局内存中，在所述三维矩阵中的任意元素记为dosei；dosei表示在三维空间中第i个位置上的剂量；

步骤4.2、获取gpu中同一个线程包内的有效线程，并将所要更新的全局内存地址相同的有效线程从所述同一线程包中筛选出来；

步骤4.3、将筛选出来的有效线程所对应的全局内存中所存储的剂量进行累加，得到的剂量结果存入序号最小的有效线程中；

步骤4.4、利用“卡汉求和”以及“比较-交换”算法将所有属于同一线程包内序号最小的有效线程中的剂量结果累加到gpu的全局内存中，从而得到人体模型中放射源粒子的总剂量；

步骤4.5、将所述总剂量除以放射源粒子的总数n，从而得到归一化的剂量结果，并将其返回到cpu内存。

与现有技术相比，本发明的有益效果在于：

1、本发明利用图像数据构建人体模型；并调用cudaruntimeapi里面的函数cudafuncattributes()获得gpu中每个线程所需寄存器的数目r，进而计算得到gpu的最优线程分配方案；将n个放射源粒子的模拟任务在gpu上划分为t个批次，每个批次同步执行；执行时基于排序后的宏观截面数据σ，利用蒙特卡洛算法计算每个批次在磁场作用下的光子和电子辐射剂量，对磁场作用下的带电粒子，需经修正运动方向后再输运；最后采用一种新的gpu快速原子加法统计剂量结果。从而解决了4个主要问题，包括：现有治疗计划系统中采用解析法无法准确计算放射源在遇到人体模型中的非均匀介质时沉积的剂量问题；现有治疗计划系统中未考虑磁场对带电粒子的影响，进一步影响沉积的剂量问题；现有治疗计划系统中若采用蒙特卡洛方法，计算剂量速度慢的问题；基于kepler构架和maxwell构架的gpu在统计剂量结果时缺乏固有双精度原子加法运算的硬件支持问题。因而本发明的算法可优化目前治疗计划系统的剂量计算，尤其在磁场作用下基于gpu的蒙特卡洛粒子输运过程可改进目前治疗计划系统中的解析法处理粒子输运过程。根据该方法可快速且准确地计算出无磁场和有磁场情况下的辐射剂量，极大提高了临床工作效率，改善了放射治疗的效果，从而减少了癌症患者的痛苦，提升了癌症患者接受放射治疗后的生活质量。

2、本发明中利用gpu的最优线程分配法完成n个放射源粒子的模拟任务，该方法利用cudaruntimeapi里面的函数cudafuncattributes()获得gpu中每个线程所需寄存器的数目r，通过r进一步获得gpu中每个流多处理器工作在满载状态的最小线程个数t，最后计算出gpu工作在满载状态时所需线程总数t，解决了gpu上线程不能充分合理利用的问题，从而能够最大化gpu的并行能力，提高了模拟的速度。

3、本发明在gpu上输运粒子之前，对宏观截面数据σ按照能量排序，从而提高了实际模拟粒子输运任务时查找宏观截面数据的速度，进而提高了模拟的速度。

4、本发明在gpu上模拟粒子的输运任务时，利用蒙特卡洛算法计算每个批次在磁场作用下的光子和电子辐射剂量，尤其针对磁场会导致剂量的变化，在输运过程中做了带电粒子运动方向的修正，优化了原始的蒙特卡洛方法，使得新算法适用范围更广；从而解决了实际中用蒙特卡洛方法无法快速处理磁场存在时的剂量计算问题，将本发明应用于核磁引导的治疗计划系统中的剂量计算，提高了剂量结果的准确性与计算速度。

5、本方法中提出一种新的gpu快速原子加法统计剂量结果，采用卡汉求和算法、nvidia的韦斯福算法和nvidia的“比较-交换”算法从软件上效仿双精度原子加法运算，解决了基于kepler构架和maxwell构架的gpu缺乏固有双精度原子加法运算的硬件支持问题，并且可快速准确实现任意类型数据的计算，还减少了线程包内部的线程竞争，提升了单精度原子加法运算的精度和双精度原子加法运算的速度，因而提高了统计剂量结果的准确性与速度。

附图说明

图1为本发明的整体流程图。

具体实施方式

本实施例中，一种基于gpu蒙特卡洛算法的磁场下光子和电子剂量计算方法，是应用于mrigrt领域，即在mri磁场环境下，放射源通过与人体模型之间复杂的相互作用，将能量沉积在人体模型的各个器官中的实验环境中，具体的说，如图1所示，按如下步骤进行：

步骤1：采集数据

步骤1.1、获取放射治疗加速器的照射源数据并进行处理，得到放射源信息其中，e表示源能量，表示源位置，表示源发射方向；

步骤1.2、获取ct机或核磁共振仪mri或正电子发射型计算机断层显像机器pet的人体解剖结构的图像数据并重建人体模型；获取核磁共振仪的磁场强度数据本实施例中采集的图像数据支持更多种格式，解决了现有技术中无法通过ct机获取人体模型中软组织、运动器官准确位置的问题，并解决了经ct扫描获得人体模型数据时会受到放射源带来的辐射剂量问题，支持mri和pet格式的数据提高了在人体模型中所包含的肿瘤区域的界限精准度，减少了癌症患者为获得人体模型数据需要进行体内植入标记物的痛苦；

步骤1.3、获取光子和电子分别与物质发生反应的核数据并进行处理，得到核数据的宏观截面数据σ，对宏观截面数据按照能量的高低进行降序排序，得到排序后的光子宏观截面数据σ^p和电子宏观截面数据σ^e；

反应的核数据来自于开源的核数据库。核数据库中的微观截面数据σ^p与σ^e经处理才能获得光子截面数据σ^p和电子截面数据σ^e。因为光子截面数据与电子截面数据计算方法相同，故将微观截面用σ表示，宏观截面数据用σ表示，σ具体计算是按如下步骤进行的：

计算粒子发生反应的宏观截面，其中n是单位体积内的原子核数：

1.3.1对于单元素材料：σ＝nσ，n0是阿伏伽德罗常数，n0＝6.0221367×10²³mol^-1，ρ为材料的密度；a为该元素的原子量，根据ct数据获得材料密度，a。

1.3.2对于混合物：σ＝σiniσi，wi为i元素在混合物中的重量百分比。

1.3.3对于化合物：σ＝σiniσi，αi为每个分子中含i种元素的原子数目。

宏观截面数据计算完成后，按照能量顺序存储，便于后面反应类型的抽取。

步骤2、确定gpu的最优线程数和输运任务的批次；

步骤2.1、利用runtimeattribute程序接口获得gpu中每个线程所需寄存器的数目r；则gpu中每个流多处理器工作在满载状态的最小线程个数为r表示每个流多处理器上的寄存器个数；从而得到gpu工作在满载状态所需线程总数为t＝mt，m表示gpu中流多处理器的个数；

步骤2.2、设置放射源粒子的数目为n，并将n个放射源粒子的输运任务划分为t个批次，则在每个批次上以串行地方式待计算放射源粒子的个数为

步骤3、利用蒙特卡洛算法计算每个批次在磁场作用下的光子和电子辐射剂量；

步骤3.1、定义每一批抽取放射源粒子的次数为w，并初始化w＝1；

步骤3.2、定义第w次抽取放射源粒子时的输运次数为u，并初始化u＝0；

步骤3.3、利用随机数生成器从放射源信息s中第w次抽取第w个放射源粒子sw；第w个放射源粒子sw第u次输运的状态为：并判断第w次抽取第w个放射源粒子的类型，若为光子记为则执行步骤3.4；若为电子记为则执行步骤3.5；

步骤3.4、基于排序后的光子宏观截面数据σ^p抽取第u+1次输运时的运动步长和运动方向再执行步骤3.7；

步骤3.5、基于排序后的电子宏观截面数据σ^e抽取第u+1次输运时的运动步长和运动方向

步骤3.6、判断电子是否处于人体模型的磁场区。若是，先将第w个电子沿运动方向移动距离再将第w个电子沿着式(1)修正的运动方向移动否则，仍然采用运动方向

式(1)中，为电子第u+1次输运时修正后的新方向；norm{}为归一化算符；为第u+1次对第w个电子进行抽样所得到的步长；q为电子电荷数；c为真空中光速，m为电子静质量，为第w个电子在第u次输运时的能量；为电子第u+1次输运时未修正的运动方向；为第w个电子在人体模型中第u+1次输运时所处位置的磁场强度；

根据原先的蒙特卡洛方法并不能计算磁场下的电子剂量问题，本方法考虑了磁场对电子的影响，修正了电子的运动方向，进而获得较为准确的磁场下电子的运动轨迹，有助于优化原先的蒙特卡洛方法，拓宽它的应用范围。本算法可用于核磁引导的治疗计划系统的剂量计算。

步骤3.7、基于排序后的光子宏观截面数据σ^p和电子宏观截面数据σ^e、第w个放射源粒子sw第u次输运的状态运动步长、运动方向或修正的运动方向，对第w个放射源粒子sw的反应类型进行抽样，得到第w个放射源粒子sw在人体模型中进行第u+1次输运的状态

在gpu中利用蒙特卡洛方法计算每个批次在磁场作用下的光子和电子辐射剂量，其中的步骤3.7进行反应类型的抽样及后续详细的步骤如下：

利用随机数生成器生成一个随机数，随机数的范围属于四大反应类型中某一反应类型对应的宏观截面范围，则第w次抽取第w个放射源粒子发生这一类型的反应。这四大类型反应分别是光电效应、康普顿散射、瑞利散射和电子对效应。每种类型产生的次级粒子包括光子和电子或正电子，将次级粒子以“栈”的形式存储。

对于光电效应，首先跳转到步骤4用快速原子加法算法将光子能量e累加到所有线程公用的剂量矩阵，结束当前粒子输运。然后判断次级粒子栈中是否有剩余粒子，有则回到步骤3.3继续模拟下一个次级粒子。反之，回到步骤3.1继续模拟下一个初级粒子，直到当前线程所有n个粒子全部模拟完毕。

对于康普顿散射，利用随机数发生器随机抽取散射光子的能量enew。用快速原子加法算法将损失的能量e-enew累加到所有线程公用的剂量矩阵，康普顿散射光子的状态则跳转到步骤3.4；

对于瑞利散射，瑞利散射光子的状态则跳转到步骤3.4；

对于电子对效应，生成一对电子和正电子，利用随机数发生器采样他们的能量(enew,1.022mev-enew)，将电子-正电子放到次级粒子栈中，如果次级粒子栈中有剩余粒子，回到步骤3.3继续模拟下一个次级粒子。否则，回到步骤3.1继续模拟下一个初级粒子，直到当前线程所有n个粒子全部模拟完毕。

步骤3.8、计算第w次抽取的第w个放射源粒子sw人体模型中进行第u+1次输运时沉积的剂量

步骤3.9、将u+1赋值给u，并返回步骤3.4执行，直到低于截止能量或第w个放射源粒子sw超出人体模型的边界为止；从而统计得到第w次抽取的第w个放射源粒子sw在人体模型中沉积的剂量dosew；

步骤3.10、将w+1赋值给w，并返回步骤3.1执行，直到w＞n为止；从而统计得到每个批次n次抽取的n个放射源粒子在人体模型中沉积的剂量dose；

步骤4、基于gpu快速原子加法统计剂量结果，累加到所有线程公用的剂量矩阵。

步骤4.1、将t个批次所获得沉积的剂量以三维矩阵的形式存入gpu的全局内存中，在三维矩阵中的任意元素记为dosei；dosei表示在三维空间中第i个位置上的剂量；

步骤4.2、获取gpu中同一个线程包内的有效线程，并将所要更新的全局内存地址相同的有效线程从同一线程包中筛选出来；对于第k个有效线程，筛选出同一个线程包内与第k个线程所要更新的全局内存地址相同的其他所有有效线程。

有效线程(activethread)指的是对于程序的分支结构正处于当前分支内的线程。比如对于一块if-else代码，当程序正在执行else时，线程包里的32个线程有些可能本应执行if分支，这些线程会被硬件屏蔽掉，成为无效线程，剩下就是有效线程。

对于第k个有效线程，用固有函数__ballot()筛选出其他所有有效线程，得到一个32比特的标记mk，第i个比特表示第i个线程是否是有效线程，1为有效，0为无效。

将第j个有效线程所要更新的全局内存地址aj用固有函数__shuffle()传到第k个有效线程，再用__ballot()筛选出线程包内满足ak＝aj的所有有效线程，将结果存储在标记mk中：如果ak＝aj，则令mk中第j个比特为1，否则令其为0。如果所有有效线程均不满足ak＝aj，则对于第j+1个有效线程重复本步操作，直到所有有效线程都处理完毕。

步骤4.3、将筛选出来的有效线程所对应的全局内存中所存储的剂量进行累加，得到的剂量结果存入序号最小的有效线程中；

步骤4.4、利用“卡汉求和”以及“比较-交换”算法将所有属于同一线程包内序号最小的有效线程中的剂量结果累加到gpu的全局内存中，从而得到人体模型中放射源粒子的总剂量；

步骤4.5、将总剂量除以放射源粒子的总数n，从而得到归一化的剂量结果，并将其返回到cpu内存。

以上步骤，从软件的角度实现了基于kepler构架和maxwell构架的gpu支持双精度原子加法运算，解决了基于kepler构架和maxwell构架的gpu硬件上不支持双精度原子加法运算的问题。具体的是解决线程包内部的线程竞争，实现高精度与低精度的混合计算时不丢失高精度的数据信息，进而提升单精度原子加法运算的精度和双精度原子加法运算的速度。

本算法的整个流程就是利用gpu硬件平台，采用改进的蒙特卡洛方法快速准确地模拟计算出病人在磁场环境下做图像引导放射治疗时沉积的剂量。