引 言

图像重建算法主要分为解析法和迭代法两大类。解析法中传统的滤波反投影(FBP)^{[ 1 ]}算法仍然被广泛使用,它可以在投影数据采集完备时重建出较为准确的图像,但当投影数据采集不完备时,重建出的图像会出现严重的伪影。这时一般采用迭代法来解决,迭代重建法的经典方法是Gorden等提出的ART(algebraic reconstruction technique)算法^{[ 2 ]},该算法可以在投影数据比较少的情况下重建出精确的图像,其原理是将重建问题转换成用迭代方法求解线性方程组的问题。ART算法的特点是:适合于不完备投影数据的图像重建,重建质量好,图像的密度分辨率高,同时对空间分辨率也有很好的保证;其次,重建算法计算简单,不同形式的采样数据重建都适用;另外还可以结合一些已知的先验知识得到更准确的解。由于在ART算法的迭代过程中,每次投影都要计算修正值,但这些修正值并不是完全相同的,穿过同一个像素网格 时,图像的误差修正将会给重建区域带来严重的噪声,且要想得到较好的重建效果,算法需要进行多次迭代,从而导致计算量大,重建时间长,因而人们一直关注的是如何提高该算法的重建效率。ART算法中花费时间最长的是加权因子的计算,本文首先研究了ART算法的3种加权因子计算模型,然后通过理论和仿真实验对用3种模型来重建图像的速度和质量进行比较分析,得出了实用的结论。

1 原 理 1.1 Radon变换

在一个平面内对函数f(x,y)沿不同的直线做线积分,就得到了函数f(x,y)的Radon变换^{[ 3 ]}。也就是说,平面上每个点的像函数值都有一个原始函数的线积分 值与之对应,如图 1所示。

图中每条射线的方程用xcosθ+ysinθ=t来表示,夹角为θ的射线的投影值可以表示为

图像重建的任务就是由各个方向不同位置的投影值P(t)恢复出原图像的分布函数f(x,y)。

1.2 ART算法原理

重建图像的方法有很多种,本文着重介绍和研究ART算法。将待重建的区域离散成一个N=n×n的正方形小网格^{[ 4 ]},其中每个小网格内的值f_j是指所在的第j个像素的平均值。根据图像成像的物理过程和其相应的数学模型,可以用一个代数方程组来表示重建的图像向量和投影数据之间的关系,这样图像重建问题就变成了求解下列线性方程组的问题:

其中,M为所有的射线总数,f_j为上述所讲的投影值,加权因子w_ij表示的是第j个像素对第i条射线的贡献值。如果把p_i看做一个超平面,那么f必须为平面上的一点,如果至少能够得到相同数量的投影值和图像像素,那么所有的超平面的交点就是最终的唯一解。在实际情况中,M、N都是很大的数,一般的矩阵理论很难求出这个方程组的解,所以通常选用迭代的方法来求解。ART算法的迭代公式为

其中λ是松弛因子,一般取值都很小,在0～2之间。由式(3)可以看出,迭代过程中要反复用到加权系数矩阵 w _i,所以在利用式(3)时计算、存储和访问 w _i 是所要面临的最大困难。通常像素的分布位置和射线的几何结构是已知的,可以利用这些条件预先计算出所有的加权系数矩阵 w _i ,然后将它们全部保存在硬盘中,等到进行迭代的时候再调出使用;但是这样频繁的硬盘访问会严重的影响重建速度,所以最有效的办法是在迭代过程中实时地计算更新 w _i ^{[ 5,6 ]}。由此可见用ART算法迭代时加权因子 w _i 的计算直接关系到重建的速度和质量。以下讨论3种计算加权因子的不同模型对图像重建的影响^{[ 7 ]}。

2 重建模型分析

假设在ART重建模型中,每条射线的宽度均为d,相邻射线之间的间距为l ,第j个网格内的像素值大小用f_j来表示。以下对加权因子的三种不同模型^{[ 8 ]}进行分析。

2.1 模型一

如图 2所示,模型一中射线的宽度d与像素网格同宽,而相邻射线之间的距离l=0,并且将像素值定义在每个像素网格的中心点处。判断加权因子w_ij值的依据是每个小像素网格的中心点在哪条射线中,即如果第j个像素的中心点包含在第i条射线内,则加权因子w_ij定义为1,否则就定义为0。

