一种用于X射线光栅相位衬度成像的背景扣除方法

摘要

本发明公开了一种用于X射线光栅相位衬度成像的背景扣除方法，成像系统包括X光机(1)、源光栅(2)、分束光栅(3)、样品室(4)、分析光栅(5)和X射线探测器(6)。本发明利用相位步进法采集到的图像，通过循环移动样品图像和背景图像的顺序，而使各像素的初始相位远离相位跳变区，从而可以使样品产生的相位没有跨过跳变区，进而可以避免相位纠缠现象的发生。同时基于本发明专利的背景扣除方法可以灵活调节相位测量的范围，如调节为(-π+3，π+3]，而传统方法的相位测量范围是固定在(-π，π]区间内的，本发明专利因此可以在某些应用情况下获得优势。

说明书

技术领域

本发明涉及X射线相位衬度成像技术领域，具体涉及一种用于X射 线光栅相位衬度成像的背景扣除方法。

背景技术

对于那些由轻元素组成的物质，硬X射线相位衬度成像技术能够提供 比传统吸收成像高达上千倍的图像衬度和测量灵敏度，该技术因而在生物 软组织成像方面有巨大的应用前景(参考文件[1，2])。在目前众多的X射 线相位衬度成像技术中，基于三块光栅的X射线相位衬度成像技术是最有 可能获得推广和实际运用的(参考文件[3-5])，这是因为该方法可以利用 普通的X光机进行X射线相位衬度成像，从而使X射线相位衬度成像技 术摆脱了同步辐射光源和微焦点光源的约束。

在基于光栅干涉仪的X射线相位衬度成像技术中，为了得到纯样品的 折射信息，背景扣除是一个必不可少的步骤(参考文件[6，7])。在实际试 验中，往往是首先把样品移入成像视场，利用相位步进方法采集一套图像， 然后把样品移出成像视场，采集另一套数据，最后利用这两套数据进行背 景扣除运算。目前普遍采用的背景扣除方法主要有两种，第一种方法称为 Arg(S)-Arg(B)，该方法的相位测量范围是(-2π，2π]，缺点是在样品产生的 相位没有超出(-π，π]范围的时候，也很有可能会产生相位纠缠现象(参考 文件[6，7])；第二种方法称为Arg(S/B)法(参考文件[8])，该方法的相位测 量范围是固定在区间(-π，π]内的，缺点是在样品产生的相位超出(-π，π] 范围的时候，该方法会产生相位缠绕现象(参考文件[7])，而解除相位缠 绕现象需要非常复杂的算法(参考文件[9，10])。

参考文献：

[1]A.Momose，″Recent advances in X-ray phase imaging，″Japanese  Journal ofApplied Physics，vol.44，p.6355，2005.

[2]A.Bravin，P.Coan，and P.Suortti，″X-ray phase-contrast imaging： from pre-clinical applications towards clinics，″Physics in medicine and  biology，vol.58，p.R1，2013.

[3]F.Pfeiffer，T.Weitkamp，O.Bunk，and C.David，″Phase retrieval  and differential phase-contrast imaging with low-brilliance X-ray sources，″ Nature Physics，vol.2，pp.258-261，Apr 2006.

[4]F.Pfeiffer，C.Kottler，O.Bunk，and C.David，″Hard x-ray phase  tomography with low-brilliance sources，″Physical review letters，vol.98，p. 108105，2007.

[5]F.Pfeiffer，M.Bech，O.Bunk，P.Kraft，E.F.Eikenberry，C. et al.，″Hard-X-ray dark-field imaging using a grating  interferometer，″Nature materials，vol.7，pp.134-137，2008.

[6]A.Momose，W.Yashiro，Y.Takeda，Y.Suzuki，and T.Hattori， ″Phase tomography by X-ray Talbot interferometry for biological imaging，″ Japanese journal of applied physics，vol.45，p.5254，2006.

[7]I.Zanette，″Interférométrie X à féseaux pour l′imagerie et l′analyse  de front d′ondes au synchrotron，″Université de Grenoble，2011，Page 44.

[8]S.Wang，R.Hu，Z.Wang，K.Gao，K.Zhang，A.Momose，et al.， ″Experimental research on the feature of Talbot-Lau interferometer vs.tube  accelerating voltage，″arXiv preprint arXiv：1410.7655，2014.

