自适应步长非局部全变分约束迭代图像重建算法

王文杰，乔志伟，牛 蕾，席雅睿

(山西大学 计算机与信息技术学院，太原030006)

0 引言

计算机断层成像(Computed Tomography, CT)在医学和工业生产中有非常重要和广泛的应用。绝大部分的商用CT机使用传统的滤波反投影(Filtered Back Projection, FBP)算法来重建图像已经有将近30年的历史[1]；然而，大剂量的辐射会增加患者患癌的风险，低剂量的扫描，即数据不完备的条件下，FBP算法难以重建出精确的图像。为了解决这个问题，有两种常用思路：一是通过插值方法对不完备的投影数据进行补足；二是通过在迭代过程中对重建图像施加某些约束条件，即对模型增加正则化项的手段来解决。一些基于迭代的算法诸如代数迭代重建(Algebraic Reconstruction Technique, ART)[2]和期望最大化(Expectation Maximization, EM)算法[3]已经被应用于图像重建中[4]。基于优化的迭代算法可以利用压缩感知、低秩矩阵等稀疏优化技术对图像进行精确重建[5]。其中，全变分(Total Variation, TV)最小化约束的迭代重建算法已经在各种CT中广泛的应用。

但是，TV模型的主要缺点是它不能恢复和保留图像中的精细结构和纹理，并倾向于产生阶梯效应(staircase artifacts)。因为传统的TV模型仅考虑图像的局部像素信息，且假设图像在单个对象内是平滑的，只在不同对象的边界处具有不连续的灰度跳跃[6]。

Buades等[7]利用兼顾了图像中更多全局信息的非局部均值(NonLocal Means, NLM)算法来进行图像去噪。该算法通过对一个像素点周围像素的加权平均值来恢复该像素，因此，相较于传统的局部算子能更好地恢复和保留图像的细节。Gilboa等[8]将传统的变分框架和非局部算子进行了有机的结合，提出了非局部变分(NonLocal TV, NLTV)框架。相较于局部变分框架,非局部变分框架在基于优化的CT重建算法中能得到更加精确的结果是可以合理预期的。

Sidky等[9]提出了自适应最速下降投影到凸集(Adaptive Steepest Descent Projection Onto Convex Sets, ASD-POCS)算法来解决图像重建中基于约束的TV最小化问题。其中，自适应最速下降(ASD)即算法采用梯度下降法来使图像的全变分最小化，同时步长因子是自适应的。投影到凸集(POCS)即为ART与正约束的合称，ART保证重建结果满足投影数据保真项的约束，正约束是因为默认图像不存在负值像素点。ASD-POCS算法可以有效和鲁棒地从不完备数据中重建出相当精确的结果。它交替地运行POCS过程和TV下降过程，从而使结果逐步地接近TV最小化解，同时满足数据保真项和正则项的要求。

然而，ASD-POCS算法使用梯度下降算法来求解TV最小化，直接使用梯度下降法解决NLTV最小化的约束问题是非常复杂和不实用的。Goldstein等[10]提出了分离Bregman (Split Bregman, SB)算法，可以有效求解很多的L1正则最优化问题。Zhang等[11]把这个方法推广到了非局部最优化问题中。Lou等[12]将非局部算子应用到了图像的反卷积和重建中。实际上，ASD-POCS算法中的TV下降过程可以看作是对POCS步骤后得到的图像的去噪过程，因此，本文把ASD-POCS算法作为基础的重建框架，以期利用其参数自适应的特点，同时，结合SB算法求解L1正则化问题的高效性来解决模型中NLTV最小化的问题，进而提出自适应非局部全变分最小化投影到凸集(Adaptive Non-local TV Minimization Projection Onto Convex Sets, ANLTVM-POCS)算法。

1 相关工作

1.1 CT模型

非约束的、基于优化的CT数学模型可以表示为如下的形式：

(1)

其中：向量g是大小为M的离散投影数据；大小为N的列向量u表示离散图像；A是一个大小为M*N的系统矩阵，这里表示平行束CT的二维Radon变换。系统矩阵A可以通过射线驱动法(Siddon算法[13])和距离驱动算法等来求取[14]。J(u)是正则项，λ是正则化参数。向量u*是最优化模型的解，即被重建图像。

那么，核心问题便是如何选择J(u)，在基于TV的重建型中，J(u)是梯度图像的L1范数：

J(u)=‖u‖TV=

(2)

