基于拟蒙特卡罗和强制探测的校正CT散射信号的方法

技术领域

本发明涉及计算机断层成像技术(CT)领域，特别是关于一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法。

背景技术

计算机断层成像技术(Computed Tomography，CT)以其无损检测、图像不重叠和高分辨率等优势，广泛应用于材料科学、无损检测、医学等诸多领域，尤其在医学诊断，彻底改变了放射诊断学。

CT图像重建是指通过测量X射线穿过被测模体M后的衰减数据，重建出被测模体M的线性衰减系数μ

I(l)＝I

得到的μ

散射信号会在测量得到的μ

因此，为了能够提高线性衰减系数μ

发明内容

本发明的目的在于提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，来克服或至少减轻现有技术的上述缺陷中的至少一个。

为实现上述目的，本发明提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，包括：

步骤10，模拟光子从X射线源S通过被测物体M到达探测器D的散射信号；包括：

步骤11，通过沃克混叠方法，用预先生成的4n维Sobol'点u＝(u

步骤12，记初始入射光到探测器D的立体角为Ω

步骤13，在确定第i-1阶相互作用点A

其中

光子以权重

步骤14，基于下式确定下一个相互作用点A

其中s

步骤15，确定光子在A

其中，

步骤16，获得光子沿路径l

P(l

光子沿路径S→A

P(l

其中，光子运动路径S→A

步骤17，获得光子在被测物体M内发生n次碰撞后，固定探测器像素D

所述方法还包括步骤20，将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。

在一种实施方式中，所述线性衰减系数μ

本发明由于采取以上技术方案，其具有以下优点：

通过计算光子理论上发生散射的概率，把QMC方法的优势和强制固定探测(FFD)技术结合起来，能够准确地得到探测器探测到的光子的散射信号，从而实现对探测器探测到的X射线的散射信号的校正。

附图说明

图1为本发明实施例提供的基于拟蒙特卡罗和强制探测的校正CT散射信号的方法的流程示意图。

图2为本发明实施例提供的光子从X射线源S射出经过被侧物体M到达探测器D的过程示意图。

图3示出本发明实施例提供的光子的能谱分布示意图。

图4示出均匀铝模体和人头模体的几何示意图。

图5示出对均匀铝模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。

图6示出对人头模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。

具体实施方式

在附图中，使用相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面结合附图对本发明的实施例进行详细说明。

在本发明的描述中，术语“中心”、“纵向”、“横向”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明保护范围的限制。

在不冲突的情况下，本发明各实施例及各实施方式中的技术特征可以相互组合，并不局限于该技术特征所在的实施例或实施方式中。

下面结合附图以及具体实施例对本发明做进一步的说明，需要指出的是，下面仅以一种最优化的技术方案对本发明的技术方案以及设计原理进行详细阐述，但本发明的保护范围并不仅限于此。

本文涉及下列术语，为便于理解，对其含义说明如下。本领域技术人员应当理解，下列术语也可能有其它名称，但在不脱离其含义的情形下，其它任何名称都应当被认为与本文所列术语一致。

本发明实施例提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，图1示出该方法的流程示意图，包括：

步骤10，模拟光子从X射线源S通过被测物体M到达探测器D的散射信号。

图2是光子从X射线源S射出经过被侧物体M到达探测器D的过程示意图。其中，该散射信号的模拟可以是针对信号，也可以是针对信号的其它参数，本文对此不做限制。

该图示仅为了便于理解本发明，不用于其它限定。

上述步骤10具体包括下述步骤11-17。

步骤11，通过沃克混叠方法，用预先生成的4n维Sobol'点u＝(u

其中，能谱分布φ(E

步骤12，记初始入射光到探测器D的立体角为Ω

本步骤中，设Ω

步骤13，在确定第i-1阶相互作用点A

其中，

光子以权重

本步骤中，用线性衰减系数μ

记光子从其当前位置A

其中

步骤14，基于下式确定下一个相互作用点A

其中u

容易理解，在确定下一个相互作用点后，可以再根据步骤13的方式确定光子从该下一个相互作用点碰撞后逃出被测物体M的概率。

步骤15，确定光子在A

其中，

本步骤中，光子在A

如图2所示，光子在A

步骤16，获得光子沿路径l

P(l

光子沿路径S→A

P(l

其中，光子运动路径S→A

步骤17，获得光子在被测物体M内发生n次碰撞后，固定探测器像素D

因此，f

其中

步骤20，将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。

校正投影数据中的散射信号，并重建CT图像的过程如下：1，读取投影数据；FDK算法重建得到体数据图像2，使用发明的基于GPU的拟蒙特卡罗方法和强制探测器技术相结合的散射模拟方法，快速模拟出散射信号；3，把每个探测器像素接收到的信号减去估计到的散射信号；4，光束硬化校正；5，使用校正散射污染的投影数据，根据FDK算法重建CT图像。6，重复步骤2-5，直到图像收敛。

在计算光子散射信号所需要的数据(例如几何参数、X射线光谱和物理数据)初始化仿真。所有的这些数据传输到GPU内存中。具体来说，由于GPU的缓存功能，体素化的模体数据和衰减系数被存储在GPU的纹理存储器中以实现快速访问。DCS数据存储在GPU全局内存中，其余变量存储在GPU常量内存中。所有图像计数器都设置为零。之后，以批处理方式进行仿真，进而使用批处理之间的结果来估计不确定性。

本发明实施例提供的上述方法可以简化为如下步骤：

1.赋n具体值，预先生成4n维Sobol’点u＝(u

2.用u

3.如果i≥n，跳到步骤5；否则令i＝i+1，用u

4.用u

5.令

6.返回P(l

为了说明本发明实施例提供的校正方法的效果，下面提供实验结果。

分别计算X射线穿过均匀铝模体和人头模体后到达探测器D的散射信号。实验所用的铝模体体积为160×28×160mm

通过计算相对差异(RD)和有效提高因子(EIF)来说明gQMCFFD算法的有效性和稳健性。

RD＝||r-t||

其中r是由gMCFRD算法计算探测器D接收到的散射信号，t可代表由gQMCFFD算法、gMCFFD算法或gQMCFRD算法计算探测器D接收到的散射信号。与gQMCFFD算法比较，gMCFFD算法用随机序列代替Sobol’序列进行抽样。gMCFRD算法基于GPU的MC方法和FRD技术相结合、gQMCFRD算法是基于GPU的QMC方法和FRD技术相结合,详见参考文献[14]。

其中

下述表1和表2分别显示了对均匀铝模体和人头模体，gMCFRD运行2

表1示出关于均匀铝模体，gQMCFRD、gMCFFD、gQMCFFD方法计算的散射信号与gMCFRD比较，其中，N表示模拟次数，T表示运行时间，RD表示其余三种方法的计算结果与gMCFRD方法的计算结果的相对差异，FOM等于运行时间和方差乘积的倒数，EIF为有效提高因子。

表2示出关于人头模体，gQMCFRD、gMCFFD、gQMCFFD方法计算的散射信号与gMCFRD比较，其中，N表示模拟次数，T表示运行时间，RD表示其余三种方法的计算结果与gMCFRD方法的计算结果的相对差异，FOM等于运行时间和方差乘积的倒数，EIF为有效提高因子。

表1

表2

图5和图6分别展示了对均匀铝模体和人头模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。图像视觉上的高度一致性说明了gQMCFFD方法的准确性。图5(a)-(d)分别示出gMCFRD、gQMCFRD、gMCFFD和gQMCFFD对Al模体计算的总散射信号结果，(e)对应上述4种方法，探测器测器中心且平行于x轴位置接收到的散射信号曲线图；(f)示出主信号。图6(a)-(d)分别示出gMCFRD、gQMCFRD、gMCFFD和gQMCFFD对人头模体计算的总散射信号结果，(e)对应上述4种方法，探测器测器中心且平行于x轴位置接收到的散射信号曲线图，(f)示出主信号。

通过采用本发明实施例提供的方法，计算光子理论上发生散射的概率，用公式(4)抽样下一阶相互作用点，这不同于已有的基于MC的散射校正方法。gQMCFFD算法使用Sobol’序列抽样，并在光子与被测物体发生碰撞后，计算其到达每个固定探测器像素的理论上的概率，把QMC方法的优势和强制固定探测(FFD)技术结合起来，能够准确地得到探测器探测到的光子的散射信号。

最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。