CN104323790B - 同轴相衬成像方法及系统和相衬ct方法及系统 - Google Patents

Abstract

本发明提供了一种同轴相衬成像方法，通过微焦斑X射线源产生具有空间相干性的射线，所述射线穿透固定位置的旋转台上的成像物体；通过探测器接收所述射线透过所述成像物体后的图像，具体是分别对所述微焦斑X射线源工作在两个电压下产生的不同射线透过所述成像物体后的图像进行曝光成像，得到双能谱图像；分别测试在所述两个电压下所述成像物体的入射面处的光强分布和谱密度分布；根据所述双能谱图像、光强分布和谱密度分布进行相位恢复。此方法不需要进行高难度的多幅图像精确配准并且成本低。此外还提供了一种同轴相衬成像系统。同时基于此同轴相衬成像方法提供了一种相衬CT方法，可以恢复成像物体的多种信息，此外还提供了一种相衬CT系统。

Description

同轴相衬成像方法及系统和相衬CT方法及系统

技术领域

本发明涉及相衬成像领域，特别是涉及一种同轴相衬成像方法及系统和相衬CT方法及系统。

背景技术

传统的X射线成像基于吸收机理，在医学、环境、材料、工业、安检等方面具有重要的应用，然而，对于一些主要由C、H、O等轻元素组成的物质而言，如生物组织、纤维、聚合物材料，由于几乎没有吸收或只有很少的吸收，使得图像对比度很差。

X射线相位衬度成像(相衬成像)是近年来发展的另一种X射线成像技术，利用X射线通过物体后发生的相移变化成像。对轻元素而言，典型的X射线相移灵敏度比吸收灵敏度大一千倍以上，从而可明显改善弱吸收物体的图像对比度而成为了国内外研究的一个热点。

目前为止，主要有五种实现X射线相衬成像的技术，干涉相衬成像、衍射增强相衬成像、光栅剪切相衬成像、同轴相衬成像与编码孔径相衬成像。同轴相衬成像因为最容易实现而被特别关注，系统无需附加任何精密的光学元器件，可以基于实验室普通的微焦斑X射线源，利用X射线通过样品后在自由空间传播一定距离即能将相位信息转化为强度调制的菲涅耳衍射原理成像。

X射线相位恢复是实现同轴相衬成像的难点与关键。为了实现对于复合组分物体进行精确相位恢复，如肿瘤血管，它由肿瘤组织、肌肉组织、血管外壁、血管内壁和血液构成的，需要在不同位置进行多次曝光，使得在相位恢复前需要对多幅图像进行精确配准，且配准需要到亚像素水平，对基于微焦斑源的球面波入射，不同位置的图像分辨率、视场也不同，更增加了配准的难度。也能采用单个位置成像、多个波长成像的方式，但需要用到准单色的同步辐射源或者昂贵的光子计数探测器，难于推广应用。

通过相衬成像获得的图像是二维图像，国内外研究人员将相衬技术与CT(Computed Tomography)即电子计算机断层扫描理论相结合发展了基于相位信息的CT技术，通过多视角扫描得到相衬成像再经CT重建可获得物体三维图像，即相衬CT。传统的相衬CT仅能获得物体的结构信息。

发明内容

基于此，为了克服上述技术问题，有必要提供一种不需要进行高难度的多幅图像精确配准并且低成本的同轴相衬成像方法及系统，同时提供一种能获得物体多种信息的相衬CT方法及系统。

一种同轴相衬成像方法，所述方法包括：

通过微焦斑X射线源产生具有空间相干性的射线，所述射线穿透固定位置的旋转台上的成像物体；

通过探测器接收透过所述成像物体后的射线进行成像，具体是分别对所述微焦斑X射线源工作在两个电压下产生的不同射线透过所述成像物体后的射线进行曝光成像，得到双能谱图像；

分别测试在所述两个电压下所述成像物体的入射面处的光强分布和谱密度分布；

根据所述双能谱图像、光强分布和谱密度分布进行相位恢复。

在其中一个实施例中，所述根据所述双能谱图像、光强分布和谱密度分布进行相位恢复的步骤包括：

根据所述双能谱图像计算所述成像物体的出射面处的光强分布；

根据所述入射面处的光强分布和所述出射面处的光强分布计算投影图像，通过阈值分割获得所述成像物体所在的数据区域；

获取参考波长，根据所述参考波长、入射面处的谱密度分布计算所述数据区域对应的所述成像物体的出射面处的谱密度分布；

获取有效传播距离，并计算入射面处的光程分布；

根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

在其中一个实施例中，所述根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布的步骤包括：

根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数；

