CT图像重建中的一种新型滤波器

杨坪坪，冯汉升，许继伟，宋云涛，

1.中国科学技术大学，安徽合肥230026；2.中国科学院等离子体物理研究所，安徽合肥230031；3.中科离子医学技术装备有限公司，安徽合肥230088

前言

随着计算机断层成像技术（CT）在医学成像、工业无损检测、地球勘探等领域应用愈加广泛，人们对CT 成像质量要求越来越高［1‐2］。重建算法［3‐7］一直是CT技术的核心，其主要分为解析算法和迭代算法，其中解析算法以滤波反投影（Filtered Back Projection,FBP）算法为主，相较迭代算法，解析算法耗时明显更少，因此商业应用更加广泛。就解析算法中的FBP算法而言，其实现主要分为滤波、反投影重建两步，其中滤波器的设计对重建图像质量起着关键作用。常见的滤波器有RL 滤波器和SL 滤波器等［8‐11］，但均各有不足。RL滤波器采用矩形窗对斜坡滤波器进行截取，有明显的Gibbs 现象，表现为严重的振荡效应。SL 滤波器采用sinc 窗函数对斜坡滤波器进行截取，振荡现象减弱，但因其低频段偏离斜坡函数，因此低频段重建质量不如RL 滤波器。此外文献［12］指出，通过加权因子的选择，一次指数滤波器能达到RL 滤波器、SL 滤波器相似的效果，其性能效果是RL滤波器、SL滤波器的折中。本文通过对一次指数滤波器的研究引入高斯滤波器，并进一步设计了一种新型滤波器——改进高斯滤波器，从滤波器设计原理出发，通过仿真实验结果对比分析验证了该滤波器在有无噪声的情况下，都能很好地抑制图像灰度波动，同时能有效改善重建图像质量，重建性能表现优越。

1 原理介绍

1.1 FBP算法

FBP 重建算法可分为二维图像重建与三维图像重建，但其滤波器设计是通用的。FBP算法中的滤波步骤主要是为了消除反投影图像中存在的星状伪影。研究指出星状伪影可以在反投影重建之前进行滤波消除，其滤波实现较先反投影后滤波更为简单，称为FBP 算法［2］。FBP 算法具体实现步骤（以二维图像重建为例）为：

（1）对某一个角度下的投影数据pθ(t)作一维傅里叶变换，记为Pθ(w)：

（2）对Pθ(w)进行一维斜坡滤波器H(w)=|w|滤波，记为P'θ(w)：

（3）对P'θ(w)进行傅里叶逆变换，记为p'θ(t)：

（4）对滤波后的投影图像（0°到180°的投影图像）作反投影累加加权，得到重建图像f(x,y)：

1.2 滤波器设计

理论上滤波器H(w)=|w|是一个无限频带的滤波函数，无法实现［13‐14］。但在实际成像过程中，投影数据高频分量比较小，且投影过程（X射线成像过程）中有平均作用，相当于低通滤波，因此可以采用滤波窗函数对斜坡函数在截止频率（，其中d为探测器间距）范围内进行加窗［15］，如下：

因此，滤波器设计实际上是窗函数选取的过程。判断窗函数的优劣一般有3个指标：主瓣、近旁瓣、远处旁瓣。窗函数的选择一般要求主瓣高而窄，空间分辨率越高，最大旁瓣相对主瓣尽可能幅度小，数值精度越高，密度分辨率越高，同时远旁瓣幅度与幅值小，曲线过渡平滑，以防止严重的Gibbs现象［16‐17］。

1.3 一次指数滤波器

文献［12］指出，一次指数滤波器能通过加权方式达到与RL 滤波器、SL 滤波器相似的效果，其形式如式（6）所示。

其中，a为加权因子，B为截止频率。通过选择不同的加权因子，一次指数滤波器能达到不同的效果。

文献指出，当a取0.64 时，其效果类似于SL 滤波器，其频域曲线［18］、空域曲线［19］如图1和图2所示。

图1 一次指数函数a=0.64，高斯函数k = 1.8与SL、RL函数频域对比Fig.1 First-order exponential function a=0.64,Gaussian function k=1.8,compared with SL and RL functions in the frequency domain

图2 一次指数函数a=0.64，高斯函数k = 1.8与SL、RL函数空域对比Fig.2 First-order exponential function a=0.64,Gaussian function k=1.8,compared with SL and RL functions in the space domain

从图1、图2中可明显看到，当一次指数滤波器权值a=0.64 时，其频域曲线近似直线，与RL 滤波器更相似，其空域曲线存在明显的振荡，频域曲线截止频率附近没有明显压低，因此会存在明显的Gibbs 效应，并且对高频噪声十分敏感。据此，引入高斯滤波器，其形式如下：

式中，k为加权因子。当k=1.8 时，其频域曲线与空域曲线见图1、图2。可以明显看出，高斯滤波器通过加权能更加逼近于SL滤波器，空域曲线光滑，没有出现明显的振荡现象，频域曲线截止频率附近较一次指数滤波器明显压低，能有效抑制Gibbs现象。且当k=0时，高斯滤波器等同于RL 滤波器，因此高斯滤波器实际应用效果是RL、SL 滤波器的折中。为了进一步验证高斯滤波器相比一次指数滤波器的优良性，后文将会采用仿真实验对比验证。

2 改进高斯滤波器

为了进一步改进高斯滤波器，本文设计了新型滤波器——改进高斯滤波器，其形式如下：

式中，a1、a2为加权因子。可以看出，当a1= 0时，改进高斯滤波器退化为高斯滤波器。其中a1加权因子的引入是为了实现对高斯滤波器截止频率附近的进一步压低以减小Gibbs 效应以及噪声影响，其加权因子对比频域曲线见图3。

