CN103460251A

CN103460251A - 利用对具有光栅布置的相衬成像的约束优化的图像积分的方法和系统

Info

Publication number: CN103460251A
Application number: CN2011800600912A
Authority: CN
Inventors: M·施坦帕诺尼; T·蒂林
Original assignee: Scherrer Paul Institut
Current assignee: Scherrer Paul Institut
Priority date: 2010-12-13
Filing date: 2011-12-09
Publication date: 2013-12-18
Also published as: WO2012080125A1; US20130279659A1; JP5818909B2; JP2014506352A; EP2652708A1; CA2821145A1; EP2652708B1; AU2011344365A1

Abstract

本发明利用一种光栅布置从物体中产生高质量、无伪像的相衬图像。新方法抑制需要直接图像积分，并显著提高了相衬图像的质量。相比已有技术，本方法不需要附加的校准工作以及增加的曝光时间。另一方面，其传送了射线照相实验中有关相位投影的出色的、直接可解释的信息。由于在使用上其通用可应用性和其简易性，所建议的发明可期望成为利用光栅布置的多种2D成像应用的标准方法，特别是在医学扫描仪(例如乳房X射线照相术)、工业生产线的检查、无损检测以及国土安全中。

Description

利用对具有光栅布置的相衬成像的约束优化的图像积分的方法和系统

本发明涉及一种用于利用约束优化从差分图像(differential image)中恢复积分图像(integrated image)的方法和系统。

众所周知，与传统的可见光学器件不同，X射线光学器件的折射率非常接近1并小于1。在第一种近似中，对于介质中小且微不足道的各向异性，表征组织的光学性质的折射率可表达为-包括X射线吸收-具有其复杂形式：n＝l-δ-iβ，其中δ为折射率实部的减量，表征相移性质，而虚部β描述样本的吸收性质。在基于传统吸收的射线照相术中，X射线的相移信息通常不直接用于图像重构。然而，在大于10keV的光子能量中并且对于轻质材料(由低原子序数元素组成)，相移项比衰减项扮演了更重要的角色，因为δ的量级典型地比β大三级。因而，相比于传统的基于吸收的成像，相衬模式可生成明显较大的图像对比度。

此外，远离吸收边沿，δ与X射线能量的平方成反比，而β降低为能量的四次方。这种机制的明显结果是，获取的相位信号可具有比吸收更低的剂量沉积，必须考虑当辐射损坏时的一个非常重要的问题，例如在生物样本或在生态系统中。

已经开发了某些方法以便记录该相位信号。它们可分为干涉方法(利用晶体)、相位传播方法、基于分析器晶体或X射线光栅的技术。所描述的发明特别以最后一种技术为背景。

基于X射线成像装置的光栅实质上检测X射线在物体中的偏转。这样的偏转可由导致差分相衬(DPC)的物体中的相移梯度的折射，或由导致所谓暗场图像(DFI)对比度的样本中的非均匀性散射而引起。该DPC图像信号可用于通过图像处理程序获取相衬(PC)图像。

如图1所示，具有两个光栅(G1和G2)或三个光栅(G0、G1和G2)的装置可应用于记录X射线的偏转。在两光栅装置的情况下，源需要满足关于其空间相干性的某种需求，而在三光栅装置中，不需要空间相干性。当光源的尺寸大于p2*l/d时，需要光栅G0，其中p2为G2的周期，l为源与G1之间的距离，而d为G1与G2之间的距离。因此，三光栅装置适用于非相干X射线源，尤其是X射线管。

为使传统衰减对比度(AC)与DPC和DFI对比度分离，应用了相位步进方法。光栅之一横向移位至入射光束，同时获取多个图像。检测器平面中每个像素的强度信号振荡，作为位移的函数。振荡的均值表示衰减对比度(AC)。振荡相位可直接关联于波前相位分布的一阶导数，因而关联于DPC信号。振荡幅度取决于X射线在物体中的散射，并且因而产生DFI信号。

对于所述(两个或三个)光栅，已经建议并应用了某些方法。光栅G0(若需要)是距离光源最近的一个光栅。其通常由具有周期p0的吸收线的透射光栅组成。其可由仅从具有相同周期的线中发出射线的源所代替。光栅G1进一步设置于源的下游。其由具有周期pl的线组成。光栅G2处于装置的最下游。其通常由具有周期p2的吸收线的透射光栅组成。其能够可替换地由具有相同周期的类似光栅敏感度的检测器系统所代替。

