CN106596594B - 一种基于成像系统特性的x射线相位成像方法 - Google Patents

Abstract

本发明公开一种基于成像系统特性的X射线相位成像方法，包括：设置X射线类同轴相衬成像参数和数字放射成像系统的曝光参数，在待成像的物平面位置放置刀口器具，采集获得刀口图像，从图像中获取不同位置的刀口截面曲线，计算二维点扩散函数。放置成像物体，对物体成像。根据物体组成的先验信息以及所得图像的噪声情况，计算松弛松弛因子的初始值。计算物体的相位信息。根据所得结果的对比度，微调松弛因子，直到获得理想结果。

Description

一种基于成像系统特性的X射线相位成像方法

技术领域

本发明属于生物医学工程及医学影像学领域，涉及一种基于成像系统特性的X射线相位成像的方法。

背景技术

乳腺癌是当前妇女第一大“杀手”，据2015年最新统计数据，全世界乳腺癌患者人数和死亡人数在2012年分别增加至170万和52万，各占所有女性癌症病例的25％和癌症死亡数的15％[2]。我国人口众多，近年来乳腺癌患病率呈明显上升趋势，乳腺癌的诊治已经成为日益沉重并亟需解决的社会问题，而实现其早期诊断是解决这一社会问题、提高患者生存率和生活质量的关键。

当前乳腺常规检查的主要手段为钼(铑)靶X射线乳房成像术，然而长期临床实践表明，该技术在灵敏度，特异性，安全性和舒适性等方面均存在重大缺欠：一方面该技术存在高达10-15％的漏检率；另一方面由该技术诊断为阳性而最终的活检确诊率为25-29％。尤其严重的是，由于年轻女性的乳房过于致密，导致诊断准确性严重降低。

直到上世纪末，X射线相位衬度成像理论(X-ray phase contrast imaging,XPCI)的提出，打破了传统的X射线成像理念，为实现理想的早期乳腺癌解剖成像诊断技术带来了新的曙光。研究表明，在相同辐射剂量下，相衬成像的衬度分辨率较之传统X线吸收衬度成像提高10倍左右，显著提高了软组织成像的图像可见度。在当前可用的X射线相衬成像技术中，类同轴相衬成像技术由于不需要引入额外的光学装置，成像光路设计相对简单稳定，被认为是当前条件下最适合实现临床医学应用转化的显微成像技术之一。

在类同轴相衬成像测量中，入射的X射线经过物体后，波面发生畸变，当畸变后的波面在继续传播到一定距离后，将和未发生畸变的波面重叠而发生干涉效应。这样，X光通过薄层弱吸收弱相移样本后，经过一定距离的自由传播，就能将相位改变信息转化成经调制后的强度信息通过像平面显示出来。类同轴相衬成像成像采用的采用菲涅耳近场衍射的系统参数设置，此时在像面上获得的衍射强度分布正比于相位改变量的Laplacian变换，并非直接反映的相位改变二维分布，所以要通过数学方法从得到的光强信息中获取相位分布信息，这种由强度测量来恢复物体相位信息问题的逆过程就是相位抽取。到目前为止，已经有很多学者开展了理想情况下类同轴相衬成像的相位抽取研究，取得了有价值的成果。

但迄今为止，X射线类同轴相衬成像在乳腺癌临床诊断上仍没有得到应用，大部分的研究是基于同步辐射源的类同轴相衬成像展开的，同步辐射源投资巨大，占地面积巨大，不利于临床应用推广。而实际工程可行的基于微焦点源的X射线相衬成像系统，由于其系统非理想性，相衬成像与理论结果存在差异，导致最终相位抽取结果并没有达到预期的结果。目前采用的各类类同轴相衬成像的相位抽取模型是建立在理想的成像系统和成像过程上，而对于实际系统中的引起成像质量恶化的因素没有充分考虑，更没有采取措施对成像结果中的恶化效应进行排除。虽然通过采用高亮度微焦斑X射线源以及高量子检测效率探测器可望在较大程度上改善上述问题，但在目前条件下，高性能X射线源及探测器尚处于研发阶段，即使这些设备得到了应用，成像质量恶化因素仍然无法完全排除。另外，如果引入高性能设备，势必造成系统成本提高，导致类同轴相衬成像的推广普及优势的极大削弱。

