CN103871056B - 基于各向异性光流场与去偏移场的脑部mr图像配准方法 - Google Patents

Abstract

本发明基于各向异性光流场与去偏移场的脑部MR图像配准方法，提出了基于光流场模型的图像配准与偏移场恢复耦合模型，针对偏移场易导致传统配准模型配准精度不高的缺点，将去偏场与求光流场两者相结合，纳入统一变分框架中，使得两者相辅相成，针对Horn模型全局性正则项因缺少图像信息不能准确指导光流运动不足，引入图像结构信息来正则光流场，从而得到光滑准确的光流信息。本发明在配准同时可以恢复图像灰度不均匀场，从而降低了偏移场的影响；并且引入图像结构信息来降低配准结果模糊程度以及保留了图像结构信息，保证了边界结构信息的完整性，有利于进一步恢复真实图像。本发明大大提升了配准结果的精度,鲁棒性好。

Description

基于各向异性光流场与去偏移场的脑部MR图像配准方法

技术领域

本发明属于图像配准技术领域，尤其是涉及一种改进的脑部MR图像配准方法。

背景技术

图像配准在遥感图像处理、计算机视觉、医学图像分析等领域有着广泛的应用。在医学图像领域，它是不同医学图像信息融合的基础。配准旨在对于指定一幅图像寻求一种空间变换，使它与另一幅医学图像上的对应点达到空间上的一致，这种一致是指人体上的同一解剖点在两张配准图像上有相同的空间位置。图像配准方法分为刚性和非刚性，由于人体组织的变化是非刚体的、非线性的，因此非刚性配准方法具有更好的配准效果。

目前，非刚性配准方法主要分为两大类:基于像素的方法和基于特征的方法。基于像素的方法中包括以互信息为相似性测度的配准方法以及基于光流场模型的配准方法。基于特征方法主要依赖特征点选择，常用于刚性配准。由于基于像素的方法不需要对图像进行分割等预处理，从而提高了配准精度，受到越来越多的重视。

互信息方法因其具有较高精度和鲁棒性，是目前最广泛用于医学图像配准方法之一。该方法假设两幅共同解剖结构图获得最佳配准时，其对应像素灰度互信息值最大。由于使用互信息，使得该类方法可以配准多模态图像。然而，该方法计算复杂度较高，难以达到医学分析实时要求，且由于仅仅使用邻域熵信息，导致相似性度量不能满足全局凸性质，从而易陷入局部最优。

鉴于基于配准所求解的位移场与光流场模型所求解的速度场之间的相似性，Palos、Pierre等人将光流场模型引入到了图像配准中，他们提出的方法都基于Horn模型:

由于图像的配准过程可以认为是从待配准图像流动到基准图像的过程，即配准所要求解的位移场可以看作光流场所求解的速度场，因此，可以通过求解光流场进行图像配准：

其中，(I_x,I_y,I_t)是图像I沿x，y和t方向的导数，(u,v)^T为所求光流场。该方程求解存在孔径问题，针对这点，Horn等人根据同一运动物体引起的光流场是连续的和平滑的特点，最先在光流场基本约束上引入光滑性约束[3，7]，并将光流场的求解问题转化为求以下能量求极小值问题：

其中，α是光滑系数，其值越大，演化过程中对图像的平滑效果越明显。通过变分，利用梯度下降流方法，得到光流场u和ν迭代方程：

相对于Horn的全局化模型，Lucas等人提出了局部化窗口模型其思想是将每一象素速度矢量的求解放到其邻域中，由邻域内各象素的信息共同约束该点速度矢量。窗口化操作中，每一象素只受其周围邻域象素的影响，可通过限定小的窗口来减少对图像造成的模糊。

其中，W是一个符合统计特性的近邻象素比远距离象素对中心象素影响大的窗口函数。

