CN106778899A - 一种基于统计相关的快速互信息图像匹配方法 - Google Patents

Abstract

本发明公开了一种基于统计相关的快速互信息图像匹配方法，先用传统方法计算实测图和第一个位置基准子图各自的直方图、熵以及两图的联合直方图、联合熵、互信息并存储；接着在遍历搜索过程中，利用相邻基准子图之间的相关性，以前一位置基准子图的直方图、熵以及与实测图的联合熵为基准，通过差量法计算当前位置基准子图的直方图、熵以及与实测图的联合熵，并计算互信息；最后，在所有搜索位置中找出互信息值最大的位置作为最终的匹配定位结果；本发明除基准图第一位置基准子图外，其它位置子图的直方图、熵及联合熵均基于差量法依次得到，在保证匹配精度的同时提高了互信息匹配的速度；此外，还可与已有方法相结合，通过多种途径来加快匹配速度。

Description

一种基于统计相关的快速互信息图像匹配方法

技术领域

本发明属于图像匹配、图像定位技术领域，特别涉及一种基于统计相关的快速互信息图像匹配方法。

背景技术

图像匹配技术是在航天技术、卫星应用技术、传感器技术、计算机技术、图像处理及模式识别的基础上发展起来的一门新技术,它在飞机辅助导航、远程武器或精确制导武器系统如巡航导弹的末制导、图像目标的搜索与跟踪等军事领域具有重要的应用价值。目前，匹配技术在工业控制、医学图像处理、遥感图像处理、制图学、目标识别等领域有着广泛应用。在基于视觉的飞行器导航与制导、遥感卫星灾害监控与环境监测、医学图像分析等应用中，常常需要对不同成像传感器获取的异源图像进行匹配。异源图像最大的特点是异源图像同名点之间的灰度值一般不具有直接联系。因此，同源图像的匹配方法多数难以直接应用于异源图像匹配。基于互信息测度的匹配方法是一种完全基于图像灰度统计概率的匹配方法，不需要对原图像间的灰度关系作任何假设,比传统的相似性度量具有更强的鲁棒性，但互信息准则最大的缺陷是计算量大，匹配耗时较长。为此，研究者们就围绕减少非匹配点相关计算量和改进搜索策略减少计算量两个方面提出了相应的加速方法。例如文献安如，王慧麟等，16阶归一化互信息和改进PSO算法的快速图像匹配，吉林大学学报(工学版)，2013年3月第43卷357-364，采用灰度压缩和粒子群优化算法相结合加快匹配速度。文献杨猛,潘泉等.基于定量定性互信息的多层次特征图像匹配算法，中国图象图形学报2010年9月第15卷第9期1376-1383，通过提取多层次特征点减少互信息测度的计算量。文献A Dame，E Marchand.Second-order optimization of mutual information for real-timeimage registration.IEEE Transactions on Image Processing A Publication of theIEEE Signal Processing Society 2012,21(9):4190-203，利用粒子群优化算法提高搜索速度。文献Y Zhuang，K Gao，X Miu，L Han，X Gong.Infrared and visual imageregistration based on mutual information with a combined particle swarmoptimization–Powell search algorithm.Optik-International Journal for Lightand Electron Optics.将粒子群优化方法与Powell优化方法相结合提高匹配速度。上述方法中，灰度压缩方法和特征提取方法会损失图像部分灰度信息，影响匹配精度，而粒子群优化方法又容易陷入局部极值，影响匹配精度，Powell方法对初始点的要求较高，初始点的选择会直接影响到最后的搜索最优值。总之，传统的快速互信息匹配算法大多是通过降灰度阶或减少参与互信息测度计算的像素数来加速匹配速度，一定程度上会对匹配精度造成影响。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供一种基于统计相关的快速互信息图像匹配方法，与已有的以灰度压缩、特征提取或采取优化搜索策略等快速互信息匹配方法不同，本发明利用匹配过程中相邻基准子图间的相关性，通过差量法减少每一个匹配位置互信息的计算量来加快匹配速度，同时，本发明方法可以与灰度压缩、特征提取等方法相结合，通过多种途径进一步加快匹配速度。该方法能有效地提高互信息匹配的实时性，可用于对实时性要求比较高的异源图像匹配或图像配准应用系统中。

