光学学报, 2019, 39 (7): 0711003, 网络出版: 2019-07-16   

基于加权总差分最小化的中子稀疏投影计算机断层重建方法 下载: 840次

Neutron Computed Tomography Reconstruction Method Using Sparse Projections Based on Weighted Total Difference Minimization
作者单位
1 北京航空航天大学机械工程及自动化学院, 北京 100191
2 中国工程物理研究院核物理与化学研究所, 四川 绵阳 621900
摘要
为提升高噪声稀疏角度投影条件下中子计算机断层扫描(CT)质量,提出同时迭代重建方法(SIRT)与加权总差分最小化(WTDM)相结合的迭代重建方法(SIRT-WTDM)。在有无噪声情况下比较代数重建算法、联合代数重建算法及同时迭代重建算法的重建图像,证明了SIRT迭代重建具有较高的图像重建精度与较强的抗噪声性能,因此将SIRT作为高噪声中子投影图像CT迭代重建算法的保真项。考虑到对图像梯度稀疏性与连续性的约束,中子CT迭代重建方法的正则化约束项采用WTDM方法。由Shepp-Logan模体与真实冷中子层析扫描数据验证可知,在极端稀疏角度投影条件下,SIRT-WTDM可获得较好的重建效果。
Abstract
Aim

ing at improving the quality of the neutron computed tomography (CT) reconstructed from high noise and sparse angle projection data, an iterative reconstruction method (SIRT-WTDM) combined the simultaneous iterative reconstruction technique (SIRT) and weighted total difference minimization (WTDM) is successfully proposed. The reconstructed images obtained by algebraic reconstruction technique, simultaneous algebraic reconstruction technique, and SIRT are compared with or without the random noise in the projections, from which the SIRT method is proved to have higher reconstruction accuracy and stronger anti-noise ability. Therefore, the SIRT method is adopted as the fidelity term of the neutron CT iterative reconstruction method with high-noise projections. Considering the constraint to the sparsity and the continuity of the image gradient, the WTDM method is adopted as the regularization term of the neutron CT iterative reconstruction method. Under the condition of extreme sparse angle projections, the SIRT-WTDM can obtain the better reconstruction images, which has been proved by the Shepp-Logan simulated data and cold neutron CT scanning data.

1 引言

近年来,中子成像作为一项重要的无损检测技术,在航空航天、材料学、核物理学、地质学、生物医学、电子、考古等领域的无损检测中发挥着越来越重要的作用[1-4]。由于低中子产额、低转换效率、高随机噪声水平等因素的影响,在工程实践中为了得到满足工程检测需求的理想重建图像,最直接有效的措施是降低探测器帧频(等同于增大积分时间)和增大扫描投影幅数,但是其直接后果是导致扫描时间成倍增加,这在工程中难以接受。另外,长时间的中子辐射会进一步加大样品被活化的风险。

文献[ 5]指出,采用稀疏角度投影计算机断层扫描(CT)技术可有效解决扫描时间过长和样品辐照活化问题,并使得投影数减少、投影噪声增加。如何利用较少的投影数据重建出理想的断层图像是CT领域的研究热点之一。常用的CT重建算法一般可分为解析重建算法和迭代重建算法两类。解析重建算法的思想清晰、计算速度快,已广泛应用于各类商业CT系统,是应用最广的重建算法,但其对投影数据要求较高,如投影数目必须足够多,以满足Shannon/Nyquist采样定理且投影角度序列必须在0°~180°(或0°~360°)范围内均匀分布等。国际上采用的标准中子CT模式通常在0°~180°(或0°~360°)范围内以较小的角度间隔均匀采集数百个甚至数千个投影。迭代重建算法需要反复迭代,计算量庞大,执行速度缓慢一直成为限制其发展的重要原因。随着计算机硬件技术的快速发展及各种加速技术的出现,迭代重建算法逐渐成为研究热点。与解析重建算法相比,迭代重建算法要求的投影数据变少,在每次迭代重建的基础上可加入不同的先验知识和约束条件,因而迭代重建算法的适用范围更广泛。

稀疏角度投影下图像重建主要采用迭代重建算法[6],如迭代滤波反投影算法[7]、代数重建算法(ART)[8]、同时迭代重建算法(SIRT)[9]及介于两者之间的联合代数重建算法(SART)[10]等。ART、SIRT及SART均采用迭代更新的策略更新重建图像,相比于解析重建算法,迭代重建算法可设置不同的迭代次数与重建参数,以完成对重建图像的多次校正,并提高图像的重建精度。3种迭代重建算法的主要区别为:ART采取ray-by-ray的更新模式,每次仅更新一条射线穿过的像素,图像更新频率最快,因此图像重建收敛速度最快;SART采取view-by-view的更新模式,每次更新每个投影角度下所有射线经过的像素点;SIRT采用point-by-point的更新模式,每个图像点的更新均需计算所有投影角度下射线对当前点的贡献值,并由所有经过该点的射线贡献值取平均获得,因而重建质量最好。

