一种基于任务驱动的蒙特卡洛散射光子模拟方法与流程

本发明涉及散射光子模拟技术领域，特别是涉及一种基于任务驱动的蒙特卡洛散射光子模拟方法。

背景技术：

锥形束计算机断层扫描(cone-beamcomputedtomography,cbct)系统由于其体积小，重量轻，扫描速度快在医院众多临床科室里有广泛应用，特别是在牙科ct，乳腺ct和图像引导放疗中。由于cbct系统采用的探测器为平板探测器，无法采用准直器技术遮挡散射线，所以重建的cbct图像中存在严重的散射伪影。

为使cbct图像在临床中深入应用，国内外研究人员对散射伪影的校正开展了大量的研究。boonejm将其归纳为基于硬件和软件的两种散射校正技术类型。其中基于硬件的散射校正方法会带来其他各种相关难题和缺陷，最后得不偿失。另一类方法是基于软件的散射校正方法，通过计算实现散射分布的估计进而实现散射校正。其中以设计散射核函数进行反卷积的方法为代表。该方法计算速度快，具有良好的临床应用价值，但散射核设计难度导致散射估计不够精确，目前只能利用骨组织形态和轮廓进行摆位确认和校正(sunm,star-lackj:improvedscattercorrectionusingadaptivescatterkernelsuperposition.physicsinmedicineandbiology2010,55(22):6695-6720.)。

此外还有另一类更为精确的方法是通过蒙特卡罗(montecarlo,mc)模拟散射光子物理效应获取散射分量实现散射校正。基于mc模拟粒子输运的计算方法是肿瘤放射治疗领域有关粒子物理输运问题的金标准。常规mc模拟方法如brownfb开发的mcnp(brownfb:mcnp–ageneralmontecarlon-particletransportcode,version5.losalamosnationallaboratory,oakridge,tn2003.)和xunjia开发的gmcdrr(jiax,yanh,l,folkertsm,jiangsb:agputoolforefficient,accurate,andrealisticsimulationofconebeamctprojections.medicalphysics2012)，相比于penelope(schmidtkej,engelw.penelope:analgorithmformontecarlosimulationofthepenetrationandenergylossofelectronsandpositronsinmatter[j].nuclearinstruments&methodsinphysicsresearch,1995,100(1):31-46.)这些经典mc模拟，在模拟效率和精确性上有了提升。但无论是penelope，mcnp还是gmcdrr，它们的模拟都是被动式的而非主动式的。被动式模拟是指在光子mc模拟中，依次对光子的射出方向，第一个散射点发生位置、类型和角度，第二个散射点发生的位置、类型和散射角度等独立抽样，直到光子离开扫描物体，最后判断光子是否到达探测器以决定是否沉积该路径光子的能量。被动式光子散射独立抽样策略思路比较自然且代码结构简单，但缺点是只能在模拟所有散射阶数和成分的光子输运后，被动式地从结果中挑选出沉积到探测器区域的散射分量。而大量未能到达探测器的光子则被丢弃，导致花费了大量时间模拟了对最终散射信号无贡献的光子。所以常规mc模拟方法由于无法对光子整体路径的概率进行计算，也就忽视光子整体路径的重要性分布，导致模拟效率过低。

为解决mc模拟的收敛速度慢和效率过低的瓶颈问题，常见的做法是引入各种减方差的技巧(haghighata,wagnerjc:montecarlovariancereductionwithdeterministicimportancefunctions.progressinnuclearenergy2003,42(1):25-53.)和使用图形处理单元(graphicsprocessingunit,gpu)并行处理技术。但同样的，这两种改善的方法均没有解决被动式模拟的根源问题。再而，减方差的技巧只有用户给出合适的粒子作用点重要性分布，才有提高计算效率的可能，否则结果会变更差，而给出合适的重要性分布通常是最为困难的。正是因为被动式模拟的根源问题没有得到解决，对如乳腺cbct和脑出血cbct等对散射精度要求更高的散射模拟任务是不可能的。

