CT 图像重建算法的比较研究

陈宗桂，董晓军，张英俊，管海辰

（湖南医药学院 医学院，湖南怀化 418000）

计算机断层扫描(Computerized Tomography，CT)在肺部疾病和脑出血检查方面优于其他设备[1-2]，是X 射线影像诊断的重要组成部分。CT 扫描过程中X 射线的生物学效应、图像质量和X 射线剂量三者之间的关系，一直是人们关注的问题[3-5]。若患者检查时X 射线的辐射剂量尽可能低，则采集图像信噪比低，分辨率差[6-7]。为此，该研究对比采用滤波反投影、直接反投影和迭代重建算法得到CT 图像，以峰值信噪比、结构相似性和平均梯度评估重建图像的质量，总结出最合适的算法，以辅助临床医师提高病灶检出率。

1 算法介绍

1.1 直接反投影算法

设f(x,y) 为被检测物体截面线性吸收系数的分布函数，一束强度为I0的X 射线从被检测物体透射出后X 射线强度降为I。X 线强度I满足比尔定理，即：

式中，θ是X 线束的投影方向与直角坐标系中X轴的夹角，表示螺旋扫描中X 线束的投影方向；R是单束X 线束在投射方向上的一维函数。令gθ(R)=ln(I0/I)得到式（2）：

gθ(R)是穿过物体后X 线束的投影值。从式（2）可以看出投影值是对沿着投照方向对被检体截面线性吸收系数的分布函数f(x,y)求线积分值，而图像重建是从每个方向上获得不同的投影值，通过反投影重建以求得物体组织密度的分布函数f(x,y)。

CT 图像重构的理论基础[8]是傅里叶投影切片定理，它指出f(x,y)在某一方向上的投影函数gθ(R)的一维傅里叶变换函数gθ(ρ)是f(x,y)的二维傅里叶变换函数F(p,θ)(极坐标形式)在(p,θ)平面上过原点的一条直线，如图1 所示。

图1 中心切片定理示意图

因此，为了从不同投照方向获得足够的数据，并对每一个方向的数据进行傅里叶变换，将变换后的数据填满频域空间。对完整频域数据进行傅里叶逆变换就可得到物体断面组织密度的图像f(x,y)，如式（3）：

式中，F(u,v)表示投影数据的二维频谱图。

1.2 滤波反投影算法

为了抑制CT 图像上的“星芒状伪影”，需要对采集的原始数据进行滤波处理。因此，令u=ρcosθ，v=ρsinθ，则式（3）可进一步变形为：

根据傅里叶变换的性质，频域范围内的图像滤波可以通过空域范围内的图像进行卷积操作来完成，因此式（4）又可以写成式（5）：

式中，h(R) 为滤波函数 |ρ|的空域形式，因而gθ(R)h(R)δ(xcosθ+ysinθ-R)表示对空域图像进行卷积得到初始图像即滤波反投影方法。

在滤波反投影算法中，滤波器的设计是影响图像精度至关重要的因素。理想滤波器是具有无限频带和积分结果不收敛。根据Paley–Wiener 准则，理想滤波器是不存在的[9]。但是，在实际成像过程中，对采集的数据增加一个窗口函数就可以得到理想的数据精度。因为CT 图像分辨率有限，使得探测器采集数据能量分布都集中在低频部分。所以，理想滤波函数 |ρ|必须加窗，换言之，仅保留频域空间上的低频数据。

该研究采用一种比较平滑的窗函数sinc(ρ/2ρ0)rect(ρ/2ρ0)来约束滤波函数 |ρ|，从而得到了S-L 滤波函数。其频域公式为：

对应的空域离散形式为：

S-L 滤波函数重建图像的波动幅值较低，抗噪声性能强和细节分辨率高。S-L 滤波函数对低频段数据响应效果好，而对于高频段数据响应就不如R-L 滤波函数。

1.3 迭代重建算法

迭代重建算法将图像数据视为未知图像矩阵，将在任何投影角处收集的投影数据与图像初始猜测解进行比较，通过比较结果更新模拟图像，直至模拟图像逼近原始图像。算法的原理是首先设置一组模拟图像矩阵，从不同角度采集数据，并与实际采集的投影数据进行比较。通过反向投影比较结果来更新模拟图像，并且可以重复获得重构图像，直到满足迭代结束条件。图2是CT图像迭代重建的系统框图。每次迭代使用一个方程组，即每次迭代都只考虑一个辐射的影响和贡献，通过一定次数的迭代逐次近似所需图像。

