CN108022261B - 一种基于改进光流场模型的非刚性图像配准方法 - Google Patents

Abstract

本发明公开了一种基于改进光流场模型的非刚性图像配准方法：采用局部互信息法求解参考图像和浮动图像之间初始位移场；构建参考图像和浮动图像的高斯金字塔图像层，在金字塔顶层图像，利用初始位移场矫正浮动图像层，最小化能量函数，求解参考图像层和浮动图像层之间运动位移场；将上一层运动位移场作为下一层初始位移场，利用光流方程计算下一层运动位移场；直到原始图像层，将获得的运动位移场作为最终运动位移场；根据最终运动位移场，对浮动图像校正，得到配准图像。本发明可以提高光流模型对非刚性图像大位移形变的估计精度，避免光流估计中的过平滑现象，提高光流估计的准确性，保留图像边缘特征，适用于物体的运动估计、目标识别、视频压缩等领域。

Description

一种基于改进光流场模型的非刚性图像配准方法

技术领域

本发明属于计算机视觉领域，更具体的说，是涉及一种基于改进光流场模型的非刚性图像配准方法。

背景技术

光流指的是空间运动物体在观察成像平面上的像素运动的瞬时速度，是通过图像序列中像素在时间域上的变化以及相邻帧之间的相关性，来寻找上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。光流不仅包含目标物体的大量运动信息，反映了图像的变化情况，还含有目标物体的三维结构等信息，被广泛地运用于运动估计、目标识别、视频检测、图像配准等领域。光流场模型算法主要分为4类：基于梯度的方法、基于特征匹配的方法、基于能量的方法、基于相位的方法。基于梯度的方法又称为微分法，它是利用时变图像灰度的梯度函数来计算像素的速度矢量，由于该方法计算简单，同时拥有较好的估计效果，因此被广泛采用。

基于梯度方法的典型代表是Horn-Schunck光流方法，该方法在光流基本约束方程的基础上增加了全局平滑假设，进而计算出图像的光流场。值得注意的是：采用Horn-Schunck光流算法对图像进行配准时，能量函数必须满足亮度守恒和全局平滑一致性假设。然而，对于非刚性形变有时并不能够满足这两个假设条件。例如，当图像的光照变化不均匀时，采用亮度守恒估计的光流位移场的误差较大，不能很好的配准浮动图像；同时，当物体发生较大形变时，两幅图像之间的运动位移场可能会较大。此时，根据全局平滑一致性的假设，通过最小化能量函数，不能获得可靠的位移场，这会在很大程度上降低配准的准确度。另一方面，较强的全局平滑一致性假设也会导致光流的扩散在图像的边缘出现过平滑现象，不利于保持图像的边缘信息，因此，在图像的梯度信息变化较明显的地方，估计的运动位移场容易存在较大偏差。

发明内容

本发明的目的是为了克服现有技术中的不足，提供了一种基于改进光流场模型的非刚性图像配准方法，利用初始位移场预先矫正浮动图像，可以提高光流模型对非刚性图像大位移形变的估计精度，同时，增加能量函数正则项的权重函数，避免光流估计中的过平滑现象，提高光流估计的准确性，保留图像的边缘特征。

本发明的目的是通过以下技术方案实现的。

一种基于改进光流场模型的非刚性图像配准方法，包括以下步骤：

步骤一，采用局部互信息法求解参考图像和浮动图像之间的初始位移场；

步骤二，构建参考图像和浮动图像的高斯金字塔图像层，在金字塔的顶层图像，利用得到的初始位移场预先矫正浮动图像层，然后最小化光流模型的能量函数，求解参考图像层和浮动图像层之间的运动位移场，向下传递位移场值；

步骤三，将上一层得到的运动位移场作为下一层的初始位移场，再利用步骤二中采用的光流方程计算下一层的运动位移场，滤波后向下传递；

步骤四，重复步骤三，直到原始图像层，将该层获得的运动位移场作为最终的运动位移场；

步骤五，根据最终得到的运动位移场，对浮动图像进行校正，得到配准图像。

