一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法

摘要

本发明公开了一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，该方法包括：步骤1，获取被测物体在多个能谱下的多色投影数据；步骤2，获取各种基材料密度图像估计值；步骤3，对基材料密度图像估计值进行正投影，得到各种基材料的投影估计值及被测物体的多色投影估计值；步骤4，计算多色投影估计值与多色投影数据之间的误差，利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差；步骤5，将投影残差进行反投影，得到各种基材料的残差图像，更新被测物体的基材料密度图像估计值；步骤6，判断终止条件是否满足，若满足则终止迭代，否则转向步骤3。本发明能够由采集到的多个能谱的多色投影数据快速重建出被测物体的多种基材料密度图像。

说明书

技术领域

本发明涉及X射线CT成像技术领域，特别是关于一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法。

背景技术

X射线计算机断层成像技术(X ray Computed Tomography,简称X射线CT)可以在不破坏或损伤物体的情况下呈现物体的内部细节，已广泛应用于医学、生物、工业、材料、古化石和航天等众多领域。传统CT成像理论假设X射线由单一能量的光子组成，忽略了X射线的多色性，因此利用传统单能量CT重建算法重建实采数据时会产生射束硬化伪影，如杯状伪影和条状伪影，严重影响成像质量。

多能谱CT成像系统利用多个不同能谱下的X射线扫描被测物体，获得被测物体在多个不同能谱下的投影数据。利用这些投影数据可以重建被测物体的基材料密度图像，或者是等效原子序数和电子密度图像。与单能谱CT相比，多能谱CT获取了更多的被测物体的信息，具有更好的物质区分能力，在硬化伪影去除、骨密度测量、PET衰减校正和伪单能图像计算等方面有广泛的应用前景。

现有的多能谱CT多色投影数据获取方式大致有两类，第一类方式是利用X射线源获取两个或多个X射线能谱，分别对被测物体进行扫描，获取两组或多组多色投影数据。代表性的技术有“全扫描”模式、双源双探扫描模式和快速电压切换扫描模式等。其中，“全扫描”模式是利用传统CT设备，在不同的管电压和管电流下对被测物体分别进行多次扫描，这种方法可在传统CT设备上完成，不需要额外添加硬件设备。第二类方式是只使用一个X射线能谱，使用三明治式探测器或者光子计数探测器获取多组多色投影数据。通常情况下，第一类方式获取的多组多色投影数据是射线路径几何不一致的，第二类方式获取的多组多色投影数据是射线路径几何一致的。

重建算法大致可以分为四类：图像域重建方法、投影域重建方法、基于深度学习的重建方法和迭代重建方法。图像域重建方法是对采集到的多组多色投影数据分别用传统CT重建算法进行重建，然后将重建图像进行线性组合得到基材料的密度图像。该方法重建的图像被认为是对真实图像的一阶近似，无法准确描述基材料密度图像与多色投影之间的非线性关系，通常在重建图像中还会存在硬化伪影，且受噪声的影响较大。为了提高图像域重建算法的图像质量，研究者提出了一些基于传统低通滤波器或者基于统计先验滤波器的改进算法，这些算法可以在一定程度上抑制噪声对重建结果的影响，然而对于图像分解精度的改进有限。投影域的重建方法首先对多组多色投影进行组合，得到多种基材料密度图像的投影数据，然后对多种基材料密度图像的投影数据分别利用传统CT重建方法进行重建，得到基材料密度图像。通常，利用投影域重建方法重建的图像要优于利用图像域重建方法重建的图像。但是该方法要求不同能谱下采集的投影数据是几何一致的，即要求每条X射线路径下的所有能谱的投影数据都要采集到。然而实际的多能谱CT系统采集的数据，并不能保证这一要求，比如双源双探的能谱CT扫描模式。近年来，利用深度学习技术，研究者在有完备训练集的情况下，在多能谱CT重建领域得到了高质量的图像重建结果。然而很多情况下，本实施例无法获得充足的训练样本。考虑到多能谱CT问题的非线性性，理论上迭代重建算法更适合这类问题求解，其利用数值方法或者优化方法构造迭代结构，通过对图像重建结果逐步修正，可以得到高精度的被测物体的各基材料密度图像信息。由于能谱CT重建问题的非线性性以及病态性，使得现有的多能谱CT重建算法无法快速重建出高质量的基材料密度图像。

发明内容

本发明的目的在于提供一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其能够由采集到的多个能谱的多色投影数据快速重建出被测物体的多种基材料密度图像。

为实现上述目的，本发明提供一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，该方法包括：步骤1，利用多能谱CT系统扫描被测物体，获得被测物体在多个能谱下的真实多色投影数据；步骤2，为被测物体的各种基材料密度图像赋初值，作为各种基材料密度图像估计值；步骤3，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值，进而利用X射线能谱信息和基材料的质量衰减系数信息得到被测物体在各个能谱下的多色投影估计值；步骤4，计算多色投影估计值与真实多色投影数据之间的误差，并利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差A；步骤5，将各种基材料的投影残差进行反投影操作，得到各种基材料的残差图像，进而更新被测物体的各种基材料密度图像估计值；步骤6，判断终止条件是否满足，若满足则终止迭代，否则转向步骤3。

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

