一种基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法

技术领域

本发明属于CT技术图像重建算法设计领域，具体涉及一种基于不完全投影数据的CT图像迭代重建算法。

背景技术

计算机断层成像CT(Computed Tomography)技术是利用X射线测量目标物体各角度下的投影数据以获取目标剖面信息的成像技术。近年来，CT技术在医学和工业上得到越来越广泛的应用。在医疗诊断方面，为临床诊断病症提供了一种无损检测方法。在工业方面，通过CT技术能够检测到物体内部的一系列物理特征，因此其被广泛地应用于材料检测、地质勘探、航天航空等领域。CT重建算法作为CT技术的核心部分，是利用图形重建算法还原出检测目标的断层图像，但由于成像硬件、扫描几何或电离辐射暴露等实际限制，使采集到的投影数据高度稀疏，导致图像重建过程发生不适定问题，从而降低重建图像的质量。不完全投影数据重建所反映的重建质量问题给图像重建算法带来挑战，如何从有限的投影数据中重构高质量CT图像成为图像重建的研究热点。

CT图像重建中的稀疏视图采样是指全扫描范围内的等间隔采样。由于稀疏视图投影数据不满足数据完整性条件，传统的解析重建方法在投影不完全的情况下是无法实现的。以最常用的滤波反投影(Filter Back Projection,FBP)方法为例，当反投影步骤不满足0-180°范围内进行不同角度的高密度投影时，重构后的图像会出现一些角度受限的伪影和一些边缘被破坏的情况。由马继明等人通过对滤波反投影原理及其重建图像与理想CT图像差值关系的分析，构造了以FBP为基础的迭代循环，改进的TV-iFBP算法在抑制受限角度伪影方面优于经典的FBP算法，可以提高重构图像的质量，但在稀疏投影角度小于50°时改进的TV-iFBP算法重建中会观察到几何畸变和条带伪影。

迭代重建法是不完全投影数据重建中另一种经典方法，它可以利用已有的投影进行连续迭代来逼近真实图像。经典的代数重建算法(Algebraic ReconstructionTechnique,ART)、联合代数重建算法(Simultaneous Algebraic ReconstructionTechnique,SART)对不完全投影数据的容错性较好，即使投影数据不满足数据完整性条件也能有效提高重建质量。实际上迭代图像重建过程为求解线性方程组的形式，随着投影缺失角范围的增大，由于缺乏约束方程，使得求解的方程集变得不适定，重构解位于过大的解域中的随机点上。因此当投影视图高度稀疏，重建图像没有引入先验信息时，重建结果会急剧下降，从而无法得到准确的重构结果。

Candes等人提出了压缩感知(Compressed Sensing,CS)理论，将其引入到CT图像重建领域。该理论认为，稀疏图像可以以远低于Shannon/Nyquist的采样频率进行精确重构，为解决不完全投影数据重建提供了新思路，有效地解决稀疏视图CT重建等不适定逆问题。传统的TV正则化方法通过结合图像梯度幅值的稀疏性，在稀疏视图图像重建中取得良好的效果。然而，传统的TV算法的梯度算子具有各向同性和梯度的单一方向的特点，容易导致稀疏视图CT重构图像中边缘和细节信息的丢失。为了从有限的投影数据中重建高质量的图像，利用CS理论开发了各种图像重建方法，例如，基于传统的TV正则化方法导致重构图像中边缘过于平滑的问题，提出了应用相邻图像体素间的各向异性边缘特性的自适应加权TV(AwTV)最小化算法，该AwTV算法可以获得具有明显增益的图像。基于TV最小化的TVM-POCS，提出了结合总变差和L1范数的稀疏投影图像重建，该方法能有效地提高图像的重建质量，并能消除因最小化重构图像的TV范数而产生的块状或斑片状伪影。