根据所述系数解泊松方程获得所述参考波长的成像物体平面处的相位分布；

根据所述参考波长的成像物体平面处的相位分布计算所述参考波长的成像物体平面处的吸收分布；

根据所述参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

一种同轴相衬成像系统，所述系统包括：

发射模块，用于通过微焦斑X射线源产生具有空间相干性的射线，所述射线穿透固定位置的旋转台上的成像物体；

成像模块，用于通过探测器接收透过所述成像物体后的射线进行成像，具体是分别对所述微焦斑X射线源工作在两个电压下产生的不同射线透过所述成像物体后的射线进行曝光成像，得到双能谱图像；

测试模块，用于分别测试在所述两个电压下所述成像物体的入射面处的光强分布和谱密度分布；

相位恢复模块，用于根据所述双能谱图像、光强分布和谱密度分布进行相位恢复。

在其中一个实施例中，所述相位恢复模块包括：

光强分布计算单元，用于根据所述双能谱图像计算所述成像物体的出射面处的光强分布；

数据区域计算单元，用于根据所述入射面处的光强分布和所述出射面处的光强分布计算投影图像，通过阈值分割获得所述成像物体所在的数据区域；

谱密度分布计算单元，用于获取参考波长，根据所述参考波长、入射面处的谱密度分布计算所述数据区域对应的所述成像物体的出射面处的谱密度分布；

光程分布计算单元，用于获取有效传播距离，并计算入射面处的光程分布；

相位恢复单元，用于根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

在其中一个实施例中，所述相位恢复单元包括：

系数计算单元，用于根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数；

相位分布计算单元，用于根据所述系数解泊松方程获得所述参考波长的成像物体平面处的相位分布；

吸收分布计算单元，用于根据所述参考波长的成像物体平面处的相位分布计算所述参考波长的成像物体平面处的吸收分布；

其它波长恢复单元，用于根据所述参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

一种相衬CT方法，所述方法包括：

上述任一实施例所提供的同轴相衬成像方法；

旋转所述旋转台，通过所述同轴相衬成像方法得到不同投影角度下的关于波长的相位分布和吸收分布；

基于波长和能量的关系得到关于能量的相位分布与吸收分布；

根据所述关于能量的相位分布与吸收分布采用重建算法重建所述成像物体的信息。

在其中一个实施例中，所述根据所述关于能量的相位分布与吸收分布采用重建算法重建成像物体的信息的步骤包括：

根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布；

根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布；

计算电子散射截面，根据所述电子散射截面计算等效原子序数分布。

一种相衬CT系统，所述系统包括：

上述任一实施例所提供的同轴相衬成像系统；

多角度恢复模块，用于旋转所述旋转台，通过所述同轴相衬成像系统得到不同投影角度下的关于波长的相位分布和吸收分布；

关于能量恢复模块，用于基于波长和能量的关系得到关于能量的相位分布与吸收分布；

重建模块，用于根据所述关于能量的相位分布与吸收分布采用重建算法重建所述成像物体的信息。

在其中一个实施例中，所述重建模块包括：

系数分布计算单元，用于根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布；

电子密度分布计算单元，用于根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布；

等效原子序数分布计算单元，用于计算电子散射截面，根据所述电子散射截面计算等效原子序数分布。

上述同轴相衬成像方法及系统，通过在一个固定的位置成像，使用微焦斑X射线源在不同电压下的两次曝光成像得到双能谱图像，进行相位恢复得到相衬成像，克服了两个或多个位置成像图像配准的难题；同时系统只用到了普通的微焦斑X射线源和普通的探测器，无需大型的同步辐射源或昂贵的光子计数探测器，降低了成像成本。

上述相衬CT方法及系统，通过改变成像物体投影角度，获取双能谱图像得到不同角度下的相衬成像，并进行重建，不仅获得成像物体结构信息还能得到成像物体成分信息同时灵敏度更高。

附图说明

图1为一个实施例中同轴相衬成像方法和相衬CT方法的应用环境图；

图2为一个实施例中同轴相衬成像方法的流程图；

图3为一个实施例中根据双能谱图像、光强分布和谱密度分布进行相位恢复的流程图；

图4为一个实施例中根据有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布的流程图；

图5为一个实施例中同轴相衬成像系统的结构框图；

图6为一个实施例中恢复模块的结构框图；

图7为一个实施例中相位恢复单元的结构框图；

图8为一个实施例中相衬CT方法的流程图；

图9为一个实施例中根据关于能量的相位分布与吸收分布采用重建算法重建成像物体的信息的流程图；

图10为一个实施例中相衬CT系统的结构框图；

图11为一个实施例中重建模块的结构框图。

