[发明专利]基于GPU蒙特卡洛算法的磁场下质子和重离子剂量计算方法

摘要

本发明公开了一种基于GPU蒙特卡洛算法的磁场下质子和重离子剂量计算方法，包括：1、采集数据；2、确定GPU的最优线程数和输运任务批次；3、利用蒙特卡洛算法计算每个批次在磁场作用下的质子和重离子辐射剂量；4、利用蒙特卡洛算法计算每个批次在磁场作用下的次级粒子的辐射剂量；5、基于GPU快速原子加法统计剂量结果。本发明能快速且准确的计算出磁场作用下质子和重离子的辐射剂量，可用于磁共振实时引导质子和重离子治疗计划系统的剂量计算，进而提高治疗计划系统中剂量计算的准确性与速度，改善放射治疗的效果。

主权项

1.一种基于GPU蒙特卡洛算法的磁场下质子和重离子剂量计算方法，其特征是按如下步骤进行：步骤1：采集数据；步骤1.1、获取放射治疗质子或重离子加速器的束流源信息其中，E表示源能量，表示源位置，表示源发射方向，σ是高斯分布的标准差；步骤1.2、获取人体解剖结构的图像数据并重建人体模型；获取核磁共振仪的磁场强度数据步骤1.3、获取质子和重离子分别与物质发生反应的核数据并进行处理，得到所述核数据的宏观截面数据∑，对所述宏观截面数据∑按照能量的高低进行降序排序，得到排序后的质子宏观截面数据∑^p、重离子宏观截面数据∑^Ion、次级粒子宏观截面数据∑^psec＝{∑^e,∑^D,∑^α}；∑^e为次级电子宏观截面数据；∑^D为氘核宏观截面数据；∑^α为α粒子宏观截面数据；步骤2、确定GPU的最优线程数和输运任务的批次；步骤2.1、利用runtime attribute程序接口获得GPU中每个线程所需寄存器的数目r；则GPU中每个流多处理器工作在满载状态的最小线程个数为R表示每个流多处理器上的寄存器个数；从而得到GPU工作在满载状态所需线程总数为T＝Mt，M表示GPU中流多处理器的个数；步骤2.2、设置放射源粒子的数目为N，并将N个放射源粒子的输运任务划分为T个批次，使得每个批次上以串行地方式待计算放射源粒子的个数为步骤3、利用蒙特卡洛算法计算每个批次在磁场作用下的初级质子和重离子辐射剂量；步骤3.1、定义每一批抽取初级放射源粒子的次数为w，并初始化w＝1；定义次级粒子的次数为w′，并初始化w′＝1；步骤3.2、定义第w次抽取放射源粒子时的输运次数为u，并初始化u＝0；定义第w′次抽取次级粒子时的输运次数为u′，并初始化u′＝0；步骤3.3、利用随机数生成器从放射源信息S中第w次抽取第w个放射源粒子S_w；所述第w个放射源粒子S_w第u次输运的状态为：步骤3.4、判断第w次抽取第w个放射源粒子S_w的类型，若为质子记为则基于排序后的质子宏观截面数据Σ^p抽取第u+1次输运时对应粒子类型的运动步长和运动方向若为重离子记为则基于排序后的重离子宏观截面数据Σ^Ion抽取第u+1次输运时对应粒子类型的运动步长和运动方向步骤3.5、判断所述放射源粒子S_w是否处于人体模型的磁场区，若是，先将第w个放射源粒子S_w沿运动方向移动的距离，再将第w个放射源粒子S_w沿着式(1)修正的运动方向移动的距离；否则，仍然采用所述运动方向式(1)中，为质子或重离子第u+1次输运时修正后的新方向；norm{}为归一化算符；d_u+1为第u+1次对第w个质子或重离子进行抽样所得到的步长；Q为质子或重离子的电荷数；c为真空中光速，m为质子或重离子质量，为第w个质子或重离子在第u次输运时的能量；为质子或重离子第u+1次输运时未修正的运动方向；为第w个质子或重离子在人体模型中第u+1次输运时所处位置的磁场强度；步骤3.6、基于排序后的质子宏观截面数据Σ^p和重离子宏观截面数据Σ^Ion、所述第w个放射源粒子S_w第u次输运的状态运动步长、运动方向或修正的运动方向，对第w个放射源粒子S_w的反应类型进行抽样，得到第w个放射源粒子S_w在人体模型中进行第u+1次输运的状态以及第u次输运过程产生的个次级粒子并存储；步骤3.7、计算第w次抽取的第w个放射源粒子S_w在人体模型中进行第u+1次输运时沉积的剂量步骤3.8、将u+1赋值给u，并返回步骤3.4执行，直到第w个放射源粒子S_w第u次输运的能量低于截止能量或第w个放射源粒子S_w超出人体模型的边界为止，记第w个放射源粒子S_w的输运总次数为U，从而统计得到第w次抽取的第w个放射源粒子S_w在人体模型中沉积的剂量Dose_w；步骤3.9、将w+1赋值给w，并返回步骤3.1执行，直到w＞n为止；从而统计得到每个批次n次抽取的n个放射源粒子在人体模型中沉积的剂量Dose；步骤3.10、统计得到所有T个批次的n个放射源粒子输运U次产生的次级粒子的数目为并将N′个次级粒子的输运任务划分为T个批次，使得每个批次上以串行地方式待计算次级粒子的个数为步骤4、利用蒙特卡洛算法计算每个批次在磁场作用下的次级粒子的辐射剂量；步骤4.1、判断N′个次级粒子是否是在质子输运过程产生的，若是，则表示N′个次级粒子是由次级质子、次级电子、氘核和α粒子组成，并执行步骤4.2；否则，表示N′个次级粒子是由次级¹²C粒子、α粒子、次级质子和次级¹¹C粒子组成，并执行步骤4.3；步骤4.2、对次级质子、次级电子、氘核和α粒子分别按如下方式进行输运：若为次级质子，则按照步骤3.4‑步骤3.9进行输运；若为次级电子，则基于次级电子宏观截面数据∑^e，按照步骤3.4‑步骤3.9进行输运；若为氘核，则基于氘核宏观截面数据∑^D，按照步骤3.4‑步骤3.9进行输运；若为α粒子，则基于α粒子宏观截面数据∑^α，按照步骤3.4‑步骤3.9进行输运；步骤4.3、对次级¹²C粒子、α粒子、次级质子和次级¹¹C粒子分别按如下方式进行输运：若为次级¹²C粒子，则按照步骤3.4‑步骤3.9进行输运；若为α粒子，则基于α粒子宏观截面数据∑^α，按照步骤3.4‑步骤3.9进行输运；若为次级质子，则按照步骤3.4‑步骤3.9进行输运；若为次级¹¹C粒子，则按照步骤3.4‑步骤3.9进行输运；步骤5、基于GPU快速原子加法统计剂量结果：步骤5.1、将T个批次所获得的放射源粒子和次级粒子的沉积剂量以三维矩阵的形式存入GPU的全局内存中，在所述三维矩阵中的任意元素记为dose_i；dose_i表示在三维空间中第i个位置上的剂量；步骤5.2、获取GPU中同一个线程包内的有效线程，并将所要更新的全局内存地址相同的有效线程从所述同一个线程包中筛选出来；步骤5.3、将筛选出来的有效线程所对应的全局内存中所存储的剂量进行累加，得到的剂量结果存入序号最小的有效线程中；步骤5.4、利用“卡汉求和”以及“比较‑交换”算法将所有属于同一线程包内序号最小的有效线程中的剂量结果累加到GPU的全局内存中，从而得到人体模型中放射源粒子和次级粒子的总剂量；步骤5.5、将所述总剂量除以放射源粒子的总数N，从而得到归一化的剂量结果，并将其返回到CPU内存。