微分相位衬度计算机层析成像的感兴趣区域重建方法*

张敬娜 张慧滔2)† 徐文峰 朱溢佞2) 邓世沃 朱佩平3)

1) (首都师范大学数学科学学院， 北京成像理论与技术高精尖创新中心， 北京 100048)

2) (琶洲实验室， 广州 510335)

3) (中国科学院高能物理研究所， 北京 100049)

基于光栅干涉仪系统的X射线微分相位衬度计算机层析成像， 不仅可以重建物体的线性衰减系数， 还可以重建物体的相移系数和线性散射系数.在实际应用时， 大面积光栅不易获得， 常常遇到样品大于光栅的情况.当用小于样品的光栅对样品进行扫描时， 样品超出光栅成像视野的部分会导致微分相位投影信息被截断.本文针对微分相位衬度计算机层析成像提出了一种相移系数的感兴趣区域重建方法.该方法利用物体相移系数和线性衰减系数(即折射率实部减小量和折射率虚部)之间的近似线性关系; 通过重建相移系数的Lambda函数和线性衰减系数的Lambda逆函数的多项式组合， 近似重建物体感兴趣区域的相移系数.数值模拟实验依据菲涅耳衍射积分理论， 进行计算机仿真X射线的传播过程和光栅成像过程.实际实验利用上海同步辐射BL13W1站的Talbot光栅干涉仪系统， 分别对标准模体和生物样品进行光栅微分相位衬度计算机层析成像.数值模拟和实际实验结果都验证了该方法的有效性.

1 引 言

基于吸收的传统计算机层析成像(computed tomography， CT)重建的是物体内部线性衰减系数的分布， 对吸收较弱的低原子序数样品难以获得高衬度CT图像.X射线穿过样品时， 不仅有振幅衰减， 还有相位改变， 即相移.对于吸收较弱的低原子序数物体， 穿过物体后的X射线相移更明显[1，2]， 因此探测相移的X射线相位衬度成像方法与CT技术相结合， 可以获得高衬度三维CT图像.

X射线相位衬度成像方法中被认为最有希望推广到临床应用的是基于光栅的微分相位衬度成像[3-5].利用光栅微分相位衬度成像可以获得样品的吸收图像、折射图像和散射图像.当样品相对X射线旋转时， 可以获得各个角度的吸收投影图、折射投影图和散射投影图， 重建出被测样品的线性衰减系数、相移系数和线性散射系数的三维空间分布.

光栅微分相位衬度成像的视野是限制该方法广泛应用的主要原因之一.由于光栅制作工艺的限制， 目前制作的光栅面积在几十毫米左右.对于大于光栅成像视野的样品进行扫描时， 在超出视野范围部分， 微分相位投影信息被截断.此时， 虽然难以对整个样品进行微分相位衬度计算机层析成像(differential phase contrast computed tomography， DPC-CT)， 但是仍然有希望应用DPC-CT对样品的感兴趣区域进行成像.这是DPC-CT的一个重要研究方向.

传统的吸收CT感兴趣区域重建可以分为精确重建方法和近似重建方法.反投影滤波(backprojection filter， BPF)[6，7]类的感兴趣区域重建属于精确重建.该方法通过对感兴趣区域的投影数据微分， 然后反投影到感兴趣区域， 最后沿着感兴趣区域的PI线进行有限区间的逆Hilbert变换， 重建感兴趣区域图像.但是有限区间的逆Hilbert变换时需要满足沿Hilbert变换方向上至少包含一个端点， 所以BPF的方法不适用于不包含图像边界的内部感兴趣区域重建问题.Lambda感兴趣区域重建[8-10]属于近似重建.该方法通过重建样品线性衰减系数的Lambda函数和Lambda逆函数的线性组合来近似重建出感兴趣区域图像.

光栅微分相位衬度成像获得的是被测样品的吸收、微分相位和散射投影.要从微分相位投影重建感兴趣区域的相移系数， 一些吸收CT的感兴趣重建方法难以直接应用.Anastasio等[11]虽然证明了BPF方法可用于物体局部相移系数的重建， 但是BPF方法应用的限制条件依然存在.迭代类算法也被用于DPC-CT的感兴趣区域重建， Cong等[12]在折射率n呈分段多项式分布的前提下， 提出了一种可用于相移系数重建的迭代算法; Pascal等[13]利用感兴趣区域上的先验信息， 提出了一种可以消除截断伪影的POCS算法; Yang等[14]提出了一种基于样条函数和希尔伯特变换的迭代算法， 迭代方法的计算比较复杂， 用途有限.Felsner等[15]在2018年提出了用吸收投影数据补全相位缺失投影的方法， 先根据线性衰减系数的重建图像进行材料分解， 再通过感兴趣区域的相位重建图像估计出每种材料的折射率实部的减小量， 然后结合每种材料对应的吸收投影图与截断的微分相位投影图外推出全部的微分相位投影图， 最后进行重建; 2020年Felsner等[16]对他们在2018年所提出的这个重建算法进行了改进， 直接根据感兴趣区域的微分相位投影图估计每种材料的折射率实部的减小量， 这种方法避免了由于截断的相位重建图像造成的误差.

