CN106682689A - 基于多尺度傅里叶梅林变换的图像匹配方法 - Google Patents

Abstract

本发明公开了一种基于多尺度傅里叶梅林变换的图像匹配方法，思路为：分别获取基准图像I₁和待匹配图像I₂，确定I₂中像素位置(x,y)处的高斯核，计算I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)；进而得到待匹配图像I₂的尺度空间图像L(t)，对L(t)分别进行双线性插值，得到双线性插值后的模拟成像S(t)；设定与L(t)尺寸相同的空白图像，将S(t)中每一个像素位置处的灰度值分别嵌入至所述空白图像中，得到多尺度图像，并分别计算I₁的傅里叶梅林不变量图像和多尺度图像S'(t)的傅里叶梅林不变量图像；计算I₁与I₂的尺度空间图像L(t)之间的相位相关匹配图像，分别计算I₁中像素位置缩放的尺度估计值和I₁中像素位置旋转的角度估计值，进而得到I₂的最优匹配图像。

Description

基于多尺度傅里叶梅林变换的图像匹配方法

技术领域

本发明涉及图像处理技术领域，特别涉及一种基于傅里叶梅林变换的频率域图像匹配方法，适用于目标跟踪、目标识别、视频监控、视频会议、生物医学、信息安全、遥感遥测、航天航空或视频多媒体中的频率域图像匹配。

背景技术

目前，随着计算机视觉技术的不断普及，图像匹配技术被广泛应用于各个领域，随之对图像匹配技术的要求也在不断提高，传统的图像匹配算法，如平均绝对差(MAD)算法是针对空间域逐像素进行处理，虽然最终匹配策略的描述方式上各不相同，但算法使用前提通常是目标尺寸不发生较大变化，目标内各像素空间相对位置不发生较大变化，因此，很难适应目标缩放及目标旋转的情况；为解决上述问题，现有技术提出了宽基线图像匹配方法，宽基线图像匹配是指由两台距离较远或3维视角差异较大的摄像机拍摄的图像之间的匹配。

常提取大量的对尺度和旋转具有不变性的高维局部特征，基于所述高维局部特征进行图像匹配，如尺度不变性特征变换(SIFT)、最稳定极值区域(MSER)、快速鲁棒性特征(SURF)，具有所述高维局部特征的算法对图像质量要求较高，对纹理不丰富目标的匹配效果不理想，且运算量较大；还有一类算法是通过监督学习或训练分类器实现快速稳定的图像匹配，但训练分类器需要较长时间和较多的样本，在要求系统有较高的实时性和可靠性的情况下时，该类算法也不适用。

传统的基于傅里叶梅林变换(Fourier-Mell in Transform)的图像匹配是一种基于频率域的图像匹配算法，是将图像傅里叶变换的幅度进行对数极坐标变换的一种变换方式；利用傅里叶梅林变换将图像尺度缩放和旋转转化为该图像傅里叶变换的幅度在对数极坐标系中的平移，并根据傅里叶变换位移理论，以及基于频率域匹配方法进行参数求解，得到最终的尺度及旋转变换参数，从而适应了目标缩放和旋转的情况，并对目标平移、尺度、旋转变换具有鲁棒性，且算法运算量较小。

然而，通过对大量图像的仿真分析发现，基于传统傅里叶梅林变换的图像匹配仅对较小范围内的尺度变化具有鲁棒性，当两图像之间尺度差异较大时，匹配失效；因此，传统的傅里叶梅林变换方法不能适应尺度差异较大的情况。

发明内容

针对上述现有技术存在的不足，本发明的目的在于提出一种基于多尺度傅里叶梅林变换的图像匹配方法，该种基于多尺度傅里叶梅林变换的图像匹配方法是对传统基于傅里叶梅林变换的图像匹配方法的改进，并基于高斯尺度空间理论，能够扩大傅里叶梅林变换在尺度差异方面的适应范围。

为达到上述技术目的，本发明采用如下技术方案予以实现。

一种基于多尺度傅里叶梅林变换的图像匹配方法，包括以下步骤：

步骤1，分别获取基准图像I₁和待匹配图像I₂，所述基准图像I₁和待匹配图像I₂分别为灰度图像，且基准图像I₁的尺寸为M×N，待匹配图像I₂的尺寸为P×Q，M＝P，N＝Q，M、N、P、Q分别为大于0的自然数；

然后确定待匹配图像I₂中像素位置(x,y)处的高斯核为G(x,y,t)，并计算得到待匹配图像I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)，(x,y)∈I₂；进而得到待匹配图像I₂的尺度空间图像L(t)，t表示尺度参数，t>0；