该类模型假设同一运动物体引起光流场是连续的、平滑的，并且提出了相近像素点有相同的运动速度的光滑约束条件，较适用图像灰度一致性假设，但是针对脑部MR图像存在偏移场现象，即图像灰度呈现缓慢变化，使得目标图像中各组织之间区分度降低，往往导致配准误差较大。该光滑性约束虽能够得到较致密光流场，但不能保持图像的不连续性,导致配准图像在演化过程中出现模糊现象。医学图像各个组织之间以及角点区域往往不连续，该模型配准结果会丢失感兴趣细节，成为限制光流场方法应用于医学图像配准的主要问题。且对于脑部MR图像，由于设备成像机制的影响，图像含有灰度不均匀场（常称为偏移场）现象，导致灰度不均匀，而Horn模型假设图像灰度变化不大，从而易降低该模型配准精度。

发明内容

为解决上述问题，本发明提出了基于光流场模型的图像配准与偏移场恢复耦合模型，主要是由数据项能量项和正则项能量项两部分组成。针对偏移场易导致传统配准模型配准精度不高的缺点，本发明将去偏场与求光流场两者相结合，纳入统一变分框架中，使得两者相辅相成，针对Horn模型全局性正则项因缺少图像信息不能准确指导光流运动不足，本发明引入图像结构信息来正则光流场，从而得到光滑准确的光流信息。

为了达到上述目的，本发明提供如下技术方案：

一种基于各向异性光流场与去偏移场的脑部MR图像配准方法，包括如下步骤：

步骤1计算光流场(u,ν)以及偏移场B，根据求解的光流场(u,ν)作为配准的位移场，拟合得到配准后图像，

光流场u通过下式运算：

光流场ν通过下式运算：

偏移场通过下式运算：

配准后图像的运算公式为：

其中，I(x,y,t)为基准图像，(x,y)为像素坐标信息,t为时间信息，新的结构张量新的结构张量对应的特征值如下：

其中，0＜α＜1,ε＞0，为常数；

步骤2将作为新的待配准图像，基准图像不变，计算光流场(u,ν)以及偏移场B，相应得到更新后的配准图像

步骤3判断两次配准图像I₃和的差值小于预先设定的阀值时，得到经步骤2更新后的光流场(u,ν)为最终配准所求位移场，得到I₃为光流场模型配准结果，根据偏移场B从而恢复真实目标图像，当两次配准图像I₃和的差值小于预先设定的阀值时，将I₃作为下一步迭代的待配准图像，执行步骤2。

作为本发明的一种优选方案，所述光流场u和光流场ν的公式通过如下步骤得到：

步骤101将去偏移场与光流场求解两者过程相结合，纳入统一变分框架中，提出以下能量泛函：

其中，I(x,y,t)为基准图像，I(x+Δx,y+Δy,t+Δt)为待配准图像，B为待恢复的偏移场，

步骤102对图像采用对数变换得到下式：

步骤103对采用泰勒展开式：

步骤104令综合式（7）和式（9）式，将（8）可以改写如下形式：

步骤105采用梯度下降流方法，得到偏移场B表达式如下：

步骤106采用低通滤波对偏移场进行光滑，即：

作为本发明的一种优选方案，所述步骤2按照预先设定的迭代次数重复执行。

本发明将去偏移场与配准两者过程相结合，在配准同时可以恢复图像灰度不均匀场，从而降低了偏移场的影响；并且引入图像结构信息来降低配准结果模糊程度以及保留了图像结构信息，保证了边界结构信息的完整性，有利于进一步恢复真实图像。本发明大大提升了配准结果的精度,鲁棒性好。

附图说明

图1为本发明提供的脑部MR图像配准方法的步骤流程示意图；

图2为采用本发明方法和现有方法进行虚拟脑MR图像配准的结果比较图，其中虚拟脑MR图像含有低对比度区域以及细长拓扑结构组织；

图3为采用本发明方法和现有方法进行虚拟脑MR图像配准的结果比较图，其中虚拟脑MR图像含有较强的灰度不均匀现象；

图4为采用本发明方法和现有方法进行虚拟脑MR图像配准的结果比较图，其中虚拟脑MR图像不仅含有较强的灰度不均匀现象而且含有较强噪声；

图5为采用本发明方法和现有方法进行真实脑部图像配准的结果比较图。

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

