基于有序子集加速拆分算法的三维CT图像重建

谌湘倩 马绍惠 须文波

摘 要： 针对计算机断层扫描（CT）重建过程中统计方法计算时间较长的问题，提出一种利用有序子集加速拆分算法的三维CT图像重建方法。该方法充分利用线性约束凸优化问题的增广拉格朗日（AL）方法在较弱条件下的收敛速度快的优势；同时针对内部最小二乘问题，使用AL方法的线性变形求解加权正则化最小二乘问题，该方法使用可分离二次型代理函数代替缩放增广拉格朗日中的二次型AL惩罚项，得到一种简单有序子集（OS）加速型拆分算法（OS?ASA），避免了繁琐的参数调整，可快速收敛。实验结果表明，该文算法显著加快了CT图像重建的收敛速度，当使用子集较多时，CT图像重建可以减少OS伪影。

关键词： 三维图像重建； 计算机断层扫描； 拉格朗日； 有序子集； 最小二乘

中图分类号： TN919?34； TP391 文献标识码： A 文章编号： 1004?373X（2016）06?0104?04

Three?dimensional CT image reconstruction based on accelerated splitting algorithm

of ordered subsets

CHEN Xiangqian1， MA Shaohui1， XU Wenbo2

（1. Department of Computer Science and Technology， Henan Institute of Technology， Xinxiang 453002， China；

2. College of Information Engineering， Jiangnan University， Wuxi 210034， China）

Abstract： Since the computing time of statistical method in CT （computed tomography） reconstruction process is long， a three?dimensional CT image reconstruction method based on accelerated splitting algorithm of ordered subsets （OS） is proposed. The method takes full advantage of the fast convergence speed when the augmented Lagrangian （AL） method of the linear constraint convex optimization is in weak condition. The weighted and regularized least square problems are solved with linear variant of AL method. The separable quadratic surrogate function is used in this method to replace the quadratic AL penalty term scaled in the AL to obtain a simple accelerated splitting algorithm of ordered?subsets （OS?ASA）， which can avoid the tedious parameter tuning and speed up the convergence rate. The experimental results show that the algorithm significantly accelerates the convergence speed of CT image reconstruction. The OS artifacts can be reduced by CT image reconstruction when the more subsets are used.

Keywords： three?dimensional image reconstruction； computed tomography； Lagrangian； ordered subset； least square

0 引 言

由于获取较低剂量X射线的CT扫描具有保持图像质量的潜力，图像重建的统计方法已广泛应用于计算机断层扫描（Computed Tomography，CT）[1]，然而，统计方法计算时间长是实践应用中的最大瓶颈。为了加速统计方法，已经有很多学者研究其优化技术，例如，文献[2]的增广拉格朗日（Augmented Lagrangian，AL）方法（包括交替方向变型）使用两个拆分求解正则化逆问题，AL方法通过引入辅助变量可以分离非光滑[l1]正则项，得到简单惩罚型最小二乘内部问题，再使用快速傅里叶变换（FFT）算法和软阈值[l1]范数的近端映射能够有效解决该问题。然而，类似于X射线CT图像重建的应用程序中，统计加权的巨大动态范围造成了Hessian矩阵高度平移，导致内部最小二乘问题非常具有挑战性[1]。为了解决该问题，文献[3]引入了分离平移变量和加权二次数据拟合项的平移不变量，使得预处理共轭梯度（Preconditioned Conjugate Gradient，PCG）方法能有效解决宽松条件下的最小二乘问题。实验结果表明，该方法显著加速了2?D CT[3]。然而，在具有锥形束几何形状的3?D CT中，为内部最小二乘问题构建预处理器更加困难。有序子集（Ordered?subsets，OS）算法[4]是一种具有对角预处理器的一阶方法，使用较小的步长，但非常简单，适用于3?D CT。通过集合投影为满足“子集平衡条件”的M个有序子集，并使用M个子集梯度增量更新图像，与标准梯度下降方法的每次外部迭代的图像更新一样多，从而导致标准梯度下降法早期迭代M次加速。当在约束条件下随机选择子集使子集梯度无偏置且具有有限方差时，也可以称为随机梯度方法[5]。当M增加时，快速OS算法有“更大的”限制周期，并在重建图像中出现伪影，但它可以通过使用松弛动量减小，即增长的对角优化器（或等价地，递减步长），以较慢的收敛速度为代价[6]去除伪影。文献[7]的研究表明，加速近端梯度方法在梯度误差和近侧映射计算中对误差更敏感。