本发明方法用于多能谱CT快速迭代重建，适用于多种常用的多能谱CT扫描模式。与现有方法相比，本发明方法能够在保证多能谱CT图像重建质量的同时，显著加快重建图像的收敛速度。

附图说明

图1为本发明提出的一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法的流程图。

图2a为本发明一个实施例用作测试模体的FORBILD胸腔模型图像。

图2b为本发明一个实施例用作测试模体的水基材料密度图像。

图2c为本发明一个实施例用作测试模体的骨基材料密度图像。

图3为本发明一个实施例所用的X射线源分别在140kV和80kV管电压下发出的高能谱和低能谱示意图。

图4a为本发明一个实施例的140kV能谱下采集到的测试模体添加初始光子数为10

图4b为本发明一个实施例的80kV能谱下采集到的测试模体添加初始光子数为10

图5a为一个实施例采用本发明方法的骨基材料密度图像重建的结果图像。

图5b为一个实施例采用本发明方法的水基材料密度图像重建的结果图像。

图5c为本发明一个实施例所重建的关于能量为70keV的光子的线性衰减系数图像。

图6为一个实施例，用以说明本发明方法与已公开的相关专利文献方法的不同。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

参阅图1，为本发明实施例提供的基于倾斜投影修正技术的多能谱CT快速迭代重建方法包括以下步骤：

步骤1，利用多能谱CT系统扫描被测物体，获得被测物体在N个能谱下的N个真实多色投影数据，其中，一个能谱下相应获得一个多色投影数据。

步骤2，为被测物体的各种基材料密度图像赋初值，作为各种基材料密度图像估计值。

步骤3，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值，进而利用X射线能谱信息和基材料的质量衰减系数信息得到被测物体在各个能谱下的多色投影估计值。

步骤4，计算多色投影估计值与真实多色投影数据之间的误差，并利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差。

步骤5，将各种基材料的投影残差进行反投影操作，得到各种基材料的残差图像，进而更新被测物体的各种基材料密度图像估计值。

步骤6，重复步骤3至步骤5，直到满足终止条件。

在一个实施例中，步骤1中的“多能谱CT系统”大致有两类。

第一类方式是利用X射线源获取两个或多个X射线能谱，分别对被测物体进行扫描，获取两组或多组多色投影数据。代表性的技术有“全扫描”模式、双源双探扫描模式，快速电压切换扫描模式等。其中，“全扫描”模式是利用传统CT设备，在不同的管电压和管电流下对被测物体分别进行多次扫描，这种方法可在传统CT设备上完成，不需要额外添加硬件设备。

第二类方式是只使用一个X射线能谱，使用三明治式探测器或者光子计数探测器获取多组多色投影数据。通常情况下，第一类方式获取的多组多色投影数据是射线路径几何不一致的，第二类方式获取的多组多色投影数据是射线路径几何一致的。

在一个实施例中，步骤1中，根据多能谱CT系统中所用探测器不同，得到的X射线多色投影数据维度也不同。

当探测器为包含W个探测器单元的线阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×V个X射线投影数据。其中，每个角度的投影包括W个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅二维X射线多色投影图。

当探测器为包含W×H个探测器单元的面阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×H×V个X射线投影数据。其中，每个角度的投影包括W×H个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅三维X射线多色投影图。

在一个实施例中，步骤2中的“为被测物体的各种基材料密度图像赋初值”是以0作为初始图像，即图像中的每个像素均为0。

在一个实施例中，以

式中，l表示一条X射线扫描路径；R

在一个实施例中，以p

式中，S

通常，若多能谱CT系统获取的多色投影数据是几何一致的，那么在ξ

在一个实施例中，以多色投影射线路径几何不一致为例，说明步骤4中的“利用倾斜投影技术对步骤3中的基材料的投影估计值

步骤4.1，将式(2)表示的各个能谱下的多色投影估计值p

