基于Wiener滤波的反投影图像重建算法

马 敏 赵 亮 姬晶晶

（中国民航大学航空自动化学院 天津 300300）

计算机层析成像(CT)是通过对物体进行不同角度的射线投影测量而获取物体横截面信息的成像技术。CT技术的核心是投影重建图像，其实质是由扫描所得的投影数据反求出成像平面上每个点的衰减系数值。图像重建方法大致分为[1]：变换法和级数展开法。滤波反投影算法是目前应用最广泛的基于变换法的图像重建算法[2,3]，具有重建速度快、空间和密度分辨率高等优点，缺点是对投影数据的完备性要求较高，否则会产生伪影。本文分析了滤波反投影重建算法产生伪影的原因，提出一种新的图像重建方法，将 传 统滤波反投影算法与Wiener滤波相结合。实验结果表明，该方法较好地克服了图像伪影现象，图像质量更好。

1 重建图像伪影

1.1 伪影概述[4]

CT图像伪影又称伪像，是指图像重建过程中，不同类型的图像互相干扰和其他各种非随机干扰在图像上的表现，使部分影像根本不反映被测物体的对应区域。CT投影是整个CT过程的输入，投影的好坏直接影响重建图像的质量。若投影数据先天不足，则无论何种算法都不能得到完美的重建图像，得到的只是带有伪影的图像，且投影数据截断越大，伪影越多。

造成投影质量差的原因主要为[5]几何误差、探测器缺陷、射线源不理想以及各种噪声。几何误差因扫描方式而异，对二维CT图像中心旋转中心不重合，会造成“拖尾”伪影；探测器缺陷和探测器阵列元素响应不一致，会造成环形伪影；射线源的多色性，会导致射线穿透试件时的硬化效应，从而造成杯状伪影；射线源能量歧离也会产生条纹状伪影；噪声主要包括随机噪声、散射噪声等。在射线穿过物体时，由于像素对射线采样的不足，通道信号会产生混叠失真，出现伪影。而滤波反投影重建的本质是把取自有限物体空间的射线投影均匀地反投影到射线所及的无限空间的各点之上，因此产生的是星状伪影。

1.2 伪影原理

本文以反投影重建算法为例来阐述伪影产生机理[6,7]，该算法可等效为以原图像为输入，重建后图像为输出的成像系统(图1)。

图1 反投影重建的等效成像系统Fig.1 Equivalent imaging system of back projection reconstruction.

该系统的点扩展函数(PSF)可如下求得。设位于坐标原点的点源δ(x,y)为x-y断面中唯一像点，并设扫描方式为平移/旋转。待重建图像为f(x,y)，射线方向由旋转坐标轴xr与固定坐标轴x的夹角确定，射线透射物体后的投影为pφ(xr)。

扫描坐标系统见图2，虚线代表射线投影方向，其方向始终与yr平行，当φ角从0°旋转到90°时，得到一系列投影pφi(xr)xr=rcos(θ–φ)。

图2 平移/旋转扫描方式的坐标系统Fig.2 Coordinate system of tran s lation/rotation scanning mode.

设输入图像为点源，利用 反 投 影重建算法，重建图像即为该系统的点扩展函数[5]。

式中，φ0为α(φ)=0的唯一解，令α(φ)=rcos(θ–φ)，则α(φ0)=rcos(θ–φ0)=0。

因r≠0，故sin(θ–φ0)=1。所以

可见，相应于反投影重建算法的系统，它的点扩展函数不是δ-函数，虽在r=0处能反映原图是点源的情况，但在r≠0处，像素值不为 0 ，这是此种算法存在伪影的本质。

2 Wiener滤波概述

