一种用于CT偏置扫描重建的指数型加权函数
杜乙, 王贤刚, 向新程
清华大学 核能与新能源技术研究院, 北京 100084
通讯作者:王贤刚,副研究员,wangxiangang@mail.tsinghua.edu.cn

作者简介:杜乙(1987—),男(汉),河南,博士研究生。

摘要

针对计算机层析成像(CT)偏置扫描滤波反投影(FBP)重建算法中加权函数形式复杂、抑噪性能不足的问题,该文提出了一种形式简洁、连续光滑的Sigmoid指数型加权函数。通过对偏置重建过程进行Fourier变换,推导得到加权函数对白噪声的频域响应余项,相较Parker、Wang这2种加权函数,该文函数频域响应余项的高频部分幅度更低。该文分别进行了扇形束与锥形束几何在无噪声和有噪声情况下的偏置扫描仿真实验,重建图像采取主客观相结合的方法进行评价。结果表明:在理想无噪声条件下,3种加权函数重建图像无显著差异;在有噪声条件下,该文函数重建图像整体平滑、边缘锐利、对比度突出;F-范数、均方差(MSE)、对比度噪声比(CNR)这3项指标均有改善。相较Parker与Wang函数,该文函数形式简洁且抑噪性能更为优异。

关键词: 计算机层析成像(CT); 偏置扫描; 重建; 加权函数; 抑噪
中图分类号:TP391 文献标志码:A 文章编号:1000-0054(2015)01-0115-07
An exponential-type weighting function for CT reconstruction with a displaced detector array
Yi DU, Xiangang WANG, Xincheng XIANG
Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China
Abstract

For computed tomography (CT) with a displaced detector array, projections should be weighted by a weighting function to fit the filtered-backprojection (FBP) image reconstruction algorithm. Existing weighting functions have complicated piecewise forms and do not give good noise suppression. This study uses Sigmoid-like exponential weighting function, which is continuous, differentiable, and simple. Fourier transforms are used for the weighted FBP algorithm to obtain the residual terms of the frequency response for white noise. The residuals of this function are lower than those of wary function at high frequencies. The fan-beam and cone-beam geometries with displaced detector arrays were simulated with the images reconstructed using the FBP weighted by the Parker function, Wang function, and the present function. Subjective and objective evaluations show that for ideal noise-free conditions, the image qualities for the three functions are quite close; while with noisy conditions, the images given by the present function are smoother with sharper edges and higher contrast with improved F-norm, and mean squared error (MSE), contrast noise ratio (CNR) indices. Thus, this function is simpler and gives better noise suppression than the existing Parker and Wang functions.

Keyword: computed tomography (CT); displaced detector scanning; reconstruction; weighting function; noise suppression

计算机层析成像( CT)技术作为一种典型的非接触检测方法,在工业无损检测中应用广泛。近年来,随着探测器技术的发展,基于平板探测器的工业锥形束 CT得到广泛应用。然而,检测对象常常在材料、形状、尺寸等方面有较大的不确定性,当检测对象尺寸大于平板探测器尺寸时,往往会发生由于视野( field of view, FOV)无法完全覆盖被测物而出现横向数据截断问题,不满足滤波反投影重建算法的数据要求,此时,必须对扫描方案进行重新设计。

采用多次扫描并进行数据拼接或重排是一种常见的应对方案[1,2,3],然而这种方法往往要求对扫描方案进行大幅调整,如修改扫描轨迹,同时引入大量插值运算,存在工程实践复杂、数值误差大的明显不足。

偏置扫描,作为另一种常用方案,需对投影数据加权后再进行滤波反投影重建。

本文提出一种指数型加权函数,通过频域分析与仿真实验验证,对其重建抑噪性能进行研究。

1 偏置扫描

偏置扫描[4,5,6,7],作为一种将扫描几何微调与重建算法改进相结合的方法,充分利用了偏置扫描的几何对称性与采集数据中的完备性,无需增加额外的硬件和时间开销,具有操作便捷、计算效率高的优势,能够很好地应对工程中一般性大尺寸被测物扫描问题。

1.1 偏置扫描几何

基于平板探测器的锥形束CT设备,光源-探测器围绕旋转中心旋转进行360°圆轨道等角度间隔采样,每一个角度位置下沿探测器方向等距采样获取投影数据。常规的圆轨道锥形束扫描几何中, FOV必须覆盖整个被测物,光源-探测器围绕旋转中心旋转扫描一周,即可获取完整的正弦图,即Radon变换,描述如下