为了实现上述目的，本发明采用的技术方案是：

一种基于统计相关的快速互信息图像匹配方法，包括如下步骤：

步骤1.获取实测图像A和基准图像I，将两幅图像的灰度值调整到同一灰度区间；

步骤2.统计实测图A各灰度值出现次数并存入灰度直方图矩阵h_A，计算A的信息熵H_A并存储，H_A的计算公式为：

其中，T为实测图总像素数。

步骤3.从基准图I的左上角(1,1)点截取与实测图大小相等的第1幅基准子图S_ini，统计S_ini各灰度值出现次数并存入灰度直方图矩阵h_{S_ini}，计算S_ini的信息熵H_{S_ini}并存储，H_{S_ini}的计算公式为：

其中，T为基准子图总像素数。

步骤4.匹配过程中，通过比较当前基准子图和前一基准子图对应行或列像素变化情况，按照差量法将前一基准子图的直方图矩阵h_S0更新为当前子图的直方图矩阵h_S。

步骤5.找出当前基准子图直方图矩阵h_S与前一基准子图直方图矩阵h_S0的不同元素，计算这些元素对应的熵值求和项，通过差量法将前一基准子图的信息熵H_S0更新为当前基准子图的信息熵H_S。

步骤6.找出当前基准子图联合直方图矩阵h_AS和前一基准子图联合直方图矩阵h_AS0的不同元素，计算这些元素对应的熵值求和项，通过差量法将前一基准子图的联合熵H_AS0更新为当前基准子图的联合熵H_AS。

步骤7.计算实测图与当前基准子图的归一化互信息NMI，计算公式为：

NMI＝(H_A+H_S)/H_AS

步骤8.遍历搜索整个基准图，以互信息最大值对应的位置作为最终匹配位置。

步骤4的“按照差量法将前一基准子图的直方图矩阵h_S0更新为当前子图的直方图矩阵h_S”，按如下步骤进行：

步骤4-1在列方向，通过统计当前基准子图最后一行和前一基准子图第一行对应像素灰度值出现次数，以前一基准子图的直方图矩阵h_S0为基准，减去前一基准子图第一行对应像素灰度值出现次数，加上当前基准子图最后一行对应像素灰度值出现次数即得到当前子图的直方图矩阵h_S；

步骤4-2在行方向，通过统计当前基准子图最后一列和前一基准子图第一列对应像素灰度值出现次数，以前一基准子图的直方图矩阵h_S0为基准，减去前一基准子图第一列对应像素灰度值出现次数，加上当前基准子图最后一列对应像素灰度值出现次数即得到当前子图的直方图矩阵h_S；

步骤5的“通过差量法将前一基准子图的联合熵H_S0更新为当前基准子图的联合熵H_S”，按如下步骤进行：

步骤5-1找出h_S与h_S0的不同元素，分别组成矩阵h_ΔS和h_ΔS0，计算h_ΔS和h_ΔS0对应的信息熵Δ_S及Δ_S0；

Δ_S计算公式为：

Δ_S0计算公式为：

其中，T为基准子图总像素数，norm_s0(s)为与h_ΔS0中元素对应的前一基准子图熵值计算过程中的求和项，此次只计算这些求和项的和，无需对各个求和项再进行重复计算。

步骤5-2以前一基准子图的信息熵H_S0为基准，通过差量法将H_S0更新为当前基准子图的信息熵H_S；

H_S对应的计算公式为：

H_S＝H_S0-Δ_S0+Δ_S