步骤2，对待匹配图像I₂的尺度空间图像L(t)分别进行双线性插值，得到双线性插值后的模拟成像S(t)；

步骤3，设定与待匹配图像I₂的尺度空间图像L(t)尺寸相同的空白图像，并将双线性插值后的模拟成像S(t)中每一个像素位置处的灰度值分别嵌入至所述空白图像中，得到多尺度图像S'(t)；所述多尺度图像S'(t)包含若干个尺度图像，且所述多尺度图像S'(t)的尺寸为X×Y；X、Y分别为大于0的自然数；

步骤4，对所述基准图像I₁和多尺度图像S'(t)分别进行快速傅里叶变换，分别得到快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)，进而分别计算得到基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)；

步骤5，根据基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)，计算得到基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)；

步骤6，根据所述基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)，分别计算得到基准图像I₁中像素位置缩放的尺度估计值和基准图像I₁中像素位置旋转的角度估计值其中所述基准图像I₁中像素位置对应待匹配图像I₂中像素位置(x,y)；将所述基准图像I₁中每个像素分别进行尺度估计值的缩放，同时将所述基准图像I₁中每个像素位置进行角度估计值的旋转后，得到待匹配图像I₂的最优匹配图像。

本发明与现有的技术相比具有以下优点：基于传统的傅里叶梅林变换的图像匹配仅能适应图像尺度差异在[0.7，1.4]之间的情况，当尺度差异超出此范围时，该算法无法正确估计尺度变换参数；本发明基于高斯尺度空间构造了一种新的基于多尺度傅里叶梅林变换的图像变换参数估计策略，即基于多尺度傅里叶梅林变换的图像匹配方法，所述方法能够将尺度差异范围扩展到[0.35，2.8]，提高了对图像尺度差异的鲁棒性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供一种基于多尺度傅里叶梅林变换的图像匹配方法流程图；

图2为本发明实施例提供的基准图像示意图；其中，该基准图像为经典的Lena图，尺寸为128×128；

图3为本发明实施例提供的待匹配图像示意图；

图4为本发明实施例提供的傅里叶梅林不变量生成图的源图像；

图5为图4的傅里叶变换幅度谱示例图；

图6为图4的傅里叶梅林不变量示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1，为本发明实施例提供一种基于多尺度傅里叶梅林变换的图像匹配方法流程图；所述基于多尺度傅里叶梅林变换的图像匹配方法，包括以下步骤：

步骤1，分别获取基准图像I₁和待匹配图像I₂，所述基准图像I₁和待匹配图像I₂分别为灰度图像，且基准图像I₁的尺寸为M×N，待匹配图像I₂的尺寸为P×Q，M＝P，N＝Q，M、N、P、Q分别为大于0的自然数。

然后确定待匹配图像I₂中像素位置(x,y)处的高斯核为G(x,y,t)，并计算得到待匹配图像I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)，(x,y)∈I₂；进而得到待匹配图像I₂的尺度空间图像L(t)，t表示尺度参数，t>0。

具体地，分别获取基准图像I₁和待匹配图像I₂，所述基准图像I₁和待匹配图像I₂分别为灰度图像，且基准图像I₁的尺寸为M×N，待匹配图像I₂的尺寸为P×Q，M＝P，N＝Q，M、N、P、Q分别为大于0的自然数。

将基准图像I₁中像素位置为(x',y')的像素值记为I₁(x',y')，且该基准图像I₁中像素位置为(x',y')的像素值I₁(x',y')对应待匹配图像I₂中像素位置为(x,y)的像素值I₂(x,y)，其对应关系为：

I₂(x,y)＝I₁(α(x'cosφ₀+y'sinφ₀)+x₀,α(-x'sinφ₀+y'cosφ₀)+y₀)

其中，(x,y)∈I₂，(x',y')∈I₁，x表示待匹配图像I₂的像素位置在水平方向的坐标，y表示待匹配图像I₂的像素位置在垂直方向的坐标，x'表示基准图像I₁的像素位置在水平方向的坐标，y'表示基准图像I₁的像素位置在垂直方向的坐标，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量。

确定待匹配图像I₂中像素位置(x,y)处的高斯核为G(x,y,t)，

并计算得到待匹配图像I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)，其表达式为：