中国专利“CN104240210A一种基于压缩感知的CT图像迭代重建方法”的技术方案要点是：首先对投影数据进行维纳滤波降噪处理和小波稀疏变换，将经过初始化的投影数据进行滤波反投影FBP重建获得先验图像，基于图像总变差TV最小化对初始化后的投影数据和重建先验图像进行代数迭代ART重建，反复迭代直至满足收敛条件。该方法能够减少传统滤波反投影FBP及代数迭代ART算法在欠采样投影数据下重建图像的锯齿伪影，同时能够在重建图像中恢复所有重要形态特征，确保重建图像质量。但传统的ART算法每次迭代只用到一条射线投影，测量噪声很容易被引入，并且需要较大的迭代次数，重建效率不高。

中国专利“CN116977470A一种基于压缩感知方向全变分的CT重建方法”的技术方案要点是：采集CT图像的投影数据，利用联合代数重建算法SART更新联合迭代图像，利用梯度下降法TV更新方向迭代修正图像，判断非负约束后的重建图像与迭代初始图像的差值，反复迭代直至满足终止条件。该方法通过将联合代数重建算法SART与方向全变分TV方法结合，解决了有限角度CT重建方法中局部边缘信息的调整问题，在对全变分TV的求解过程中引入了四个方向的信息，能更好的描述图像方向性信息的特征。但是对于高度稀疏的投影数据，使用传统的迭代重建算法重构会发生初始数据与最终重构结果偏差较大，更多的迭代次数会引入噪声导致重建结果质量下降。

中国专利“CN112381904A一种基于DTw-SART-TV迭代过程的有限角度CT图像重建方法”的技术方案要点是：将SART用作每个内部迭代的基本值，并使用TV梯度下降对其进行约束，基于IRS的DTw-SART-TV使用前两个迭代重建结果的线性组合调整下一次迭代的值，将TV的有效降噪性能和两步迭代加速的有效组合来加快CT图像重建的收敛速度，采用动态松弛因子进一步提高了整体收敛速度。但是传统的TV算法的梯度算子具有各向同性和梯度的单一方向的特点，容易导致稀疏视图CT重构图像中边缘和细节信息的丢失。

本发明的目的是在现有图像重建算法基础上提出新的改进方法，基于解析重建算法FBP重构的特点改进重构结果，并将该重构结果作为联合代数重建算法SART的迭代初始条件，再通过改进传统单一方向TV梯度算子对重构切片进行校正，然后循环迭代改进的FBP和SART-TV过程以提高重构图像的质量。

发明内容

为解决上述问题，本发明的目的是提供一种基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法，本发明的技术方案为：

步骤1：使用x射线探测器进行投影数据采集，获得不完全投影数据p；

步骤2：通过FBP算法对不完全投影数据重建，获得预处理的重建图像f；

步骤3：重建图像f作为CFBP过程的初始值，利用重建图像再投影数据与不完全投影数据的差值，以及Ram-Lak滤波器在重建过程中通常会放大投影数据中的高频噪声的特点校正差值投影数据；对校正后的差值投影数据执行滤波反投影重建算法，重建边缘信息来补充先验图像f；

步骤4：对CFBP过程校正的先验信息做矩阵变换后，将其作为SART迭代初始条件重构图像，使SART迭代过程中求解的不适定方程集的适定性得到增强；

步骤5：引入八方向二分均值梯度的FTV项梯度算子执行SART-FTV迭代重构，更好地保留缺失方向上投影数据的边缘和细节；

步骤6：改进的CFBP过程与SART-FTV迭代过程相互优化迭代，使重建的切片逐步逼近真实的图像。

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤1具体为：

步骤1.1：使用x射线探测器在全扫描范围内等间隔采集CT图像的投影数据，组成不完全投影数据集。

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤2具体为：

步骤2.1：通过FBP算法对不完全投影数据重建，滤波反投影算法是把投影数据先做滤波处理，得到修正的投影数据，然后再对修正后的投影数据进行反投影累加处理，获得预处理的先验图像f(x,y)，重建图像在(x,y)处的解析表达式为：