步骤6的“通过差量法将前一基准子图的联合熵H_AS0更新为当前基准子图的联合熵H_AS”，按如下步骤进行：

步骤6-1找出h_AS与h_AS0的不同元素，分别组成矩阵h_ΔAS和h_ΔAS0，计算h_ΔAS和h_ΔAS0对应的信息熵Δ_AS及Δ_AS0；

Δ_AS计算公式为：

Δ_AS0计算公式为：

其中，T为基准子图总像素数，norm_AS(s)为与h_ΔAS0中元素对应的前一基准子图联合熵计算过程中的求和项，此次只计算这些求和项的和，无需对各个求和项再进行重复计算。

步骤6-2以前一基准子图的信息熵H_AS0为基准，通过差量法将H_AS0更新为当前基准子图的信息熵H_AS；

H_AS对应的计算公式为：

H_AS＝H_AS0-Δ_AS0+Δ_AS

步骤7.计算实测图与当前基准子图的归一化互信息NMI，计算公式为：

NMI＝(H_A+H_S)/H_AS

步骤8.遍历整个基准图，以互信息最大值对应的位置作为最终匹配点。

与现有技术相比，本发明的有益效果是：

1.本发明利用了匹配过程中相邻基准子图之间具有相关性的特点，在进行互信息测度计算时，只用传统方法计算实测图和第1幅基准子图的灰度直方图矩阵、熵、联合熵以及实测图与各基准子图的联合直方图矩阵，遍历搜索过程中其它所有基准子图的灰度直方图、熵以及联合熵都采用差量法完成计算。因为互信息测度的计算量主要集中在熵值计算过程中各个求和项的计算，采用差量法只需更新熵值计算过程中发生变化的求和项，与传统的独立计算每个匹配位置的互信息匹配方法相比，计算量大为减少。设实测图大小为m*n，本发明方法基准子图直方图的计算量为传统方法的或多次实验表明，每个基准子图熵计算量与传统方法熵计算量相比平均减少约50％，联合熵计算量减少约25％；当图像灰度级别越多,像素灰度值越分散,图像熵值越大的情况本发明的提速效果越好。

2.本发明采用遍历搜索，能够避免一般优化搜索方法可能收敛到局部极值的缺陷，保证匹配结果的可靠性。

3.本发明方法能够与已有的灰度压缩、特征提取等方法相结合，进一步提高匹配的实时性。

附图说明

图1是本发明应用系统主要组件。

图2是本发明互信息匹配实施过程的流程图。

图3是本发明应用于SAR与可见光匹配的仿真结果示意图。

图4是本发明应用于红外与可见光匹配的仿真结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图对本发明作进一步地详细描述。

参照附图1，本发明的应用系统主要包括5大组件：

组件1：匹配处理器CPU模块1，主要用来对实时拍摄的实测图和基准图进行匹配计算。

组件2：内部存储器模块2，主要用来暂时存储匹配过程中的原始数据、中间结果和最终匹配结果。

组件3：实时图像传感器模块3，主要用来实时地获取实测图像；

组件4：数字地图存储器模块4，主要用来存储预先拍摄好的基准图像；

组件5：基准图像选取模块5，主要根据外界其它导航信息在数字地图库中选取适用于本次匹配的基准图。

参照附图2，组件1采用基于互信息的快速匹配方法进行匹配，实现过程包括以下步骤：

步骤1.计算实测图A的灰度直方图矩阵h_A和信息熵H_A；

按照传统方法，通过遍历整个实测图统计实测图各像素灰度值的出现次数存入直方图矩阵h_A，设实测图总像素数为T，按照下面的求和公式用传统算法计算实测图熵H_A:

步骤2.计算基准图第一幅基准子图S_ini的灰度直方图矩阵h_{S_ini}和信息熵H_{S_ini}；

从基准图I的左上角(1,1)点截取第1幅基准子图S_ini，与实施步骤2的过程类似，用传统方法获得S_ini的灰度直方图矩阵h_{S_ini}，按照下面公式用传统算法计算S_ini的信息熵H_{S_ini}:

其中，norm_{_ini}为熵计算的求和项。

步骤3.从基准图I的(1,2)点开始，按照从左到右、从上至下的方法在基准图每一个点截取与实测图大小相同的基准子图，以前一基准子图的直方图、熵、联合熵为基准，按照差量法计算当前基准子图的直方图、熵及联合熵，并计算与实测图的互信息，遍历完整个基准图后，找出互信息最大值对应的位置的坐标作为最终匹配点。设当前基准子图左上角位置为(i,j)，下面结合图1和图2对本步骤作进一步地详细描述。

步骤3-1计算每行第一列处基准子图的直方图

设当前基准子图所在行号为i，该行第一列(i,1)点处基准子图的直方图为h_S，前一行即i-1行第一列(i-1,1)处基准子图直方图h_S0，则h_S是由h_S0通过差量法更新得到的。具体做法为：以(i-1,1)处基准子图的直方图则h_s0为基准，减去该子图第一行像素灰度值的出现次数，加上当前(i,1)处基准子图最后一行像素灰度出现次数即可得到当前基准子图直方图h_S。

步骤3-2计算每行第一列处基准子图的熵

设当前基准子图所在行号为i，该行第一列(i,1)点处基准子图的熵H_S，则H_S是由(i-1,1)处基准子图熵H_S0通过差量法更新得到的。具体做法为：在步骤4-1的基础上，找出h_S与h_S0的不同元素，分别组成矩阵h_ΔS和h_ΔS0，计算h_ΔS和h_ΔS0对应的信息熵Δ_S及Δ_S0；

Δ_S计算公式为：

Δ_S0计算公式为：

其中，T为基准子图总像素数,norm_s0(s)为与h_ΔS0中元素对应的前一基准子图熵值计算过程中的求和项，此次不需再另行计算。

再根据差量法将H_S0更新为H_S，对应的公式为：

H_S＝H_S0-Δ_S0+Δ_S

计算实测图与当前基准子图的归一化互信息NMI，计算公式为：

NMI＝(H_A+H_S)/H_AS

步骤3-3计算每行第一列处基准子图与实测图的联合熵

设当前基准子图所在行号为i，实测图与该行第一列(i,1)点处基准子图的联合熵H_AS是由实测图与(i-1,1)处基准子图的联合熵H_AS0通过差量法更新得到的。设实测图与(i-1,1)处基准子图的联合直方图为h_AS0，具体做法为：按照传统遍历法统计得到实测图与(i,1)处基准子图的联合直方图h_AS，找出h_AS与h_AS0的不同元素，分别组成矩阵h_ΔAS和h_ΔAS0，先计算h_ΔAS和h_ΔAS0对应的信息熵Δ_AS及Δ_AS0；

Δ_AS计算公式为：

Δ_AS0计算公式为：

其中，T为基准子图总像素数，norm_AS(s)为与h_ΔAS0中元素对应的前一基准子图联合熵计算过程中的求和项，此次不需再另行计算。

再根据差量法将H_AS0更新为H_AS，对应的公式为：

H_AS＝H_AS0-Δ_AS0+Δ_AS

计算实测图与当前基准子图的归一化互信息NMI，计算公式为：

NMI＝(H_A+H_S)/H_AS

步骤3-4在行方向，以每行第一列(i,1)处基准子图的直方图、熵、联合熵为基准，从左到右，依次通过差量法计算该行其它基准子图的直方图、熵及与实测图的联合熵，并计算与实测图的互信息。

设(i,j-1)处基准子图的直方图为h_S0、熵为H_S0、与实测图的联合熵为H_AS0；(i,j)处基准子图的直方图为h_S，熵为H_S、与实测图的联合熵为H_AS，具体做法为：以(i,j-1)处基准子图的直方图h_S0为基准，减去该子图第一列像素灰度值的出现次数，加上当前(i,j)处基准子图最后一列像素灰度出现次数即可得到当前基准子图直方图h_S；当前基准子图熵H_S、以及与实测图的联合熵H_AS的计算方法与步骤3-2及步骤3-3计算方法类似。