其中us,t表示图像单个像素点。显然，传统的TV模型仅仅考虑了一个像素毗邻的另外两个像素，相较于非局部算法框架来说，难以保留图像的细节和纹理。

1.2 非局部变分框架

接下来，简要回顾在文献[8]中介绍的非局部算子。令Ω⊂R2，x∈Ω，并且u(x):Ω→R是实函数。设w:Ω×Ω→R是关于参考图像的非负对称实函数，即：w(x,y)=w(y,x)。关于一像素点的非局部导数可以定义为：

(3)

则关于像素点x的非局部梯度算子▽NLu(x):Ω×Ω→Ω定义为其所有偏导数▽NLu(x,·)组成的向量：

(4)

对于向量u在点x∈Ω的非局部梯度L2范数：

(5)

所有像素点的非局部梯度L2范数也可以组成一幅图像，即加权非局部梯度图像。

对于向量v1=v1(x,y)，v2=v2(x,y)，向量内积定义为：

(6)

根据向量内积的定义及散度和梯度的共轭关系，

〈▽NLu,v〉=-〈u,divNLv〉

(7)

非局部散度divNLv(x):Ω×Ω→Ω定义为：

(8)

拉普拉斯算子可以写为：

(9)

权函数w(x,y)：

(10)

其中：Ga是标准差为a的高斯核函数，参数h是滤波参数。式(10)显示了权函数只有在关于点x和点y的比较像素块相似度大时具有较大的取值。

为了便于计算机进行离散运算，分别给出非局部梯度、非局部梯度L2范数、非局部散度和拉普拉斯算子的离散形式：

(11)

(12)

(13)

(14)

综上所述，非局部全变分可以看作加权非局部梯度图像的L1范数：

(15)

其中：wi, j是w(x,y)的离散形式，i和j是图像的像素索引。

对于非局部框架，其核心问题是对三个参数的选取：比较块的尺寸(patch size)、搜索窗的尺寸(search window size)和式(10)中滤波参数h的取值；然而，这三个参数之间的耦合影响非常严重，作出一个恰当的选择是一项比较艰巨的任务。本文在大量实验的基础上，尝试了几种参数组合，然后给出一些可以得到相对精确结果的经验参数选择。

理论上，搜索窗的尺寸应该是整幅图像，然而实际上，一个过大的搜索窗会导致计算开销的几何倍增加。许多文献给出了从11×11到21×21之间的不同选择，一个好的搜索窗尺寸应该是在兼顾精确结果和计算开销的中间选择。在本次工作中，发现11×11大小的搜索窗已经足够得出相对精确的结果且仅需要相对较少的计算开销。

在一些关于NLTV图像去噪的研究[15]中提到，大的比较块尺寸可以增强算法对于噪声的鲁棒性，较小的比较块尺寸可以保留更多的低对比度图像细节。在式(10)中，可以发现，结果像素点的取值还和邻域像素点的高斯均值相关。对于不同的成像模型以及不同的噪声类型，应该根据情况作出恰当的选择。在本文的实验中，大小为3×3的比较块已经可以得到较好的视觉结果图像。

滤波参数h可以说在非局部全变分模型中扮演了灵魂角色。该参数控制着比较块的相似度大小，而相似度的准确性又是非局部框架能得到精确结果的基本保证。尽管如此，很少有文献提到在图像重建中关于参数h应该如何选择。在许多关于非局部框架图像去噪的文献中参数h一般取图像噪声方差σ大小的若干倍。在文献[12]中，对预处理图像进行小波变换；然后使用小波系数的绝对中位差来确定参数h的大小。许多图像的噪声估计算法可以在这里用来对参数h的选择作出参考，本文不再详细展开。

1.3 分离Bregman算法

在引言中曾提到，ASD-POCS算法交替运行POCS步骤和TV最小化步骤，而TV最小化可以看作是对POCS步骤后得到的图像的去噪过程，因此，可以把TV最小化表述为如下最优化模型：

(16)

其中u0是POCS步骤后得到的图像。

标准ASD-POCS算法使用梯度下降法来最小化TV。对于NLTV，使用分离Bregman算法来求解这个最优化问题。引入一个辅助变量d:Ω×Ω→R来替换式(16)中的▽NLu，然后可以得出如下的约束最优化问题：

(17)

s.t.d=▽NLu

式(17)的约束问题可以转化为如下无约束形式：

(18)

通过引入分离Bregman算法框架，式(18)可以转化为如下Bregman迭代形式：