随着压缩感知理论的发展,基于总变分(TV)的代数迭代重建理论在稀疏角度投影重建方面具有巨大的潜力[11-12],Rudin等[13]提出了非线性总变分降噪模型,将总变分作为正则约束,并用于图像降噪,在降噪的同时较好地保留图像边缘和细节信息。Sidky等[14-15]提出了TVM-POCS (total variation minimization-projection onto convex sets)和ASD-POCS (adaptive steepest descent-projection onto convex sets)算法,并在稀疏投影重建中取得了良好的效果。张海娇等[16]针对能谱CT重建,利用不同能量通道下重建图像具有结构相似性,提出一种基于结构先验的加权非局部全变分(NLTV)重建算法,该方法在复杂模型和高噪声模型的能谱重建方面具有明显优势。余维等[17-18]针对总差分最小化(TDM)理论研究了加权总差分最小化(WTDM)重建算法,该方法在不完备数据重建中具有良好的效果。

本文针对冷中子稀疏角度CT扫描投影数据量较少与随机噪声严重的特征,讨论冷中子断层成像质量提升问题,对比ART、SART和SIRT在有无噪声情况下的重建图像质量,选取得到的最优迭代重建方法;考虑图像梯度的稀疏性与连续性约束,将加权总差分测度引入到中子稀疏角度投影下的CT重建,有效保护了图像边缘不同方向的梯度信息。模拟数据与真实冷中子稀疏角度投影数据验证了SIRT-WTDM重建算法在高噪声稀疏角度投影下中子CT重建效果。

2 SIRT迭代重建算法

在迭代重建算法中,待重建图像被离散化为M个像素点的图像,每个像素点的值为uj'(j'=1,2,…,M),探测器单元记录的投影数据为N,投影值记为pi'(i'=1,2,…,N),其中i'j'为序号。当一束单能射线穿过物质时,其强度呈指数规律衰减,遵循朗伯-比尔(Lambert-Beer)定律。通过探测器对射线强度进行定量测量可得[19]

p=μ(x,y)dL=-ln(I/I0),(1)

式中:I为探测器单元接收到的衰减后射线辐射强度;I0为探测器单元接收到的衰减前射线辐射强度;p为物质衰减系数函数μ(x,y)在积分路径L上的积分值。SIRT迭代重建方法的方程为

uj'(k+1)=uj'(k)+λ(k)i'=1Nai'j'pi'-m'=1Mai'm'um'm'=1Mai'm'i'=1Nai'j',(2)

式中:uj'为待重建图像的像素值;k为SIRT算法迭代次数;ai'j'为系统矩阵元,是重建像素j'对第i'条射线投影值pi'的贡献度;λ为松弛因子(0<λ<2);m'为重建图像像素序号。在本文算法中,将每一条射线束看作一条无限窄的直线,射线i'穿过重建像素j'的长度定义为ai'j',因此0≤ai'j'2(1≤i'N,1≤j'M)。

针对Shepp-Logan模型[20],设计了平行束CT模拟程序,在0°~180°范围内均匀采集Shepp-Logan模型360幅投影,验证ART、SART及SIRT的重建效果。ART、SART与SIRT的投影数据访问顺序默认为顺序访问,即按照CT扫描的实际顺序,从第一个投影顺序访问至最后一个投影。当松弛因子λ=1.5、迭代次数为300时各个算法的重建图像如图1所示。受中子源强度、探测器转换效率、辐射本底、信号噪声和记录系统随机噪声等因素影响,中子照相与CT的信噪比较低,比X射线图像信噪比低1~2个数量级。因此,中子CT重建算法应具有较强的抗噪声性能,以获取更多的重建图像细节信息。为验证ART、SART及SIRT的抗噪声特征,在Shepp-Logan模型的归一化投影正弦图中加入平均值为0、标准差为0.01的高斯噪声,迭代重建参数保持不变,最终重建图像如图2所示。由图2可知,SIRT的重建图像与标准Shepp-Logan模拟模体的总体差距非常小,且像素灰度值更加平稳,具有较高的稳定性。

图 1. Shepp-Logan重建图像。(a)理想图像;(b) ART算法;(c) SART算法;(d) SIRT算法

Fig. 1. Reconstructed images of Shepp-Logan. (a) Ideal image; (b) ART algorithm; (c) SART algorithm; (d) SIRT algorithm

下载图片 查看所有图片

图 2. 含噪声投影数据下Shepp-Logan重建图像。(a)理想图像;(b) ART算法;(c) SART算法;(d) SIRT算法

Fig. 2. Reconstructed images of Shepp-Logan using noisy projections. (a) Ideal image; (b) ART algorithm; (c) SART algorithm; (d) SIRT algorithm

下载图片 查看所有图片

为定量描述有无噪声情况下3种重建方法的重建图像质量,采用均方误差(MSE)与峰值信噪比(PSNR)来衡量图像的重建质量。MSE与PSNR的数学表达式分别为

