总变差约束的数据分离最小图像重建模型及其Chambolle-Pock求解算法∗

乔志伟

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

(2018年4月27日收到；2018年7月2日收到修改稿)

1 引 言

在计算机断层成像(computed tomography,CT)图像重建中,以滤波反投影(filtered backprojection,FBP)法为代表的解析重建方法在商用CT中占据主流[1].其缺点是需要完备的投影数据.如稀疏视角下的FBP重建算法会引起条状伪影[2].随着压缩感知等稀疏优化理论的提出,基于优化的迭代算法成为了十余年来的研究热点[3,4].其中,总变差最小(total variation minimization,TV)算法是最成功的优化算法之一,其已经在平板CT[5]、偏置扫描平板CT[6]、短扫描平板CT[7],C-arm CT[8]、低剂量平板CT[9]及机载成像仪(on-board imager,OBI)平板CT等[10]三维CT中得到了成功的应用,展示了其高精度的稀疏重建能力.

TV算法可以高精度地稀疏重建的原因可以从多个方面来解释.TV算法体现了压缩感知的思想.压缩感知认为,如果物体在某个变换域是稀疏的,那么可以仅仅使用稀疏采集到的原始数据来高精度地重建物体.其实现策略是在数据保真的约束下,使得物体的稀疏变换的ℓ1范数最小[11].TV算法使用的稀疏变换是梯度大小变换(gradient magnitude transform,GMT)[12].TV算法体现了先验约束的思想.在稀疏重建时,迭代法构建的线性方程组是欠定的,使得该线性方程组有无穷多个解.TV算法可以被认为加入了一种先验知识,即物体的GMT是稀疏的.这样,TV算法是从无穷多个解中选取TV最小的解.TV算法也体现了低通滤波的作用.TV算法最初被当作去噪方法应用到图像处理中,其可以在去噪的同时保持图像的边缘信息[13].FBP的稀疏重建引入的条状伪影是一种高频噪声.TV正是通过对这种噪声的低通滤波,消除了条状伪影,实现了高精度稀疏重建.

TV最优化模型有约束形式和非约束形式两种类型.非约束形式的模型参数是一个平衡因子,用来平衡数据保真力度和正则力度,其没有明确的物理意义,难以选择.约束形式的TV模型,其模型参数有物理意义,更容易选择[9].常用的约束的TV模型,以数据保真项为约束,TV正则项为目标函数,其缺点是数据保真项对应的数据容差限这一参数对重建质量敏感,选取方法鲁棒性较差[12].芝加哥大学Pan研究组从2014年陆续提出了TV约束的、数据分离最小(TV constrained,data divergence minimization,TVcDM)新型TV模型,并将之应用到了扇形束ROI重建[14],C-arm CT[8]、短扫描平板CT[7]及正电子发射断层扫描(positron emission computed tomography,PET)[15]中.该模型的模型参数是TV限,其在较大范围的变化对图像质量影响较小,故易于选取,选取方法鲁棒性较好.

TVcDM模型是一种约束的最优化模型,传统的非约束优化模型求解算法,如交替方向乘子法(alternating direction method of multipliers,ADMM)和分裂布莱格曼法(split Bregman,SB)[16],不适宜用来求解该模型.Chambolle-Pock(CP)算法框架[17−20]是一种优秀的最优化算法框架,不但适用于非约束的优化模型,而且适用于约束的优化模型.它有诸多优点:其一,算法参数均可以计算得到而不需人为设定；其二,可以解所有凸最优化模型,无论模型平滑还是非平滑；其三,算法的每个子最优化问题一般均有闭合解.CP算法框架可以用来推导求解TVcDM模型的CP算法实例.

Pan研究组[20]于2012年推导了一系列CP算法实例.这些实例分别用来求解最小二乘模型、非约束的ℓ2-TV模型、非约束的ℓ1-TV模型、非约束的KL-TV模型以及数据分离约束的、TV最小(data divergence constrained,TV minimization,DDcTV)模型,但是当时尚未提出TVcDM模型,所以并没有给出TVcDM模型的CP算法实例的推导.2014—2017年,该组的一系列使用TVcDM模型的CT和PET图像重建的文章,主要聚焦在完全重建、有限角度或者投影被截断情况下的实际重建问题.2018年,本文作者与Pan研究组及芝加哥大学EPRI中心研究人员将这一模型应用于电子顺磁共振成像(electron paramagnetic resonance imaging,EPRI)的稀疏重建中[21].这些研究聚焦在模型和算法在实际数据中的评估,而没有使用仿真数据,即存在金标准的情况下,从定量的角度展开收敛规律、稀疏重建能力评估和参数选择等理论和方法层面的研究.