P(u,v)=Rnf(ξ)=∫ξf( r)·d m( r),- Ld2 uLd2.(1)

其中: ξ表示投影角度; Ld表示探测器长度; f( r)表示待测样品; d m( r)表示沿 r方向的积分弦长; P( u, v)表示以( u, v)为横纵坐标的2维投影。而偏置扫描几何中,受探测器宽度限制FOV仅覆盖投影区域的一部分,沿探测器方向存在数据截断(如图1所示),此时旋转扫描得到的投影数据,仅为非偏扫描中所得数据的一部分,可用式(2)描述。

P( u, v) = ξf( r)·d m( r), -Lo u≤( Ld-Lo).(2)

其中: u表示探测器横向坐标; Lo表示偏置扫描时探测器的重叠长度(见图1所示)。

图1 等距扇形束偏置扫描示意图

1.2 偏置数据近似完备性

Li等[7,8,9]证明了圆轨道扇形束偏置扫描条件下投影数据完备性问题,现作简述。圆轨道扇形束滤波反投影重建算法由平行束重建算法发展得到的,平行束重建以中心切片定理为基础,要求精确重建 至少需要的投影数据为

{P(α,s)|α∈[0,π], s∈[-s m,s m]}. (3)

其中: P( α, s)为平行束投影数据; α为采集投影的角度范围; s为探测器单元位置坐标; sm为探测器半宽度。扇形束全采样数据为

{ G( β, γ)∈[0, π], γ∈[m, γm]}. (4)

其中 γm为扇形束半扇角。根据采样几何关系 β=α-arcsin sR, γ=arcsin sR,可将 G( β, γ)重排为对应的 P( α, s),即

P(α,s)=Gα-arcsinsR,arcsinsR.(5)

根据平行束共轭采样原理 P( α, s) =P( α±π, -s),平行束采样投影数据存在冗余性,即沿探测器方向仅需一半的投影数据就可精确重建物体,对应扇形束偏置扫描,当重叠距离大于等于零时(即 Lo≥0),投影数据是完备的。

对于基于平板探测器的锥形束全扫描,在中间层面退化为扇形束,投影共轭关系为

P( β, u, v) =P(2 γm +π+β, -u, v). (6)

式(6)表明,偏置扫描条件下可精确重建被测物体的中间层面。在非中间层面,理论上不存在共轭采样,考虑到投影的连续性和小锥角的事实,借鉴Feldkamp-Davis-Kress(FDK)重建算法的近似思想,可认为偏置扫描采集几何存在近似的共轭关系,即投影数据近似完备,可误差较小地完成非中间层面重建。

1.3 加权重建及常用加权函数

偏置扫描可扩大扫描和重建视野,理论上当重叠区域Lo=0时可获得极限 FOV,此时获取的截断正弦图中无冗余数据,故无需特殊处理,沿用经典滤波反投影法即可完成被测物重建。然而,实际中平板探测器像素单元极小, Lo=0条件难以在工程中得到实现,一旦调整不当,即会发生因中心区域数据缺失而导致的错误重建。为避免此情况发生,偏置扫描实际操作中须保证部分探测区域发生重叠(即Lo>0),如此也减小了射术扇角/锥角,利于降低近似共轭重建的模型误差。

然而重叠区域的存在,意味着此区域内投影数据冗余,若仍沿用原有的滤波反投影算法,则在反投影重建时区域内体素值会发生多次叠加,导致重建图像中心出现明亮的盘状伪像。

为消除重叠区域影响,必须对重建算法进行改进,即在原有的经典滤波反投影算法基础上,引入加权函数,通过对投影数据进行加权,使共轭采样数据的权值之和为1,消除数据冗余区域内共轭采样的影响。

常见的方法主要有2种:预加权法与后加权法[7]。预加权法中,仅需在数据缺失部分补零,并对采集到的投影数据进行加权处理即可,处理简捷,计算十分高效,但补零会造成一定程度的 CT数下降;后加权法针对补零方法的不足,通过外推、迭代等方法对缺失数据进行估计,在滤波环节结束后再进行加权,方法准确度较预加权精确高,但依赖于数据估计方法,且牺牲了计算效率,应用中往往需在2者间进行折中。鉴于数据估计方法研究超出本研究范围,本文选择便捷高效的预加权法进行函数性能对比研究。

2 指数形式的加权函数
2.1 函数形式与特点