(20)

式(19)可以分离为分别对u和d的递归迭代：

(21)

(22)

利用Euler-Lagrange公式，可以得出式(21)的最优条件：

λ(u-u0)-γdivNL(dk-▽NLu-bk)=0

(23)

根据式(23)，使用Gauss-Seidel迭代算法来求取u的近似解(给出关于一个像素点在第k+1次迭代的离散形式，同d和b)：

(24)

关于子问题d，可以通过收缩算子(shrinkage operators)[10]来求取：

(25)

最后，b可以通过直接迭代求解：

(26)

2 本文方法

本文方法主要受到文献[9,11-12]的启发。2008年Sidky等[9]提出的ASD-POCS算法，应用正则化方法将目标函数定义为优化准则及惩罚项二者之和的形式，以罚函数的形式增强解的稳定性，以罚函数的形式增强了解的稳定性[16]。2008年的方法改进自他们2006年提出的重建算法，2006的重建算法首次利用医学图像梯度变换图像具有稀疏性这一先验知识，结合凸集投影，无论是缺失角度投影数据(e.g. 0～45°)，还是大间隔角度投影数据(e.g. 0～180°间隔6°采样)，相较于传统的FBP解析类重建算法，ART和EM等迭代类算法都表现出了极其精确高效的稀疏重建能力。2008年的ASD-POCS算法考虑了投影数据中存在多种不一致性的问题，新增了数据容差控制因子，大幅度提高了算法的鲁棒性，改善了2006算法重建轮廓附近存在灰度渐变伪影的缺点，并且该算法在迭代过程中步长因子自适应的特点，也改善了其他迭代方法容易陷入局部最优解的问题，能在较少的迭代次数下得到更加精确的解。

本文方法主要基于ASD-POCS算法，并且针对ASD-POCS算法中TV模型的缺点，改用NLTV模型来约束最终解，提高算法保留重建结果中精细结构和纹理的特性。

且针对NLTV传统迭代方法(如梯度下降法)难以求解的问题，利用分离Bregman方法求解L1正则化问题的优势来求取NLTV最小化。

下面给出提出算法的伪代码。

1)

Initialization:β=1.0;βred=0.995;α=0.2;αred=0.95;n=2;rmax=0.95;u=0;;

2)

Repeat main loop

3)

up=u;

4)

fori=1 toNd

5)

u=u+βAi(g-Ai·u)/(Ai·Ai)

6)

end for

7)

fori=1 toNi

8)

ifui<0 thenui=0

9)

end for

10)

dd=‖g-g0‖2

11)

dp=‖u-up‖2

12)

if (first iteration)

13)

dtvg=α*dp

14)

end if

15)

up=u;

16)

Initialization:u0=u;d0=0;b0=0;λ=1;γ=dtvg;

17)

fork=1 ton

18)

updateuk+1by equation (24)

19)

updatedk+1by equation (25)

20)

updatebk+1by equation (26)

21)

end for

22)

dg=‖u-up‖2

23)

ifdg>rmax*dpanddd>

24)

dtvg=dtvg*αred

25)

end if

26)

β=β*βred

27)

Until (stopping criteria)

行4)的Nd是投影的总个数，行7)的Ni是总像素个数。算法中有6个参数(β;βred;α;αred;n;rmax)控制整个算法的迭代过程，在文献[9]中有详细的叙述，本文不再展开。u是算法的初始化图像，在这里设定为零。行3)和15)的up是临时变量，用来存储算法当前步骤的结果并参与步长自适应调整的计算。dtvg是算法最小化子问题17)到21)行的步长因子。算法的另一个重要参数是数据容差，是多种数据不一致的简化，包括简化模型带来的误差、投影数据中的噪声以及X射线的散射等，这使得不可能总是重建出与数据完全一致的结果图像[9]。参数确定了重建生成图像的投影与实际投影数据的欧氏距离(L2范数)。在结果部分讨论了关于的选择对于重建结果的影响。

把标准ASD-POCS算法的TV梯度下降过程替换为基于NLTV的分离Bregman迭代过程。这里有两个重要参数(λ,γ)在分离Bregman的迭代过程中占据主导性地位，进而影响最终重建结果的精度。在下一章详细展开讨论。

3 实验结果与讨论