鉴于如上原因,本文在总结出CP算法实例推导方法的基础上,详细推导了TVcDM-CP算法实例,并给出算法的详细解释；以Shepp-Logan和FORBILD仿真模体[22]的CT重建为例,实现该算法,验证模型和算法的正确性；通过分析多个度量标准的迭代行为,揭示其收敛规律；通过稀疏重建,评估该模型的高精度稀疏重建能力；通过对含噪投影数据的重建,研究了模型参数对重建质量的影响.最后研究了算法参数对收敛速率的影响.

2 方 法

连续-连续(continuous-to-continuous,CC)成像模型把投影和图像均看成连续函数.而离散-离散(discrete-to-discrete,DD)成像模型把投影和图像均看成离散函数.这样,DD模型就是一个线性方程组,基于DD模型的重建就是求方程组的解.一般情况下,该线性方程组是一个大规模的、病态的且往往欠定的(如稀疏重建时)方程组,使得求解成为需要重点攻克的问题.为了得到所需要的有用解,往往需要将此成像模型进一步建模为一个最优化问题,使其可以将稀疏先验等先验知识引入模型.其后,需要设计求解算法,以解最优化模型.

基于如上脉络,以CT重建为背景,以TVcDM模型为研究对象,在该部分将按照成像系统建模-最优化模型-CP算法-参数选择这一研究链展开叙述.因为本文旨在探索新颖的TV模型及其求解算法,不失一般性,以平行束CT作为研究对象.别的成像模态或者CT成像的其他扫描模式遵循同样的规律,其研究方法只需要做相应的类推或扩展.

2.1 成像模型

平行束CT的DD成像模型是一个线性方程组,如(1)式所示:

这里,g是一个大小为M的列向量,表示离散的投影数据；u是一个大小为N的列向量,表示离散图像；A是一个大小为M×N的矩阵,这里代表2D Radon变换；Am,n表示第n个像素对第m个投影测量值的贡献.对2D Radon变换来说,Am,n是第n个像素(正方形)与第m条射线(直线)的相交长度.

设图像的大小为nx×ny,将其一列一列串行化,就可以得到图像的向量表示,其大小为N=nx×ny.设投影集大小为np×na,即共有na个角度下的投影,每个投影上有np个点,将其一列一列串行化,就可以得到投影数据的向量表示,其大小为M=np×na.

实际上,求取Am,n的方法就是所谓的投影方法,因为投影是系统矩阵与图像向量的乘积.投影方法包括像素驱动法[23],Siddon射线驱动法[24],Joseph射线驱动法[25]和距离驱动法[26]等.传统的像素驱动法会引入高频震荡伪影,已经设计了一种精确的像素驱动法[23]；本文使用了精确的像素驱动法来生成投影,即求取系统矩阵.

2.2 最优化模型

由于(1)式所示的线性方程组的病态、大规模及欠定等因素,需要进一步将之建模为一个最优化模型.根据重建的需要,可以将压缩感知、低秩矩阵等理论及各种先验知识耦合,以设计最优化模型.本文研究一种新颖的基于压缩感知的TV模型——TVcDM.

TVcDM,即TV约束的、数据分离最小模型,可以表示为

这里,u?是最优化模型的解,即被重建的图像；∥.∥表示ℓ2范数的平方；图像的TV范数∥u∥TV是图像梯度大小变换|Du|mag的ℓ1范数:

这里D是一个大小为2N ×N的矩阵

D1和D2均为N×N的矩阵,分别表示沿着x和y轴方向的离散梯度变换.

令ux,y,x∈[1,nx],y∈[1,ny]表示二维图像的每个像素,则D1变换可以表示为

D2变换可以表示为