假设经过n次迭代后，已有基材料密度图像估计值(f

其中，

式中，

其中，S

在几何不一致的情况下，在某些X射线路径下只已知一个多色投影数据，假定当前已知的第1个能谱对应的多色投影估计值p

其中：

i表示第i个x射线能谱对应的数据，i＝1…N；j表示第j种基材料密度图像对应的数据，j＝1…N；f

步骤4.2，获取当前迭代点向超平面做正交投影的投影方向：利用线性多色投影方程(6.1)，计算当前迭代点(R

式中，dir21表示当前迭代点(R

取式(8)的单位向量作为当前射线路径下已知能谱下的线性多色投影方程(6.1)对应的超平面做正交投影的投影方向dir2即如式(9)：

其中，||·||为向量取模运算符；dir2表示当前射线路径下下已知能谱对应的投影方程(6.1)所表示的超平面的法方向；dir21表示式(8)的结果。

步骤4.3，计算由线性多色投影方程组(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向：

由线性多色投影方程(6.2)～(6.N)可知，向量(a

其中：|·|为行列式；e

由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向可以表示为式(11)：

dir11＝(D

dir11表和dir12均表示由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向，dir11和dir12为求解过程中的中间变量，最终使用的是dir1。

取式(11)所表示的两个方向中与dir2成锐角的方向，作为由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向dir1，即式(12)：

其中，<，>表示两个向量的内积。

步骤4.4，利用式(13)，将dir1和dir2两个方向加权求和，得到从当前迭代点(R

dir＝λ

为从当前迭代点(R

沿该方向dir在当前射线路径下已知的多色投影方程上的投影点(R

(R

式中，A为各种基材料的投影残差，即当前迭代点(R

在一个实施例中，步骤5中，根据式(16)获得被测物体的各基材料密度图像的估计值：

(f

式中，T表示转置。

得到本发明实施例更新后的基材料密度图像估计值f

需要说明的是，在多能谱CT重建问题中，通常多色投影射线路径几何不一致时有效的方法，可以直接推广到几何一致时使用，所以，几何一致下的步骤和几何不一致下的步骤是一样的。

在一个实施例中，步骤6所述的终止条件，包括但不限于最大迭代次数，得到的被测物体的各种基材料密度图像估计值是否符合主观判断或所需要求，两次迭代得到的各种基材料密度图像估计值的差小于设定的阈值。

在一个实施例中，无论是步骤3所述将各种基材料密度图像估计值进行正投影得到各种基材料的投影估计值，还是步骤4所述由多色投影估计值与真实多色投影数据之间的误差利用倾斜投影修正技术得到各种基材料的投影残差，或是步骤5所述将各种基材料的投影残差进行反投影操作得到各种基材料的图像残差，均可采用并行方法计算，并基于硬件并行计算平台加速实现。

以下通过一个具体实施例来说明本发明提出的一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法的具体实现过程。

图2a为本具体实施例中用作测试模体的FORBILD胸腔模型图像。图2b为本具体实施例中用作测试模体的水基材料密度图像。图2c为本具体实施例中用作测试模体的骨基材料密度图像。在本实施例中使用的高低能X射线能谱用开源软件使用开源软件SpectrumGUI，仿真了GE Maxiray 125球管在80kVp和140kVp下的能谱，并对能谱做了归一化处理。实验管电压分别为80kV和140kV。在模拟实验中，模体假定为不同密度的水和骨的组合。水和骨也是重建图像的基材料。材料的质量衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。在模拟多色投影时，X射线能谱和材料的质量衰减系数的采样间隔均为1kV。

实验扫描参数设置如下：射线源到转台中心的距离为490mm，射线源到探测器的距离为880mm。线探测器由512个探测器单元组成，每个探测器单元的尺寸为0.2mm。

根据方程(2)模拟多色投影，某一射线路径下X射线与基材料密度图像的交线长由光线投射算法计算。无噪声模拟实验和含噪声模拟实验均以快速电压切换的形式采集得到，即偶数多色投影是在某一种电压下模拟扫描得到，奇数多色投影是在另一种电压下模拟扫描得到，因此得到的投影为几何不一致的投影。两个电压下的多色投影，均在CT系统旋转一周的情况下采集了720个角度的投影数据。对投影数据添加初始光子数为10

具体实施步骤如下：

步骤1：为基材料密度图像f

步骤2：假定经过n次迭代后，有基材料密度图像估计值f

步骤3：首先根据公式(8)和(9)计算出dir2；然后根据公式(10)、(11)和(12)计算出dir1；最后根据公式(13)计算出本文算法的倾斜投影方向dir；

步骤4：根据公式(14)计算当前路径下更新后的基材料密度图像投影R

步骤5：判断终止条件是否满足，若满足则终止迭代，否则转向步骤2。

选取参数λ

以下通过一个具体实施例来说明本发明方法与已公开的相关专利文献的不同。

使用开源软件Spectrum GUI生成，仿真了GE Maxiray 125球管140kVp下的能谱，并对能谱做了归一化处理。将140kV下的能谱分为1-70kV和71-140kV两段，能谱的采样间隔δ为1kV，并分别归一化，分别作为本实施例中使用的X射线能谱

根据(2)式，本实施例构建在一条特定射线路径下的二元非线性方程组

其中，

图6为用E-ART方法和本发明方法求解该方程组时的迭代轨迹图，黑色实线和黑色虚线分别为S

可以看出，本发明方法与上述已公开的相关专利文献方法的不同之处在于E-ART方法以正交投影的方式更新迭代点，迭代点收敛到真解的速度很慢，迭代50次时仍不能收敛至真解附近。而本发明方法采用倾斜投影技术，大大加快了迭代点到真解的收敛速度，20次左右即收敛到真解。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

本领域普通技术人员可以理解：本实施例采用快速电压切换模式获得所需数据，但是本发明保护范围并不限于此扫描模式。

本领域普通技术人员可以理解：本实施例获得的两组数据是几何不一致的，但是本发明可以对几何一致的数据进行，获得相应的基材料密度图像。

本领域普通技术人员可以理解：本发明可以对三个或更多个能谱下的扫描数据进行，获得三个或更多种基材料的密度图像。

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