偏置重建中加权函数形式并不唯一,但应满足式(7)条件,即近似共轭采样对权值之和相加等于1。不仅如此,加权函数还须具有较好的光滑性[10],以避免阶跃变化或奇异点引入的频谱泄露效应而导致的重建伪像。

k-u,v)+k(u,v)=1,-LouLo;k(u,v)=1,Lo<uLd-Lo.(7)

加权函数形式多样,经典的 Parker函数出现最早,应用也最多[5,6,7],如式(8)所示。

k p(u,v)= sin2π4·u+LoLo,-Lou<0;cos2π4·u-LoLo,0u<Lo;1,LouLd-Lo.(8)

此后 Wang[11,12]从扇形束等角扫描滤波反投影重建公式出发,通过变量替换关系,提出了一种以正弦函数为架构的加权函数,如式(9)所示,并取得了很好的仿真实验结果。

k w(u,v)= 12sinπ2·arctanuRarctanLoR+1,-Lou<Lo;1,LouLd-Lo.(9)

其中R表示光源—旋转中心距离。

然而这2种函数均以分段形式给出,形式较为复杂,不够简洁。

在此本文提出了一种 Sigmoid型的指数函数形式的加权函数,如式(10)所示。

k d(u,v)= 11+exp-T·arctanuRarctanLoR,-Lo≤u≤Ld-Lo. (10)

其中T表示边界权重变换因子,由给定的边界权重α代入式(11)求得。

T=- ln(1/α-1). (11)

相较 Parker Wang函数,此加权函数借鉴了 Wang函数中权值随R变化的特点,且在整个取值区间内以唯一表达式给出,形式简单,具有无限可导性,连续光滑,同时权值在靠近数据截断一侧较小,远离侧较大,使加权后向数据完备侧集中(如图2所示)。本函数的不足之处在于,在使用时需设定重叠边界点权值,并在非重叠区域上权值逼近1,而不等于1,然而由此引起的数值误差在计算精度要求范围内可基本忽略。

图2 文中3种加权函数对比示意图

2.2 抑噪特性分析

实际应用中,方法的抑噪性能至关重要,偏置扫描中投影数据发生截断,数据量相应减少,加权函数的抑噪性能更为凸显。为简化书写,现以极坐标下等角采样扇形束偏置扫描为例,对加权重建过程中的噪声传递进行分析。

实际测量得到的投影数据,可认为由2部分相加构成,即

P(β,γ)=P(β,γ)+Noise(β,γ).(12)

其中:P(β,γ)表示投影真值; Noise(β,γ)表示噪声。

根据中心切片定理和滤波反投影算法,重建图像可表示为

(13)

其中: 分别表示Fourier与其逆变换; k( γ)表示加权函数; H( ω)表示滤波函数。将式(12)代入(13),得到噪声引起的误差分布图像