令|.|mag表示对一个二维向量求模,可以是ℓ1范数模,也可以是ℓ2范数模.在本文中,采用ℓ2范数模,即

这样,(3)式的各向同性形式可以表示为

为了提高收敛速度,引入λ和ν两个算法参数以调整TV范数和数据保真项这两个凸函数的相对大小,如此得到的新的TVcDM模型为[7,8,14]

这两个参数虽然出现在最优化模型中,但是这里称其为算法参数,因为它们不能决定最优化模型的解,但是可以影响解的收敛路径和速度.

2.3 Chambolle-Pock算法及其解释

表1给出了用于求解最优化模型(9)的CP算法伪代码.

在表1中,∥.∥sv是矩阵的最大奇异值,其求法见文献[20]中的算法8；向量u,u¯及s的大小为N；向量p的大小为M；向量q和a的大小为2N；6.3中的向量1I是一个大小为N 的“1”向量；6.2中的投影操作符ProjectOntoℓ1Ballνt1的求法见表2的伪代码；uconv是达到设定的实用收敛条件后的收敛解.

在表2中,x表示一个长度为N的向量；a是ℓ1球的半径；m是一个长度为N的向量,其每个元素对应x中相应元素的绝对值；sign(x)是符号函数:正数的函数值是1,0的函数值是0,负数的函数值是−1.

表1 用于求解最优化模型(9)的CP算法伪代码[7,8,14]Table 1.Pseudocode of the CP algorithm for solving the optimization model shown in(9)[7,8,14].

表2 ProjectOntoℓ1Ball的求解算法伪代码[14]Table 2. Pseudocode of the solving algorithm for ProjectOntoℓ1Ball[14].

2.4 CP算法的推导步骤

本节首先从CP算法框架总结出算法实例的推导方法,然后推导TVcDM模型的CP算法实例.

2.4.1 CP算法框架及其推导方法

CP算法框架可以解形如(10)式所示的最优化模型.

在此模型中,x和y是分别属于向量空间X和Y的多维向量；K是一个矩阵,表示一个线性变换,定义了一个从X到Y的线性映射；G和F分别是关于x和y的多元函数,分别定义了从X和Y到实数的一个映射.CP框架要求这两个函数是凸函数,但不要求平滑.

CP算法框架见表3.

表3 CP算法框架(N步迭代)[17]Table 3.The CP algorithm framework(N steps of iterations)[17].

表3中,∥K∥SV是指矩阵K的最大奇异值；F∗是指凸函数F 的共轭凸函数；proxσ[F∗]和proxτ[G]是最邻近映射操作；T表示转置.

任意一个凸函数H(z)的共轭凸函数定义为

这里,需要注意,max操作符不是求取使得目标函数最大的自变量,而是求取当c为最大化值时目标函数的值.

prox最邻近映射的表达式为

该最邻近映射本质上是一个最优化问题,但是因为z是一个变量,所以其结果是一个关于z的函数.

现在,可以总结出使用CP算法框架推导CP算法实例的步骤为:

1)构造如(10)式所示的最优化问题；

2)求取F∗；

3)求 取proxσ[F∗]和proxτ[G]两 个 最 邻 近映射；

4)代入表3的第4行和第5行形成算法实例.

在如上步骤中,关键问题是基于(11)式的凸共轭函数的求取和基于(12)式的最邻近映射的求取.CP算法的优势之一是每个子最优化问题都有闭合解.这就要求在(12)式中,函数H(Z)要足够简单,以使得可以求得闭合形式的解.现有研究表明,在图像重建中使用的大部分凸函数所对应的最邻近映射均有闭合解[20].

2.4.2 TVcDM-CP算法实例的推导

(9)式所对应的最优化问题,可以通过如下变换形成(10)式的形式.

其中,δℓ1Ball是一个指示函数,表示如果一个n维向量在特定半径的ℓ1球内部,则函数值为0,否则,值为∞,其表达式为:该式表示,如果向量x在半径为a的ℓ1球内部,则函数值为0,否则值为∞.

现在,根据如下所示的一系列变换,将(12)与(10)式对应起来.

这样,

根据CP算法实例的推导步骤,现在的任务是求取两个凸函数F1(p)和F2(q)的共轭凸函数,以及求取FF和G三个函数的prox最邻近映射.

