随着计算机断层成像(computed tomography, CT)在医疗诊断、工业无损探伤和食品安全检测等领域广泛的应用^[1-5],其对图片质量的要求也越来越高,进而对软硬件特别是算法的要求也越来越高.在众多算法中,FDK算法^[6]应用最为普遍.FDK算法是由Feldkamp、Davis和Kress提出的一种基于圆轨道扫描的近似重建算法,其本质是滤波反投影(filtered back projection,FBP),它对图像质量的影响非常明显.传统的滤波函数能保持较高的空间分辨率,但同时引起明显的Gibbs现象(即有明显的振荡效应),致使重建效果较差.因此,在保持高图像分辨率的情况下设计新的滤波函数,降低Gibbs现象就尤为重要.为此,研究人员相继提出NEW滤波函数^[7]、R-L-S-L、R-L-NEW和R-L-MS-L等混合滤波函数^[8-10]来消除Gibbs现象.虽然这些滤波函数都具有保持高空间分辨率的同时减小Gibbs现象的作用,但无法完全消除Gibbs现象,主要原因是与之结合的R-L滤波函数的近旁瓣突出,远旁瓣的幅度和宽度较大^[9-10],而NEW滤波函数和R-L滤波函数一样,由理想滤波器推导出,所以其近旁瓣也较为突出^[7].本研究基于FDK算法,设计了一种新型滤波函数进行图像重建,在保持较高空间分辨率的情况下大幅减小Gibbs现象,提高了密度分辨率,进而改进了重建图像质量.

FDK算法因其高效、简便以及易于图像处理器(graphics processing unit,GPU)加速,至今仍被大量使用.FDK算法本质上是FBP,把得到的锥束投影数据进行滤波,然后利用扇形束近似重建而得.FDK算法的计算公式^[11]如式(1)和式(2).

对投影数据P的加权滤波为

P'=(R/((R²+a²+b²)^1/2)P(β, a, b))*G(a)(1)

其中, R为光源到物体中心的距离; β为源绕中心旋转轴z轴的旋转角度; a为旋转中心的横坐标; b为旋转中心的纵坐标; P(β, a, b)为投影数据; G(a)为滤波函数; 符号*表示卷积运算.

反投影重建结果为

f(x, y, z)=∫₀^2π(R²)/(U(x, y, β))P'(β, a, b)d β(2)

其中, U(x, y, β)=R+xcos β+ysin β为加权函数,这里x和y为探测器坐标; P'(β, a, b)为加权滤波后的投影数据; 其他变量定义如式(1).

影响重建结果好坏的关键在于滤波函数的选择、重建过程中插值的方式,以及锥角的大小等.从式(2)可直观地看出,滤波函数能直接影响到反投影重建的结果.所以,合适的滤波器是实现高质量重建图像的关键因素之一.

基于FDK算法,分别用NEW、R-L-S-L、R-L-N、R-L-MS-L和NEW-MS-L滤波函数对尺寸为256×256×256像素的Shepp-Logan三维头模型^[16]进行重建,然后观察中心面处的截面及其纵向第128行的灰度曲线图,结果如图4.

图5(a)至图5(f)分别为采用NEW、R-L-S-L、

图4 原始模型图像与用NEW-MS-L滤波函数的重建图像
Fig.4 Original model image and reconstruction image by NEW-MS-L filter function

图5 原始模型及使用不同滤波函数重建结果的灰度曲线
Fig.5 The profile pictures of original model and reconstructed results by different filter functions

R-L-NEW、NEW-MS-L和R-L-MS-L滤波函数重建的纵向第128行的灰度曲线图.灰度曲线图的波动大小反应与原始模型的差异,波动越大说明差异大,表现为图像越粗糙.由图5可见,NEW-MS-L较其他滤波函数能够更好的还原原始模型.

仅通过灰度曲线图的对比来说明NEW-MS-L滤波函数的优越性是不够的,还需通过计算NEW、R-L-S-L、R-L-N、R-L-MS-L和NEW-MS-L滤波函数的归一化均方误差d和归一化平均绝对误差r, 如式(7)和式(8),进行定量对比.

d=[(∑^N_i=1∑^N_j=1(f ^rec_ij-f_ij)²)/(∑^N_i=1∑^N_j=1(f ^rec_ij-f ^-)²)]^1/2(7)

r=(∑^N_i=1∑^N_j=1|f ^rec_ij-f_ij|)/(∑^N_i=1∑^N_j=1f_ij)(8)

其中, i和j分别为图像的横坐标和纵坐标; N为横纵坐标的点数; f_ij^rec为重建结果的灰度值; f_ij为原始模型重建结果的灰度值; f ^-为原始模型重建结果的灰度值的平均值.

d反映少数点的大误差, d=0表示重建后图像完全再现了原始模型图像; d越小表示两者的偏差越小. r反映多数点的小误差, r=0表示重建图像与头模型原始图像无误差; r越小说明误差越小.

表1为分别采用NEW、R-L-S-L、R-L-NEW、R-L-MS-L和NEW-MS-L滤波函数重建图像后的d值和r值对比.

表1 重建结果的误差值对比
Table 1 Reconstruction error value

由表1可见,采用NEW-MS-L滤波函数后的d值与采用NEW滤波函数后的d值相当,且所得r值较其他滤波函数的都小,说明采用NEW-MS-L滤波函数所得的重建效果更佳.

上述结果还可通过调整k₁和k₂值来进一步优化^[10]. d和r随k₁的变化如表2,由表2可知,当k₁=0.8时, d值最小; 但当k₁=0.7、 0.8和0.9时, d值虽然相差无几,但r值随k₁的增加而变大,且变化较明显.所以,本研究权衡d和r的值,选取k₁=0.7.

表2 误差值d与r随k1的变化
Table 2 Variation of error value d and r with k1

当k₁=0.7, k₂=0.3时,可得到d=0.315 5, r=0.737 3. 虽然r值增大了,但d值进一步减小.将NEW-MS-L滤波函数的d值和r值与表1的其他滤波函数的d值和r值进行对比,结果反映新滤波函数性能比其他滤波函数要好.

基于混合滤波和加权平均理论,构建新型滤波函数NEW-MS-L.其空域主瓣分布图具有主瓣高而窄、近旁瓣小、远旁瓣衰减快等特征,能很好地抑制ibbs现象,保持较高的空间分辨率.仿真结果表明,通过对比NEW、R-L-S-L、R-L-NEW、R-L-MS-L滤波函数重建结果与灰度曲线图,并分析归一化均方误差和归一化平均绝对误差,发现NEW-MS-L滤波函数可在保持较高空间分辨率的情况下更有效地抑制Gibbs现象,密度分辨率更高,重建结果也更平滑.此外,根据实际情况适当调整k₁和k₂, 能使图像效果达到符合预期的要求.

NEW-MS-L滤波函数在本次仿真实验中主要针对传统基于吸收效应的CT,今后可进一步探讨NEW-MS-L滤波函数在相衬CT^[17-22]中的应用.