[9]T.R.Judge and P.Bryanston-Cross，″A review of phase  unwrapping techniques in fringe analysis，″Optics and Lasers in Engineering， vol.21，pp.199-239，1994.

[10]K.Itoh，″Analysis of the phase unwrapping algorithm，″Applied  Optics，vol.21，pp.2470-2470，1982.

发明内容

(一)要解决的技术问题

本发明指在提出一种新的背景扣除方法，用于X射线光栅相位衬度成 像中，一方面消除样品引起相位没有超出(-π，π]时候产生的相位纠缠现象； 另一方面调节测量相位的范围，如调节测量范围到(-π+3，π+3]，从而可以 适应于不同的应用环境。

(二)技术方案

本发明提出一种用于X射线光栅相位衬度成像的背景扣除方法，应用 于X射线相位衬度成像系统，该系统包括X光机(1)、源光栅(2)、分 束光栅(3)、样品室(4)、分析光栅(5)和X射线探测器(6)。

所述方法包括如下步骤：

S1、在垂直于光路的横向平面上，沿垂直于分析光栅栅条的方向上使 分析光栅在一个光栅周期内逐步移动，在原始位置和每步移动后的位置均 分别采集样品图像和背景图像，得到的样品图像记为S_j，得到的背景图像 记为B_j，j表示图像序号并有j＝0，1，……，M-1；

S2、利用如下公式计算所述样品图像和背景图像的折射角矩阵Φ：

其中(x，y)是所述样品图像或背景图像的像素坐标，是折射角矩 阵Φ的元素，表示像素坐标(x，y)处的折射角，z_T是样品和分析光栅之间的 距离，k是整数，I_k(x，y)代表所述像素坐标(x，y)处图像的灰度值，简写为 I_k，arg[]表示对[]中的复数求复角运算，其中：

当I₀，I₁，…，I_M-1分别为样品图像S₀，S₁，…，S_M-1的像素的灰度值时，得到样 品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别为背 景图像B₀，B₁，…，B_M-1的像素的灰度值时，得到背景图像的折射角矩阵其内部元素为

当I₀，I₁，…，I_M-1分别为样品图像S_M-1，S₀，S₁，…，S_M-2的像素的灰度值时，得 到样品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别 为背景图像B_M-1，B₀，B₁，…，B_M-2的像素的灰度值时，得到背景图像的折射角 矩阵其内部元素为

……；

当I₀，I₁，…，I_M-1分别为样品图像S₁，S₂，…，S_M-1，S₀的像素的灰度值时，得到 样品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别为 背景图像B₁，B₂，…，B_M-1，B₀的像素的灰度值时，得到背景图像的折射角矩阵 其内部元素为

S3、对于背影图像的折射角矩阵中像素坐标(1，1)的各元 素比较其绝对值的大小，取其中绝对值最小的 元素的从而确定其对应的样品图像的折射角矩阵中的元素i为自然数且不大于M，通过公式计算得到扣除背景后的 样品图像的像素坐标(1，1)处的折射角

S4、对背景图像和样品图像中的每一个像素(x，y)，重复步骤S3，从 而得到样品的折射角矩阵Φ_i的每一个元素

根据本发明的优选实施方式，所述步骤S1中，每次移动分析光栅的 距离均为d/M，d是所述分析光栅的周期，M是上述过程中总共采集图像 的数量

根据本发明的优选实施方式，M大于或等于3。

(三)有益效果

(1)本发明通过循环移位原始图像的顺序而使各像素的初始相位远 离相位跳变区，从而可以使样品产生的相位没有跨过跳变区，进而可以避 免在样品产生的相位没有超出(-π，π]范围的时候产生的相位纠缠现象。

(2)本发明可以在一定程度上调节测量相位的区间，(传统Arg(S/B) 方法测量的相位区间是(-π，π]，是固定不变的)。假设样品的折射引起的相 位是π+0.1，那么根据Arg(S/B)方法计算得到的值是-π+0.1，此时就会引起 计算错误，而利用本发明提出的方法，通过合理移位，即可保证测量的相 位区间是(-π+2，π+2]，从而可以准确地测量出样品的相位π+0.1。

附图说明

图1是X射线光栅相位衬度成像系统光学结构图；

图2是本发明的背景扣除新方法的原理示意图，其中(a)图是位移 曲线示意图，(b)图是测量得到的相位的示意图；