根据(11)式,可得

根据(12)式,可得

将(15),(19)—(21)代入表3,可以得到表1所示的TVcDM模型的CP算法实例.

2.5 重建参数

最优化模型及其算法构成了一套完整的重建方案,其中所有的参数构成了重建参数.与模型相关的参数称为模型参数；与算法相关的参数称为算法参数.

模型参数包括TV限t1、像素大小、投影方法等.在本文中,像素大小和探元大小都归一化为1,投影方法采用精确的像素驱动法,故将在3.3节重点研究TV限对重建的影响.

算法参数包括σ,τ,λ,ν和θ,因为σ,τ和θ参数均可以计算得到,而不需要人为设定,所以不需要关注.λ和ν对算法的收敛速度和路径有影响,故将在3.4节研究其对收敛速率的影响.

3 结 果

在本节中,设计4个研究点:1)重建方法的验证及CP算法的收敛行为分析；2)重建方法的稀疏重建能力评估；3)TV限对重建精度的影响；4)λ和ν的选择对收敛速率的影响.

3.1 逆犯罪研究——重建方法的验证及CP算法的收敛行为

从数学的角度讲,图像重建是一个反问题.对于平行束CT重建,解析法,如滤波反投影算法,是反Radon变换问题；迭代法是线性方程组的求解问题.逆犯罪是一种评估反问题是否正确的验证方法.当仿真数据理想且充分的情况下,如果求得的解的误差在计算机精度范围内可以充分小,则认为逆犯罪成功,表明设计的整套重建方案及其实现——从成像系统建模、最优化问题建模、算法设计到计算机实现,都是正确的.

在本部分中,采用Shepp-Logan和FORBILD模型为仿真模型；仿真生成理想而充分的投影数据；使用第二部分的重建方法进行重建；通过逆犯罪研究验证整套重建方法的正确性.

3.1.1 仿真模体及重建条件

分别采用Shepp-Logan模体和FORBILD模体为仿真模型,图像大小为256×256,每个像素被看成单位1大小的像素.旋转中心定位于128×128.探测器长度为256,每个探元大小为单位1.均匀采集[0,π]范围360个角度下的投影,所以投影数据的大小为256×360.采用精确的像素驱动法定义系统矩阵A,并生成投影数据.如此构建的线性方程组,有256×256个未知数,有256×360个方程,属于过定方程.但是因为方程组是相容的,所以方程组有惟一的精确解.

在CP算法中,设定λ=1,b=1即ν=νA,其余算法参数按照表1中计算得到.TV限设定为仿真模型的TV值.

CP迭代算法的收敛条件定义为RMSE(un,utruth)6 10−4.RMSE的定义见(22)和(23)式.

3.1.2 重建结果的定性及定量分析

重建结果采用定性观察和定量分析的方法.一方面视觉检查图像质量,一方面比较图像中线位置的波形,以进行定性观察.定量分析,则采用均方根误差(root-mean-square-error,RMSE)作为度量标准,如下式所示:

式中x是待评估向量,s是向量真值,n是向量的长度.

对Shepp-Logan重建来说,CP算法在迭代到2273次时,达到了设定的收敛条件.重建结果见图1. 图1(a)是Shepp-Logan仿真模型图像；图1(b)是收敛态的重建图像.比较这两个图像可以发现它们几乎完全一样,用肉眼难以分辨彼此.图1(c)是重建图像和仿真模型图像的中心线波形的比较,从图中可以看出,两条波形几乎完全重合.图像和波形的主观评价表明取得了高精度重建.

对FORBILD重建来说,CP算法在迭代到1712次时,达到了设定的收敛条件.重建结果见图2.图2(a)是FORBILD仿真模型图像；图2(b)是收敛态的重建图像.比较这两个图像可以发现,它们几乎完全一样,用肉眼难以分辨彼此；图2(c)是重建图像和仿真模型图像的中心线波形比较图.从图中可以看出,两条波形几乎完全重合.图像和波形的主观评价表明取得了高精度重建.

两种模体的重建图像的RMSE均达到了10−4,一方面表示算法的收敛条件达到了,另一方面表明逆犯罪成功,因为此反问题被高精度地求解了出来.