图2 迭代重建算法的迭代过程

迭代重建(Algebraic Reconstruction Technique，ART)算法先设定模拟图像矩阵投影值qi。当投影角度为θ时，计算模拟图像矩阵的投影值qi和实际投影值pi。通过误差Δp=pi-qi对投影角度θ下的所有像素点进行校正。迭代重建的公式为:

2 CT图像的重建结果

在该研究中，采用三种不同的算法重构CT图像，如图3 所示。图3（b）是采用同步迭代算法重建的图像。在实际应用中，由于图像尺寸大，为了求解图像上每一个像素值需要建立多个方程组，计算量大，复杂程度高，所以提出了迭代重建方法。图3（b）是采用式（8）进行迭代重建CT 图像的。该算法重建图像时对噪声不敏感，它与其他代数迭代算法不同的地方在于每次迭代不止对一条投影方向上的像素进行修正，还修正了与该直线平行的所有投影方向上的投影数据，这是该算法能有效抑制图像重建过程中的噪声的主要原因。随着迭代次数的无限增加，该算法得到的图像无限逼近原始图像。但是，算法收敛速度慢，在重构图像的边缘和细节存在不同程度的模糊。图3（c）是采用直接反投影算法得到的图像。由于在反投影不断叠加的过程中，中间低频信号不断叠加，而周围高频信号不断缺失，使得图像中感兴趣区域的边缘出现了“星芒状伪影”。这种情况在投影数据少的时候更明显。图3（d）是采用滤波反投影算法得到的图像。探测器采集的是0～180°方向的投影信号，滤波反投影算法是分别对每一个方向的投影数据进行一维傅里叶变换后填满频谱图。所有投影信号在频域中滤波，即乘以加权因子r。通过傅里叶逆变换将滤波后的投影数据恢复到时域，得到原始图像。每一次滤波后投影信号被反投影且最后重叠。该算法采用的滤波器是S-L滤波器。S-L 滤波器不是在频域中加窗截断噪声和干扰信号，而是通过一些相对平滑的窗口函数对空域图像卷积。因此，S-L 滤波后图像振荡幅值更小，图像的结构特征显示更清晰。

图3 三种不同算法重建得到颅脑CT图像

3 图像质量的评价

该研究采用三种算法对CT 图像进行重建，如果想知道原始图像的质量是否有所提高，则需要采用一系列的定量指标进行评价。图像信息越详细，图像边缘灰度变化越快，图像越清晰，人类视觉分辨率越高，结构特征越明显。该研究采用峰值信噪比、结构相似性、平均梯度对重构后的图像边缘结构特征、清晰度和图像相似性进行评估。

3.1 图像重建时间

为了提高X 射线的检查效率，缩短设备的成像时间，改进图像重建算法的运行速度就显得尤其重要。为此，选择的五组图片的尺寸分别是128 pixel×128 pixel、256 pixel×256 pixel、512 pixel×512 pixel、1 024 pixel×1 024 pixel 和2 048 pixel×2 048 pixel。采用不同的算法重建相同尺寸的图像，比较这三种算法的运行时间，如图4 所示。通过图4 可知，针对相同尺寸的图像，直接反投影和滤波反投影的重建时间基本一致。但是，迭代重建算法的运行时间明显大于另外两种算法。另外，图像矩阵从128 pixel×128 pixel 增至2 048 pixel×2 048 pixel，三种重建算法的运行时间都变长。

图4 三种算法重建CT图像所需的运行时间

3.2 基于峰值信噪比的图像质量评价