其中，*表示卷积操作，(x,y)∈I₂；进而得到待匹配图像I₂的尺度空间图像L(t)，t表示尺度参数，t>0，t＝σ²，σ表示待匹配图像I₂中像素位置(x,y)处的高斯核G(x,y,t)的标准差；所述待匹配图像I₂的尺度空间图像L(t)包含若干个尺度空间图像，每一个尺度参数值对应一个尺度空间图像；本实施例中尺度参数t的取值分别为1.0、2.0和4.0。

随着待匹配图像I₂中每个像素位置处的高斯核的标准差的增大，待匹配图像I₂的尺度空间图像L(t)中的尺度参数t也随之逐渐增大，进而使得待匹配图像I₂的尺度空间图像L(t)中小于待匹配图像I₂中每个像素位置处的高斯核的标准差σ的对应尺度空间图像的纹理信息将被抑制，对应的每一个尺度空间图像的模糊程度也随之增大，最终模拟出人在距离目标由近到远时目标在人眼中的成像。

步骤2，对待匹配图像I₂的尺度空间图像L(t)分别进行双线性插值，得到双线性插值后的模拟成像S(t)。

具体地，将待匹配图像I₂的尺度空间图像L(t)的尺寸设定为X×Y，X＝P，Y＝Q，X、Y分别为大于0的自然数；则双线性插值后的模拟成像S(t)的尺寸为(X/s)×(Y/s)，s表示由待匹配图像I₂的尺度空间图像L(t)到双线性插值后的模拟成像S(t)的缩放比率，t表示尺度参数，t>0，并将双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)处的灰度值记为S(m,n,t)，其表达式为：

x″＝m×s；y″＝n×s；

x₂＝x₁+1；

y₂＝y₁+1；

其中，(x″,y″)表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置，x″表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置在水平方向的坐标，y″表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置在垂直方向的坐标，(x″,y″)∈I₂；与双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应的待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置(x″,y″)最接近的4个像素位置分别为(x₁,y₁)、(x₂,y₁)、(x₁,y₂)和(x₂,y₂)，L(x₁,y₁,t)表示待匹配图像I₂中像素位置(x₁,y₁)处的尺度空间图像像素值，L(x₂,y₁,t)表示待匹配图像I₂中像素位置(x₂,y₁)处的尺度空间图像像素值，L(x₁,y₂,t)表示待匹配图像I₂中像素位置(x₁,y₂)处的尺度空间图像像素值，L(x₂,y₂,t)表示待匹配图像I₂中像素位置(x₂,y₂)处的尺度空间图像像素值，表示向下取整操作，t表示尺度参数，t>0。

步骤3，设定与待匹配图像I₂的尺度空间图像L(t)尺寸相同的空白图像，并将双线性插值后的模拟成像S(t)中每一个像素位置处的灰度值分别嵌入至所述空白图像中，得到多尺度图像S'(t)，t表示尺度参数，t>0；所述多尺度图像S'(t)包含若干个尺度图像，且所述多尺度图像S'(t)的尺寸为X×Y；X、Y分别为大于0的自然数。

基准图像I₁中像素位置(x',y')的像素值I₁(x',y')对应多尺度图像S'(t)中像素位置处的尺度图像其对应关系为：

其中，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0；本实施例中尺度参数t的取值分别为1.0、2.0和4.0。

步骤4，对基准图像I₁和多尺度图像S'(t)分别进行快速傅里叶变换，分别得到快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)，进而分别计算得到基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)，t表示尺度参数，t>0。

步骤4的子步骤为：

4.1对基准图像I₁和多尺度图像S'(t)分别进行快速傅里叶变换，分别得到快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)；将快速傅里叶变换后的基准图像频谱F₁中频率位置为的频谱记为且该快速傅里叶变换后的基准图像频谱F₁中频率位置为的频谱对应快速傅里叶变换后的多尺度图像频谱F₂(t)中频率位置(u,v)处的频谱F₂(u,v,t)，其对应关系为：

其中， 表示快速傅里叶变换后的基准图像频谱F₁的频率位置在水平方向的坐标，表示快速傅里叶变换后的基准图像频谱F₁的频率位置在垂直方向的坐标；(u,v)∈F₂(t)，u表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在水平方向的坐标，v表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在垂直方向的坐标；α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，·表示点乘，t表示尺度参数，t>0。

4.2分别提取快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)的幅度部分，得到快速傅里叶变换后基准图像频谱的幅度图像|F₁|和快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|。

将快速傅里叶变换后基准图像频谱的幅度图像|F₁|中频率位置为的频谱幅度记为且该快速傅里叶变换后基准图像频谱的幅度图像|F₁|中频率位置为的频谱幅度对应快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|中频率位置(u,v)处的频谱幅度|F₂(u,v,t)|，其对应关系为：

