CN105608719B - 一种基于两阶段投影调整的快速ct图像重建方法 - Google Patents

Abstract

本发明提供了一种基于两阶段投影调整的快速CT图像重建方法。针对传统代数重建法无法精确求解非一致线性方程组最小二乘解的缺点，本发明将非一致线性方程组的求解转化为两个一致线性方程组的求解，并利用基于投影调整的快速代数重建算法分别求解两个一致线性方程组，最终实现非一致线性方程组最小二乘解的快速精确求解。为进一步减小算法计算量，本发明采用平方方程误差描述第一个一致方程组迭代解的逼近程度，据此将算法分为前后两个收敛阶段。在第一阶段，只更新第一个一致方程组的解；在第二阶段，分别更新两个一致方程组的解。与传统代数重建法相比，本发明算法能够快速精确求解最小二乘解，可用于不完全投影情况下的快速CT图像重建。

Description

一种基于两阶段投影调整的快速CT图像重建方法

技术领域

本发明涉及一种基于两阶段投影调整的快速CT图像重建方法，属于生物医学成像、无损检测等技术领域。

背景技术

计算机断层扫描(CT)成像是一种重要的无损检测技术，经过几十年的发展，已经从医学领域扩展到工业等其他领域。随着临床应用与工业应用范围的日益扩大与深入，受检测环境、成本、时间、人员健康等诸多因素限制,不完全投影数据下的CT图像重建已不可避免。在现今主流商用CT机中，以滤波反投影算法(FBP)为代表的解析法是最常用的重建算法。该类算法的优点是计算量小，重建速度快，但要求获得充足的投影数据。在不完全投影情况下,现有的滤波反投影类算法难以重建出高质量的CT图像。

与解析法不同，代数重建算法将CT图像重建问题转化为大规模线性方程组求解问题，通过结合一些正则化先验知识，在投影数据不完全时仍可重建高质量图像。另外，代数重建算法不易受所处理问题模型的限制并且可以引入被检物体的先验信息，容易扩展到CT应用的其他领域，展示出了很大的发展空间与应用潜力。

虽然世界上第一台CT机所采用的图像重建算法就是代数重建法，但计算量大、重建速度慢的缺点，在很长一段时间内限制了代数重建法在CT成像系统的大规模应用。目前随着并行计算理论及计算机技术的快速发展，代数重建法的技术瓶颈已得到很大程度的缓解，其优势更加凸显，近年来代数重建技术重新引起了众多研究者的注意。如何在不完全投影情况下提高代数重建法的收敛速度与精度，成为了代数重建算法能否在实际成像系统中得到普遍应用的关键。

发明内容

本发明提供一种基于两阶段投影调整的快速CT图像重建方法，达到收敛速度快，能够精确求解非一致线性方程组的最小二乘解，可用于不完全投影情况下的快速CT图像重建的目的。

本发明的理论分析为：

代数重建法(ART)将图像重建问题转化为如下大规模线性方程组求解问题：

Ax＝b (1)

其中，x即为待求的图像，为N×1向量，由对二维或三维图像按一维重新排列得到，x的每个元素x_i对应原图像的一个像素。A是M×N系统矩阵，每个元素a_ij描述第j个像素对第i条射线的贡献。b为M×1测量数据向量，其第i个元素b_i对应第i条射线的衰减。

最早的ART方法由Kaczmarz提出，也称Kaczmarz方法(简称KA)。其基本思想是：将各方程看作N维空间的一个超平面，从初始解向量x⁰出发，依次在各超平面上投影，逐步逼近方程组的真实解，可用如下数学公式表示：

x^k+1＝P(b,x^k)＝P_Mο…οP₂οP₁(b,x^k) (2)

x^k,i＝P_i(b,x^k,j)＝x^k,j+λ^i,ja_i (3)

x^k，0＝x^k,x^k+1＝x^k,M (5)

其中，a_i为矩阵A的第i行向量，〈·〉与||·||分别表示内积与欧拉范数，a_i ^T表示对a_i的转置。P_i(b,x^k,j)表示将向量x^k,j投影到第i个超平面的操作，x^k,j表示KA第k次迭代时投影到第j个超平面的向量。每次迭代，对向量x^k执行操作P(b,x^k)，即依次在各超平面上进行一轮投影，得到下一代的解向量x^k+1＝P(b,x^k)＝x^k,M。

为了提高Kaczmarz方法的效率，我们提出了基于投影调整的快速代数重建算法。在每次迭代对KA投影到超平面i的投影向量x^k,i沿其与前代该超平面投影向量x^k-1,i连线方向进行进一步的调整优化，使得调整后的新向量x^c在x^k,i-1与x^k,i连线方向上最优(即到真实解向量x^*的平方距离或误差最小)。基于投影调整的快速代数重建算法每次迭代的数学公式如下：

x^e＝P_i(b,x^k,j)-x^k-1,i＝x^k,i-x^k-1,i (8)

式(7)中的参数β根据如下表达式计算：

△d＝d^k-1,i-d^k,i (10)