fMSE=1m·ni=0m-1j=0n-1u(i,j)-u^(i,j)2,(3)fPSNR=20·lgMIfMSE,(4)

式中:mn为图像的尺寸;ij分别为图像的横、纵坐标;u为重建后图像; u^为理论图像。MSE可衡量两张图像之间的平均像素差距,MSE越小,表示重建图像与真实图像的像素灰度差距越小,重建质量越高。PSNR为信号最大可能功率与影响它的表示精度的破坏性噪声功率的比值,MI为图像像素的最大取值,对于16位图像,MI=65535,同理对于8位图像,MI=255。PSNR越大,表示图像的失真越小。在有无噪声的情况下,ART、SART及SIRT重建图像的MSE与PSNR指标如表1所示。在有无噪声情况下SIRT的MSE均为最小,且PSNR均为最大,可以证明SIRT在有无噪声的情况下重建质量最高。此外,比较了ART、SART及SIRT 3种方法在含随机噪声时相比于无噪声时重建图像的MSE、PSNR的上升率与下降率。如表1所示,在有噪声情况下,SIRT的MSE与PSNR的上升率与下降率分别为6.52%和0.36%,均为3种方法中的最小值,因此证明SIRT比ART与SART具有较好的抗噪声性能。

表 1. Shepp-Logan重建图像的MSE和PSNR

Table 1. MSE and PSNR of reconstructed images of Shepp-Logan

MethodWithout noise in the projectionsWith noise in the projectionsImage quality reduction rate
MSEPSNR /dBMSEPSNR /dBIncreasing rate of MSE /%Reduction rate of PSNR /%
ART0.004771.45160.007869.183265.963.17
SART0.004771.42080.007769.250863.833.04
SIRT0.004671.51670.004971.25716.520.36

查看所有表

3 稀疏角度投影SIRT图像重建

采用稀疏角度投影CT扫描时,采集得到的投影数据量不能满足Shannon/Nyquist采样定理,滤波反投影重建算法(FBP)的重建图像会出现欠采样条状伪影。随着投影数据量的减少,重建图像中伪影情况大幅度加剧。图3所示为0°~180°范围内稀疏程度不同的投影条件下Shepp-Logan模型的FBP以及SIRT重建结果。对比图3,在投影数据比较完备的条件下(180幅投影),采用FBP与SIRT均能精确地恢复图像的轮廓信息。当投影幅数为20幅时,FBP与SIRT重建图像中像素灰度值出现不符合实际模型的不连续性,像素灰度分布不均匀,图像细节信息大量丢失,视觉品质大幅降低,相比较而言,SIRT的重建图像中没有出现明显条状伪影,总体图像重建质量略优于FBP。为此,本研究将基于加权总差分最小化的正则化约束加入到稀疏角度投影SIRT重建中,在每次SIRT迭代重建后,均采用加权总差分最小化的方法对图像各个方向的梯度进行最小化约束,以得到稀疏投影角度下最优CT重建质量。

图 3. 不同稀疏角度投影FBP与SIRT重建图像。(a) 180幅投影FBP重建图像;(b) 20幅投影FBP重建图像;(c) 180幅投影SIRT重建图像;(d) 20幅投影SIRT重建图像

Fig. 3. Reconstructed images obtained by FBP and SIRT with different sparse angle projections. (a) FBP reconstructed image with 180 projections; (b) FBP reconstructed image with 20 projections; (c) SIRT reconstructed image with 180 projections; (d) SIRT reconstructed image with 20 projections

下载图片 查看所有图片

4 基于加权总差分最小化的SIRT图像重建

基于总差分最小化[21]约束的代数重建算法针对角度稀疏型CT重建问题取得了很好的效果。然而总差分正则化仅考虑了图像梯度的稀疏性,对物体边缘结构信息的保护稍显不足。余维等[17-18]研究了基于加权总差分最小化的代数重建算法。加权总差分最小化不仅考虑了图像梯度的稀疏性约束,还考虑了图像梯度的连续性约束,能对图像边缘不同方向的梯度信息进行有效保护。本研究将加权总差分最小化测度作为高噪声压缩采样中子图像迭代重建的正则化约束,进一步研究WTDM正则化方法在中子稀疏角度投影重建中的应用。

投影数据迭代重建问题,可以等价于优化问题的求解过程,即

u'=argminu'Au'-g22+R(u'),(5)

