CN104240230A - 一种提高相位相关算法匹配精度的方法 - Google Patents

Abstract

本发明公开了一种提高相位相关算法匹配精度的方法，包括如下步骤：（1）将参考图像与目标图像乘以改进的窗口函数；（2）对（1）所得的乘积结果进行傅里叶变换，得到频谱，在频谱上叠加权重函数；（3）利用相位相关算法，获得参考图像与目标图像之间的整数位移矢量；（4）利用最小二乘法，对整数峰值周围3×3邻域进行二次曲面拟合；（5）对二次曲面拟合求导，获得小数位移矢量；（6）整数位移矢量加上小数位移矢量，所得和为目标图像与参考图像之间的亚像素级位移矢量。本发明方法可以提高图像间的匹配精度，在图像识别、图像匹配、稳定成像领域具有重大意义。

Description

一种提高相位相关算法匹配精度的方法

技术领域

本发明涉及图像处理技术领域，特别是基于相位相关算法的匹配精度提高方法。 

背景技术

亚像素级别的图像配准是一项重要的图像预处理技术，广泛应用于如遥感、高精度3D重建、视觉定位、医学图像等领域中的许多关键问题。相位相关法是一种常用的图像配准方法，由于它对噪声鲁棒性强，实现简单，测量精度高，计算量小，代价-精度比最高而广受关注。传统的相位相关法通过改进可达到亚像素配准精度，其目前发展方向大致有三个：一是寻找变换参数新的频域对应特征；二是发展各种加权、频率选择算法，减少图像的离散化和噪声的影响，进一步提高精度；三是与其他方法结合，如频率求解的扩展相位相关法与解最优化问题法互补，可克服后者对噪声敏感的弱点。 

扩展相位相关分为空域求解法和频域求解法。空域求解的相位相关扩展具有配准精度高，计算量小的优点，需要先用像素级别的配准算法达到像素级别精度后才能使用，遵循从粗到精的配准过程。其中，影响亚像素配准精度的因素包括由傅里叶变换的周期性引起的边缘效应，和采样过程中尤其是对于非带限信号所引入的频谱混叠。边缘效应的处理方法是对输入图像加窗口函数。为了抑制混叠效应，可对频谱进行加权，强调信噪比更高的低通分量。 

发明内容

本发明的目的在于提供一种基于相位相关算法提高匹配精度的方法，即提高图像间的位移矢量精度的方法。本发明提出了一种平顶窗口函数，将平顶窗口函数叠加至输入图像，并对输入图像的功率谱叠加频谱权重函 数，对所得的相关峰值模型的二维的数值矩阵用移动最小二乘法进行曲面拟合，得到两幅图像之间的位移矢量。 

本发明通过如下技术方案实现。 

一种提高相位相关算法匹配精度的方法，其包括如下步骤： 

(1)将参考图像与目标图像乘以改进的窗口函数； 

(2)对(1)所得的乘积结果进行傅里叶变换，得到频谱，在频谱上叠加权重函数； 

(3)利用相位相关算法，获得参考图像与目标图像之间的整数位移矢量； 

(4)利用最小二乘法，对整数峰值周围3×3邻域进行二次曲面拟合； 

(5)对二次曲面拟合求导，获得小数位移矢量； 

(6)整数位移矢量加上小数位移矢量，所得和为目标图像与参考图像之间的亚像素位移矢量。 

进一步地，步骤(1)中所述的改进的窗口函数为平顶窗口函数。所述平顶窗口函数为其中  $\mathrm{\Γ}(n_1,n_2)=k\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_1}{M}\right)\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_2}{N}\right)),$ M，N为图像的尺寸，k为窗口函数的拉伸系数，n₁，n₂代表图像的两个方向坐标。 

进一步地，步骤(2)中所述的权重函数为高斯函数或者理想低通滤波器。 

进一步地，步骤(4)中所述的最小二乘法为移动最小二乘法。 