因此，当前制约同轴相衬成像技术在临床推广应用的关键问题主要体现在像系统方面，图像系统自身存在着缺陷，比如X光源并非理想点源，探测器性能受到自身分辨率及点扩散函数等因素的限制，系统存在各类有害噪声等。

由于基于微焦点源的类同轴相衬成像是更适合实现临床推广应用的成像技术，因此，针对该成像技术实现相位信息抽取的准确度，获取表征物体内部结构信息的相位图像具有显著的研究和应用价值。为了在成像系统存在缺陷情况下，实现物体结构真实相位的获取，必须充分考虑到成像系统自身的特性，发展出一种基于成像系统特性的X射线相位成像方法，这对于类同轴相衬成像系统的发展应用具有十分重要的意义。

发明内容

本发明的主旨是提出一种基于成像系统特性的X射线相位成像方法，以此解决当前工程技术条件下，基于类同轴相衬成像的相位信息抽取所面临的关键问题：由于当前工程条件下，X光源，探测器的非理想性以及系统噪声等引起了相衬图像质量恶化，相衬度降低，基于传统的相位抽取方法无法保证相位信息抽取结果的准确性。本发明的基于成像系统特性的X射线相位成像方法，保证了非理想成像系统下准确获取被成像物体的相位信息，为类同轴相衬成像系统及技术的临床推广应用提供有力条件，从而为实现乳腺癌的早期诊断提供技术支持。技术方案如下：

一种基于成像系统特性的X射线相位成像方法，包括下列步骤：

①设置X射线类同轴相衬成像参数：设定光源到物体的距离，并设定物体到探测器的距离d。

②设置数字放射成像系统的曝光参数，在待成像的物平面位置放置刀口器具，采集获得刀口图像，从图像中获取不同位置的刀口截面曲线，而后将刀口截面曲线进行平均，再对平均曲线求导数，获得对应的表征系统特性的线扩散函数，将该线扩散函数在直角坐标系中旋转一周，计算得到二维点扩散函数PSF(x,y),(x,y)是探测器平面上的坐标位置。

③放置成像物体，保持成像参数不变，对物体成像，获得成像结果I_d(x,y)。

④根据物体组成的先验信息以及所得图像的噪声情况，计算松弛松弛因子α的初始值其中λ为X射线的波长，δ为被成像物体的复折射率n＝1-δ+jβ中的相移部分，j为虚数单位，m_u与复折射率n中的吸收因子β相关，SNR为实际获得信号的信噪比。

⑤计算g(I)函数，其中I_in(x,y)选择无物体部分的X射线强度平均值。

⑥利用公式计算物体的相位信息，这里u，v为探测器平面所在直角坐标系下两个方向的空间频率，FT为傅里叶变换，FT^-1为傅里叶逆变换。

⑦根据所得结果的对比度，微调松弛因子，直到获得理想结果。

附图说明

图1类同轴相衬成像系统的线扩散函数

图2仿体的真实相位剖面曲线

图3仿体在理想情况下的相衬成像结果的剖面曲线

图4仿体在实际系统下的相衬成像结果的剖面曲线

图5传统的相位抽取方法获得的相位剖面曲线

图6本发明相位抽取方法获得的相位剖面曲线

具体实施方式

下面结合附图和实施例对本发明进行说明。

1数字X射线成像系统

实验成像系统为Pixarray 100小动物数字放射成像系统，由美国BIOPTICS公司制造。该系统的探测器为1024×1024的CCD阵列，像素大小为50μm×50μm，14级灰度。横向及纵向的空间分辨率均为每毫米20像素。X射线管的焦斑尺寸为50μm。实验中，X射线源的工作电压为33kVp，工作电流为0.5mA。成像物体采用直径分别为100μm，50μm，20μm和10μm的聚乙烯纤维构建。实验设置X射线源到物体的距离为100cm，对应的物体到探测器的距离为100cm。在以上设置下，光源在探测器上成的焦斑像为50μm。由于成像系统的点扩散函数是探测器点扩散函数和光源点扩散函数的卷积，在上述实验设置情况下，需要考虑探测器和光源焦斑点扩散函数对类同轴相衬成像结果的恶化效应以及系统噪声问题，最终针对程序系统特性，实现X射线相位成像方法。