因此，针对现有技术不足，提供一种基于任务驱动的蒙特卡洛散射光子模拟方法以解决现有技术不足甚为必要。

技术实现要素：

本发明的目的在于避免现有技术的不足之处而提供一种基于任务驱动的蒙特卡洛散射光子模拟方法，该基于任务驱动的蒙特卡洛散射光子模拟方法，从抽样原理的层面提升光子散射模拟效率，主动把模拟任务需求纳入可控制的路径变异空间中，通过自动抽样模型对光子路径进行变异抽样，然后比较前后抽样路径的重要性，选择对模拟任务相对重要的路径作为当前能量沉积路径，从而实现光子路径的自动重要性采样。任务驱动的路径抽样方法摒弃了光子的独立抽样策略，且能主动式调控散射任务，使得模拟效率大大提高。

本发明的上述目的通过如下技术手段实现。

提供一种基于任务驱动的蒙特卡洛散射光子模拟方法，具体步骤如下：

第一步：对光子散射参数进行预设和初始化，生成一条原始光子路径；

第二步：采用均匀抽样算法和随机游走抽样算法对散射点进行位置抽样，产生一条模拟光子路径；

第三步：比对原始光子路径和模拟光子路径的概率；

第四步：应用路径抽样算法对模拟光子路径进行自动抽样；

第五步：将模拟光子能量沉积到对应的探测器中，并判断沉积的光子路径数量是否满足预设的光子路径数量，若是，则结束操作，否则返回第二步。

具体而言的，光子散射的参数包括光子散射能量、类型、阶数、路径概率和光子路径数量，且光子散射路径的总概率为所有片段中每一个片段概率的乘积。

优选的，第三步光子路径概率的算法具体如下：

设光子从源点至da探测器像素中，历经a1、…、ai、…、an的n个散射点的散射作用，经过n+1个片段，每一个片段记为lk(k＝1,2,…,n+1)，对应片段的长度为sk，散射点ai对应瑞利散射和康普顿散射中的任意一种散射类型tj(j＝1或2)，t1为瑞利散射，t2为康普顿散射，散射点ai在xi位置发生康普顿散射或瑞利散射的线性衰减系数为散射点ai在xi位置发生康普顿散射和瑞利散射的总线性衰减系数为μ(xi)；

当k＝1时，光子经过初始片段l1的概率函数p1的求解方法如下：

其中，s1为片段l1的长度，光子分布的方向概率函数f(x)，根据郎伯比尔定律，光子到达x1点的概率函数为光子在x1点发生tj类型的散射作用的概率函数为累乘三种概率函数即为初始片段l1的概率函数，如式(1)所示；

当1＜k≤n时，光子经过中间片段l2～ln中的任意一个片段的概率函数如下：

其中，为ai、ai+1两点所成直线方向的散射微分截面系数，且tm∈tj，光子从散射点ai运动到ai+1且发生tj类型相互作用的概率为

当k＝n+1时，光子经过结尾片段ln+1的概率函数pn+1的求解方法如下：

其中，α为ln+1片段的输运方向和探测器的法线方向夹角，光子在an点发生散射偏转到da探测器元的微分散射截面系数为探测器像素da单位元对应an点的立体角为光子从an到达探测器da位置的概率为

光子完整路径的总概率函数如下所示：

具体而言的，第四步的具体操作如下：

s1、根据式(4)分别计算原始光子路径的概率函数ptotal和模拟光子路径的概率函数p’total；

s2、根据式(5)判断接受概率paccept的值，若paccept≥1，则概率函数p’total所在路径为最终路径，并进入第五步；

若paccept＜1，则进入步骤s3；

s3、从[0，1]中抽样一个随机数ξ，若ξ＜paccept，则概率函数p’total所在路径为最终路径并进入第五步；

否则，定义概率函数p’total所在路径为原始光子路径，并进入第五步。

进一步的，概率函数paccept的计算步骤如下：