本发明基本原理如下： 

已知两张连续2D图像的空间采样尺寸为M×N的图像f(n₁,n₂)和g(n₁,n₂)，若g(n₁,n₂)相对f(n₁,n₂)分别存在沿n₁、n₂方向的亚像素位移(包含整数位移矢量+小数位移矢量，如3.5像素，3为整数位移矢量，0.5为小数位移矢量)δ₁、δ₂，则g(n₁,n₂)＝f(n₁-δ₁,n₂-δ₂)。它们的二维离散傅里叶变换为F(k₁,k₂)和G(k₁,k₂)，k₁,k₂为频谱空间的变量。考虑到傅里叶变换定义在连续空间和离散空间中的性质区别，F(k₁,k₂)和G(k₁,k₂)的关系式可近似写作 

$g(k_1,k_2)\≅F(k_1,k_2)\·e^{-j\frac{2\mathrm{\π}}{M}{k}_1{\mathrm{\δ}}_1}{e}^{-j\frac{2\mathrm{\π}}{N}{k}_2{\mathrm{\δ}}_2}---\left(1\right)$

F(k₁,k₂)和G(k₁,k₂)的交叉功率谱R(k1,k2)可表达为： 

$R(k_1,k_2)=F(k_1,k_2)\stackrel{\‾}{G(k_1,k_2)}---\left(2\right)$

表示为G(k₁,k₂)的共轭项,交叉相位谱是归一化的交叉功率谱的，定义为： 

$\hat{R}(k_1,k_2)=\frac{F(k_1,k_2)\stackrel{\‾}{G(k_1,k_2)}}{\left|F(k_1,k_2)\stackrel{\‾}{G(k_1,k_2)}\right|}\≅e^{j\frac{2\mathrm{\π}}{M}{k}_1{\mathrm{\δ}}_1}{e}^{j\frac{2\mathrm{\π}}{N}{k}_2{\mathrm{\δ}}_2}$ (3)相位相关(POC)函数就是的二维离散傅里叶逆变换，可表达为： 

$\begin{array}{c}\hat{r}(n_1,n_2)=\frac{1}{\mathrm{MN}}\underset{k_1{k}_2}{\mathrm{\Σ}}\hat{R}(k_1,k_2)e^{j\frac{2\mathrm{\π}}{M}{k}_1{n}_1}{e}^{j\frac{2\mathrm{\π}}{N}{k}_2{n}_2}\\ \≅\frac{\mathrm{\α}}{\mathrm{MN}}\frac{\sin \left\{\mathrm{\π}(n_1,{\mathrm{\δ}}_1)\right\}}{\sin \left\}\right\{\frac{\mathrm{\π}}{M}(n_1+{\mathrm{\δ}}_1)\}}\frac{\sin \left\{\mathrm{\π}(n_2+{\mathrm{\δ}}_2)\right\}}{\sin \left\{\frac{\mathrm{\π}}{N}(n_2+{\mathrm{\δ}}_2)\right\}}\end{array}---\left(3\right)$

当原图像中无噪声时，α＝1。当有小噪声分量掺入到原图像中时，峰值α将会减小，但不会改变函数原本的形状，因此，α≤1。式(4)代表了一般图像相互间有亚像素级别位移时POC函数的峰值形状，其峰值的位置对应两图像之间的位移量。 

基于POC的图像位移矢量测量的精度会受到边缘效应的影响，这是由于DFT(离散傅里叶变换)的周期性会使得图像因为每条边界的不同而出现不连续性，人为地在图像频谱引入高频干扰，而且两幅图像因微小位移而在边界部分引入的新信息也会使频谱也发生畸变。为了抑制DFT的边缘效应，可对输入图像上加窗口函数。 

