一种应用于ct成像的快速代数重建方法
【技术领域】
[0001] 本发明涉及一种应用于CT成像的快速代数重建方法，属于生物医学成像、无损检 测等技术领域。
【背景技术】
[0002] 计算机断层扫描(CT)成像作为一种重要的无损检测技术，在医学领域及工业领域 都有广泛应用。目前常用的CT图像重建算法主要有两大类:一类是以滤波反投影算法(FBP) 为代表的解析法;另一类是以代数重建法(ART)为代表的迭代算法。解析法的优点是计算量 小，重建速度快，当投影数据充足时能获得良好的重建图像质量，因此现今主流商业CT仍采 用解析类方法重建图像。
[0003] 在实际医学应用和工业应用中，受检测环境、成本、时间、人员健康等因素限制，经 常遇到投影数据不完备的情形。在不完全投影情况下，现有的滤波反投影类算法难以重建 出完整的CT图像。与解析法不同，代数重建算法将图像重建问题转化为大规模线性方程组 求解问题，通过结合一些正则化先验知识，在投影数据不完全时仍可重建高质量图像。另 外，代数重建算法不易受所处理问题模型的限制并且可以引入被检物体的先验信息，容易 扩展到CT应用的其他领域，展示出了很大的发展空间与应用潜力。
[0004]代数重建法的缺点是计算量大，重建速度慢，这在很长一段时间内限制了它的应 用。但随着并行计算理论及计算机硬件技术的发展，此类技术瓶颈已得到很大程度的缓解， 代数重建法的优势更加凸显，近年来代数重建技术重新引起众多学者的注意。如何在不完 全投影情况下提高代数重建法的收敛速度与精度，成为了代数重建算法能否在实际成像系 统中得到普遍应用的关键。
【发明内容】
[0005] 本发明提供一种可应用于CT成像的快速代数重建方法。相对于传统的代数重建方 法(ART)，该算法收敛速度快，可以用于不完全投影情况下的快速CT图像重建。
[0006] 本发明的理论分析为：
[0007] 代数重建法(ART)将图像重建问题转化为如下大规模线性方程组求解问题：
[0008] Ax = b (1)
[0009] 其中，X即为待求的图像，为NX 1向量，由对二维或三维图像按一维重新排列得到。 X的每个元素 Xl对应原图像的一个像素。A是Μ X N系统矩阵，每个元素描述第j个像素对第 i条射线的贡献。b为MX 1测量数据向量，其第i个元素 h对应第i条射线的衰减。
[00?0] 最早的ART方法由Kaczmarz提出，也称Kaczmarz方法（简称KA)。其基本思想是：将 各方程看作N维空间的一个超平面，从初始解向量X*3出发，依次在各超平面上投影，逐步逼 近方程组的真实解，可用如下数学公式表示：
[0011] xk+1 = P(b,xk) =PM〇---〇P2°Pi(b,xk) (2)
[0012] xk'1 = Pi(b，xk'j)=xk'j+A1' :iai (3)
[0014] xk'° = xk,xk+1 = xk'M (5)
[0015] 其中，ai为矩阵A的第i行向量，〈·>与Μ · I I分别表示内积与欧拉范数，aiT表示对 的转置表示将向量投影到第i个超平面的操作，表示KA第k次迭代时投 影到第j个超平面的向量。每次迭代，对向量xki行操作P(b，xk)，即依次在各超平面上进行 一轮投影，得到下一代的解向气
[0016] 为了提高Kaczmarz方法的效率，本发明在每次迭代对KA投影到第i个超平面的投 影向量xk〃沿其与前代该超平面投影向量X 1-1，1连线方向进行进一步的调整优化，使得调整 后的新向量在X k-1与Xk>1连线方向上最优(到真实解向量，的平方距离或误差最小）。真 实解向量，为所有超平面的交点，调整优化后得到的向量X e应满足：（x*-xe)丄（χ?-χΗ'1)。 几何示意图如附图图1所示。
[0017] 本发明所提快速代数重建方法的迭代数学公式描述如下：
[0020] xe = Pi(b，xk'j)-xk-hLxb-xk-U (8)
[0021] 该算法的关键在于如何求得最优的β。利用的上述垂直关系，有：
[0022] dM=||x*_xM||2=||x*_ xC||2+i32||xe||2 ⑶
[0023] dk-14= | |x*-xk-14| |2= | |x*-xc| |2+(β+ι)2| |xe| |2 (1〇)
[0024] Ad = dk-14-(1Μ=(2β+1)| |xe| I2 (11)
[0025]
[0026] 其中，d15，1表示向量x15，1到真实解向量x*的平方误差(或距离）。根据式(3)、（4)，KA 中将向量xk'」投影到第i个平面时，得到的xk"满足：以"-分"丄以"-林八因此有：
[0027] dk>i = dkJ-| |xk>i-xkJ| |2 = dkJ-| I2 (13)
[0028] 这样由一个初始值d*3开始= 根据式（13)依次迭代，即可求出相邻两次迭 代dM相对于dk~的减小值：
[0029] Ad = dk-14-(1Μ (14)
[0030] 在式（14)中，d*3会被消去。所以虽然d*3的值未知，但在编程实现时可将d*3设为任意 一个常量C，式（14)的结果不会发生改变，也不影响算法的迭代计算。将式（14)计算得到的 Ad代入式（12)即可求出β，再将求出的β代入式(7)求出。然后根据式（15)、（16)更新d，1 与dM:
[0031] xk'i = xc (15)
[0032] dk>i = dk>i-e2| |xe| I2 (16)
[0033] 需要指出的是，一般情况下编程设定的W关I |#-/| |2,这样算法迭代得到的dk'1 与| |Χ*-ΧΜ| |2相差一个常量，但Ad=dk-|x*-xk-Ul |2-| |χ*-χΜ| |2成立。
[0034] 本发明所提的快速代数重建法中，每次迭代可以采用不同的超平面投影顺序。仿 真与初步分析表明，在每次迭代采用随机的投影顺序，在许多情况下可以明显提高算法的 效率。
[0035] 本发明所提的算法，每个超平面i需要记录前一次迭代的投影向量X1-1，1。对于CT成 像应用，图像的尺寸（即向量A 1'1的维数)及超平面的个数（即测量数据的个数)很大。这种 情况下，记录每个超平面的前一次迭代向量3+ 1需要很大的存储空间。为了降低存储要求， 可以推广本算法，只选择少量超平面对上面的投影向量进行进一步调整优化。这样可以 在算法加速程度与所需额外存储空间进行平衡，在增加少量存储空间的条件下加快算法的 收敛速度。幸运的是，对于本发明算法，这种方式很容易实现，只需事先选定一个进行加速 调整操作的超平面子集Ω即可，迭代更新公式不变。
[0036] 本发明所提的快速代数重建方法的技术方案包括以下步骤：
[0037] 步骤一:初始化。设定k = 0,设定初始解ALA (1^ = ^设置进行加速调整操作 的超平面子集Ω。
[0038] 步骤二:迭代过程。重复以下步骤直到算法满足收敛条件：
[0039] 1.在第k次迭代，首先设定超平面投影顺序，用idx(l)表示第1个投影超平面的索 引下标。
[0040] 2.在第k次迭代，对每一个超平面i = idx(l)，1 = 1，2，. . .M，依次执行如下操作：
[0041] 1)将当前向量投影到超平面i，得到向量士1。具体操作如下：
[0042] a)j = idx(l_l)
[0044] c)xk,1 = xk,:i+Aai
[0045] d)dk>i = dkJ-| |Aai| I2
[0046] 2)如果k>0并且ie Ω，则对超平面i上的向量xk，i执行如下调整操作：
[0047] a)Ad = dk"1>i-dk>i
[0048] b)xe = xk>i-xk"1>i
[0050] d)xk>i = xk>i+0xe
[0051] e)dk>i = dk>i-e2| |xe| I2
[0052] 3.更新向量，迭代次数加1。具体操作如下：
[0053] xk+1 = xk'M,dk+1 = dk,M
[0054] xk+1'0 = xk+1，dk+1'0 = dk+1
[0055] k = k+l
[0056] 本发明可以选择较少的超平面对其投影向量执行加速调整操作，并且在每次迭代 采用随机投影顺序，这样可减小额外存储量，缩短收敛时间。
【附图说明】
[0057]图1投影向量调整几何示意图。
[0058] 表1迭代过程的平方向量距离(或误差）
【具体实施方式】
[0059] 下面将结合附图和【具体实施方式】，对本发明进行进一步的说明。
[0060] 假设CT成像问题转化得到的大规模线性方程组为Ax = b。其中，A是MXN系统矩阵， Μ为测量数据的维数，N = nXn为图像空间的维数。在不完全投影情况下，M〈N。为了重建高质 量图像，还需要引入一些先验知识进行正则化。一种有效且广泛使用的正则化技术是全变 分(Total variance,TV)正则化技术。
[0061] 本发明提供一种类似TV的正则化技术，正则化后CT图像重建问题转化为如下目标 函数的最小化问题：