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

摘要

本发明公开了一种基于各向异性全变分的有限角度CT重建方法，其通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。

权利要求

1.一种基于各向异性全变分的有限角度CT重建方法，包括如下步骤：(1)利用探测器在不同角度方向采集获取CT图像的投影数据，组成投影数据集 所述投影数据集 由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；

(2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；

低剂量模式下的CT图像重建方程表达式为：

稀疏视角模式下的CT图像重建方程表达式为：其中：为CT图像数据集且维度为k，k为CT图像的像素点总数，中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| ||ATV表示各向异性全变分；

(3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量 和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；

(4)根据实际情况选择对应模式下的目标函数，并采用交替投影接近算法对其进行优化求解以重建得到CT图像。

2.根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：首先，引入变量 将低剂量模式下的CT图像重建方程改写为如下形式：其中：表示由CT图像数据集 计算得到的正向投影数据；

然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：式(1.3)对应的拉格朗日方程 如下：

其中：为拉格朗日乘子向量,T表示转置；

则将式(1.3)转换成鞍点问题，具体表达形式如下：最后，通过极小化式(1.4)将变量 去除，从而得到对应的目标函数如下：其中：·表示内积运算。

3.根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：式(2.2)对应的拉格朗日方程 如下：

则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：其中：为拉格朗日乘子向量,T表示转置，β为给定的权重系数。

4.根据权利要求2所述的有限角度CT重建方法，其特征在于：对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：P＝F-1R(ω)F

其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。

5.根据权利要求3所述的有限角度CT重建方法，其特征在于：对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：P＝F-1R(ω)F

其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。

说明书

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

技术领域

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

背景技术

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

[0003] 对于大多数的自然图像而言，尤其是医学图像，快速变化的区域只存在于一些结构的边界处，图像的大部分区域都是局部平滑的，所以即使一幅图像其本身并不稀疏，它的梯度图像却很有可能是稀疏的。全变分(Total Variation，TV)极小化就利用了图像的这种稀疏性质，常被应用于CT图像重建领域；一幅图像的TV就是其梯度图像的L1范数，对其进行极小化可以作为由CT投影得到的数据保真项的一个约束条件，这种方法可以得到较好的重建结果，伪影较少，且边界也恢复的较好。

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

[0005] 当前已存在很多图像重建算法，通常是基于对一些经典的迭代算法如梯度下降法(Gradient Descent，GD)、代数重建算法(Algebraic Reconstruction Technique，ART)、同步迭代重建技术(Simultaneous Iterative Reconstruction Technique，SIRT)、同步代数重建算法(Simultaneous Algebraic Reconstruction Technique，SART)等进行改进得到的，常常还会结合有序子集法(Ordered-Subsets Technique，OS)来进行加速。这些算法中的大部分主要存在两个问题：一是在从投影数据重建图像时精度不高，二是收敛速度慢导致重建所需要的时间很长。Kudo等人在2016年提出了一种快速迭代重建算法，用来求解低剂量CT和稀疏视角CT重建问题，该方法通过极小化数据保真项来进行迭代重建，并用TV惩罚项来对此过程进行约束，其求解思路如下：首先，通过使用拉格朗日对偶方法，将最初的极小化问题转换成鞍点(原始-对偶)问题，然后使用一阶原始-对偶方法中的交替投影接近法(Alternating Projection Proximal，APP)来对其进行求解；由于这种方法的整体结构和同时迭代重建算法如SIRT和SART算法相同，收敛速度很慢，于是引入了滤波反投影(Filtered Back-Projection，FBP)重建算法中的斜坡滤波器来对迭代方程进行预处理，这种处理方式并不会影响最终的迭代结果，其值与原始问题的解完全相同，最终得到的算法可以理解为通过对斜坡滤波器的使用来达到利用FBP型预处理以加速一阶原始-对偶方法的目的。

发明内容

[0006] 鉴于上述，本发明提供了一种基于各向异性全变分的有限角度CT重建方法，该方法除了能够解决现有CT重建算法存在的部分边界模糊的问题，提高CT图像重建的质量之外，还能够以很快的收敛速度精确求解目标方程。

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

[0008] (1)利用探测器在不同角度方向采集获取CT图像的投影数据，组成投影数据集所述投影数据集 由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；

[0009] (2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；

[0010] 低剂量模式下的CT图像重建方程表达式为：

[0011]

[0012] 稀疏视角模式下的CT图像重建方程表达式为：

[0013]

[0014] 其中：为CT图像数据集且维度为k，k为CT图像的像素点总数， 中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| ||ATV表示各向异性全变分；

[0015] (3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量 和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；

[0016] (4)根据实际情况选择对应模式下的目标函数，并对其优化求解以重建得到CT图像。

[0017] 进一步地，所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：

[0018] 首先，引入变量 将低剂量模式下的CT图像重建方程改写为如下形式：

[0019]

[0020] 其中：表示由CT图像数据集 计算得到的正向投影数据；

[0021] 然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：

[0022]

[0023] 式(1.3)对应的拉格朗日方程 如下：

[0024]

[0025] 其中： 为拉格朗日乘子向量,T表示转置；

[0026] 则将式(1.3)转换成鞍点问题，具体表达形式如下：

[0027]

[0028] 最后，通过极小化式(1.4)将变量 去除，从而得到对应的目标函数如下：

[0029]

[0030] 其中：·表示内积运算。

[0031] 进一步地，所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：

[0032] 首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：

[0033]

[0034] 式(2.2)对应的拉格朗日方程 如下：

[0035]

[0036] 则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：

[0037]

[0038] 其中：为拉格朗日乘子向量,T表示转置，β为给定的权重系数。

[0039] 对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：

[0040] P＝F-1R(ω)F

[0041]

[0042] 对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：

[0043] P＝F-1R(ω)F

[0044]

[0045] 其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。

[0046] 进一步地，所述步骤(4)中采用交替投影接近算法对目标函数进行优化求解以重建得到CT图像，该算法中涉及到的各项参数包括总迭代次数Niter、TV去噪迭代次数NATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ、ATV权重因子γh和γv；其中，总迭代次数Niter的取值范围为1～50000，TV去噪迭代次数NATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围为10-10～1，TV项权重系数β的取值范围为10-6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10-3～103，参数σ的取值范围为10-3～103，ATV权重因子γh和γv的取值范围为10-6～1。

[0047] 本发明CT图像重建方法通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。

附图说明

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

[0049] 图2为Shepp&Logan phantom模型的真值图像，显示灰度范围为[1.00,1.05]。

[0050] 图3为Shepp&Loganphantom模型真值图像对应的灰度数值图。

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

[0052] 图5为重建过程中两种方法(本发明方法和各向异性全变分方法)得到的重建图像相比于真值图像的均方根误差(rootmean square error，RMSE)随迭代次数的变化示意图，投影数据的角度范围为120°，角度数为120。

[0053] 图6(a)为Shepp&Loganphantom模型采用本发明方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。

[0054] 图6(b)为Shepp&Loganphantom模型采用各向异性全变分方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。

[0055] 图7为两种方法(本发明方法和各向异性全变分方法)得到的重建图像的水平中线的剖线与真值剖线的对比示意图，迭代次数为3000次。

具体实施方式

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

[0057] 如图1所示，本发明基于各向异性全变分的有限角度CT快速迭代重建算法，包括如下步骤：

[0058] (1)采集探测器在不同方向测得的CT图像的投影数据，组成投影数据集 投影数据集 由对应各个角度方向下所采集得到的投影数据向量组成，投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数。

[0059] (2)分别建立低剂量和稀疏视角CT的图像重建方程式：

[0060] 低剂量CT：

[0061] 稀疏视角CT：其中：为CT图像数据集且维度为k，k为CT图像的像素点总数，中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A是系统矩阵，它适用于平行光束投影、扇形光束投影及锥形光束投影等各种形式；β为权重系数； 为各向异性全变分项。

[0062] 在低剂量CT情况下，我们假定nm＞k并且 含有噪声；在稀疏视角CT情况下，我们假定nm＜k并且 不含噪声，m为角度方向数量。

[0063] (3)对公式(1)和(2)分别进行预处理：

[0064] 对于低剂量CT，引入变量 则式(1)可改写成如下形式：

[0065]

[0066] 其中： 可以看作是从图像 计算得到的正向投影。引入非奇异矩阵P1/2，则上式可进一步改写为：

[0067]

[0068] 上式对应的拉格朗日方程可写成如下形式：

[0069]

[0070] 其中：是拉格朗日乘子向量，则式(3)式可转换成鞍点问题：

[0071]

[0072] 因为 可以通过极小化上式去除，于是我们可以得到：

[0073]

[0074] P＝F-1R(ω)F

[0075]

[0076] 其中：·表示内积运算，F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量。

[0077] 对于稀疏视角CT，同样引入非奇异矩阵P1/2，则式(2)可改写成如下形式：

[0078]

[0079] 同样的可写出上式对应的拉格朗日方程：

[0080]

[0081] 则式(5)可转换成如下形式：

[0082]

[0083] P＝F-1R(ω)F

[0084]

[0085] (4)使用交替投影接近算法APP来对式(4)和(6)进行求解：

[0086] 4-1初始化 并设置各项参数的值：总迭代次数Niter、TV去噪迭代次数NATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ，ATV权重因子γh和γv；其中，总迭代次数Niter的取值范围为1～50000，TV去噪迭代次数NATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围均为10-10～1，TV项权重系数β取值范围为10-6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10-3～103，参数σ的取值范围为10-3～103，ATV权重因子γh和γv取值范围均为10-6～1。

[0087] 4-2初始令迭代计数k＝1；

[0088] 4-3对 进行更新，当k为1时：

[0089]

[0090] 当k不为1时：

[0091] 低剂量CT：

[0092] 稀疏视角CT：

[0093] 4-4对 进行更新：

[0094]

[0095] 可用Chambolle的TV去噪算法来对上式进行求解，由于 这一项在这一步中是常量，为方便表述引入 则式(7)的求解方程为：

[0096]

[0097]

[0098] 其中： 为梯度算子，假定重建图像的大小为N×N像素，在本实施方式中，其定义如下：

[0099]

[0100]

[0101]

[0102] 则各项异性全变分项 的计算方式如下：

[0103]

[0104] 对每一个 div算子的计算方式如下：

[0105]

[0106] 计算 若∈c＞∈且重复次数未达到NATV，重复计算式(8)～式(9)。

[0107] 4-5再次对 进行更新：

[0108] 低剂量CT：

[0109] 稀疏视角CT：

[0110] 4-6判断重建图像 与重建图像 的差值图像的所有像素值的平方和是否小于迭代终止阈值δ或k是否大于Niter：若是，执行步骤4-7；若否，令k＝k+1，执行步骤4-3～步骤4-5。

[0111] 4-7终止迭代，得到最终重建图像

[0112] 以下我们通过对Shepp&Loganphantom模型的正弦图进行重建来验证实施方法的实用性和可靠性。Shepp&Loganphantom模型的真值图像如图2所示(由于其大部分结构的灰度值非常接近，为了清楚的显示其结构，这里选择灰度范围[1.00,1.05]进行显示)，其对应的灰度数值图如图3所示，正弦图(角度范围为180°，角度数为180)如图4所示。采用本发明方法(proposed algorithm)和各向异性全变分方法(ART+ATV)进行重建，图5给出了重建过程中两种方法得到的重建图像相比于真值图的均方根误差随迭代次数的变化，可以看到两种方法都随着迭代次数的增加而逐渐收敛，本发明方法的重建结果均方根误差更小，且其收敛速度非常的快。我们选取最大迭代次数(即3000次)时的重建图像加以显示，图6(a)为采用各向异性全变分方法重建后的CT图像，图6(b)为采用本发明方法重建后的CT图像，可看到本发明方法重建的图像伪影要小得多，能更好的恢复图像的细节信息，重建效果显著优于各向异性全变分方法。为了更直观的证明本发明方法的优越性，图7给出了图6(a)和图6(b)所示的重建图像的水平中线的剖线与真值剖线的对比，可以看到，本发明方法的剖线几乎和真值完全重合，即其能够得到更加接近于真值的重建结果。

[0113] 上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。