具体实施方式

本发明实施例中的同轴相衬成像方法和相衬CT方法可应用于如图1所示的应用环境中，此环境包括微焦斑X射线源、固定位置的旋转台和探测器。其中，微焦斑X射线源可以产生具有空间相干性的射线，旋转台可以360度旋转并放置成像物体，探测器接收射线透过成像物体后的图像。

在一个实施例中，如图2所示，提出了一种同轴相衬成像方法，该方法包括：

步骤S110，通过微焦斑X射线源产生具有空间相干性的射线，射线穿透固定位置的旋转台上的成像物体。

本实施例中，微焦斑X射线源可以工作在不同的电压下，产生不同能谱的具有空间相干性的射线，成像物体放在固定位置的旋转台上并且保持角度不变，射线穿透成像物体。

步骤S120，通过探测器接收射线透过成像物体后的图像，具体是分别对微焦斑X射线源工作在两个电压下产生的不同射线透过成像物体后的图像进行曝光成像，得到双能谱图像。

本实施例中，探测器采用数字化探测器，并设置数字化探测器在全分辨工作模式。设置微焦斑X射线源分别工作在高、低管电压kVp₁、kVp₂下，产生不同的射线透过成像物体，探测器分别接收射线透过成像物体后的图像进行曝光成像，得到双能谱图像I₁和I₂。

步骤S130，分别测试在两个电压下成像物体的入射面处的光强分布和谱密度分布。

本实施例中，成像物体入射面如图1所示，为成像物体靠近微焦斑X射线源的一面，分别对成像物体进行高、低管电压kVp₁、kVp₂下入射面处光强分布测试，得到成像物体入射面处的光强分布I_in1和I_in2。分别对成像物体进行高、低管电压kVp₁、kVp₂下入射面处谱密度分布测试，得到成像物体入射面处的谱密度分布和

步骤S140，根据双能谱图像、光强分布和谱密度分布进行相位恢复。

从探测器探测的强度图像即双能谱图像I₁和I₂反解相位分布为非线性问题，为了获得解析解，有不同的线性近似方法，也就产生了不同应用条件下不同的相位恢复方法。衬度传输函数法，适用于弱吸收与相位缓变物体；强度传播方程法，适合于近场成像；基于傅里叶变换的迭代算法相对解析算法，更精确、稳定，但计算效率低。本实施例中，可以根据双能谱图像I₁和I₂、光强分布I_in1和I_in2和谱密度分布和采用强度传播方程法进行相位恢复。

在近轴菲涅尔衍射与近场条件下，物平面与像平面的光强差与物平面相位导数的关系可以用强度传播方程表示。假设M₀为对应参考波长λ₀射线穿过物体后的吸收变化，可表示为

M₀＝∫μ(r_⊥,z,λ₀)dz (1)

其中，μ为线性吸收系数，r_⊥为横向坐标，z为光轴方向坐标。

在横向吸收变化M′₀(r_⊥)接近0的条件即近似纯相位物体，球面波入射物平面处相位分布满足以下泊松方程：

$-\frac{R^{\′}{\mathrm{\λ}}_0}{2\mathrm{\π}}{\▿}_{\⊥}^2{\mathrm{\φ}}_0=\frac{-I\left(R_2\right)+{\∫S}^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}}{\∫S^{\mathrm{out}}\left(E\right){\mathrm{\σ}}^2\mathrm{dE}}---\left(2\right)$

其中，R′为有效传播距离R′＝R₁R₂/(R₁+R₂)，其中R₁为如图1所示X射线源与成像物体的距离，R₂为成像物体与探测器的距离，φ₀为物平面处在参考波长时的相位分布，为拉普拉斯算子，I(R₂)为像平面处测得光强分布，S^out(E)为成像物体出射面处谱分布，ψⁱⁿ(E)为成像物体入射面处光程分布，E为光子能量，且

$S^{\mathrm{out}}\left(E\right)=S^{\mathrm{in}}\left(E\right)\exp (-{\mathrm{\σ}}^3{\stackrel{\‾}{M}}_0)---\left(3\right)$

Sⁱⁿ(E)为入射面处谱密度分布，为样品的平均吸收，由成像物体入射面处总射线强度与图像的总强度决定，σ＝E₀/E＝λ/λ₀，E₀＝hc/λ₀为对应于参考波长的光子能量，h为普朗克常量，c为真空中的光速。方程(2)中仅有一个未知数φ₀，故只需一个测试图像I(R₂)即可通过求解泊松方程获得相位分布φ₀。其它波长如果远离成像物体的吸收边界，则对应的相位分布为φ(λ)＝σ³φ₀。I(R₂)就对应于本实施例中的能谱图像I₁或I₂，可见对于成像物体是近似纯相位物体只需要一个能谱图像就能进行相位恢复。