该装置的两种状态可区别为：所谓“近场状态”和“塔尔波特状态”。在“近场状态”中，光栅周期p、光栅距离d以及X射线波长λ被选择，以使衍射效应可忽略。在这种情况下，所有光栅需要包括吸收线。在“塔尔波特状态”中，光栅结构的衍射显著。两个状态之间的明显差别并未简单给出，因为该确切标准取决于光栅结构的占空因数，以及光栅是否正在吸收或相移。例如，对于具有吸收线和占空因数0.5的光栅，“近场状态”的条件为d≥p²/2λ。此处，G1应由正在吸收，或者优选为正在相移的光栅线组成。可有若干数量的相移，优选为π/2或其倍数。该光栅周期必须与光栅间的相对距离相匹配。在“塔尔波特状态”中的装置的情况下，需要考虑塔尔波特效应，以获取好的对比度。用于光栅周期和距离的公式为现有技术所知。

该样本主要置于G0和G1之间(或在两光栅装置的情况下位于G1的上游)，然而将其置于G1和G2之间是有利的。

本发明与上述所有情况相关，即两光栅和三光栅的情况、“近场状态”和“塔尔波特状态”的情况，以及对于置于G1上游或下游的样本。

此外，所提出的本发明还结合基于扫描的系统来工作，用于利用平面光栅平行和准平行几何结构，或用于利用圆柱或球形曲面光栅的紧密的扇形束或锥形束的几何结构。

具有上述系统的一种成像实验对X射线折射敏感。折射角α(相对于光轴)与利用系统所测量的差分相移之间的关系由下式给出：

其中d为光栅间的距离，并且p₂为G2的周期。X轴对应于相位步进方向。此外，α(x)与样本后波形的相位分布

之间的关系由[12]给出：

α (x) = \frac{λ}{2 π} \frac{&PartialD; φ (x)}{&PartialD; x} - - - (2)

将公式(2)代入公式(1)，产生

因此，

与φ(x)的一阶导数成比例。根据该公式，φ(x)的重构需要x方向上的DPC测量的积分，

总体上，噪声信号的积分将噪声误差和噪声方差累加。结果，具有增加幅度的水平条纹伪像在积分的方向产生。图2示出了噪声DPC图像的积分，其由修正的SheppLogan幻象所生成，其中在图2c中，条纹伪像清晰可见。

因而本发明的一个目标在于，提供一种系统和方法，以提高积分的相衬图像质量。

根据本发明，该目标通过独立权利要求1和19的特征来实现。本发明的优选实施例分别在从属权利要求2至18和20至36中给出。

用于利用约束优化从差分图像中恢复积分图像的创造性的系统和创造性的方法包括：

a)用于标识图像的装置，其将成本函数

最小化，服从约束

其写成：

最小化

服从

b)其中f为图像矢量，T为变换算子矩阵，D_x为差分算子矩阵，W为对角矩阵，包含可逆噪声标准偏差，为测量的图像矢量而ε为噪声功率的边界。

表示矢量的p-范数，p为正数，优选为间隔[0，2]中的整数。

当差分图像获取自用于X射线，尤其是硬性X射线的装置时实现本发明优选的实施例，用于从样本中获取定量的X射线图像，包括：

a.X射线源；

b.三个或至少两个光栅，分别命名为G0、G1与G2，或G1与G2；

c.位置灵敏检测器，具有含有多个独立像素的空间调制的检测灵敏度；

d.用于记录检测器图像的装置；

e.用于评估图像序列中每个像素强度，以便将用于每个独立像素的物体特征标识为吸收控制(dominate)的像素和/或差分相衬控制的像素和/或X射线散射控制的像素的装置；

f.其中该图像的序列通过连续或逐步地从0到π或2π来旋转样本或相对于样本旋转该装置与源被收集。

优选地，该系统和方法可操作于所谓“近场状态”，或“塔尔波特状态”中。

用于光栅G1的优选的实施例将G1提供为线光栅，其为吸收光栅或相位光栅，其中该相位光栅为低吸收光栅但生成相当多的X射线相移，后者优选为∏/2或其倍数。

因此，光栅G2可实现为线光栅，其具有高X射线吸收对比度，而其周期与G1的自身图像的周期相同；G2位于检测器的前方，其线与G1的线平行。

该系统和方法考虑了一定的操作自由度，其中操作可在平行束、准平行束、扇形束或锥形束模式中选择，而G0、G1和G2分别具有相应的平面、圆柱形或球形。因而，该操作可在具有二维光栅的全场模式中或具有一维光栅的扫描模式中选择。