窗口函数相当于在空间域内的滤波，本质是通过中心加权的方法降低图像边缘部分的影响，从而减小边缘信息对变换后图像频域的干扰。加在图像上的窗口函数应该是平滑连续的，且具有良好的中心选择性，这样才能在既削弱边缘干扰的同时又不引入新的截断误差。图像配准常用的窗口 函数有Hanning窗，Hamming窗，Blackman窗。但是对于背景较单一、目标物体较小而且信息集中的图像，上述几种经典窗口对配准精度的提高是不利的，因为窗口函数通带内不均的加权可能会大大改变原图像信息。此外，单一背景的图像边缘效应本身不显著，上述窗口的通带范围不可调节则不能对图像信息进行有效的截取。 

为此，本发明构造一个改进的平顶窗口函数，它的中间部分相当于一个矩形窗，而在边缘附近过渡为汉宁窗(Hanning窗)，为了调节窗口函数的通带范围引入参数k，表达式可写作： 

(4)其中  $\mathrm{\Γ}(n_1,n_2)=k\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_1}{M}\right)\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_2}{N}\right)),$ M，N为图像的尺寸，k为窗口函数的拉伸系数。 

改进窗口函数既保证了图像中心信息的不失真，又继承了Hanning窗在边缘连续性的过渡而控制了截断误差。拉伸系数k的引入使得该窗口函数的通带范围可调，适应性更广，通过实验改变k可调整窗口函数的性能，最优化k值使得测量精度更理想。 

在频域上添加加权函数主要有两个功能：一是滤除高频畸变分量，二是改变用于拟合的POC函数的峰值模型。自然图像的大部分能量普遍集中在低频率分量中。图像因边缘效应造成的高频畸变，由加窗或者非带限信号采样带来的高频混叠都会使得高频分量相比低频分量具有更低的信噪比。同时，在频域乘上一个加权函数，相当于在空间域内做卷积运算，这样会相应地改变峰值模型。因此，应选择可能实现高精度的函数拟合的峰值模型z(x,y)，以及能有效滤除高频信号的滤波器H(k₁,k₂)。此外，在选择了z(x,y)和H(k₁,k₂)之后，需要最优化用于函数拟合的数据点的数量和加 权函数H(k₁,k₂)的通带带宽。 

利用移动最小二乘法(MLSQ)来拟合二次曲面峰值模型z(x,y)＝a₂₀x²+a₁₁xy+a₀₂y²+a₁₀x+a₀₁y+a₀₀的算法如下： 

在拟合窗口内，取函数f(x,y)为函数z(x,y)的逼近函数： 

$f(x,y)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_i(x,y)p_i(x,y)=p^T(x,y)\mathrm{\α}(x,y)---\left(6\right)$

其中α(x,y)＝[α₁(x,y),α₂(x,y),Λ,α_m(x,y)]^T为待定系数，它是坐标(x,y)的函数。p(x,y)＝[p₁(x,y),p₂(x,y),Λ,p_m(x,y)]^T称为基函数，它是一个k阶完备多项式，m是基函数的个数。对于二维曲面的二次基，p＝[1,x,y,x²,xy,y²]^T，m＝6。 

α(x,y)由加权最小二乘法确定，它应使用f(x,y)局部逼近z(x,y)误差最小，即它们的n维的加权离散L₂范式最小，可写作： 

$\begin{array}{c}J=\underset{I=1}{\overset{n}{\mathrm{\Σ}}}{w}_I\left(s\right){[f(x,y)-z(x_I,y_I)]}^2\\ =\underset{I=1}{\overset{n}{\mathrm{\Σ}}}{w}_I\left(s\right){[p^T(x_I,y_I)\mathrm{\α}(x,y)-z(x_I,y_I)]}^2\end{array}---\left(7\right)$

其中n是紧支域内节点的数目，(x_I,y_I)是紧支域内的节点，w_I(s)是节点(x_I,y_I)的权函数， $s=\sqrt{{(x-x_I)}^2+{(y-y_I)}^2}.$

假设构造矩阵P和W(s)，它们的表达式为： 

