Fast Mutual Information Based Matching Algorithm in Terms of Gray-level Statistical Correlation of Adjacent Reference Sub-images

  • FU Yanjun 1 ,
  • SUN Kaifeng 2
Expand
  • 1 Xi’an International University, Xi’an 710077, China
  • 2 Xi’an Precision Machinery Research Institute, Xi’an 710075, China

Received date: 2019-08-30

  Online published: 2025-05-30

Abstract

Aiming at the poor real-time performance of the matching algorithm caused by the large amount of mutual information (mutual information, MI for short) computation, a fast MI computation method is proposed in this paper. In the matching process, the gray statistics of the previous matching position reference sub-image is updated to the current position by using the correlation between the adjacent reference sub-images, which avoids a lot of repeated calculation. Experimental results show that, compared with the standard method, the average matching time consumed is reduced by more than 20% without affecting the matching probability and matching accuracy. Furthermore, the proposed method can be combined with the existed gray-level compression based method or multi-resolution based method to further improve the real-time performance of matching.

Cite this article

FU Yanjun , SUN Kaifeng . Fast Mutual Information Based Matching Algorithm in Terms of Gray-level Statistical Correlation of Adjacent Reference Sub-images[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2020 , 40(6) : 5 -10 . DOI: 10.15892/j.cnki.djzdxb.2020.06.002

0 引言

在基于视觉的飞行器导航与制导、遥感卫星灾害监控与环境监测、医学图像分析等应用中, 常常需要对不同成像传感器获取的异源图像进行匹配。
基于互信息测度的匹配方法是一种完全基于图像灰度统计概率的匹配方法,不需要对原图像间的灰度关系作任何假设,非常适合于异源图像的匹配,但互信息测度最大的缺陷是计算量大, 匹配耗时较长[1]。为此,研究者们就围绕减少互信息测度计算量以及改进搜索策略等方面提出了相应的加速方法。文献[2-5]分别采用模拟退火、粒子群等优化搜索算法提高匹配速度;文献[6]采用灰度压缩和粒子群优化算法相结合加快匹配速度;文献[7]通过提取多层次特征点减少了互信息测度的计算量;文献[8-10]采用由粗到精的多层匹配方法对算法进行加速;文献[11-15]将分层匹配与各种优化算法相结合提高匹配速度。
文中利用匹配过程中相邻基准子图间灰度统计特征之间的相关性,通过差量法减少每一个匹配位置互信息的计算量来加快匹配速度。与已有的以灰度压缩、特征提取或采取优化搜索策略等快速互信息匹配方法相比,该方法既没有减少参与匹配的像素数,也无需对图像灰度等级进行处理,不会对匹配精度造成影响。

1 图像熵及互信息匹配

图像熵描述了图像信源的平均信息量,反映了图像灰度的统计信息,相似的两幅图像其图像熵也相近。利用互信息进行图像匹配的实质是:当两幅图像在空间位置配准时,其重叠部分所对应像素对的灰度互信息达到最大值。 两幅图像 U V的互信息 I ( U , V )可以用下式表示:
I ( U , V ) = H ( U ) + H ( V ) - H ( U , V )
其中:H(U)、H(V) 分别表示图像UV的熵; H ( U , V )表示两幅图像的联合熵。熵通常由变量的概率密度来表示。针对灰度图像,可以分别用两幅图像的灰度直方图和联合直方图进行估计。 计算公式分别为:
H ( U ) = - i = 0 255 h U ( i ) T l o g h U ( i ) T
H ( V ) = - i = 0 255 h V ( i ) T l o g h V ( i ) T
H ( U , V ) = - j = 0 255 i = 0 255 h U V ( i , j ) T l o g h U V ( i , j ) T
式(2)~式(4)中, h U ( i )表示图像U中灰度值为 i的像素出现的次数; h V ( i )表示图像 V中灰度值为 i的像素出现的次数; h U V ( i , j )表示同一像素位置上图像U中灰度值为 i、图像 V中灰度值为j的像素对出现的次数; T为图像的总像素数。
研究发现[16],互信息本身的大小与待配准两图像间的重叠度有一定的关联性,为了消除这种关联关系,文献[17]提出了利用归一化互信息(normalized mutual information, NMI)作为相似性测度。NMI能够减少对图像重叠部分的敏感性,实验表明它比标准的互信息方法更具鲁棒性。NMI表达式为:
N M I ( U , V ) = ( H ( U ) + H ( V ) ) / H ( U , V )
在飞行器导航中,匹配的目的是确定飞行器的当前位置,通过计算实测图与基准图中每一个匹配位置对应基准子图的互信息并找出最大互信息值对应位置,即可确定基准图与实测图的相对位置, 实现对飞行器的定位。

2 快速互信息匹配算法

与基于特征的图像匹配算法相比,基于最大归一化互信息的图像匹配最大缺点是计算速度慢。 每计算一次互信息测度,需要遍历图像中的所有像素,这使得计算量大大增加。由式(5)可以看出,归一化互信息的计算量主要集中在两幅图像各自的熵及其联合熵的计算。从式(2)~式(4)可以看出,图像熵及联合熵均为求和运算,其计算量主要集中于每一个求和项的计算,因此,如果能减少每次匹配过程中需要计算的求和项的项数, 则可以加快匹配速度。因为图像分布具有块状结构,在遍历式搜索匹配过程中,各相邻基准子图的灰度统计量具有很强的相关性,基准图像 ( i , j )位置处对应子图的灰度直方图与其四邻域位置处对应子图的灰度直方图非常相近。
图1所示,图1(a)为原始基准图像,图1(b)图1(a)的一个子图,图1(c)图1(b)在原基准图上右移一个像素对应的子图,图1(d)图1(e)分别为图1(b)图1(c)的灰度直方图。从图中可以看出,相邻两个子图的直方图非常相似,也即两相邻子图存在大量的灰度值其出现的频率没有发生变化,相应的,在式(2)~式(4)的求和项中分别有很多项数没有发生变化,因此, 当前匹配位置处基准子图的熵以及与实测图的联合熵的计算可以前一匹配位置处基准子图的熵及与实测图的联合熵的计算为基准,通过比较相邻两个匹配位置处基准子图对应的直方图及联合直方图的变化情况,用差量法进行计算。即在整个匹配过程中,除了要完整计算基准图上第一匹配位置(1,1)处对应基准子图的熵以及与实测图的联合熵外,其它所有匹配位置的基准子图的熵及与实测图的联合熵的计算主要集中在与前一匹配位置相比灰度直方图矩阵和联合直方图矩阵中发生变化的那些元素对应的求和项上。据此, 文中互信息匹配算法的流程图如图2所示。具体步骤如下:
图1 相邻基准子图及其灰度直方图
图2 文中匹配算法流程图
1)获取实测图像A和基准图像R,将两幅图像的灰度值调整到同一灰度区间。
2)按照式(2)或式(3)计算A的信息熵 H ( A )并存储。
3)从基准图像R的左上角(1,1)点截取与实测图像A大小相等的第1幅基准子图 S 1,统计 S 1各灰度值出现次数并存入初始灰度直方图矩阵 h S 1,按照式(2)或式(3)计算 S 1的信息熵 H ( S 1 )并存储,最后将 h S 1 H ( S 1 )分别赋值给 h S 0 H ( S 0 )
4)统计实测图A与第1幅基准子图 S 1对应像素位置的灰度对出现次数并存入初始联合直方图矩阵 h A S 1,按照式(4)计算初始联合熵 H ( A S 1 )并保存,将 h A S 1 H ( A S 1 )分别赋值给 h A S 0 H ( A S 0 );再按照式(5)计算实测图与基准子图的归一化互信息值并存储。
5)从基准图上(1,2)点位置开始,逐行逐像素按照“Z”字形进行迭代匹配。匹配过程中, 当前匹配位置记为 ( i ,   j ),当 j 1时,当前匹配位置的前一匹配位置应为 ( i ,   j - 1 );当 j = 1时,则当前匹配位置的前一匹配位置应为 ( i - 1 ,   j ),则当前位置 ( i ,   j )处的归一化互信息可按如下步骤计算:
①通过比较当前匹配位置 ( i ,   j )处基准子图和前一匹配位置 ( i ,   j - 1 ) ( i - 1 ,   j )处基准子图对应列或行像素变化情况,按照差量法将前一基准子图的直方图矩阵 h S 0更新为当前子图的直方图矩阵 h S。具体做法为: 在行方向匹配时,当前匹配位置的前一匹配位置为 ( i ,   j - 1 ),则以 h S 0为基准,减去前一基准子图第一列对应像素灰度值出现次数,加上当前基准子图最后一列对应像素灰度值出现次数,即得到 h S;在列方向匹配时,当前匹配位置的前一匹配位置为 ( i - 1 ,   j ),则以 h S 0为基准,减去前一基准子图第一行对应像素灰度值出现次数,加上当前基准子图最后一行对应像素灰度值出现次数,即得到 h S
②找出 h S h S 0的不同元素,分别组成矩阵 h Δ S h Δ S 0,按照式(6)、式(7)分别计算 h Δ S h Δ S 0对应的信息熵 Δ S Δ S 0;其中, n o r m S 0 ( S )为与 h Δ S 0中元素对应的前一基准子图熵值计算过程中的求和项,因为这些求和项在前一匹配位置已经计算过,此次无需再重复计算,而只需计算这些项的和即可。
Δ S = - S h Δ S ( S ) T l o g h Δ S ( S ) T
Δ S 0 = - S n o r m S 0 ( S )
③以前一基准子图的信息熵H(S0)为基准,按照式(8)通过差量法将H(S0)更新为当前基准子图的信息熵H(S)。
H ( S ) = H ( S 0 ) - Δ S 0 + Δ S
④计算实测图与当前基准子图的联合直方图矩阵 h A S,找出 h A S和前一基准子图联合直方图矩阵 h A S 0的不同元素,分别组成矩阵 h Δ A S h Δ A S 0,按照式(8)、式(9)分别计算 h Δ A S h Δ A S 0对应的联合熵 Δ A S Δ A S 0;其中, n o r m A S 0为与 h Δ A S 0中元素对应的前一基准子图联合熵计算过程中的求和项,因为这些求和项在前一匹配位置已经计算过,此次无需再重复计算,而只需计算这些项的和即可。
Δ A S = - S h Δ A S ( S ) T l o g h Δ A S ( S ) T
Δ A S 0 = - S n o r m A S 0 ( S )
⑤以前一基准子图的联合熵 H ( A S 0 )为基准,按照式(11)通过差量法将 H ( A S 0 )更新为当前基准子图的联合熵 H ( A S )
H ( A S ) = H ( A S 0 ) - Δ A S 0 + Δ A S
⑥按照式(12)计算实测图与当前基准子图的归一化互信息。
N M I ( A , S ) = ( H ( A ) + H ( S ) ) / H ( A , S )
6) 遍历整个基准图,选取归一化互信息最大值对应的位置作为最终匹配点。

