一种改进的中子外部CT图像伪影抑制方法

林强; 杨民; 张晓敏; 唐彬; 刘斌; 霍合勇

doi:doi:10.3788/AOS202040.2210001

光学学报, 2020, 40 (22): 2210001, 网络出版: 2020-10-25

一种改进的中子外部CT图像伪影抑制方法下载： 1029次封面文章

Improved Image Artifacts Suppression Method for Neutron External CT

论文大纲

林强 ¹杨民 ^1,*张晓敏 ¹唐彬 ²刘斌 ²霍合勇 ²

作者单位

¹ 北京航空航天大学机械工程及自动化学院, 北京 100191

² 中国工程物理研究院核物理与化学研究所, 四川绵阳 621900

AI 词云图 AI语音精读 AI语音超短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

径向边缘的伪影是外部CT图像明显存在的缺点。为了抑制外部CT图像径向边缘的伪影,采用加权方向全变差(WDTV)算法计算沿径向和切向方向的局部方向差分,引入两个权重参数对这两个局部方向差分进行加权求和。WDTV算法可以更好地描述外部CT图像梯度模的稀疏性,有效提升重建图像的质量。针对高噪声中子外部CT检测的需求,在WDTV重建模型的框架下,对沿着径向和切向两个方向附近的一定角度范围内分别施加不同的全变差(TV)最小化约束。改进的WDTV算法中TV最小化的作用更明显,具有更强地抑制径向边缘伪影与抗噪声的性能。计算机仿真结果和齿轮冷中子实验结果表明,改进的WDTV重建模型能够更有效地抑制重建图像的噪声,提高重建图像的质量。

Abstract

The artifacts of radial edges are obvious shortcoming of external CT image. In order to suppress the artifacts of the radial edge of external CT images the weighted directional total variation (WDTV) algorithm is used to calculate the local directional differences along the radial and tangential directions, and two weighted parameters are introduced to carry out the weighted sum of these two local directional differences. The WDTV algorithm can better describe the sparsity of external CT image gradient magnitude and effectively improve the quality of reconstructed image. In order to meet the needs of high-noise neutron external CT inspection, different total variation (TV) minimization constraints are applied in certain angle ranges near the radial and tangential directions under the framework of WDTV reconstruction model. In the improved WDTV algorithm, TV minimization plays a more obvious role and has stronger performance of suppressing radial edge artifacts and anti-noise. The results of computer simulation and gear cold neutron experiment show that the improved WDTV reconstruction model can suppress the noise of reconstructed images more effectively and improve the quality of reconstructed images.

1 引言

中子与X射线作为两种重要的辐射源,在材料学、生物医学、航空航天、汽车工业、电子、考古和安检等领域的无损检测中发挥着越来越重要的作用^[1-5]。众所周知,X射线和中子与物质的作用机理不同。X射线光子作用于物质原子中的核外电子,其作用截面与核素的原子序数有确定的函数关系,高原子序数的物质具有较大的X射线作用截面,则其在X射线图像中具有较高的对比度。中子则是直接与原子核相互作用,散射截面的尺寸与原子序数无关,同一能量的中子对于不同核素的散射截面也有较大的差异。由此可见,中子成像与X射线成像存在着很好的互补性,中子成像在检测含较轻元素(氢、锂和硼等)的材料、原子序数相近的元素、同位素和放射性材料等方面有X射线成像无法具备的优势^[6]。

外部计算机断层扫描(CT)成像多应用于工业无损检测领域,但受到待检测样品的尺寸和密度、射线能量以及探测器尺寸等限制,所以其只能对物体边缘进行照射以获取投影数据^[7]。对于火箭装药筒内壁的包覆层检测、金属管道内壁的损伤检测和压力容器内壁的保温层检测等,使用合适的外部CT策略以及相应的重建算法就可以获得较好的成像效果。当使用中子束外部CT系统扫描样品时,在2π角度的旋转过程中,射线束与探测器构成的透照场仅覆盖了样品的外部区域,感兴趣重构区域为透照场扫过的环形区域^[8],而外部轮廓区域的CT图像由专用的重建算法得到^[9-10]。

外部轮廓区域的研究在医学领域中备受关注,医生采用了外部CT成像技术实现对心脏外围组织的重建,以减少辐照区域、辐照剂量和心脏跳动带来的伪影^[11]。目前,外部CT成像技术已在工业CT的无损检测领域中发挥着重要的作用。现有的外重建算法主要有奇异值分解(SVD)算法、LT(Lambda Tomography)法和迭代法。1963年,Cormack^[12]提供了外部问题的反演公式。受到Cormack^[12]的启发,Quinto^[13-15]提出了一种SVD算法来解决外部问题。SVD算法通过构造正交基对待重建图像进行正交分解,利用已知投影和先验信息分别重建相互正交的零空间及其正交补空间,最后将二者合成为最终的重建结果^[16]。LT算法属于局部成像算法的一种,该算法将偏微分方程中的Λ算子引入CT成像中,通过求解密度函数的Λ相关函数以重建出Λf的分布图像,这不同于传统的CT图像^[17];Quinto^[11]对LT算法进行了改进,采用有界外推的方法对缺失数据进行补偿,获得了火箭发动机模型的外部重建结果。迭代重建算法将待重建图像的先验信息转换为约束条件或优化准则,将CT重建问题转换为具有约束条件的优化问题,采用迭代算法来获得逼近理想图像的最优解,如近几年备受热议和关注的TVM(Total Variation Minimization)类算法,其被认为是解决外重建问题的有效算法^[18-20];Zeng等^[21]针对外部问题提出了一种子区域平均TVM-POCS(Subregion-Averaged-TVM-Projection Onto Convex Sets)算法,数值仿真结果表明了该算法的有效性。Guo等^[22]在TVM迭代过程中加入RSF(Region Scalable Fitting)模型,证明该模型在抑制渐变伪影和硬化伪影两个方面的表现更胜一筹。Guo等^[23]基于外部CT重建图像含有伪影的特点,提出了加权方向全变差(WDTV)正则化模型以减少重建图像中在与径向相切边缘附近的伪影,计算机仿真结果和实验结果证明,基于WDTV的正则化模型更适用于解决外部问题。然而上述重建算法是基于X射线扇束、锥束和螺旋CT扫描场景而设计的,本文根据高噪声中子平行束外部CT检测的需求,对WDTV正则化模型进行改进,这可以有效提高外部CT成像的质量,进而获得冷中子外部CT图像。