在横向吸收变化M₀′(r_⊥)不可忽略时，即吸收与相位相互混合，满足下述方程：

$a-b{M}_0^{\′}-c{\▿}_{\⊥}^2{\mathrm{\φ}}_0+d{M}_0^{\′}{\▿}_{\⊥}^2{\mathrm{\φ}}_0=0---\left(4\right)$

其中， $a=I\left(R_2\right)-{\∫S}^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}---\left(5\right)$

$b=-{\∫S}^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}---\left(6\right)$

c＝R′∫S^out(E)σ²dE (7)

d＝R′∫S^out(E)σ⁵dE (8)

方程(4)中有两个未知数M₀与φ₀，故需要两个不同的多色谱(双能谱)入射进行两次光强测试求解。物平面处相位分布满足以下泊松方程：

$-\frac{R^{\′}{\mathrm{\λ}}_0}{2\mathrm{\π}}{\▿}_{\⊥}^2{\mathrm{\φ}}_0=\frac{a_1{b}_2-a_2{b}_1}{b_1{c}_2-b_2{c}_1+a_2(d_1-d_2{b}_1/b_2)}---\left(9\right)$

其中，a、b、c、d的定义如公式(5)(6)(7)(8)所示，下标1、2分别表示第一、第二次光强测试，如： $a_j=I_j\left(R_2\right)-\∫S_j^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}_j^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}$ ,j＝1,2对应于本实施例中的对微焦斑X射线源工作在两个电压下进行的两次测试，I₁(R₂)和I₂(R₂)分别对应于双能谱图像I₁和I₂，求解泊松方程(9)可获得相位分布φ₀，代回方程(4)可求得吸收分布M₀。其它波长如果远离成像物体的吸收边界，则对应的吸收分布为M(λ)＝σ³M₀，对应的相位分布为φ(λ)＝σ³φ₀。可见对于成像物体是吸收与相位相互混合时需要两个能谱图像才能进行相位恢复。进行相位恢复得到相衬成像，可以有效显示成像物体的内部结构得到成像物体的二维图像。

本实施例中，通过在一个固定的位置成像，使用微焦斑X射线源在不同电压下的两次曝光成像得到双能谱图像，进行相位恢复得到相衬成像，克服了两个或多个位置成像图像配准的难题；同时系统只用到了普通的微焦斑X射线源和普通的探测器，无需大型的同步辐射源或昂贵的光子计数探测器，降低了成像成本。

在一个实施例中，成像物体为吸收与相位混合的多组份物体，如图3所示，步骤S140包括：

步骤S141，根据双能谱图像计算成像物体的出射面处的光强分布。

本实施例中，成像物体出射面如图1所示，为成像物体靠近探测器的一面，根据双能谱图像I₁和I₂，基于光强与距离平方成反比的关系，由I₁、I₂，计算在紧贴成像物体的出射面处的光强分布I_out1、I_out2，公式如下：

$I_{\mathrm{outj}}={\left(\frac{R_1+R_2}{R_1}\right)}^2{I}_j$

其中j＝1，2，其中R₁为如图1所示X射线源与成像物体的距离，R₂为成像物体与探测器的距离。

步骤S142，根据入射面处的光强分布和出射面处的光强分布计算投影图像，通过阈值分割获得成像物体所在的数据区域。

本实施例中，根据入射面处的光强分布I_in1和I_in2以及出射面处的光强分布I_out1和I_out2代入ln(I_inj/I_outj)分别计算2个电压下的投影图像，其中j＝1，2。然后设定一个阈值，根据ln(I_inj/I_outj)的值把小于阈值的像素去掉，即排除图像中无吸收的边缘空白区域，得到成像物体所在的数据区域。

步骤S143，获取参考波长，根据参考波长、入射面处的谱密度分布计算数据区域对应的成像物体的出射面处的谱密度分布。

本实施例中，选择一个参考波长λ₀，根据λ₀、成像物体入射面处的谱密度分布和代入 $S_j^{\mathrm{out}}\left(E\right)=S_j^{\mathrm{in}}\left(E\right)\exp (-{\mathrm{\σ}}^3{\stackrel{\‾}{M}}_{0j})$ 计算数据区域对应的成像物体的出射面处的谱密度分布其中j＝1，2，σ＝E₀/E＝λ/λ₀，为成像物体的平均吸收，n为图像的像素索引，N为图像中成像物体所在的数据区域的总像素数。

步骤S144，获取有效传播距离，并计算入射面处的光程分布。