式中g(t,θ)表示经过p(t,θ)滤波修正后的投影数据，h(t)是传递函数H(w)＝|w|的傅里叶逆变换，先验图像f(x,y)为(x,y)点的所有滤波投影的累加。

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤3具体为：

步骤3.1：对不完全投影数据p(t,θ)进行滤波反投影重建，得到重构图像f(x,y)；

步骤3.2：对重构图像f(x,y)执行正向投影求radon变换，得到预估投影数据p

步骤3.3：计算预估投影数据和初始不完全投影数据的差值，得到差值投影数据p

步骤3.4：对差值投影数据执行Ram-Lak滤波器对投影数据滤波处理，然后校正差值投影数据p

步骤3.5：对校正后的差值投影数据执行滤波反投影得到边缘图像f

步骤3.6：用求得的边缘图像f

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤4具体为：

步骤4.1：初始化待重构图像

步骤4.2：计算W的行和向量，估计某一投影角度下通过像素的所有射线的投影值

步骤4.3：计算某一投影角度下实际投影和估计投影的误差

步骤4.4：计算W的列和向量，对第j个像素点的修正值

步骤4.5：对第j个像素点的值进行修正

步骤4.6：重复步骤4.2～4.5，直到达到收敛条件的要求或者指定的迭代次数。

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤5具体为：

步骤5.1：SART过程的初始化

步骤5.2：借助SART算法完成投影数据迭代

步骤5.3：非负约束

步骤5.4：FTV过程的初始化

步骤5.5：计算增量因子

增量因子提供了数据迭代前的图像估计与非负迭代约束执行后的估计之前的差值度量。

步骤5.6：计算全变分梯度以及梯度方向

其中，n为全变分最小化过程的迭代序号，n＝1,2,…,N，N为FTV过程最大迭代次数。

步骤5.7：沿全变分梯度下降的方向迭代修正图像

其中，a为调节因子。反复计算全变分梯度并迭代修正图像，直到n＝N为止。梯度下降过程通过指定该调节因子来控制增量因子沿其方向递增N所执行的梯度下降步骤的总数。

步骤5.8：初始化下一轮迭代

SART-FTV迭代过程依赖SART过程和梯度下降之前的平衡，只要梯度下降引起的图像总变化不超过SART迭代过程引起的图像变化，整个迭代步骤将使图像估计更接近成像线性系统的解空间。

在本发明的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法中，所述步骤6具体为：

步骤6.1：重复步骤3～6，直到达到收敛条件的要求或者指定的迭代次数。

与现有技术相比，本发明提出的技术方案具有以下有益效果：

1、现有的解析重建算法FBP在重建过程中会产生一个低信噪比的重构图，对该重构图正向投影生成低信噪比投影数据，结合Ram-Lak滤波器会放大投影数据中高频噪声的特点，提出改进的FBP算法优化重构结果。

2、将改进的FBP方法优化的重建结果作为SART迭代初始条件应用到SART算法迭代重建的过程中，解决传统的迭代重建算法由于投影数据不完全导致最终重构结果偏差较大的问题，使SART迭代过程中求解的不适定方程集的适定性得到增强。

3、提出一种新的TV项梯度算子，解决传统TV算法由于梯度单一方向导致稀疏视图CT重构图像中边缘和细节信息丢失的问题，加强正则化迭代更新过程使图像估计更接近成像线性系统的解空间。

4、提出一种新的不完全投影CT重建方法，将改进的FBP与改进的SART-TV过程相互优化迭代，使重建图像收敛于原始图像。

附图说明

图1CFBP-SART-FTV算法流程图；

图2稀疏角度投影数据正弦图；

图3稀疏角度投影数据的滤波反投影FBP算法重构结果；

图4FBP算法和改进的CFBP算法重构对比结果；

图5SART算法和改进的CFBP-SART算法重构对比结果；