式中:g为实测投影数据;u'为待重建断层图像;A为当前扫描模式下对应的系统矩阵;‖Au'-g22为数据保真项,代表迭代重建过程中断层模拟投影与实测投影之间的误差项;R(u')为重建图像所要满足的图像先验信息及准则项,是迭代重建算法的正则项。

加权总差分要实现对图像梯度的稀疏性和连续性的双重约束。对于梯度的稀疏性,约束任意一个像素点与其邻域内像素点像素保持一致;对于梯度的连续性,则需要约束任意一个像素位置与邻域内其他像素点的局部梯度保持一致[22]DhDv分别为水平和竖直方向的梯度算子;DvhDhv为对角方向上的梯度算子;DhhDh沿竖轴方向的导数;DvvDv沿横轴方向的导数。图4(a)为图像任意像素点Pi,j的四邻域示意图,图4(b)为图像任意像素点Pi,j各方向的局部梯度。对图像梯度的约束可等价为对各个像素点水平垂直差分和对角差分的联合约束。

图 4. WTDM方法在不同方向的差分测度[17]。(a)像素点四邻域示意图;(b)像素点各方向的局部梯度

Fig. 4. Difference measures of WTDM method in different directions[17]. (a) Diagram of four neighbors of pixel point; (b) partial gradients of pixel point in different directions

下载图片 查看所有图片

对局部梯度稀疏性和局部梯度连续性的约束差分可分别表示为

TD(Pi,j)1=DhPi,j1+DvPi,j1,(6)TD(Pi,j)2=DhvPi,j1+DvhPi,j1+DvvPi,j1+DhhPi,j1,(7)

其中

DhPi,j1=Pi+1,j-Pi,jDvPi,j1=Pi,j+1-Pi,jDvhPi,j1=Pi,j+1-Pi+1,jDhvPi,j1=Pi+1,j+1-Pi,j,(8)DhhPi,j1=Pi+1,j+1+Pi,j-Pi+1,j-Pi,j+1=DvvPi,j1,(9)

则对局部梯度连续性约束差分可简化为

TD(Pi,j)2=DhvPi,j1+DvhPi,j1+2·DhhPi,j1(10)

加权求和得到加权总差分,实现对梯度连续性和稀疏性的双重约束:

WTD(Pi,j)=DhPi,j1+DvPi,j1+α(DhvPi,j1+DvhPi,j1+2·DhhPi,j1),(11)

式中:α为平衡梯度连续性与梯度稀疏性约束之间的权值。进一步分析发现,差分约束的目的是实现图像整体梯度的最大稀疏效果。在极限情况下, DhPi,j1+ DvPi,j1+ DhvPi,j1+ DvhPi,j1已取得最小化,即 DhPi,j1DvPi,j1DhvPi,j1DvhPi,j1均为零时, DhhPi,j1DvvPi,j1必定达到最小化,因此可忽略对于 DhhPi,j1DvvPi,j1的约束,简化后的加权总差分为

WTD(Pi,j)=DhPi,j1+DvPi,j1+α(DhvPi,j1+DvhPi,j1)(12)

α=0 时,则可实现对图像梯度稀疏性的约束;当α=1时,则可以实现对梯度稀疏性和连续性同等程度的约束(α默认为1)。将WTDM正则项约束代入投影数据迭代重建问题[(5)式]中,则图像的迭代重建模型转化为[17-18]

u'=argminu'Au'-g22+2ω·ij[|ui+1,j-ui,j|+|ui,j+1-ui,j|+α(|ui+1,j+1-ui,j|+|ui,j+1-ui+1,j|)],(13)

式中:ui,j为重建图像(i,j)坐标下的像素值;ω为加权总差分最小化的作用程度。ω越大,加权总差分最小化作用程度越强,重建图像梯度的变化幅度较大;ω越小,加权总差分最小化作用程度越弱,重建图像梯度的变化幅度较小。因此对于加权总差分最小化CT重建模型的求解表示为(13)式的求解。(13)式的求解有多种方法,其中软阈值滤波方法的收敛性和有效性在理论上已得到证明,并成功应用于CT重建[23]。本研究选用软阈值滤波方法实现(13)式的求解。在软阈值滤波的框架下对基于加权总差分最小化重建模型进行求解,需要构造一个伪逆[17-18,21]:

ui,jn+1=14+4·α×{f(ω,u~i,jn+1,u~i+1,jn+1)+f(ω,u~i,jn+1,u~i,j+1n+1)+f(ω,u~i,jn+1,u~i,j-1n+1)+f(ω,u~i,jn+1,u~i-1,jn+1)+α·[f(ω,u~i,jn+1,u~i+1,j+1n+1)+f(ω,u~i,jn+1,u~i+1,j-1n+1)+f(ω,u~i,jn+1,u~i-1,j-1n+1)+f(ω,u~i,jn+1,u~i-1,j+1n+1)]},(14)

式中: u~i,jn+1为迭代重建图像上(i,j)坐标位置处灰度值;f为软阈值函数。

f(ω,y,z)=(y+z)/2,|y-z|<ωy-ω/2,(y-z)ωy+ω/2,(y-z)-ω,(15)

式中:yz为软阈值函数f的输入参量。

余维等[17-18]针对模拟模体的稀疏角度投影数据,采用交替最小化的方式实现对(13)式的求解:首先利用SART进行重建以缩小原始投影数据与模拟投影数据之间的差异,在每次SART重建后进行一次软阈值滤波,从而减小重建图像的加权总差分,最终得到理想的重建图像,因此其算法中一次SART-WTDM算法主循环主要包括一次SART重建与一次WTDM过程。本研究针对中子投影数据噪声高、图像质量差的特点,提出了SIRT-WTDM重建方法。SIRT-WTDM重建方法主要对余维等[17-18]提出的SART-WTDM方法进行两个方面的调整:1)将SIRT作为基于WTDM重建方法的保真项,以增强算法的抗噪声性能;2)设置SIRT-WTDM算法主循环中WTDM的次数为NTD。这样可将一次SIRT-WTDM算法主循环中WTDM的次数由1调整为NTD,进一步恢复欠采样投影重建图像的细节信息,提升中子稀疏角度投影CT质量。因此一次SIRT-WTDM算法主循环可表示为一次SIRT重建与NTD次WTDM过程。

