基于加权总差分最小化的中子稀疏投影计算机断层重建方法 下载: 840次
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迭代重建算法
在迭代重建算法中,待重建图像被离散化为
式中:
式中:
针对Shepp-Logan模型[20],设计了平行束CT模拟程序,在0°~180°范围内均匀采集Shepp-Logan模型360幅投影,验证ART、SART及SIRT的重建效果。ART、SART与SIRT的投影数据访问顺序默认为顺序访问,即按照CT扫描的实际顺序,从第一个投影顺序访问至最后一个投影。当松弛因子
图 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的数学表达式分别为
式中:
表 1. Shepp-Logan重建图像的MSE和PSNR
Table 1. MSE and PSNR of reconstructed images of Shepp-Logan
|
3 稀疏角度投影SIRT图像重建
采用稀疏角度投影CT扫描时,采集得到的投影数据量不能满足Shannon/Nyquist采样定理,滤波反投影重建算法(FBP)的重建图像会出现欠采样条状伪影。随着投影数据量的减少,重建图像中伪影情况大幅度加剧。
图 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正则化方法在中子稀疏角度投影重建中的应用。
投影数据迭代重建问题,可以等价于优化问题的求解过程,即
式中:
加权总差分要实现对图像梯度的稀疏性和连续性的双重约束。对于梯度的稀疏性,约束任意一个像素点与其邻域内像素点像素保持一致;对于梯度的连续性,则需要约束任意一个像素位置与邻域内其他像素点的局部梯度保持一致[22]。
图 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
对局部梯度稀疏性和局部梯度连续性的约束差分可分别表示为
其中
则对局部梯度连续性约束差分可简化为
加权求和得到加权总差分,实现对梯度连续性和稀疏性的双重约束:
式中:
当
式中:
式中:
式中:
余维等[17-18]针对模拟模体的稀疏角度投影数据,采用交替最小化的方式实现对(13)式的求解:首先利用SART进行重建以缩小原始投影数据与模拟投影数据之间的差异,在每次SART重建后进行一次软阈值滤波,从而减小重建图像的加权总差分,最终得到理想的重建图像,因此其算法中一次SART-WTDM算法主循环主要包括一次SART重建与一次WTDM过程。本研究针对中子投影数据噪声高、图像质量差的特点,提出了SIRT-WTDM重建方法。SIRT-WTDM重建方法主要对余维等[17-18]提出的SART-WTDM方法进行两个方面的调整:1)将SIRT作为基于WTDM重建方法的保真项,以增强算法的抗噪声性能;2)设置SIRT-WTDM算法主循环中WTDM的次数为
5 实验验证
5.1 模拟数据的实验验证
针对Shepp-Logan模型,本研究在0°~180°范围内均匀采集90幅投影,比较FBP、SIRT、SART-WTDM及SIRT-WTDM算法在稀疏角度投影时的CT重建图像质量,Shepp-Logan模型的CT条件如
表 2. 基于Shepp-Logan模拟的CT条件
Table 2. Simulative CT conditions based on Shepp-Logan
|
SIRT的松弛因子
图 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算法的
图 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
由
为了更好地对
表 3. Shepp-Logan的90幅投影重建图像的MSE和PSNR
Table 3. MSE and PSNR of reconstructed images of Shepp-Logan with 90 projections
|
5.2 冷中子数据的实验验证
基于中国工程物理研究院的CMRR(China Mianyang Research Reactor)反应堆,获得了编码管样品冷中子层析扫描数据。中子CT扫描条件如
在上述采集条件下,对投影正弦图左右两端冗余部分进行裁切并进行均匀采样,获得75幅投影正弦图,对75幅稀疏投影条件下的编码管样品分别进行FBP、SIRT、SART-WTDM及SIRT-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
|
图 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算法的
图 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
将
表 5. 编码管样品重建图像的MSE和PSNR
Table 5. MSE and PSNR of reconstructed images of code tube sample
|
为了进一步比较FBP、SIRT、SART-WTDM与SIRT-WTDM算法对稀疏角度中子投影数据的CT重建效果,采用FBP、SIRT、SART-WTDM及SIRT-WTDM算法分别对
本实验中设置SIRT-WTDM算法的
对比
图 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领域的应用,因此未来的研究方向主要集中于重建算法加速与重建参数的自适应调整方面的研究。
[4] 魏国海, 韩松柏, 陈东风, 等. 中子照相技术在核燃料元件无损检测中的应用[J]. 核技术, 2012, 35(11): 821-826.
[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.
[11] 王林元, 刘宏奎, 李磊, 等. 基于稀疏优化的计算机断层成像图像不完全角度重建综述[J]. 物理学报, 2014, 63(20): 208702.
[12] 蔺鲁萍, 王永革. 不完全角度CT图像重建的模型与算法[J]. 北京航空航天大学学报, 2017, 43(4): 823-830.
[16] 张海娇, 孔慧华, 孙永刚. 基于结构先验的加权NLTV能谱CT重建算法[J]. 光学学报, 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.
[19] 闫镔, 李磊. CT图像重建算法[M]. 北京: 科学出版社, 2014: 5- 14.
YanB, LiL. CT image reconstruction algorithm[M]. Beijing: Science Press, 2014: 5- 14.
[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.
Article Outline
林强, 杨民, 唐彬, 刘斌, 霍合勇, 刘家伟. 基于加权总差分最小化的中子稀疏投影计算机断层重建方法[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.