图6CFBP-SART算法和单次迭代的CFBP-SART-FTV算法重构对比结果；

图7单次迭代的CFBP-SART-FTV算法和CFBP-SART-FTV算法重构对比结果；

图8不完全投影数据基于Shepp-Logan模型的重建算法对比结果；

图9不完全投影数据重构图像和Shepp-Logan模型图像的绝对差值结果；

图10不完全投影数据基于Shepp-Logan模型的重建算法质量评价。

具体实施方式

如图1所示，本发明的一种基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法，包括：

步骤1：使用x射线探测器进行投影数据采集，获得不完全投影数据p；

步骤1.1：使用x射线探测器在全扫描范围内等间隔采集CT图像的投影数据，组成不完全投影数据集。在稀疏视图重构中，投影数据量是影响重建结果的关键因素。本实验在0-180°范围进行等间隔采样，最终使用的稀疏投影分别为30°、60°和90°，采样间隔分别为6°、3°和2°，采样结果如图2所示。

步骤2：通过FBP算法对不完全投影数据重建，获得预处理的重建图像f；

步骤2.1：通过FBP算法对不完全投影数据重建，滤波反投影算法是把投影数据先做滤波处理，得到修正的投影数据，然后再对修正后的投影数据进行反投影累加处理，获得预处理的先验图像。图3为滤波反投影FBP算法对30°、60°和90°稀疏角度的投影数据的重建结果。

步骤3：重建图像f作为CFBP过程的初始值，利用重建图像再投影数据与不完全投影数据的差值，以及Ram-Lak滤波器在重建过程中通常会放大投影数据中的高频噪声的特点校正差值投影数据；对校正后的差值投影数据执行滤波反投影重建算法，重建边缘信息来补充先验图像f；

步骤3.1：对不完全投影数据p(t,θ)进行滤波反投影重建，得到重构图像f(x,y)；

步骤3.2：对重构图像f(x,y)执行正向投影求radon变换，得到预估投影数据p

步骤3.3：计算预估投影数据和初始不完全投影数据的差值，得到差值投影数据p

步骤3.4：对差值投影数据执行Ram-Lak滤波器对投影数据滤波处理，然后校正差值投影数据p

步骤3.5：对校正后的差值投影数据执行滤波反投影得到边缘图像f

步骤3.6：用求得的边缘图像f

FBP算法和改进的CFBP算法重构对比结果如图4所示。在图4中，比较了FBP和改进的CFBP两种方法对不同稀疏角度的重建结果，其中从左到右的列分别为稀疏角数为30°、60°和90°稀疏角度的重构结果，第一行是FBP解析重构算法的结果，第二行是改进的CFBP方法的结果，第三行是原始图像与FBP重构图像的绝对差值，第四行是原始图像与CFBP方法重构图像的绝对差值。

步骤4：对CFBP过程校正的先验信息做矩阵变换后，将其作为SART迭代初始条件重构图像，使SART迭代过程中求解的不适定方程集的适定性得到增强；

步骤4.1：初始化待重构图像

步骤4.2：计算W的行和向量，估计某一投影角度下通过像素的所有射线的投影值

步骤4.3：计算某一投影角度下实际投影和估计投影的误差

步骤4.4：计算W的列和向量，对第j个像素点的修正值

步骤4.5：对第j个像素点的值进行修正

步骤4.6：重复步骤4.2～4.5，直到达到收敛条件的要求或者指定的迭代次数。

SART算法和改进的CFBP-SART算法重构对比结果如图5所示。在图5中，比较了SART算法和改进了迭代初始条件的CFBP-SART这两种方法对不同稀疏角度的重建结果，其中从左到右的列分别为稀疏角数为30°、60°和90°稀疏角度的重构结果，第一行是SART算法的结果，第二行是改进的CFBP-SART方法的结果，第三行是原始图像与SART算法重构图像的绝对差值，第四行是原始图像与改进的CFBP-SART方法重构图像的绝对差值。

