CN109447913B - 一种应用于不完备数据成像的快速图像重建方法 - Google Patents

Abstract

本发明公开一种应用于不完备数据成像的快速图像重建方法，基于成像系统离散数学模型，利用图像全变分(TV)范数作为先验信息，建立成像问题的最优化目标函数；结合自适应代数重建技术、邻近映射技术及单调外插加速技术迭代求解最优化问题，利用自适应代数重建技术提高数据保真度，利用邻近映射技术降低对应图像的TV范数，实现稀疏或不完备数据条件下重建高质量图像。

Description

一种应用于不完备数据成像的快速图像重建方法

技术领域

本发明属于图像重建技术领域，特别是涉及一种应用于不完备数据成像的快速图像重建方法。

背景技术

图像包含的信息量十分丰富，是人类获取信息的主要来源。成像系统已在生物医学工程、航空航天、无损检测、公共安全等众多领域获得了广泛的应用。比如，计算机断层扫描成像(CT)作为一种重要的无损检测技术，以其无损、精确、三维可视化等优点广泛应用于工业、医疗、安检等领域。成像系统由硬件和软件两部分构成，成像算法是软件部分的核心，在很大程度上决定了成像系统的性能，是目前成像技术研究的热点。

对于许多线性成像系统，可采用的成像算法有解析法与迭代求解方法两个大类。一般而言，解析法具有计算量小、重建速度快的优点，当投影数据完备与理想时，可获得良好的重建图像。但在实际检测中，由于各种原因(如成本、操作条件、人体健康等)，经常出现获取的测量数据不完备的情况。当测量数据不完备或稀疏时，解析法通常难以得到较理想的重建图像，无法满足实际应用需求。

与解析法相比，迭代求解方法的应用范围更广，适用于不同方式的采样数据，并可灵活结合先验知识进行重建，实现在稀疏或不完备数据条件下重建高质量图像。但迭代求解法也存在计算量大、重建速度慢的缺点，这些缺点在很长一段时间内限制了其在成像领域的推广与应用。近年来随着计算机与硬件技术的快速发展，迭代求解方法重新受到人们的青睐。因此，结合线性成像系统数学模型与先验信息，研究开发简单高效的快速图像重建算法，对实现不完备投影数据情况下快速精确成像，拓展成像系统应用领域具有重要的意义。

代数重建法(ART)是一种广泛应用于图像重建(特别是CT成像)的迭代求解方法，它将图像重建问题转化为求解线性方程组，当投影数据不完备时(对应欠定方程组)，仍可结合先验信息进行求解。在噪声条件下，ART中的松弛参数对算法性能有重要的影响。总体上，在大噪声条件下应使用较小的松弛参数，在小噪声条件下应使用较大的松弛参数。但实际中松弛参数的选择取决于诸多因素，因此若能在算法收敛过程中自适应地调整松弛参数的值对于提高算法收敛速度与成像质量将有良好效果。在不完备数据条件下，由于对应的方程组是欠定的，仅采用代数重建法求解线性方程组无法获得满意的成像结果，还需要结合先验信息进行优化。由于实际图像的全变分(total variation,TV)范数一般都比较小，因此常被作为先验信息用于图像重建。在实际应用中，图像的TV范数常作为待优化问题的正则项，本发明方法亦将采用TV范数正则化技术。这样，待重建的图像需要满足两类约束：(1)数据保真度约束(重建的图像要满足测量数据的约束，这里主要指线性方程组约束)；(2)先验信息约束(重建的图像要与先验信息吻合，这里主要指图像的TV范数应该比较小)。在设计基于迭代求解的快速图像重建算法时，就要兼顾这两方面的约束。同时为了降低计算复杂度，需要探求计算量低且能快速收敛的高效算法。

发明内容

为了解决上述问题，本发明提供一种应用于不完备数据成像的快速图像重建方法，该方法结合自适应代数重建技术与邻近算子(proximal operator)求解TV范数约束的图像重建问题。在每次迭代，包含两类操作：第一类操作执行自适应代数重建法的算子来提高数据保真度，松弛参数根据当前估计向量到对应超平面的相对距离自适应调整；第二类操作采用邻近算子在第一类操作得到的向量附近寻找一个新的向量以降低图像的TV范数。在此基础上，还采用加速技术以实现算法的快速收敛。

本发明采用的技术方案是：

一种应用于不完备数据成像的快速图像重建方法，其特征在于，包括以下步骤：