要实现本发明提供的更为精确的图像配准，需要首先建立新的图像配准模型。基于变分方法光流场配准模型主要是由去偏变分框架和正则项能量项两部分组成。针对偏移场易导致传统配准模型配准精度不高的缺点，本发明将去偏场与求光流场两者相结合，纳入统一变分框架中，使得两者相辅相成。针对Horn模型全局性正则项因缺少图像信息不能准确指导光流运动不足，本发明将引入图像结构信息来正则光流场，使之得到光滑准确的光流信息。

因此，首先我们必须建立基于光流与去偏变分框架：

由于受到射频场不均匀性、设备本身、组织之间差异等影响，使得脑部MR图像同一组织像素灰度沿空间呈缓慢变化，即图像带有偏移场现象。设观察得到图像为I，真实图像为J，偏移场为B，噪声为N，则有，

I＝J×B+N （5）

本发明中暂不考虑噪声影响。为了恢复图像中的偏移场B，本发明借鉴多数学者采用的将偏移场B与真实图像J相分离方法，则对式（5）等式两边取对数，可得：

logI＝logB+logJ （6）

为方便起见，我们令则（6）可写成当图像含有较强偏移场时，图像目标和背景之间区别度将会降低。因此，直接使用原始光流场模型配准，会使得配准感兴趣区域分析带来较大误差。为了降低偏移场对配准影响，本发明将去偏移场与光流场求解两者过程相结合，纳入统一变分框架中，提出以下能量泛函：

其中，I(x,y,t)为基准图像，I(x+Δx,y+Δy,t+Δt)为待配准图像，B为待恢复的偏移场。为了方便求出偏移场B，对图像采用对数变换得到下式：

对采用泰勒展开式，不考虑高次阶导数影响，可得：

令综合式（7）和式（9）式，可将（8）可以改写如下形式：

极小化式（10），采用梯度下降流方法，得到偏移场B表达式如下：

考虑到图像偏移场具有光滑性特点，本发明采用低通滤波对偏移场进行光滑，即：

其中，*为卷积算子，G_σ是标准差为σ的高斯滤波器。

本发明将光流场求解与图像去偏移场过程纳入统一变分框架下，相比传统微分光流场模型，通过去偏移场，降低图像灰度不均匀性对光流计算影响，来提高配准结果精度，并且使得待配准图像中包含更多的目标图像特征信息，能够进一步恢复真实图像。通过式（8）可以得知，准确的光流场可以更好地拟合偏移场，从而在图像演化过程中，使得两者相辅相成。

我们知道，为了解决光流场孔径问题，Horn等人在光流场整体平滑基础上，提出了全局光滑约束能量项：

对式（13）进行变分，得到光流场u和ν扩散方程：

由于div(▽u)＝Δu，div(▽ν)＝Δν，全局性光流正则项其实质是各向同性的热扩散方程[10]，当该正则项作用于光流场形成过程中时，易导致光流场无方向性扩散。另外，该正则项仅仅依靠光流场信息，没有利用图像本身信息，进而缺少图像信息指引其扩散行为，尤其在配准图像特征结构处出现模糊现象，会使得配准误差较大。

由于拉普拉斯变化等效于高斯卷积，进一步分析式（14）可知：

u＝G_σ*u，ν＝G_σ*ν （15）

其中，G_σ是标准差为的高斯核。虽然采用类似高斯光滑能够得到致密光流场，但是全局正则项没有考虑到图像信息在光流场中的限定作用。为有效利用图像结构信息，Weickert等人[10]首次将扩散张量引入到扩散PDE中，这一方法的思想利用结构张量作为局部结构描绘子，通过其构造扩散张量，使PDE的扩散速度在图像边缘方向较大，而在垂直边缘的方向上相对较小。利用这个性质，我们将其改进并构造各向异性光流正则项。传统的结构张量定义如下：

其中，向量μ₁表示垂直图像边缘方向，向量μ₂表示图像边缘方向（对应D的特征向量）。λ_1,2为两方向上的量级（对应D的特征值）。计算如下：