4.3对快速傅里叶变换后基准图像频谱的幅度图像|F₁|和快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|分别进行对数极坐标变换，即分别令对数坐标系中的距离参数为ρ，令对数极坐标系中的角度参数为θ，且分别满足u＝e^ρcosθ，v＝e^ρsinθ，分别得到快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像和快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像进而得到快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像中频率位置为的频谱幅度对应快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像中频率位置(u,v)处的频谱幅度其对应关系为：

其中，t表示尺度参数，t>0；α表示基准图像I₁中像素位置(x',y')缩放的尺度。

分别将快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像作为基准图像I₁的傅里叶梅林不变量图像F_1LP；将快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像作为多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)，t表示尺度参数，t>0。

将多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)中对数坐标系中的距离参数ρ、对数极坐标系中的角度参数θ处的像素位置记为S'_LP(ρ,θ,t)，该多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)中距离参数ρ、对数极坐标系中的角度参数θ处的像素位置S'_LP(ρ,θ,t)对应基准图像I₁的傅里叶梅林不变量图像F_1LP中距离参数ρ'、角度参数θ'处的像素位置F_1LP(ρ',θ')，θ'＝θ-φ₀。

其中，ln表示以e为底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0；本实施例中尺度参数t的取值分别为1.0、2.0和4.0。

步骤5，根据基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)，计算得到基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)。

步骤5的子步骤为：

5.1对基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)分别进行快速傅里叶变换，分别得快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP和快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP。

快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像尺度图像频谱与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)，满足下述关系：

其中，ln表示以e为底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0。

5.2计算得到快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像频谱与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)的互功率谱C(u,v,t)，其表达式为：

其中，*表示求复共轭操作，| |表示取幅值操作。

进而计算得到快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像频谱Γ_1LP(u,v)，与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)的相位相关互功率谱其表达式为：

5.3对所述相位相关互功率谱进行快速傅里叶反变换，得到基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)，t表示尺度参数，t>0，所述相位相关匹配图像I_MATCH(t)中包含若干个相位相关匹配图像，每一个尺度参数t对应一个相位相关匹配图像。

将相位相关匹配图像I_MATCH(t)中像素位置为(i,j)处的相位相关匹配图像记为I_MATCH(i,j,t)，其表达式为：

其中，δ()表示冲击函数，ln表示以为e底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，t表示尺度参数，t>0。

分别计算相位相关匹配图像I_MATCH(t)中最大相位相关匹配图像幅度值对应的像素位置在水平方向上的坐标i_max和相位相关匹配图像I_MATCH(t)中最大相位相关匹配图像幅度值对应的像素位置在垂直方向上的坐标j_max，其表达式分别为：

j_max＝φ₀；ln表示以为e底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，t表示尺度参数，t>0；本实施例中尺度参数t的取值分别为1.0、2.0和4.0。

步骤6，根据基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)，分别计算得到基准图像I₁中像素位置缩放的尺度估计值和基准图像I₁中像素位置旋转的角度估计值其中所述基准图像I₁中像素位置对应待匹配图像I₂中像素位置(x,y)；将所述基准图像I₁中每个像素分别进行尺度估计值的缩放，同时将所述基准图像I₁中每个像素位置进行角度估计值的旋转后，得到待匹配图像I₂的最优匹配图像。

具体地，对相位相关匹配图像I_MATCH(t)中每一个尺度参数t对应的相位相关匹配图像分别进行灰度最大值检测，分别得到相位相关匹配图像的最大值I_Max(t)，以及相位相关匹配图像的最大值I_Max(t)在水平方向的坐标和相位相关匹配图像的最大值I_Max(t)在垂直方向的坐标然后分别获取相位相关匹配图像的最大值I_Max(t)中的极大值I_Best、极大值I_Best对应的尺度参数极大值t_best、尺度参数极大值t_best对应的水平方向坐标和尺度参数极大值t_best对应的垂直方向坐标然后分别计算基准图像I₁中像素位置(x',y')缩放的尺度估计值和基准图像I₁中像素位置(x',y')旋转的角度估计值其表达式分别为：

将所述基准图像I₁中每个像素分别进行尺度估计值的缩放，同时将所述基准图像I₁中每个像素位置进行角度估计值的旋转后，得到待匹配图像I₂的最优匹配图像。

通过下述仿真实验对本发明效果作进一步验证说明。