随着现代医疗影像设备的迅速发展，其无论是在影像成像方面还是在图像后处理方面，都有着广泛的应用前景，新的成像模式和图像处理技术一直在提高和完善。因此，医学图像质量评价已成为制定治疗方案、评估外科手术、预后疗效评估的重要参考。高质图像能如实反映人体内部的解剖结构，为影像诊断和治疗提供丰富的信息。在精准影像服务中，图像质量在保证图像准确度方面起着重要作用。目前，评价图像质量的方法有主观评价方法和客观评价法[10-13]。主观评价方法以观察者为主体，对被测图像的优缺点进行评估。客观评价法是采用相应算法计算图像的一些指标，通过这些指标来评价图像优劣。该研究采用均方误差（Mean Squared Error，MSE）和峰值信噪比（Peak Signal to Noise Rate，PSNR）评价重构图像[14]。均方误差法先计算初始图像和重建图像的差值，再除以行与列的乘积，如式（9）。均方值反映了该算法重建得到图像与初始图像的逼近程度。该方法的优点是简单、快速，可以将误差量化表达且具有一定的参考价值。根据式（10）可求出不同算法下重构图像的峰值信噪比。若重建后图像的PSNR 值越高，就代表该图的微细结构特征反映好，失真少。该研究通过比较峰值信噪比来判断三种算法重建得到的CT 图像的质量，如图5 所示。针对五组不同大小的图像，每一种算法重建得到图像的PSNR 变化基本一致。但是，和另外两种算法相比，滤波反投影算法的峰值信噪比最好，失真最小。直接反向投影算法重建图像的峰值信噪比最小，微细结构特征反映差。

图5 比较三种算法重建得到的CT图像的峰值信噪比

3.3 基于结构相似性的图像质量评价

结构相似度(SSIM)是对亮度相似性、对比度相似性以及结构相似性这三个参量的联合度量，是评价两幅图像特征结构相似性的有效度量。通过比较这三部分的信息，从而获得两幅图像的整体度量[15]。最后采用求均值的方法计算得重建前后两幅图像的结构相似信息。其表达式为：

式中，μx、μy分别表示图像X 和Y 的均值，σx、σy分别表示图像X 和Y 的标准差，分别表示图像X 和Y 的方差。常数C1、C2用来防止分母可能等于零时，分式出现不稳定现象。该研究采用SSIM比较三种不同算法重建得到图像的结构相似性以评价算法的优缺点，如图6 所示。从图6 可知，针对相同尺寸的图像，滤波反投影重建得到图像最接近原始图像。而直接反投影未对采集数据进行预处理，使得重建图像存在“星芒状伪影”与原始图像偏离程度最大。

图6 比较三种算法重建得CT图像的结构相似性

3.4 基于平均梯度的图像质量评价

平均梯度（Average Gradient）是指图像边界点附近的灰度变化差异，可用来表示图像清晰度[16]。平均梯度反映图像中感兴趣区域的细节差异与纹理特征变化[17-18]。一般来说图像平均梯度越大，表达信息越全面，图像清晰度越高。因此，其值越大越好。图像平均梯度计算公式如下：

式中，M×N表示图像矩阵，ΔIx与ΔIy分別表示采用水平方向与垂直方向的结构元素对图像形态学处理后得到的梯度信息。通过平均梯度评价三种不同算法重建后的CT 图像的质量，如图7 所示。从图7 可知，滤波反投影和迭代重建算法在不同图像矩阵中的平均梯度最好，这说明滤波反投影法重建得到图像能够较好反映微小细节反差与边缘特征变化，而直接反投影算法得到的CT 图像的平均梯度最低，说明该算法对于边缘特征变化的反应较差。

图7 比较三种算法重建得CT图像的平均梯度

4 结束语

该研究采用了三种不同算法重建CT 图像，并比较重建图像的质量差异。从运行时间分析可知，对于相同尺寸的图像，滤波反投影算法的重建时间短，重建过程中图像损失的细节和失真都较少，且接近原始图像。直接反投影算法和滤波反投影算法的运行时间相近，但是前者重建的图像效果较差。迭代重建算法与其他两种算法比较，耗时最长。从算法的复杂程度和收敛程度进行分析可知，三种算法都收敛。但是，迭代重建算法为了使重建图像收敛至初始图像，迭代次数必须接近无穷。这必然会增加计算机工作的复杂性和运行时间。从抗噪性分析可知，直接反投影算法的峰值信噪比较低，说明该算法重建得到的图像鲁棒性较差。虽然迭代重建算法的峰值信噪比比直接反投影算法高但是其收敛慢，图像边缘和细节存在不同程度的模糊。从结构相似性和平均梯度分析可知，滤波反投影算法的结构相似性和平均梯度比另外两种算法好，说明该算法重建得到CT 图像的细节分辨最清楚，重建得到CT 图像与原始图像最接近。因而在实际应用中，针对相同尺寸的CT 图像重建，滤波反投影算法的运行时间最短，重建图像失真小，抗噪性强，分辨率高。对于图像重建算法的未来研究重点应该是提高重建算法的抗噪性，缩短运行时间，同时降低算法重建过程的复杂性，加快算法收敛。