由于D为对称正定矩阵，且秩为1，因此切线方向上的特征值为0。为了避免此现象，通常会对D的各个通道作高斯卷积。

根据脑图像的特点，要求配准过程具备以下几方面能力：一方面，脑图像含有的丰富纹理从而图像边界较多，这就要求配准过程要尽可能减少对运动边界造成的模糊；另一方面，图像演化要能够保持图像的不连续性，即要求扩散过程具有边缘保持能力。然而，从(16)式可以看出，由于仅仅考虑梯度信息，从而导致该结构张量易受噪声影响，针对这点，我们引入图像非局部结构信息来构造μ₁，μ₂

其中，*为卷积操作，若给定图像I＝{I{i},i∈Ω}，ω(i,j)是该图像非局部核权重，是图像像素点i邻域和像素点j的邻域相似度量，定义如下：

为归一化因子，N(i)表示中心点在像素点i的邻域，表示标准差为σ的高斯加权欧式距离。可用计算获得。

与传统的结构方向信息相比，本发明采用非局部核平滑传统法线方向。由于非局部核使用了图像结构信息，对图像细节信息以及纹理特征保持较好，并且对图像噪声有较强鲁棒性，使得所得到的结构方向更稳定。

下面分析如何定义特征向量对应的特征值，以构造适合脑图像配准的结构张量。特征值表征沿特征向量方向上扩散量的大小。为了在图像演化过程中尽量保持图像边缘，图像演化应尽可能保持沿图像边缘的方向上的扩散，而抑制垂直于边缘的方向上的扩散。依上述分析，定义新的结构张量对应的特征值如下：

其中，0＜α＜1,ε＞0，为常数。依据(18)式，当λ₁≈λ₂≈0时，表征该点位于图像同质区域，从而进行同质扩散，因而令否则，λ₁＞＞λ₂且λ₁越大表示图像亮度变化越大，所在区域越接近边缘，应抑制该处垂直于边缘方向上的扩散量。特征值之差(λ₁-λ₂)反映了图像区域的一致性，差值越小表示一致性越强；根据（19）式，当(λ₁-λ₂)→0时，表示该区域具有很强的一致性，因而为了保持各方向扩散量的一致，令；否则，为了增强图像结构的一致性，令一致性方向μ₂上的扩散量随(λ₁-λ₂)的增大而增大，当λ₁＞＞λ₂时，只沿一致性方向扩散，即

确定了特征向量和特征值，根据计算公式（16），就可以计算出新的结构张量：

将构造的各向异性结构张量引入到扩散方程（14）中，得到新的扩散方程如下：

结合（2）（8）与（21）式即可得到本发明提出的各向异性配准与偏移场耦合模型：

在采用上述改进光流场模型进行脑部MR图像配准时，需将基准图像和每次迭代配准后图像作为运动前后相邻的两帧图像，来求解模型光流场。设基准图像为I₁和待配准图像I₂，(u,ν)为进行计算所得光流场。基于光流场配准的方法如图1所示，包括如下步骤:

（1）初始化光流场u和ν以及偏移场：根据式（22）和式（23）以及式（12），计算光流场模型的光流场(u,ν)以及偏移场B。然后，根据求解的光流场(u,ν)作为配准的位移场，拟合得到配准后图像

（2）将作为新的待配准图像，基准图像不变，利用本文提出的关于ｕ，ｖ，Ｂ的公式计算光流场(u,ν)以及偏移场B，相应得到更新配准图像

（3）设置一定迭代次数或者判断两次配准图像I₃和的差值小于设置阀值，所求光流场(u,ν)为最终配准所求位移场，I₃为本发明提供的光流场模型配准结果，根据偏移场B进一步恢复真实目标图像。否则，将I₃作为下一步迭代的待配准图像，进入步骤（2），直到满足步骤（3）。

为了进一步说明本发明的优越性，以下将列举实例，将本方法与一些现有技术进行对比，实验环境为Dell2.0GHz1GB RAM计算机上Matlab7.0，。实验对象为虚拟脑数据以及真实脑数据。虚拟脑数据来自于mcgill库，该脑图像库可以提供不同噪声水平、灰度不均匀水平（intensity inhomogeneity,INU）的数据，是目前常用的脑MR图像分析标准库之一。