对于该系统和方法的装置，两种类型的操作可定义如下：a)对于近场状态操作，光栅之间的距离在状态内自由选择，以及b)对于塔尔波特状态，根据下式选择：

D_{n, sph} = \frac{L \cdot D_{n}}{L - D_{n}} = \frac{L \cdot n \cdot p_{1}^{2} / 2 η^{2} λ}{L - n \cdot p_{1}^{2} / 2 η^{2} λ}

其中n＝1，3，5，......，并且

其中l＝1，2，3，......，D_n为当应用平行X射线束时的奇分数塔尔波特距离，而D_n，sph为当应用扇形或锥形X射线束时的距离，L为源与G1之间的距离。

为了受益于该差分相衬方法的全部优点，该系统和方法通过以下方式执行，即相位步进通过光栅G0、G1和G2之一相对于其它光栅的机械移动来实现。

关于光栅的结构，光栅结构可有利地通过平面技术来制造。

为了将图像中的差分相位信息进行分离，如来自CT、MRI或超声设备的医学图像，可根据欧洲专利申请EP10167569.2来获取差分相位信息，其通过引用合并于此。

在本发明进一步优选的实施例中，G1和G2之间的相位关系可精确地对应于使强度曲线通过一阶泰勒级数可展开的值，并且该差分相位信息可优选地根据国际专利申请PCT/EP2010/051291(WO2010/089319)来获取，其通过引用合并于此。

在数学成本函数的求解期间，算子D_x可设计为光栅的相位步进方向上的任意阶差分算子。因而，变换算子T可选择为垂直于光栅步进方向上的任意阶差分算子。进一步地，加权算子W可选择为对角加权矩阵，包含每个像素中DPC图像的逆向标准偏差1/σ_DPC。

优选地，该约束优化问题可通过将其重新计算为二阶锥形规划(SOCP)来求解。可替换地，约束优化问题可通过将其根据下式计算为可能无约束的形式来求解：

最小化

可替换地，该无约束优化问题可利用梯度下降或(非线性)共轭梯度算法来求解。

因此，该新发明解决条纹伪像对噪声DPC图像的直接积分的问题。基本思想是通过求解约束优化问题来抑制图像垂直方向上的变化。该相位图像通过最小化成本函数来检索，同时保持与测量数据的一致性。应用于噪声DPC测量的情况，其促使该积分生成图像中较低的变化，并且因此提高了图像质量。

本发明的优选实例在下文中更详细地描述。

下文中的图像，如医学图像或来自结构分析或材料测试等诸如此类的图像，将表示为矢量。图像矢量f(i)获取自图像I(x，y)像素值的按列提取：

其中n_x×n_y为图像尺寸。图像矢量的尺寸为l×n，n＝n_x·n_y。图像变换由算子(矩阵)表示，其可应用于图像矢量。算子矩阵的尺寸为m×n，其中m为变换矢量的尺寸。在多数情况下，m＝n成立。

新方法基于标准线性回归模型，其中DPC图像的测量由下式给出：

D_x为测量算子，对

的差分测量进行建模，而w为随机矢量，对DPC图像的噪声进行建模。D_x可利用

在x方向上的有限差分变换来实现。

关于噪声模型，DPC图像像素的噪声方差由下式给出：

σ_{DPC}^{2} &Proportional; \frac{1}{N V^{2}} - - - (7)

其中N为光子数，并且V为该像素上的平均条纹可见度。N和V可利用该测量的AC和DFI图像来计算。

对于来自DPC测量的相位检索，约束优化问题可定义为：

最小化

服从：

在该优化问题中，f为图像矢量，T为变换算子，为测量的DPC图像矢量，以及ε为用于噪声功率的边界。W为考虑噪声模型的对角加权矩阵，

表示矢量的p-范数，其被定义为：

(8)中表达的问题对任意f寻找矢量Tf最小的p-范数，其满足数据一致性约束

在(8)的最小项中，图像可变换为任一线性变换域，由矩阵T(例如，傅里叶变换、小波、有限差分等)表示。其使数据约束的利用非常灵活。在DPC测量的情况下，如下有限差分变换是优选的：

T＝D_y (10)

因为积分的图像由垂直于相位步进的方向上高强度的变化(水平条纹)而失真。

有多种可能来求解问题(8)。在凸面情况下(p≥1)，其可被重新计算并求解为二阶锥形规划。一个替换是利用拉格朗日乘子将(8)重新计算为无约束的形式，