本实施例中，根据R₁即如图1所示X射线源与成像物体的距离，R₂即成像物体与探测器的距离，计算有效传播距离R′＝R₁R₂/(R₁+R₂)，计算入射面处光程分布其中r_⊥为横向坐标，j＝1，2。

步骤S145，根据双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

本实施例中，根据双能谱图像I₁和I₂、有效传播距离R′，出射面处的谱密度分布入射面处的光程分布其中j＝1，2先解泊松方程得到对应参考波长λ₀的物平面处相位分布和吸收公布，再根据其它波长如果远离成像物体的吸收边界时满足φ(λ)＝σ³φ₀、M(λ)＝σ³M₀计算其它波长的相位分布和吸收分布，其中σ＝λ/λ₀。

在一个实施例中，如图4所示，步骤S145包括：

步骤S145a，根据双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数。

本实施例中，将双能谱图像I₁和I₂、有效传播距离R′、出射面处的谱密度分布和入射面处的光程分布其中j＝1，2代入以下公式计算系数a_j、b_j、c_j、d_j其中j＝1，2代表2个不同电压下的2套系数。

$a_j=I_j-S_j^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}_j^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}$

$b_j=-\∫S_j^{\mathrm{out}}\left(E\right)[1-R^{\′}{\▿}_{\⊥}^2{\mathrm{\ψ}}_j^{\mathrm{in}}\left(E\right)-\frac{2R^{\′}}{R_1}]\mathrm{dE}$

$c_j=R^{\′}\∫S_j^{\mathrm{out}}\left(E\right){\mathrm{\σ}}^2\mathrm{dE}$

$d_j=R^{\′}\∫S_j^{\mathrm{out}}\left(E\right){\mathrm{\σ}}^5\mathrm{dE}$

步骤S145b，根据得到的系数解泊松方程获得参考波长的成像物体平面处的相位分布。

本实施例中，将步骤S145a计算出来的a₁、b₁、c₁、d₁和a₂、b₂、c₂、d₂代入如下方程中解此泊松方程得到对应参考波长λ₀的成像物体平面处的相位分布φ₀。

$-\frac{R^{\′}{\mathrm{\λ}}_0}{2\mathrm{\π}}{\▿}_{\⊥}^2{\mathrm{\φ}}_0=\frac{a_1{b}_2-a_2{b}_1}{b_1{c}_2-b_2{c}_1+a_2(d_1-d_2{b}_1/b_2)}$

其中为拉普拉斯算子，λ₀为参考波长，R′为有效传播距离。

步骤S145c，根据参考波长的成像物体平面处的相位分布计算参考波长的成像物体平面处的吸收分布。

本实施例中，将a₁、b₁、c₁、d₁和对应参考波长λ₀的成像物体平面处的相位分布φ₀代入如下方程中得到对应参考波长λ₀的成像物体平面处的吸收分布M₀。 $a_0-b_1{M}_0^{\′}-c_1{\▿}_{\⊥}^2{\mathrm{\φ}}_0+d_1{M}_0^{\′}{\▿}_{\⊥}^2{\mathrm{\φ}}_0=0$

步骤S145d，根据参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

本实施例中，其它波长如果远离成像物体的吸收边界，则对应的相位分布为φ(λ)＝σ³φ₀，对应的吸收分布为M(λ)＝σ³M₀，代入σ＝λ/λ₀，得到

$M\left(\mathrm{\λ}\right)=\frac{{\mathrm{\λ}}^3}{{\mathrm{\λ}}_0^3}{M}_0.$

本实施例中，针对吸收与相位混合的多组份物体根据双能谱图像采用强度传播方程法进行相位恢复，此方法通过解泊松方程可对成像物体进行相位恢复简单方便，效率高，也不需要对图像进行配准。

如图5所示，在一个实施例中，提供了一种同轴相衬成像系统，包括：

发射模块210，用于通过微焦斑X射线源产生具有空间相干性的射线，射线穿透固定位置的旋转台上的成像物体。

成像模块220，用于通过探测器接收透过成像物体后的射线进行成像，具体是分别对微焦斑X射线源工作在两个电压下产生的不同射线透过成像物体后的射线进行曝光成像，得到双能谱图像。

测试模块230，用于分别测试在两个电压下成像物体的入射面处的光强分布和谱密度分布。

相位恢复模块240，用于根据双能谱图像、光强分布和谱密度分布进行相位恢复。

在另一个实施例中，如图6所示，上述相位恢复模块240包括：

光强分布计算单元241，用于根据双能谱图像计算成像物体的出射面处的光强分布。