图2为本发明实施例提供的基准图像示意图；其中，该基准图像为经典的Lena图，尺寸为128×128；图3为本发明实施例提供的待匹配图像示意图，是对图1中每个像素分别进行尺度为2.4的缩放、角度为-45°的旋转变换后切取的部分图像，尺寸为128×128；根据本发明方法对图2分别进行尺度和旋转参数估计后，分别得到图2中每个像素位置缩放的尺度估计值为2.392，图2中每个像素位置旋转的角度估计值为-45°；而根据原始的傅里叶梅林变换的图像匹配方法进行参数估计后，得到图2中每个像素位置缩放的尺度估计值为1.106，图2中每个像素位置旋转的角度估计值为4.2°；其中，图4为本发明实施例提供的傅里叶梅林不变量生成图的源图像，图5为图4的傅里叶变换幅度谱示例图，图6为图4的傅里叶梅林不变量示意图。

结果表明，本发明方法在鲁棒性上明显优于原始的基于傅里叶梅林变换的图像匹配方法。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (10)

1.一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，包括以下步骤：

步骤1，分别获取基准图像I₁和待匹配图像I₂，所述基准图像I₁和待匹配图像I₂分别为灰度图像，且基准图像I₁的尺寸为M×N，待匹配图像I₂的尺寸为P×Q，M＝P，N＝Q，M、N、P、Q分别为大于0的自然数；

然后确定待匹配图像I₂中像素位置(x,y)处的高斯核为G(x,y,t)，并计算得到待匹配图像I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)，(x,y)∈I₂；进而得到待匹配图像I₂的尺度空间图像L(t)，t表示尺度参数，t>0；

步骤2，对待匹配图像I₂的尺度空间图像L(t)分别进行双线性插值，得到双线性插值后的模拟成像S(t)；

步骤3，设定与待匹配图像I₂的尺度空间图像L(t)尺寸相同的空白图像，并将双线性插值后的模拟成像S(t)中每一个像素位置处的灰度值分别嵌入至所述空白图像中，得到多尺度图像S'(t)；所述多尺度图像S'(t)包含若干个尺度图像，且所述多尺度图像S'(t)的尺寸为X×Y；X、Y分别为大于0的自然数；

步骤4，对基准图像I₁和多尺度图像S'(t)分别进行快速傅里叶变换，分别得到快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)，进而分别计算得到基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)；

步骤5，根据基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)，计算得到基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)；

步骤6，根据基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)，分别计算得到基准图像I₁中像素位置缩放的尺度估计值和基准图像I₁中像素位置旋转的角度估计值其中所述基准图像I₁中像素位置对应待匹配图像I₂中像素位置(x,y)；将所述基准图像I₁中每个像素分别进行尺度估计值的缩放，同时将所述基准图像I₁中每个像素位置进行角度估计值的旋转后，得到待匹配图像I₂的最优匹配图像。

2.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，在步骤1中，所述得到待匹配图像I₂的尺度空间图像L(t)，其过程为：

将基准图像I₁中像素位置为(x',y')的像素值记为I₁(x',y')，且该基准图像I₁中像素位置为(x',y')的像素值I₁(x',y')对应待匹配图像I₂中像素位置为(x,y)的像素值I₂(x,y)，其对应关系为：

I₂(x,y)＝I₁(α(x'cosφ₀+y'sinφ₀)+x₀,α(-x'sinφ₀+y'cosφ₀)+y₀)

其中，(x,y)∈I₂，(x',y')∈I₁，x表示待匹配图像I₂的像素位置在水平方向的坐标，y表示待匹配图像I₂的像素位置在垂直方向的坐标，x'表示基准图像I₁的像素位置在水平方向的坐标，y'表示基准图像I₁的像素位置在垂直方向的坐标，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量；

确定待匹配图像I₂中像素位置(x,y)处的高斯核为G(x,y,t)，

并计算得到待匹配图像I₂中像素位置(x,y)处的尺度空间图像像素值L(x,y,t)，其表达式为：