2 传统的解析重建与迭代重建算法

中子束外部CT扫描原理如图1所示。为了验证不同重建算法对于平行束外部的重建效果,设计了外部CT扫描模体,该模体包含不同灰度的圆孔、不同方向的方孔以及横向与纵向的线状结构,如图2(a)所示,ROI-t表示不同区域,t=1,2,3,4,5,6,7,8。外部CT扫描参数如表1所示,CT扫描过程以及投影正弦图像如图2(b)和图2(c)所示。从图2(c)可以看到,正弦图像左侧的数据未发生截断,而正弦图像右侧的数据发生明显截断,模拟的模体尺寸为512 pixel×512 pixel,而正弦图像的宽度为128 pixel,缺失待扫描样品129~512 pixel宽的投影数据,投影数据量为完整数据量的1/4。

图 1. 中子外部CT扫描原理图

Fig. 1. Schematic of neutron exterior CT scan

下载图片查看所有图片

表 1. 外部CT扫描参数

Table 1. Parameters of external CT scan phantom

Scanning mode	Unit size of detector /pixel	Detector width /pixel	Scanning angular range /(°)	Number of projected views	Phantom size /(pixel×pixel)
2D parallel beam	1	128	0-360	720	512×512

查看所有表

为了实现FBP(Filtered Back Projection)算法的重建,首先根据平行束数据的对称性原理将正弦图像右侧的投影数据进行补充,然后采用FBP算法对其进行CT重建,FBP算法重建的图像如图3(a)所示。

图 2. 外部CT扫描模体的示意图。(a)模体;(b)扫描过程;(c)正弦图像

Fig. 2. Schematic of external CT scan phantom. (a) Phantom; (b) scanning process; (c) sinusoidal image

下载图片查看所有图片

采用SART(Simultaneous Algebraic Reconstruction Technique)算法对模体进行外部CT重建,设置SART算法的松弛因子λ=0.05,迭代次数N_iter为500次,SART算法重建的图像如图3(b)所示。从图3可以看到,FBP算法与SART算法重建的图像中存在严重的结构缺失状伪影,得到与圆周切向(tangential direction)相切的边缘较为清晰,与径向(radial direction)相切的边缘被伪影扭曲。

图 3. 不同算法重建的图像。(a) FBP算法;(b) SART算法

Fig. 3. Images reconstructed by different algorithms. (a) FBP algorithm; (b) SART algorithm

下载图片查看所有图片

3 WDTV最小化重建方法

由于外部CT扫描图像中投影数据的严重缺失,传统的FBP算法与SART算法的重建图像均表现出严重的结构缺失状伪影。为了保留重建图像在切向方向与径向方向的边缘细节,首先采用WDTV最小化模型对模体进行CT重建,该模型是计算沿重建图像的切向方向和径向方向全变差的加权和,然后采用SART-WDTV算法来寻找具有最小WDTV的求解方案,在减小图像伪影的同时保留图像的边缘细节信息,以提高重建图像的质量^[23]。

3.1 局部方向梯度

对于单位向量α和离散图像f₍_i_,_j₎,其中(i,j)为离散图像中的像素点, $\nabla_{α}$ f₍_i_,_j₎为f₍_i_,_j₎沿着向量α方向的局部梯度向量,表达式为

\nabla_{α} f_{(i, j)} = [\nabla f_{(i, j)} \cdot α] α \approx D_{α} [f_{(i, j)}] α, (1)

式中:Ñf₍_i_,_j₎表示f₍_i_,_j₎的梯度向量;D_α[f₍_i_,_j₎]表示f₍_i_,_j₎沿着α方向的局部差分。采用线性插值的方法来计算离散的D_α[f₍_i_,_j₎]。

待重建图像如图4(a)所示,O点为重建图像的中心点,γ₍_i_,_j₎为O点和点(i,j)的连线与I轴之间的角度,点(i,j)的径向和切向用虚箭头来标出,即[cos γ₍_i_,_j₎,sin γ₍_i_,_j₎]为径向,{cos[γ₍_i_,_j₎+π/2],sin [γ₍_i_,_j₎+π/2]}为切向,这两个方向构成的局部正交坐标系,如图4(b)所示。对于给定的单位向量α=(cos α,sin α),则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = \frac{f_{α} - f_{(i, j)}}{ρ_{α}}, (2)

