1.Beijing Advanced Innovation Center for Imaging Theory and Technology, School of Mathematical Sciences, Capital Normal University, Beijing 100048, China 2.Pazhou Lab, Guangzhou 510335, China 3.Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
Fund Project:Project supported by the National Natural Science Foundation of China (Grant Nos. 61671311, 61827809), the National Key Research and Development Program of China (Grant No. 2020YFA0712200), and the National Defense Technology Foundation Project (Grant No. JSZL2018208C003)
Received Date:22 December 2020
Accepted Date:31 January 2021
Available Online:28 May 2021
Published Online:05 June 2021
Abstract:X-ray differential phase contrast computed tomography imaging based on grating interferometer system can reconstruct not only the linear attenuation coefficient, but also the phase shift coefficient and the linear scattering coefficient of the object. In practical application, it is very difficult to make a large area grating, so the sample is often larger than the grating. When the sample is scanned with a grating smaller than the sample, the part of the sample beyond the field of view of the grating will cause the differential phase projection information to be truncated. In this paper, a method of reconstructing the region of interest for differential phase contrast computed tomography is proposed. The method is based on the approximate linear relation between the phase shift coefficient of the object and the linear attenuation coefficient (i.e. the decrement in the real part of the refractive index and the imaginary part of the refractive index), the phase shift coefficient of the region of interest is approximately reconstructed by the polynomial of Lambda function of the phase shift coefficient and Lambda inverse function of linear attenuation coefficient. In this paper, according to the Fresnel diffraction theory and differential phase grating phase step-by-step method of imaging a simulation experiment is performed. In the experiment, conducted is the approximate reconstruction by using the first order polynomial and quadratic polynomial of Lambda function of the phase shift coefficient and Lambda inverse function of linear attenuation coefficient. The sample size is five times of grating imaging field, and the results show that this method can approximately reconstruct the region of interest for the sample image. We also carry out the actual data experiment. The actual data are obtained by the Talbot grating interferometer system of Shanghai synchrotron radiation BL13W1 station, and the standard model and biological sample are imaged. The method of reconstructing the region of interest is proposed in this paper. This method can be applied to the multi-material samples with a similar relationship between the decrement in the real part of the refractive index and the decrement in the imaginary part of the refractive index, and also to single-material samples. The comparison between the numerical simulations and the actual experimental results verifies the effectiveness of the proposed method. Keywords:differential phase contrast/ computed tomography/ reconstruction of region of interest/ phase shift coefficient
本文依据菲涅耳衍射理论, 模拟了X射线光栅干涉仪相位步进法的成像过程. 首先采用单色平行束光源照射样品, 模拟了穿过物质后X射线的传播过程; 然后通过移动吸收光栅获得每一步所探测到的光强信息; 最后利用相位步进法提取得到样品的微分相位投影和吸收投影[23-26]. 模拟实验中, 所用的X射线能量为30 keV, 探测器为$205$个探元的线阵探测器, 每个探元为0.08 mm. 成像示意图如图1所示, 样品的尺寸远大于光栅尺寸. 两个光栅的尺寸均为16.4 mm, 相位光栅的周期为0.016 mm, 吸收光栅的周期为0.008 mm, 两个光栅之间的距离为$7.743 \times {10^2}\;{\rm{mm}}$. 实验样品的直径为82 mm, 由大小不同的圆和椭圆组成, 模拟时离散的图像尺寸为$2048 \times 2048$, 每个像素大小为0.04 mm; 黑色部分的材质为PMMA, 折射率实部减小量$\delta $为$2.968 \times {10^{ - 7}}$, 折射率虚部$\beta $为$1.02 \times {10^{ - 10}}$, 其余白色部分的$\delta $和$\beta $均为$0$; 红色虚线内为光栅成像能够覆盖的感兴趣区域(region of interest, ROI), 直径为16.4 mm. 图2(a)为信息提取所得到的吸收投影的正弦图, 图2(b)为微分相位投影的正弦图, 扫描角度为360°, 扫描角度数为720个. 图 1 模拟实验成像示意图 Figure1. Imaging schematic diagram of simulation experiment.
图 2 (a) 感兴趣区域吸收投影的正弦图; (b) 感兴趣区域微分相位投影的正弦图 Figure2. (a) Sinogram of absorption projection for the ROI; (b) sinogram of differential phase projection for the ROI.
对感兴趣区域正弦图采用本文方法进行重建, 重建结果见图3. 图3(a)为采用一次多项式近似重建的结果图, 图3(b)为采用二次多项式近似重建的结果图. 多项式各项的系数采用样品自身的相移系数值, 通过(16)式和(17)式获得, 一次多项式为${P_1} = {{1}}{{.29}}\varLambda \delta +{{ 0}}{{.28}}{\varLambda ^{ - 1}}\beta $, 二次多项式为${P_2} = $$- \; {{0}}{{.56}}\varLambda \delta - {{1}}{{.26}}{\varLambda ^{ - 1}}\beta +{{ 0}}{{.28\;(}}\varLambda \delta {)^2} +{{ 1}}{{.61}}{\varLambda ^{ - 1}}\beta \varLambda \delta\; +$$ {{ 2}}{{.05(}}{\varLambda ^{ - 1}}\beta {)^2} $. 图3(c)为图3(a)和图3(b)在绿色虚线位置处的剖线图, 蓝色曲线表示采用一次多项式近似重建的结果, 红色曲线表示采用二次多项式近似重建的结果, 黑色曲线表示理论值. 从图3(c)可以看出采用二次多项式近似重建的结果更接近理论值. 从表1的均方误差(mean-square error, MSE)和峰值信噪比(peak signal-to-noise ratio, PSNR)的值也可以看出二次多项式近似重建结果优于一次多项式近似重建结果. 图 3 相移系数重建结果 (a) 采用一次多项式近似的重建图像; (b)采用二次多项式近似的重建图像; (c) 图3(a)和图3(b)在绿色虚线位置处的剖线图 Figure3. Reconstruction results of phase shift coefficient: (a) The reconstruction image using a first order polynomial approximation; (b) the reconstruction image using a second order polynomial approximation; (c) Fig.3 (a) and Fig.3 (b) in the green dotted line location profile chart.
方法
MSE
PSNR
一次多项式
0.0796
10.9894
二次多项式
0.0271
15.6647
表1一次多项式和二次多项式重建结果的MSE和PSNR Table1.MSE and PSNR of reconstruction results of the first order polynomial and the second order polynomial.
23.2.实际实验 -->
3.2.实际实验
本文实验采用上海同步辐射BL13 W1站的Talbot光栅干涉仪系统, 数据采集中所用射线的光子能量为20 keV, 探测器为2048 × 2048个探元的面阵探测器, 探元大小为6.5 μm. 实验所用相位光栅的周期为2.396 μm, 成像面积直径为20 mm左右, 吸收光栅的周期为2.4 μm, 成像面积直径为20 mm左右, 两块光栅之间的距离为46.38 mm. 为了更好的验证本文方法重建的效果, 实验采用对完整的投影数据截断的方式, 获得感兴趣区域投影数据. 第一组实验采用的模体由LDPE, PMMA, PTFE三种成分不同的圆柱以及水组成, 将这三个圆柱和水放在一个外直径为10.2 mm的聚乙烯塑料试管内. 图4(a)是模体的照片. LDPE, PMMA, PTFE三种成分和水的折射率实部减小量$\delta $如表2所示. 实验扫描角度为180°, 角度采样数540个. 采用8步相位步进法, 通过信息提取可获得样品的吸收投影数据和微分相位投影数据. 取一个断层的投影数据, 两点合并后的探测器单元个数为1024个, 截取的覆盖感兴趣区域的探测器单元个数为512个. 图 4 相移系数重建结果 (a)实验模体; (b)全局数据重建图像, 红色虚线内为感兴趣区域图像; (c)截断数据采用一次多项式近似的重建图像; (d)截断数据采用二次多项式近似的重建图像; (e) 图4(c)和图4(d)在绿色虚线位置处的剖线图 Figure4. Reconstruction results of phase shift coefficient: (a) Experimental modle; (b) the reconstruction image of global data, the ROI image is in the red dotted line; (c) the reconstruction image of the truncated data using a first order polynomial approximation; (d) the reconstruction image of the truncated data using a second order polynomial approximation; (e) Fig.4 (c) and Fig.4 (d) in the green dotted line location profile chart.
材料
水(H2O)
PTFE (C2F4)
PMMA (C5O2H8)
LDPE (C2H4)
$\delta $/10–7
5.26
9.65
6.30
5.46
表2水、PTFE、PMMA、LDPE的折射率实部减小量$\delta $ Table2.The decrement of the real part of the refractive index of water, PTFE, PMMA, and LDPE.
截断数据的感兴趣区域重建和全局数据重建的结果见图4. 全局数据重建采用的是将FBP重建算法中的滤波器改为Hilbert滤波器的方式[20], 重建结果见图4(b), 图中圆环为聚乙烯塑料试管, 内部不同灰度的小圆为不同成分的圆柱截面, 其余灰色部分为水, 红色虚线内为截断的感兴趣区域图像. 图4(c)和图4(d)是采用本文方法对截断数据的重建结果, 图4(c)为采用一次多项式近似重建的结果图, 图4(d)为采用二次多项式近似重建的结果图. 多项式系数利用感兴趣区域内已知材质的相移系数值, 通过(16)式和(17)式获得, 一次多项式为${P_1} = 0.08\varLambda \delta +0.79{\varLambda ^{ - 1}}\beta $, 二次多项式为${P_2} = $$0.76\varLambda \delta \;+\; 0.35{\varLambda ^{ - 1}}\beta\; -\; 0.15\;{{\rm{(}}\varLambda \delta )^2} - 0.56{\varLambda ^{ - 1}}\beta \varLambda \delta\; +$$0.61{{\rm{(}}{\varLambda ^{ - 1}}\beta )^2} $. 图4(e)为图4(c)和图4(d)在绿色虚线位置处的剖线图, 其中蓝色曲线为一次多项式近似重建的结果, 绿色曲线为二次多项式近似重建的结果, 黑色曲线为理论值. 对比图4(b)—图4(d)的重建结果, 可以看出本文方法可以很好地重建出感兴趣区域图像, 与全局数据重建的感兴趣区域图像相比, 对比度略低. 由于LDPE材质的$\delta $为$5.46 \times {10^{ - 7}}$, 水的$\delta $为$5.26 \times $$ {10^{ - 7}}$二者之间相差不大, 图4(b)中LDPE和水无法区分开, 在图4(c)和图4(d)中LDPE和水稍有区别. 从图4(e)可以看出采用二次多项式近似重建结果比一次多项式近似重建结果更接近理论值. 第二组实验样品是一只仓鼠的前趾, 长度为5 mm左右, 实验参数与第一组实验相同. 选取光栅成像视野的1/4作为感兴趣区域, 探测器单元个数为256. 图5为全局吸收投影和微分相位投影的正弦图, 其中两条红色虚线间为截取的感兴趣区域正弦图, 图5(a)为吸收投影的正弦图, 图5(b)为微分相位投影的正弦图. 图 5 (a) 全局吸收投影的正弦图; (b)全局微分相位投影的正弦图. 两红色虚线间为截断的感兴趣区域正弦图 Figure5. (a) Sinogram of global absorption projection; (b) sinogram of the global differential phase projection. Between the two red dotted line for ROI of truncation sinogram.
截断数据的感兴趣区域重建和全局数据重建的结果见图6. 图6(a)为全局数据重建的结果图, 图6(b)为全局数据的感兴趣区域重建图像; 图6(c)为截断数据的感兴趣区域重建结果图, 重建过程中采用的是一次多项式${P_1} = \varLambda \delta + 0.58{\varLambda ^{ - 1}}\beta $, 由于生物软组织的主要元素为C, H, O, 其折射率实部减小量与水的比较接近, 因此多项式系数利用实验一中水的相移系数值, 通过(16)式和(17)式获得. 由于本文方法利用了微分相位投影和吸收投影两种信息, 对比图6(b)和图6(c)可以看出本文方法重建图像比全局的微分相位重建图像细节更加丰富. 图 6 相移系数重建图像 (a) 全局数据重建图像, 红色虚线内为感兴趣区域图像; (b) 全局数据的感兴趣区域重建图像; (c) 截断数据采用一次多项式近似的重建图像 Figure6. Reconstruction image of phase shift coefficient: (a) The reconstruction image of global data, the ROI image is in the red dotted line; (b) the reconstruction image of the ROI from the global data; (c) the reconstruction image of the truncated data using a first order polynomial approximation.