计算机断层扫描(Computed Tomography,CT)^{[1, 2]}是计算机与X线检查技术相结合的产物,它能够得到被检测物体的断层灰度图像并且不损伤原物体,然后根据这些灰度值检测出断层面的内部结构。实际应用中常用滤波反投影算法(FBP)和代数重建算法(ART)这两大类方法来实现CT图像重建。其中滤波反投影算法具有重建速度快,重建质量好的优点,但它的局限在于重建前必须有完备的投影数据,可完备的投影数据在实际操作中却往往不容易得到。ART算法却能很好地克服这一缺点,它是运用迭代的方法解线性方程组,实现了在投影数据较少的情况下同样重建出高质量的图像的目的。虽然ART算法重建的速度较慢,但随着计算机技术的发展,计算速度不再是需要考虑的问题,只需把精力集中在如何提高图像的重建质量上。

本文简单地介绍了ART算法的基本原理,分析了引入松弛因子的原因和松弛因子的有无对重建图像质量的影响,借助仿真实验来研究不同投影数下选择松弛因子的标准,尽可能用最短的时间得到最优的重建图像,为实际重建时选择合适的松弛因子提供理论上的参考依据。

1 ART算法理论

假设被测物体某断层面的线型衰减系数的分布函数为f(x,y)。如果单色射线以I₀的入射强度穿过物体时,探测器上检测到的射线强度用I来表示,那么根据比尔定理这个过程满足

令穿过物体后射线的投影值为p,则p可以表示为p=ln(I₀/I)=∫f(x,y)ds。可见投影值p对应的就是f(x,y)沿每一条射线的线积分。用各个方向不同位置的投影值来重新得到图像的分布函数f(x,y)^{[3, 4]}的过程就是图像重建的过程。

ART算法是一种在离散域^{[5, 6]}中进行运算的级数展开法。首先把待重建图像所在的区域分解成N=n×n个矩形小方格,所有角度下穿过图像的射线总条数为M,待求解的图像分布函数为f(x,y),根据射线衰减的物理过程就可以构建出如下数学方程

ART算法的第一步是先给函数f(x,y)定义一个初始值,然后进行迭代运算,其迭代公式为:

事实上,式(3)只有在理想的情况下才成立,在实际采集投影数据时必然会存在误差,求解时就不再是求解一个方程组而是求解一个不等式组,那么也就不能直接利用式(3)进行迭代了,需要在式(3)中添加一个松弛因子λ。此时的迭代公式就变为:

所以实际上用ART算法重建图像的过程就变成了解一个带有松弛因子的线性方程系统的过程,因此,提高算法的重建速度和重建质量都与如何选择合适的松弛因子有很大的关系,特别是在实际的重建过程中。但是目前只有一些基本的标准来规定如何选择松弛因子,完整规范的理论依据还有待进一步研究。

2 松弛因子对重建图像质量的影响 2.1 松弛因子的引入

一般松弛因子λ∈(0,2),因为只有在这个区间内ART算法才能按照一定的标准保持其收敛性。图 1用几何的方式解释了松弛因子在ART算法中起到的作用。

图 1

Fig. 1

图 1 松弛因子的几何解释 Fig. 1 Geometric interpretation of relaxation factor

假设只有两个变量,图 1中的L₁和L₂表示理想情况下的投影方程,L₃和L₄表示带有误差的投影方程,X₀表示随机给定的一个迭代初始值。由于存在误差的关系,如果按照式(3)进行迭代,迭代过程中X₁和X₃会投影到L₃上,X₂投影到L₄上,迭代至最后得到的是L₃和L₄交点,但这并不是想要得到的结果。如果用加了松弛因子后的式(4)进行迭代,那么迭代过程中X₁和X₃投影到L₁上,X₂投影到L₂上,最终得到的解就是L₁和L₂的交点,这才正是所期望的结果。