本文使用线性AL方法（Linear Augmented Lagrangian Method，LALM）求解正则化（加权）最小二乘问题，优化了二次AL惩罚项，代替了平滑数据拟合项，在缩放增广拉格朗日中，使用固定对角线优化器，产生更加简单的OS加速型拆分算法（OS?ASA）。为了进一步加速，使用二阶递归系统分析设计了确定性向下连续化方法，避免了繁琐的参数调整，可快速收敛。实验结果表明，提出的算法在早期迭代中显著加速了X射线CT图像重建。

1 本文提出的OS?加速拆分算法

首先考虑一个普通复合凸优化问题：

[x∈argminx{g（Ax）+h（x）}] （1）

式（1）的等效约束最小化问题为：

[（x，u）∈argminx，u{g（u）+h（x）} s.t. u=Ax] （2）

式中：[g]和[h]是封闭且恰当的凸函数；CT中，[A]表示系统（投影）矩阵；[x]表示待构建图像；[g]表示一个加权二次数据拟合项；[h]表示边缘保持正则项。求解约束最小化问题（2）的一种方式是使用AL方法（交替方向），交替最小化缩放后的增广拉格朗日：

[LA（x，u，d；ρ）?g（u）+h（x）+ρ2Ax-u-d22] （3）

相对于[x]和[u]，按照[d]的梯度增加，产生下列AL迭代[2]：

[x（k+1）∈argminxh（x）+ρ2Ax-u（k）-d（k）22u（k+1）∈argminug（u）+ρ2Ax（k+1）-u-d（k）22d（k+1）=d（k）-Ax（k+1）+u（k+1）] （4）

式中：[d]为拆分变量[u]缩放后的拉格朗日乘数；[ρ>0]为对应的AL惩罚参数。

线性AL方法（LALM）[8]在式（4）的x更新中替代二次AL惩罚项：

[θk（x）?ρ2Ax-u（k）-d（k）22] （5）

可将二次代理（SQS）函数分离为：[θk（x；x（k））?θk（x（k））+?θk（x（k）），x-x（k）+ρL2x-x（k）22 =ρ2tx-（x（k）-tA′（Ax（k）-u（k）-d（k）））22+ （constant independent of x）] （6）

式中，[L>A22=λmax（A′A）]确保[LI-A′A>0]和[t?1L]，该函数满足“优化”条件：

[θk（x；x）≥θk（x）， ?x，x∈Dom θkθk（x；x）=θk（x）， ?x∈Dom θk] （7）

可很容易地将[L]泛化为对称半正定矩阵[S+]，例如，基于OS的算法中所用的对角矩阵[7][DL]，仍能确保式（7）。当[S+=A′A]时，LALM则为标准AL方法，优化对角矩阵产生更加简单的[x]?更新，对应的LALM迭代如下[8]：

[x（k+1）∈argminx?k（x）?h（x）+θk（x；x（k））u（k+1）∈argminug（u）+ρ2Ax（k+1）-u-d（k）22d（k+1）=d（k）-Ax（k+1）+u（k+1）] （8）

将[g]限制为二次损失函数，即[g（u）?（12）y-u22]，则最小化问题（1）变为正则化最小二乘问题：

[x∈argminxΨ（x）?12y-Ax22+h（x）] （9）

令[L（x）?g（Ax）]为式（9）的二次数据拟合项。假设[L]适合利用OS加速，即[L]可分解为[M]个较小的满足“子集平衡条件”[7]的二次型函数[L1，L2，…，LM]：

[?L（x）≈ML1（x）≈ML2（x）≈…≈MLM（x）] （10）

有助于子集梯度近似[L]的梯度。由于[g]为二次型，其近端映射为线性，故简单闭合形式解可表示为：

[u（k+1）=ρρ+1（Ax（k+1）-d（k））+1ρ+1y] （11）

[d]的更新如下：

[u（k+1）+ρd（k+1）=y] （12）

若初始化[d]为[d（0）=ρ-1（y-u（0））]。令[u?u-y]为拆分剩余差，则简化的LALM迭代如下所示（求x的迭代）：

[s（k+1）=A′（ρ（Ax（k）-y）+（1-ρ）u（k））x（k+1）∈prox（ρ-1t）h（x（k）-（ρ-1t）s（k+1））u（k+1）=ρρ+1（Ax（k+1）-y）+1ρ+1u（k）] （13）

每次迭代式（13），其计算复杂度减少，为一个[A]的乘法、一个[A′]的乘法和一个[h]的近端映射，可非迭代求解或不使用[A]，[A′]迭代求解。由于[L]的梯度为[A′（Ax-y）]，令[g?A′u]（拆分剩余的逆投影）表示拆分梯度，式（13）可重写为：