步骤一中采用局部互信息法求解参考图像和浮动图像之间的初始位移场的主要流程为：

1)对参考图像和浮动图像进行分块；

2)对参考图像中的图像块N₁，通过最大化局部互信息值，在浮动图像中找到最优匹配块N₂，将它们之间的位移作为图像块N₂中心点的位移值；

3)重复步骤2)，直到得到浮动图像中所有图像块中心点的位移值；利用局部互信息法得到的光流位移场为稀疏光流场，通过数据插值方法，获得最终的初始位移场。

局部互信息定义为：

其中，N₁和N₂分别代表参考图像I₁和浮动图像I₂的图像块；W₀＝(u₀,v₀)^T是图像块N₁和N₂之间的位移场；u₀代表水平位移，v₀代表垂直位移；H(N₁)和H(N₂+W₀)分别代表图像块N₁和图像块N₂+W₀的独立熵；H(N₁,N₂+W₀)代表图像块N₁和图像块N₂+W₀之间的联合熵；

代表图像块N₁中像素值a的独立概率分布；

代表图像块N₂+W₀中像素值b′的独立概率分布；

代表图像块N₁中像素值a和图像块N₂+W₀中像素值b′之间的联合概率分布；

通过最大化局部互信息值得到稀疏光流场，具体公式如下所示：

其中，W₀ ^*＝(u₀ ^*,v₀ ^*)^T代表图像块N₁和N₂之间的最优位移场。

步骤二中最小化能量函数的过程为：

改进的正则项如下式所示：

边缘特征匹配项定义为：

根据改进的正则项和边缘特征匹配项构建的能量函数为：

E(W)＝E_D(W)+αE_S(W)+λE_match(W)

其中，

将能量函数E(W)分别对u、v求导，并另其导数为0，得Euler-Lagrange方程，如下式所示：

其中，

其中，Ψ(s²)是惩罚函数，

ε＝0.001，s＝(x,y)^T关于横纵坐标x、y的变量；Ψ′(s²)是Ψ(s²)的一阶偏导；div(·)为散度；

是空间梯度算子，u、v分表代表水平位移和垂直位移；W＝(u,v)^T代表参考图像和浮动图像之间的位移场；X＝(x,y)^T代表图像区域Ω中的一点；W_f＝(u_f,v_f)^T是边缘匹配的位移场，u_f、v_f分表代表水平位移和垂直位移；γ是梯度守恒项和亮度守恒项之间的平衡系数，取0.3；α是正则项的权重系数，取0.8；λ是特征匹配项的权重系数，取0.3；I₁、I₂分别代表参考图像和浮动图像；u₁、v₁分表代表水平位移和垂直位移。

是对x求偏导，

是对y求偏导，

是先对x求偏导再对y求偏导，

是对y求两次偏导，

是对x求两次偏导。

与现有技术相比，本发明的技术方案所带来的有益效果是：

(1)本发明在传统光流场的模型中，加入初始位移场，采用初始位移场预先矫正浮动图像，提高了光流估计的精度，有利于矫正非刚性图像的大位移形变，解决了传统光流场模型无法对发生大位移形变的非刚性图像进行有效配准的问题；

(2)本发明将改进的正则项加入到传统光流模型的能量函数中，可以有效避免光流估计中的过平滑现象，有利于保留图像的边缘特征，避免光流过平滑现象，有着广泛的应用前景；

(3)本发明将边缘特征匹配项加入到传统光流模型的能量函数中，可以较好地保留图像的内部细节特征。

附图说明

图1为本发明基于改进光流场模型的非刚性图像配准方法的流程图。

图2为本发明中基于局部互信息法计算初始位移场的流程图。

图3～图5为本发明针对几种非刚性图像测试图片的配准结果对比。