d^k,i＝||x*-x^k,i||²,d^k-1,i＝||x*-x^k-1,i||² (11)

其中，d^k,i表示向量x^k,i到真实解向量x^*的平方误差(或距离)。在KA中，将向量x^k,j投影到第i个平面得到x^k,i，对应的d^k,i满足如下关系：

d^k,i＝d^k,j-||x^k,i-x^k,j||²＝d^k,j-||λ^i,ja_i||² (12)

这样由一个初始值d⁰开始(d^k,0＝d^k)，根据式(12)依次迭代计算d^k-1,i与d^k,i，再将d^k-1,i与d^k,i代入式(10)，可消去d⁰，求得Δd。将式(10)计算得到的△d代入式(9)求出β，再将求出的β代入式(7)即可求出调整后的x^c。然后根据式(13)、(14)更新x^k,i与d^k,i：

x^k,i＝x^c (13)

d^k,i＝d^k,i-β²||x^e||² (14)

上述快速代数重建算法中每次迭代可以采用不同的超平面投影顺序。在每次迭代采用随机的投影顺序，许多情况下可以明显提高算法的效率。本发明后面的例子中，也将采用随机投影顺序的快速代数重建算法。

理论分析表明，如果待求解线性方程组(1)是一致线性方程组(即方程组存在真实解)，则KA将收敛到方程组的解。若初始解x⁰属于A^T的值域R(A^T)，KA收敛的解是方程组的最小二范数解。若方程组(1)是非一致线性方程组(即方程组不存在解)，则KA最终只能收敛到一个区域，而不是某个特定解。基于投影调整的快速代数重建法在求解一致线性方程组的收敛性与KA相同，在非一致方程组情况下也不能收敛到某个特定解。

由于噪声及正则化的影响，在很多情况下基于代数重建法的CT图像重建问题被转化为非一致线性方程组求解问题。由于非一致线性方程组Ax＝b无解，我们转而求该方程组均方误差||Ax-b||²最小的解，即最小二乘解。为此，本发明将非一致线性方程组的求解转化为两个一致线性方程组的求解，并利用基于投影调整的快速代数重建算法分别求解两个一致线性方程组，最终实现非一致线性方程组最小二乘解的精确求解与CT图像的快速重建。其理论依据与实施方法如下。

若方程组(1)为非一致线性方程组，则将方程组右边的向量b分解为互相垂直的两个分量：

其中，b_R(A)为b在矩阵A的列空间或值域空间R(A)的投影分量，为b在矩阵A^T的零空间N(A^T)的投影分量。假设x^*为非一致线性方程组的最小二乘解，则x^*满足：

Ax*＝b_R(A) (16)

而实际是如下方程的解：

A^Ty＝0 (17)

其中，y为M×1向量。方程组(16)与(17)都是一致线性方程组，都可以用基于投影调整的快速代数重建算法求精确解。因此我们可以以b作为y的初始解y⁰，先求出式(17)的解然后利用求出b_R(A)，最后再解方程组(16)求出最小二乘解x^*。为加速求解过程，我们交互地迭代求解y与x。每次迭代的数学公式如下：

b^k+1＝b-y^k+1 (19)

在算法起始收敛阶段，由式(20)迭代产生的x^k+1实际上是收敛到方程组(21)的解而不是方程组(16)对应的最小二乘解x^*。

Ax＝b^k+1 (21)

但是，随着y^k+1逐渐收敛到b^k+1将逐渐收敛到b_R(A)。这样方程组(21)就变为了(16)，因此x^k最终仍可以收敛到与(16)对应的最小二乘解x^*。但是在算法起始收敛阶段，x^k的收敛方向可能偏离x^*。

为进一步减小算法计算量与加快算法收敛速度，本发明将算法分为前期与后期两个阶段。在y^k离距离较远的前期(或第一)阶段，我们只更新y^k；在y^k离距离较近时的后期(或第二)阶段，我们再分别更新y^k与x^k。由于该策略在算法前期阶段只更新y^k而不需要更新x^k，减小了计算量。同时，对于在算法开始阶段x^k的收敛方向较严重偏离x^*的情况，该策略还可以提高算法的收敛速度。

为了判断算法当前迭代属于第一还是第二收敛阶段，本发明采用式(22)定义的平方方程误差ε^k来描述y^k的收敛程度：

上式中，将方程组的平方误差除以向量b的长度M是为了让这个定义可以适用于不同规模的场景。事先定义一个阈值ε_T，若当前迭代的方程组平方误差ε^k>ε_T，则认为算法当前属于第一阶段，只更新y^k；若ε^k≤ε_T，则认为算法当前属于第二阶段，分别更新y^k与x^k。

本发明所提的适用于非一致线性方程组最小二乘解精确求解的两阶段投影调整快速代数重建算法的技术方案包括以下步骤：

步骤一：初始化：设定k＝0，设定初始解x⁰，y⁰＝b及d⁰，ε_T。

步骤二：迭代过程：重复以下步骤直到算法收敛：

1)k＝k+1；

2)根据式(18)更新y^k。

3)根据式(22)计算ε^k

4)如果ε^k≤ε_T，执行以下操作：