最小化

或一般地，利用任意数量的数据约束，

最小化

该问题公式还称为正则化，其主要应用于不适定问题的反转。拉格朗日乘子(或正则化参数)λ_i控制正则化项相比于数据一致性项

的权重。

问题(12)具有强大性质以允许任意数量的正则化项，其在更多关于物体先验知识可用的情况下尤其有用。

正则化项中范数参数p的选择取决于所应用的变换算子T。典型选择为p＝2，因为在这种情况中，问题(11)/(12)是线性的，并存在显式解。对于T＝D_y的情况，-范数可导致模糊，并且因此降低了图像分辨率。另一方面，众所周知

-范数(p＝1)的最小化以保留图像中的边沿，并且因此可期望产生比

-范数更好的图像分辨率[1]。

-范数的最小化致使其中不存在显式解的非线性优化问题。

此处，迭代算法用于求解该非线性优化问题。其基于非线性共轭梯度(NLCG)方法，特征在于对大规模线性系统反转的快速收敛[2]。

参考文献

[1]L.Rudin，S.Osher，and E.Fatemi，″Nonlinear total variationbased noise removal algorithms，″Phys.D Nonlinear Phenom.60，259？268(1992)

[2]J.Nocedal and S.Wright，Numerical optimization(Springerverlag，1999)

Claims

1.一种用于利用约束优化从差分图像中恢复积分图像的方法，包括如下步骤：

标识一图像，其将成本函数

最小化，服从约束

其写成：

最小化

服从

其中f为图像矢量，T为变换算子矩阵，D_x为差分算子矩阵，W为对角矩阵，包含可逆噪声标准偏差，

为测量的图像矢量而ε为用于噪声功率的边界。

表示矢量的p-范数，其中p为正数，优选为间隔[0，2]中的整数。

2.根据权利要求1的方法，其中该差分数据获取自用于X射线，尤其是硬性X射线的装置，用于从样本中获取定量的X射线图像，包括：

a.X射线源；

b.三个或至少两个光栅，命名为G0、G1与G2，或G1与G2；

c.具有多个独立像素的位置灵敏检测器，具有空间调制的检测灵敏度；

d.用于记录检测器图像的装置；

e.用于评估图像序列中每个像素强度，以便将每个独立像素的物体特征标识为吸收控制的像素和/或差分相衬控制的像素和/或X射线散射控制的像素的装置；

f.其中该图像序列通过连续或逐步地从0到π或2π来旋转样本或相对于样本旋转该装置与源来收集。

3.根据权利要求1或2的方法，操作于所谓“近场状态”或“塔尔波特状态”中。

4.根据前述权利要求1至3任一的方法，其中G1为线光栅(G1)，吸收光栅或相位光栅，其为低吸收光栅但生成相当多的X射线相移，后者优选为π/2或其倍数。

5.根据前述权利要求1至4任一的方法，其中G2为线光栅，其具有高X射线吸收对比度，而其周期与G1的自身图像的周期相同；G2位于检测器的前方，其线与G1的线平行。

6.根据前述权利要求1至5任一的方法，其中操作可在平行束、准平行束、扇形束或锥形束模式中选择，而G0、G1和G2分别具有相应的平面、圆柱形和球形。

7.根据前述权利要求1至6任一的方法，其中该操作可在具有二维光栅的全场模式或具有一维光栅的扫描模式中选择。

8.根据前述权利要求1至7任一的方法，其中对于近场状态操作，光栅之间的距离在状态中自由选择，且对于塔尔波特状态，根据下式选择：

D_{n, sph} = \frac{L \cdot D_{n}}{L - D_{n}} = \frac{L \cdot n \cdot p_{1}^{2} / 2 η^{2} λ}{L - n \cdot p_{1}^{2} / 2 η^{2} λ}

其中n＝1，3，5，......，并且

9.根据前述权利要求1至8任一的方法，其中相位步进通过光栅G0、G1和G2之一相对于其它光栅的机械移动来执行。

10.根据前述权利要求1至9任一的方法，其中光栅结构通过平面技术来制造。

11.根据前述权利要求1至10任一的方法，其中差分相位信息根据欧洲专利申请EP10167569.2获得。

12.根据前述权利要求1至11任一的方法，其中G1和G2之间的相位关系精确地响应于使强度曲线通过一阶泰勒级数可展开的值，并且该差分相位信息优选地根据国际专利申请PCT/EP2010/051291(WO2010/089319)来获得。