3 算法时间复杂性分析

文中算法与传统方法相比,计算量的差别主要体现在每一个匹配位置处基准子图各统计量的计算上,为此只分析基准子图各统计量计算的时间复杂度。
设基准图尺寸为 N × N,实测图尺寸为 n × n,则基准子图大小与实测图大小一致均为 n × n
1)基准子图直方图计算
每一个匹配位置处,传统算法需要遍历整个基准子图来计算该子图的直方图,而文中方法只须遍历基准子图的第一行(或列)和最后一行(或列)即可获得该子图直方图。对于整个匹配过程来说,传统方法其时间复杂度为 O ( ( N - n ) 2 n 2 ),文中方法时间复杂度为 O ( ( N - n ) 2 n )
2)基准子图熵及联合熵计算
由于相邻基准子图其直方图发生变化情况跟图像内容有关,很难从理论上推导两子图间有多少个灰度值其出现频率会发生变化,也很难推导出其联合直方图中有多少个灰度值对的频率会发生变化。以直方图为例,最坏情况时有 2 n( n为基准子图各边长所含像素个数)个灰度值频率发生变化,且 2 n不超过基准图灰度等级数。
为此,文中在标准图像集上进行了大量的仿真实验,通过对仿真结果进行统计得出:相邻子图间灰度直方图中发生变化的元素平均为50%左右,联合灰度直方图中发生变化的元素平均为75%左右,这就意味着文中方法相比传统方法,其基准子图熵值的计算量减少了近50%,联合熵计算量减少了约25%,时间复杂度或时间频度有所降低。部分实验结果列于表1~表4
表1 SAR与可见光图像在列方向匹配时部分位置计算量统计结果
实测图尺寸/像素 基准子图左上角点位置 文中差量法 传统方法
M M1 M2 M M1 M2
150×150 (2,1) 302 123 10 597 22 801 249 14 326
(3,1) 302 118 10 455 22 801 249 14 492
(4,1) 302 114 10 557 22 801 249 14 392
(5,1) 302 120 10 301 22 801 250 14 449
(6,1) 302 119 10 288 22 801 249 14 455
100×90 (30,80) 182 105 5 490 9 191 249 7 278
(30,81) 182 102 5 518 9 191 249 7 154
(30,82) 182 101 5 534 9 191 248 7 235
(30,83) 182 99 5 485 9 191 248 7 251
(30,84) 182 118 5 512 9 191 247 7 210
表2 SAR与可见光图像在行方向匹配时部分位置计算量统计结果
实测图尺寸/像素 基准子图左上角点位置 文中差量法 传统方法
M M1 M2 M M1 M2
150×150 (2,1) 302 128 10 805 22 801 251 14 650
(3,1) 302 124 10 648 22 801 250 14 485
(4,1) 302 119 10 778 22 801 250 14 691
(5,1) 302 118 10 585 22 801 249 14 479
(6,1) 302 129 10 789 22 801 251 14 576
100×90 (30,80) 202 114 5 568 9 191 251 7 324
(30,81) 202 95 5 602 9 191 251 7 305
(30,82) 202 106 5 565 9 191 251 7 363
(30,83) 202 100 5 461 9 191 250 7 252
(30,84) 202 100 5 561 9 191 251 7 330
表3 差量法与传统标准法在SAR与可见光图像上的匹配结果
实测图尺寸/像素 真实匹配位置 文中差量法 传统方法
仿真匹配位置 匹配耗时/s 仿真匹配位置 匹配耗时/s
150×150 (180,150) (180,150) 193.121 4 (180,150) 273.313 8
(100,50) (100,50) 196.565 8 (100,50) 273.844 6
100×90 (90,30) (90,30) 164.165 2 (90,30) 216.169 2
(60,120) (60,120) 164.018 3 (60,120) 213.864 2
表4 差量法结合多分辨率分层算法与单纯多分辨率分层算法在SAR与可见光图像上的匹配结果
实测图尺寸/像素 真实匹配位置 差量法与多分辨率分层匹配相结合算法 单纯多分辨率分层匹配算法
仿真匹配位置 匹配耗时/s 仿真匹配位置 匹配耗时/s
150×150 (180,150) (180,150) 63.615 (180,150) 92.153
(100,50) (100,50) 60.352 (100,50) 89.623
100×90 (90,30) (90,30) 47.796 (90,30) 69.092
(60,120) (60,119) 47.593 (60,119) 69.103
另一方面,因文中方法需要存储前一匹配位置基准子图的多个统计量,算法的空间复杂度有所增加,但由于硬件技术的发展大大提高了计算机的存储容量,使得存储容量的局限性对于算法的影响大大降低。