2类同轴相衬成像的相位抽取原理

在X射线类同轴相衬成像中，物体对入射X射线的传递函数T(x,y)可以表示为

其中(x,y)是物平面上的坐标位置，j为虚数单位，μ(x,y)和分别对应物体对入射X射线的吸收作用和相移效应。

X射线穿过物体后，相关光在继续传播过程中发生Fresnel衍射效应，X射线波面表示为：

W(x,y)＝T(x,y)*P_d(x,y) (2)

其中下标d表示物体到探测器的距离，*为线性卷积算子，P_d(x,y)是Fresnel衍射算子：

这样，由最终探测器采集到的图像为(假设入射X射线强度为均匀分布)

I_d(x,y)＝|T(x,y)*P_d(x,y)|² (4)

将类同轴相衬成像的近轴条件代入上式，可以得到探测器上光强和物体相位信息的关系式：

实际成像时，由于系统探测器光源的非理想性以及系统噪声，在探测器平面上获得的真实相衬成像结果需要将这些因素考虑进来，即像面上光强为

这里PSF(x,y)是光源和探测器点扩散函数(Point Spread Function,PSF)的综合效应，即光源焦斑投影到探测器平面上的点扩散函数与探测器自身点扩散函数的卷积结果，N(x,y)为系统噪声。

为了推导出相位信息表达式，将上式修改为

将上式两端进行傅里叶变换后整理，得到频域表达式

这里FT表示傅里叶变换算子。

令利用拉普拉斯算子项的求解原理，引入滤波项H_p(u,v)((u,v)为探测器平面对应的直角坐标系空间频率坐标)，可以得到相位信息表达式

这里FT-1为傅里叶逆变换，由于H_p(u,v)在零频率位置造成无穷大，而FT[PSF(x,y)]也存在某些频率位置幅度为0的可能，为了解决这个问题，需要引入松弛因子α，而松弛因子的引入，可以将系统原有的噪声问题也实现缓解，则可以将g(I)函数的噪声项忽略，设置新的g(I)为：

对应的相位求解的最终计算公式为：

这里u，v为探测器平面所在直角坐标系下两个方向的空间频率，FT^-1为傅里叶逆变换，松弛因子α的初始值可以按照下列公式计算：

其中λ为X射线的波长，δ为被成像物体的复折射率(n＝1-δ+jβ，j为虚数单位)中的相移部分，m_u与复折射率n＝1-δ+jβ中的吸收因子β相关，SNR为实际获得信号的信噪比。最佳α值的选择取决于物体的吸收及图像的噪声：物体吸收强时需要用大的α来补偿；大的α值能够更好地抑制噪声但是会导致信号失真，而小的α值虽然能得到更准确地信号，但是其代价是引入了更多噪声。

3本发明的一种基于成像系统特性的X射线相位成像方法的流程描述如下：

1)X射线类同轴相衬成像参数设置：本发明中，设置光源到物体的距离为100cm，与此对应的物体到探测器的距离为100cm。

2)设置数字放射成像系统的曝光参数，在待成像的物平面位置放置刀口器具，采集获得刀口图像，从图像中获取不同位置的刀口截面曲线，而后将刀口截面曲线进行平均，再对平均曲线求导数，获得对应的表征系统特性的线扩散函数，将该线扩散函数在直角坐标系中旋转一周，计算得到二维点扩散函数PSF(x,y)，(x,y)是探测器平面上的坐标位置。

3)放置成像物体，本发明中采用直径分别为100μm，50μm，20μm和10μm的聚乙烯纤维构建。保持成像参数不变，对物体成像，获得成像结果I_d(x,y)。