13.根据前述权利要求1至12任一的方法，其中算子D_x为光栅的相位步进方向上任意阶的差分算子。

14.根据前述权利要求1至13任一的方法，其中变换算子T为垂直于光栅相位步进方向上的任意阶的差分算子。

15.根据前述权利要求1至14任一的方法，其中加权算子W为对角加权矩阵，包含每个像素中DPC图像的逆向标准偏差1/σ_DPC。

16.根据前述权利要求1至15任一的方法，其中该约束优化问题通过将其重新计算为二阶锥形规划(SOCP)来求解。

17.根据前述权利要求1至15任一的方法，其中该约束优化问题通过将其根据下式计算为可能无约束的形式来求解：

最小化

18.根据前述权利要求1至17任一的方法，其中该无约束优化问题利用梯度下降或(非线性)共轭梯度算法或其它来求解。

19.一种用于利用约束优化从差分图像中恢复积分图像的系统，包括：

a)用于标识图像的装置，其将成本函数最小化，服从约束

其写成：

最小化

服从

b)其中f为图像矢量，T为变形算子矩阵，D_x为差分算子矩阵，W为对角矩阵，包含可逆噪声标准偏差，为测量的图像矢量而ε为用于噪声功率的边界。

20.根据权利要求19的系统，其中该差分图像获取自用于X射线，尤其是硬性X射线的装置，用于从样本中获取定量的X射线图像，包括：

a.X射线源；

b.三个或至少两个光栅，命名为G0、G1与G2，或G1与G2；

d.用于记录检测器图像的装置；

f.其中该图像序列通过连续或逐步地从0到π或2π来旋转样本或相对于样本旋转该装置与源收集。

21.根据权利要求19或20的系统，操作于所谓“近场状态”或“塔尔波特状态”中。

22.根据前述权利要求19至21任一的系统，其中G1为线光栅(G1)，为吸收光栅或为相位光栅，其中相位光栅为低吸收光栅但生成相当多的X射线相移，后者优选为π/2或其倍数。

23.根据前述权利要求19至22任一的系统，其中G2为线光栅，其具有高X射线吸收对比度，而其周期与G1的自身图像的周期相同；G2位于检测器的前方，其线与G1的线平行。

24.根据前述权利要求19至23任一的系统，其中操作可在平行束、准平行束、扇形束或锥形束模式中选择，而G0、G1和G2分别具有相应的平面、圆柱形或球形。

25.根据前述权利要求19至24任一的系统，其中该操作可在具有二维光栅的全场模式或具有一维光栅的扫描模式中选择。

26.根据前述权利要求19至25任一的系统，其中对于近场状态操作，光栅之间的距离在状态中自由选择，且对于塔尔波特状态，根据下式选择：

D_{n, sph} = \frac{L \cdot D_{n}}{L - D_{n}} = \frac{L \cdot n \cdot p_{1}^{2} / 2 η^{2} λ}{L - n \cdot p_{1}^{2} / 2 η^{2} λ}

其中n＝1，3，5，......，并且

27.根据前述权利要求19至26任一的系统，其中相位步进通过光栅G0、G1和G2之一相对于其它光栅的机械移动来执行。

28.根据前述权利要求19至27任一的系统，其中光栅结构通过平面技术来制造。

29.根据前述权利要求19至28任一的系统，其中差分相位信息根据欧洲专利申请EP10167569.2获得。

30.根据前述权利要求19至29任一的系统，其中G1和G2之间的相位关系精确地响应于使强度曲线通过一阶泰勒级数可展开的值，并且该差分相位信息优选地根据国际专利申请PCT/EP2010/051291(WO2010/089319)来获得。

31.根据前述权利要求19至30任一的系统，其中算子D_x为光栅的相位步进方向上任意阶差分算子。

32.根据前述权利要求19至31任一的系统，其中变换算子T为垂直于光栅相位步进方向上的任意阶差分算子。

33.根据前述权利要求19至32任一的系统，其中加权算子W为对角加权矩阵，包含每个像素中DPC图像的逆向标准偏差1/σ_DPC。

34.根据前述权利要求19至33任一的系统，其中该约束优化问题通过将其重新计算为二阶锥形规划(SOCP)来求解。

35.根据前述权利要求19至33任一的系统，其中该约束优化问题通过将其根据下式计算为可能无约束的形式来求解：

最小化

36.根据前述权利要求19至35任一的系统，其中该无约束优化问题利用梯度下降或(非线性)共轭梯度算法来求解。