4 仿真实验

为了验证算法的有效性,将文中方法应用于SAR与可见光图像匹配,以Matlab R2015为仿真平台,进行两类实验:一类是通过比较每个匹配位置处文中方法与传统遍历式互信息匹配方法在计算量方面的差别以及总耗时来说明文中方法的有效性;另一类是将文中方法与多分辨率分层匹配相结合,与单纯的多分辨率方法进行匹配耗时及匹配精度的比较,其中的多分辨率分层匹配采用d3小波基,小波分解级数为1。限于篇幅,文中只列出部分结果。
匹配过程中设: 基准子图直方图计算中需要统计的像素灰度值的次数为 M,基准子图熵值计算中需要计算的求和项的项数为 M 1,基准子图与实测图联合熵计算中需要计算的求和项项数为 M 2
选取可见光与SAR图像进行匹配,如图3所示,其中图3(a)为可见光基准图,大小为422×358;图3(b)为SAR实测图像1,大小为150×150,图3(c)为SAR实测图像2,大小为100×90,图3(d)为基于文中差量法匹配时实测图在基准图上的匹配定位结果;表1为差量法与传统标准算法列方向匹配时部分匹配位置计算量统计结果的比较,表2为差量法与传统标准算法行方向匹配时部分匹配位置计算量统计结果的比较,表3为文中差量法与传统标准方法匹配结果和匹配耗时的比较,表4为差量法结合多分辨率分层匹配法与单纯多分辨率分层匹配法的匹配结果和匹配耗时的比较。
图3 SAR与可见光图像匹配
图3可以看出,采用归一化互信息作为匹配测度,可以实现异源图像之间的匹配;由表1~表3可以看出,相比传统方法,文中差量法能够有效减少除第一个匹配位置外其它所有匹配位置互信息测度的计算量,总匹配耗时平均减少20%以上;由表4可以看出,文中差量法可以与多分辨率分层匹配技术相结合,其匹配精度与单纯的多分辨率分层匹配算法精度相同(在有些情况下匹配精度有所下降),但其匹配耗时平均下降20%左右。