其中，图3～图5中图(a)为参考图像及浮动图像；图(b)为采用Brox算法获得的配准结果及其差值图像；图(c)是利用LDOF算法获得的配准结果及差值图像；图(d)是利用Sun算法获得的配准结果及差值图像；图(e)是利用本发明算法获得的配准结果及差值图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步描述。为了解决背景技术中存在的问题，本发明尝试利用局部互信息法求解初始位移场，对浮动图像做预处理矫正，以提高光流场算法对大位移形变的配准精度；同时利用增加的改进的正则项对光流估计进行约束，保留图像的边缘特征，避免过平滑现象。

本发明基于改进光流场模型的非刚性图像配准方法，针对传统光流场模型无法对大位移形变进行有效配准，以及光流估计中容易出现过平滑问题进行了改进，具体过程如下：

(一)采用局部互信息法求解参考图像和浮动图像之间的初始位移场，以矫正浮动图像的大位移形变，即通过估计初始位移场来替代原来由粗到细框架中的零初始位移场。

采用局部互信息法求解初始位移场的主要流程为：

1)对参考图像和浮动图像进行分块；为了采用线性模型近似地模拟图像中的非线性变形，同时兼顾求解的运算速度，本发明中可设图像块的大小为16×16

2)对参考图像中的图像块N₁，通过最大化局部互信息值，在浮动图像中找到最优匹配块N₂，将它们之间的位移作为图像块N₂中心点的位移值；

3)重复步骤2)，直到得到浮动图像中所有图像块中心点的位移值；利用局部互信息法得到的光流位移场为稀疏光流场，通过数据插值方法，获得最终稠密的初始位移场。

本发明通过最优化参考图像和浮动图像之间的局部互信息值，获得参考图像和浮动图像之间的初始位移场。其中局部互信息定义为：

其中，N₁和N₂分别代表参考图像I₁和浮动图像I₂的图像块；W₀＝(u₀,v₀)^T是图像块N₁和N₂之间的位移场；u₀代表水平位移，v₀代表垂直位移；H(N₁)和H(N₂+W₀)分别代表图像块N₁和图像块N₂+W₀的独立熵；H(N₁,N₂+W₀)代表图像块N₁和图像块N₂+W₀之间的联合熵；

代表图像块N₁中像素值a的独立概率分布；

代表图像块N₂+W₀中像素值b′的独立概率分布；

代表图像块N₁中像素值a和图像块N₂+W₀中像素值b′之间的联合概率分布；

通过最大化局部互信息值得到稀疏光流场，具体公式如下所示：

其中，W₀ ^*＝(u₀ ^*,v₀ ^*)^T代表图像块N₁和N₂之间的最优位移场。

为了提高光流位移场估计的精度，在图像块位移的求解过程中，只将获得的位移场值作为中心点的位移。对于非刚性图像来说，不可能每一点的位移都是相同的，因此，通过对稀疏位移场进行数据插值来获得稠密位移场更具有实际意义；同时，可以有效避免由于图像块中所有像素点的位移都相等而产生的块效应。

(二)基于初始位移的由粗到细的高斯金字塔分层迭代策略求解参考图像层和浮动图像层之间的运动位移场：对参考图像和浮动图像分别构建高斯金字塔图像层，金字塔的图像层数可设置为4；在金字塔的顶层图像，利用得到的初始位移场预先矫正浮动图像层；最小化光流模型的能量函数，求解参考图像层和浮动图像层之间的运动位移场，向下传递位移场值。

传统的光流场模型算法中的正则项采用全局平滑一致性假设，即假设图像的位移变化是平缓的，并采用L2范数进行度量，L2范数对图像的约束作用较强，能够较好的保留图像的边缘信息。然而，当图像发生较大变形时，容易出现比较严重的块效应。实际非刚性图像的形变一般较大，采用L1范数对正则项进行约束，可以允许图像出现较大的形变，但L1范数约束力较弱，容易出现光流过平滑的现象，导致图像的边缘细节特征丢失。为解决上述问题，本发明提出一个改进的正则项，将权重函数加入到传统光流模型的正则项中，完成本发明能量函数正则项的构建，如下式所示：

其中，Ψ(s²)是惩罚函数，

ε＝0.001，s＝(x,y)^T关于横纵坐标x、y的变量；