图1 Shepp-Logan仿真模体的TVcDM重建结果 (a)Shepp-Logan仿真模体图像；(b)重建图像；(c)模型图像和重建图像左右中心线波形的比较:红色线是真实波形,蓝色线是重建波形,二者几乎完全重合；图像的显示窗口是[0,1]Fig.1.The TVcDM reconstruction images of the Shepp-Logan simulation phantom:(a)The Shepp-Logan phantom image；(b)the reconstructed image；(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

图2 FORBILD仿真模体的TVcDM重建结果 (a)FORBILD仿真模体图像；(b)重建图像；(c)模体图像和重建图像左右中心线波形的比较图:红色线是真实波形,蓝色线是重建波形,二者几乎完全重合；图像的显示窗口是[0,1]Fig.2.The TVcDM reconstruction images of the fORBILD simulation phantom:(a)the fORBILD phantom image；(b)the reconstructed image；(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

3.1.3 CP算法的收敛行为

为了描述收敛行为,引入3个度量标准,通过观察3个度量标准的迭代走势,揭示CP算法的迭代规律.这些度量标准如下式所示:

(23)式描述的是重建的图像与仿真模型图像的RMSE距离；(24)式描述的是重建图像生成的投影与原始投影之间的ℓ2范数距离；(25)式描述重建图像的TV值与TV真值的相对误差.在此逆犯罪研究中,因为模型的投影和图像是完全一致的、数据是充分且精确的,所以这些度量标准在收敛态均可以充分小,表示取得了高精度重建.

图3和图4分别展示了Shepp-Logan和FORBILD模体的CP算法的收敛行为,图3(a)—(c)分别展示了(23)—(25)式所示的度量标准的迭代走势.

从图3和图4的(a)可以看出,图像误差不断下降,达到了设定的收敛条件,而且可以看到其有继续下降的趋势.如果继续推进迭代,图像误差甚至可以达到10−6.灰度图像的灰度只有256个等级,即从0到255,理论上说,当图像误差达到1/256≈3.9‰时,显示器已经无法分辨存在此误差水平的两幅图像.所以,我们设定逆犯罪标志,即收敛条件为图像误差小于等于10−4.

图3和图4的(b)显示的是数据误差,即数据残差的走势.数据残差是TVcDM模型的目标函数.可以看出,数据残差不断下降,而且在当前的停止点有继续向下的趋势.

图3 TVcDM-CP算法重建Shepp-Logan模体时的收敛行为 (a)重建误差M1(n)的迭代走势；(b)数据误差M2(n)的迭代走势；(c)重建图像的TV相对误差M3(n)的迭代走势.Fig.3.The convergence behavior of the TVcDM-CP algorithm for reconstructing the Shepp-Logan phantom:(a)The iteration trend of the reconstruction error M1(n)；(b)the iteration trend of the data error M2(n)；(c)the iteration trend of the relative TV error M3(n).

图4 TVcDM-CP算法重建FORBILD模体时的收敛行为 (a)重建误差M1(n)的迭代走势；(b)数据误差M2(n)的迭代走势；(c)重建图像的TV相对误差M3(n)的迭代走势Fig.4.The convergence behavior of the TVcDM-CP algorithm for reconstructing the fORBILD phantom:(a)The iteration trend of the reconstruction error M1(n)；(b)the iteration trend of the data error M2(n)；(c)the iteration trend of the relative TV error M3(n).

图3和图4的(c)显示了TV值相对于TV限的相对误差,其整体趋势向下,但是在两个重建案例中,当迭代到后期,出现了较大抖动,相对误差变大.这种现象是该算法迭代行为的一个特点,在我们各种各样的重建实践中经常出现,这并不影响算法的收敛性.虽然TV值相对误差出现了向上抖动,但是图像误差和数据误差是整体不断下降的,表明最优化模型的解仍然在向收敛态逼近.其实,可以观察到TV值下降到了一个很小的值,比如10−7之后产生了向上的抖动,其达到一定值之后,又会继续向下.此后,即使再有抖动,其上限也会保持在一个很小的值附近,比如10−4.而TV值是指图像梯度变换的ℓ1范数,TV值在此极小范围的波动,分配到图像每个像素上(在本节案例中,图像包含256×256=65536个像素),图像基本没有变化.

从图3和图4的所有子图可以看出,所有的迭代曲线均存在抖动现象,这是TVcDM-CP算法迭代行为的一个特点.图像误差和数据残差(目标函数)整体会不断下降,直到计算机浮点数精度(本文因为投影与反投影操作使用GPU加速,所以均采用了单精度数)所允许的收敛态.

3.2 重建方法的稀疏重建能力评估

为了评估TVcDM模型的稀疏重建能力,在仿真实验中,分别从30,60,90和120个角度重建图像,分析投影个数对重建的影响.重建条件与3.1节相同,除了Shepp-Logan模体重建的收敛条件变为迭代至1000次结束,而FORBILD模体重建的收敛条件为迭代至2000次结束.

图5是Shepp-Logan模体的TVcDM的重建结果与FBP算法重建结果的比较.由图可以看出随着投影个数从120个减少到30个,FBP重建的条状伪影越来越严重,30个角度下的FBP重建引入了明显的伪影；同时看到,各个角度下的TVcDM重建图像几乎完全相同,看不出明显的因为投影个数减少而产生的质量退化现象.30个角度下的投影已经可以实现TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

图5 Shepp-Logan模体重建中TVcDM与FBP算法的重建结果比较 图像上面的数字表示投影个数,左边的文字表示所使用的算法；显示窗口为[0,1]Fig.5.Comparison of the reconstructed Shepp-Logan images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

图6是FORBILD模体的TVcDM的重建结果与FBP算法重建结果的比较.由图可以看出,随着投影个数从120个减少到30个,FBP重建图像的条状伪影越来越严重,60个角度下的FBP重建引入了明显的伪影而30个角度下的伪影已经异常严重；同时看到,各个角度下的TVcDM重建图像几乎完全相同,看不出明显的因为投影个数减少而产生的质量退化现象.30个角度下的投影已经可以实现TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

为了定量比较TVcDM与FBP算法在不同投影个数情形下的重建精度,图7绘制了重建图像的RMSE相对投影个数的走势图.图7(a)对应Shepp-Logan模体重建；图7(b)对应FORBILD模体重建.两个模体重建对应的RMSE走势规律是相同的.随着投影个数从120减小到30,FBP重建结果的RMSE逐渐增大,表示FBP算法对投影个数的变化是敏感的,其没有高精度稀疏重建能力,当投影个数为30个时,其RMSE过大,图像精度过低,导致重建图像失去有效性.而TVcDM算法对投影个数变化不敏感,投影个数为60,90和120时,RMSE非常接近.投影个数为30时,因为投影视角过于稀疏,重建误差增大.但即使TVcDM算法只使用30个角度下的投影,其重建图像的精度也远远高于FBP算法从120个角度下的投影重建图像的精度.

图6 FORBILD模体重建中TVcDM与FBP算法的重建结果比较 图像上面的数字表示投影个数,左边的文字表示所使用的算法；显示窗口为[0,1]Fig.6.Comparison of the reconstructed FORBILD images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

图7 重建图像RMSE随投影个数变化的走势图 (a)对应Shepp-Logan模体重建；(b)对应FORBILD模体重建Fig.7.Plot of RMSE of the reconstructed image as function of projection number with(a)corresponding to Shepp-Logan phantom reconstruction and(b)FORBILD phantom reconstruction.

图5和图6的定性主观视觉检查和图7的定量误差分析均表明:TVcDM-CP算法有从稀疏投影中高精度重建图像的能力.图7(b)中30个角度对应的重建RMSE比图7(a)中30个角度对应的重建RMSE大.此现象表明:不同的重建对象,因为复杂程度不同,所需的最少投影个数是不同的.FORBILD图像比Shepp-Logan图像复杂,所以要达到特定的重建精度,需要的投影个数更多.

3.3 TV限对重建精度的影响

在仿真实验中,给投影加45.0 dB的高斯白噪声,进行重建.投影个数为30个；TV限分别设定为TV真值的0.6,0.8,1.0,1.2和1.4倍；其余重建条件相同.

图8和图9是Shepp-Logan模体重建结果.图8是不同TV限约束下的对含噪投影的重建图像.图9是不同TV限约束下的对含噪投影的重建图像与仿真模体(真值)的差分绝对值图像.从图8(c)和图9(c)可以看出,当TV限选取正确,即选取TV真值时,可以获得高精度的重建结果,TVcDM模型不但压制了因为稀疏采集(只有30个投影)可能引起的条状伪影,而且压制了因为高斯白噪声污染而产生的图像噪声.从图8和图9的(a)和(b)可以看出,当TV限选取过小时,TV范数的平滑作用过大,使得重建图像产生了低频过平滑伪影.而从图8和图9的(d)和(e)可以看出,当TV限选取过大时,TV范数的平滑作用减小,使得重建图像产生了高频噪声.可见,TV限的选择对重建有重要影响,其决定了最优化模型的收敛解.只有在平滑和细节保持之间,寻找到平衡点,选取出最适合的TV限,才可以高精度地重建图像.

图10(a)是Shepp-Logan模体重建中RMSE随TV限变化的走势图.可以看到TV限过大或者过小,都导致了更大的误差,选取合适的TV限会取得较小误差,实现高精度重建.

图11和图12是FORBILD模体重建结果.结合图10(b),可以看出,过小或者过大的TV限均导致了更大的误差,合适的TV限可以取得较小误差,取得高精度重建.

图8 TV限的不同选择对重建精度的影响 (a)—(e)分别是TV限为TV真值的0.6,0.8,1.0,1.2和1.4倍时对应的重建图像；为了突出显示,采用了伪彩色显示方式；显示窗口是[0,1]Fig.8.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

图9 TV限的不同选择对重建精度的影响 (a)—(e)分别是TV限为TV真值的0.6,0.8,1.0,1.2和1.4倍时对应的重建图像与仿真模体的差分绝对值图像；为了突出显示,采用了伪彩色显示方式；显示窗口是[0,0.2]Fig.9.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

两种模体复杂程度不同,图像特征不同,但是在此研究中,TV限对重建精度的影响规律一致.实验结果和定性分析一致表明,TV限对重建精度有重要影响,只有选择最佳TV限才可以在图像的平滑度和细节保持之间找到最佳平衡点,实现高精度重建.

图10 RMSE随TV限变化的走势图 (a)和(b)分别对应Shepp-Logan重建和FORBILD重建Fig.10.Plot of RMSE as function of TV tolerance with(a)corresponding to Shepp-Logan reconstruction and(b)FORBILD reconstruction.

图11 TV限的不同选择对重建精度的影响 (a)—(e)分别是TV限为TV真值的0.6,0.8,1.0,1.2和1.4倍时对应的重建图像；为了突出显示,采用了伪彩色显示方式；显示窗口是[0,1]Fig.11.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

图12 TV限的不同选择对重建精度的影响 (a)—(e)分别是TV限为TV真值的0.6,0.8,1.0,1.2和1.4倍时对应的重建图像与仿真模体的差分绝对值图像；为了突出显示,采用了伪彩色显示方式；显示窗口是[0,0.2]Fig.12.The impact of the selection of TV tolerance on reconstruction accuracy.(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

3.4 λ和ν的选择对收敛速率的影响

仍然使用3.1.1中设定的重建条件,只是将投影个数变为30.采用Shepp-Logan和FORBILD两个模体分别进行实验.为了分析λ和ν的不同选择对重建的影响,设计5种不同的组合形式,使算法均迭代至1000次,观察RMSE的下降态势,比较其对重建速率的影响.

图13是收敛速率比较图,其中图13(a)对应Shepp-Logan重建；图13(b)对应FORBILD重建.

从图13可以看出,不同的λ和ν对应不同的收敛速率.对于两个模体重建来说,λ=1,ν=0.1νA均为收敛速率最快的一个组合.但是如果按照从慢到快排序,对两个模体重建而言,顺序不完全相同.λ和ν的最佳选择是依赖于特定重建对象的.即使对于同一个重建对象,当投影个数发生变化时,λ和ν的最快组合也可能发生变化.对于一个特定的重建任务,要想找到最快的组合,需要使用不同的组合重建,通过观察迭代收敛行为,选择最快的一个组合.

图13 不同的λ和ν的组合下RMSE在1000次迭代范围内的收敛速率比较 (a)和(b)分别对应Shepp-Logan重建和FORBILD重建Fig.13. Comparison of the convergence rates of RMSE of reconstructions with different λ and ν in the iteration range of 1 to 1000 with(a)and(b)corresponding to Shepp-Logan and FORBILD reconstruction,respectively.

4 讨 论

4.1 模型设计相关问题

为了提高TV模型的精确重建能力,人们提出了各种TV范数的变种:广义的各向异性TV范数,保边TV范数,适用于有限角度重建的各向异性TV范数,高阶TV范数；TpV(p∈[0,1])范数等.类似的最优化模型还包括小波变换最小、基于字典学习的稀疏模型以及非局部均值等.所有这些模型均可以用一个统一优化模型来描述.在统一优化模型具体化时,其一要考虑选取何种正则模型,其二要考虑选取何种数据保真模型,其三则是考虑选择约束形式还是非约束形式.

最优化模型是数学模型,它是对物理模型的一种近似,因此各种各样的最优化模型没有优劣之分,只有在特定重建场合下的适用性问题.对于图像重建的一个特定任务,应用驱动地设计最适合的最优化模型是基于优化的重建方案设计的最重要的工作.

最优化模型决定了最终的解,所以模型参数对重建有重要影响.当前,尽管人们已经系统研究了模型参数对重建的影响,但是模型参数的自动选择问题应该是今后的研究重点.

4.2 关于重建算法

本文所用模型的求解是非常具有挑战性的:1)规模巨大,系统矩阵很大,设计的算法中不能出现矩阵求逆操作；2)往往欠定,方程个数远小于未知量个数；3)病态,因为噪声和数据的不一致性,造成方程组是不相容的；4)TV范数是一个不平滑的凸函数.