S1、建立线性成像系统的离散化数学模型，并结合图像先验信息构造目标函数，将成像问题转化为最优化问题；

所述先验信息为图像的全变分范数，构造的目标函数为：

其中，x为待求图像的向量表示，为N×1向量，由对二维或三维图像按一维重新排列得到，A是M×N系统矩阵，用向量a_i表示矩阵A第i行向量的转置，b为M×1测量数据向量，用b_i表示向量b的第i个元素；||x||_TV表示x所对应图像的全变分范数，γ为平衡数据保真度约束(||Ax-b||²)与全变分范数约束(||x||_TV)的平衡系数，为正实数，取值范围一般在[0.001,10]之间；

S2、对S1得到的目标函数进行迭代求解，所述迭代求解的过程包括以下步骤：

S21、初始化：设置参数，所述参数包括容差阈值ε_i，最小松弛参数ρ_s；设定初始解向量x⁰，初始辅助向量y¹＝x⁰；并令迭代数k＝1，初始辅助标量t₁＝1；

S22、迭代过程：

S221、对辅助向量y^k执行自适应代数重建技术算子操作，得到临时估计向量u^k，所述临时估计向量u^k的表达式为：

其中，“°”表示复合函数符号，令i＝1,2,…,M，M为测量数据的个数，则

为对第i个超平面的投影算子，

为

的复合算子，表示在所有超平面的一轮依次循环投影操作；

S222、对临时估计向量u^k执行邻近映射，得到更新的向量x^k，所述向量x^k的表达式为：

上式中，prox表示实现邻近映射的邻近算子；

S223、更新辅助向量y^k，具体更新步骤如下：

S2231、按如下式子更新辅助标量t_k：

S2232、计算

S2233、如果

则y^k+1＝x^k；

其中，

由下式定义：

S224、进入下一次迭代，即k←k+1；

S23、重复步骤S221-S224，直到满足终止条件，迭代求解完成，得到重建图像。

所述步骤S221中的复合算子

是非扩展的(nonexpansive)；所述步骤S221与S222中从y^k到x^k的更新算子

也是非扩展的。非扩展性(nonexpansive)是数学上的一种定义。如果算子

是非扩展的，则对任意的x₁与x₂，满足

所述步骤S2232中，辅助向量y^k+1实际上是通过对x^k与x^k-1的外插操作得到的，发明方法中引入辅助向量y^k是为了加速算法，实现快速收敛。为了保证引入“外插操作”的加速效果，避免产生可能的发散操作，所述步骤S2233中在得到y^k+1后比较了

与

的值。如果

优于

(即

)，则保留加速操作(即外插操作)得到的y^k+1，若

没有优于

(即

)，则摒弃外插操作，令y^k+1＝x^k。

进一步的是，所述步骤S21中，设定的初始解向量x⁰可以是满足图像基本特征(如所有元素非负)的任意向量，一般情况下可简单地设为全零向量；最小松弛参数ρ_s为正实数，取值范围一般在(0,0.1)之间；容差阈值ε_i为正实数，是关于第i个测量数据偏差水平的估计，即假设第i个测量数据的偏差小于ε_i，当测量数据的偏差满足独立高斯分布且已知待求图像x对应的||Ax-b||²≤ε(ε为已知的测量误差阈值，取值范围一般为0.001M～M)时，ε_i的取值方式为：

其中M为测量数据的个数。

进一步的是，所述步骤S221中的投影算子

定义如下：

其中，

v＝b_i-<a_i,x