是空间梯度算子，u、v分表代表水平位移和垂直位移；W＝(u,v)^T代表参考图像和浮动图像之间的位移场；X＝(x,y)^T代表图像区域Ω中的一点。

在图像内部时，梯度信息变化较弱，正则函数

是L2范数约束，对图像有较强的约束力，可以保证图像内部的平滑性；在图像的边缘时，梯度信息变化较大，

是L1范数约束，允许图像发生较大的位移，可以实现较大形变的配准。

传统光流模型及其由粗到细的金字塔图像求解策略，在光流迭代的过程中容易丢失图像的细节特征。鉴于边缘特征能够处理两幅图像间的旋转、平移、缩放等变换，进行图像配准时具有较强的鲁棒性。因此，本发明在传统的光流场算法中，加入边缘特征守恒项，以保持图像的细节特征信息；同时，由于边缘匹配的限制，还有助于防止图像边缘出现过平滑的现象。边缘特征匹配项定义下式所示：

E_match(W)＝∫_ΩΨ|W-W_f|²dX (7)

其中，W_f＝(u_f,v_f)^T是边缘匹配的位移场，u_f、v_f分表代表水平位移和垂直位移。

将边缘特征匹配项加入到传统的能量函数中，并且增加改进的正则项，构成本发明的能量函数，如下式所示：

E(W)＝E_D(W)+αE_S(W)+λE_match(W)

(8)

其中，

式中，γ是梯度守恒项和亮度守恒项之间的平衡系数，本发明中取0.3；α是正则项的权重系数，本发明中取0.8；λ是特征匹配项的权重系数，本发明中取0.3；I₁、I₂分别代表参考图像和浮动图像。

采用基于初始位移场的改进光流场算法对非刚性图像进行配准，该问题可以归结为最小化上式的能量函数。将能量函数E(W)分别对u、v求导，并另其导数为0，得Euler-Lagrange方程，如下式所示：

其中，

式中，Ψ′(s²)是Ψ(s²)的一阶偏导；div(·)为散度。u₁、v₁分表代表水平位移和垂直位移。

是对x求偏导，

是对y求偏导，

是先对x求偏导再对y求偏导，

是先对y求偏导再对x求偏导，

是对y求两次偏导，

是对x求两次偏导。I_x代表图像I₂(X+W)对x求偏导的结果，I_y代表图像I₂(X+W)对y求偏导的结果，I_z是图像I₂(X+W)和I₁(X)之间的差值图像，I_xy是图像I₂(X+W)先对x求偏导再对y求偏导的结果，I_yy是图像I₂(X+W)对y求两次偏导的结果，I_xz是图像I_z对x求偏导的结果，I_xx是图像I₂(X+W)对x求两次偏导的结果，I_yz是图像I_z对y求偏导的结果。

(三)将上一层金字塔图像得到的运动位移场作为下一层的初始位移场，再继续光流迭代，利用步骤(二)中采用的光流方程公式(6)-(11)计算下一层的运动位移场，滤波后向下传递。

也就是说，将上一层的计算结果插值作为初始值传递给下一层图像，利用这个位移初始值对下一层的浮动图像进行变换，得到变换后的浮动图像，利用参考图像层和变换后的浮动图像层，求解下一层的运动位移场。

(四)重复步骤(三)，在金字塔图像层之间不断的传递位移值，并且变换浮动图像层，直到金字塔的最底层，即原始图像层，并将这一层获得的运动位移场作为最终的运动位移场。

(五)根据最终得到的运动位移场，对浮动图像进行校正，得到配准图像。

具体实施例：

下面结合具体的实例对本发明的技术方案做进一步详细描述。

实验结果均为本发明在CPU为Intel i5-4590，3.3GHz，内存为8G的台式电脑上运行所得，操作系统为Windows 7，仿真软件为64位Matlab R2015b。图3～图5为本发明针对非刚性图像测试图片的配准结果对比示意图。主要参数选择为：γ＝0.3，α＝0.8，λ＝0.3。光流迭代次数为60次。