4)根据物体组成的先验信息以及所得图像的噪声情况，计算松弛松弛因子α的初始值其中λ为X射线的波长，δ为被成像物体的复折射率(n＝1-δ+jβ，j为虚数单位)中的相移部分，m_u与复折射率n＝1-δ+jβ中的吸收因子β相关，SNR为实际获得信号的信噪比。

5)计算g(I)函数，其中I_in(x,y)选择无物体部分的X射线强度平均值。

6)利用公式计算物体的相位信息，这里u，v为探测器平面所在直角坐标系下两个方向的空间频率，FT^-1为傅里叶逆变换。根据所得结果的对比度，微调松弛因子，直到获得理想结果。

本实施例采用美国BIOPTICS公司生产的Pixarray 100小动物数字放射成像装置构建类同轴相衬成像系统。首先通过刀口装置获取系统探测器的传递函数。图1给出了通过刀口法测量获得的系统线扩散函数，由此测量获得的系统线函数曲线的半高宽是80微米。将该线扩散函数在直角坐标系中旋转一周，可以获得对应的点扩散函数。

根据仿体中采用纤维的物质组成以及其直径，可以推算出仿体的真实相位，图2给出了仿体中4根不同直径纤维真实相位的剖面曲线。

采用数值模拟方法，可以得到仿体在理想成像系统下对应的相衬成像结果，图3为理想相衬成像结果的剖面曲线。

当考虑光源和探测器的非理想性以及系统噪声时，我们通过所构建的成像系统，获得的相衬成像结果并取对应的剖面曲线，结果如图4所示。

对实际相衬成像结果采用传统的相位抽取方法，即忽略成像系统的非理想性，可以得到图5所示的相位抽取结果。可以看到对直径为100μm，50μm和20μm的纤维，抽取的相位信息低于真实相位，而对于直径10μm的纤维，已经很难获得准确的相位信息。

采用本发明提出的基于成像系统特性的X射线相位成像方法，考虑成像系统光源探测器非理想性以及系统噪声问题，通过微调松弛因子，可以获得如图6所示的相位成像结果，逼近图2推算得到的仿体真实相位，优于传统的相位抽取方法。

最终结果表明，针对当前工程条件下，类同轴相衬成像系统成像结果的恶化效应，采用本发明的基于成像系统特性的X射线相位成像方法，可以有效提高非理想系统情况下，相位抽取结果的准确性。由于本发明实现采用的物体为聚乙烯材料，其折射率与人体乳腺组织非常接近，因此本发明的研究成果可以进一步有效的一指到实际临床上用于乳腺癌早期诊断的同轴相衬成像系统中。该方法的应用，将为有效实现早期乳腺癌微小病变组织的诊断提供技术支持，为深入开展乳腺癌的同轴相衬成像的临床实践和研究提供有力支持。

Claims (1)

1.一种基于成像系统特性的X射线相位成像方法，包括下列步骤：

①设置X射线类同轴相衬成像参数：设定光源到物体的距离，并设定物体到探测器的距离d；

②设置数字放射成像系统的曝光参数，在待成像的物平面位置放置刀口器具，采集获得刀口图像，从图像中获取不同位置的刀口截面曲线，而后将刀口截面曲线进行平均，再对平均曲线求导数，获得对应的表征系统特性的线扩散函数，将该线扩散函数在直角坐标系中旋转一周，计算得到二维点扩散函数PSF(x,y),(x,y)是探测器平面上的坐标位置；

③放置成像物体，保持成像参数不变，对物体成像，获得成像结果I_d(x,y)；

④根据物体组成的先验信息以及所得图像的噪声情况，计算松弛松弛因子α的初始值其中λ为X射线的波长，δ为被成像物体的复折射率n＝1-δ+jβ中的相移部分，j为虚数单位，m_u与复折射率n中的吸收因子β相关，SNR为实际获得信号的信噪比；

⑤计算g(I)函数，其中I_in(x,y)选择无物体部分的X射线强度平均值；

⑥利用公式计算物体的相位信息，这里u，v为探测器平面所在直角坐标系下两个方向的空间频率，FT为傅里叶变换，FT^-1为傅里叶逆变换；

⑦根据所得结果的对比度，微调松弛因子，直到获得理想结果。