其中，q(x)和q(x)’分别对应ptotal和p’total所在光子路径的路径转移概率函数。

具体而言的，均匀概率密度函数或高斯概率密度函数为对称函数，所以式(5)可简化为：

进一步的，第五步分别记录康普顿散射和瑞利散射信号的光子路径模拟数。

优选的，在第二步中，根据散射点位置不同选取对应抽样算法进行随机抽样操作，具体如下：

p1、判断散射点的位置；

若为一阶散射点，则采用均匀抽样算法进行随机选取；

否则，进入步骤p2；

p2、若为多阶散射点，则判断散射点发生的散射类型；

若发生的是康普顿散射，则采用均匀抽样算法进行随机选取；

若发生的是瑞利散射，则采用rws算法进行随机选取。

进一步的，对于多阶散射点中，采用均匀抽样算法的方法具体如下：

在360度范围内随机抽取一个散射角θc，在散射角θc的开口方向上选取一点作为康普顿二阶散射点。

具体而言的，对于多阶散射点中，采用rws算法的方法具体如下：

k1、判断概率函数p’total所在路径为最终路径；

若概率函数p’total所在路径为最终路径，则通过在当前散射角上进行小幅度的高斯偏移，获得新的散射角；

否则，进入步骤k2；

k2、随机一个散射角θr，以σ²为方差进行高斯抽样，获得新的散射角θr’，在θr’方向上均匀选取一点作为新的散射点。

本发明从抽样原理的层面提升光子散射模拟效率，主动把模拟任务需求纳入可控制的路径变异空间中，通过自动抽样模型对光子路径进行变异抽样，然后比较前后抽样路径的重要性，选择对模拟任务相对重要的路径作为当前能量沉积路径，从而实现光子路径的自动重要性采样。任务驱动的路径抽样方法摒弃了光子的独立抽样策略，且能主动式调控散射任务，使得模拟效率大大提高。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1是本发明一种基于任务驱动的蒙特卡洛散射光子模拟方法的系统框图。

图2为圆轨道cbct扫描示意图。

图3为分段显示的光子路径图。

图4为经过路径变异后产生的三条光子路径的示意图。

图5为自动路径抽样算法的部分代码示意图。

图6为实施例2中采用本发明方法模拟出的光子路径与gmcdrr方法模拟出的光子路径的对比图及实施例2中采用本发明方法和gmcdrr方法模拟出的效果对比图。

具体实施方式

结合以下实施例对本发明作进一步描述。

实施例1。

如图1-5所示，一种基于任务驱动的蒙特卡洛散射光子模拟方法，具体步骤如下：

第一步：对光子散射参数进行预设和初始化，生成一条原始光子路径。

其中，光子散射的参数包括光子散射能量、类型、阶数、路径概率和光子路径数量，且光子散射路径的总概率为所有片段中每一个片段概率的乘积。

第二步：采用均匀抽样算法和随机游走抽样算法对散射点进行位置抽样，产生一条模拟光子路径。

在第二步中，根据散射点位置不同选取对应抽样算法进行随机抽样操作，具体如下：

p1、判断散射点的位置。

若为一阶散射点，则采用均匀抽样算法进行随机选取。

否则，进入步骤p2。

p2、若为多阶散射点，则判断散射点发生的散射类型。

若发生的是康普顿散射，则采用均匀抽样算法进行随机选取。

对于多阶散射点中，采用均匀抽样算法的方法具体如下：

在360度范围内随机抽取一个散射角θc，在散射角θc的开口方向上选取一点作为康普顿二阶散射点。

若发生的是瑞利散射，则采用rws算法进行随机选取。

对于多阶散射点中，采用rws算法的方法具体如下：

k1、判断概率函数p’total所在路径为最终路径。

若概率函数p’total所在路径为最终路径，则通过在当前散射角上进行小幅度的高斯偏移，获得新的散射角。

否则，进入步骤k2。