为了评估和检验本文的算法，设计了两组实验。实验采用大小为128×128的Sheep-Logan和Forbild头部仿真模体(http://www.imp.uni-erlangen.de/forbild/english/results/index.htm)，成像采样旋转中心为第64×64号像素点。探测器个数为128，单个探元大小为1。采样角度范围为0到180°，所有采样均匀间隔。系统矩阵采用Siddon算法生成。图像第一组实验是通过理想的不含噪声的欠采样数据来进行运算。使用两组不同投影角度个数(20和30)生成的投影数据来检验本文算法的稀疏重建能力。第二组实验，在模拟的投影数据中加入高斯噪声(强度通过高斯噪声方差衡量)来评估算法的抗噪能力。如图1所示。

本文通过重建图像剖线比较的方式来评估算法的重建特性。通过均方根误差(Root Mean Squared Error, RMSE)来定量地评估重建结果的精度：

其中：u和f是分别是重建结果图像和真值图像；N是图像的总像素个数。

图1 仿真模体Fig. 1 Simulation phantoms

3.1 不含噪声数据重建

至于算法的停止条件(stopping criteria)，可以当计算结果达到一个特定的条件(如‖uk-uk+1‖/‖uk‖<10-3)时结束循环，但是为了更直观地描述本文的算法，在本组实验中使算法直接循环500次，以便更全面地通过算法迭代过程的RMSE曲线来分析算法的收敛特性。

参数γ的恰当选择能使分离Bregman算法的子问题u和d收敛的速度更快。在标准分离Bregman算法中，为了使算法更好地收敛，γ既不能过大也不能过小。文献[10]中提到，当γ=2λ时，分离Bregman算法通常可以取得较好的收敛结果。在本文的算法中，参数γ可以说是连接数据保真项和正则化项之间的桥梁。本文直接设定γ等于标准ASD-POCS算法框架中的自适应步长参数dtvg，以使ASD-POCS算法步长参数自适应的特点与分离Bregman迭代过程有机地结合起来。

对于参数λ，它影响了算法在保留更多原始投影数据内容和去除噪声之间的平衡。如果λ过大，数据保真项的权重会增加，换句话说，结果可能将会过拟合给定的投影数据，相应的数据中的噪声也会更多地被保留；相反，较小的λ将会增大正则化项的权重，也就是说，算法不仅会去除更多的噪声信息，同时也会平滑模糊掉图像中更多的细节。一个好的选择应该是在去噪和尽可能保留图像细节之间的某种折中。因为参数γ在本文的算法中是自适应的，所以直接设定λ为1。在文献[17]中，对这两个参数(γ和λ)的选择进行了更加详细的讨论。

关于滤波参数h的选择，经过大量实验后，对于Sheep-Logan模体和Forbild头部模体，分别取0.02和0.03可以得出较为精确的结果。

图2是ASD-POCS算法和本文算法在20和30个投影角度下的Sheep-Logan模体的重建结果图像。图2的实验结果表明：传统TV算法重建的结果图像中含有更多的噪点，并且图像细节也被平滑掉；相对来说，NLTV重建的结果图像保留了更多的细节，图像结构边缘清晰。

为了更清晰地定量比较图2四幅重建结果的图像，在图3中给出了四幅图像的垂直于图像中心的切面垂直中心剖线。对于20个角度的重建条件，NLTV的垂直中心剖线拥有更少的波动。在30个投影角度的情况下，NLTV的垂直中心剖线与真值图像的垂直中心剖线基本重合。结果表明，相对于原TV算法，NLTV算法可以获得更高的精度。

图2 Sheep-Logan模体不含噪声投影数据重建图像Fig. 2 Reconstruction results of Sheep-Logan phantom by noise-free projection data

图3 Sheep-Logan模体不含噪投影数据重建结果图像垂直中心剖线Fig. 3 Vertical profile of Sheep-Logan phantom reconstructed by noise-free projection data

图4给出了算法迭代过程的RMSE趋势曲线。显然，NLTV模型可以在更少的迭代次数下重建出更好的结果图像，且在收敛点处拥有更高的精度。在20个角度投影数据重建的条件下，TV模型最终的收敛均方误差是0.011，NLTV模型的最终收敛精度是0.003；在30个角度的条件下，TV模型最终精度是0.002，NLTV模型则得到了十万分之一精度(5.3×10-5)。

图4 Sheep-Logan模体不含噪重建RMSE收敛趋势曲线Fig. 4 RMSE convergence curve of reconstructing Sheep-Logan phantom by noise-free projection data

图5显示了30个和40个投影角度下Forbild头部模体的重建结果。Forbild模体的右侧有蜂窝状结构，在较小的区域有密集的灰度变化，对于重建算法的要求较高。可以看到，因为TV模型对于单一结构内灰度均匀这一先验条件的假设，TV模型对于蜂窝状结构的重建结果较差，有明显的放射状伪影，且TV模型在增加投影角度数量后，对于蜂窝状结构的重建结果并没有明显的改善；而NLTV模型对于蜂窝状结构的恢复更加精确。

图5 Forbild头部模体不含噪声投影数据重建结果Fig. 5 Reconstruction results of Forbild head phantom by noise-free projection data

表1显示了6组实验的均方根误差结果。在Sheep-Logan模体的仿真实验中，TV模型在50个投影角度下的重建精度与NLTV模型20个角度下的精度在一个数量级。

表1 6组实验结果RMSE对比 Tab. 1 RMSE comparison of 6 sets of experimental results

3.2 含噪声数据重建

本节比较了TV和NLTV模型在含噪声Sheep-Logan模体投影数据的重建情况。首先，生成了50个角度的投影数据；然后，在投影数据中加入了方差为0.01的高斯白噪声。

重建结果如图6所示。可以看出，TV模型和NLTV模型都显示出了对于噪声的抑制特性，但是，NLTV模型图像结构之间有着更清晰的边缘间隔，且整幅图像看起来比TV模型更加锐利清晰。

图6 Sheep-Logan模体含噪声投影数据重建结果Fig. 6 Construction results of Sheep-Logan phantom by noise-added projection data

图7显示的垂直中心剖线结果表明，NLTV模型比TV模型拥有更高的精度，且在像素灰度值跳跃较大的地方更好地与真值重合。这也进一步证实了重建结果图像中显示的，NLTV模型在重建图像的不同结构之间有着更锐利的边缘区隔。

图7 Sheep-Logan模体含噪重建结果图像垂直中心剖线Fig. 7 Vertical profile of Sheep-Logan phantom reconstructed by noise-added projection data

图8显示的含噪声重建过程的RMSE收敛曲线进一步表明了，NLTV相比传统TV模型能在更少的迭代次数下达到更高的收敛精度。NLTV模型在迭代到30次时已经基本接近最终的收敛精度0.002 2，而TV模型则需要迭代到150次时，才基本接近最终的收敛精度0.005 5。

图8 Sheep-Logan 模体含噪重建RMSE 收敛趋势曲线Fig. 8 RMSE convergence curve of reconstructing Sheep-Logan phantom by noise-added projection data

3.3 参数对重建结果的影响

在这组实验中，并不计划去寻找能使重建结果最优的参数组合，仅仅通过选定测试参数的不同值来讨论其对于重建结果的影响。

(28)

图9 数据容差对重建结果的影响Fig. 9 Effect of data tolerance on reconstruction result

3.3.2 参数λ

接下来，简单讨论平衡参数λ对于重建结果的影响。为了使对比更加明显，把噪声方差增加到0.03。

图10的结果显示，平衡参数λ对于重建结果的影响与数据容差有某种相似性：较大的λ可以增加数据保真项的权重，但会造成重建图像中含有较多的噪声，这是因为没有很好地抑制在原始投影数据中所含的噪声内容；相反地，较小的λ可以对原始投影数据中的噪声内容进行很好的过滤，但是会使结果丢失细节。

图10 平衡参数λ对重建结果的影响Fig. 10 Effect of balance parameter λ on reconstruction result

4 结语

本文以ASD-POCS算法为基础，结合非局部全变分模型的优点，采用分离Bregman算法求解算法中的L1约束问题，提出了ANLTVM-POCS算法。该算法结合了ASD-POCS算法参数自适应的优点和非局部全变分模型能更好保存图像精细结构及纹理的特点，在一定程度上改善了ASD-POCS算法采用TV模型不能很好恢复图像细节的弱点，同时该算法也可以简单地扩展到诸如三维束重建和螺旋重建中。实验结果表明，该算法可以在较少的重建次数下达到较高的重建精度，且相比ASD-POCS算法达到收敛点所需迭代次数更少。由于非局部框架的引入，对于平衡参数λ和滤波参数h的选择，本文只是在大量实验结果的基础上给出的经验参数。合适的参数是能重建出精确结果的保证，能根据模型的特点开发出更加自动精确的参数选择算法是未来的工作重点之一。