则式(7)可用矩阵形式表达为： 

J＝(Pα-z)^TW(s)(Pα-z)    (10) 

其中z＝(z₁,z₂,...,z_n)^T，z_n＝z(x_n,y_n)。 

为了得到α(x,y)，应对J取极值，也就是： 

$\frac{\∂J}{\∂\mathrm{\α}}=A(x,y)\mathrm{\α}(x,y)-B(x,y)z=0$

A(x,y)α(x,y)＝B(x,y)z         (11) 

其中， 

A(x,y)＝P^TW(s)P           (12) 

B(x,y)＝P^TW(s)            (13) 

根据式(11)进行矩阵运算后，可得： 

α(x,y)＝A^-1(x,y)B(x,y)z            (14) 

那么，拟合函数f(x,y)可表达为： 

$f(x,y)=\mathrm{\Φ}(x,y)z=\underset{I=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Φ}}_I(x,y)z_I---\left(15\right)$

其中，Φ(x,y)为形函数，表达式为： 

Φ(x,y)＝[Φ₁(x,y),Φ₂(x,y),...Φ_n(x,y)]＝p^T(x,y)A^-1(x,y)B(x,y)。 

权函数取三次样条权函数， 

$w\left(\stackrel{\‾}{s}\right)=\left\{\begin{array}{cc}\frac{2}{3}-4{\stackrel{\‾}{s}}^2& (\stackrel{\‾}{s}\≤\frac{1}{2})\\ \frac{4}{3}-4\stackrel{\‾}{s}+4{\stackrel{\‾}{s}}^2-\frac{4}{3}{\stackrel{\‾}{s}}^3& (\frac{1}{2}\≤\stackrel{\‾}{s}\≤1)\\ 0& (\stackrel{\‾}{s}>1)\end{array}\right.---\left(16\right)$

其中 $\stackrel{\‾}{s},\frac{s}{s_{\max }},s=\sqrt{{(x-x_I)}^2+{(y-y_I)}^2},$ 这里取圆形影响域，其半径记为s_max。 

与现有技术相比，本发明具有如下优点和技术效果：本发明提出了一个改进窗口函数，结合空间域加窗口函数和频域加权的扩展相位相关法，对所得的相关峰值模型的二维的数值矩阵进行移动最小二乘法曲面拟合，并通过变化参数值使得精度最优化。 

附图说明

图1为实施方式中的方法流程图； 

图2a图2b分别为实施方式中的参考图像与目标图像； 

图3为实施方式中的改进平顶窗口函数图； 

图4a、图4b分别为改进的移动最小二乘法(MLSQ)拟合原理图和流程图； 

图5为实施方式中得到的整数像素值的相关峰图； 

图6实施方式中利用移动最小二乘法拟合峰值点获得的峰值曲面。 

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明，但本发明的实施不限于。需指出的是，以下若有未特别详细说明的参量或过程，均是本领域技术人员可以参照现有技术实现的。 

本发明的流程图如图1所示，具体步骤如下： 

(1)、输入参考图像与目标图像，所用的参考图像如图2a所示，目标图像如图2b所示,两幅图像之间存在已知的位移矢量。将参考图像与目标图像乘以改进的平顶窗口函数，平顶窗口函数为 

其中 $\mathrm{\Γ}(n_1,n_2)=k\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_1}{M}\right)\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_2}{N}\right)),$ M，N为图像的尺寸，k为窗口函数的拉伸系数。本实例中k＝2，窗口函数如图3所示。 

(2)、对(1)所得的乘积结果分别进行离散傅里叶变换，得到其各自的频谱，在各自频谱上叠加权重函数，对参数进行全局优化。采用的权重函数为高斯函数H₂(k₁,k₂)＝exp(-(k₁ ²+k₂ ²)/2σ²)，其中σ是高斯函数的参数，k₁,k₂代表图像的两个方向坐标，控制函数的宽度。本实例中，最优化σ取值为28； 