图3是本发明的一个实施例的X射线光栅相位衬度成像系统的结构图；

图4是本发明的一个实施例的原始图像示意图，s①，s②，s③，s④ 和s⑤是五张样品图像，b①，b②，b③，b④和b⑤是五张背景图像；

图5是本发明的一个实施例的三种背景扣除方法得到的折射图像的示 意图，其中(a)图是利用传统Arg(S)-Arg(B)方法计算的结果，(b)图是 利用传统Arg(S/B)方法计算得到的结果，(c)是利用本发明计算得到的结 果；

图6是利用传统Arg(S/B)方法以及本发明提出的方法得到的实验结果 的定量比较图。

具体实施方式

本发明的总体构思是：充分利用相位步进法采集到的图像的特点，通 过循环移位原始样品图像和背景图像的顺序，而使各像素的初始相位远离 相位跳变区，从而可以使样品产生的相位没有跨过跳变区，进而可以避免 相位纠缠现象的发生。

具体来说，本发明的用于X射线光栅相位衬度成像的背景扣除方法应 用于X射线相位衬度成像系统，X射线相位衬度成像系统至少包括有X 光机、源光栅、分束光栅、样品室、分析光栅和X射线探测器。本发明的 方法包括如下步骤：

S1、在垂直于光路的横向平面上，沿垂直于分析光栅栅条的方向上使 分析光栅在一个光栅周期内逐步移动，在原始位置和每步移动后的位置均 分别采集样品图像和背景图像，得到的样品图像记为S_j，得到的背景图像 记为B_j，j表示图像序号并有j＝0，1，……，M-1。一般，每次移动分析光栅的 距离均为d/M，d是所述分析光栅的周期，M是上述过程中总共采集图像 的数量。M通常大于或等于3。

S2、利用如下公式计算所述样品图像和背景图像的折射角矩阵Φ：

其中(x，y)是所述样品图像或背景图像的像素坐标，是折射角矩 阵Φ的元素，表示像素坐标(x，y)处的折射角，z_T是样品和分析光栅之间的 距离，k是整数，I_k(x，y)代表所述像素坐标(x，y)处图像的灰度值，简写为 I_k，arg[]表示对[]中的复数求复角运算，其中：

当I₀，I₁，…，I_M-1分别为样品图像S₀，S₁，…，S_M-1的像素的灰度值时，得到样 品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别为背 景图像B₀，B₁，…，B_M-1的像素的灰度值时，得到背景图像的折射角矩阵其内部元素为

当I₀，I₁，…，I_M-1分别为样品图像S_M-1，S₀，S₁，…，S_M-2的像素的灰度值时，得 到样品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别 为背景图像B_M-1，B₀，B₁，…，B_M-2的像素的灰度值时，得到背景图像的折射角 矩阵其内部元素为

……；

当I₀，I₁，…，I_M-1分别为样品图像S₁，S₂，…，S_M-1，S₀的像素的灰度值时，得到 样品图像的折射角矩阵其内部元素为当I₀，I₁，…，I_M-1分别为 背景图像B₁，B₂，…，B_M-1，B₀的像素的灰度值时，得到背景图像的折射角矩阵 其内部元素为

S3、对于背影图像的折射角矩阵中像素坐标(1，1)的各元 素比较其绝对值的大小，取其中绝对值最小的 元素的从而确定其对应的样品图像的折射角矩阵中的元素i为自然数且不大于M，通过公式计算得到扣除背景后的 样品图像的像素坐标(1，1)处的折射角

S4、对图像中的每一个像素(x，y)，重复步骤S3，重复步骤S3，从而 得到样品的折射角矩阵Φ_i的每一个元素

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体五 次图像取样的实施例，对本发明作进一步的详细说明。也就是说，在该实 施例中，M＝5。

步骤S1、图像采集过程。图1是该实施例的X射线光栅相位衬度成 像系统的光学结构图。如图1所示，主要包括X光机1，源光栅2，分束 光栅3，样品室4，分析光栅5和X射线探测器6。在图像采样时，首先 对X射线探测器进行暗场和明场校正，然后以莫尔条纹为参考对象，尽量 对准三块光栅，使视场中莫尔条纹的周期尽量无穷大。然后，采集第一张 图像；接着在垂直于光栅栅条方向移动分析光栅1μm，光栅运动完成后， 采集第二张图像；继续移动分析光栅1μm，光栅运动完成后，采集第三张 图像；然后继续移动分析光栅1μm，光栅运动完成后，采集第四张图像； 最后移动分析光栅1μm，光栅运动完成后，采集第五张图像。然后把样品 移入视场，把分析光栅回归原始位置，按照上述同样的步骤，采集五张样 品图像。