k2、随机一个散射角θr，以σ²为方差进行高斯抽样，获得新的散射角θr’，在θr’方向上均匀选取一点作为新的散射点。

第三步：比对原始光子路径和模拟光子路径的概率。

第三步光子路径概率的算法具体如下：

设光子从源点至da探测器像素中，历经a1、…、ai、…、an的n个散射点的散射作用，经过n+1个片段，每一个片段记为lk(k＝1,2,…,n+1)，对应片段的长度为sk，散射点ai对应瑞利散射和康普顿散射中的任意一种散射类型tj(j＝1或2)，t1为瑞利散射，t2为康普顿散射，散射点ai在xi位置发生康普顿散射或瑞利散射的线性衰减系数为散射点ai在xi位置发生康普顿散射和瑞利散射的总线性衰减系数为μ(xi)。

当k＝1时，光子经过初始片段l1的概率函数p1的求解方法如下：

其中，s1为片段l1的长度，光子分布的方向概率函数f(x)，根据郎伯比尔定律，光子到达x1点的概率函数为光子在x1点发生tj类型的散射作用的概率函数为累乘三种概率函数即为初始片段l1的概率函数，如式(1)所示。

当1＜k≤n时，光子经过中间片段l2～ln中的任意一个片段的概率函数如下：

其中，为ai、ai+1两点所成直线方向的散射微分截面系数，且tm∈tj，光子从散射点ai运动到ai+1且发生tj类型相互作用的概率为

当k＝n+1时，光子经过结尾片段ln+1的概率函数pn+1的求解方法如下：

其中，α为ln+1片段的输运方向和探测器的法线方向夹角，光子在an点发生散射偏转到da探测器元的微分散射截面系数为探测器像素da单位元对应an点的立体角为光子从an到达探测器da位置的概率为

光子完整路径的总概率函数如下所示：

第四步：应用路径抽样算法对模拟光子路径进行自动抽样。

第四步的具体操作如下：

s1、根据式(4)分别计算原始光子路径的概率函数ptotal和模拟光子路径的概率函数p’total。

s2、根据式(5)判断接受概率paccept的值，若paccept≥1，则概率函数p’total所在路径为最终路径，并进入第五步。

若paccept＜1，则进入步骤s3。

s3、从[0，1]中抽样一个随机数ξ，若ξ＜paccept，则概率函数p’total所在路径为最终路径并进入第五步。

否则，定义概率函数p’total所在路径为原始光子路径，并进入第五步。

其中，概率函数paccept的计算步骤如下：

其中，q(x)和q(x)’分别对应ptotal和p’total所在光子路径的路径转移概率函数。

均匀概率密度函数或高斯概率密度函数为对称函数，所以式(5)可简化为：

具体算法如图(5)所示。

第五步：将模拟光子能量沉积到对应的探测器中，并判断沉积的光子路径数量是否满足预设的光子路径数量，若是，则结束操作，否则返回第二步。

第五步分别记录康普顿散射和瑞利散射信号的光子路径模拟数。

从抽样原理的层面提升光子散射模拟效率，主动把模拟任务需求纳入可控制的路径变异空间中，通过自动抽样模型对光子路径进行变异抽样，然后比较前后抽样路径的重要性，选择对模拟任务相对重要的路径作为当前能量沉积路径，从而实现光子路径的自动重要性采样。任务驱动的路径抽样方法摒弃了光子的独立抽样策略，且能主动式调控散射任务，使得模拟效率大大提高。

实施例2。

一种基于任务驱动的蒙特卡洛散射光子模拟方法，其它特征与实施例1相同，不同之处在于：如图6所示，该实施例为实际模拟操作，以下所采用和获取的数据均为实际操作所得，具有真实性。

该模拟操作的所在平台为ubuntu-12.04.4系统gpu架构，显卡设备为nvidiageforcegtxtitanz，x射线源能量为60kvp，所使用的均匀体模大小为10×10×2.8cm³，矩阵大小为250×250×280，探测器的大小为40cm×30cm，矩阵大小为512×384，射线源到旋转中心距离和旋转中心到探测器的距离分别为15.59cm和49.41cm，以下为具体的模拟步骤：