ρ＝max(ρ_s,min(1.5ρ'₁,1))。

进一步的是，所述步骤S222中，使用快速投影梯度算法来实现邻近映射。

进一步的是，在步骤S23中采用的终止条件为，若对于1≤i≤5，连续5次迭代，下式都成立：

则算法终止，迭代求解完成。

上述的终止条件可保证成像效果，但对于本领域普通技术人员来说，作为公知常识，理论上，在算法收敛的情况下，迭代次数越多越好，但实际情况一般不会无限制的进行算法运行，因此，当满足上述终止条件时迭代次数k过大时(即长时间运行都不能满足上述终止条件时)，可以设定迭代次数k达到最大迭代次数k_max，时算法终止，最大迭代次数k_max的合适值一般与待解决的具体问题相关，可根据实际情况选择或通过实验确定，只要重建的图像满足需要即可。

本发明的理论分析为：

代数重建法(ART)将图像重建问题转化为求解线性方程组Ax≈b，式子中用“≈”是因为考虑到实际测量数据存在偏差。设含误差(或噪声)的测量数据表示为：

b＝b*+e (1)

其中，e为M×1的噪声向量，b*为M×1的理想测量数据向量，即b*满足：

Ax＝b* (2)

代数重建法将方程组Ax＝b的各方程<a_i,x>＝b_i(1≤i≤M)看作N维空间的一个超平面，从初始解向量x⁰出发，在每次迭代，按照一定的投影顺序，依次向各超平面进行一轮循环投影操作，逐步逼近方程组的真实解。以顺序投影(依次投影的超平面序号为1,2,…,M)为例，在每次迭代，由当前估计解向量(简称当前向量)x出发，按超平面投影顺序i＝1,2,…,M，依次执行如下操作，不断更新估计向量x，更新的向量记为x′：

其中，

v＝b_i-<a_i,x> (4)

上式中的<·>为内积符号，ρ为ART算法的松弛参数，应用于噪声环境，取值范围一般在(0，1]之间。当ρ＝1时，由式(3)得到的更新向量x′是x在超平面i(<a_i,x>＝b_i)上的投影点。若测量数据不存在噪声(即e_i＝0，<a_i,x>＝b_i*)，ρ＝1时x′是由式(3)得到的离最优解x*最近的向量，从这个角度来看此时选ρ＝1时最优的。但是，由于测量数据b中含有噪声(e_i≠0)，我们希望将x投影到超平面<a_i,x>＝b_i*，而不是超平面<a_i,x>＝b_i，此时最优的松弛参数应为：

由于在实际中测量数据的偏差(或噪声)e_i是未知的，我们只能选与ρ′尽量接近的松弛参数。在假设|e_i|≤ε_i的条件下，经过数学推导，我们选用的自适应松弛参数为

ρ＝max(ρ_s,min(1.5ρ'₁,1)) (6)

其中，ρ_s为事先设定的最小松弛参数,

由于式(4)中的变量v是根据当前估计向量x与超平面i(<a_i,x>＝b_i)的距离实时计算更新的，因此，本发明方法中ART的松弛参数ρ在迭代过程是自适应调整的，并且是|v|的连续函数。

每次迭代自适应代数重建技术(ART)的算子可表示为：

其中，“°”表示复合函数符号，

为对第i个超平面的投影算子，

为

的复合算子，表示在所有超平面的一轮依次循环投影操作，i＝1,2,…,M，M为测量数据的个数；

我们将松弛参数(ρ)自适应调整的代数重建技术(ART)称为自适应代数重建技术，执行自适应代数重建技术每次迭代操作的算子称为自适应代数重建技术(ART)算子。

可以证明本发明方法中自适应ART的算子

是非扩展的(nonexpansive)，这个性质对保证本发明算法的收敛性十分重要。

由于实现邻近映射的邻近算子(proximal operator)是非扩展的，因此本发明方法中从y^k到x^k的更新算子也是非扩展的，结合步骤S2233中的单调性更新操作，保证了本发明算法的收敛性(至少不会发散)。

采用本技术方案的有益效果：

1.本发明方法结合自适应ART与邻近算子操作，并采用单调外插加速技术，显著提高了重建算法的收敛速度。

2.本发明方法结合了自适应代数重建技术，可根据测量数据的偏差水平，在算法收敛过程自适应地调整ART的松弛参数，应用灵活。

3.合理平衡测量数据信息与图像先验信息，实现不完备数据条件下的快速图像重建。

附图说明

下面结合附图对本发明作进一步详细说明。

图1为本发明方法流程示意图；

图2为实施例1中本发明的算法与对比成像算法的平方向量误差收敛曲线图；

图3为实施例1的原始图像及重建图像图；

图4为实施例2中本发明的算法与对比成像算法的平方向量误差收敛曲线图；

图5为实施例2的原始图像及重建图像图；

图6为实施例3中本发明的算法与对比成像算法的平方向量误差收敛曲线图；

图7为实施例3的原始图像及重建图像图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明作进一步阐述。

实施例1：

在实施例1中，参见附图1-3所示，应用本发明方法实现平行波束CT成像系统的快速图像重建。待重建CT图像是一个128×128的二维图像，即图3中的原始图像，对应表示成16384×1的向量x^*。一种用于平行波束CT成像系统不完备数据成像的快速图像重建方法，参见图1所示，

S1、首先建立线性成像系统的离散化数学模型Ax＝b，其中，A是维度为7200×16384的系统矩阵，即M＝7200，N＝16384，16384×1向量a_i表示矩阵A第i行向量的转置，b＝Ax*+e为7200×1的测量数据向量，b_i表示向量b的第i个元素。e为7200×1零均值高斯白噪声向量，标准方差σ＝0.1。选用图像的全变分(total variation,TV)范数作为先验信息，结合线性成像系统的离散化数学模型Ax＝b与图像先验信息(图像TV范数)建立的最小化目标函数为：

其中，平衡测量数据信息与图像先验信息的平衡系数γ＝0.003。

S2、对上述目标函数的迭代求解过程包括以下步骤：

S21、初始化。本实施例中，初始化步骤设定的参数为：容差阈值ε_i＝3σ(相当于取ε＝7200σ²)，最小松弛参数ρ_s＝0.05，初始解向量x⁰取全零向量，初始辅助向量y¹＝x⁰，辅助标量t₁＝1，并令迭代数k＝1。

S22、迭代过程：对第k次迭代重复以下步骤直到满足终止条件：

S221、对辅助向量y^k执行自适应代数重建技术算子操作，得到临时估计向量u^k，以提高估计向量对应的数据保真度：

其中，“°”表示复合操作符号，

为对第i个超平面的投影算子，

为

(i＝1,2,…,M)的复合算子，表示在所有超平面(i＝1,2,…,M)的一轮依次循环投影操作。

投影算子

定义如下：

其中，

v＝b_i-<a_i,x>.

ρ＝max(ρ_s,min(1.5ρ'₁,1))。

S222、对临时估计向量u^k执行邻近映射(proximal mapping)，得到更新的向量x^k，以减少估计向量所对应图像的TV范数：

上式中，prox表示实现邻近映射的邻近算子；使用快速投影梯度(fast projectedgradient,FGP)

算法来实现邻近映射(proximal mapping)。

S223、利用单调外插操作更新辅助向量y^k。

S2231、按如下式子更新辅助标量t_k：

S2232、计算

S2233、如果

则y^k+1＝x^k；

其中，

由下式定义：

本发明方法中引入辅助向量y^k是为了加速算法，实现快速收敛。步骤S2232中，辅助向量y^k+1实际上是通过对x^k与x^k-1的外插操作得到的。为了保证引入“外插操作”的加速效果，避免产生可能的发散操作，步骤S2233中在得到y^k+1后比较了

与

的值。如果

优于

(即

)，则保留加速操作(即外插操作)得到的y^k+1，若

没有优于

(即

)，则摒弃外插操作，令y^k+1＝x^k。

S224、进入下一次迭代，即k←k+1。

S23、重复步骤S221-S224，直到满足终止条件，迭代求解完成，得到重建图像。

这里为了满足实施例实验的实际需要，设定最大迭代次数400为终止条件(实际实验测试能满足成像效果)，即若迭代次数k达到最大迭代次数400时算法终止；但若迭代次数k达到最大迭代次数400前，对于1≤i≤5，下式都成立：

则算法也停止。

图2比较了本发明方法与其它几种经典快速成像算法的平方向量误差(10*lg(||x^k-x^*||²))收敛曲线，被比较的算法包括迭代收缩阈值算法(ISTA-FGP)，单调快速迭代收缩阈值算法(MFISTA-FGP),自适应最陡梯度凸集投影法(ASD-POCS)等。为了验证本发明算法引入辅助向量y^k的单调加速效果，图2还包含了本发明方法去掉单调外插操作更新辅助向量y^k(即将本发明方法的步骤S223内容替换为y^k+1＝x^k)的版本(标注为PAART-FGP)。

从图2可以看出，本发明方法的收敛速度明显快于其它几种成像算法。原始图像及含本发明方法在内的几种成像算法得到的重建图像如图3所示，从中可以看出本发明方法与自适应最陡梯度凸集投影法(ASD-POCS)的重建效果优于迭代收缩阈值算法(ISTA-FGP)与单调快速迭代收缩阈值算法(MFISTA-FGP)的重建效果。

实施例2：

实施例2与实施例1相比，不同之处在于测量数据偏差e的标准方差σ＝0.2，其它部分都相同。本实施例发明方法与其它几种经典快速成像算法的平方向量误差(10*lg(||x^k-x^*||²))收敛曲线如附图4所示，原始图像及几种成像算法得到的重建图像如图5所示。由图中可知，对于与实施例1不同的噪声环境，本发明方法采用与实施例1相同的参数设置，仍然可以获得良好的成像效果，收敛速度明显优于其它几种待比较的成像算法。

实施例3：

实施例3与实施例2的不同之处在于实施例3中应用本发明方法实现扇形波束(而不是平行波束)CT成像系统的快速图像重建。待重建CT图像与实施例1相同，也是128×128的二维图像；线性成像系统离散化数学模型Ax＝b中的A，b数值与实施例1、2不同，但A，b的维度与实施例1、2相同，即A是7200×16384的系统矩阵，b＝Ax*+e为7200×1的测量数据向量；e为7200×1零均值高斯白噪声向量，标准方差σ＝0.2。

本实施例发明方法与其它几种经典快速成像算法的平方向量误差(10*lg(||x^k-x^*||²))收敛曲线如附图6所示，原始图像及几种成像算法得到的重建图像如图7所示。由图中可知，对于扇形波束CT成像系统，本发明方法采用与实施例1平行波束成像系统相同的参数设置，仍然可以获得良好的成像效果，收敛速度明显优于其它几种待比较的成像算法。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (4)

1.一种应用于不完备数据成像的快速图像重建方法，其特征在于，包括以下步骤：

S1、建立线性成像系统的离散化数学模型，并结合图像先验信息构造目标函数，将成像问题转化为最优化问题；

所述先验信息为图像的全变分范数，构造的目标函数为：

其中，x为待求图像的向量表示，为N×1向量，由对二维或三维图像按一维重新排列得到，A是M×N系统矩阵，用向量a_i表示矩阵A第i行向量的转置，b为M×1测量数据向量，用b_i表示向量b的第i个元素；||x||_TV表示x所对应图像的全变分范数，γ为平衡数据保真度约束与全变分范数约束的平衡系数，为正实数，取值范围：0.001≤γ≤10；

S2、对S1得到的目标函数进行迭代求解，所述迭代求解的过程包括以下步骤：

S21、初始化：设置参数，所述参数包括容差阈值ε_i，最小松弛参数ρ_s；设定初始解向量x⁰，初始辅助向量y¹＝x⁰；并令迭代数k＝1，初始辅助标量t₁＝1；

S22、迭代过程：

S221、对辅助向量y^k执行自适应代数重建技术算子操作，得到临时估计向量u^k，所述临时估计向量u^k的表达式为：

其中，

表示复合函数符号，令i＝1,2,…,M，M为测量数据的个数，

则

为对第i个超平面的投影算子，

为

的复合算子，表示在所有超平面的一轮依次循环投影操作；

S222、对临时估计向量u^k执行邻近映射，得到更新的向量x^k，所述向量x^k的表达式为：

上式中，prox表示实现邻近映射的邻近算子；

S223、更新辅助向量y^k，具体更新步骤如下：

S2231、按如下式子更新辅助标量t_k：

S2232、计算

S2233、如果

则y^k+1＝x^k；

其中，

由下式定义：

S224、进入下一次迭代，即k←k+1；

S23、重复步骤S221-S224，直到满足终止条件，迭代求解完成，得到重建图像；

在步骤S23中采用的终止条件为：若对于1≤i≤5，连续5次迭代，下式都成立：

则算法终止，迭代求解完成。

2.根据权利要求1所述的一种应用于不完备数据成像的快速图像重建方法，其特征在于，所述步骤S21中，设置参数的方法为：初始解向量x⁰设为全零向量；最小松弛参数ρ_s为正实数，取值范围：0＜ρ_s＜0.1；容差阈值ε_i为正实数，是关于第i个测量数据偏差水平的估计，将测量数据偏差看作满足独立高斯分布，则ε_i的取值方式为：

其中ε为已知的测量误差阈值，取值范围为0.001M～M。

3.根据权利要求1所述的一种应用于不完备数据成像的快速图像重建方法，其特征在于，所述步骤S221中的投影算子

定义如下：

其中，

v＝b_i-<a_i,x>

ρ＝max(ρ_s,min(1.5ρ'₁,1))。

4.根据权利要求1所述的一种应用于不完备数据成像的快速图像重建方法，其特征在于，所述步骤S222中，使用快速投影梯度算法来实现邻近映射。

Images