CN103559698A - 一种基于混合迭代的同轴相衬成像相位恢复方法及系统 - Google Patents

Abstract

本发明公开了一种基于混合迭代的同轴相衬成像相位恢复方法及系统，该方法包括：采集两个位置的物体图像和背景图像；根据物体图像和背景图像，将两个位置的物体图像进行归一化处理，选择归一化后的一个位置的物体图像为参照物体图像、另一个位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；根据第一迭代算法、第三物体图像和参照物体图像，计算参照物体图像对应位置的初始收敛相位分布；根据第二迭代算法、第三物体图像、参照物体图像以及参照物体图像对应位置的初始收敛相位分布，计算参照物体图像对应位置的最终收敛相位分布。实施本发明实施例，可以增强相位恢复的精确度、效率和稳定性。

Description

一种基于混合迭代的同轴相衬成像相位恢复方法及系统

技术领域

本发明涉及光学技术领域，具体涉及一种基于混合迭代的同轴相衬成像相位恢复方法及系统。

背景技术

X射线相位衬度成像（X-ray phase contrast imaging，XPCI）是一种X射线成像技术，利用X射线通过物体后发生的相移变化成像。基于不同的成像原理，目前有五种实现X射线相位衬度成像的技术，分别为：干涉法、衍射增强法、光栅剪切法、同轴相衬法以及编码孔径相衬成像法。同轴相衬成像装置最为简单，无需精密的光学元器件，可以基于实验室普通的微焦斑X射线源，利用X射线通过样品后在自由空间传播，基于菲涅耳衍射成像原理将相位信息转化为强度信息。

同轴相衬成像可以增强图像的边缘亮度，从而提高了图像的对比度，但这种亮度的增强不是线性的，可能会呈现出错误的样品厚度或密度。然而，可以通过算法从获得的强度图中恢复出相位分布图，正确反应出样品的真实结构与组织特性，但对基于普通微焦斑X射线源的同轴相衬成像，想实现精确相位恢复非常困难。目前，同轴相衬成像的相位恢复方法有：解析算法（或线性近似算法）和迭代算法。解析算法是将非线性方程进行线性近似得到方程的解，计算效率高，但相位的解可能不稳定。此外，由于解析算法基于一定的假设和简化，不同的方法适用的范围不同、受限于成像物体组成或者受限于成像距离。迭代算法由于近似少，适用范围更广，且简单、灵活、结果稳定、精确。其中，基于傅里叶变换（Fourier transform，FT）的迭代算法的效率高，但精度低，而基于盖师贝格-撒克斯通（Gerchberg and Saxton）算法的精度高，但效率低。

发明内容

本发明公开了一种基于混合迭代的同轴相衬成像相位恢复方法及系统，用于增强相位恢复的精确度、效率和稳定性。

本发明第一方面公开了一种基于混合迭代的同轴相衬成像相位恢复方法，包括：

采集第一位置的物体图像以及所述第一位置的背景图像；

采集第二位置的物体图像以及所述第二位置的背景图像；

根据所述第一位置的物体图像和所述第一位置的背景图像，以及所述第二位置的物体图像和所述第二位置的背景图像，将采集的所述第一位置的物体图像和所述第二位置的物体图像进行归一化处理，选择归一化处理后的所述第一位置的物体图像为参照物体图像，选择归一化处理后的所述第二位置的物体图像为对照物体图像，将所述参照物体图像和所述对照物体图像进行图像精确配准，获得第三物体图像；

根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布；

根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布。

本发明第二方面公开了一种基于混合迭代的同轴相衬成像相位恢复系统，包括：

第一单元，用于采集第一位置的物体图像以及所述第一位置的背景图像；

所述第一单元，还用于采集第二位置的物体图像以及所述第二位置的背景图像；

第二单元，用于根据所述第一位置的物体图像和所述第一位置的背景图像，以及所述第二位置的物体图像和所述第二位置的背景图像，将采集的所述第一位置的物体图像和所述第二位置的物体图像进行归一化处理，选择归一化处理后的所述第一位置的物体图像为参照物体图像，选择归一化处理后的所述第二位置的物体图像为对照物体图像，将所述参照物体图像和所述对照物体图像进行图像精确配准，获得第三物体图像；

第三单元，用于根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布；

第四单元，用于根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布。

本发明实施例中，采集第一位置的物体图像以及第一位置的背景图像；并采集第二位置的物体图像以及第二位置的背景图像；根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布；根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。本发明实施例中，在进行迭代计算前，先对采集的物体图像进行了图像精确配准，在进行迭代计算时，先用第一种迭代算法进行计算收敛相位分布，然后用该收敛相位分布为初始值继续用第二种迭代算法进行计算得到最终的收敛相位分布，可以增强相位恢复的精确度、效率和稳定性。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明第一实施例公开的一种基于混合迭代的同轴相衬成像相位恢复方法的流程图；