5 实验验证

5.1 模拟数据的实验验证

针对Shepp-Logan模型,本研究在0°~180°范围内均匀采集90幅投影,比较FBP、SIRT、SART-WTDM及SIRT-WTDM算法在稀疏角度投影时的CT重建图像质量,Shepp-Logan模型的CT条件如表2所示。

表 2. 基于Shepp-Logan模拟的CT条件

Table 2. Simulative CT conditions based on Shepp-Logan

Scanning modeReconstructing image size /(pixel×pixel)Range of scanning angle /(°)Angle step /(°)
Parallel beam256×2561802

查看所有表

SIRT的松弛因子λ=1.5;设置SART-WTDM算法中SART的松弛因子λ=0.1,NTD=1,以保证与文献[ 18]参数一致;设置本文SIRT-WTDM算法中SIRT的λ=1.5,NTD=1,ω=0.00035;设置SIRT、SART-WTDM及SIRT-WTDM算法主循环次数为700。FBP、SIRT、SART-WTDM及SIRT-WTDM算法整体重建图像以及对应局部放大图像如图5所示。为验证4种重建算法的抗噪声性能,

图 5. 90幅投影Shepp-Logan的重建图像。(a) FBP重建;(b) SIRT重建;(c) SART-WTDM重建;(d) SIRT-WTDM重建

Fig. 5. Reconstructed images of Shepp-Logan using 90 projections. (a) FBP reconstruction; (b) SIRT reconstruction; (c) SART-WTDM reconstruction; (d) SIRT-WTDM reconstruction

下载图片 查看所有图片

在Shepp-Logan模型的归一化投影正弦图中加入平均值为0、标准差为0.01高斯噪声,设置SIRT-WTDM算法的NTD=2,ω=0.0005,其余重建参数保持不变。含噪声投影条件下FBP、SIRT、SART-WTDM及SIRT-WTDM算法整体重建图像以及对应局部放大图像如图6所示。

图 6. 90幅含噪声投影Shepp-Logan的重建图像。(a) FBP重建;(b) SIRT重建;(c) SART-WTDM重建;(d) SIRT-WTDM重建

Fig. 6. Reconstructed images of Shepp-Logan using 90 noisy projections. (a) FBP reconstruction; (b) SIRT reconstruction; (c) SART-WTDM reconstruction; (d) SIRT-WTDM reconstruction

下载图片 查看所有图片

图5图6可知,稀疏角度投影对重建图像的影响主要表现为重建图像中的条状伪影以及图像细节不清晰。相比于FBP与SIRT,SART-WTDM、SIRT-WTDM重建结果中条状伪影基本消除,图像中像素变化平稳,细节信息较为清晰,重建结果更优。从图6可以看出, SIRT-WTDM重建图像像素灰度值更加平滑,说明SIRT-WTDM方法对噪声的抑制更为显著。

为了更好地对图5图6中重建图像质量进行表征,采用MSE、PSNR对图像质量进行定量评价,结果如表3所示。可以看出:FBP、SIRT、SART-WTDM、SIRT-WTDM方法的MSE依次降低,PSNR依次增加。因此可以得出FBP、SIRT、SART-WTDM、SIRT-WTDM方法的重建图像质量依次递增。图5图6中4幅图像的第128列数据绘制灰度值分布曲线如图7所示,可得SIRT-WTDM方法获得图像灰度与真实模型的差距最小,并且像素灰度值更加平稳,验证了SIRT-WTDM方法可得到较高的图像质量。

表 3. Shepp-Logan的90幅投影重建图像的MSE和PSNR

Table 3. MSE and PSNR of reconstructed images of Shepp-Logan with 90 projections

MethodWithout noise in projectionsWith noise in projections
MSEPSNR /dBMSEPSNR /dB
FBP0.010068.13740.011867.4080
SIRT0.004671.45930.006470.0438
SART-WTDM0.004571.58420.006270.1929
SIRT-WTDM0.004571.58500.006270.2200

查看所有表

5.2 冷中子数据的实验验证

基于中国工程物理研究院的CMRR(China Mianyang Research Reactor)反应堆,获得了编码管样品冷中子层析扫描数据。中子CT扫描条件如表4所示,对面阵探测器采集的投影数据进行重排获得样品位于探测器第776层的投影正弦图。