人们已经为各种各样的最优化模型设计了求解算法,如ASD-POCS算法[12],ADMM算法以及本文使用的CP算法等.需要注意的是,模型决定解,求解算法仅仅决定解的收敛路径和速度.所以,各种各样的求解算法的重建精度是没有可比性的,惟一可以比较的是其收敛速度.本文采用CP算法的原因是该算法保证收敛、算法参数不需要调节且每个子最优化问题均有闭合解.

5 结 论

本文设计了TVcDM-CP重建方法.将图像重建问题建模为一个TVcDM最优化模型.在该约束的最优化模型中,TV正则项是约束项,数据保真项是目标函数.基于CP算法框架,详细推导了TVcDM模型的CP算法实例；经仿真模型重建,验证了模型及算法的正确性并分析了算法的收敛行为；评估了模型的稀疏重建能力；分析了模型参数的选择对重建的影响.

TVcDM模型相较于DDcTV模型,一个可能的优势是模型参数比较容易选择.DDcTV模型的模型参数是数据容差限[12],其大小决定了重建图像的投影和实际的测量投影的距离.在实际重建中,投影数据包含噪声和不一致性,只有选择合适的数据容差限才可以有效抑制噪声和不一致性.选取过小,则噪声被引入；选取过大则重建图像过于平滑,细节丢失.TVcDM的模型参数是TV限,合理的TV限估计会实现高精度重建,而选取过小则图像过于平滑,选取过大则图像被引入噪声.TV限在小范围内的变化对重建影响不大(参考文献[8]的图7),所以该参数的选择比较鲁棒.

本文研究表明,TVcDM-CP重建方法可以高精度地从稀疏角度下的投影中重建图像；在CP算法的收敛过程中,会出现抖动现象,但是整体呈现下降趋势；TV限对重建质量有重要影响,合理选取该参数可以取得高精度的重建质量；λ和ν的不同取值组合会导致不同的收敛速率.

本文所设计的TVcDM-CP算法是一种通用的图像重建模型,可以被应用到平行束CT、扇形束CT、锥形束CT中,也可以被应用到PET,MRI及EPRI中.依据本文的研究链,即成像系统建模-最优化问题建模-CP算法设计-算法实现及验证-实际重建,各种重建模态均可以实现高精度稀疏重建.在此研究链中,需要做的改变仅仅是将本文的成像模型中的系统矩阵变为相应的成像模态对应的系统矩阵,然后任务驱动地选择模型参数和算法参数.