$L(x,y,t)=\left\{\begin{array}{cc}G(x,y,t)*I_2(x,y)& t>0\\ I_2(x,y)& t=0\end{array}\right.$

其中，*表示卷积操作，(x,y)∈I₂；进而得到待匹配图像I₂的尺度空间图像L(t)，t表示尺度参数，t>0，t＝σ²，σ表示待匹配图像I₂中像素位置(x,y)处的高斯核G(x,y,t)的标准差；所述待匹配图像I₂的尺度空间图像L(t)包含若干个尺度空间图像，每一个尺度参数值对应一个尺度空间图像。

3.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，在步骤2中，所述双线性插值后的模拟成像S(t)，还包括：

所述双线性插值后的模拟成像S(t)的尺寸为(X/s)×(Y/s)，待匹配图像I₂的尺度空间图像L(t)的尺寸为X×Y，X＝P，Y＝Q，X、Y分别为大于0的自然数；s表示由待匹配图像I₂的尺度空间图像L(t)到双线性插值后的模拟成像S(t)的缩放比率，t表示尺度参数，t>0；并将双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)处的灰度值记为S(m,n,t)，其表达式为：

$\begin{array}{c}S(m,n,t)=\frac{(x_2-x^{\′\′})(y_2-y^{\′\′})}{(x_2-x_1)(y_2-y_1)}L(x_1,y_1,t)+\frac{(x^{\′\′}-x_1)(y_2-y^{\′\′})}{(x_2-x_1)(y_2-y_1)}L(x_2,y_1,t)\\ +\frac{(x_2-x^{\′\′})(y^{\′\′}-y_1)}{(x_2-x_1)(y_2-y_1)}L(x_1,y_2,t)+\frac{(x^{\′\′}-x_1)(y^{\′\′}-y_1)}{(x_2-x_1)(y_2-y_1)}L(x_2,y_2,t)\end{array}$

$s=\sqrt{t};$

x”＝m×s；y”＝n×s；

其中，(x”,y”)表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置，x”表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置在水平方向的坐标，y”表示双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置在垂直方向的坐标，(x”,y”)∈I₂；与双线性插值后的模拟成像S(t)中任意一个像素位置(m,n)对应的待匹配图像I₂的尺度空间图像L(t)中尺度空间图像位置(x”,y”)最接近的4个像素位置分别为(x₁,y₁)、(x₂,y₁)、(x₁,y₂)和(x₂,y₂)，

L(x₁,y₁,t)表示待匹配图像I₂中像素位置(x₁,y₁)处的尺度空间图像像素值，L(x₂,y₁,t)表示待匹配图像I₂中像素位置(x₂,y₁)处的尺度空间图像像素值，L(x₁,y₂,t)表示待匹配图像I₂中像素位置(x₁,y₂)处的尺度空间图像像素值，L(x₂,y₂,t)表示待匹配图像I₂中像素位置(x₂,y₂)处的尺度空间图像像素值，表示向下取整操作。

4.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，在步骤3中，所述多尺度图像S'(t)，还包括：

基准图像I₁中像素位置(x',y')的像素值I₁(x',y')对应多尺度图像S'(t)中像素位置处的尺度图像其对应关系为：

$S^{\′}(\hat{x},\hat{y},t)=I_1\left(\frac{\mathrm{\α}(x^{\′}{\mathrm{cos\φ}}_0+y^{\′}{\mathrm{sin\φ}}_0)+x_0}{\sqrt{t}},\frac{\mathrm{\α}(-x^{\′}{\mathrm{sin\φ}}_0+y^{\′}{\mathrm{cos\φ}}_0)+y_0}{\sqrt{t}}\right)$

其中，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0。

5.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，步骤4的子步骤为：

4.1对基准图像I₁和多尺度图像S'(t)分别进行快速傅里叶变换，分别得到快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)；

4.2分别提取快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)的幅度部分，得到快速傅里叶变换后基准图像频谱的幅度图像|F₁|和快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|，进而分别得到快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像和快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像

4.3分别将快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像作为基准图像I₁的傅里叶梅林不变量图像F_1LP；将快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像作为多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)；

其中，将多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)中对数坐标系中的距离参数ρ、对数极坐标系中的角度参数θ处的像素位置记为S'_LP(ρ,θ,t)，该多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)中距离参数ρ、对数极坐标系中的角度参数θ处的像素位置S'_LP(ρ,θ,t)对应基准图像I₁的傅里叶梅林不变量图像F_1LP中距离参数ρ'、角度参数θ'处的像素位置F_1LP(ρ',θ')，θ'＝θ-φ₀；ln表示以e为底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度；φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0。

6.如权利要求5所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，所述快速傅里叶变换后的基准图像频谱F₁和快速傅里叶变换后的多尺度图像频谱F₂(t)，还包括：

将快速傅里叶变换后的基准图像频谱F₁中频率位置为的频谱记为且该快速傅里叶变换后的基准图像频谱F₁中频率位置为的频谱对应快速傅里叶变换后的多尺度图像频谱F₂(t)中频率位置(u,v)处的频谱F₂(u,v,t)，其对应关系为：