图2是本发明第二实施例公开的另一种基于混合迭代的同轴相衬成像相位恢复方法的流程图；

图3是本发明第三实施例公开的一种基于混合迭代的同轴相衬成像相位恢复系统的结构图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供了一种基于混合迭代的同轴相衬成像相位恢复方法及系统，用于增强相位恢复的精确度、效率和稳定性。以下分别进行详细说明。

请参阅图1，图1是本发明第一实施例公开的一种基于混合迭代的同轴相衬成像相位恢复方法的流程图。其中，图1所示的基于混合迭代的同轴相衬成像相位恢复方法适用于X射线相衬成像系统。如图1所示，该基于混合迭代的同轴相衬成像相位恢复方法可以包括以下步骤。

S101、采集第一位置的物体图像以及第一位置的背景图像。

本发明实施例中，X射线相衬成像系统采集第一位置的物体图像以及第一位置的背景图像。

本发明实施例中，采集的背景图像为无物体放置的亮场图像。

S102、采集第二位置的物体图像以及第二位置的背景图像。

本发明实施例中，X射线相衬成像系统采集第二位置的物体图像以及第二位置的背景图像。

本发明实施例中，第二位置与第一位置是两个不同的位置，在采集图像时，X射线图像探测器与物体间的距离不同。

S103、根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像。

本发明实施例中，X射线相衬成像系统根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像。

本发明实施例中，一般需要采集两幅物体图像，一幅X射线图像探测器与物体近贴的吸收图像、一幅X射线图像探测器与物体具有一定间隔的相衬图像，或者直接采集两幅相衬图像。

本发明实施例中，采集了两个不同位置的物体图像和相应位置的背景图像后，首先需要对物体图像进行归一化处理，即用采集的物体图像除以物体图像对应位置的背景图像，得到归一化处理后的物体图像。

本发明实施例中，在对物体图像进行迭代计算以前，还需要对物体图像进行图像精确配准预处理，这是由于在不同位置采集的图像其几何放大倍数不同，而且系统自身也存在一定的机械精度，应用中采集的两幅物体图像存在缩放、平移和旋转，在进行相位恢复前对物体图像的精确配准，可以提高相位恢复精度。采用基于傅氏变换的方法完成图像精确配准，相对基于灰度或基于特征的配准算法，更方便处理采集的物体图像。因为图像的比例、旋转和平移变换均能在傅里叶变换频域中反映出来，且在频域内对噪声干扰有一定的抵抗能力，同时，傅里叶变换可以采用快速傅里叶变换（FFT）方法提高处理速度，并且有成熟的快速算法，易于硬件实现。实现图像精确配准的步骤如下：

1）选择归一化后的物体图像A为参照物体图像，归一化后的物体图像B为对照物体图像，对参照物体图像A与对照物体图像B分别做傅里叶变换，得到新的参照物体图像

与新的对照物体图像

2）在傅里叶空间对新的参照物体图像

与新的对照物体图像

分别进行对数极坐标变换；

3）再在对数极坐标下利用相位相关求得对照物体图像B相对于参照物体图像A的缩放因子a和旋转角度θ₀；

4）根据缩放因子a和旋转角度θ₀，对对照物体图像B进行角度和缩放比例补偿，补偿以后得到的物体图像B₂与参照物体图像A之间仅存在平移量的差别；

5）利用直角坐标下的相位相关求出对照物体图像B相对于参照物体图像A的平移量；

6）根据平移量，对物体图像B₂进行平移补偿，得到新的物体图像C。

S104、根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布。

本发明实施例中，X射线相衬成像系统根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布。

本发明实施例中，第一迭代算法为傅里叶变换算法。傅里叶变换（Fouriertransform,FT）迭代算法的原理：

基于Wu和Liu提出的一般性相衬成像公式[Xizeng Wu and Hong Liu,Ageneral theoretical formalism for X-rayphase contrast imaging,2003]，公式可写为：

