CN104535595A - 一种用于x射线光栅相位衬度成像的背景扣除方法 - Google Patents

Abstract

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

Description

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

技术领域

本发明涉及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，″JapaneseJournal 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 andbiology，vol.58，p.R1，2013.

[3]F.Pfeiffer，T.Weitkamp，O.Bunk，and C.David，″Phase retrievaland 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 phasetomography 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 gratinginterferometer，″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′analysede 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.tubeaccelerating voltage，″arXiv preprint arXiv：1410.7655，2014.

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

[10]K.Itoh，″Analysis of the phase unwrapping algorithm，″AppliedOptics，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。

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

接着，将背景图像和样品图像改变排列序列分别为⑤①②③④、④⑤①②③、③④⑤①②、②③④⑤①，再次依上述公式进行计算，得到样品图像的折射角矩阵背景图像的折射角矩阵 ${\mathrm{\Φ}}_3^B,{\mathrm{\Φ}}_4^B,{\mathrm{\Φ}}_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，工作管电压为40KV，管电流为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。需要指出的是基于本发明的这个优势在这个实施例上没有得到体现，但本领域的技术人员可以理解。

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

Claims (3)

1.一种用于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的每一个元素

2.如权利要求1所述的用于X射线光栅相位衬度成像的背景扣除方法，其特征在于，所述步骤S1中，每次移动分析光栅的距离均为d/M，d是所述分析光栅的周期，M是上述过程中总共采集图像的数量。

3.如权利要求1所述的用于X射线光栅相位衬度成像的背景扣除方法，其特征在于，M大于或等于3。