[发明专利]X射线蒙特卡洛模拟的路径积分方法

说明书

技术领域

本发明涉及医学图像处理技术，尤其涉及一种X射线蒙特卡洛模拟的路径积分方法。 

背景技术

对简单几何模型的散射数值模拟是X-ray(X射线)散射矫正技术的基础，数值模拟仿真是检验散射矫正技术优劣的必不可少的手段，高效率的数值模拟仿真算法不仅可以加快散射矫正技术的研发工作，本身也可以直接发展为一种散射矫正技术。 

在已有的蒙特卡洛X射线散射模拟中，人们采用一系列伪随机数或准随机序列确定光子的出射扇角、张角和透射深度，并跟踪每一束光子到达检测器格点的概率并加以积分。如文献“An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT”，G.Poludniowski，P.M.Evans，V.N.Hansen and S.Webb，Phys.Med.Biol.54(2009)3847-3864.描述的CRFD方法。 

在这一方法中，随机产生的出射角度和透射深度不能确保每个随机光子的散射点落在模体内部，这些无效散射点导致一定的额外计算开销，并且在模体边界处被积函数不连续导致一定的误差，不能充分利用准随机序列的低差异性质。这些误差需要提供更多数量的积分点加以弥补，从而导致计算量增大。 

因此，有必要提供一种X射线蒙特卡洛模拟的路径积分方法，使得积分区域等同于模体的有效区域，被积函数在积分区域上连续，消除边界上不连续导致的误差，减少计算量。 

发明内容

本发明所要解决的技术问题是提供一种X射线蒙特卡洛模拟的路径积分方法，使得积分区域等同于模体的有效区域，被积函数在积分区域上连续，消除了边界上不连续导致的误差，减少计算量。 

本发明为解决上述技术问题而采用的技术方案是提供一种X射线蒙特卡洛模拟 的路径积分方法，包括如下步骤：a)产生[0，1]³ⁿ空间上长度为m的低差异序列，其中n为散射次数，m为样本个数；b)通过将所述低差异序列进行同构映射变换生成与模体的几何形状匹配的散射点序列；c)由该些散射点序列对每个检测器网格点构成m条光子路径；d)确定经由所述光子路径到达检测器网格点的概率密度；e)将所述概率密度加权求和，得到每个检测器网格点上的散射强度。 

进一步地，所述散射强度的计算公式为 其中dΩ为三维空间的积分微元， 是第i条散射路径，V为模体的体积，m为样本个数，n为散射次数。 

进一步地，所述低差异序列为Halton序列，生成的准随机数范围在0～1之间。 

进一步地，所述模体的几何形状为椭球体或椭圆柱或四面体或平行六面体。 

进一步地，所述散射次数n为0～∞的整数，每个检测器网格点上的总散射强度为所有散射次数分别计算后得到的散射强度之和。 

进一步地，所述每个检测器网格点上的总散射强度为述散射次数n分别取值1、2、3计算后得到的散射强度之和。 

进一步地，所述散射次数n取值为2时的散射强度为两个散射点确定的光子路径概率在R⁶空间上的积分。 

进一步地，所述光子经一次散射到达检测器网格点的概率密度 按如下公式计算： 

其中， 是散射角为θ时的概率密度，P_Compton是发生Compton散射的分概率，L₁是X射线源到散射点的距离，L₂是散射点到检测器格点的距离，D₁是X射线到达散射点之前在模体内穿过的距离，D₂是X射线经过散射之后在模体内穿过的距离，φ_d是X射线在检测器上的入射角度，μ是模体的X射线衰减系数，s_d是检测器的有效面积。 

进一步地，所述光子经n次散射到达检测器网络点的概率密度 按如下公式计算： 

其中，n是散射次数，n为大于1的正整数，θ_i是第i次散射的散射角， 是散射角为θ_i时的概率密度，P_Compton是发生Compton散射的分概率，L_j是光子路径上第j段折线段的长度，D_j是光子路径上第j段折线段在模体中经过的长度，φ_d是X射线在检测器上的入射角度，μ是模体的X射线衰减系数，s_d是检测器的有效面积。