步骤S2、具体计算过程。把5张背景图像依次编号为①、②、③、④、 ⑤，5张样品图像也编号为①、②、③、④和⑤，背景图像和样品图像均 排成序列①②③④⑤，利用如下公式分别计算背景图像和样品图像的折射 角矩阵Φ：

公式中的参数含义与前相同，且M＝5，d＝5μm。

该排列序号下得到的样品图像的折射角矩阵背景图像的折射角 矩阵

接着，将背景图像和样品图像改变排列序列分别为⑤①②③④、 ④⑤①②③、③④⑤①②、②③④⑤①，再次依上述公式进行计算，得到 样品图像的折射角矩阵背景图像的折射角矩阵${\Phi }_3^B,{\Phi }_4^B,{\Phi }_5^B$

步骤S3、对对于背影图像的折射角矩阵的第一个元素，即x＝1、y＝1， 根据计算得到的5个背景图像在x＝1，y＝1像素处的5个折射角，即和寻找这五个折射角中绝对值最小的 那个数所处的位置序号，记为i，然后在对应的5个样品图像的折射角 [和]中寻找第i个数，利用 计算出第一个像素处的折射角

S4、对背景图像和样品图像中的每一个像素(x，y)，重复步骤S3，从 而得到样品的折射角矩阵Φ_i的每一个元素

上述计算方法的原理如图2所示。图2的(a)图中的黑色曲线是某 一个像素的光强随分析光栅横向运动而变化的位移曲线，一般来说，该曲 线是一个标准的余弦曲线。图像中的①、②、③、④和⑤代表相位步进方 法采集的五张图像的位置，这五个位置均匀的把位移曲线的一个周期分成 5等份。图2的(b)图中的黑色曲线代表提取出来的相位，该曲线与(a) 图中的曲线的一一对应，对应关系描述如下，假设相位步进采集的5张图 像是①、②、③、④和⑤，那么提取出来的相位就是位置①垂直投影在图 2(b)中曲线的点的纵坐标的值。在样品引起相位没有超出(-π，π]时候， Arg(S)-Arg(B)背景扣除方法会产生相位纠缠现象的原因是：没有样品的时 候对应的位置在图2的(b)图曲线的跳变区附近，加上样品以后， 位移曲线产生移动，因此计算出来的相位可能会越过图2的(b) 图曲线的跳变区，此时在利用公式进行背景扣除的时候， 就会产生一个2π或者-2π的增量，这就是错误产生的原因。本发明的核心 思想就是使每一个像素的初始相位都尽可能的远离跳变区，远离的办法是 通过循环移位原始相位步进图像来实现的，如把①②③④⑤移位成 ⑤①②③④，那么根据余弦曲线的性质，我们知道其等效于⑤¹①②③④， 如图2的(a)图所示，这里认为点⑤¹，④¹，③¹，②¹在位移曲线中的相 位和点⑤，④，③，②在位移曲线中的相位相同。那么根据序列⑤¹①②③④ 计算得到的相位就是如图2的(b)图中的⑤①②③④所示。同样的办法， ④¹⑤¹①②③等效于④¹⑤¹①②③；③¹④¹⑤¹①②等效于③¹④¹⑤¹①②； ②¹③¹④¹⑤¹①等效于②¹③¹④¹⑤¹①，我们可以分别求出序列④⑤①②③、 ③④⑤①②、②③④⑤①的相位。计算发现，这五种图像序列计算得到的 相位均匀地分布在图2中曲线的一个周期上。在这5个相位中，找到最远 离跳变区，也就是绝对值最小的那个值。那么在加上样品的时候，由于样 品折射引起的相位变化就有最大的可能不会跨过跳变区，因此也不会引起 计算错误。

以下公式证明了样品图像和背景图像经过同样的循环移位后，纯样品 的折射信息固定不变，从而证实了本发明提出的背景扣除方法的可行性。

需要指出的是在上述公式中，sample表示样品图像，background表示 背景图像，我们假定背景的初始相位稍微偏离跳变边界，同时样品引起的 相位移动非常小，不会跨过跳变区。