在上述采集条件下,对投影正弦图左右两端冗余部分进行裁切并进行均匀采样,获得75幅投影正弦图,对75幅稀疏投影条件下的编码管样品分别进行FBP、SIRT、SART-WTDM及SIRT-WTDM重建,设置SIRT-WTDM算法的NTD=1,ω=0.00003,其余重建参数与5.1节均一致。整体重建图像以及对应局部放大图像如图8所示。受中子CT系统噪声影响,图8中FBP重建图像中出现明显条状伪影与像素突变情况;SIRT、SART-WTDM与SIRT-WTDM重建图像中条形伪影有所减轻。相比于FBP、SIRT及SART-WTDM方法,SIRT-WTDM重建结果中图像轮廓还原度更高,且噪声得到了很好的抑制,条形伪影进一步削弱,像素灰度值更加符合真实情况。

图 7. Shepp-Logan的90幅投影重建图像的灰度曲线。(a)投影数据无噪声;(b)投影数据含噪声

Fig. 7. Gray curves of reconstructed images of Shepp-Logan with 90 projections. (a) Without noise in projections; (b) with noise in projections

下载图片 查看所有图片

表 4. 编码管样品CT扫描条件

Table 4. CT scanning conditions of code tube sample

Neutron beam typeDetector typeNumber of detector pixelsRange of scanning angle /(°)Total projection number
Cold neutronAndor iKon-L 9362048×2048180450

查看所有表

图 8. 编码管样品75幅投影重建图像。(a)完整数据FBP重建;(b) 75幅投影FBP重建;(c) 75幅投影SIRT重建;(d) 75幅投影SART-WTDM重建;(e) 75幅投影SIRT-WTDM重建

Fig. 8. Reconstructed images of code tube sample with 75 projections. (a) FBP reconstruction with complete projections; (b) FBP reconstruction with 75 projections; (c) SIRT reconstruction with 75 projections; (d) SART-WTDM reconstruction with 75 projections; (e) SIRT-WTDM reconstruction with 75 projections

下载图片 查看所有图片

为了验证SIRT-WTDM算法在严重稀疏投影情况下的重建图像质量,采取更加苛刻的投影幅数对SIRT-WTDM重建算法的图像信息还原能力进行验证。对编码管样品第313层在0°~180°范围内进行均匀采样,获得18幅投影。设置SIRT-WTDM算法的NTD=2,ω=0.0001,其余重建参数与5.1节均一致。获得的重建图像如图9所示。对比图9所示的重建结果可得:FBP的重建图像受投影数据严重稀疏的影响,重建图像条状伪影非常明显,物体内部结构的识别受到严重影响。与FBP、SIRT、SART-WTDM方法相比,SIRT-WTDM方法得到的重建图像质量最佳,重建图像中物体内外部轮廓结构更加清晰,且基本没有明显条状伪影。从SIRT-WTDM方法的重建图像中可清晰地看出该样品断层图像内部具有26个柱状结构。

图 9. 编码管样品18幅投影重建图像。(a)完整数据FBP重建;(b) 18幅投影FBP重建;(c) 18幅投影SIRT重建;(d) 18幅投影SART-WTDM重建;(e) 18幅投影SIRT-WTDM重建

Fig. 9. Reconstructed images of code tube sample with 18 projections. (a) FBP reconstruction with complete projections; (b) FBP reconstruction with 18 projections; (c) SIRT reconstruction with 18 projections; (d) SART-WTDM reconstruction with 18 projections; (e) SIRT-WTDM reconstruction with 18 projections

下载图片 查看所有图片

图8图9中所有重建图像的灰度值进行归一化,把完整数据FBP重建图像当作理想图像,并以理想图像为基准采用MSE、PSNR对各个图像的重建精度进行定量评价。图8图9所示重建图像的MSE、PSNR如表5所示,可以发现SIRT-WTDM方法重建图像的MSE、PSNR均为最优,进而证明基于SIRT-WTDM的中子稀疏角度投影CT重建方法可以得到更优的重建结果。

表 5. 编码管样品重建图像的MSE和PSNR

Table 5. MSE and PSNR of reconstructed images of code tube sample

Method776th slice image313rd slice image
MSEPSNR /dBMSEPSNR /dB
FBP9.8263×10-478.20690.013166.9626
SIRT9.4837×10-478.36100.009068.6035
SART-WTDM7.6070×10-479.31870.004271.8748
SIRT-WTDM4.0869×10-482.01690.003772.4440

查看所有表

为了进一步比较FBP、SIRT、SART-WTDM与SIRT-WTDM算法对稀疏角度中子投影数据的CT重建效果,采用FBP、SIRT、SART-WTDM及SIRT-WTDM算法分别对图10(a)所示的编码管样品感兴趣区域 (ROI)进行多层重建以构建编码管样品的三维体数据,编码管样品采集投影数为75幅,重建图像大小为600 pixel×600 pixel,总重建层数为561层。