步骤3-5按上述步骤遍历完整个基准图后，找出互信息值最大的点作为最终匹配位置。

通过以上步骤，就可以实现基于互信息的快速匹配。

下面以两个实施例为例，通过比较每个匹配位置处本文方法与传统方法在互信息计算量方面的差别以及总耗时来说明本发明方法的有效性。具体实施过程为：以MatlabR2015为仿真平台，首先按照具体实施方案中的步骤3-1、步骤3-2以及步骤3-3分别统计列方向各基准子图的直方图、熵以及联合熵的计算量，并与传统计算方法进行比较；然后再按照具体实施方式中的步骤3-4分别统计行方向各基准子图的直方图、熵以及联合熵的计算量，并与传统计算方法进行比较；最后通过完整的匹配过程比较本发明方法与传统方法的匹配耗时。

匹配过程中设：基准子图直方图计算中像素灰度值统计次数为N，基准子图熵值计算中求和项的项数为N_H_S，基准子图与实测图联合熵计算中的求和项项数为N_H_AS。

实施例1：将本发明的遍历式互信息匹配方法和传统的遍历式互信息匹配方法分别应用于可见光与SAR图像的匹配，如图3所示，其中图3(a)为可见光基准图I，大小为422*358；图3(b)为SAR实测图像A1，大小为150*150，图3(c)为SAR实测图像A2，大小为100*90，图3(d)为实测图在基准图上的匹配定位结果；表1为本实施例列方向匹配时部分位置计算量统计结果，表2为本实施例行方向匹配时部分位置计算量统计结果。表3为本实施例的总匹配耗时。

表1实施例1在列方向匹配时部分位置计算量统计结果

表2实施例1在行方向匹配时部分位置计算量统计结果

表3实施例1的匹配耗时比较

实施例2：将本发明的遍历式互信息匹配方法和传统的遍历式互信息匹配方法分别应用于可见光与红外图像的匹配，如图4所示，其中图4(a)为可见光基准图I，大小为473*734；图4(b)为红外实测图像A1，大小为80*80，图4(c)为红外实测图像A2，大小为120*120，图4(d)为实测图在基准图上的匹配定位结果；表4为本实施例列方向匹配时部分位置计算量统计结果，表5为本实施例行方向匹配时部分位置计算量统计结果。表6为本实施例的匹配结果和匹配耗时。

表4实施例2在列方向匹配时部分位置计算量统计结果

表5实施例2在行方向匹配时部分位置计算量统计结果

表6实施例2的匹配耗时比较

Claims (4)

1.一种基于统计相关的快速互信息图像匹配方法，其特征在于，包括如下步骤：

步骤1.获取实测图像A和基准图像I，将两幅图像的灰度值调整到同一灰度区间；

步骤2.统计实测图像A各灰度值出现次数并存入灰度直方图矩阵h_A，计算A的信息熵H_A并存储，H_A的计算公式为：

$H_A=-\underset{a}{\Σ}\frac{h_A}{T}\left(a\right).log\frac{h_A}{T}\left(a\right)$

其中，T为实测图的像素总数，a表示像素值；

步骤3.从基准图像I的左上角(1,1)点截取与实测图像A大小相等的第1幅基准子图S_ini，统计S_ini各灰度值出现次数并存入灰度直方图矩阵h_{S_ini}，计算S_ini的信息熵H_{S_ini}并存储，H_{S_ini}的计算公式为：

$H_{S\_ini}=-\underset{s}{\Σ}\frac{h_{S\_ini}}{T}\left(s\right).log\frac{h_{S\_ini}}{T}\left(s\right)$

其中，T为基准子图的像素总数，s表示像素值；