本文先给出了平行束扫描模式下Lambda函数和Lambda逆函数的重建公式， 并利用相移系数和线性衰减系数(即折射率实部减小量和折射率虚部)之间的近似关系， 提出了一种针对微分相位衬度CT的感兴趣区域重建方法.这种方法同时利用了微分相位投影和吸收投影两种信息， 通过重建相移系数的Lambda函数和线性衰减系数的Lambda逆函数的多项式组合来近似重建样品的相移系数.其中， 相移系数的Lambda函数由微分相位投影进行一阶差分重建， 线性衰减系数的Lambda逆函数由吸收投影进行反投影重建.

2 微分相位衬度CT的感兴趣区域方法

2.1 数学准备

利用光栅微分相位衬度成像可以获得样品的吸收投影、微分相位投影和散射投影.吸收投影的表达式为

其中:μ为样品的线性衰减系数，μ=2kβ;l表示射线路径;k为波数，k=2π/λ，λ为波长;β为样品的折射率虚部.

微分相位投影的表达式为

其中，pdpc表示相位投影，r表示射线的法向，δ为样品的折射率实部减小量， 一般我们重建的是样品的相移系数kδ.散射投影就是散射角的二阶矩σ2，其表达式为

其中，α为样品的线性散射系数.

本文所用的Λ算子是负拉普拉斯算子的平方根[9]， 即

其中， (x1，x2) 表示二维直角坐标系下的坐标.

δ(x1，x2) 表示点x=(x1，x2) 处的折射率实部减小量.δ(x1，x2) 的二维Fourier变换可以表示为

其 中，ω1=ωcosφ，ω2=ωsinφ，ω1，ω2是 二 维Fourier变换后的频域直角坐标，ω，φ是极坐标.

的二维逆Fourier变换可以表示为

Λδ的二维Fourier变换为

Λ-1β的二维Fourier变换为

因为光栅相位衬度成像获得的是样品的吸收投影和散射投影， 所以可以直接利用吸收CT的重建算法[17-19]对样品的线性衰减系数和线性散射系数进行重建.但光栅相位衬度成像获得的微分相位投影是相位投影的偏导数， 所以在对样品的相移系数进行重建时， 不能直接应用吸收CT的重建方法.对于样品相移系数的重建， 一种方法是先将微分相位投影信息积分， 得到相位投影， 之后采用线性衰减系数的重建方法重建样品的相移系数; 另一种方法是直接利用微分相位投影信息重建， 将吸收重建的滤波反投影(filter back-projection， FBP)算法中所用的斜坡滤波器， 改为希尔伯特滤波器即可进行滤波反投影重建[20].对微分相位投影积分和希尔伯特滤波都要求投影数据没有截断， 因此这两种方法都无法直接应用于截断数据的感兴趣区域重建.

2.2 感兴趣区域重建方法

为了简明起见， 我们只讨论平行束扫描模式下的DPC-CT的重建方法， 其结论也可推广到扇束和锥束CT.在平行束条件下，l表示射线路径，l:x·φ=r， 其中x=(x1，x2)，φ=(cosφ，sinφ) ，r为射线的法向，φ为射线旋转角度，φ∈[0，π].沿射线l的相位投影可以表示为

pdpc(r，φ)沿r方向的一维Fourier变换为

根据FBP重建公式

微分相位衬度成像Λδ的重建公式为

β(x1，x2)表示点x处的折射率虚部， 沿射线l的吸收投影可以表示为

微分相位衬度成像Λ-1β的重建公式为

根据折射率实部减小量和折射率虚部之间存在一种线性关系，δ=εβ[21，22]， 其中ε为一个与材质有关的常数.因此，

相移系数kδ则可以由Λδ和Λ-1β的多项式组合Pn近似重建.

其中，n表示多项式组合的最高次数，aji为各个项的系数.

多项式系数aji可以通过扫描已知材质模体进行求解， 假设模体包含M种材质， 每个材质的相移系数值为kδm， 所在区域为Ωm，m=1，2，···，M.多项式的系数可以通过求解优化问题

得到， 其中a=(a10，a11，···，aji).

3 实 验

3.1 模拟实验

本文依据菲涅耳衍射理论， 模拟了X射线光栅干涉仪相位步进法的成像过程.首先采用单色平行束光源照射样品， 模拟了穿过物质后X射线的传播过程; 然后通过移动吸收光栅获得每一步所探测到的光强信息; 最后利用相位步进法提取得到样品的微分相位投影和吸收投影[23-26].