本实验中设置SIRT-WTDM算法的NTD=2,ω=0.00003;设置SIRT、SART-WTDM与SIRT-WTDM算法主循环次数为500,其余重建参数均与5.1节一致,最终得到各个方法重建的体数据结果。在三维体数据可视化之前,首先将体素灰度值线性拉伸至0~65535以满足三维可视化软件的需要。为了更清晰地比较各个方法的重建结果,采用伪彩渲染的方式显示三维体数据,其中灰度窗口统一设置为1000~19000,如图10(b)所示。编码管样品完整投影数据FBP重建三维可视化结果如图10(b)所示。75幅投影数据的FBP、SIRT、SART-WTDM及SIRT-WTDM重建三维可视化结果如图10(c)~(f)所示。

对比图10(c)~(f)所示结果,可以发现:FBP的三维重建结果噪声最大,三维可视化图像不够清晰;由于中子投影数据中存在严重噪声,SART-WTDM方法的三维重建结果中出现了与真实灰度值不符的像素灰度突变,如图10(e)中箭头所示区域;相比于FBP、SIRT及SART-WTDM方法的重建图像,SIRT-WTDM方法能够较好地抑制重建图像噪声,重建得到的样品轮廓更加清晰。观察图10(c)~(f)虚线框标记的局部区域可得:SIRT-WTDM方法得到重建图像的像素灰度值更加平稳,能够更加清晰地展现样品的内部三维结构,其三维可视化结果与图10(b)结果最为相似。因此可证明SIRT-WTDM方法能够较好地抑制中子投影数据中的随机噪声,并能够有效地减轻由投影数据缺失导致的重建伪影,最终得到较好的中子CT重建结果。

图 10. 编码管样品的三维可视化结果。(a)编码管样品的重建部位;(b)完整数据FBP重建;(c) 75幅投影数据FBP重建;(d) 75幅投影数据SIRT重建;(e) 75幅投影数据SART-WTDM重建;(f) 75幅投影数据SIRT-WTDM重建

Fig. 10. 3D visualization results of code tube sample. (a) Reconstructed region of code tube sample; (b) FBP reconstructed result with complete projections; (c) FBP reconstructed result with 75 projections; (d) SIRT reconstructed result with 75 projections; (e) SART-WTDM reconstructed result with 75 projections; (f) SIRT-WTDM reconstructed result with 75 projections

下载图片 查看所有图片

6 结论

稀疏角度投影CT是一种能够有效解决中子CT时间过长与被检样品辐照活化问题的有效方法。由于中子投影数据噪声高,稀疏角度投影CT所造成的投影数据量不足,传统的CT重建方法难以恢复被检样品的真实结构信息。为提升高噪声稀疏角度投影下中子断层成像质量,提出了采用SIRT与WTDM相结合的迭代重建方法。与ART与SART相比,SIRT具有对投影噪声敏感程度低的特点。在同等投影噪声强度情况下,SIRT重建图像的像素灰度值更加平稳,图像质量更优。考虑到图像梯度的稀疏性与连续性的约束,迭代重建方法的正则化项采用WTDM方法,并采用软阈值滤波方法实现对图像重建模型的求解。所提出的SIRT-WTDM方法实现了对编码管样品冷中子稀疏角度投影数据的三维重建,其三维可视化结果证明:相比于FBP、SIRT及SART-WTDM方法,SIRT-WTDM方法能较好地抑制中子投影数据中的随机噪声,并能有效地减轻由投影数据缺失导致的重建伪影。但是SIRT-WTDM迭代重建方法计算量大、重建速度慢、所需调整参数多的特点限制了该算法在中子CT领域的应用,因此未来的研究方向主要集中于重建算法加速与重建参数的自适应调整方面的研究。

参考文献