(3)、根据公式(2)获得参考图像与目标图像的交叉功率谱，根据公式(3)获得其归一化联合功率谱，根据公式(4)对归一化联合功率谱进行逆向离散傅里叶变换，获得参考图像与目标图像之间的整数像素级别位移矢量，如图5； 

(4)、对整数像素峰值坐标周围3×3邻域，根据相关峰周围的9个点，利用移动最小二乘法对整数像素相关峰进行曲面拟合，移动最小二乘 法算法及流程图如图4a、图4b所示，得到二次曲面如图6所示； 

下面具体对步骤(4)作进一步说明： 

我们令第一次拟合时的影响域圆心是由式(4)计算出的像素级别的匹配点，根据描述的移动最小二乘法原理式(7)-式(14)求出系数矩阵α_I(x,y)。把α_I(x,y)看作是常系数代入二次曲面方程式z(x,y)＝a₂₀x²+a₁₁xy+a₀₂y²+a₁₀x+a₀₁y+a₀₀，按照传统的LSQ法求出峰值点坐标，设为(x_I,y_I)。 

为进一步提高计算精度，以(x_I,y_I)作为为影响域的新圆心，半径不变，再进行一次MLSQ曲面拟合，得到新系数矩阵α_II(x,y)，再重复上面的步骤可得到新的亚像素极值点坐标(x_II,y_II)，以此类推。可以设定某一阈值T，当前后两次求得的极值点坐标间的距离小于T时停止下一次的拟合。本实例设定阈值T＝0.01。若i表示拟合的次数，则相邻两次拟合所得亚像素匹配点间的距离可表达为： 

$\mathrm{\Δd}=\sqrt{{(x_i-x_{i+1})}^2+{(y_i-y_{i+1})}^2}.$

(5)、对拟合的二次曲面求导，获得目标图像与参考图像间的亚像素位移精度； 

(6)、整数位移矢量加上小数位移矢量，所得和为目标图像与参考图像之间的亚像素级位移矢量，进而完成目标图像与参考图像之间的亚像素匹配。 

本发明中应用了各种图对具体实施方式进行了阐述，以上所述仅为本发明较佳可行的实施例子而已。对于本领域的技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改善之处。综上所述，本说明书内容不应理解为对本发明的限制。 

Claims (5)

1.一种提高相位相关算法匹配精度的方法，其特征在于，包括如下步骤：

(1)将参考图像与目标图像乘以改进的窗口函数；

(2)对(1)所得的乘积结果进行傅里叶变换，得到频谱，在频谱上叠加权重函数；

(3)利用相位相关算法，获得参考图像与目标图像之间的整数位移矢量；

(4)利用最小二乘法，对整数峰值周围3×3邻域进行二次曲面拟合；

(5)对二次曲面拟合求导，获得小数位移矢量；

(6)整数位移矢量加上小数位移矢量，所得和为目标图像与参考图像之间的亚像素级位移矢量。

2.如权利要求1所述的一种提高相位相关算法匹配精度的方法，其特征在于，步骤(1)中所述的改进的窗口函数为平顶窗口函数。

3.如权利要求2所述的一种提高相位相关算法匹配精度的方法，其特征在于，所述平顶窗口函数为其中 $\mathrm{\Γ}(n_1,n_2)=k\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_1}{M}\right)\·0.5(1-\cos \left(2\mathrm{\π}\frac{n_2}{N}\right)),$ M，N为图像的尺寸，k为平顶窗口函数的拉伸系数，n₁，n₂代表图像的两个方向坐标。

4.如权利要求1所述的一种提高相位相关算法匹配精度的方法，其特征在于，步骤(2)中所述的权重函数为高斯函数或者理想低通滤波器。

5.如权利要求1所述的一种提高相位相关算法匹配精度的方法，其特征在于，步骤(4)中所述的最小二乘法为移动最小二乘法。