$F_2(u,v,t)=\frac{t}{{\mathrm{\α}}^2}{e}^{-j\frac{2\mathrm{\π}}{\sqrt{t}}(x_{0}u+y_{0}v)}\·F_1(\stackrel{\‾}{x},\stackrel{\‾}{y})$

$\stackrel{\‾}{x}=\frac{\sqrt{t}}{\mathrm{\α}}(u{\mathrm{cos\φ}}_0+v{\mathrm{sin\φ}}_0),\stackrel{\‾}{y}=\frac{\sqrt{t}}{\mathrm{\α}}(-u{\mathrm{sin\φ}}_0+v{\mathrm{cos\φ}}_0)$

其中， 表示快速傅里叶变换后的基准图像频谱F₁的频率位置在水平方向的坐标，表示快速傅里叶变换后的基准图像频谱F₁的频率位置在垂直方向的坐标；(u,v)∈F₂(t)，u表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在水平方向的坐标，v表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在垂直方向的坐标；α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，·表示点乘，t表示尺度参数，t>0。

7.如权利要求5所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，所述快速傅里叶变换后基准图像频谱的幅度图像|F₁|和快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|，还包括：

将快速傅里叶变换后基准图像频谱的幅度图像|F₁|中频率位置为的频谱幅度记为且该快速傅里叶变换后基准图像频谱的幅度图像|F₁|中频率位置为的频谱幅度对应快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|中频率位置(u,v)处的频谱幅度|F₂(u,v,t)|，其对应关系为：

$\left|F_2(u,v,t)\right|=\frac{t}{{\mathrm{\α}}^2}\left|F_1(\stackrel{\‾}{x},\stackrel{\‾}{y})\right|,\stackrel{\‾}{x}=\frac{\sqrt{t}}{\mathrm{\α}}(u{\mathrm{cos\φ}}_0+v{\mathrm{sin\φ}}_0),\stackrel{\‾}{y}=\frac{\sqrt{t}}{\mathrm{\α}}(-u{\mathrm{sin\φ}}_0+v{\mathrm{cos\φ}}_0)$

其中， 表示快速傅里叶变换后的基准图像频谱F₁的频率位置在水平方向的坐标，表示快速傅里叶变换后的基准图像频谱F₁的频率位置在垂直方向的坐标；(u,v)∈F₂(t)，u表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在水平方向的坐标，v表示快速傅里叶变换后的多尺度图像频谱F₂(t)的频率位置在垂直方向的坐标；α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，·表示点乘，t表示尺度参数，t>0。

8.如权利要求5所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，所述对快速傅里叶变换后基准图像频谱的幅度图像|F₁|和快速傅里叶变换后多尺度图像频谱的幅度图像|F₂(t)|分别进行对数极坐标变换，具体为：

分别令对数坐标系中的距离参数为ρ，令对数极坐标系中的角度参数为θ，且分别满足u＝e^ρcosθ，v＝e^ρsinθ，分别得到快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像和快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像进而得到快速傅里叶变换后基准图像频谱的幅度图像在对数极坐标下的幅度图像中频率位置为的频谱幅度对应快速傅里叶变换后多尺度图像频谱的幅度图像在对数极坐标下的图像中频率位置(u,v)处的频谱幅度其对应关系为：

$\left|{\tilde{F}}_2(e^{\mathrm{\ρ}}cos\mathrm{\θ},e^{\mathrm{\ρ}}sin\mathrm{\θ},t)\right|=\frac{t}{{\mathrm{\α}}^2}\left|{\tilde{F}}_1\left(e^{\mathrm{\ρ}-\ln \frac{\mathrm{\α}}{\sqrt{t}}}cos\right(\mathrm{\θ}-{\mathrm{\φ}}_0),e^{\mathrm{\ρ}-\ln \frac{\mathrm{\α}}{\sqrt{t}}}sin(\mathrm{\θ}-{\mathrm{\φ}}_0\left)\right)\right|$

其中，t表示尺度参数，t>0；α表示基准图像I₁中像素位置(x',y')缩放的尺度。

9.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，步骤5的子步骤为：

5.1对基准图像I₁的傅里叶梅林不变量图像F_1LP和多尺度图像S'(t)的傅里叶梅林不变量图像S'_LP(t)分别进行快速傅里叶变换，分别得快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP和快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP；

快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像尺度图像频谱与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)，满足下述关系：