e( r) = 02π F-1[ F[Noise( β, γ k( γ)]· H( ω)]d β.(14)

对式(14)左右两端取Fourier变换,根据卷积定理,可得

E( r) = 02π[ N( β, ω) *K( ω)]· H( ω)d β.(15)

其中: N( β, ω)表示噪声的Fourier变换; K( ω)表示加权函数的Fourier变换。 CT系统中存在光子统计噪声、探测器暗电流、量子化噪声等多个噪声源,一般认为系统总的噪声服从白噪声模型,此时 N( β, ω)为常数,故可得

E(r)02πK(ω)·H(ω)dβ.(16)

其中 H( ω)为Ramp滤波函数,具有显著的高通特性,故噪声水平取决于 K( ω)的带通性能,高频响应越低,噪声通过份额越少,抑噪性能越好。

鉴于Parker函数与Wang函数相似度较高,仅选取Wang函数与本文函数进行抑噪性能对比。将Wang函数与本文函数在等角间隔条件下进行重新表述,并略去投影角度变量 β,分别如式(17)与式(18)所示。

kw( γ) = 12sinπ2·βγ0+1,-γoγ<γo;1,γoγγm.(17)

kd( γ) = 11+exp-T·γγ0,o u γm.(18)

其中: γo表示重叠角范围; γm表示扇形束半扇角。

对式(17)取Fourier变换,得

Kw( ω) = - kw( γ)·e -jωγd γ= 1 + δ(0) +Rw( ω).(19)

式(19)中余项为

Rw( ω) =cos( γ0 ω2ωjπ24γ02-ω2.(20)

对式(18)取Fourier变换,得:

Kd( ω) = - kd( γ)·e -jωγd γ= 1 + δ(0) +Rd( ω).(21)

式(21)中余项为

Rd( ω) = 2γ0ωsin( γ0 ω) + j2γ0ωcos(γ0ω)-1γ0ωsin(γ0ω).(22)

由于 H( ω)的直流分量为0,故将 Kw( ω)、 Kd( ω)与分别 H( ω)相乘结果等于各自余项 Rw( ω)、 Rd( ω)与 H( ω)相乘。现将 Rw( ω H( ω)与 Rd( ω H( ω)的频率响应进行幅值归一化处理,即得图3。从图中可看出,本文函数的高频响应明显低于Wang函数,故抑噪性能应较之突出。

3 仿真实验研究
3.1 实验条件与参数

为研究本文提出的加权函数偏置重建效果,进行了如下2组Shepp-Logan头模仿真实验: 1) 无噪声条件下扇形束与锥形束偏置扫描重建; 2) 有噪声条件下扇形束偏置扫描重建。

这2组仿真实验扫描与重建几何参数一致:光源-旋转中心-探测器距离分别为600与300(常用的2∶1比例),头模尺寸为256 × 256(× 256),探测器尺寸为313(× 512),水平重叠区域长57,角度采样数为360。其中,在实验2中,选取Gauss噪声模型,Gauss噪声的均值为0,方差为(0.01)2

3.2 实验结果

3.2.1 无噪声条件下扇形束与锥形束重建结果

本组实验用于研究本文提出的加权函数在无噪声理想情况下重建效果,分别进行扇形束与锥形束偏置扫描重建,相应重建结果与对比分如图4图5所示。

图4 扇形束全扫描与偏置扫描重建结果对比

图5 锥形束全扫描与偏置扫描重建结果对比

从扇形束重建结果(图4a—图4d)中可看出,本文提出的加权函数重建结果图像清晰,边缘明显,未出现严重伪像,剖面线对比(图4e)表明,与全扫描、Parker函数、Wang函数基本一致。图5为锥形束第90层(锥角约为4°)重建结果,结论与扇形束一致。

本组实验结果表明,本文提出的加权函数在无噪声条件下可得到理想的偏置扫描重建结果。

3.2.2 噪声条件下扇形束重建结果

本组实验在3.2.1基础上引入Gauss噪声,以扇形束偏置扫描为例,对噪声干扰下3种加权函数的噪声抑制效果进行比较研究,用于验证2.2中理论分析。3.2.1结果表明,相应结论可推广至锥形束偏置扫描。

3种不同加权函数的锥束偏置重建结果分别如图6图7图8 所示。为便于重建图像质量对比,选取2块方形感兴趣区域(region of interest,ROI)分别进行主观与客观评价。

作为主观评价,对图6a、图7a、图8a中白线区域进行局部放大分别得到图6b、图7b、图8b,可看出:在相同投影噪声条件下, Parker函数与Wang函数重建结果噪声颗粒明显,同时区域内低对比度物体边界模糊,导致目标难以分辨,而本文函数重建结果噪声较之降低,图像更加平滑,边缘保持优于Parker函数与Wang函数。

作为客观评价,选用误差F-范数、均方误差(MSE)、对比度—噪声比(CNR)这3个常用的图像质量评价指标,分别如式(23)、式(24)和式(25)定义。

Error F= i,j=1(fi,j-fi,j)2, (23)

MSE= i,j=1(fi,j-fi,j)2N2, (24)

CNR= contrastσnoise2. (25)

其中:F-范数与MSE指标的数值越小,表示图像噪声水平越低,进而说明加权函数抑噪性能越好, CNR则相反。

作为对比,选取在图6a、图7a、图8a中上部白色框内区域作为ROI分别对F-范数、MSE和CNR进行统计,结果如表1所列。

表1 F-范数、均方误差、对比度—噪声比计算结果

表1可看出,在相同噪声条件下,Parker函数与Wang函数性能接近,而本文函数3项指标均表现优异,尤其在反映抑噪和对比度保持的CNR指标上最为突出。

综合主观与客观评价结果,本文加权函数在抑噪性能上要好于Parker与Wang函数。

本文函数也存在一定局限:所使用的加权函数以 e为指数,无法对重叠区域内的共轭数据进行相对权重调整,当改变指数时,将得到一组加权函数簇,若指数趋于无穷大,即为1/2加权[11],说明指数选择对重建性能存在影响,这种影响因素需做进一步探讨,同时也为偏置扫描重建问题中的加权函数选择和设计提供一个新的参考。

4 结论

本文利用偏置扫描方法来解决基于平板探测器的大尺寸被测物CT检测问题,在原有的加权重建算法基础上,提出了一种新型加权函数,通过理论推导和实验验证,得到以下主要结论:

1) 本文所提出的新型加权函数,具有形式简单、连续光滑、计算效率高的特点,经仿真实验验证后表明,在无噪声或正常测量条件下,其重建图像无明显伪像,准确度较高。

2) 与已有的2种常用加权函数对比研究表明,本文所提出的加权函数的高频响应余项较低,可有效抑制高频噪声影响,仿真实验研究表明:其重建图像更为平滑,边缘较为清晰,对比度表现突出,具有更好的抑噪性能。