一些学者认为可以用实验的方法选择出可靠的松弛因子^[7],首先用实验模拟仿真出类似重建任务的环境,然后选用不同的松弛因子对图像进行多次重建,在最后得到的重建结果中选择出重建效果最佳的松弛因子,就选择这个松弛因子作为最佳值运用到实际的重建任务中。在投影数据含有很大的噪声时,可以直接选择较小的松弛因子^[8],一般λ∈(0.025,0.25)。在这个范围内选择松弛因子进行重建,得到的图像会比较平滑,噪声也相对较小。目前针对这种现象比较合理的解释是:加入小的松弛因子相当于加入了一个低通滤波器,它有效地滤掉了重建图像中含有很多噪声的高频成分,也就是说较小的松弛因子可以有效地抑制迭代过程中的高频噪声。

上述提及到的方法虽然能够达到不错的效果,但该方法却存在很大的局限性,因为不断变化的实际环境在实验室中是不可能模拟出来的,所以需要寻找一种新的选取松弛因子的算法。从ART算法的迭代式(4)可以看出,投影数据值和图像被修正的程度存在着正比例的关系,即投影数据值越大,图像向量被修正的程度也就越大,而且由式(1)可以知道投影值与物体强度衰减的大小也是成正比例关系的。通常选用的常数松弛因子相当于对所有的投影值都采用了同一个低通滤波器,重建图像的边缘看起来比较模糊,因为低通滤波器在抑制高频噪声的同时也抑制了图像本身的高频成分。基于上述所说的正比例关系,文献^[9]提出了一种能够反映投影数据变化特点的动态选择松弛因子的方法。假设λ₀为通常重建图像时选用的常数松弛因子,文献^[9]提出的动态松弛因子的表达式是:

3 实验分析 3.1 常数形式的松弛因子实验分析

实验中选用经典的Shepp-Logan^{[10, 11]}头模型作为重建对象,改变投影角数和松弛因子对其进行迭代,然后分析结果,研究松弛因子λ对重建图像质量的影响。该实验中选取的重建图像大小为128×128。分别选取λ=1,λ=0.2和λ =1.5在投影角数θ为60个,90个和180个进行一次迭代,所得到的重建图像如图 2所示。

图 2

Fig. 2

图 2 不同的松弛因子在不同投影角数下的重建图像 Fig. 2 Reconstructed image of different relaxation factor and the number of projection angle

为了更加客观地评价重建图像和原图像之间存在的误差,实验还计算了重建图像的归一化平均绝对距离判据r。实验中将迭代初值F取为0,松弛因子分别取0.02，0.08，0.2，0.5，1.0和1.5进行实验仿真。ART算法迭代后的重建图像与原始图像的误差如图 3所示。观察图 3可以得出结论:ART重建算法进行迭代时,如果投影数比较多时,可以选取稍小一点的松弛因子;随着投影数的减少,就要逐渐选取稍大的松弛因子,但一般不会超过1。如果松弛因子选取过大,则对图像向量的修正程度就会偏大,重建出的图像与原图像有比较大的偏差。如图中所示当投影数θ为180个时,松弛因子选择0.2最佳。一般在有松弛因子的情况下,仅需要4～6次迭代就可以得到比较满意的重建图像了。

图 3

Fig. 3

图 3 选取不同投影时的重建误差 Fig. 3 Reconstruction error of different projections

从实验结果中可以看出越小的松弛因子重建的图像越平滑,伪影越少。但是加入松弛因子后,重建图像的边缘却趋于模糊了。

用256×256的Shepp-Logan头模型来研究动态松弛因子对图像重建的影响,式(5)中的λ₀通常选择0.2。选用动态因子和常数因子为0.2的图像重建结果如图 4所示。

图 4

Fig. 4

图 4 两种松弛因子的比较图 Fig. 4 The comparison of two kinds of relaxation factor

由图 4的实验结果表明:按照文献^[9]所提出的方法选择动态松弛因子进行重建,与常数的松弛因子进行重建结果相对比,动态松弛因子重建的图像更清晰,边界效果更好。

4 结 论

如上述实验结果所示,在ART算法中如何选择松弛因子会直接影响到图像的重建质量。判断选择的松弛因子是否合适时还需要考虑以下两个因素:投影数据的采集方式和测量环境的噪声类型。在用ART算法进行图像重建时,选择合适的松弛因子可以达到用比较少的迭代次数得到同等质量的图像的目的。如果用动态的方法来选择松弛因子则可以使重建图像的边界效果比较好。今后将进一步实验研究出更多有效的松弛因子的选取方法,这对工业CT上实现不完全投影重建具有重要的意义。