[s（k+1）=ρ?L（x（k））+（1-ρ）g（k））x（k+1）∈prox（ρ-1t）h（x（k）-（ρ-1t）s（k+1））g（k+1）=ρρ+1?L（x（k+1））+1ρ+1g（k）] （14）

式（14）称为基于梯度的LALM，因为仅[L]的梯度用于执行更新，每次迭代式（14）的计算复杂度变为[L]的梯度评估和[h]的近端映射。

基于梯度的LALM为正则化最小二乘成本函数[ψ]的广义近端梯度下降法，其中，步长为[ρ-1t]、搜索方向为[s（k+1）]，即[L]的梯度和拆分梯度的加权平均。[ρ]越小，产生的步长越大，当[ρ=1]时，式（14）变为近端梯度方法或迭代收缩/阈值算法（ISTA）。即利用LALM，通过减小[ρ]可任意增加近端梯度方法的步长。假设所有更新均准确，即对于所有[k]，[εk=0]。由式（12）得：

[-ρd（k）=u（k）-y→Ax-y=z，k→∞（-ρd（0））-z=u（0）-Ax]

合理的初始化，即[u（0）=Ax（0）]，因此，[g（0）=?L（x（0））]，可重写为[ρ]的函数：

[C（ρ）=x（0）-x22ρ-1t+A（x（0）-x）22ρ] （15）

当：

[ρopt=A（x（0）-x）2x（0）-x2≤1] （16）

该常量获得最小值；当[L?A22]，[ρopt?1]且满足第一项控制[C]，使[ρopt<ρ≤1]。由于[ρρopt?1]，原始对偶差的上界变为：

[Ω（zk，x）-Ω（z，xk）≤C2k≈L2x（0）-x22ρ-1k] （17）

即相比近端梯度方法（[ρ=1]），利用[ρ-1]的因子提升了式（14）的收敛速度（边界），[ρopt<ρ≤1]。

最后，由于式（14）仅要求[L]的梯度执行更新，从而产生本文提出的OS?加速LALM（OS?ASA）：

[s（k，m+1）=ρM?Lm（x（k，m））+（1-ρ）g（k，m））x（k，m+1）∈prox（ρ-1t）h（x（k，m）-（ρ-1t）s（k，m+1））g（k，m+1）=ρρ+1M?Lm+1（x（k，m+1））+1ρ+1g（k，m）c（k，m+1）=c（k+1）=c（k+1，1）， c∈s，x，g] （18）

当[LM+1=L1]时，类似于典型OS算法，当[M=1]时该算法收敛，当[M>1]时，无法确保收敛。

2 OS?ASA重建CT图像

X?ray CT图像重建问题是一种约束正则化加权最小二乘问题，使用下式置换求解式（14）：

[A←W12Ay←W12yh←R+lΩ] （19）

式中：[lΩ]指凸集[C]的特征函数，因此，式（14）中内部极小化问题转变为一种约束去噪问题。使用快速迭代收缩/阈值算法（FISTA）[9]求解该内部约束去噪问题，以上一次更新作为热启动的开始。

通常，FISTA迭代次数越多，本文算法收敛速度越快，但是，当[n]较大时，迭代内更新的开销不可忽略。在典型X?ray CT图像重建问题中，优化一般较容易，因为统计学加权[W]的巨大动态范围。因此，大部分情况下[t?1]，极大地消除了约束去噪问题中的正则化功率。实践中，约束去噪问题可解决，一至两次迭代后，正则化问题就可满足容差条件。

本文使用具有对角Hessian的SQS函数[DL?diag{A′WA1}]优化加权二次型数据拟合项[8][L]，还使用具有Hessian [DL]的SQS函数优化缩放增广拉格朗日中的二次型惩罚参数，并使用具有位反转顺序的子集梯度增量更新图像。

3 实验结果

针对具有各种几何形状的真实CT扫描的3?D X?ray CT图像，对几种基于OS算法的重建结果进行比较，算法包括：OS?SQS?M，具有[M]个子集的标准OS算法[7]；OS?Nes05?M：基于[M]个子集的Nesterov快速梯度方法，即OS+momentum算法[10]；OS?ASA?M?[ρ]?[n]，利用[M]个子集和[n]次FISTA迭代，提出固定的AL惩罚参数用于求解内部约束去噪问题。

OS?SQS是层析重建的标准迭代方法，OS?Nes05是使用Nesterov技术的X?ray CT图像快速重建方法。不同于其他几种基于OS算法，由于迭代内更新，本文算法具有额外开销，然而，当[n=1]时，即约束去噪问题具有单梯度降落，每次迭代一个向前/向后投影对和[M]个正则化梯度评估，故几种算法具有相同的计算复杂度。

