CN106846427B - 一种基于重加权各向异性全变分的有限角度ct重建方法 - Google Patents

Abstract

本发明公开了一种基于重加权各向异性全变分的有限角度CT重建方法，包括根据CT成像原理建立CT图像重建方程式，初始化各项参数，迭代求解数据保真项和重加权各向异性全变分极小化项，直至满足迭代停止条件。本发明通过将各向异性的思想融入重加权的方法，有效的提高了有限角度CT重建算法的重建恢复效果，更好的还原了真实的CT图像。与现有重建方法的实验比较表明，本发明方法不仅能够得到更好的重建结果，而且重建所需的时间也更少。

Description

一种基于重加权各向异性全变分的有限角度CT重建方法

技术领域

本发明属于CT成像技术领域，具体涉及一种基于重加权各向异性全变分的有限角度CT重建方法。

背景技术

X光CT在很多应用领域中都有很重要的应用，如用于诊断和治疗的医学成像领域、安全检查领域以及产品质量检测控制领域。除了精确和鲁棒的图像重建，有限角度问题也受到了广泛的关注。一方面为了病人的健康考虑，希望能够尽量减少病人受到的X光照射，另一方面病人本身很难长时间的保持不动配合X光照射，所以希望能够尽可能的缩短扫描的时间。这些一方面可以通过降低X光剂量来实现，另一方面可以通过减少X光照射时间来实现，包括以更加稀疏的视角或者限制投影视角的大小来获取投影数据。还有一些实际场合也会有扫描角度的限制，比如X光显微镜对生物样品进行成像时会受到样品固定器的限制而无法进行全角度的成像。有限角度CT重建是一个严重的病态问题，因为投影数据的角度范围小于精确重建的理论要求。为了解决这个问题，一些先验信息常被用来作为问题的约束条件，如图像的非负性，轮廓或边界，图像的稀疏性等等。

传统的解析算法如滤波反投影算法(filtered back-projection，FBP)在投影数据角度缺失时无法得到较为精确的重建效果，相比而言迭代图像重建算法能得到稍好一些的结果。其主要分为两类，一类是统计图像重建算法(statistical imagereconstruction，SIR)，其利用基本的物理特性(如光子服从泊松分布的特性)建立模型，然后通过最大化重建图像和测量数据之间的匹配程度来估计最佳重建图像。另一类是代数重建算法(algebraic reconstruction technique，ART)，其假定目标横截面由一系列未知的像素/体素组成，然后建立代数方程，利用测量的投影数据和系统矩阵来进行重建。在CT重建领域较为常用的是ART算法。对于这两种算法，通常都还需要利用一些先验信息作为约束条件。

近年来压缩感知(compressed sensing，CS)理论的发展大大推动了欠采样信号重建算法的研究。该理论说明，在一定的条件下，即使是采样不满足奈奎斯特采样定理的信号，也可以被精确的重建和恢复。压缩感知的关键就在于寻找一个合适的转换空间来对信号进行稀疏的表达。在图像重建领域经常使用的转换空间是离散梯度空间和小波空间。对于大多数自然图像而言，尤其是医学图像，快速变化的区域只存在于一些结构的边界处，很多区域都是局部平滑的。所以即使一幅图像其自身不是稀疏的，它的梯度图像却很有可能是稀疏的。压缩感知的一个特例是全变分(total variation，TV)极小化，其利用了图像的这种稀疏的性质，常被应用于稀疏视角CT图像重建领域。一幅图像的TV就是其梯度图像的l₁范数，对其进行极小化可以作为由CT投影得到的数据保真项的一个约束条件。这种方法可以得到较好的重建结果，伪影较少，且边界也恢复的较好。

但由于传统的TV方法假定图像是局部平滑的，所以重建结果的边界区域往往会有过平滑的问题。为了解决这个问题，很多人在这方面做了很多工作，提出了很多新的算法，然而效果并不显著。于是很多研究者就开始探究是否可以对l₁正则化进行进一步的优化。实际上，l₀范数才是最直接的稀疏表达，它通过计算图像中非零元素的个数得到，所以它应该是利用稀疏先验的更好的选择。然而，由于求解l₀范数极小化问题是一个非凸问题，所以还没有算法可以对其进行直接求解。数学上来说，l₀范数可以看成是l_p范数在p趋近于0时的极限情况。当0＜p＜1时，l_p范数是非凸的，有很多个局部极值点。为了解决这个问题，一种迭代重加权的方法被提出，有很多人都在这方面做了很多研究，验证了它的有效性。这种迭代重加权方法的基本原理是将非凸的l_p范数极小化问题转换成一系列重加权的TV极小化问题，虽然整个算法是一个非凸的问题，但在每一次迭代过程中，需要解决的是一个凸问题。通过范数的定义可以知道，l₀范数和l₁范数的区别就在于对于信号幅值的利用方式，l₀范数计算非负元素的个数，而l₁范数计算的是所有非负元素的值的绝对值之和。重加权方法的目的就是尽可能的矫正l₀范数和l₁范数之间这一关键性的差异。经常使用的方法是交替的进行图像重建和权重系数更新，直到迭代终止条件被满足。