需要指出的是，在上述本发明的方法介绍中，我们均以5步作为实例 说明其操作过程，在实际应用中，任何一个大于或等于3的整数都是可以 的。另需要说明的是，本发明提出的背景扣除方法，同样可以运用在基于 两块光栅的相位衬度成像技术和基于晶体干涉仪的X射线相位衬度成像 技术中。

图3是该发明的具体实施方式的成像系统的结构图，X光机603、源 光栅604、分束光栅612、样品室613和分析光栅618均通过光学精密位 移台602固定在光学平台601上，在分析光栅618横向运动方向，配置有 超精密的压电陶瓷电机617，用于高精密的相位步进扫描，X射线探测器 619通过机械部件固定。X光机603的焦点大小为1mm，工作管电压为40 KV，管电流为45mA。源光栅604周期为22.7μm，金厚度为70μm，占空 比为1∶1，面积为2×2cm²。分束光栅612周期为4.36μm，金厚度为2.43μm， 占空比为1∶1，面积为5×5cm²。分析光栅618周期为5.4μm，金厚度为65μm， 占空比为1∶1，面积为5×5cm²。X射线探测器的像素尺寸为18×18μm²，探 测器619的活动区域面积为6.84×6.84cm²。源光栅604距离X光机603出 光点的距离为8cm，源光栅604和分束光栅612的距离为106.9cm，样品 室613位于分束光栅612后面，分析光栅618和分束光栅612的距离为 25.6cm，X射线探测器619紧贴着分析光栅618放置。

三块光栅经过精确对准后，固定好样品，该实施例中样品是一根直径 为10mm的PMMA玻璃棒、一根直径为5mm的POM玻璃棒、一根直径 为5mm的PMMA玻璃棒和一根直径为10mm的POM玻璃棒。首先采集 第一张样品图像①；然后在垂直于光栅栅条方向移动分析光栅1μm，光栅 运动完成后，采集第二张样品图像②；接着移动分析光栅1μm，光栅运动 完成后，采集第三张样品图像③；然后继续移动分析光栅1μm，光栅运动 完成后，采集第四张样品图像④；最后移动分析光栅1μm，光栅运动完成 后，采集第五张样品图像⑤。然后把样品移出视场，按照同样的步骤，采 集背景图像①、②、③、④和⑤。采集的5张样品图像额5张背景图像分 别如图4所示。

图5是基于上述五张背景图像(b①，b②，b③，b④和b⑤)和五张 样品图像(s①，s②，s③，s④和s⑤)，分别采用2种传统背景扣除方法 以及本发明提出的背景扣除方法计算得到的样品的折射图像。(a)图是利 用传统Arg(S)-Arg(B)方法计算得到的结果，(b)图是利用传统Arg(S/B) 方法计算得到的结果，(c)图是利用本发明计算得到的结果。可以看出利 用传统Arg(S)-Arg(B)方法计算得到的实验结果中存在明显的坏点(相位缠 绕现象)，而利用传统Arg(S/B)方法以及本发明提出的方法则完全消除了 相位缠绕现象。

图6是利用传统Arg(S/B)方法以及本发明提出的方法得到的实验结果 的定量比较图。(a)图是利用本发明计算得到的折射图像，(b)图是利用 传统Arg(S/B)方法计算得到的样品的折射图像，(c)图是两种方法得到的 实验结果的断面图，(c)图中的黑色实线是本发明的断面图，断面选取如 (a)中的黑色横线所示，(c)图中的虚线是传统Arg(S/B)方法的断面图， 断面选取如图(b)图中的白色横线所示。从图中可以看出，两条曲线完 全重合，本发明的方法可以得到与传统Arg(S/B)方法一模一样的数值。

但是传统Arg(S/B)方法测量的相位区间是(-π，π]，是固定不变的。假 设样品的折射引起的相位是π+0.1，那么根据传统Arg(S/B)方法计算得到 的相位值是-π+0.1，此时就会引起错误。而利用本发明提出的新方法，通 过合理移位，即可保证测量的相位区间是(-π+2，π+2]，从而可以准确测量 相位π+0.1。需要指出的是基于本发明的这个优势在这个实施例上没有得 到体现，但本领域的技术人员可以理解。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行 了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已， 并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、 等同替换、改进等，均应包含在本发明的保护范围之内。