Wiener滤波属于反卷积算法，是从噪声中提取信号的一种滤波方法，由Wiener[8]提出，在一、二维信号处理和图像复原领域得到广泛应用。在滤波反投影算法(Filtered Back Projection, FBP)对图像进行重建时，滤波函数的选择对消除反投影所产生的星形伪影及图像重建质量影响很大，目前常用R-L滤波函数和S-L滤波函数。R-L滤波函数优点是形式简单、实用，重建图像轮廓清晰，缺点为重建图像有明显的振荡现象，且当投影数据含有噪音时，重建质量比较差。S-L滤波函数优点是对投影数据的高频成分有抑制作用，重建图像的震荡响应小，缺点是图像分辨率不高。相比而言，Wiener滤波法频域形式简单，计算效率高，复原效果良好，且抗噪性能优良，在重建图像的精度和速度方面都有明显的优势。

以一线性系统为例来阐述Wiener滤波原理[9]，若它的单位样本响应为h(n)，Wiener滤波器的输入-输出关系如图3所示。

图3 Wiener滤波原理示意图Fig.3 Schematics of the principle of multistage Wiener filter.

由图3，x(n)为随机输入信号，且

其中，s(n)表示信号，v(n)表示噪声，则输出y(n)为

y(n)为s(n)的估计值，用 ŝ(n)表示，即

用e(n)表示真实值与估计值之间的误差，即

显然，e(n)是一随机变量，可能为正或负。因此，选择用它的均方误差来表达误差，所谓均方误差最小即它的平方的统计期望最小：

3 Wiener−Hopf方程的求解[10]

按最小均方误差准则确定 Wiener滤波器的冲激响应h(n)，令ξ(n)对h(j)的导数等于零，即得式中，Rxs(m)=E(x(n)s(n+m))是s(n)与x(n)的互相关函数，Rxx(m) =E(x(n)x(n+m))是x(n)的自相关函数。式(8)所列为 W iener滤波器的标准方程或Wiener-Hopf方程。若已知Rxs(m)和Rxx(m)，则解此方程即可求出Wiener滤波器的冲激响应。

设滤波器冲激响应序列的长度为N[11]，冲激响应矢量为

滤波器输入数据矢量为

则滤波器的输出为

Wiener-Hopf 方程也可写成

其中

式(13)为s(n)与x(n)的互相关函数，是N维列矢量；R=E﴾x(n)xT(n)﴿ 是x(n)的自相关函数，是N阶方阵。利用求逆矩阵的方法直接求解式( 1 2 )，得

式中，opt表示“最佳”，即FIR Wiener滤波器的冲激响应。

利用Wiener滤波器以平滑图像边缘，它对于图像变化越明显的区域所起的作用越强。实验结果显示，在程序中加入Wiener滤波后得到图像更清晰。

4 仿真及实验结果

4.1 数值仿真

为验证该法的有效性，进行了数值仿真实验。本文重建图像分辨率为 560×420。对三种气液两相流流型进行仿真。程序在 Intel(R) Core(TM)2 Duo CPU、1 G内存的PC上运行。操作系统为Windows XP系统，程序运行平台为MAT LAB7.0。由投影数据分别对泡状流、核心流和层流采用不同图像重建算法进行图像重建，其中滤波反投影重建算法采用的是R-L滤波函数，两种方法迭代次数均为20。仿真重建图像如图4所示。

图4 FBP与改进的FBP算法重建的图像Fig.4 Images reconstructed using the FBP and improved FBP algorithm.

4.2 实验结果

所用多相流CT实验装置的射线源为1.11×1010Bq(300 mCi)241Am源，γ射线能量为59.5 keV，探测器阵列由17个CdZnTe探测器组成，单个探测器的尺寸为3 mm×7 mm×3 mm。将直径为29、22.4、13.4 mm的有机玻璃棒插入充满水的流体管道中，模拟三个泡状体的泡状流；将Φ60 mm有机玻璃管套在Φ82 mm流体管内，将水充入内管道模拟核心流；用Φ82 mm有机玻璃管模拟流体管道，在该管道中充入一半水模拟层流。从7个角度采集实验数据，在 M atlab7.0环境下分别用滤波反投影重建算法和加Wiener滤波的反投影重建算法作图像重建。

