利用CAD模型的不完全扫描CT图像重建 
下载： 763次
1 引言
工业上,对大尺寸和高密度工件进行X射线计算机断层成像(CT)时,会出现穿不透、扫描角度有限等情况,导致扫描得到的投影数据不完备[1]。例如工业上对在役管道进行检测时,由于受检测环境(如靠附于墙壁)的限制,只能以有限的旋转角度进行旋转扫描检测。此时采用基于CT图像重建的解析重建算法将不能够精确地重建出被扫描物体的断层面信息,且重建图像中会有明显的滑坡伪影[2]。对于有限角CT重建中出现的滑坡伪影,Frikel等[3]利用微局部分析进行了说明:如果图像特征的边界与投影X射线相切,则该边界较容易从有限角投影数据中重建出来,反之,该边界较难重建出来。
在投影数据有缺失的情况下,迭代算法可以获得比解析算法更高质量的图像。迭代算法包括代数重建算法(ART)[4]、联合代数重建算法(SART)[5]、最大似然期望最大化算法(MLEM)[6]、基于有序子集的最大似然期望最大化算法(OSEM)[7],但是迭代算法的速度慢,且是近似重建。正则化方法[8]可利用不同的先验信息约束重建模型的解,来加强重建图像与先验信息的一致性。由于大多数CT图像在梯度域是稀疏或近似稀疏的,因此基于全变差(TV)[9]、各向异性TV(ATV)[10]、广义TV(TGV)[11]、分数阶TV[12]、自适应权重TV[13]、乘性[14]正则化及图像梯度L0范数[15]正则化的有限角CT重建模型能够重建更高质量的图像。
基于图像域先验的方法是解决有限角度问题的主流手段,利用先验信息构建约束条件,可以有效改善有限角度问题的病态性。如将形状或结构相似的先验图像作为约束,应用先验图像约束压缩感知(PICCS)[16]来提升不完全角度下的图像重建质量;利用结构相似的先验图像构建字典[17]来减少有限角度伪影;用TV正则化与迭代边缘检测相结合的边缘引导TV最小化方法(EGTVM)[18]重建优质图像;基于小波紧框架和先验图像[19]、基于先验图像诱导的非局部均值[20-21]的有限角CT重建方法也是对重建图像进行约束,进一步提高了图像的质量。最近几年流行的基于深度学习的重建方法[22-23]也有被应用在有限角投影数据的CT图像重建当中。
上述介绍的基于图像域先验的方法主要针对医学CT图像重建问题,它们的CT先验图像容易获得。而工业CT的有限角是CT系统本身的局限性造成的,其重建算法的CT先验图像无法获取。为了提高不完全扫描CT图像重建的质量,本文提出一种解决方法:将工件的CAD设计模型作为先验知识,从中提取需要的先验图像,通过提出的优化迭代重建算法(SART+TVM+PRIOR)来实现该类对象的有限角CT重建。实验结果表明,所提方法能有效减少或消除边缘伪影,使得重建图像的边缘结构清晰可辨。
2 方法原理
2.1 有限角CT扫描模型
2.2 SART+TVM图像重建算法
基于图像TV的CT重建模型为
式中:图像f∈RN,N为图像大小;A为投影系数矩阵;gδ为探测器得到的投影数据的向量形式,gδ∈RM,M为射线总数。图像f的TV范数为
式中:fs,t为第s行、第t列的图像灰度值;τ为一个很小的扰动参数。
求解未知的CT图像f,即求解(1)式模型,结合投影到凸集方法和梯度下降方法的SART+TVM算法具体包括以下步骤。
1) 初始化重建参数,SART重建算法的松弛因子ωn和停止迭代次数Ncount,TVM步骤的松弛因子α和每轮的停止迭代次数NTV,初始图像f(0)为0。
2) 采用SART重建算法,
3) 对得到的中间结果
4) 构造TV最小化的梯度下降方向:
5) 转入TVM最小化的迭代步骤的子循环,l=1,2,…,NTV,l为停止迭代次数,得到
6) 如满足停止迭代的条件,则输出f(TV-GRAD),否则f(n)=f(TV-GRAD),n=n+1,转到步骤2)。
从SART+TVM重建算法迭代步骤中可以看出:步骤2)SART重建算法用来确保数据的一致性;步骤3)是为了让重建图像在约束条件f≥0内;步骤4)和5)迭代步骤使得目标函数‖f‖TV最小化(也即TVM)。
2.3 PICCS重建算法
PICCS算法最早应用于医学CT图像重建,其先验图像来自于前期的医学CT扫描,然后利用先验图像作为约束,来提升不完全角度下的图像重建质量,模型为
式中:fpri为先验图像;λ为正则化参数。针对(7)式,PICCS算法也是采用两步交替迭代(SART和TVM交替迭代)的方式进行求解,与SART+TVM算法不同的是步骤5):
也即是先验图像fpri在TV最小化的迭代步骤的子循环中,其他的重建参数和重建步骤与SART+TVM算法相同。
2.4 所提图像重建算法
为了突破工业CT有限角扫描无法获取先验图像的难点,本文提出一种基于CAD模型的有限角CT重建算法,主要步骤是:读取工件的CAD三维模型图,根据CT扫描断层位置获取先验图像fpri,再对扫描投影数据进行带先验图像的SART+TVM+PRIOR算法重建。
2.4.1 先验图像的获取
首先需要利用AutoCAD软件打开被扫描工件的CAD三维模型图M,根据实际扫描状态确定工件的基准面,建立三维坐标系。设基准面的分层位置为O,基准面上的坐标为X、Y轴坐标,则Z轴坐标代表垂直于基准面的分层位置P。
1) 由工件CT扫描的断层位置计算其CAD模型对应的分层位置P,利用AutoCAD软件自带的“平面截面”功能,即可得到位置P处的截面(剖切面),根据实际域到像素域的转换系数,将得到的截面转换成像素截面。
2) 对得到的像素截面进行灰度赋值,任一点像素的灰度值和该点物质的衰减系数相等。对于多材质工件,需要根据CT系统X射线能量和模型M中各部分材质来确定一张衰减系数表,具体的X射线能量和不同的材质对应一个衰减系数,然后再确定每一个像素的衰减系数,从而生成一幅衰减系数图像fs;对于单材质工件,由于工件材料单一,工件内衰减系数相同,可通过计算像素截面工件部分的像素总个数Nnum和有限角度的投影数据的平均值Aaverg(如有限角度[1°,120°]内有166个投影角度,所有的投影数据累加再除以166即为Aaverg),然后将Aaverg/Nnum值(衰减系数)赋给像素截面工件部分的每一个像素,从而生成衰减系数图像fs。
3) 在实际工件CT扫描后,对得到的原始投影数据进行SART+TVM算法的数次迭代重建,设重建图像大小为W×W,确定图像中工件的具体区域S;创建W×W大小的零灰度图像,将fs配准到该图像的区域S内,最后得到先验图像fpri。
2.4.2 SART+TVM+PRIOR重建算法
在得到先验图像fpri的情况下,在SART+TVM重建算法的步骤1)和2)中加入fpri,变成SART+TVM+PRIOR算法,具体步骤如下。
1) 初始化重建参数:SART重建算法的松弛因子ωn和停止迭代次数Ncount,TVM步骤的松弛因子α和每轮的停止迭代次数NTV,将初始图像f(0)设为先验图像fpri。
2) 采用(3)式的SART重建算法,
后续步骤与SART+TVM算法相同。
3 仿真与实际扫描实验
利用仿真截面和实际扫描工件的投影数据实验对所提方法的可行性进行论证。本实验的计算机配置为Windows 7 64位系统,NVIDIA Quadro K620显卡,2.8 GHz CPU,8.0 GB内存。对有限角CT重建图像采取主观和客观的质量评价方法,经常使用的几种量化指标有均方根误差(RMSE)、峰值性噪比(PSNR)和结构相似性(SSIM)指标。它们的表达式分别为
式中:f为被重建的图像;fref为参考图像;
3.1 仿真实验
3.1.1 先验图像的获取
现有一多材质工件的CAD设计模型,如
图 2. 工件CAD设计模型及剖切面示意图。(a) CAD模型;(b)剖切面
Fig. 2. CAD model of object and schematic of section. (a) CAD model; (b) section
第一步:由工件CT扫描的断层位置可得其CAD模型对应的分层位置P=40 mm,利用AutoCAD软件自带的“平面截面”功能得到位置P处的截面;根据实际域到像素域的转换系数,将得到的截面转换成像素截面,如
第二步:可根据CT系统X射线能量和CAD模型的各部分材质确定多材质工件每一像素的衰减系数,生成一幅衰减系数图像fs。如
图 3. 获取过程图像。(a)像素截面;(b)衰减系数图像;(c)先验图像
Fig. 3. Image of acquisition process. (a) Pixel section; (b) attenuation coefficient image; (c) prior image
第三步:创建512×512大小的零灰度图像,将衰减系数图像fs配准到该图像的中心区域,最后得到先验图像fpri,如
3.1.2 不完备投影数据
对工件在Z轴坐标40 mm的分层进行有限角度的CT扫描。实际的工件内部常常会有缺陷裂纹等损伤,在先验图像的基础上,利用大小不一的黑色孔洞来代表孔隙,长度不规则的黑色细线来代表裂纹,以含各类缺陷的仿真图作为实际工件在Z轴坐标40 mm的分层,如
图 4. 仿真图。(a)含缺陷;(b)含噪声和缺陷
Fig. 4. Simulated images. (a) Include defects; (b) include defects and noise
然后对该分层也即是含噪声和缺陷的仿真图先进行归一化处理(将0~255的灰度值线性映射到0~1),再进行有限角度120°(θ∈[1°,120°])的CT扫描,角度增量为1°,获得的不完备投影数据如
3.1.3 重建参数的选择
对工件模型的有限角度的不完备投影数据进行CT图像重建。大量文献选取的SART的松弛因子ωn=0.5和TVM步骤的松弛因子α=0.1对于SART+TVM算法有较好效果,接下来通过仿真实验来比较不同NTV对重建图像的影响。
在停止迭代次数Ncount=100的情况下,依次改变NTV值为10,20,30,40,50。加入先验图像fpri的SART+TVM+PRIOR算法和未加入fpri的SART+TVM算法的重建图像如
图 6. 不同NTV重建的图像。(a) 10;(b) 20;(c) 30;(d) 40
Fig. 6. Reconstructed images with different NTV. (a) 10; (b) 20; (c) 30; (d) 40
将含缺陷的仿真图[
表 1. 不同NTV重建图像的量化指标
Table 1. Quantitative index of reconstructed images with different NTV
|
图 7. SART+TVM+PRIOR算法不同NTV重建图像的PSNR和SSIM
Fig. 7. PSNR and SSIM of images with different NTVreconstructed by SART+TVM+PRIOR algorithm
3.1.4 与PICCS算法的比较
为公正地比较所提SART+TVM+PRIOR算法与PICCS算法的重建效果,将它们对不完备投影数据进行重建时的参数大小设置一致。PICCS的重建参数初始化:初始图像值f(0)=0,SART的松弛因子ωn=0.5,TVM步骤每轮的停止迭代次数NTV=30、松弛因子α=0.1,先验图像fpri相同,正则化参数λ=0.1,λ越小,先验信息越主要。
图 8. 不同迭代次数下的重建图像。(a) 100次;(b) 300次;(c) 500次
Fig. 8. Reconstructed images under different iterations. (a) 100; (b) 300; (c) 500
为比较SART+TVM+PRIOR、PICCS和SART+TVM三种算法的优劣,将
图 9. PSNR、SSIM与迭代次数的变化曲线。(a) PSNR;(b) SSIM
Fig. 9. PSNR and SSIM versus number of iterations. (a) PSNR; (b) SSIM
记图像a为SART+TVM+PRIOR算法迭代100次的重建图像,记图像b为PICCS算法迭代500次的重建图像。
表 2. 图像a和b的质量指标
Table 2. Quantitative index of image a and b
|
3.2 扫描工件实验
3.2.1 工件制作
利用AutoCAD软件绘制工件的三维立体模型,制作的两个尺寸大小相同的模型如
图 10. AutoCAD制作的三维模型。(a)模型1;(b)模型2
Fig. 10. 3D model made by AutoCAD software. (a) Model 1; (b) model 2
制作好的含各种缺陷的三维模型文件在AutoCAD软件中以.dwg格式保存,软件自带的格式转换功能可将其转换成.stl格式。利用MakerBot软件打开.stl格式的工件模型2,该软件工作界面如
图 11. 3D打印的软硬件及工件实物。(a)软件界面;(b)打印机;(c)工件
Fig. 11. Software and hardware of 3D printing and real object. (a) Software interface; (b) printer; (c) object
3.2.2 工件扫描
工件扫描实验在CD-130BX/UCT型号的微焦CT扫描仪上进行,如
表 3. 扫描参数
Table 3. Parameters of the scanning
|
3.2.3 先验图像的获取
先验图像的获取过程如
3.2.4 不完备投影数据重建
为了模拟有限角度重建,抽取[1°,150°]、[1°,120°]、[30°,150°]和[30°,120°]内的投影数据进行CT图像重建,重建算法的参数设置:ωn=0.5,NTV=10(考虑到噪声不大),α=0.1,重建图像大小为512×512。在实际的CT图像重建中参考图像很难获得或者没有参考图像,因此在工件扫描实验中没有图像的量化指标,使用重建图像来对比加入和未加入先验图像的重建方法性能。
在有限角度为[1°,150°]和[1°,120°]、投影角度数分别为208和166时,加入和未加入先验图像的SART和SART+TVM算法分别重建的图像如
图 14. 有限角度为[1°,150°]和[1°,120°]时重建的图像,显示窗口为[0,1]。(a) SART;(b) SART+PRIOR;(c) SART+TVM;(d) SART+TVM+PRIOR
Fig. 14. CT images reconstructed in the range of [1°,150°] and [1°,120°], display window is [0,1]. (a) SART; (b) SART+PRIOR; (c) SART+TVM; (d) SART+TVM+PRIOR
在有限角度为[30°,150°]和[30°,120°]、投影角度数分别为166和126时,各算法重建的图像如
图 15. 有限角度为[30°,150°]和[30°,120°]时重建的图像,显示窗口为[0,1]。(a) SART; (b) SART+PRIOR; (c) SART+TVM; (d) SART+TVM+PRIOR
Fig. 15. CT images reconstructed in the range of [1°, 150°] and [1°, 120°], display window is [0, 1]. (a) SART; (b) SART+PRIOR; (c) SART+TVM; (d) SART+TVM+PRIOR
4 结论
将工件的CAD设计模型作为先验知识,从中提取需要的先验图像,利用优化迭代重建算法SART+TVM+PRIOR来实现该类对象的有限角CT重建。在多材质工件的仿真实验中,选择好适当的重建参数后,加入先验图像重建的图像质量要远远优于未加入的重建质量。在单材质工件的实际扫描实验中,在所提重建算法中加入工件的先验图像能巨大地改善图像质量,工件的边缘结构更加清晰,减弱或消除了滑坡伪影,图像更加平滑明亮,这有利于判断和确定缺陷的具体位置。然而,这种算法不能改善大量投影数据缺失造成的缺陷伪影。对于多材质工件的像素截面灰度的赋值,需根据CT系统X射线能量和模型中各部分材质来确定,对于先验图像的获取,需要更快更好的配准方法,这都需要进行下一步的研究。
[1] 张瀚铭. CT不完全数据重建算法研究[D]. 北京: 信息工程大学, 2017.
Zhang HM. Computed tomography image reconstruction with incomplete projection data[D]. Beijing: University of Information Engineering, 2017.
[2] 王成祥. 有限角CT的正则化图像重建算法研究[D]. 重庆: 重庆大学, 2016.
Wang CX. Regularization based image reconstruction algorithms for limited-angle CT[D]. Chongqing: Chongqing University, 2016.
[4] Jiang M, Wang G. Development of iterative algorithms for image reconstruction[J]. Journal of X-Ray Science and Technology, 2002, 10: 77-86.
[6] Dempster A P, Laird N M, Rubin D B. Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society: Series B (Methodological), 1977, 39(1): 1-22.
[7] Browne J, de Pierro A B. A row-action alternative to the EM algorithm for maximizing likelihood in emission tomography[J]. IEEE Transactions on Medical Imaging, 1996, 15(5): 687-699.
[8] 马继明, 张建奇, 宋顾周, 等. 全变分约束迭代滤波反投影CT重建[J]. 光学学报, 2015, 35(2): 0234002.
[11] Bredies K, Kunisch K, Pock T. Total generalized variation[J]. SIAM Journal on Imaging Sciences, 2010, 3(3): 492-526.
[12] Zhang Y, Zhang W, Lei Y, et al. Few-view image reconstruction with fractional-order total variation[J]. Journal of the Optical Society of America. A, Optics, Image Science, and Vision, 2014, 31(5): 981-995.
[14] 卢孝强, 孙怡. 基于乘性正则化的有限角度CT重建算法[J]. 光学学报, 2010, 30(5): 1285-1290.
[15] Yu W, Wang C X, Huang M. Edge-preserving reconstruction from sparse projections of limited-angle computed tomography using ℓ0-regularized gradient prior[J]. The Review of Scientific Instruments, 2017, 88(4): 043703.
[16] Chen G H, Tang J, Leng S. Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[J]. Medical Physics, 2008, 35(2): 660-663.
[17] 刘进, 亢艳芹, 顾云波, 等. 稀疏张量约束的低剂量CT图像重建[J]. 光学学报, 2019, 39(8): 0811004.
[19] 张伶俐. 有限角及低信噪比CT图像重建的优化模型与算法研究[D]. 重庆: 重庆大学, 2018.
Zhang LL. Optimization model and algorithm analysis for limited-angle and low SNR CT image reconstruction[D]. Chongqing: Chongqing University, 2018.
[21] 张海娇, 孔慧华, 孙永刚. 基于结构先验的加权NLTV能谱CT重建算法[J]. 光学学报, 2018, 38(8): 0811003.
[23] 朱斯琪, 王珏, 蔡玉芳. 基于改进型生成对抗网络的低剂量CT去噪算法[J]. 光学学报, 2020, 40(22): 2210002.
Article Outline
余浩松, 邹永宁, 张智斌, 姚功杰, 周日峰. 利用CAD模型的不完全扫描CT图像重建[J]. 光学学报, 2021, 41(6): 0611002. Haosong Yu, Yongning Zou, Zhibin Zhang, Gongjie Yao, Rifeng Zhou. Image Reconstruction of Incomplete CT Scanning Using a CAD Model[J]. Acta Optica Sinica, 2021, 41(6): 0611002.
PDF全文