式中:f_α表示点(i,j)沿着α方向与网格的交点;ρ_α表示f_α与点(i,j)之间的距离。f_α的灰度值是采用线性插值法来计算的。计算过程中,给定点(i,j)和α=(cos α,sin α),因α代表不同的角度,则D_α[f₍_i_,_j₎]具有不同的表示形式。实验过程中,设置相邻像素之间的距离为1,即(i,j)和(i+1,j)之间的距离为1。对于0<α≤π/4,则D_α[f₍_i_,_j₎]可以表示为

\begin{array}{l} D_{α} [f_{(i, j)}] = \frac{f_{α} - f_{(i, j)}}{ρ_{α}} = \frac{[\tan α \cdot f_{(i + 1, j + 1)} + (1 - \tan α) \cdot f_{(i + 1, j)}] - f_{(i, j)}}{\frac{1}{\cos α}} = \\ (\cos α - \sin α) \cdot f_{(i + 1, j)} + \sin α \cdot f_{(i + 1, j + 1)} - \cos α \cdot f_{(i, j)} 。 (3) \end{array}

类似的计算过程可用于其他情况。

1) 对于0<α≤π/4,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\cos α - \sin α) \cdot f_{(i + 1, j)} + \sin α \cdot f_{(i + 1, j + 1)} - \cos α \cdot f_{(i, j)} 。 (4)

2) 对于π/4<α≤π/2,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\sin α - \cos α) \cdot f_{(i, j + 1)} + \cos α \cdot f_{(i + 1, j + 1)} - \sin α \cdot f_{(i, j)} 。 (5)

3) 对于π/2<α≤3π/4,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\sin α + \cos α) \cdot f_{(i, j + 1)} - \cos α \cdot f_{(i - 1, j + 1)} - \sin α \cdot f_{(i, j)} 。 (6)

4) 对于3π/4<α≤π,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (- \sin α - \cos α) \cdot f_{(i - 1, j)} + \sin α \cdot f_{(i - 1, j + 1)} + \cos α \cdot f_{(i, j)} 。 (7)

5) 对于π<α≤5π/4,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\sin α - \cos α) \cdot f_{(i - 1, j)} - \sin α \cdot f_{(i - 1, j - 1)} + \cos α \cdot f_{(i, j)} 。 (8)

6) 对于5π/4<α≤3π/2,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\cos α - \sin α) \cdot f_{(i, j - 1)} - \cos α \cdot f_{(i - 1, j - 1)} + \sin α \cdot f_{(i, j)} 。 (9)

7) 对于3π/2<α≤7π/4,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (- \cos α - \sin α) \cdot f_{(i, j - 1)} + \cos α \cdot f_{(i + 1, j - 1)} + \sin α \cdot f_{(i, j)} 。 (10)

8) 对于7π/4<α≤2π,则D_α[f₍_i_,_j₎]可以表示为

D_{α} [f_{(i, j)}] = (\cos α + \sin α) \cdot f_{(i + 1, j)} - \sin α \cdot f_{(i + 1, j - 1)} - \cos α \cdot f_{(i, j)} 。 (11)

设α₍_i_,_j₎=γ₍_i_,_j₎+π/2,那么

\begin{array}{l} S_{γ TV} = \sum_{(i, j)} | D_{γ_{(i, j)}} [f_{(i, j)}] |, (12) \\ S_{α TV} = \sum_{(i, j)} | D_{α_{(i, j)}} [f_{(i, j)}] |, (13) \\ S_{WDTV} = w_{1} S_{γ TV} + w_{2} S_{α TV}, (14) \end{array}

式中:S_γ_TV和S_α_TV分别表示沿径向和切向方向的局部差分和;S_WDTV表示S_γ_TV和S_α_TV的加权和; $D_{γ_{(i, j)}}$ 表示沿着径向方向的局部差分; $D_{α_{(i, j)}}$ 表示沿着切向方向的局部差分;ω₁和ω₂分别表示WDTV中径向和切向局部差分的权值。 $D_{γ_{(i, j)}}$ [f₍_i_,_j₎]和 $D_{α_{(i, j)}}$ [f₍_i_,_j₎]可以通过(4)~(11)式来计算求得。基于WDTV最小化的重建模型可以表示为

\min S_{WDTV} s.t. ‖ AS - p ‖_{2} < ε, S \geq 0, (15)

式中:A表示当前扫描模式下对应的系统矩阵;p表示实测投影数据矩阵;ε表示由噪声、散射和射束硬化等引起数据不一致的上限。综上所述,外部问题的图像重建可以转换为(15)式最优化问题。

3.2 SART-WDTV重建算法

采用SART-WDTV算法来求解(15)式的最优化问题。每次迭代的过程中,分别执行SART迭代重建步骤和WDTV最小化步骤以寻找最优的求解结果。SART算法更新重建图像后,采用非负约束来修正重建结果。WDTV最小化步骤中,采用梯度下降法来实现WDTV最小化的过程,该过程中最重要的是计算 $D_{γ_{(i, j)}}$ 和 $D_{α_{(i, j)}}$ 关于f₍_i_,_j₎的偏导数。设当前遍历的是f₍_i_,_j₎,并且0<γ_i_,_j≤π/4则π/2<α_i_,_j≤3π/4,根据(4)~(11)式可以得到