数据区域计算单元242，用于根据入射面处的光强分布和出射面处的光强分布计算投影图像，通过阈值分割获得成像物体所在的数据区域。

谱密度分布计算单元243，用于获取参考波长，根据参考波长、入射面处的谱密度分布计算数据区域对应的成像物体的出射面处的谱密度分布。

光程分布计算单元244，用于获取有效传播距离，并计算入射面处的光程分布。

相位恢复单元245，用于根据双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

在一个实施例中，如图7所示，上述相位恢复单元245包括：

系数计算单元245a，用于根据双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数。

相位分布计算单元245b，用于根据系数解泊松方程获得参考波长的成像物体平面处的相位分布。

吸收分布计算单元245c，用于根据参考波长的成像物体平面处的相位分布计算参考波长的成像物体平面处的吸收分布。

其它波长恢复单元245d，用于根据参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

在一个实施例中，如图8所示，提出了一种相衬CT方法，也应用于如图1所示的环境中，该方法除了包括上述各实施例所提供的同轴相衬成像方法外，还包括以下步骤：

步骤S320，旋转旋转台，通过同轴相衬成像方法得到不同投影角度下的关于波长的相位分布和吸收分布。

本实施例中，旋转旋转台的角度和同轴相衬成像的次数可根据需要进行选择，进行n次同轴相衬成像，得到不同投影角度下的关于波长的相位分布和吸收分布，即其中j＝0-n代表不同角度，从而得到不同角度下的成像物体的内部结构的二维图像。

步骤S330，基于波长和能量的关系得到关于能量的相位分布与吸收分布。

本实施例中，基于波长和能量的关系E＝hc/λ，其中E为光子能量，h为普朗克常量，c为真空中的光速，h为波长。代入和得到 ${\mathrm{\φ}}_j\left(E\right)=\frac{h^3{c}^3}{E^3{\mathrm{\λ}}_0^3}{\mathrm{\φ}}_j$ 和 $M_j\left(E\right)=\frac{h^3{c}^3}{E^3{\mathrm{\λ}}_0^3}{M}_j$ 其中j＝0-n。

步骤S340，根据关于能量的相位分布与吸收分布采用重建算法重建成像物体的信息。

本实施例中，X射线与物体作用决定于物体折射率，同时考虑光的吸收与相位变化时，物体折射率可表示为一复数n(E)＝1-δ(E)+iβ(E)，其中n(E)为能量相关的复折射率，β(E)为复折射率的虚部，与物体吸收相关，δ(E)为复折射率的实部，与物体引起相移相关。关于能量的相位分布φ(E)与吸收分布M(E)可表示为，

$\mathrm{\φ}\left(E\right)=-\mathrm{\λ}\left(E\right)r_e\∫{\mathrm{\ρ}}_e\mathrm{dl}=-\frac{2\mathrm{\π}}{\mathrm{\λ}\left(E\right)}\∫\mathrm{\δ}\left(E\right)\mathrm{dl}$

$M\left(E\right)=\∫\mathrm{\μ}\left(E\right)\mathrm{dl}=\frac{4\mathrm{\π}}{\mathrm{\λ}\left(E\right)}\∫\mathrm{\β}\left(E\right)\mathrm{dl}$

其中，μ(E)为线性衰减系数，r_e为经典电子半径，ρ_e为电子密度分布。通过重建可以得到电子密度分布ρ_e或δ(E)、线性衰减系数μ(E)或β(E)。并且衰减系数μ(E)可进一步分解为光电吸收与散射两部分，如下

$\mathrm{\μ}\left(E\right)=\mathrm{NK}\frac{Z^5}{e^3\left(E\right)}+{\mathrm{NZ\σ}}_{\mathrm{KN}}\left(E\right)=K\frac{{\mathrm{\ρ}}_e{Z}^4}{e^3\left(E\right)}+{\mathrm{\ρ}}_e{\mathrm{\σ}}_{\mathrm{KN}}\left(E\right)$

其中，N为原子密度，K为无量纲常数，Z为等效原子序数，e(E)为相对单个电子能量归一化的光子能量，σ_KN(E)为单个电子的散射截面，可由如下Klein-Nishina(克莱因-仁科)方程求得

${\mathrm{\σ}}_{\mathrm{KN}}\left(E\right)=2{\mathrm{\πr}}_e^2\left\{\frac{1+\mathrm{\η}}{{\mathrm{\η}}^2}\right[\frac{2(1+\mathrm{\η})}{1+2\mathrm{\η}}-\frac{1}{\mathrm{\η}}\ln (1+2\mathrm{\η})]+\frac{1}{2\mathrm{\η}}\ln (1+2\mathrm{\η})-\frac{1+3\mathrm{\η}}{{(1+2\mathrm{\η})}^2}\}$