因为我们所需要研究的是有限角度重建问题，所以除了图像的稀疏性这一先验信息之外，我们还有另一个先验信息可以加以利用，那就是角度范围信息。1988年的时候，Quinto从理论上分析了从有限角度投影数据重建得到的图像的性质特点，给出了一个结论：与投影方向相切的方向的边界和细节信息更容易被恢复，而在另一些特定的角度则会出现一些伪影和模糊。TV方法在进行图像重建时无法探测到某些特定角度的较为模糊的边界。此外，由于TV是一幅图像中所有像素点的梯度幅值之和，是各向同性的，它同时也会包含那些模糊边界的贡献，导致对不模糊的边界的保边能力也会下降。为了更好的解决这个问题，几年前各向异性全变分(anisotropic total variation，ATV)的概念被提出，除了图像的稀疏性之外，它将投影的角度范围作为另一个先验信息。这个方法的主要思想就是尽可能的减少模糊边界的信息对边界探测的影响，从而可以得到更好的图像重建结果。

发明内容

鉴于上述，本发明提供了一种基于重加权各向异性全变分的有限角度CT重建方法，该方法解决了现有CT重建算法存在的过平滑和部分边界模糊的问题，能够提高CT图像重建的质量。

一种基于重加权各向异性全变分的有限角度CT重建方法，包括如下步骤：

(1)采集探测器在不同方向测得的CT图像的投影数据，组成投影数据集p；

(2)根据CT成像原理和投影数据集p，建立CT图像重建方程式：

min_f||f||_RwATV， s.t. Wf＝p，f_x，y≥0(1)

公式(1)中Wf＝p为数据保真项，min_f||f||_RwATV为极小化项；

其中，W是系统矩阵，W中的每一个元素是投影光束与重建图像f中像素点的相交长度，f_x，y为重建图像f在位置(x，y)处的像素值，其大小对应于待重建断层图像不同像素处的X光线性吸收系数，x为每一个像素值所处位置的横向坐标，且x∈1～N_width，y为每一个像素值所处位置的纵向坐标，且y∈1～N_height，N_width和N_height分别表示图像f的宽度和高度；

(3)对公式(1)进行迭代优化求解，得到最终重建图像。

所述公式(1)中，

其中，R为权重矩阵，R中的每一个元素r_x，y是的权重值，为f的各向异性全变分，其计算方式如下：

其中，A和B分别为式(3)中差分项(f_x，y-f_x-1，y)²和(f_x，y-f_x，y-1)²的权值系数，均为正实数；f_x-1，y为重建图像f在位置(x-1，y)处的像素值，f_x，y-1为重建图像f在位置(x，y-1)处的像素值。

所述的投影光束为平行光束、扇形光束或锥形光束。

步骤(2)中，为了减少计算复杂度并节约程序运行时间，作为优选，系统矩阵W的每一个元素是投影光束与重建图像f中像素点的相交长度。

所述步骤(3)的求解过程为：

(3-1)初始化各项参数：总迭代次数N_iter、极小化项求解迭代次数N_TV、极小化项求解权重系数a、参数ε、参数ξ、迭代终止阈值δ；

(3-2)初始令迭代计数k＝1，

(3-3)采用代数重建算法求解公式(1)中的数据保真项，得到重建图像f_k，0，具体计算公式为：

其中，p_i为第i条投影光线的投影数据值，i的取值范围为1～N_view，N_view为投影总数目，fⁱ为对第i个投影数据重建得到的图像，f^i-1为对第i-1个投影数据重建得到的图像，w_i是系统矩阵W的一个子矩阵，w_i中的每一个元素是第i条投影光线与图像fⁱ中的像素点的相交长度，为fⁱ在位置(x，y)处的像素值；

(3-4)采用梯度下降法求解公式(1)中的极小化项，得到重建图像f_k，1；

(3-5)判断重建图像f_k，1与重建图像f_k-1，1的差值图像的像素值是否小于迭代终止阈值δ或k是否大于N_iter，若是，执行步骤(3-6)，若否，令k＝k+1，执行步骤(3-3)～步骤(3-5)；

(3-6)终止迭代，得到最终重建图像f_k，1。