[1] Sinha V, Srivastava A, Koo Lee H. A novel method for NDT applications using NXCT system at the Missouri University of Science & Technology[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2014, 750: 43-55.

[2] Tamaki M. Conceptual monochromatic digital neutron radiography using continuous cold neutron beam[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2005, 542(1/2/3): 32-37.

[3] Yasuda R, Matsubayashi M, Nakata M, et al. Application of neutron imaging plate and neutron CT methods on nuclear fuels and materials[J]. IEEE Transactions on Nuclear Science, 2005, 52(1): 313-316.

[4] 魏国海, 韩松柏, 陈东风, 等. 中子照相技术在核燃料元件无损检测中的应用[J]. 核技术, 2012, 35(11): 821-826.

    Wei G H, Han S B, Chen D F, et al. Application of neutron radiography for non-destructive testing nuclear fuel elements[J]. Nuclear Techniques, 2012, 35(11): 821-826.

[5] Strobl M, Manke I, Kardjilov N, et al. Advances in neutron radiography and tomography[J]. Journal of Physics D: Applied Physics, 2009, 42(24): 243001.

[6] 司凯. 稀疏投影角度下的CT迭代图像重建算法应用研究[D]. 济南: 山东大学, 2015: 2- 14.

    SiK. Application research on CT iterative image reconstruction algorithm of sparse angles projection[D]. Jinan: Shandong University, 2015: 2- 14.

[7] 马继明, 张建奇, 宋顾周, 等. 全变分约束迭代滤波反投影CT重建[J]. 光学学报, 2015, 35(2): 0234002.

    Ma J M, Zhang J Q, Song G Z, et al. Total variation constrained iterative filtered back projection CT reconstruction method[J]. Acta Optica Sinica, 2015, 35(2): 0234002.

[8] Gordon R, Bender R, Herman G T. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[J]. Journal of Theoretical Biology, 1970, 29(3): 471-481.

[9] Gilbert P. Iterative methods for the three-dimensional reconstruction of an object from projections[J]. Journal of Theoretical Biology, 1972, 36(1): 105-117.

[10] Andersen A H, Kak A C. Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm[J]. Ultrasonic Imaging, 1984, 6(1): 81-94.

[11] 王林元, 刘宏奎, 李磊, 等. 基于稀疏优化的计算机断层成像图像不完全角度重建综述[J]. 物理学报, 2014, 63(20): 208702.

    Wang L Y, Liu H K, Li L, et al. Review of sparse optimization-based computed tomography image reconstruction from few-view projections[J]. Acta Physica Sinica, 2014, 63(20): 208702.

[12] 蔺鲁萍, 王永革. 不完全角度CT图像重建的模型与算法[J]. 北京航空航天大学学报, 2017, 43(4): 823-830.

    Lin L P, Wang Y G. CT image reconstruction model and algorithm from few views[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(4): 823-830.

[13] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms[J]. Physica D: Nonlinear Phenomena, 1992, 60: 259-268.

[14] Sidky E Y, Kao C M, Pan X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J]. Journal of X-ray Science and Technology, 2006, 14(2): 119-139.

[15] Sidky E Y, Pan X C. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Physics in Medicine and Biology, 2008, 53(17): 4777-4807.

[16] 张海娇, 孔慧华, 孙永刚. 基于结构先验的加权NLTV能谱CT重建算法[J]. 光学学报, 2018, 38(8): 0811003.

    Zhang H J, Kong H H, Sun Y G. Weighted NLTV reconstruction algorithm based on structural prior information for spectral CT[J]. Acta Optica Sinica, 2018, 38(8): 0811003.

[17] 余维. 不完备投影数据的CT重建算法研究[D]. 重庆: 重庆大学, 2014: 24- 35.

    YuW. Investigation of reconstruction algorithms for incomplete CT projections[D]. Chongqing: Chongqing University, 2014: 24- 35.

[18] Yu W, Zeng L. A novel weighted total difference based image reconstruction algorithm for few-view computed tomography[J]. PLoS One, 2014, 9(10): e109345.

[19] 闫镔, 李磊. CT图像重建算法[M]. 北京: 科学出版社, 2014: 5- 14.

    YanB, LiL. CT image reconstruction algorithm[M]. Beijing: Science Press, 2014: 5- 14.

[20] Shepp L A, Logan B F. The Fourier reconstruction of a head section[J]. IEEE Transactions on Nuclear Science, 1974, 21(3): 21-43.

[21] Yu H Y, Wang G. A soft-threshold filtering approach for reconstruction from a limited number of projections[J]. Physics in Medicine and Biology, 2010, 55(13): 3905-3916.

[22] Shu XB, AhujaN. Hybrid compressive sampling via a new total variation TVL1[M] // Daniilidis K, Maragos P, Paragios N. Lecture notes in computer science. Berlin, Heidelberg: Springer, 2010, 6316: 393- 404.

[23] Liu B D, Wang G, Ritman E L, et al. Image reconstruction from limited angle projections collected by multisource interior X-ray imaging systems[J]. Physics in Medicine and Biology, 2011, 56(19): 6337-6357.

林强, 杨民, 唐彬, 刘斌, 霍合勇, 刘家伟. 基于加权总差分最小化的中子稀疏投影计算机断层重建方法[J]. 光学学报, 2019, 39(7): 0711003. Qiang Lin, Min Yang, Bin Tang, Bin Liu, Heyong Huo, Jiawei Liu. Neutron Computed Tomography Reconstruction Method Using Sparse Projections Based on Weighted Total Difference Minimization[J]. Acta Optica Sinica, 2019, 39(7): 0711003.

本文已被 1 篇论文引用
被引统计数据来源于中国光学期刊网
引用该论文: TXT   |   EndNote

相关论文

加载中...

关于本站 Cookie 的使用提示

中国光学期刊网使用基于 cookie 的技术来更好地为您提供各项服务,点击此处了解我们的隐私策略。 如您需继续使用本网站,请您授权我们使用本地 cookie 来保存部分信息。
全站搜索
您最值得信赖的光电行业旗舰网络服务平台!