其中，η＝E/m_ec²，m_ec²是静止电子能量511keV。可以看出通过关于能量的相位与吸收分布，采用CT重建算法，可以重建成像物体多种信息，包括：电子密度分布ρ_e或δ(E)、线性衰减系数μ(E)或β(E)、等效原子序数Z。重建可以采用经典的滤波反投影算法(filterback-projection,FBP)，也可采用迭代重建算法，但考虑到投影数越多，相位恢复时间越长，辐射剂量也更大的问题，本实施例中选用可少角度采样的迭代重建算法。

常规CT仅能重建成像物体线性衰减系数μ，且不考虑能谱，本实施例中采用双能谱相衬成像后再进行CT重建可得到成像物体的多种信息，重建得到的线性衰减系数为μ(E)，考虑了能谱分布，可以减小射束硬化效应。同时对弱吸收成像物体，重建得到的电子密度分布ρ_e比线性衰减系数μ重建灵敏度更高。重建得到的等效原子序数Z可以反应成像物体的成分信息，突破传统吸收CT仅能获得成像物体结构信息的限制。

在一个实施例中，如图9所示，步骤S340包括：

步骤S341，根据关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布。

本实施例中，关于能量的吸收分布M(E)可表示为M(E)＝∫μ(E)dl其中，μ(E)为线性衰减系数。采用代数迭代算法(Algebraic Reconstruction Technique，ART)，从μ(E)的初始估计值出发，通过对估计值的进行迭代修正使其逐渐朝着μ(E)的真实值逼近，直到满足所需要求，则结束迭代，得到线性衰减系数分布μ(E)。

步骤S342，根据关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布。

本实施例中，关于能量的相位分布φ(E)可表示为φ(E)＝-λ(E)r_e∫ρ_edl其中，λ₀为参考波长，E₀＝hc/λ₀为对应于参考波长的光子能量，h为普朗克常量，c为真空中的光速，r_e为经典电子半径，ρ_e为电子密度分布。采用代数迭代算法(AlgebraicReconstruction Technique，ART)，从ρ_e的初始估计值出发，通过对估计值的进行迭代修正使其逐渐朝着ρ_e的真实值逼近，直到满足所需要求，则结束迭代，得到线性衰减系数分布ρ_e。

步骤S343，计算电子散射截面，根据电子散射截面计算等效原子序数分布。

本实施例中，σ_KN(E)为单个电子的散射截面，可由如下Klein-Nishina(克莱因-仁科)方程求得

${\mathrm{\σ}}_{\mathrm{KN}}\left(E\right)=2{\mathrm{\πr}}_e^2\left\{\frac{1+\mathrm{\η}}{{\mathrm{\η}}^2}\right[\frac{2(1+\mathrm{\η})}{1+2\mathrm{\η}}-\frac{1}{\mathrm{\η}}\ln (1+2\mathrm{\η})]+\frac{1}{2\mathrm{\η}}\ln (1+2\mathrm{\η})-\frac{1+3\mathrm{\η}}{{(1+2\mathrm{\η})}^2}\}$

其中，η＝E/m_ec²，m_ec²是静止电子能量511keV。

然后根据计算得到等效原子序数分布Z，其中，N为原子密度，K为无量纲常数，e(E)为相对单个电子能量归一化的光子能量

本实施例中采用代数迭代算法重建成像物体的信息，通过少角度稀疏采样，降低针对所有投影图像需要进行相位恢复的总时间，也可降低辐射剂量。

如图10所示，在一个实施例中，提供了一种相衬CT系统，除了包括上述各实施例所提供的同轴相衬成像系统外，还包括：

多角度恢复模块420，用于旋转旋转台，通过同轴相衬成像系统得到不同投影角度下的关于波长的相位分布和吸收分布。

关于能量恢复模块430，用于基于波长和能量的关系得到关于能量的相位分布与吸收分布。

重建模块440，用于根据关于能量的相位分布与吸收分布采用重建算法重建成像物体的信息。

如图11所示，在一个实施例中，重建模块440包括：

系数分布计算单元441，用于根据关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布。

电子密度分布计算单元442，用于根据关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布

等效原子序数分布计算单元443，用于计算电子散射截面，根据电子散射截面计算等效原子序数分布。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (8)

1.一种同轴相衬成像方法，所述方法包括：

通过微焦斑X射线源产生具有空间相干性的射线，所述射线穿透固定位置的旋转台上的成像物体；

通过探测器接收所述射线透过所述成像物体后的图像，具体是分别对所述微焦斑X射线源工作在两个电压下产生的不同射线透过所述成像物体后的图像进行曝光成像，得到双能谱图像；