步骤5：引入八方向二分均值梯度的FTV项梯度算子执行SART-FTV迭代重构，更好地保留缺失方向上投影数据的边缘和细节；

步骤5.1：SART过程的初始化

步骤5.2：借助SART算法完成投影数据迭代

步骤5.3：非负约束

步骤5.4：FTV过程的初始化

步骤5.5：计算增量因子

增量因子提供了数据迭代前的图像估计与非负迭代约束执行后的估计之前的差值度量。

步骤5.6：计算全变分梯度以及梯度方向

其中，n为全变分最小化过程的迭代序号，n＝1,2,…,N，N为FTV过程最大迭代次数。

步骤5.7：沿全变分梯度下降的方向迭代修正图像

其中，a为调节因子。反复计算全变分梯度并迭代修正图像，直到n＝N为止。梯度下降过程通过指定该调节因子来控制增量因子沿其方向递增N所执行的梯度下降步骤的总数。

步骤5.8：初始化下一轮迭代

SART-FTV迭代过程依赖SART过程和梯度下降之前的平衡，只要梯度下降引起的图像总变化不超过SART迭代过程引起的图像变化，整个迭代步骤将使图像估计更接近成像线性系统的解空间。

CFBP-SART算法和单次迭代的CFBP-SART-FTV算法重构对比结果如图6所示。在图6中，比较CFBP-SART算法和单次迭代的CFBP-SART-FTV两种方法对不同稀疏角度的重建结果，其中从左到右的列分别为稀疏角数为30°、60°和90°稀疏角度的重构结果，第一行是CFBP-SART方法的结果，第二行是单次迭代的CFBP-SART-FTV算法的结果，第三行是原始图像与CFBP-SART方法重构图像的绝对差值，第四行是原始图像与单次迭代的CFBP-SART-FTV算法重构图像的绝对差值。

步骤6：改进的CFBP过程与SART-FTV迭代过程相互优化迭代，使重建的切片逐步逼近真实的图像。

步骤6.1：重复步骤3～6，直到达到收敛条件的要求或者指定的迭代次数。

单次迭代的CFBP-SART-FTV算法和CFBP-SART-FTV算法重构对比结果如图7所示。在图7中，比较单次迭代的CFBP-SART-FTV方法和提出的CFBP-SART-FTV方法对不同稀疏角度的重建结果，其中从左到右的列分别为稀疏角数为30°、60°和90°稀疏角度的重构结果，第一行是单次迭代的CFBP-SART-FTV方法的结果，第二行是提出的CFBP-SART-FTV方法的结果，第三行是原始图像与单次迭代的CFBP-SART-FTV方法重构图像的绝对差值，第四行是原始图像与提出的CFBP-SART-FTV方法重构图像的绝对差值。

为了证明本发明提出的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法能够在稀疏投影数据下重建高质量图像，使用不同的图像重建算法如FBP算法、TV-iFBP算法、TVM-POCS算法以及SART-TV算法的重构结果作为参照进行对比实验，来评估新算法的有效性。

实验应当遵循以下规则，使用相同的稀疏投影数据对照组，重构过程中的参数设置一致。图像的质量评价使用主观评价和客观评价，本发明采用的主观评价是观察者对不同重建方法的重建结果和参考CT图像进行比较，从重建图像是否含有噪声和伪影，图像的边缘结构是否保留完整，有无过度平滑等做出主观判断，本发明采用均方误差(MSE)、峰值信噪比(PSNR)和结构相似性(SSIM)作为客观评价指标对重建结果进行定量评估。

主观评价：

本发明通过对各重建算法的重建结果以及各重建结果与原始图像的绝对差值进行主观评价。基于稀疏投影数据的FBP算法、TV-iFBP算法、TVM-POCS算法、SART-TV算法以及本发明提出的CFBP-SART-FTV重建算法的重构结果如图8所示，图9显示了原始图像与图8所示各算法重构图像的绝对差值。