步骤(3-1)中，所述的总迭代次数N_iter的取值范围为3～30000；所述的极小化项求解迭代次数N_TV的取值范围为1～1000；所述的极小化项求解权重系数a的取值范围为0.01～10；所述的参数ε和ξ的取值范围都为10^-10～1。

步骤(3-4)的具体步骤为：

(3-4-1)更新权重矩阵R，更新公式为：

其中，为求解极小化项时第n次迭代的重建图像在位置(x，y)处的像素值，为第n次迭代的重建图像在位置(x-1，y)处的像素值，为第n次迭代的重建图像在位置(x，y-1)处的像素值，为的权重值，n的取值范围为1～N_TV，参数ξ为防止出现奇异值而设定的一个数值较小的实数；

(3-4-2)对重加权各项异性全变分项||f||_RwATV进行求导，具体为：

其中，ε为防止出现奇异值而设定的一个数值较小的实数；u_x，y为位置(x，y)处的导数值，为的权重值，为的权重值；

(3-4-3)利用公式(6)对待重建图像进行更新，得到重建图像f_k,1；

其中，u为u_x,y组成的矩阵。

本发明的CT图像重建方法通过利用有约束的CT图像重建模型，将迭代重加权的方法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的有条状伪影和部分区域模糊的问题，在对该模型的求解过程中交替采用代数重建算法和梯度下降法，直到迭代终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。

附图说明

图1为本发明CT图像重建算法的流程示意图；

图2为Shepp&Logan phantom模型的真值图像；

图3为Shepp&Logan phantom模型的真值图像对应的灰度数值图；

图4为Shepp&Logan phantom模型的真值图像对应的正弦图(即不同角度的投影数据组成的图，角度范围为180°)；

图5为重建过程中两种方法得到的重建图像相比于真值图像的均方根误差(rootmean square error，RMSE)随迭代次数的变化，投影数据的角度范围为120°；

图6为两种方法重建后的CT图像，其中，图6(a)为Shepp&Loganphantom模型采用本发明方法重建后的CT图像，投影数据角度范围为120°，迭代次数为10000次；图6(b)为Shepp&Logan phantom模型采用各向异性全变分方法重建后的CT图像，投影数据角度范围为120°，迭代次数为10000次；

图7为两种方法得到的重建图像的垂直中线的剖线图与真值的对比，迭代次数为10000次。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，本发明基于重加权各向异性全变分的有限角度CT重建方法，包括如下步骤：

S1.根据CT成像原理和采集的投影数据建立CT图像重建方程式：

min_f||f||_RwATV， s.t. Wf＝p，f_x，y≥0(1)

其中，Wf＝p为数据保真项，min_f||f||_RwATV为极小化项，且：

S2.初始化各项参数：总迭代次数N_iter、极小化项求解迭代次数N_TV、极小化项求解权重系数a、参数ε、参数ξ、迭代终止阈值δ；

S3.用代数重建算法来对数据保真项进行求解，得到重建图像f_k，0，具体公式为：

S4.用梯度下降法来对极小化项进行求解，得到重建图像f_k，1；

S5.判断重建图像f_k，1与重建图像f_k-1，1的差值图像的像素值是否小于迭代终止阈值δ或k是否大于N_iter，若否，返回到步骤S3继续执行，若是，迭代停止，且保存相关变量。

实施例1

本实施例中，采用平行光束照射CT图像，投影总数目为256×80，设定总迭代次数N_iter为10000，极小化项求解迭代次数N_TV为20，极小化项求解权重系数a为0.15，参数ε为10^-8，参数ξ为为10^-8，迭代终止阈值δ为10^-4。

通过对Shepp&Logan phantom模型的正弦图进行重建来验证本实施方式的实用性和可靠性。Shepp&Logan phantom模型的真值图像如图2所示，其对应的灰度数值图如图3所示，正弦图(180°全角度)如图4所示。采用本发明方法(ART+RwATV，参数的选择为：A＝1，B＝0.001)和各向异性全变分方法(ART+ATV，参数的选择为：A＝1，B＝0.001)对投影数据(采用的投影数据角度范围为120°)进行重建，图5给出了重建过程中两种方法得到的重建图像相比于真值图像的均方根误差随迭代次数的变化。可以看到两种方法都随着迭代次数的增加而逐渐收敛，本发明方法的重建结果均方根误差更小，其收敛速度也更快。我们选取最大迭代次数(即10000次)时的重建图像加以显示，图6(a)为采用各向异性全变分方法重建后的CT图像，图6(b)为采用本发明方法重建后的CT图像，可看到本发明方法重建的图像伪影较小，能更好的恢复图像的细节信息，重建效果优于各向异性全变分方法。为了更直观的证明本发明方法的优越性，图7给出了图6(a)和6(b)所示的重建图像的剖线图与真值的剖线图的对比，可以看出，本发明方法能得到更加接近于真值的重建结果。对算法的相关参数及重建结果相关信息进行总结，如表1所示，可以看到本发明方法得到的重建结果均方根误差更小，且程序运行时间也小于各向异性全变分算法。