$\begin{array}{c}\mathrm{FT}\left[M^{2}I(M\stackrel{\→}{\mathrm{\ρ}};R_2)\right]=I_0\·{\mathrm{OTF}}_{G.U.}\left(\frac{\stackrel{\→}{u}}{M}\right)\·{\mathrm{OTF}}_{\det }\left(\frac{\stackrel{\→}{u}}{M}\right)\\ \left\{\cos \left({\mathrm{\αu}}^2\right)\mathrm{FT}\right[A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]+2\sin \left({\mathrm{\αu}}^2\right)\mathrm{FT}[A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)\mathrm{\φ}\left(\stackrel{\→}{\mathrm{\ρ}}\right)]\\ +i\frac{{\mathrm{\λR}}_2}{2M}\sin \left({\mathrm{\αu}}^2\right)\stackrel{\→}{u}\·\mathrm{FT}[\▿A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]+i\frac{{\mathrm{\λR}}_2}{M}\cos \left({\mathrm{\αu}}^2\right)\stackrel{\→}{u}\·\mathrm{FT}[\mathrm{\φ}\left(\stackrel{\→}{\mathrm{\ρ}}\right)\▿A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]\}\end{array}---\left(1\right)$

其中分别表示物平面空间向量与空间频率向量，α＝πλR₂/M，FT[·]表示傅立叶变换，I₀表示入射到物平面的光强，

表示有限焦斑大小引起的几何模糊光学传输函数，表示探测器的空间频率响应；M表示几何放大率(R₁＋R₂)/R₁，R₁，R₂分别表示源物距与物像距，λ表示波长，A表示光波振幅，φ表示相位。把上面大括号中四项分别简称为T₁、T₂、T₃与T₄项。T₃很小可以忽略，把T₄视为微扰项，得到迭代公式：

$\mathrm{FT}\left[A^2\right]={\tilde{I}}_1$

$\mathrm{FT}\left[A^2{\mathrm{\φ}}^{(n+1)}\right]=\frac{{\tilde{I}}_2-\cos \left({\mathrm{\αu}}^2\right)\mathrm{FT}\left[A^2\right]-T4(A^2,{\mathrm{\φ}}^{\left(n\right)},R_2)}{2\sin \left({\mathrm{\αu}}^2\right)}---\left(2\right)$

其中

表示空间频域中归一化与去卷积后的图像，OTF表示焦斑几何模糊与探测器响应两者总的光学传输函数。式(2)即为基于傅立叶变换(FT)的迭代公式。

假设X射线源焦斑为直径为f的均匀圆分布，则焦斑几何模糊光学传输函数OTF_G.U.可写为

$\left|{\mathrm{OTF}}_{G.U.}\left(\frac{u}{M}\right)\right|=\frac{{2J}_1\left[\mathrm{\πf}(M-1)\right|u|/M]}{\mathrm{\πf}(M-1)\left|u\right|/M}---\left(3\right)$

其中J₁(x)为第一类贝塞尔函数。

目前，普遍应用的探测器是数字化成像探测器，其响应为矩形窗函数，光学传输函数OTF_det可写为

$\left|{\mathrm{OTF}}_{\det }\left(\frac{u}{M}\right)\right|=\sin c\left(\frac{p_{\mathrm{ax}}{u}_x}{M}\right)\sin c\left(\frac{p_{\mathrm{ay}}{u}_y}{M}\right)---\left(4\right)$

其中sinc(u)是sinc函数，表示矩形窗函数的傅里叶变换，p_ax,y是有效像素大小，决定于探测器像素周期与填充因子。

在实际应用中，经常由于系统几何或图像分辨要求的限制，无法采集近贴的吸收图像，而需要用不同放大率下的两个相衬像I₁与I₂进行相位恢复，假设R₁>R₂，只需将两幅图像强度分别代入公式(2)，即可得到R₁位置处的振幅与相位分布，再根据光传播的可逆性，用菲涅耳传播因子去卷积即得到物平面处的振幅与相位分布。公式(5)即为利用两幅相衬图像的相位恢复迭代公式。

$[A^2(\stackrel{\→}{\mathrm{\ρ}};R_1),\mathrm{\φ}(\stackrel{\→}{\mathrm{\ρ}};R_1)]=\mathrm{Iter}(I_1,I_2),$

$A(\stackrel{\→}{\mathrm{\ρ}};0)\exp \left(\mathrm{i\φ}(\stackrel{\→}{\mathrm{\ρ}};0)\right)=\mathrm{IFT}\left[\frac{\mathrm{FT}\left[A(\stackrel{\→}{\mathrm{\ρ}};R_1)\exp \left(\mathrm{i\φ}(\stackrel{\→}{\mathrm{\ρ}};R_1)\right)\right]}{\mathrm{FT}[\exp (-\mathrm{i\π}{\mathrm{\ρ}}^2)/\mathrm{\λ}{R}_1]}\right]---\left(5\right)$

其中Iter表示采用公式（2）的迭代运算，IFT[〃]表示反傅里叶变换。

S105、根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。

本发明实施例中，X射线相衬成像系统根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。

本发明实施例中，第二迭代算法为盖师贝格-撒克斯通算法或泊松分布算法。其中，盖师贝格-撒克斯通（Gerchberg and Saxton，GS）算法的原理：