图3是第一组大位移形变非刚性图像配准效果对比，其中，图(a)为参考图像及浮动图像；图(b)为采用Brox算法获得的配准结果及差值图像；图(c)是利用LDOF算法获得的配准结果及差值图像；图(d)是利用Sun算法获得的配准结果及差值图像；图(e)是利用本发明算法获得的配准结果及差值图像。从图3中可以看出，Brox算法可以对非刚性图像中较小的形变位移进行矫正，例如：在浮动图像的内部区域，其配准结果和参考图像之间的差值较小；但是对于浮动图像边缘较大的形变区域，却出现了较大的误差。LDOF算法是针对于大位移形变的图像配准算法，因此，图像的边缘区域得到了很好的矫正；但是，由于正则项的约束能力较弱，导致配准图像的内部区域出现了过平滑现象，例如：左下角区域的配准效果不够理想。Sun算法采用了约束性较强的正则项，在图像内部的小位移形变得到了很好的配准，没有出现过平滑现象；但是，在图像的边缘区域，由于形变位移较大，配准图像与参考图像之间存在较大的误差。本发明算法采用了局部互信息法，首先对初始位移场进行估计，可以有效的矫正非刚性图像的大位移形变，在图像边缘区域获得了比较理想的配准效果。此外，由于本文算法能量函数的正则项中加入了权重函数，可以有效地保证图像内部区域的平滑，同时保留图像的边缘细节特征，没有出现过平滑的现象。相较于其它算法，本文算法估计得到的位移场更加准确，配准效果更好。

图4是第二组大位移形变非刚性图像配准效果对比，其中，图(a)为参考图像及浮动图像；图(b)为采用Brox算法获得的配准结果及差值图像；图(c)是利用LDOF算法获得的配准结果及差值图像；图(d)是利用Sun算法获得的配准结果及差值图像；图(e)是利用本发明算法获得的配准结果及差值图像。从图4中可以看出，由于该组实验浮动图像和参考图像之间的位移形变较小，所以几种光流配准算法都获得了比较理想的配准效果。但是可以看出，Brox算法、LDOF算法、Sun算法在图像的边缘区域都产生了不同程度的误差。而本发明算法，由于用改进的正则项进行约束，较好的保留了图像的边缘细节特征，同时由于初始位移场的加入，浮动图像发生较大形变的边缘图像区域得到很好的矫正，与参考图像之间的误差较小，取得了较好的配准效果。

图5是第三组大位移形变非刚性图像配准效果对比，其中，图(a)为参考图像及浮动图像；图(b)为采用Brox算法获得的配准结果及差值图像；图(c)是利用LDOF算法获得的配准结果及差值图像；图(d)是利用Sun算法获得的配准结果及差值图像；图(e)是利用本发明算法获得的配准结果及差值图像。从配准结果图5中可以看出，对于浮动图像和参考图像之间的较大位移形变，Brox算法几乎不起作用，Sun算法的配准效果也很差；这是由于Brox算法和Sun算法采用全局平滑假设，因此不适合处理非刚性图像的大位移形变。LDOF算法在整体上对大位移形变进行了矫正，但是对于图像右下角的边缘，光流没有得到很好的扩散，得到配准结果与参考图像之间的误差较大；而配准结果图像的最上角出现了光流的过平滑现象，导致图像的边缘没有得到保持。与之相反，由于本发明算法加入了带有权重函数的正则项，允许光流在边缘区域扩散较快，同时在图像内部区域，保持图像的平滑性，因此，没有出现过平滑的现象，较好的保留了图像的边缘细节特征。同时，由于初始位移场的加入，对发生大位移形变的非刚性图像进行了预先的矫正，再采用由粗到细的金字塔模型进行光流迭代，得到的光流估计结果精度较高。因此相比于其它算法，本发明算法的配准效果更好，证明了本发明所提出的基于初始位移场的改进光流场算法的有效性，能够处理非刚性图像之间的大位移形变。