表1

算法	A	B	迭代次数	均方根误差	时间(s)
						ART+ATV	1	0.001	10000	1.470846E-02	1632.676000
ART+RwATV	1	0.001	10000	4.498879E-03	1398.482666

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的一个优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于重加权各向异性全变分的有限角度CT重建方法，包括如下步骤：

(1)采集探测器在不同方向测得的CT图像的投影数据，组成投影数据集p；

(2)根据CT成像原理和投影数据集p，建立CT图像重建方程式：

min_f||f||_RwATV，s.t.Wf＝p，f_x，y≥0 (1)

公式(1)中Wf＝p为数据保真项，min_f||f||_RwATV为极小化项；

其中，W是系统矩阵，W中的每一个元素是投影光束与重建图像f中像素点的相交面积，f_x，y为重建图像f在位置(x，y)处的像素值，其大小对应于待重建断层图像不同像素处的X光线性吸收系数，x为每一个像素值所处位置的横向坐标，且x∈1～N_width，y为每一个像素值所处位置的纵向坐标，且y∈1～N_height，N_width和N_height分别表示图像f的宽度和高度；

所述公式(1)中，

其中，R为权重矩阵，R中的每一个元素r_x，y是的权重值，为f的各向异性全变分，其计算方式如下：

其中，A和B分别为式(3)中差分项(f_x，y-f_x-1，y)²和(f_x，y-f_x，y-1)²的权值系数，均为正实数；f_x-1，y为重建图像f在位置(x-1，y)处的像素值，f_x，y-1为重建图像f在位置(x，y-1)处的像素值；(3)对公式(1)进行迭代优化求解，得到最终重建图像；具体求解过程为：

(3-1)初始化各项参数：总迭代次数N_iter、极小化项求解迭代次数N_TV、极小化项求解权重系数a、参数ε、参数ξ、迭代终止阈值δ；

(3-2)初始令迭代计数k＝1，

(3-3)采用代数重建算法求解公式(1)中的数据保真项，得到重建图像f_k，0，具体计算公式为：

其中，p_i为第i条投影光线的投影数据值，i的取值范围为1～N_view，N_view为投影总数目，fⁱ为对第i个投影数据重建得到的图像，f^i-1为对第i-1个投影数据重建得到的图像，w_i是系统矩阵W的一个子矩阵，w_i中的每一个元素是第i条投影光线与图像fⁱ中的像素点的相交长度；

(3-4)采用梯度下降法求解公式(1)中的极小化项，得到重建图像f_k，1；

(3-5)判断重建图像f_k，1与重建图像f_k-1，1的差值图像的像素值是否小于迭代终止阈值δ或k是否大于N_iter，若是，执行步骤(3-6)，若否，令k＝k+1，执行步骤(3-3)～步骤(3-5)；

(3-6)终止迭代，得到最终重建图像f_k，1。

2.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：所述的投影光束为平行光束、扇形光束或锥形光束。

3.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：步骤(3-1)中，所述的总迭代次数N_iter的取值范围为3～30000。

4.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：步骤(3-1)中，所述的极小化项求解迭代次数N_TV的取值范围为1～1000。

5.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：步骤(3-1)中，所述的极小化项求解权重系数a的取值范围为0.01～10。

6.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：步骤(3-1)中，所述的参数ε和ξ的取值范围都为10^-10～1。

7.根据权利要求1所述基于重加权各向异性全变分的有限角度CT重建方法，其特征在于：步骤(3-4)的具体步骤为：

(3-4-1)更新权重矩阵R，更新公式为：

其中，为求解极小化项时第n次迭代的重建图像在位置(x，y)处的像素值，为第n次迭代的重建图像在位置(x-1，y)处的像素值，为第n次迭代的重建图像在位置(x，y-1)处的像素值，为的权重值，n的取值范围为1～N_TV，参数ξ为防止出现奇异值而设定的一个数值较小的实数；

(3-4-2)对重加权各项异性全变分项||f||_RwATV进行求导，具体为：

其中，ε为防止出现奇异值而设定的一个数值较小的实数；u_x，y为位置(x，y)处的导数值，为的权重值，为的权重值；

(3-4-3)利用公式(6)对待重建图像进行更新，得到重建图像f_k，1；

其中，u为u_x，y组成的矩阵。