5 结语

与已有的以灰度压缩、特征提取或采取优化搜索策略等快速互信息匹配方法不同,文中基于图像的块状结构特性以及匹配过程中相邻子图统计特性的相关性, 通过差量法减少互信息测度的计算量来加快匹配速度, 其互信息测度计算精度及匹配精度与原始的基于标准互信息计算方法的计算精度及匹配精度相同, 但计算量大为减少。 另外,如果参与匹配的两幅图像为高分辨率清晰图像,则可以把文中方法与灰度压缩、多分辨率分层匹配等方法相结合,即先对待匹配的图像进行灰度压缩或多分辨率分解,然后采用遍历式方法进行搜索,而在每一个搜索点处采用文中的差量法计算互信息测度,这样可以进一步减少匹配耗时。
[1]
张见威, 韩国强. 基于互信息的医学图像配准中互信息的计算[J]. 生物医学工程学杂志, 2008, 25(1):12-17.

[2]
DAME A, MARCHAND E. Second-order optimization of mutual information for real-time image registration[J]. IEEE Transactions on Image Processing, 2012, 21(9):4190-4203.

DOI PMID

[3]
ZHUANG Y W, GAO K, MIU X H, et al. Infrared and visual image registration based on mutual information with a combined particle swarm optimization-powell search algorithm[J]. Optik, 2016, 127(1): 188-191.