模拟实验中， 所用的X射线能量为30 keV，探测器为 2 05 个探元的线阵探测器， 每个探元为0.08 mm.成像示意图如图1所示， 样品的尺寸远大于光栅尺寸.两个光栅的尺寸均为16.4 mm， 相位光栅的周期为0.016 mm， 吸收光栅的周期为0.008 mm， 两个光栅之间的距离为 7.743×102mm.实验样品的直径为82 mm， 由大小不同的圆和椭圆组成， 模拟时离散的图像尺寸为 2 048×2048 ， 每个像素大小为0.04 mm; 黑色部分的材质为PMMA， 折射率实部减小量δ为 2.968×10-7， 折射率虚部β为 1.02×10-10， 其余白色部分的δ和β均为0; 红色虚线内为光栅成像能够覆盖的感兴趣区域(region of interest， ROI)， 直径为16.4 mm.图2(a)为信息提取所得到的吸收投影的正弦图， 图2(b)为微分相位投影的正弦图， 扫描角度为360°， 扫描角度数为720个.

图1 模拟实验成像示意图Fig.1.Imaging schematic diagram of simulation experiment.

图2 (a) 感兴趣区域吸收投影的正弦图; (b) 感兴趣区域微分相位投影的正弦图Fig.2.(a) Sinogram of absorption projection for the ROI;(b) sinogram of differential phase projection for the ROI.

对感兴趣区域正弦图采用本文方法进行重建，重建结果见图3.图3(a)为采用一次多项式近似重建的结果图， 图3(b)为采用二次多项式近似重建的结果图.多项式各项的系数采用样品自身的相移系数值， 通过(16)式和(17)式获得， 一次多项式为P1=1.29Λδ+0.28Λ-1β， 二次多项式为P2=-0.56Λδ-1.26Λ-1β+0.28(Λδ)2+1.61Λ-1βΛδ+2.05(Λ-1β)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)在绿色虚线位置处的剖线图Fig.3.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.

表1 一次多项式和二次多项式重建结果的MSE和PSNRTable 1.MSE and PSNR of reconstruction results of the first order polynomial and the second order polynomial.

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三种成分和水的折射率实部减小量δ如表2所示.实验扫描角度为180°， 角度采样数540个.采用8步相位步进法， 通过信息提取可获得样品的吸收投影数据和微分相位投影数据.取一个断层的投影数据， 两点合并后的探测器单元个数为1024个， 截取的覆盖感兴趣区域的探测器单元个数为512个.

图4 相移系数重建结果 (a)实验模体; (b)全局数据重建图像， 红色虚线内为感兴趣区域图像; (c)截断数据采用一次多项式近似的重建图像; (d)截断数据采用二次多项式近似的重建图像; (e) 图4(c)和图4(d)在绿色虚线位置处的剖线图Fig.4.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.

表2 水、PTFE、PMMA、LDPE的折射率实部减小量δTable 2.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)式获得， 一次多项式为P1=0.08Λδ+0.79Λ-1β， 二次多项式为P2=0.76Λδ+0.35Λ-1β-0.15(Λδ)2-0.56Λ-1βΛδ+0.61(Λ-1β)2.图4(e)为图4(c)和图4(d)在绿色虚线位置处的剖线图， 其中蓝色曲线为一次多项式近似重建的结果， 绿色曲线为二次多项式近似重建的结果， 黑色曲线为理论值.

对比图4(b)—图4(d)的重建结果， 可以看出本文方法可以很好地重建出感兴趣区域图像， 与全局数据重建的感兴趣区域图像相比， 对比度略低.由于LDPE材质的δ为 5.46×10-7， 水的δ为5.26×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)全局微分相位投影的正弦图.两红色虚线间为截断的感兴趣区域正弦图Fig.5.(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)为截断数据的感兴趣区域重建结果图， 重建过程中采用的是一次多项式P1=Λδ+0.58Λ-1β，由于生物软组织的主要元素为C， H， O， 其折射率实部减小量与水的比较接近， 因此多项式系数利用实验一中水的相移系数值， 通过(16)式和(17)式获得.由于本文方法利用了微分相位投影和吸收投影两种信息， 对比图6(b)和图6(c)可以看出本文方法重建图像比全局的微分相位重建图像细节更加丰富.

图6 相移系数重建图像 (a) 全局数据重建图像， 红色虚线内为感兴趣区域图像; (b) 全局数据的感兴趣区域重建图像; (c) 截断数据采用一次多项式近似的重建图像Fig.6.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.

4 结 论

本文通过重建相移系数的Lambda函数和线性衰减系数的Lambda逆函数的多项式组合， 来近似重建物体的相移系数分布.由于Lambda函数和Lambda逆函数的重建具有局部性， 因此该方法可用于截断数据的感兴趣区域重建.模拟实验和生物样品的重建结果都表明了该方法可以对超过光栅成像视野大小的样品进行感兴趣区域重建.