由实验结果表明，当有足够的稀疏视图时，FBP算法、TV-iFBP算法、TVM-POCS算法、SART-TV算法和CFBP-SART-FTV算法都可以获得高质量的图像。当投影数据高度稀疏时，重建问题变不适定，FBP算法重构图像中出现一些角度受限的伪影，TV-iFBP算法重建结果中会观察到条带伪影以及结构信息丢失。TVM-POCS算法、SART-TV算法重构结果在感兴趣区域处细节不够清晰。相比之下，我们的CFBP-SART-FTV算法具有更好的重建质量，并能有效地抑制稀疏视角带来的噪声。

客观评价：

本发明采用均方误差(MSE)、峰值信噪比(PSNR)和结构相似性(SSIM)作为客观评价指标，对重建结果进行定量评估和对比。

评价指标一：均方误差(MSE)主要用于衡量重建图像和参考图像间各像素的差异，重建图像与参考图像像素值误差的平方和与总个数的比值越小，重建图像与参考图像也就越接近。均方误差(MSE)定义为

其中，X是重建图像，X是参考图像，x

评价指标二：峰值信噪比(PSNR)作为评估图像的常用客观指标，单位为dB，PSNR值越大，表明图像失真越少，重建图像与参考图像也就越接近，峰值信噪比(PSNR)定义

其中MAX为图片可能的最大灰度值。

评价指标三：结构相似性(SSIM)综合考虑了重建图像和参考图像的亮度、对比度、结构的相似性，是通过图像结构信息的变化来衡量图像相似性的，它的取值范围为[0,1]，SSIM的值越接近1，表示重建图像与标准图像结构越相似。结构相似性(SSIM)定义为

其中，μ

基于稀疏投影数据的FBP算法、TV-iFBP算法、TVM-POCS算法、SART-TV算法以及本发明提出的CFBP-SART-FTV迭代重建方法的重构结果质量的定量指标比较如图10所示。

由图中数据结果可知：①在各稀疏投影角度下，CFBP-SART-FTV算法始终优于其它四种传统重建算法，即使只有30个投影(PSNR＝75.4426,MSE＝0.001706,SSIM＝0.99995)的CFBP-SART-FTV算法重构的图像也优于90个投影的(PSNR＝74.0253,MSE＝0.002577,SSIM＝0.99995)FBP算法、(PSNR＝73.9285,MSE＝0.002632,SSIM＝0.99993)TV-iFBP算法、(PSNR＝71.5399,MSE＝0.004561,SSIM＝0.99983)TVM-POCS算法、(PSNR＝71.6335,MSE＝0.004464,SSIM＝0.99984)SART-TV算法，因此新提出的CFBP-SART-FTV算法的性能是最优的。②当稀疏投影角度在30°到90°时，CFBP-SART-FTV算法对应的结构相似性(SSIM)的值总是大于0.99990，且幅度变换很小，因此CFBP-SART-FTV算法的重构质量在这个投影数范围内是稳定的。③通过比较30个投影(PSNR＝75.4426,MSE＝0.001706,SSIM＝0.99995)、60个投影(PSNR＝75.1536,MSE＝0.001985,SSIM＝0.99992)以及90个投影(PSNR＝74.2557,MSE＝0.002441,SSIM＝0.99990)的重建结果表明，CFBP-SART-FTV算法在高度稀疏角度下的重构质量更好。

综合上述对比实验的数据表现，本发明所设计的基于CFBP-SART-FTV迭代重建的不完全投影CT重建方法的重建质量都要优于其他重建方法的质量，本发明有效解决了传统FBP算法由于投影缺失角范围的增大而出现的角度伪影，解决了SART算法由于缺乏约束方程而使重构解位于过大的解域中的随机点上，解决了传统的具有单一梯度的TV算法导致重构图像中边缘和细节信息丢失的问题，证明了该方法在有限角度应用中的可行性，且具有较大的现实意义与应用价值。