[4]
杜晓刚, 党建武, 王阳萍, 等. 基于萤火虫算法的互信息医学图像配准[J]. 计算机科学, 2013, 40(7):273-276.

[5]
柏连发, 韩静, 张毅, 等. 采用改进梯度互信息和粒子群优化算法的红外与可见光图像配准算法[J]. 红外与激光工程, 2012, 41(1):248-254.

[6]
安如, 王慧麟, 王盈, 等.16 阶归一化互信息和改进PSO算法的快速图像匹配[J]. 吉林大学学报(工学版), 2013, 43(增刊1):357-364.

[7]
杨猛, 潘泉, 张绍武, 等. 基于定量定性互信息的多层次特征图像匹配算法[J]. 中国图象图形学报, 2010, 15(9):1376-1383.

[8]
索庆亮, 王硕. 最大互信息在模版匹配中的应用[J]. 电光与控制, 2013, 20(4):84-86.

[9]
吴泽鹏, 郭玲玲, 朱明超, 等. 结合图像信息熵和特征点的图像配准方法[J]. 红外与激光工程, 2013, 42(10):2846-2852.

[10]
尉强, 刘忠. 基于互信息量的小幅员目标图像匹配[J]. 飞航导弹, 2018(3):66-69.

[11]
余春超, 杨智雄, 夏宗泽, 等. 采用GPU并行架构的基于互信息和粒子群算法的异源图像配准[J]. 红外技术, 2016, 38(11):938-946.

[12]
周鸣, 朱振福. 基于互信息和二级搜索的图像配准[J]. 微型机与应用, 2013, 32(9):1-4.

[13]
马宾, 马德贵, 廖建敏. 基于多分辨率的多模态医学图像配准[J]. 量子电子学报, 2012, 29(1):15-20.

[14]
冯林, 张名举, 贺明峰, 等. 用分层互信息和薄板样条实现医学图像弹性自动配准[J]. 计算机辅助设计与图形学学报, 2005, 17(7):1492-1496.

[15]
KUTTEN K S, CHARON N, MILLER M I, et al. A large deformation diffeomorphic approach to registration of CLARITY images via mutual information[J]. Lecture Notes in Computer Science, 2017, 10433(9): 275-282.

[16]
PLUIM J P W, MAINTZ J B A, VIERGEVER M A. Mutual information based registration of medic al images: a survey[J]. IEEE Transactions on Medical Imaging, 2003, 22(8): 986-1004.

[17]
KNOPS Z F, MAINTZ J B, VIERGEVER M, et al. Normalized mutual information based registration using k-means clustering and shading correction[J]. Medical Image Analysis, 2006, 10(3): 432-439.

PMID

Outlines

/