The authors have declared that no competing interests exist.

参考文献
[1] 赵飞, 路宏年, 孙翠丽. 一种新的二维 CT 扫描方式及其重建算法[J]. 光学技术, 2006, 32(2): 284-289.
ZHAO Fei, LU Hongnian, SUN Cuili. New scan mode for 2D-CT and its reconstruction algorithm[J]. Optical Technique, 2006, 32(2): 812-817. (in Chinese) [本文引用:1] [CJCR: 0.416]
[2] 陈云斌, 李寿涛, 王远, . 转台一次偏置扫描的 ICT 重建算法[J]. 中国体视学与图像分析, 2011, 16(3): 248-253.
CHEN Yunbin, LI Shoutao, WANG Yuan, et al. An ICT reconstruction algorithm for rotation off-centered scan[J]. Chinese Journal of Stereology and Image Analysis, 2011, 16(3): 248-253. (in Chinese) [本文引用:1] [CJCR: 0.232]
[3] 傅健, 路宏年. 工业 CT 半扫描成像技术[J]. 北京航空航天大学学报, 2005, 31(9): 966-969.
FU Jian, LU Hongnian. Half-scan mode for industrial CT[J]. Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(9): 966-969. (in Chinese) [本文引用:1] [CJCR: 0.526]
[4] Cho P S, Rudd A D, Johnson R H. Cone-beam CT from width-truncated projections[J]. Computerized Medical Imaging and Graphics, 1996, 20(1): 49-57. [本文引用:1] [JCR: 1.496]
[5] Mettivier G, Russo P, Lanconelli N, et al. Cone-beam breast computed tomography with a displaced flat panel detector array[J]. Medical Physics, 2012, 39(5): 2805-2819. [本文引用:1] [JCR: 3.012]
[6] ZHANG Guozhi, Jacobs R, Nuyts J, et al. Assessment of the central artifact in cone beam CT imaging with an offset geometry [C]//Medical Imaging 2012: Physics of Medical Imaging. San Diego, USA: Society of Photo-Optical Instrumentation Engineers, 2012: 83132U1-83132U7. [本文引用:1]
[7] Schafer D, Grass M, Vav De Haar P. FBP and BPF reconstruction methods for circular X-ray tomography with off-center detector[J]. Medical Physics, 2011, 38(7): S85-S94. [本文引用:4] [JCR: 3.012]
[8] LI Liang, CHEN Zhiqiang, ZHANG Li, et al. A new cone-beam X-ray CT system with a reduced size planar detector[J]. High Energy Physics and Nuclear Physics, 2006, 30(8): 812-817. [本文引用:1] [JCR: 0.233] [CJCR: 0.281]
[9] 李亮, 陈志强, 张丽, . 一种用于小体积偏置探测器锥束 CT 系统的反投影滤波重建算法[J]. CT 理论与应用研究, 2007, 16(1): 1-9.
LI Liang, CHEN Zhiqiang, ZHANG Li, et al. A BPF-type reconstruction algorithm for cone-beam CT system with a reduced size detector[J]. Computerized Tomography Theory and Applications, 2007, 16(1): 1-9. (in Chinese) [本文引用:1] [CJCR: 0.659]
[10] Parker D L. Optimal short scan convolution reconstruction for fan beam CT[J]. Medical Physics, 1982, 9: 254-257. [本文引用:1] [JCR: 3.012]
[11] WANG Ge. X-ray micro-CT with a displaced detector array[J]. Medical Physics, 2002, 29(7): 1634-1636. [本文引用:2] [JCR: 3.012]
[12] Zamyatin A A, Taguchi K, Silver M D. Helical cone beam CT with an asymmetrical detector[J]. Medical Physics, 2005, 32(10): 3117-3127. [本文引用:1] [JCR: 3.012]