图3 a2=1.8时，a1取0、0.2、0.3时改进高斯滤波器频域对比图Fig.3 Comparison of the frequency domain of the improved Gaussian filter when a1=0,0.2 and 0.3,with a2=1.8

更进一步，本文做了针对不同加权因子的剖线灰度曲线对比图，其参数与后文的重建参数一致。灰度曲线表征了图像的粗糙程度，曲线波动越大，图像越粗糙。从图4中可以看出，当a1 给定时，a2 越大图像越平滑，能更好地还原原始模型。

图4 a1= 0.5,a2取2.3、4.0时的灰度曲线对比Fig.4 Comparison of grayscale curve when a2=2.3 and 4.0,with a1=0.5

为了进一步验证实验结果评价图像质量，引入常用的医学图像判据，其形式如下：

（1）归一化均方距离d：

式中，tx,y、rx,y分别表示模型原图和重建图对应点的像素密度，t'表示模型原图密度的平均值，图像像素N×N个。归一化均方距离表征了少数点的大误差情况。

（2）归一化平均绝对距离r：

归一化平均绝对距离表征了多数点的小误差情况。

从表1可以看出，a2 越大图像越平滑，但相应的重建误差将会增大，其中归一化均方距离误差明显升高。因此实际重建过程中，需要根据需求灵活选择加权因子。

表1 a1= 0.5,a2取2.3、4.0时改进高斯滤波器的重建误差对比Tab.1 Comparison of reconstruction error of the improved Gaussian filter when a2 = 2.3 and 4.0,with a1=0.5

3 改进高斯滤波器仿真实验对比

为了测试改进高斯滤波器实际应用效果，对滤波器在无噪声情况与加噪声情况下的重建过程进行模拟实验验证。其中重建模型采用二维Shepp‐Logan图像［20］，重建模式选用二维平行束重建，投影角度取0°到179°，步长为1，重建图像大小为512×512，探测器采样点数512，重建像素间隔与探测器采样间隔均取1 mm。分别取无噪声与有噪声情况下的重建结果图，灰度曲线图与重建误差表（分别为归一化均方距离误差表与归一化平均绝对距离误差表）进行对比，其中灰度曲线图采用像素位置256处的重建图像横向剖线，一次指数滤波器权值取0.64，高斯滤波器权值取1.8，改进高斯滤波器权值a1=0.5，a2=2.3。

3.1 无噪声情况下不同滤波器重建结果与相应的灰度曲线图

通过图5与表2可以发现，一次指数滤波器重建效果与RL 滤波器相似，灰度曲线波动很大，图像粗糙，重建误差也很大。高斯滤波器重建效果与SL 滤波器相似。改进高斯滤波器灰度曲线相比其他滤波器明显更平滑，其重建误差明显减小，归一化均方距离与SL 滤波器差距不大，其归一化平均绝对距离却比SL滤波器小许多。因此，在无噪声情况下，改进高斯滤波器重建结果明显优于其他滤波器。

表2 不同滤波器重建误差表Tab.2 Reconstruction errors of different filters

图5 无噪声情况下不同滤波器重建结果与相应的灰度曲线图Fig.5 Reconstruction results of different filters in the case of noise-free and the corresponding grayscale curves

3.2 有噪声情况下不同滤波器重建结果与相应的灰度曲线图

从图6可以看出，在有噪声情况下，一次指数滤波器重建灰度曲线图略优于RL滤波器，且明显劣于SL滤波器、高斯滤波器与改进高斯滤波器，因此重建图像非常粗糙。高斯滤波器重建灰度曲线较SL滤波器更平滑。改进高斯滤波器重建灰度曲线较其他滤波器波动明显更小，更加趋近于原始图像。

图6 有噪声情况下不同滤波器重建结果与相应的灰度曲线图Fig.6 Reconstruction results of different filters in case of noises and the corresponding grayscale curves

为了进一步验证改进高斯滤波器的抗噪性能，采用抗噪性能明显优于RL 滤波器的SL 滤波器进行对比，做了3 组噪声对比实验，加噪噪声采用平均值为0、方差为1 的不同强度高斯噪声进行对比，产生噪声命令如下：

其中，k为参数，I为原始图像，noise为噪声图像，randn能根据图像大小产生标准正态分布序列。因此可以通过选取不同k值达到模拟不同噪声强度的目的。

从表3、表4可以看出，改进高斯滤波器在不同噪声强度下保持重建归一化均方距离误差相差不大的情况下，其归一化平均绝对距离误差明显更小，因此拥有比SL滤波器更优的抗噪性能。

表3 SL滤波器、改进高斯滤波器在不同噪声强度下的重建归一化均方距离Tab.3 Reconstruction normalized mean square distances of SL filter and improved Gaussian filter under different noise intensities

表4 SL滤波器、改进高斯滤波器在不同噪声强度下的重建归一化平均绝对距离Tab.4 Reconstruction normalized average absolute distances of SL filter and improved Gaussian filter under different noise intensities

4 结论

本文从已有对指数滤波器的研究出发，分析了一次指数滤波器的不足，引入高斯滤波器，并进一步改进提出一种新型滤波器。通过仿真实验对比分别验证了该滤波器较RL 滤波器、SL 滤波器、一次指数滤波器、高斯滤波器的优良特性，它有效抑制了重建时的灰度波动，重建图像更加平滑，且在有无噪声的情况下，均能保持更佳的重建图像质量，重建误差明显更小。此外，在改进高斯滤波器加权因子选择a1=0.5、a2=2.3 时，该滤波器能在有效保持图像平滑的同时减小重建误差。