为比较改进前后算法的仿真结果，用空泡份额及相对误差对重建图像的质量进行评估。空泡份额是指气液两相流中，气相所占截面积与总流通截面积之比。以流体管道中含有三个气泡为例计算空泡份额，并与实验模型中的实际比例进行比较。根据已知实验模拟的流体管道管径和三个模拟泡状气体管径数据，可算出三个模拟泡状体管体的截面占流体管道截面的比例为Pa=0.1251、Pb=0.0746、Pc=0.0267。用图像处理技术算出的各模拟泡状管截面所成像所占流体管道截面成像的比例用P'a、P'b、P'c表示，空泡份额及相对误差比较结果见表1。

滤波反投影重建算法与基于 Wiener滤波的反投影重建算法均采用20次迭代，对泡状流、核心流和层流进行图像重建，重建时间见表2。

表1 空泡份额及相对误差比较Table 1 Comp a r ison of void fra c tio n and relative e rrors.

表2 FBP与改uter time of the FBP算法的重建 oved FBP and impr algo进的 时间Table 2 Comp FBP ns).(20次迭代)rithm (20 iteratio

从实验结果对比可见，本文所用改进的重建算法的重建效果明显好于传统的滤波反投影的重建效果。在近似相等的时间内，滤波反投影重建图像虽然也可以重建图像，但表面有许多细条纹和亮纹，有明显的 伪 影，而基于维纳滤波的反投影重建算法重建的图像则清晰得多，表层均匀，无杂质。由表2和表3，加入维纳滤波后，重建图像明显伪影减少，更接近原始图像，因此有一定的可行性和应用性。

5 结语

本文提出了一种应用于CT成像系统图像重建算法的改进方案，针对滤波反投影算法产生伪影的问题，通过采取将Wiener滤波技术融入传统算法的方法提高重建图像的质量。此方法对有泡状流、核心流和层流等典型流型的模拟管道进行图像重建，实验结果表明，方法较好地克服了伪影，成像的边界更加清晰，获得较好的图像重建效果。

1 Gilbert B K , Welch B, Holmes D R. Rapid execution of fan beam image reconstruction algorithms using efficient computational techniques and special-purpose[J]. IEEE Trans Biomedical Engineering, 1981, 28(2): 98–115

2 刘杰, 施寅, 阮秋琦. CT快速图像重建算法研究[J]. 中国医学物理学杂志, 2003, 3: 15–17 LIU Jie, SHI Yin, RUAN Qiuqi. The study of fast image reconstruction[J]. Chinese Journal of medical Physics,2003, 3: 15–17

3 Lucas G P, Albusaidi K. An axial scanning system for investigating vertical gas-liquid flows in the slug and bubbly-slug transition regimes[J]. Flow Measurement and Instrumentation, 1998, 9: 9–35

4 Luenberger D G. Linear and nonlinear programming[J].MA: Addison-Lesley, 1989, 5: 19–25

5 庄天戈. CT原理与算法[M]. 上海: 上海交通大学出版社, 1992: 34–35,77 ZHUANG Tiange. Principles andAlgorithms of CT[M].Shanghai: Shanghai Jiaotong UniversityPress, 1992:34–35,77

6 Fossa M. Design and performance of a conductance probe for measuring the liquid fractionin two-phase gas-liquid flows[J]. Flow Measurement and Instrumentation, 1998, 9:103–109

7 Thorn R, Hamme E A, Johanen G A. Recent developments inhere-phase flow measurement[J].Measurement Science Technology, 1997, 8(7): 691–701

8 Klaus Muelle. Anti-aliased three-dimentional cone-beam reconstruction of low-constrast objects with algebraic methods[J]. IEEE Transactiol Imagining, 1999, 5(2):116–134

9Alexander D P, Zayed M R. Adaptive filtering primer with MATLAB[J]. Chinese Journal of medical Physics,2006, 10: 50–75

10 Johansen G A, Froystein T, Hjertaker B T. A dual sensor flow imaging topographic system[J]. Measurement Sci Tech, 1996, 7(2): 297–307

11 Jiang Hsieh. Computer tomography principle, design,artifacts and recent advances[J]. SPIE of America, 2003,5(10): 10–12