步骤4.匹配过程中，通过比较当前基准子图和前一基准子图对应行或列像素变化情况，按照差量法将前一基准子图的直方图矩阵h_S0更新为当前子图的直方图矩阵h_S；

步骤5.找出当前基准子图直方图矩阵h_S与前一基准子图直方图矩阵h_S0的不同元素，计算这些元素对应的熵值求和项，通过差量法将前一基准子图的信息熵H_S0更新为当前基准子图的信息熵H_S；

步骤6.找出当前基准子图联合直方图矩阵h_AS和前一基准子图联合直方图矩阵h_AS0的不同元素，计算这些元素对应的熵值求和项，通过差量法将前一基准子图的联合熵H_AS0更新为当前基准子图的联合熵H_AS；

步骤7.计算实测图像A与当前基准子图的归一化互信息NMI，计算公式为：

NMI＝(H_A+H_S)/H_AS

步骤8.遍历整个基准图，以互信息最大值对应的位置作为最终匹配点。

2.根据权利要求1所述基于统计相关的快速互信息图像匹配方法，其特征在于，所述步骤4中按照差量法将前一基准子图的直方图矩阵h_S0更新为当前子图的直方图矩阵h_S，按如下步骤进行：

步骤4.1在列方向，通过统计当前基准子图最后一行和前一基准子图第一行对应像素灰度值出现次数，以前一基准子图的直方图矩阵h_S0为基准，减去前一基准子图第一行对应像素灰度值出现次数，加上当前基准子图最后一行对应像素灰度值出现次数，即得到当前子图的直方图矩阵h_S；

步骤4.2在行方向，通过统计当前基准子图最后一列和前一基准子图第一列对应像素灰度值出现次数，以前一基准子图的直方图矩阵h_S0为基准，减去前一基准子图第一列对应像素灰度值出现次数，加上当前基准子图最后一列对应像素灰度值出现次数，即得到当前子图的直方图矩阵h_S。

3.根据权利要求1所述基于统计相关的快速互信息图像匹配方法，其特征在于，所述步骤5中通过差量法将前一基准子图的联合熵H_S0更新为当前基准子图的联合熵H_S，按如下步骤进行：

步骤5.1找出h_S与h_S0的不同元素，分别组成矩阵h_ΔS和h_ΔS0，计算h_ΔS和h_ΔS0对应的信息熵Δ_S及Δ_S0；

Δ_S计算公式为：

Δ_S0计算公式为：Δ_S0＝-∑_snorm_s0(s)

其中，norm_s0(s)为与h_ΔS0中元素对应的前一基准子图熵值计算过程中的求和项，此次只计算这些求和项的和，无需对各个求和项再进行重复计算；

步骤5.2以前一基准子图的信息熵H_S0为基准，通过差量法将H_S0更新为当前基准子图的信息熵H_S，H_S＝H_S0-Δ_S0+Δ_S。

4.根据权利要求1所述基于统计相关的快速互信息图像匹配方法，其特征在于，所述步骤6中通过差量法将前一基准子图的联合熵H_AS0更新为当前基准子图的联合熵H_AS，按如下步骤进行：

步骤6.1找出h_AS与h_AS0的不同元素，分别组成矩阵h_ΔAS和h_ΔAS0，计算h_ΔAS和h_ΔAS0对应的信息熵Δ_AS及Δ_AS0；

Δ_AS计算公式为：

Δ_AS0计算公式为：Δ_AS0＝-∑_snorm_AS(s)

其中，norm_AS(s)为与h_ΔAS0中元素对应的前一基准子图联合熵计算过程中的求和项，此次只计算这些求和项的和，无需对各个求和项再进行重复计算；

步骤6.2以前一基准子图的信息熵H_AS0为基准，通过差量法将H_AS0更新为当前基准子图的信息熵H_AS，H_AS＝H_AS0-Δ_AS0+Δ_AS。