2.2 模型二

如图 3所示,模型二中射线的宽度d=0,而相邻射线之间的距离l与像素网格同宽。用第i条射线经过第j个像素时,与网格相交的线段长度来定义加权因子w_ij,图 3中的加权因子w_ij=l_AB 。

2.3 模型三

如图 4所示,模型三对射线的定义与模型一相同,射线的宽度d与像素网格同宽,而相邻射线之间的距离l=0。用第i条射线和第j个像素相重叠的面积占整个像素小网格的面积的比例来定义加权因子w_ij,图 4中的加权因子w_ij=S_ABCD /d²。

3 实验结果

为了比较上述所讲到的三种模型的不同,选用了128×128的经典Sheep-Logan头模型来进行重建仿真实验。图 5是180个投影角度时分别用3种加权因子模型进行一次迭代后的重建图像,然后取投影角总数分别为60,90,180时进行迭代分析。还选用了归一化平均绝对距离r和归一化均方距离判据q来进行误差比较,用来更精确地反映重建质量。两个判据的具体定义为:

式中,x_i,j,y_i,j分别是指原图像和重建后图像中第i行、j列的像素值;x — 则是指原图像中所有像素值的平均值。其中q可以较敏感的反映出某几点产生较大误差的情况,而r则可以较敏感的反映出许多点均有的一些小误差的情况。这两个判据q和r的值越小,表示重建图像与原图像的误差越小,质量也就越好。

表 1是三种模型分别在60、90、180个投影角度下所用的重建时间,表 2和表 3则列出了用3种加权因子模型重建下的两个判据值。由上述结果可以看出模型三重建的效果最好,但是由于模型三在计算时要计算每条射线与每个像素之间相交的面积,而大部分面积却是不规则的,所以计算加权因子所需要的时间就比较长,也就是说模型三这一计算加权因子的方法没有重建速度上的优势。对于模型一,其在计算加权因子时将加权因子简化成只有1和0两个值,大大提高了重建的速度,缩短了迭代过程所需要的时间。但是这种简化加权因子的方法往往会引起相邻方程之间的矛盾,给重建结果加入了噪声,并且当某一步迭代出现误差后,误差就会直接传递到后面的迭代过程中,使得这种模型重建出的图像质量比较差,常常会出现椒盐(salt and pepper)现象 。而模型二中计算加权因子时仅仅计算了每条射线穿过每个像素的长度,计算量小并且能够通过一定的算法快速实现,另外模型二中射线束所穿过的像素数目一般比模型三中射线束所覆盖的像素数目少,重建时间也就比模型三大大的缩短了,但是重建质量却可以基本与模型三相当。

针对模型二用第i条射线经过第j个像素的长度来定义加权因子w_ij来说,最基本的计算方法是siddon算法^{[ 9 ]},但是siddon算法的计算时间比较长,效率比较低。很多学者对模型二的计算方法已经做了大量的研究,目前国内已经提出的加权因子的计算方法有:张顺利,张定华等提出的快速网格遍历算法^{[ 5 ]};刘郁纪,李公平等提出的siddon改进算法^{[ 8 ]};杨文良,魏东波等提出的改进的逐点遍历法^{[ 10 ]}等。ART迭代算法的速度在不断地提高,其重建时间上的瓶颈也在不断被打破。

4 结 论

本文对ART算法的三种加权因子计算模型进行了研究、比较和分析,得出以下结论:模型一重建时间比较短,但重建质量却很差,并且伴随着严重的椒盐噪声; 模型三重建质量最好,但花费的重建时间却相对比较长;模型二可以在重建时间大大缩短的情况下使得图像的重建质量与模型三相当。因此,要想同时保持较好的重建质量和较快的重建速度,ART迭代过程中加权因子计算模型的最佳选择就是模型二。需要指出的是,除了加权因子模型本身外,还有其他一些因素严重影响着ART算法的重建速度和重建质量,例如投影数据的访问方式^{[ 11 ]}、松弛因子^{[ 12 ]}的选择以及先验知识等。随着计算机技术的不断提高,ART算法的重建速度和重建质量都将会有更大程度上的提高。