GS算法硬件实现较为灵活，算法也比较简单，利用光传播的可逆性，在不同的像面间正反向反复迭代计算，利用采集的光强图像作为迭代限制条件。光传播为近轴菲涅尔衍射时，光波场计算公式为：

$U(x,y;z)=\frac{\mathrm{iexp}(-\mathrm{ikz})}{\mathrm{\λz}}\∫\∫U(x^{\′},y^{\′};0)\×\exp [-\mathrm{ik}\frac{{\left({x-x}^{\′}\right)}^2+{\left({y-y}^{\′}\right)}^2}{2z}]{\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′},---\left(6\right)$

其中U(x,y;z)表示光波场，k为波数。其中，主要迭代步骤如下：

1）任给一个初始相位φ⁽ⁿ⁾(x，y；0)，令 $U^{\left(n\right)}(x,y;0)\⇐\sqrt{I_1}\exp \left[{\mathrm{i\φ}}^{\left(n\right)}(x,y;0)\right];$

2）根据U⁽ⁿ⁾(x,y;0)和菲涅尔衍射公式(6)计算U⁽ⁿ⁾(x,y;R₂，取幅角得到U⁽ⁿ⁾(x,y;R₂)；

3）令 $U^{\left(n\right)}(x,y;R_2)\⇐\sqrt{I_2}\exp \left[{\mathrm{i\φ}}^{\left(n\right)}(x,y;R_2)\right];$

4）根据U⁽ⁿ⁾(x，y；R₂和菲涅尔衍射公式计算U⁽ⁿ⁺¹⁾(x,y;0)，取幅角得到φ⁽ⁿ⁺¹⁾(x,y;0)；

5）判断φ⁽ⁿ⁺¹⁾(x,y;0)的变化是否足够小，若是，则迭代结束；若否，则令

重复步骤2-5；其中，n和n+1表示迭代的次数。

在图1所示的基于混合迭代的同轴相衬成像相位恢复方法中，X射线相衬成像系统采集第一位置的物体图像以及第一位置的背景图像；并采集第二位置的物体图像以及第二位置的背景图像；根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布；根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。本发明实施例中，在进行迭代计算前，先对采集的物体图像进行了图像精确配准，在进行迭代计算时，先用第一种迭代算法进行计算收敛相位分布，然后用该收敛相位分布为初始值继续用第二种迭代算法进行计算得到最终的收敛相位分布，可以增强相位恢复的精确度、效率和稳定性。

请参阅图2，图2是本发明第二实施例公开的另一种基于混合迭代的同轴相衬成像相位恢复方法的流程图。其中，图2所示的基于混合迭代的同轴相衬成像相位恢复方法适用于X射线相衬成像系统。如图2所示，该基于混合迭代的同轴相衬成像相位恢复方法可以包括以下步骤。

S201、X射线相衬成像系统设置射线源的优化工作参数以及探测器的工作模式。

本发明实施例中，X射线成像系统在采集图像之前，首先需要设置射线源的优化工作参数以及探测器的工作模式，以保证采集的图像效果。

S202、X射线相衬成像系统采集第一位置的物体图像以及第一位置的背景图像。

本发明实施例中，采集的背景图像为无物体放置的亮场图像。

S203、X射线相衬成像系统采集第二位置的物体图像以及第二位置的背景图像。

本发明实施例中，第二位置与第一位置是两个不同的位置，在采集图像时，X射线图像探测器与物体间的距离不同。

S204、X射线相衬成像系统根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像。

本发明实施例中，一般需要采集两幅物体图像，一幅X射线图像探测器与物体近贴的吸收图像、一幅X射线图像探测器与物体具有一定间隔的相衬图像，或者直接采集两幅相衬图像。

本发明实施例中，采集了两个不同位置的物体图像和相应位置的背景图像后，首先需要对物体图像进行归一化处理，即用采集的物体图像除以物体图像对应位置的背景图像，得到归一化处理后的物体图像。

本发明实施例中，在对物体图像进行迭代计算以前，还需要对物体图像进行图像精确配准预处理，这是由于在不同位置采集的图像其几何放大倍数不同，而且系统自身也存在一定的机械精度，应用中采集的两幅物体图像存在缩放、平移和旋转，在进行相位恢复前对物体图像的精确配准，可以提高相位恢复精度。采用基于傅氏变换的方法完成图像精确配准，相对基于灰度或基于特征的配准算法，更方便处理采集的物体图像。因为图像的比例、旋转和平移变换均能在傅里叶变换频域中反映出来，且在频域内对噪声干扰有一定的抵抗能力，同时，傅里叶变换可以采用快速傅里叶变换（FFT）方法提高处理速度，并且有成熟的快速算法，易于硬件实现。实现图像精确配准的步骤如下：

1）选择归一化后的物体图像A为参照物体图像，归一化后的物体图像B为对照物体图像，对参照物体图像A与对照物体图像B分别做傅里叶变换，得到新的参照物体图像

与新的对照物体图像

2）在傅里叶空间对新的参照物体图像

与新的对照物体图像

分别进行对数极坐标变换；

3）再在对数极坐标下利用相位相关求得对照物体图像B相对于参照物体图像A的缩放因子a和旋转角度θ₀；

4）根据缩放因子a和旋转角度θ₀，对对照物体图像B进行角度和缩放比例补偿，补偿以后得到的物体图像B₂与参照物体图像A之间仅存在平移量的差别；

5）利用直角坐标下的相位相关求出对照物体图像B相对于参照物体图像A的平移量；

6）根据平移量，对物体图像B2进行平移补偿，得到新的物体图像C。

S205、X射线相衬成像系统计算焦斑几何模糊光学传输函数。

本发明实施例中，在根据迭代算法计算参照物体图像位置的收敛相位分布时，需要焦斑几何模糊光学传输函数。因此，在根据迭代算法计算参照物体图像位置的收敛相位分布之前，需要计算焦斑几何模糊光学传输函数。

S206、X射线相衬成像系统计算探测器响应的光学传输函数。

本发明实施例中，在根据迭代算法计算参照物体图像位置的收敛相位分布时，需要探测器响应的光学传输函数。因此，在根据迭代算法计算参照物体图像位置的收敛相位分布之前，需要计算探测器响应的光学传输函数。

S207、X射线相衬成像系统将参照物体图像和第三物体图像进行去卷积处理。

本发明实施例中，在根据迭代算法计算参照物体图像位置的收敛相位分布之前，需要对参照物体图像和第三物体图像进行去卷积预处理。

S208、X射线相衬成像系统设置第一位置的初始相位分布。

本发明实施例中，在根据迭代算法计算参照物体图像位置的收敛相位分布时，需要一个初始相位分布。因此，在根据迭代算法计算参照物体图像位置的收敛相位分布之前，首先给参照物体图像位置设置一个初始相位分布，如全0分布。

S209、X射线相衬成像系统根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布。

本发明实施例中，第一迭代算法为傅里叶变换算法。傅里叶变换（Fouriertransform,FT）迭代算法的原理：

基于Wu和Liu提出的一般性相衬成像公式[Xizeng Wu and Hong Liu,Ageneral theoretical formalism for X-rayphase contrast imaging,2003]，公式可写为：

$\begin{array}{c}\mathrm{FT}\left[M^{2}I(M\stackrel{\→}{\mathrm{\ρ}};R_2)\right]=I_0\·{\mathrm{OTF}}_{G.U.}\left(\frac{\stackrel{\→}{u}}{M}\right)\·{\mathrm{OTF}}_{\det }\left(\frac{\stackrel{\→}{u}}{M}\right)\\ \left\{\cos \left({\mathrm{\αu}}^2\right)\mathrm{FT}\right[A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]+2\sin \left({\mathrm{\αu}}^2\right)\mathrm{FT}[A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)\mathrm{\φ}\left(\stackrel{\→}{\mathrm{\ρ}}\right)]\\ +i\frac{{\mathrm{\λR}}_2}{2M}\sin \left({\mathrm{\αu}}^2\right)\stackrel{\→}{u}\·\mathrm{FT}[\▿A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]+i\frac{{\mathrm{\λR}}_2}{M}\cos \left({\mathrm{\αu}}^2\right)\stackrel{\→}{u}\·\mathrm{FT}[\mathrm{\φ}\left(\stackrel{\→}{\mathrm{\ρ}}\right)\▿A^2\left(\stackrel{\→}{\mathrm{\ρ}}\right)]\}\end{array}---\left(1\right)$

其中分别表示物平面空间向量与空间频率向量，α＝πλR₂/M，FT[·]表示傅立叶变换，I₀表示入射到物平面的光强，

表示有限焦斑大小引起的几何模糊光学传输函数，

表示探测器的空间频率响应；M表示几何放大率(R₁＋R₂)/R₁，R₁，R₂分别表示源物距与物像距，λ表示波长，A表示光波振幅，φ表示相位。把上面大括号中四项分别简称为T₁、T₂、T₃与T₄项。T₃很小可以忽略，把T₄视为微扰项，得到迭代公式：

$\mathrm{FT}\left[A^2\right]={\tilde{I}}_1$

$\mathrm{FT}\left[A^2{\mathrm{\φ}}^{(n+1)}\right]=\frac{{\tilde{I}}_2-\cos \left({\mathrm{\αu}}^2\right)\mathrm{FT}\left[A^2\right]-T4(A^2,{\mathrm{\φ}}^{\left(n\right)},R_2)}{2\sin \left({\mathrm{\αu}}^2\right)}---\left(2\right)$

其中表示空间频域中归一化与去卷积后的图像，OTF表示焦斑几何模糊与探测器响应两者总的光学传输函数。式(2)即为基于傅立叶变换(FT)的迭代公式。

假设X射线源焦斑为直径为f的均匀圆分布，则焦斑几何模糊光学传输函数OTF_G.U.可写为

$\left|{\mathrm{OTF}}_{G.U.}\left(\frac{u}{M}\right)\right|=\frac{{2J}_1\left[\mathrm{\πf}(M-1)\right|u|/M]}{\mathrm{\πf}(M-1)\left|u\right|/M}---\left(3\right)$

其中J₁(x)为第一类贝塞尔函数。

目前，普遍应用的探测器是数字化成像探测器，其响应为矩形窗函数，光学传输函数OTF_det可写为

$\left|{\mathrm{OTF}}_{\det }\left(\frac{u}{M}\right)\right|=\sin c\left(\frac{p_{\mathrm{ax}}{u}_x}{M}\right)\sin c\left(\frac{p_{\mathrm{ay}}{u}_y}{M}\right)---\left(4\right)$

其中sinc(u)是sinc函数，表示矩形窗函数的傅里叶变换，p_ax,y是有效像素大小，决定于探测器像素周期与填充因子。

在实际应用中，经常由于系统几何或图像分辨要求的限制，无法采集近贴的吸收图像，而需要用不同放大率下的两个相衬像I₁与I₂进行相位恢复，假设R₁>R₂，只需将两幅图像强度分别代入公式(2)，即可得到R₁位置处的振幅与相位分布，再根据光传播的可逆性，用菲涅耳传播因子去卷积即得到物平面处的振幅与相位分布。公式(5)即为利用两幅相衬图像的相位恢复迭代公式。

$[A^2(\stackrel{\→}{\mathrm{\ρ}};R_1),\mathrm{\φ}(\stackrel{\→}{\mathrm{\ρ}};R_1)]=\mathrm{Iter}(I_1,I_2),$

$A(\stackrel{\→}{\mathrm{\ρ}};0)\exp \left(\mathrm{i\φ}(\stackrel{\→}{\mathrm{\ρ}};0)\right)=\mathrm{IFT}\left[\frac{\mathrm{FT}\left[A(\stackrel{\→}{\mathrm{\ρ}};R_1)\exp \left(\mathrm{i\φ}(\stackrel{\→}{\mathrm{\ρ}};R_1)\right)\right]}{\mathrm{FT}[\exp (-\mathrm{i\π}{\mathrm{\ρ}}^2)/\mathrm{\λ}{R}_1]}\right]---\left(5\right)$

其中Iter表示采用公式（2）的迭代运算，IFT[·]表示反傅里叶变换。

S210、X射线相衬成像系统根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。

本发明实施例中，第二迭代算法为盖师贝格-撒克斯通算法或泊松分布算法。其中，盖师贝格-撒克斯通（Gerchberg and Saxton，GS）算法的原理：

GS算法硬件实现较为灵活，算法也比较简单，利用光传播的可逆性，在不同的像面间正反向反复迭代计算，利用采集的光强图像作为迭代限制条件。光传播为近轴菲涅尔衍射时，光波场计算公式为：

$U(x,y,z)=\frac{\mathrm{iexp}(-\mathrm{ikz})}{\mathrm{\λz}}\∫\∫U(x^{\′},y^{\′};0)\×\exp [-\mathrm{ik}\frac{{\left({x-x}^{\′}\right)}^2+{\left({y-y}^{\′}\right)}^2}{2z}]{\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′},---\left(6\right)$

其中U(x,y;z)表示光波场，k为波数。其中，主要迭代步骤如下：

1）任给一个初始相位φ⁽ⁿ⁾(x，y；0)令 $U^{\left(n\right)}(x,y;0)\⇐\sqrt{I_1}\exp \left[{\mathrm{i\φ}}^{\left(n\right)}(x,y;0)\right];$

2）根据U⁽ⁿ⁾(x,y;0)和菲涅尔衍射公式(6)计算U⁽ⁿ⁾(x,y;R₂，取幅角得到U⁽ⁿ⁾(x,y;R₂)；

3）令 $U^{\left(n\right)}(x,y;R_2)\⇐\sqrt{I_2}\exp \left[{\mathrm{i\φ}}^{\left(n\right)}(x,y;R_2)\right];$

4）根据U⁽ⁿ⁾(x,y;R₂和菲涅尔衍射公式计算U⁽ⁿ⁺¹⁾(x,y;0)，取幅角得到φ⁽ⁿ⁺¹⁾(x,y;0)；

5）判断φ⁽ⁿ⁺¹⁾(x,y;0)的变化是否足够小，若是，则迭代结束；若否，则令

重复步骤2-5；其中，n和n+1表示迭代的次数。

S211、X射线相衬成像系统判断第一位置的最终收敛相位分布是否在物平面，若否，则用菲涅尔传播因子对第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。

本发明实施例中，需要的收敛相位分布是物平面的收敛相位分布，因此，需要判断得到的第一位置的最终收敛相位分布是否在物平面，若否，则用菲涅尔传播因子对第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。若第一位置的最终收敛相位分布在物平面，则结束。

在图2所示的基于混合迭代的同轴相衬成像相位恢复方法中，X射线相衬成像系统首先设置射线源的优化工作参数以及探测器的工作模式；之后采集第一位置的物体图像以及第一位置的背景图像；并采集第二位置的物体图像以及第二位置的背景图像；根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；计算焦斑几何模糊光学传输函数；计算探测器响应的光学传输函数；将参照物体图像和第三物体图像进行去卷积处理；设置第一位置的初始相位分布；根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布；根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布；判断第一位置的最终收敛相位分布是否在物平面，若否，则用菲涅尔传播因子对第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。本发明实施例中，在进行迭代计算前，先对采集的物体图像进行了图像精确配准，在进行迭代计算时，先用第一种迭代算法进行计算收敛相位分布，然后用该收敛相位分布为初始值继续用第二种迭代算法进行计算得到最终的收敛相位分布，可以增强相位恢复的精确度、效率和稳定性。

请参阅图3，图3是本发明第三实施例公开的一种基于混合迭代的同轴相衬成像相位恢复系统的结构图。其中，图3所示的基于混合迭代的同轴相衬成像相位恢复系统可以是独立的系统，也可以是X射线相衬成像系统，本发明在此不作限定。如图3所示，该基于混合迭代的同轴相衬成像相位恢复系统300可以包括：

第一单元301，用于采集第一位置的物体图像以及第一位置的背景图像；

第一单元301，还用于采集第二位置的物体图像以及第二位置的背景图像；

第二单元302，用于根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；

第三单元303，用于根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布；

第四单元304，用于根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布。

作为一种可能的实施方式，基于混合迭代的同轴相衬成像相位恢复系统300还可以包括：

第五单元305，用于在第一单元301采集第一位置的物体图像以及第一位置的背景图像之前，设置射线源的优化工作参数以及探测器的工作模式。

作为一种可能的实施方式，基于混合迭代的同轴相衬成像相位恢复系统300还可以包括：

第六单元306，用于在第三单元303根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布之前，计算焦斑几何模糊光学传输函数；

第七单元307，用于计算探测器响应的光学传输函数；

第八单元308，用于将参照物体图像和第三物体图像进行去卷积处理；

第九单元309，用于设置第一位置的初始相位分布。

作为一种可能的实施方式，基于混合迭代的同轴相衬成像相位恢复系统300还可以包括：

第十单元310，用于在第四单元304根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布之后，判断第一位置的最终收敛相位分布是否在物平面；

第十一单元311，用于在第十单元310的判断结果为否时，用菲涅尔传播因子对第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。

本发明实施例中，第一迭代算法为傅里叶变换算法，第二迭代算法为盖师贝格-撒克斯通算法或泊松分布算法。

在图3所示的基于混合迭代的同轴相衬成像相位恢复系统中，首先设置射线源的优化工作参数以及探测器的工作模式；之后采集第一位置的物体图像以及第一位置的背景图像；并采集第二位置的物体图像以及第二位置的背景图像；根据第一位置的物体图像和第一位置的背景图像，以及第二位置的物体图像和第二位置的背景图像，将采集的第一位置的物体图像和第二位置的物体图像进行归一化处理，选择归一化处理后的第一位置的物体图像为参照物体图像，选择归一化处理后的第二位置的物体图像为对照物体图像，将参照物体图像和对照物体图像进行图像精确配准，获得第三物体图像；计算焦斑几何模糊光学传输函数；计算探测器响应的光学传输函数；将参照物体图像和第三物体图像进行去卷积处理；设置第一位置的初始相位分布；根据第一迭代算法、第三物体图像和参照物体图像，计算第一位置的初始收敛相位分布；根据第二迭代算法、第三物体图像、参照物体图像以及第一位置的初始收敛相位分布，计算第一位置的最终收敛相位分布；判断第一位置的最终收敛相位分布是否在物平面，若否，则用菲涅尔传播因子对第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。本发明实施例中，在进行迭代计算前，先对采集的物体图像进行了图像精确配准，在进行迭代计算时，先用第一种迭代算法进行计算收敛相位分布，然后用该收敛相位分布为初始值继续用第二种迭代算法进行计算得到最终的收敛相位分布，可以增强相位恢复的精确度、效率和稳定性。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：闪存盘、只读存储器（Read-Only Memory，ROM）、随机存取器（Random Access Memory，RAM）、磁盘或光盘等。

以上对本发明实施例所提供的基于混合迭代的同轴相衬成像相位恢复方法及系统进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于混合迭代的同轴相衬成像相位恢复方法，其特征在于，包括：

采集第一位置的物体图像以及所述第一位置的背景图像；

采集第二位置的物体图像以及所述第二位置的背景图像；

根据所述第一位置的物体图像和所述第一位置的背景图像，以及所述第二位置的物体图像和所述第二位置的背景图像，将采集的所述第一位置的物体图像和所述第二位置的物体图像进行归一化处理，选择归一化处理后的所述第一位置的物体图像为参照物体图像，选择归一化处理后的所述第二位置的物体图像为对照物体图像，将所述参照物体图像和所述对照物体图像进行图像精确配准，获得第三物体图像；

根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布；

根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布。

2.如权利要求1所述的方法，其特征在于，所述采集第一位置的物体图像以及所述第一位置的背景图像之前，还包括：

设置射线源的优化工作参数以及探测器的工作模式。

3.如权利要求1所述的方法，其特征在于，所述根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布之前，还包括：

计算焦斑几何模糊光学传输函数；

计算探测器响应的光学传输函数；

将所述参照物体图像和所述第三物体图像进行去卷积处理；

设置所述第一位置的初始相位分布。

4.如权利要求1所述的方法，其特征在于，所述第一迭代算法为傅里叶变换算法，所述第二迭代算法为盖师贝格-撒克斯通算法或泊松分布算法。

5.如权利要求1-4任一项所述的方法，其特征在于，所述根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布之后，还包括：

判断所述第一位置的最终收敛相位分布是否在物平面，若否，则用菲涅尔传播因子对所述第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。

6.一种基于混合迭代的同轴相衬成像相位恢复系统，其特征在于，包括：

第一单元，用于采集第一位置的物体图像以及所述第一位置的背景图像；

所述第一单元，还用于采集第二位置的物体图像以及所述第二位置的背景图像；

第二单元，用于根据所述第一位置的物体图像和所述第一位置的背景图像，以及所述第二位置的物体图像和所述第二位置的背景图像，将采集的所述第一位置的物体图像和所述第二位置的物体图像进行归一化处理，选择归一化处理后的所述第一位置的物体图像为参照物体图像，选择归一化处理后的所述第二位置的物体图像为对照物体图像，将所述参照物体图像和所述对照物体图像进行图像精确配准，获得第三物体图像；

第三单元，用于根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布；

第四单元，用于根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布。

7.如权利要求6所述的系统，其特征在于，所述系统还包括：

第五单元，用于在所述第一单元采集第一位置的物体图像以及所述第一位置的背景图像之前，设置射线源的优化工作参数以及探测器的工作模式。

8.如权利要求6所述的系统，其特征在于，所述系统还包括：

第六单元，用于在所述第三单元根据第一迭代算法、所述第三物体图像和所述参照物体图像，计算所述第一位置的初始收敛相位分布之前，计算焦斑几何模糊光学传输函数；

第七单元，用于计算探测器响应的光学传输函数；

第八单元，用于将所述参照物体图像和所述第三物体图像进行去卷积处理；

第九单元，用于设置所述第一位置的初始相位分布。

9.如权利要求6所述的系统，其特征在于，所述第一迭代算法为傅里叶变换算法，所述第二迭代算法为盖师贝格-撒克斯通算法或泊松分布算法。

10.如权利要求6-9任一项所述的系统，其特征在于，所述系统还包括：

第十单元，用于在所述第四单元根据第二迭代算法、所述第三物体图像、所述参照物体图像以及所述第一位置的初始收敛相位分布，计算所述第一位置的最终收敛相位分布之后，判断所述第一位置的最终收敛相位分布是否在物平面；

第十一单元，用于在所述第十单元的判断结果为否时，用菲涅尔传播因子对所述第一位置的最终收敛相位分布进行去卷积处理，获得物体在物平面的最终收敛相位分布。

Images