分别测试在所述两个电压下所述成像物体的入射面处的光强分布和谱密度分布；

根据所述双能谱图像、光强分布和谱密度分布进行相位恢复；

所述根据所述双能谱图像、光强分布和谱密度分布进行相位恢复的步骤包括：

根据所述双能谱图像计算所述成像物体的出射面处的光强分布；

根据所述入射面处的光强分布和所述出射面处的光强分布计算投影图像，通过阈值分割获得所述成像物体所在的数据区域；

获取参考波长，根据所述参考波长、入射面处的谱密度分布计算所述数据区域对应的所述成像物体的出射面处的谱密度分布；

获取有效传播距离，并计算入射面处的光程分布；

根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

2.根据权利要求1所述的方法，其特征在于，所述根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布的步骤包括：

根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数；

根据所述系数解泊松方程获得所述参考波长的成像物体平面处的相位分布；

根据所述参考波长的成像物体平面处的相位分布计算所述参考波长的成像物体平面处的吸收分布；

根据所述参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

3.一种同轴相衬成像系统，其特征在于，所述系统包括：

发射模块，用于通过微焦斑X射线源产生具有空间相干性的射线，所述射线穿透固定位置的旋转台上的成像物体；

成像模块，用于通过探测器接收所述射线透过所述成像物体后的图像，具体是分别对所述微焦斑X射线源工作在两个电压下产生的不同射线透过所述成像物体后的图像进行曝光成像，得到双能谱图像；

测试模块，用于分别测试在所述两个电压下所述成像物体的入射面处的光强分布和谱密度分布；

相位恢复模块，用于根据所述双能谱图像、光强分布和谱密度分布进行相位恢复；

所述相位恢复模块包括：

光强分布计算单元，用于根据所述双能谱图像计算所述成像物体的出射面处的光强分布；

数据区域计算单元，用于根据所述入射面处的光强分布和所述出射面处的光强分布计算投影图像，通过阈值分割获得所述成像物体所在的数据区域；

谱密度分布计算单元，用于获取参考波长，根据所述参考波长、入射面处的谱密度分布计算所述数据区域对应的所述成像物体的出射面处的谱密度分布；

光程分布计算单元，用于获取有效传播距离，并计算入射面处的光程分布；

相位恢复单元，用于根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算关于波长的相位分布和吸收分布。

4.根据权利要求3所述的系统，其特征在于，所述相位恢复单元包括：

系数计算单元，用于根据所述双能谱图像、有效传播距离、出射面处的谱密度分布、入射面处的光程分布计算解泊松方程必须的系数；

相位分布计算单元，用于根据所述系数解泊松方程获得所述参考波长的成像物体平面处的相位分布；

吸收分布计算单元，用于根据所述参考波长的成像物体平面处的相位分布计算所述参考波长的成像物体平面处的吸收分布；

其它波长恢复单元，用于根据所述参考波长的成像物体平面处的相位分布和吸收分布得到其它波长的相位分布和吸收分布。

5.一种相衬CT方法，所述方法包括：

权利要求1至2任一项所述的同轴相衬成像方法；

旋转所述旋转台，通过所述同轴相衬成像方法得到不同投影角度下的关于波长的相位分布和吸收分布；

基于波长和能量的关系得到关于能量的相位分布与吸收分布；

根据所述关于能量的相位分布与吸收分布采用重建算法重建所述成像物体的信息。

6.根据权利要求5所述的方法，其特征在于，所述根据所述关于能量的相位分布与吸收分布采用重建算法重建所述成像物体的信息的步骤包括：

根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布；

根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布；

计算电子散射截面，根据所述电子散射截面计算等效原子序数分布。

7.一种相衬CT系统，其特征在于，所述系统包括：

权利要求3至4任一项所述的同轴相衬成像系统；

多角度恢复模块，用于旋转所述旋转台，通过所述同轴相衬成像系统得到不同投影角度下的关于波长的相位分布和吸收分布；

关于能量恢复模块，用于基于波长和能量的关系得到关于能量的相位分布与吸收分布；

重建模块，用于根据所述关于能量的相位分布与吸收分布采用重建算法重建所述成像物体的信息。

8.根据权利要求7所述的系统，其特征在于，所述重建模块包括：

系数分布计算单元，用于根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算线性衰减系数分布；

电子密度分布计算单元，用于根据所述关于能量的相位分布与吸收分布采用迭代重建算法计算电子密度分布；

等效原子序数分布计算单元，用于计算电子散射截面，根据所述电子散射截面计算等效原子序数分布。