图2是虚拟脑MR图像配准结果比较图。图像噪声水平为0%，INU水平为0%。图2.a为基准图像，图2.b为待配准图像。从图中我们可以发现图中含有低对比度区域以及细长拓扑结构组织。图2.c为Horn算法迭代300次得到的结果，可以看出由于使用热扩散方程作为正则项，导致该算法得到的目标边界不光滑，且低对比区域信息易丢失。图2.d为Lucas算法配准结果，由于使用了小邻域信息，使得该算法可以得到较为光滑的边界，然而，由于所使用的邻域信息为各向同性，导致细长拓扑结构信息丢失。图2.e为本发明算法得到的结果，由于使用各向异性正则项，使得所得到的边界较为光滑，且保持了细长拓扑结构信息。

图3是虚拟脑MR图像配准结果比较图。图3.a为基准图像，图3.b为待配准图像。基准图像噪声水平为0%，INU水平为0%，待配准图像噪声水平为0%，INU水平为100%。从图中我们可以发现待配准图像含有较强的灰度不均匀现象。图3.c为Horn算法迭代300次得到的结果，可以看出由于灰度不均匀现象的影响导致该模型陷入局部最优，从而配准失败。图3.d为Lucas算法配准结果，同样由于灰度不均匀现象的影响导致该模型配准失败。图3.e为本发明算法得到的结果，由于本模型为配准与偏移场耦合模型，在配准同时可以恢复图像灰度不均匀场，从而降低了偏移场的影响。而且使用各向异性正则项，使得所得到的边界较为光滑，且保持了细长拓扑结构信息。图3.f为本发明恢复的偏移场，可以看出本发明恢复的偏移场与真实情况较为接近。

图4是虚拟脑MR图像配准结果比较图。图4.a为基准图像，图4.b为待配准图像。基准图像噪声水平为0%，INU水平为0%，待配准图像噪声水平为3%，INU水平为80%。从图中我们可以发现待配准图像不仅含有较强的灰度不均匀现象而且含有噪声。图4.c为Horn算法迭代300次得到的结果，图4.d为Lucas算法配准结果，由于灰度不均匀现象的影响导致这两种模型配准失败。图4.e为本发明算法得到的结果，由于使用patch信息构造各向异性正则项，从而降低了噪声的影响，且在配准时可恢复图像偏移场。图4.f为本发明恢复的偏移场，从图中可以看出本发明算法可以得到较理想结果。

图5是真实脑部图像配准结果。图5.a是基准图像，图5.b是待配准图像。基准图像中含有明显偏移场现象，并且图像有细长拓扑结构以及角点信息。图5.c是Horn模型迭代300次配准结果。从配准结果可以看出，由于HS模型正则项没有考虑到图像信息，配准结果在图像角点处有明显模糊现象。图5.d是Lucas算法配准结果，由于使用各向同性结构信息，导致配准结果边界处以及低对比度区域较为模糊。图5.e为本发明算法得到的结果，由于本发明设计的非局部各向异性正则项较好保留了目标图像角点等重要结构信息，更好保持图像一致性能力，大大降低了配准图像模糊程度，且在配准过程中恢复图像偏移场，从而降低偏移场影响。图5.f为本发明恢复的偏移场，本发明恢复的偏移场与真实情况较为接近。

为进一步量化本发明算法鲁棒性，本发明使用20组虚拟脑MR图像配准结果以及10组真实脑MR图像配准结果进行分析。选用平均灰度差(Mean)和均方差(Var)以及峰值信噪比(PSNR)评价三种模型配准。其定义如下：

（1）平均灰度差(Mean)

（2）均方差(Var)

（3）峰值信噪比(PSNR)

表1配准精度比较

分析结果如表1所示。由表1可以看出，本发明方法配准结果评价参数明显优于Horn模型以及Lucas模型。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (3)