尽管上面结合附图对本发明进行了描述，但本发明并不局限于上述，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，还可以做出很多形式，这些均属于本发明的保护之内。

Claims (3)

1.一种基于改进光流场模型的非刚性图像配准方法，其特征在于，包括以下步骤：

步骤一，采用局部互信息法求解参考图像和浮动图像之间的初始位移场；

步骤二，构建参考图像和浮动图像的高斯金字塔图像层，在金字塔的顶层图像，利用得到的初始位移场预先矫正浮动图像层，然后最小化光流模型的能量函数，求解参考图像层和浮动图像层之间的运动位移场，向下传递位移场值；

其中，最小化能量函数的过程为：

改进的正则项如下式所示：

边缘特征匹配项定义为：

E_match(W)＝∫_ΩΨ|W-W_f|²dX

根据改进的正则项和边缘特征匹配项构建的能量函数为：

E(W)＝E_D(W)+αE_S(W)+λE_match(W)

其中，

将能量函数E(W)分别对u、v求导，并另其导数为0，得Euler-Lagrange方程，如下式所示：

其中，

其中，Ψ(s²)是惩罚函数，

ε＝0.001，s＝(x,y)^T关于横纵坐标x、y的变量；Ψ′(s²)是Ψ(s²)的一阶偏导；div(·)为散度；

是空间梯度算子，u、v分表代表水平位移和垂直位移；W＝(u,v)^T代表参考图像和浮动图像之间的位移场；X＝(x,y)^T代表图像区域Ω中的一点；W_f＝(u_f,v_f)^T是边缘匹配的位移场，u_f、v_f分表代表水平位移和垂直位移；γ是梯度守恒项和亮度守恒项之间的平衡系数，取0.3；α是正则项的权重系数，取0.8；λ是特征匹配项的权重系数，取0.3；I₁、I₂分别代表参考图像和浮动图像；u₁、v₁分表代表水平位移和垂直位移；

是对x求偏导，

是对y求偏导，

是先对x求偏导再对y求偏导，

是对y求两次偏导，

是对x求两次偏导；

步骤三，将上一层得到的运动位移场作为下一层的初始位移场，再利用步骤二中采用的光流方程计算下一层的运动位移场，滤波后向下传递；

步骤四，重复步骤三，直到原始图像层，将该层获得的运动位移场作为最终的运动位移场；

步骤五，根据最终得到的运动位移场，对浮动图像进行校正，得到配准图像。

2.根据权利要求1所述的基于改进光流场模型的非刚性图像配准方法，其特征在于，步骤一中采用局部互信息法求解参考图像和浮动图像之间的初始位移场的主要流程为：

1)对参考图像和浮动图像进行分块；

2)对参考图像中的图像块N₁，通过最大化局部互信息值，在浮动图像中找到最优匹配块N₂，将它们之间的位移作为图像块N₂中心点的位移值；

3)重复步骤2)，直到得到浮动图像中所有图像块中心点的位移值；利用局部互信息法得到的光流位移场为稀疏光流场，通过数据插值方法，获得最终的初始位移场。

3.根据权利要求2所述的基于改进光流场模型的非刚性图像配准方法，其特征在于，局部互信息定义为：

其中，N₁和N₂分别代表参考图像I₁和浮动图像I₂的图像块；W₀＝(u₀,v₀)^T是图像块N₁和N₂之间的位移场；u₀代表水平位移，v₀代表垂直位移；H(N₁)和H(N₂+W₀)分别代表图像块N₁和图像块N₂+W₀的独立熵；H(N₁,N₂+W₀)代表图像块N₁和图像块N₂+W₀之间的联合熵；

代表图像块N₁中像素值a的独立概率分布；

代表图像块N₂+W₀中像素值b′的独立概率分布；

代表图像块N₁中像素值a和图像块N₂+W₀中像素值b′之间的联合概率分布；

通过最大化局部互信息值得到稀疏光流场，具体公式如下所示：

其中，W₀ ^*＝(u₀ ^*,v₀ ^*)^T代表图像块N₁和N₂之间的最优位移场。

Images