\begin{array}{l} \{\begin{array}{l} \frac{\partial | D_{γ_{(i, j)}} [f_{(i, j)}] |}{\partial f_{(i, j)}} \approx \frac{- D_{γ_{(i, j)}} [f_{(i, j)}] \cdot \cos γ_{(i, j)}}{\sqrt[]{{D_{γ_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \\ \frac{\partial | D_{γ_{(i, j)}} [f_{(i, j)}] |}{\partial f_{(i + 1, j)}} \approx \frac{D_{γ_{(i, j)}} [f_{(i, j)}] \cdot [\cos γ_{(i, j)} - \sin γ_{(i, j)}]}{\sqrt[]{{D_{γ_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \\ \frac{\partial |D_{γ_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i + 1, j + 1)}} \approx \frac{D_{γ_{(i, j)}} [f_{(i, j)}] \cdot \sin γ_{(i, j)}}{\sqrt[]{{D_{γ_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \end{array}, (16) \\ \{\begin{array}{l} \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i, j)}} \approx \frac{- D_{α_{(i, j)}} [f_{(i, j)}] \cdot \sin α_{(i, j)}}{\sqrt[]{{D_{α_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \\ \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i, j + 1)}} \approx \frac{D_{α_{(i, j)}} [f_{(i, j)}] \cdot [\sin α_{(i, j)} + \cos α_{(i, j)}]}{\sqrt[]{{D_{α_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \\ \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i - 1, j + 1)}} \approx \frac{- D_{α_{(i, j)}} [f_{(i, j)}] \cdot \cos α_{(i, j)}}{\sqrt[]{{D_{α_{(i, j)}} [f_{(i, j)} {]}}^{2} + ξ}} \end{array}, (17) \end{array}

式中:ξ表示一个小的正数,此处设置ξ=1×10^-8。设v¹和v²是与待重建图像f具有相同尺寸的矩阵,分量 $v_{(i, j)}^{1}$ 是 $D_{γ_{(i, j)}}$ 关于f₍_i_,_j₎的偏导数,分量 $v_{(i, j)}^{2}$ 是 $D_{α_{(i, j)}}$ 关于f₍_i_,_j₎的偏导数。对v¹和v²的分量进行更新,表达式为

\begin{array}{l} \{\begin{array}{l} v_{(i, j) update}^{1} = v_{(i, j)}^{1} + \frac{\partial |D_{γ_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i, j)}} \\ v_{(i + 1, j) update}^{1} = v_{(i + 1, j)}^{1} + \frac{\partial |D_{γ_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i + 1, j)}} \\ v_{(i + 1, j + 1) update}^{1} = v_{(i + 1, j + 1)}^{1} + \frac{\partial |D_{γ_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i + 1, j + 1)}} \end{array}, (18) \\ \{\begin{array}{l} v_{(i, j) update}^{2} = v_{(i, j)}^{2} + \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i, j)}} \\ v_{(i, j + 1) update}^{2} = v_{(i, j + 1)}^{2} + \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i, j + 1)}} \\ v_{(i - 1, j + 1) update}^{2} = v_{(i - 1, j + 1)}^{2} + \frac{\partial |D_{α_{(i, j)}} [f_{(i, j)}]|}{\partial f_{(i - 1, j + 1)}} \end{array} 。 (19) \end{array}

对于γ₍_i_,_j₎的其他角度,执行类似(16)~(19)式的计算过程。在遍历所有f₍_i_,_j₎之后获得偏导数 $v_{(i, j)}^{1}$ 和 $v_{(i, j)}^{2}$ ,然后采用梯度下降算法来更新重建图像。

重建图像的过程中,梯度下降算法的步长为β,β值太大会破坏数据的一致性,则重建图像变得均匀,β值太小会削弱S_WDTV,则重建图像中的结构缺失状伪影不能得到有效校正,因此同时调整WDTV最小化步骤中的迭代总数N_WDTV和β,便可以获得更好的效果。

图 4. 局部方向差分的计算方法^[23]。(a)局部方向差分的坐标系;(b)局部方向差分的计算方法

Fig. 4. Calculation method of local directional difference ^[23]. (a) Coordinate system of local directional difference; (b) calculation method of local directional difference

下载图片查看所有图片

3.3 改进的SART-WDTV算法

采用SART-WDTV重建算法来计算沿径向和切向方向的局部方向差分并作加权和,在重建模型的WDTV项中引入两个权重参数以控制径向边缘和切向边缘的不同响应强度,进而抑制重建图像中沿着径向方向的图像伪影,从而提升重建图像的质量。SART-WDTV算法在每次WDTV过程中仅计算并实现一次待重建像素沿径向和切向方向差分的最小化,因此提升重建效果有限。实验过程中,对沿着重建图像的径向和切向两个方向且在[-ϕ,ϕ]角度范围内分别施加不同的全变差(TV)最小化约束,以期达到进一步提升重建图像质量的目的,其中ϕ为计算局部方向差分选择的角度。

WDTV重建模型的过程中,待重建像素局部方向差分选用的方向为径向(γ角度方向)与切向(α角度方向)。改进的SART-WDTV算法的计算方式如图5所示。为了加强TV最小化作用的程度,WDTV最小化过程中待重建像素局部方向差分选用的方向为γ'₍_i_,_j₎和α'₍_i_,_j₎,可以表示为

\begin{array}{l} γ'_{(i, j)} = γ_{(i, j)} \pm (ϕ - n Δ ϕ), (20) \\ α'_{(i, j)} = α_{(i, j)} \pm (ϕ - n Δ ϕ), (21) \end{array}

式中:Δϕ表示角度间隔;n表示常数,即n=0,1,2,…,(ϕ/Δϕ)。

改进的SART-WDTV算法在每次WDTV过程中选用的方向为γ'₍_i_,_j₎和α'₍_i_,_j₎。根据(20)~(21)式,由于γ'₍_i_,_j₎和α'₍_i_,_j₎选用的角度数量为2(ϕ/Δϕ)+1,改进的SART-WDTV算法在每次WDTV最小化过程中实现2(ϕ/Δϕ)+1次的TV最小化。因此改进的SART-WDTV算法中TV最小化的作用得以加强,具有更明显的降噪与抑制外部CT图像径向边缘伪影的优点。

图 5. 改进方法中局部方向差分的计算方法。(a)局部方向差分的坐标系;(b)局部方向差分的计算方法

Fig. 5. Calculation method of local directional difference in improved method. (a) Coordinate system of local directional difference; (b) calculation method of local directional difference

下载图片查看所有图片

4 实验验证

4.1 模拟数据的实验验证

采用设计的模拟模体来验证改进WDTV正则化模型的有效性。外部CT扫描条件的设置如表1所示,采用SART、SART-TV和SART-WDTV算法对图2(c)外部CT扫描数据进行CT重建。设置SART、SART-TV和SART-WDTV算法的λ=0.05,SART-TV算法β=0.2、TV最小化步骤中的迭代总数N_TV=20,SART-WDTV算法的β=0.18、N_WDTV=80、w₁=1.00以及w₂=0.15,设置三种算法的总迭代次数为500次,得到的重建结果如图6所示。

从图6可以看到,采用SART算法重建的图像中存在明显的结构缺失条状伪影与图像噪声,所以图像的信噪比较低;采用SART-TV算法能够明显减少结构缺失条状伪影,但重建图像中的圆孔形状与方孔形状沿着切向方向拉长,所以出现严重失真的现象;采用SART-WDTV算法能够明显降低重建图像的噪声,所以图像失真的现象得到明显改善。对于宽度为128 pixel的正弦图像,设置SART-WDTV算法的β=0.05、N_WDTV=10、w₁=1.00、w₂=0、ϕ=15°和Δϕ=2.5°,500次迭代得到的重建结果如图7(b)所示。SART-WDTV算法与改进的SART-WDTV算法的重建结果如图7所示。

图 6. 不同算法的重建结果。(a) SART算法;(b) SART-TV算法;(c) SART-WDTV算法

Fig. 6. Reconstruction results of different algorithms. (a) SART algorithm; (b) SART-TV algorithm; (c) SART-WDTV algorithm

下载图片查看所有图片

图 7. 不同算法重建结果的对比。(a) SART-WDTV算法;(b)改进的SART-WDTV算法

Fig. 7. Comparison of reconstruction results of different algorithms. (a) SART-WDTV algorithm; (b) improved SART-WDTV algorithm

下载图片查看所有图片

从图6和图7可以看到,ROI-7的横条特征较为模糊,主要原因是外部CT扫描仅具有与ROI-7横条方向垂直的投影数据[图2(b)],缺失与横条方向一致的投影数据;射线垂直穿过横条特征后衰减很小,因此投影数据具有微弱的横条特征信息。迭代重建的过程中,微弱的横条特征投影信息再反投影到该条射线穿过所有重建的像素中,因此横条特征出现模糊的现象;对于ROI-1的横条特征情况正好相反,射线平行穿过ROI-1的横条特征后衰减较大,因此投影数据具有较强的特征信息,最终使ROI-1区域的重建效果更清晰。对比图7的重建结果,SART-WDTV算法仅计算并实现待重建像素沿径向和切向方向差分的最小化,导致重建图像中的细节(圆孔、方孔和线状结构边缘等)出现沿着径向或切向的失真现象。主要表现为方孔夹角边缘不够锐利,圆孔沿着径向方向的边缘不够圆滑,三个平行分布的线状结构发生向心倾斜。在改进的SART-WDTV重建模型中,SART-WDTV最小化的方向不再是径向与切向两个方向,而是径向和切向两个方向附近的角度方向,因此改进的SART-WDTV算法中TV最小化的作用更明显,具有更强的抑制图像噪声与恢复图像细节的能力。从图7(b)可以看到,改进的SART-WDTV算法重建图像的细节边缘更清晰,图像失真情况可以得到更好的改善,重建图像具有更高的质量。

采用均方误差(MSE)与峰值信噪比(PSNR)指标^[24-25]对SART、SART-TV、SART-WDTV和改进的SART-WDTV算法重建的图像进行定量评价。MSE是衡量两张图像之间的平均像素差距,MSE值越小,表示重建图像的像素灰度差距与真实图像越小,重建质量越高。PSNR值越大则表示图像的失真现象越少。得到的MSE和PSNR指标如表2所示。从表2可以看到,SART、SART-TV、SART-WDTV与改进的SART-WDTV 4种算法中,改进的SART-WDTV算法的MSE值最小且PSNR值最大,表明改进的SART-WDTV重建算法重建图像的效果最优。

表 2. 不同算法重建图像的MSE和PSNR

Table 2. MSE and PSNR of images reconstructed by different algorithms

Algorithm	MSE	PSNR
SART	1319.476	16.927
SART-TV	1034.756	17.982
SART-WDTV	874.630	18.713
Improved SART-WDTV	828.189	18.950

查看所有表

4.2 冷中子数据的实验验证

基于中国工程物理研究院的CMRR (China Mianyang Research Reactor)反应堆,得到了钢材质齿轮样件的冷中子层析扫描数据,获得的钢材质齿轮的直径约为150 mm,高度为50 mm。CT扫描序列投影图像的尺寸为512 pixel×1504 pixel,投影幅数为1200幅,CT扫描角度为360°。对第360层齿轮样品投影的正弦图像左右两端冗余部分进行裁切,并进行2 pixel×2 pixel合并采样处理。获得的第360层取对数操作后投影的正弦图像如图8(a)所示,正弦图像的尺寸为654 pixel×1200 pixel。采用SART-TV算法得到的重建图像如图8(b)所示。

为了验证不同算法的CT重建图像质量,对图8(a)的投影正弦图像进行裁切,获得图8(a)中第1 pixel至第218 pixel的投影正弦图像,如图9(c)所示。将裁切后的投影数据作为外部CT扫描投影的数据,此时投影正弦图像的尺寸为218 pixel×1200 pixel,仅为完整投影正弦图像的1/3。采用FBP、SART、SART-TV以及SART-WDTV算法对图9(c)的投影正弦图像进行CT重建。当采用FBP算法重建图像时,需要根据平行束数据的对称性原理将正弦图像右侧的投影数据进行补充,补充后的投影数据如图9(a)所示,SART算法重建参数λ=0.05。由于中子投影数据的信噪比比X射线图像低,为了保证迭代重建算法的收敛性,设置的较大迭代次数N_iter=1000。FBP和SART算法重建的图像如图9(b)和图9(d)所示。

图 8. 齿轮样品完整投影正弦图与重建图像。(a)完整投影正弦图;(b)重建图像

Fig. 8. Complete projected sinogram and reconstructed image of gear sample. (a) Complete projected sinogram; (b) reconstructed image

下载图片查看所有图片

图 9. 不同算法对齿轮样品的重建图像。(a)对称补数据后正弦图;(b)FBP算法;(c)外部CT扫描正弦图;(d) SART算法

Fig. 9. Reconstruction images of gear sample by different algorithms. (a) Sinogram after symmetrical complement data; (b) FBP algorithm; (c) sinogram obtained by exterior CT scan; (d) SART algorithm

下载图片查看所有图片

设置SART-TV算法参数β=0.2、N_TV=30和N_iter=1000,设置SART-WDTV算法参数β=0.1、N_WDTV=120、w₁=0.75、w₂=1.00和N_iter=1000。SART-TV算法与SART-WDTV算法的外部CT重建结果分别如图10(a)和图10(b)所示。采用改进的SART-WDTV算法对模体进行外部CT重建,设置β=0.05、N_WDTV=16、w₁=0.10、w₂=1.00、ϕ=20°、Δϕ=5°和N_iter=1000。改进的SART-WDTV算法的重建结果如图10(c)所示。对比图9(b)、图9(d)和图10,可以得到FBP算法与SART算法外部CT重建图像的细节十分模糊,严重影响重建图像的识别与判读;SART-TV、SART-WDTV与改进的SART-WDTV算法重建图像的质量得到有效提高,重建图像的细节得到一定的改善,重建图像能够较好地反映样品的真实结构。重建图像的局部放大结果如图11所示。从图11可以看到,相比于SART-TV算法,SART-WDTV算法与改进的SART-WDTV算法中齿轮内圈轮廓的圆度更大,而且ROI-Ⅰ和ROI-Ⅲ中槽的轮廓更符合真实情况。另外观察图11中箭头标记的齿根部位,发现改进的SART-WDTV算法中的齿根部位更清晰,与真实情况较为接近。总而言之,改进的SART-WDTV算法重建图像的噪声明显减少,重建图像的细节较为清晰,重建图像的质量更优,因此验证所提算法的有效性。另外,对不同迭代重建算法的收敛速度进行分析,利用方均根误差(RMSE)指标来评价重建算法的收敛速度^[26]。图12为采用SART-TV、SART-WDTV与改进的SART-WDTV算法的重建图像与图8(b)理想图像的RMSE随迭代次数的变化曲线。需要注意的是,当计算RMSE时,首先将重建图像与理想图像进行归一化,然后再计算两者之间的RMSE。从图12可以看到,改进的SART-WDTV重建算法经过少量的迭代次数后,重建图像就足以达到很小的RMSE,由RMSE曲线变化规律可以看到,所提算法在迭代500次就能够得到高质量的CT图像,表明改进的SART-WDTV算法比SART-TV算法和SART-WDTV算法具有更快的收敛速度。

图 10. 齿轮样品不同算法重建图像。(a) SART-TV算法;(b) SART-WDTV算法;(c)改进的SART-WDTV算法

Fig. 10. Gear sample reconstruction image with different algorithms. (a) SART-TV algorithm; (b) SART-WDTV algorithm; (c) improved SART-WDTV algorithm

下载图片查看所有图片

图 11. 不同算法重建图像的局部放大结果。(a)完整数据重建;(b) SART-TV算法;(c) SART-WDTV算法;(d)改进的SART-WDTV算法

Fig. 11. Local magnification results of reconstructed images with different algorithms. (a) Complete data reconstruction; (b) SART-TV algorithm; (c) SART-WDTV algorithm; (d) improved SART-WDTV algorithm

下载图片查看所有图片

图 12. 不同重建算法得到的RMSE曲线

Fig. 12. RMSE curves obtained by different reconstruction algorithms

下载图片查看所有图片

5 结论

外部CT是一种常用的大尺寸筒状构件层析检测方法。当外部CT扫描时,射线仅仅穿过待检样品的外部轮廓,造成投影数据严重缺失,不完全的投影数据会导致重建图像出现严重的径向边缘伪影。针对外部CT图像含有伪影的特点,采用WDTV算法计算沿径向和切向方向的局部方向差分并作加权和,引入两个权重参数以控制径向边缘和切向边缘的不同响应强度,可以有效地抑制重建图像的伪影。然而WDTV算法仅计算并实现待重建像素沿径向和切向方向差分的最小化,导致重建图像的细节出现沿着径向或切向的失真现象。为了进一步提高高噪声中子外部CT重建图像的质量,在WDTV算法重建模型的框架下,对沿着重建图像的径向和切向两个方向附近的一定角度范围内分别施加不同的TV最小化约束。改进的WDTV算法中TV最小化的作用更明显,而且具有更强地抑制径向边缘伪影与抗噪声的性能。计算机仿真结果和齿轮冷中子实验结果表明,改进的WDTV最小化重建算法能够更好地抑制重建图像的噪声,减小径向边缘的伪影,提高重建图像的质量。但是改进的WDTV最小化重建算法所需调整参数过多,这限制该算法在中子外部CT领域的推广和应用,因此未来的研究方向主要集中于重建参数的优化与自适应调整方面的研究。

参考文献

[1] Endrizzi M. X-ray phase-contrast imaging[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2018, 878: 88-98.

[2] 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.

[3] 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.

[4] 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.

[5] Bayanov B F, Belov V P, Bender E D, et al. Accelerator-based neutron source for the neutron-capture and fast neutron therapy at hospital[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 1998, 413(2/3): 397-426.

[6] Sinha V, Avachat A V, Lee H K. Design and development of a neutron/X-ray combined computed tomography system at Missouri S & T[J]. Journal of Radioanalytical and Nuclear Chemistry, 2013, 296(2): 799-806.

[7] 郭雨濛, 曾理. CT外问题的原理与算法综述[J]. 中国体视学与图像分析, 2017, 22(3): 320-326.

Guo Y M, Zeng L. Review on the principle and algorithms of the exterior problem of CT[J]. Chinese Journal of Stereology and Image Analysis, 2017, 22(3): 320-326.

[8] 倪风岳. CUDA加速CV图像分割和外部CT图像重建算法研究[D]. 重庆: 重庆大学, 2011: 45- 55.

Ni FY. Algorithm study on accelerate CV image segmentation and exterior industrial image reconstruction by CUDA[D]. Chongqing: Chongqing University, 2011: 45- 55.

[9] 陈庆贵, 卢洪义, 齐强, 等. 固体火箭发动机界面脱粘切向CT检测[J]. 固体火箭技术, 2016, 39(3): 347-352.

Chen Q G, Lu H Y, Qi Q, et al. Tangential CT inspection of interface debonding of SRM[J]. Journal of Solid Rocket Technology, 2016, 39(3): 347-352.

[10] Guo J Q, Zeng L, Liu B D. High-quality image reconstruction from exterior helical cone-beam CT data for NDE of industrial pipelines[J]. Insight, 2011, 53(10): 534-541.

[11] Quinto E T. Local algorithms in exterior tomography[J]. Journal of Computational and Applied Mathematics, 2007, 199(1): 141-148.

[12] Cormack A M. Representation of a function by its line integrals, with some radiological applications[J]. Journal of Applied Physics, 1963, 34(10): 2722-2727.

[13] Quinto E T. Tomographic reconstructions from incomplete data-numerical inversion of the exterior Radon transform[J]. Inverse Problems, 1988, 4(3): 867-876.

[14] Quinto E T. Singularities of the X-ray transform and limited data tomography in R2 and R3[J]. SIAM Journal on Mathematical Analysis, 1993, 24(5): 1215-1225.

[15] Quinto E T. Exterior and limited-angle tomography in non-destructive evaluation[J]. Inverse Problems, 1998, 14(2): 339-353.

[16] Quinto E T. Singular value decompositions and inversion methods for the exterior Radon transform and a spherical transform[J]. Journal of Mathematical Analysis and Applications, 1983, 95(2): 437-448.

[17] Vainberg E I, Kazak I A, Kurozaev V P. Reconstruction of the internal three-dimensional structure of objects based on real-time integral projections[J]. Soviet Journal of Nondestructive Testing, 1981, 17(6): 415-423.

[18] Sidky E Y, Kao C, 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.

[19] 张海娇, 孔慧华, 孙永刚. 基于结构先验的加权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.

[20] 刘进, 亢艳芹, 顾云波, 等. 稀疏张量约束的低剂量CT图像重建[J]. 光学学报, 2019, 39(8): 0811004.

Liu J, Kang Y Q, Gu Y B, et al. Low dose computed tomography image reconstruction based on sparse tensor constraint[J]. Acta Optica Sinica, 2019, 39(8): 0811004.

[21] Zeng L, Liu B D, Liu L H, et al. A new iterative reconstruction algorithm for 2D exterior fan-beam CT[J]. Journal of X-Ray Science and Technology, 2010, 18(3): 267-277.

[22] Guo Y M, Zeng L. Improved iterative image reconstruction algorithm for the exterior problem of computed tomography[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2017, 842: 96-108.

[23] Guo Y M, Zeng L, Wang C X, et al. Image reconstruction model for the exterior problem of computed tomography based on weighted directional total variation[J]. Applied Mathematical Modelling, 2017, 52: 358-377.

[24] Huynh-Thu Q, Ghanbari M. The accuracy of PSNR in predicting video quality for different video scenes and frame rates[J]. Telecommunication Systems, 2012, 49(1): 35-48.

[25] 林强, 杨民, 唐彬, 等. 基于加权总差分最小化的中子稀疏投影计算机断层重建方法[J]. 光学学报, 2019, 39(7): 0711003.

Lin Q, Yang M, Tang B, et al. Neutron computed tomography reconstruction method using sparse projections based on weighted total difference minimization[J]. Acta Optica Sinica, 2019, 39(7): 0711003.

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

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

林强, 杨民, 张晓敏, 唐彬, 刘斌, 霍合勇. 一种改进的中子外部CT图像伪影抑制方法[J]. 光学学报, 2020, 40(22): 2210001. Qiang Lin, Min Yang, Xiaomin Zhang, Bin Tang, Bin Liu, Heyong Huo. Improved Image Artifacts Suppression Method for Neutron External CT[J]. Acta Optica Sinica, 2020, 40(22): 2210001.

一种改进的中子外部CT图像伪影抑制方法 下载： 1029次封面文章

1 引言

2 传统的解析重建与迭代重建算法

图 1. 中子外部CT扫描原理图

Fig. 1. Schematic of neutron exterior CT scan

表 1. 外部CT扫描参数

Table 1. Parameters of external CT scan phantom

图 2. 外部CT扫描模体的示意图。(a)模体;(b)扫描过程;(c)正弦图像

Fig. 2. Schematic of external CT scan phantom. (a) Phantom; (b) scanning process; (c) sinusoidal image

图 3. 不同算法重建的图像。(a) FBP算法;(b) SART算法

Fig. 3. Images reconstructed by different algorithms. (a) FBP algorithm; (b) SART algorithm

3 WDTV最小化重建方法

3.1 局部方向梯度

3.2 SART-WDTV重建算法

图 4. 局部方向差分的计算方法[23]。(a)局部方向差分的坐标系;(b)局部方向差分的计算方法

Fig. 4. Calculation method of local directional difference [23]. (a) Coordinate system of local directional difference; (b) calculation method of local directional difference

3.3 改进的SART-WDTV算法

图 5. 改进方法中局部方向差分的计算方法。(a)局部方向差分的坐标系;(b)局部方向差分的计算方法

Fig. 5. Calculation method of local directional difference in improved method. (a) Coordinate system of local directional difference; (b) calculation method of local directional difference

4 实验验证

4.1 模拟数据的实验验证

图 6. 不同算法的重建结果。(a) SART算法;(b) SART-TV算法;(c) SART-WDTV算法

Fig. 6. Reconstruction results of different algorithms. (a) SART algorithm; (b) SART-TV algorithm; (c) SART-WDTV algorithm

图 7. 不同算法重建结果的对比。(a) SART-WDTV算法;(b)改进的SART-WDTV算法

Fig. 7. Comparison of reconstruction results of different algorithms. (a) SART-WDTV algorithm; (b) improved SART-WDTV algorithm

表 2. 不同算法重建图像的MSE和PSNR

Table 2. MSE and PSNR of images reconstructed by different algorithms

4.2 冷中子数据的实验验证

图 8. 齿轮样品完整投影正弦图与重建图像。(a)完整投影正弦图;(b)重建图像

Fig. 8. Complete projected sinogram and reconstructed image of gear sample. (a) Complete projected sinogram; (b) reconstructed image

图 9. 不同算法对齿轮样品的重建图像。(a)对称补数据后正弦图;(b)FBP算法;(c)外部CT扫描正弦图;(d) SART算法

Fig. 9. Reconstruction images of gear sample by different algorithms. (a) Sinogram after symmetrical complement data; (b) FBP algorithm; (c) sinogram obtained by exterior CT scan; (d) SART algorithm

图 10. 齿轮样品不同算法重建图像。(a) SART-TV算法;(b) SART-WDTV算法;(c)改进的SART-WDTV算法

Fig. 10. Gear sample reconstruction image with different algorithms. (a) SART-TV algorithm; (b) SART-WDTV algorithm; (c) improved SART-WDTV algorithm

图 11. 不同算法重建图像的局部放大结果。(a)完整数据重建;(b) SART-TV算法;(c) SART-WDTV算法;(d)改进的SART-WDTV算法

Fig. 11. Local magnification results of reconstructed images with different algorithms. (a) Complete data reconstruction; (b) SART-TV algorithm; (c) SART-WDTV algorithm; (d) improved SART-WDTV algorithm

图 12. 不同重建算法得到的RMSE曲线

Fig. 12. RMSE curves obtained by different reconstruction algorithms

5 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

一种改进的中子外部CT图像伪影抑制方法下载： 1029次封面文章

图 4. 局部方向差分的计算方法^[23]。(a)局部方向差分的坐标系;(b)局部方向差分的计算方法

Fig. 4. Calculation method of local directional difference ^[23]. (a) Coordinate system of local directional difference; (b) calculation method of local directional difference