第一步，进行光子散射能量，类型，阶数，路径概率等项目的预设和初始化，产生原始路径x。

光子能量在调试阶段设置为单能60kvp，散射类型主要为康普顿散射和瑞利散射两种，设置一个类型控制阈值controltype＝0.5。

然后在[0,1]中均匀抽取随机数ζ，若ζ<0.5，则认为该点发生康普顿散射，否则认为发生瑞利散射。

本发明中为了防止相互干扰，不同散射阶数分开模拟。

第二步：采用均匀抽样算法和随机游走抽样算法(random-walksampling，rws)进行散射点位置抽样，产生新的光子路径y。

对于一阶散射点，无论是康普顿散射还是瑞利散射，均采用均匀抽样策略。

即在均匀体模内随机选取一个作用点作为散射作用点。对于二阶散射点，康普顿散射和瑞利散射分别采取不同的抽样方法。

对康普顿散射采用均匀抽样算法，首先在360度范围内随机抽取一个散射角θc，然后再在此角度方向上均匀选取一点作为康普顿二阶散射点。

对瑞利散射采取rws，即如果接受当前路径，则在当前散射角下进行小幅度的高斯偏移，从而获得新的散射角。

具体为，随机初始化一个散射角θr，以此θr为均值，以σ²为方差进行高斯抽样(σ²的选择任意，但一般取比较小的值)，获得新的散射角θr’，并在θr’方向上均匀选取一点作为瑞利二阶散射点。

第三步：计算每个光子片段在扫描物体内的长度lk，查询线性衰减系数μ(x)以计算光子路径积分查询微分散射截面数据并计算立体角根据公式(1)、(2)、(3)计算p1、pk和pn+1。最后根据公式(4)计算原始光子路径的概率函数ptotal和模拟光子路径的概率函数p’total。

第四步：应用路径抽样算法自动抽样光子路径，由于本发明中选取的路径转移概率函数为均匀概率密度函数或高斯概率密度函数，这些概率密度函数是对称的，所以q(x)＝q(x)’，式(5)简化后如下：

判断接受概率paccept的值，若paccept≥1，则概率函数p’total所在路径为最终路径，并进入第五步。

若paccept＜1，则进入步骤s3。

s3、从[0，1]中抽样一个随机数ξ，若ξ＜paccept，则概率函数p’total所在路径为最终路径并进入第五步。

否则，定义概率函数p’total所在路径为原始光子路径，并进入第五步。

第五步：沉积路径中光子能量到对应的探测器像素点，分别记录康普顿散射和瑞利散射信号，并判断是否完成所需的路径模拟数。

若是，则结束模拟，否则返回第二步。

仿真结果如下：

图6中的图(a)-(d)分别为本方法所模拟出的一阶康普顿散射，一阶瑞利散射，二阶康普顿散射和二阶瑞利散射信号。

图6中的图(e)-(h)为gmcdrr所模拟出的一阶康普顿散射，一阶瑞利散射，二阶康普顿散射和二阶瑞利散射信号。

图6中的图(i)-(l)对应为图6中的图(a)-(h)图中黄色参考线的轮廓线图。

对比实验均是在同一个不确定水平δ(uncertaintylevel)下比较的，其中n为探测器像素数，stdi为第i个像素点的信号标准差。

从图(a)-图(i)中可见，本方法模拟的结果与gmcdrr模拟的结果基本一致，也就是说本方法可以达到与gmcdrr相同的模拟水平。

其次，还统计了两种方法的模拟时间，如表(1)所示：

表(1)

从表可见，本方法在模拟结果相同的条件下，模拟的速度明显要优于gmcdrr方法。

对于一阶散射来说，方法加速了大约208倍；对于二阶散射来说，加速了大约12倍。

这说明了从抽样原理层面提高mc散射模拟效率是可行的。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。