本文从肩区域螺旋CT扫描重建一幅512[×]512[×]109图像，其中，正弦图大小为888[×]32[×]7 146，间隙为0.5，子集最大数约为40。图1为初始FBP图像中央横断平面的裁剪图像、参考重建以及使用本文算法（OS?ASA?40?c?1）在第30次迭代时的重建图像。图（a）为肩扫描从初始FBP图像的中央轴平面裁剪的图像[x（0）]（显示从800～1 200 HU）；图（b）为参考重建x*；图（c）为使用提出的算法（OS?ASA?40?c?1）在第30次迭代处的重建图像[x（30）]。图2表示差分图像，即[x（30）-x*]，针对几种基于OS的算法，标准OS算法（20和40子集）在差分图像中展示了可见条纹伪影和结构化高频噪声。图2中肩扫描：使用基于OS的算法，从[x（30）]-x*的中央轴平面裁剪的差分图像（显示从-30～30 HU）。当M=20时，OS+Momentum算法的差分图像更加结构化和不均匀，但从整体上看，OS+Momentum算法与本文算法的差分图像相似。当M=40时，本文算法的差分图像仍然均匀，而OS+Momentum算法的差分图像带有类似噪声的OS伪影，当M增加时，M=40时，OS+Momentum算法并没有伪影恶化。显然，OS?ASA比OS+Momentum方法具有更好的梯度容差。

图1 参考图像和重建图像

图2 差分图像

针对X?ray CT 图像重建，图3给出了使用FISTA求解内部约束去噪问题。图3中肩扫描为本文提出求解内约束的去噪问题。其中，n取1，2或5以重建图像x（k）和参考图像之间x*的均方差作为迭代的函数。从图3可看出，当使用多次FISTA迭代求解内部约束去噪问题时，收敛速度没有明显提升。因此，只需1次FISTA迭代（n=1），每个子集更新即可满足快速精确的X?ray CT图像重建。

图3 内部约束去噪

4 结 语

增广拉格朗日（AL）方法和有序子集（OS）是两种分别使用分解和近似的加速优化算法，通过考虑AL方法的线性变形，结合这两种技术，本文提出了一种快速OS?加速型拆分算法（OS?ASA），解决了正则化（加权）最小二乘问题。利用OS?ASA对X射线计算机断层扫描（CT）图像进行重建，并将其与几种先进的OS算法进行比较，使用具有不同几何形状的真实CT扫描图像作为数据集。实验结果表明，OS?ASA具有较快的收敛速度和优异梯度容差。今后的研究工作是确定性向下连续化方法的收敛速度分析，以及M大于1时OS?ASA的更严格收敛分析。

参考文献

[1] 赵杰，龚硕然，王龙.脑部CT图像的三维重建及纹理特征研究[J].激光杂志，2014，35（6）：62?65.

[2] 石英，肖金生，刘春晓，等.配气凸轮优化设计的惩罚函数法和增广拉格朗日乘子法[J].武汉理工大学学报（交通科学与工程版），2002，26（3）：365?368.

[3] RAMANI S， FESSLER J A. A splitting?based iterative algorithm for accelerated statistical X?ray CT reconstruction [J]. IEEE transactions on medical imaging， 2012， 31（3）： 677?688.

[4] 林少春.基于有序子集的惩罚加权最小二乘低剂量CT优质重建算法研究[J].中国组织工程研究与临床康复，2008，12（4）：679?682.

[5] 张华军，赵金，罗慧，等.基于梯度投影法与随机优化算法的约束优化方法[J].控制与决策，2014，27（10）：1777?1782.

[6] 何玲君.基于最大似然和罚似然估计的CT统计重建算法研究[D].太原：中北大学，2011.

[7] SCHMIDT M， ROUX N L， BACH F. Convergence rates of inexact proximal?gradient methods for convex optimization [J]. Advances in neural information processing systems， 2011， 24： 1458?1466.

[8] LIN Z， LIU R， SU Z. Linearized alternating direction method with adaptive penalty for low?rank representation [J]. Advances in neural information processing systems， 2011， 24： 612?620.

[9] 孙念，胡炳樑，王爽，等.基于FISTA算法的编码孔径光谱图像压缩与复原系统[J].红外与激光工程，2014，19（1）：238?242.

[10] KIM D， RAMANI S， FESSLER J A. Accelerating X?ray CT ordered subsets image reconstruction with Nesterovs first?order methods [C]// Proceedings of the 12th International Meeting on Fully 3?D Image Reconstruction in Radiology and Nuclear Medicine. [S.l.： s.n.]， 2013： 22?25.