1.一种基于各向异性光流场与去偏移场的脑部MR图像配准方法，其特征在于，包括如下步骤：

步骤1计算光流场(u,ν)以及偏移场B，根据求解的光流场(u,ν)作为配准的位移场，拟合得到配准后图像，

光流场u通过下式运算：

$\frac{\∂u}{\∂t}=\mathrm{\α}div(D\▿u)-{\tilde{I}}_x({\tilde{I}}_{x}u+{\tilde{I}}_{y}v+{\tilde{I}}_t-\tilde{B})$

光流场ν通过下式运算：

$\frac{\∂v}{\∂t}=\mathrm{\α}div(D\▿v)-{\tilde{I}}_y({\tilde{I}}_{x}u+{\tilde{I}}_{y}v+{\tilde{I}}_t-\tilde{B})$

偏移场通过下式运算：

$\tilde{B}=G_{\mathrm{\σ}}*({\tilde{I}}_{x}u+{\tilde{I}}_{y}v+{\tilde{I}}_t)$

配准后图像的运算公式为：

$I_2^{new}=I_2(x+u,y+v)$

其中，I(x,y,t)为基准图像，(x,y)为像素坐标信息,t为时间信息，新的结构张量新的结构张量对应的特征值如下：

其中，0＜α＜1,ε＞0，为常数；

步骤2将作为新的待配准图像，基准图像不变，计算光流场(u,ν)以及偏移场B，相应得到更新后的配准图像

步骤3判断两次配准图像I₃和的差值小于预先设定的阀值时，得到经步骤2更新后的光流场(u,ν)为最终配准所求位移场，得到I₃为光流场模型配准结果，根据偏移场B恢复真实目标图像，当两次配准图像I₃和的差值小于预先设定的阀值时，将I₃作为下一步迭代的待配准图像，执行步骤2。

2.根据权利要求1所述的基于各向异性光流场与去偏移场的脑部MR图像配准方法，其特征在于，所述光流场u和光流场ν的公式通过如下步骤得到：

步骤101将去偏移场与光流场求解两者过程相结合，纳入统一变分框架中，提出以下能量泛函：

$E=\underset{\mathrm{\Ω}}{\∫}{\left(I\right(x,y,t)-B\×I(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y,t+\mathrm{\Δ}t\left)\right)}^{2}dxdy---\left(7\right)$

其中，I(x,y,t)为基准图像，I(x+Δx,y+Δy,t+Δt)为待配准图像，B为待恢复的偏移场，步骤102对图像采用对数变换得到下式：

$E=\underset{\mathrm{\Ω}}{\∫}{\left(\tilde{I}\right(x,y,t)-\tilde{B}-\tilde{I}(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y,t+\mathrm{\Δ}t\left)\right)}^{2}dxdy---\left(8\right)$

步骤103对采用泰勒展开式：

$\tilde{I}(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y,t+\mathrm{\Δ}t)-\tilde{I}(x,y,t)\≈{\tilde{I}}_x\frac{\mathrm{\Δ}x}{\mathrm{\Δ}t}+{\tilde{I}}_y\frac{\mathrm{\Δ}y}{\mathrm{\Δ}t}+{\tilde{I}}_t---\left(9\right)$

步骤104令综合式(7)和式(9)式，将(8)可以改写如下形式：

$E=\underset{\mathrm{\Ω}}{\∫}{(\tilde{B}-({\tilde{I}}_{x}u+{\tilde{I}}_{y}v+{\tilde{I}}_t\left)\right)}^{2}dxdy---\left(10\right)$

步骤105采用梯度下降流方法，得到偏移场B表达式如下：

$\tilde{B}={\tilde{I}}_{x}u+{\tilde{I}}_{y}v+{\tilde{I}}_t---\left(11\right)$

步骤106采用低通滤波对偏移场进行光滑，即：

$\tilde{B}=G_{\mathrm{\σ}}*\tilde{B}---\left(12\right).$

3.根据权利要求1或2所述的基于各向异性光流场与去偏移场的脑部MR图像配准方法，其特征在于：所述步骤2按照预先设定的迭代次数重复执行。