a)根据式(19)更新b^k。

b)根据式(20)更新x^k。

本发明提供一种基于两阶段投影调整的快速CT图像重建算法。相对于传统的代数重建方法(ART)，该算法收敛速度快，能够精确求解非一致线性方程组的最小二乘解，可用于不完全投影情况下的快速CT图像重建。

附图说明

图1待重建的原始CT图像；

图2本发明算法的平方方程误差(dB)收敛曲线；

图3本发明算法重建的CT图像。

具体实施方式

下面将结合附图和一个平行波束CT图像重建例子的具体实施方式，对本发明进行进一步的说明。

待重建CT图像是一个100×100的二维图像(见附图图1)，对应表示成10000×1的向量x*，A是6000×10000的系统矩阵，b＝Ax*+e为6000×1的测量数据向量。e为零均值高斯白噪声向量，强度为b的5％。本例中M＝6000，N＝10000，属于不完全投影情况。为了重建高质量图像，还需要引入一些先验知识进行正则化。一种有效且广泛使用的正则化技术是全变分(Total variance,TV)正则化技术。

本发明提供一种类似全变分(TV)的正则化技术，正则化后CT图像重建问题转化为如下目标函数的最小化问题：

f(x)＝||Ax-b||²+γ²||Rx||² (23)

其中，γ为正则化系数。R＝R₁+R₂为N×N矩阵。假设原100×100二维图像按列重排成N×1向量x*，令h_r(i)表示原图像第i个像素右边相邻的像素下标，v_b(i)表示第i个像素下方相邻的像素下标。矩阵R₁的元素定义如下：

矩阵R₂的元素定义如下：

上式中，mod(i,n)表示i除以n的余数。

最小化目标函数(23)等效于求如下非一致线性方程组的最小二乘解：

其中，

接下来采用本发明提出的基于两阶段投影调整的快速代数重建算法求解非一致线性方程组(26)，重建CT图像。求解步骤如下：

步骤一：初始化：设定k＝0，x⁰＝0，d⁰＝0，ε_T1＝1。

步骤二：迭代过程。分为前期与后期两个收敛阶段：

1)由附图图2可知，k≤4时，ε^k>ε_T，处于第一收敛阶段。算法只根据式(18)更新y^k，不更新x^k。因此图中的平方方程误差曲线在前期维持不变。

k>4时,ε^k≤ε_T，进入第二收敛阶段。算法分别根据式(18)、(19)、(20)更新y^k、b^k与x^k。到第8代时，算法已经收敛。本发明算法重建的CT图像如附图图3所示。

Claims (2)

1.一种基于两阶段投影调整的快速CT图像重建方法，该方法基于快速代数重建法求解非一致线性方程组Ax＝b的最小二乘解；其中x即为待求的图像，为N×1向量，由对二维或三维图像按一维重新排列得到，x的每个元素x_i对应原图像的一个像素；A是M×N系统矩阵，每个元素a_ij描述第j个像素对第i条射线的贡献；b为M×1测量数据向量，其第i个元素b_i对应第i条射线的衰减；将非一致线性方程组Ax＝b的求解转化为两个一致线性方程组的求解，并基于两阶段加速策略用快速代数重建法分别求解两个一致线性方程组；包含以下步骤：

步骤1：将向量b分解为互相垂直的两个分量：

其中，b_R(A)为b在矩阵A的列空间或值域空间R(A)的投影分量，为b在矩阵A^T的零空间N(A^T)的投影分量；

步骤2：设x为非一致线性方程组的最小二乘解，则x满足：

Ax＝b_R(A)；

则为如下方程的解：

A^Ty＝0

先求解一致方程组A^Ty＝0得到然后利用求出b_R(A)，最后再解一致方程组Ax＝b_R(A)，求得的解x即为非一致线性方程组的最小二乘解；

步骤3：以b作为y的初始解y⁰，采用快速代数重建法交互地迭代求解两个一致线性方程组，最终收敛的x即为非一致线性方程组的最小二乘解；每次迭代的数学公式为：

y^k+1＝P(0,y^k)

b^k+1＝b-y^k+1

x^k+1＝P(b^k+1,x^k)；

其中：k表示迭代数，P(0,y^k)表示采用基于投影调整的快速代数重建法求解方程组A^Ty＝0的迭代更新算子，P(b^k+1,x^k)表示采用基于投影调整的快速代数重建法求解方程组Ax＝b^k+1的迭代更新算子；

其特征在于，所述步骤3交互地迭代求解y与x的过程分为前期与后期两个阶段，其中前期阶段只更新y^k；在后期阶段分别更新y^k与x^k。

2.如权利要求1所述的一种基于两阶段投影调整的快速CT图像重建方法，其特征在于采用当前迭代线性方程组A^Ty＝0的平方方程误差ε^k来划分前期与后期阶段，ε^k的计算公式为：

根据经验事先设定阈值ε_T，若当前迭代的方程组平方误差ε^k＞ε_T，则认为算法当前属于前期阶段；若ε^k≤ε_T，则认为算法当前属于后期阶段。