${{\mathrm{\Γ}}^{\′}}_{LP}(u,v,t)=e^{-j2\mathrm{\π}(uln\frac{\mathrm{\α}}{\sqrt{t}}+{\mathrm{v\φ}}_0)}\×\frac{t}{{\mathrm{\α}}^2}{\mathrm{\Γ}}_{1LP}(\stackrel{\‾}{x},\stackrel{\‾}{y})$

其中，ln表示以e为底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，x₀表示基准图像I₁中像素位置(x',y')在水平方向的平移量，y₀表示基准图像I₁中像素位置(x',y')在垂直方向的平移量，t表示尺度参数，t>0；

5.2计算得到快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像频谱与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)的互功率谱C(u,v,t)，其表达式为：

$C(u,v,t)=\frac{{\mathrm{\Γ}}_{1LP}(\stackrel{\‾}{x},\stackrel{\‾}{y}){{\mathrm{\Γ}}^{\′}}_{LP}^{*}(u,v,t)}{\left|{\mathrm{\Γ}}_{1LP}(\stackrel{\‾}{x},\stackrel{\‾}{y}){{\mathrm{\Γ}}^{\′}}_{LP}(u,v,t)\right|}$

其中，*表示求复共轭操作，| |表示取幅值操作；

进而计算得到快速傅里叶变换后基准图像的傅里叶梅林不变量图像频谱Γ_1LP中频率位置为处的傅里叶梅林不变量图像频谱Γ_1LP(u,v)，与快速傅里叶变换后多尺度图像的傅里叶梅林不变量图像频谱Γ'_LP中频率位置(u,v)处的傅里叶梅林不变量图像频谱Γ'_LP(u,v,t)的相位相关互功率谱其表达式为：

$\hat{C}(u,v,t)=e^{j2\mathrm{\π}(u\ln \frac{\mathrm{\α}}{\sqrt{t}}+{\mathrm{v\φ}}_0)};$

5.3对所述相位相关互功率谱进行快速傅里叶反变换，得到基准图像I₁与待匹配图像I₂的尺度空间图像L(t)之间的相位相关匹配图像I_MATCH(t)，t表示尺度参数，t>0；

所述相位相关匹配图像I_MATCH(t)中包含若干个相位相关匹配图像，每一个尺度参数t对应一个相位相关匹配图像；

将相位相关匹配图像I_MATCH(t)中像素位置为(i,j)处的相位相关匹配图像记为I_MATCH(i,j,t)，其表达式为：

$I_{MATCH}(i,j,t)=\mathrm{\δ}(i-ln\frac{\mathrm{\α}}{\sqrt{t}},j-{\mathrm{\φ}}_0)$

其中，δ()表示冲击函数，ln表示以为e底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度，t表示尺度参数，t>0；

分别计算相位相关匹配图像I_MATCH(t)中最大相位相关匹配图像幅度值对应的像素位置在水平方向上的坐标i_max和相位相关匹配图像I_MATCH(t)中最大相位相关匹配图像幅度值对应的像素位置在垂直方向上的坐标j_max，其表达式分别为：

j_max＝φ₀；ln表示以为e底的对数操作，α表示基准图像I₁中像素位置(x',y')缩放的尺度，φ₀表示基准图像I₁中像素位置(x',y')旋转的角度。

10.如权利要求1所述的一种基于多尺度傅里叶梅林变换的图像匹配方法，其特征在于，在步骤6中，所述得到基准图像I₁中像素位置缩放的尺度估计值和基准图像I₁中像素位置旋转的角度估计值其过程为：

对相位相关匹配图像I_MATCH(t)中每一个尺度参数t对应的相位相关匹配图像分别进行灰度最大值检测，分别得到相位相关匹配图像的最大值I_Max(t)，以及相位相关匹配图像的最大值I_Max(t)在水平方向的坐标和相位相关匹配图像的最大值I_Max(t)在垂直方向的坐标然后分别获取相位相关匹配图像的最大值I_Max(t)中的极大值I_Best、极大值I_Best对应的尺度参数极大值t_best、尺度参数极大值t_best对应的水平方向坐标和尺度参数极大值t_best对应的垂直方向坐标然后分别计算基准图像I₁中像素位置缩放的尺度估计值和基准图像I₁中像素位置旋转的角度估计值其表达式分别为：

$\widehat{\mathrm{\α}}=\sqrt{t_{best}}{e}^{{\tilde{i}}_{max}\left(t_{best}\right)},{\widehat{\mathrm{\φ}}}_0={\tilde{j}}_{max}\left(t_{best}\right).$