Robust Multimodal Remote Sensing Image Registration Based on Local Statistical Frequency Information

Abstract

Multimodal remote sensing image registration is a prerequisite for comprehensive application of remote sensing image data. However, inconsistent imaging environment and conditions often lead to obvious geometric deformations and significant contrast differences between multimodal remote sensing images, which makes the common feature extraction extremely difficult, resulting in their registration still being a challenging task. To address this issue, a robust local statistics-based registration framework is proposed, and the constructed descriptors are invariant to contrast changes and geometric transformations induced by imaging conditions. Firstly, maximum phase congruency of local frequency information is performed by optimizing the control parameters. Then, salient feature points are located according to the phase congruency response map. Subsequently, the geometric and contrast invariant descriptors are constructed based on a joint local frequency information map that combines Log-Gabor filter responses over multiple scales and orientations. Finally, image matching is achieved by finding the corresponding descriptors; image registration is further completed by calculating the transformation between the corresponding feature points. The proposed registration framework was evaluated on four different multimodal image datasets with varying degrees of contrast differences and geometric deformations. Experimental results demonstrated that our method outperformed several state-of-the-art methods in terms of robustness and precision, confirming its effectiveness.

1. Introduction

Nowadays, advanced remote sensing technology can realize omni-directional and multi-granularity perception of the same target scene under different imaging environment and conditions [1]. In the meanwhile, rich multimodal remote sensing images can be acquired with different sensors, or in different time periods [2,3]. A prerequisite for comprehensive utilization of multimodal image information is accurate image registration, which has been widely used in many computer vision tasks such as image mosaic, image fusion, change detection, vision-based navigation, scene matching guidance [4]. Therefore, the accuracy of image registration is critical to the application effect of the above top-level tasks.

Multimodal remote sensing image registration is the process of aligning remote sensing images of the same scene taken by different sensors, at different times, or/and from distinct viewpoints [5]. Although many related methods have been proposed [6,7,8,9,10,11,12], it is still very challenging work due to the contrast inconsistency caused by differences in imaging environment (different sensors, weather, and time periods) and the large geometric deformations induced by imaging conditions (different platform attitudes and positions) (Figure 1). Recent mainstream methods have been applied successfully in the situation where the geometric changes are small [6,7,8,9] or can be greatly alleviated according to the capture information [10,13]. However, automatic multimodal remote sensing image registration has not been solved effectively in complicated environments with large geometric changes and significant contrast differences. The challenges in the accurate registration of multimodal remote sensing images are specifically analyzed as follows:

• Differences in imaging mechanisms or different weather capture conditions cause non-linear radiation changes between images, which leads to significant contrast differences, rendering traditional feature representation methods based on grayscale or gradient less effective or even invalid.
• Large geometric deformations occur between images acquired from different azimuths (viewpoints) or different platforms (airborne camera, space camera), which makes it extremely difficult to extract invariant features.
• Images acquired at different times or by different sensors contain some structural changes, resulting in poor consistency of the feature representation for the same target, making it difficult to achieve accurate registration.

To alleviate the difficulty mentioned above, this paper presents a robust multimodal remote sensing image registration framework using local statistical frequency information. The main contributions of our paper can be summarized as follows:

• The maximum phase congruency optimization method is proposed, which is a guarantee for stable structural feature localization in multimodal remote sensing images and determines the center of the feature description regions.
• To make full use of the frequency information and dig out structural features in image better, a joint local frequency information map that combines Log-Gabor filter responses over scales and orientations was constructed, which offers the main information to the feature descriptors.
• The geometric and contrast invariant descriptors were generated through the selection of feature scales and the orientation statistics of the region to be described on the joint feature map, which is critical to achieve accurate registration.

Due to the improvements above, the proposed method is robust for multimodal remote sensing image registration with large geometric variations and significant contrast differences (Figure 1). The rest of this paper is organized as follows: The related works of registration for multimodal remote sensing images are reviewed and the main differences of our proposed method from others are described in Section 2. Detailed description of the proposed registration framework using local statistical information is given in Section 3. Comparative experiments and analysis are performed in Section 4. Finally, conclusions are drawn in Section 5.

2. Related Works

Over the years, a variety of multimodal remote sensing image registration methods have been proposed, which can be classified into two categories: global intensity-based methods and local feature-based methods [14]. Global intensity-based methods obtain the optimum transformation parameters via maximizing the similarity of global intensities between input images, mainly including mutual information [6], cross correlation [15], phase correlation [16], Fourier transformation [17] and wavelet features [18]. These methods perform well for the images with high correlation in global intensity or in transform domain. However, contrast reversal, occlusion, small area overlaps, and clutters occur frequently in some regions of input images, which makes the intensity-based methods unable to achieve an accurate registration.

In contrast, local feature-based methods first achieve feature matching by extracting and comparing local features in images and then compute the transformation via the correspondence of features. The advantage of these methods is that they can deal with significant geometric deformations as well as contrast differences in images. These methods can be classified into three groups: typical feature-based methods, deep learning-based methods, and local structural feature-based methods. The main idea of typical feature-based methods is to manually design different feature extraction models according to different applied scenarios. The related typical feature extraction methods include contour-based [19], line-based [20], region-based [21], and gradient distribution-based [8,22,23,24,25]. These methods have a good performance in the case that corresponding features are obvious or their gradient distribution is consistent in input images. The localization performance of several local features was compared in [26]; the best result was obtained by Root-SIFT [25]. However, they treat all image content equally, so they are highly sensitive to structural disparities caused by insignificant structures, such as shadow, illuminance, and resolution. Therefore, it will lead to severe degradation of matching performance when large contrast differences appear in multimodal remote sensing images. In recent years, with the prosperous application of deep learning in computer vision field, deep learning-based methods have been developed for remote sensing image registration [27,28,29], which can automatically learn high-level semantic features and get better matching performance compared to typical feature-based methods in some complicated cases. Especially, superglue [30] outperforms other learned approaches and achieves state-of-the-art results on the task of pose estimation in challenging real-world indoor and outdoor environments. However, these methods cannot deal with large geometric deformation (scale and rotation), and their performance is usually affected by the size of training data. As a result, the multimodal registration problem still cannot be effectively solved by the current deep learning-based methods.

Different from the methods of the first two groups, local structural features-based methods can extract more robust common features from different modalities and are less sensitive to contrast differences. Due to these advantages, they have been successfully applied to multimodal image registration [9,10,31,32,33,34]. As a typical representative of structural features, phase congruency of local frequency information was first proposed by Morrone et al. [35]. Subsequently, to improve the robustness of phase congruency to noise and contrast, Kovesi proposed a new sensitivity measure [36] and a highly contrast invariant localized feature detector [37]. Recently, Liu et al. [9] proposed mean local phase angle and frequency spread phase congruency by using local frequency information on Log-Gabor wavelet transformation space, which improved the robustness compared with traditional multimodal matching. To extract the structural features, Ye et al. [10] developed the histogram of orientated phase congruency descriptor, which outperforms several methods in matching performance. Lately, they proposed channel features of orientated gradients (CFOG) [38] as an extension of [10] with superior performance in image matching and computational efficiency. Xie et al. [32] achieved multimodal image registration through combining phase correlation and multi-scale structural information.

Moreover, a structural information orientation histogram descriptor is constructed by concatenating the orientation of magnitude and the minimum moment [37]. Li et al. [33] proposed RIFT method for multimodal images registration via constructing a maximum index map (MIM) for feature description, which can achieve good performance on images with radiation-variation. However, the above-mentioned methods based on phase congruency are not enough to make full use of local frequency information, leading to limitations in processing large geometric deformation.

To address the registration problem of infrared and visible image, we constructed the maximally stable phase congruency (MSPC) descriptor using maximally stable extremal regions (MSER) [39] and local frequency information [31], nevertheless, in which MSER cannot obtain a higher repeatability for multimodal images. Subsequently, we developed a robust matching method for electro-optical by combining phase congruency with kernelized correlation filter (KCF) [40], which can first adjust the input images according to the platform parameters, and then get the registration results [13]. To deal with the multimodal remote sensing image registration with large geometric deformation and contrast difference, this paper extends our early works [13,31], by fully utilizing local frequency information. The main process of the proposed framework is shown in Figure 2. We first obtain the maximum phase congruency of local frequency information by optimizing the control parameters, and then extract the salient feature points according to their optimized phase congruency response. After that, we construct the geometric and contrast invariant descriptors (GCID) based on a joint local frequency information map (JLFM) that combines Log-Gabor filter responses over multiple scales and orientations. Finally, image matching and registration can be achieved by the correspondence of descriptors. In the experimental part, we performed a more thorough evaluation in terms of robustness and precision for the proposed registration framework with four different multimodal remote sensing image datasets.

3. Methodology

In this section, the robust multimodal remote sensing image registration framework using local statistical frequency information is presented. Section 3.1 introduces the calculation method of the maximum phase congruency of local frequency information through parameter optimization, and further realizes the location of salient feature points on the optimized feature map. Then, Section 3.2 presents the construction of GCID using JLFM in detail. Finally, the entire multimodal remote sensing image registration framework is given in Section 3.3.

3.1. Maximum Phase Congruency and Feature Detection

Phase congruency provides a measure that is independent of the overall magnitude of the signal making it invariant to variations in image illumination and/or contrast. To improve the insensitive to noise and provide good localization, Kovesi [36] proposed a new sensitivity measure and noise compensation method for phase congruency, which can locate the features that remain constant over scales. The new measure is calculated by the following formula:

$$PC(x,y)=\frac{{\displaystyle \sum _s{\displaystyle \sum _o{w}_o(x,y)\lfloor A_{so}(x,y)\Delta {\Phi }_{so}(x,y)-T\rfloor }}}{{\displaystyle \sum _s{\displaystyle \sum _o{A}_{so}(x,y)+\epsilon }}},$$

(1)

where $w_o$ is a factor that weights for frequency spread, $A_{so}(x,y)$ is an amplitude of the Fourier component at position $(x,y)$, and $\Delta {\Phi }_{so}$ is a phase deviation function. $\epsilon$ is a small constant to avoid the denominator being zero, and $T$ is a threshold that eliminates noise influence. The symbol $\lfloor \rfloor$ denotes that the enclosed quantity is equal to itself when its value is positive and zero otherwise. Based on this new measure, Kovesi [37] presented a highly localized feature detector whose responses are invariant to image contrast, which can be achieved as follows:

(1)

The moment analysis equations at each point are calculated as following:

$$A={\displaystyle \sum _o{\left(PC\left({\theta }_o\right)\cos \left({\theta }_o\right)\right)}^2},$$

(2)

$$B=2{\displaystyle \sum _o\left(PC\left({\theta }_o\right)\cos \left({\theta }_o\right)\right)\left(PC\left({\theta }_o\right)\sin \left({\theta }_o\right)\right)},$$

(3)

$$C={\displaystyle \sum _o{\left(PC\left({\theta }_o\right)\sin \left({\theta }_o\right)\right)}^2},$$

(4)

(2)

The maximum moment $M$ and minimum moment $m$ are given by,

$$M=\frac{1}{2}\left(A+C+\sqrt{B^2+{\left(A-C\right)}^2}\right),$$

(5)

$$m=\frac{1}{2}\left(A+C-\sqrt{B^2+{\left(A-C\right)}^2}\right).$$

(6)

A large value of $M$ indicates an edge feature point and A large value of $m$ means that point should be a corner; therefore, $M+m$ contains more features than anyone of them. To highlight features and to reduce computational complexity, we adopted the sum of phase congruency in multiple orientations as a candidate feature map (CFM) according to (7).

$$M+m={\displaystyle \sum _{o}P{C}^2\left({\theta }_o\right)}\le {\displaystyle \sum _{o}PC\left({\theta }_o\right)}$$

(7)

However, a very important factor rarely considered in the previous literature for phase congruency feature detection is the optimization and fine adjustment of control parameters. To improve robustness of phase congruency to nonlinear radiation deformations, RIFT [33] proposed MIM for feature description, nevertheless, other important structural information is lost. The literature [41] presented that the 2D-MSPC parameters can be optimally and automatically tuned by maximizing the norm of cost function formed by $M$ and $m$; however, their calculations are very complicated. In this paper, we use the similarity of CFMs extracted from two input images as the criterion of parameter optimization. The larger the similarity of structural features, the better the parameters.

There are eight parameters that mainly affect the performance of phase consistency, which can be denoted by $[N_s,N_o,{\lambda }_{\min },\eta ,\sigma ,k,C_o,g]$ and described as

Table 1. List of phase congruency parameters.

Parameters	Range	Default	Meaning
$N_s$	[3~6]	4	Number of wavelet scales.
$N_o$	[1~6]	6	Number of filter orientations.
${\lambda }_{\min }$	$\ge 3$	3	Wavelength of smallest scale filter
$\eta$	[1.3, 1.6, 2.1, 3]	1.6	Scaling factor between successive filters.
$\sigma$	[0.1~1]	0.55	Ratio of the standard deviation of the Gaussian.
$k$	[10~20]	2	Noise scaling factor.
$C_o$	(0~1]	0.5	The fractional measure of frequency spread.
$g$	[1~50]	10	Controls the sharpness of the transition in the sigmoid function.

. Among them, $N_s$ and $N_o$ are the number of filter scales and orientations respectively. Optimal values of $N_s=4$ and $N_o=6$ are verified through detailed experiments in [33]. The changes of ${\lambda }_{\min }$, $k$, $C_o$, and $g$ do not have significant effect on the phase consistency extraction results, which has been confirmed by our experiments. Therefore, the default values of those parameters are used in this paper. The values of scaling factor between successive filters $\eta$ and ratio of the standard deviation of the Gaussian $\sigma$ have a relatively large impact on the results, which can be seen from Figure 3.

The cosine similarity of CFMs extracted from the input images was used as the criterion of optimizing parameters. The specific optimization process can be described as follows:

(1)

CFMs are produced from the input images through changing the values of $\eta$ and $\sigma$, and keeping default values for other parameters.

(2)

Compute the cosine similarity of CFMs by the following formula:

$$CS(CF{M}_{ref},CF{M}_{sen})=\frac{{\stackrel{\rightharpoonup }{CFM}}_{ref}·{\stackrel{\rightharpoonup }{CFM}}_{sen}}{\Vert {\stackrel{\rightharpoonup }{CFM}}_{ref}\Vert ·\Vert {\stackrel{\rightharpoonup }{CFM}}_{sen}\Vert },$$

(8)

where $CF{M}_{ref}$ and $CF{M}_{sen}$ are the corresponding CFMs of reference image and sensed image respectively. ${\stackrel{\rightharpoonup }{CFM}}_{ref}$ and ${\stackrel{\rightharpoonup }{CFM}}_{sen}$ are their histogram statistical vectors.

(3)

The average cosine similarity of different CFMs for a group parameter is computed, and the optimal parameters can be determined according to the maximum cosine similarity of CFMs obtained by different parameter combinations.

In the algorithm above, each $\eta$ has 19 combinations by changing the value of $\sigma$ within the range of (0.1, 1) at 0.05 intervals. Therefore, the optimal parameters are determined by comparing the cosine similarity between 76 CFM pairs, those obtained from one multimodal remote sensing image pair with different parameters. In this paper, the optimal parameters for CFMs of a category image can be determined by their average cosine similarity. CFMs extracted from infrared and visible images by using optimized and default parameters are shown in Figure 3, from which we can see that the optimized CFMs (OCFMs) contained more structural information than those obtained by default parameters.

After the extraction of optimized CFMs from multimodal remote sensing image pairs, feature points will be detected from those CFMs by using FAST [42] and response ranking, which can be performed as follows:

(1)

FAST is applied on OCFMs obtained from input multimodal remote sensing image pairs to get plentiful candidate feature points.

(2)

To enhance the saliency of feature points, the extracted candidate feature points are ranked according to their response value in CFMs, and the top-k points will be selected as salient feature points.

(3)

To ensure the uniform distribution of feature points, non-maximum suppression is implemented on their $n\times n$ neighborhood.

Among the above feature point detection methods, OCFMs and response ranking can ensure the saliency of feature points, and non-maximum suppression can promote the uniformity of their distribution. Examples of feature point detection on a multimodal image pair are shown in Figure 3d and Figure 4, in which, the remote sensing images were acquired by different sensors (Infrared and Visible) or by the same sensor (Visible) but at different times, which leads to an obvious contrast difference between the input images. It can be seen from the results that the feature points detected by our proposed method had high saliency and repeatability and that their distributions were very uniform.

3.2. Construction of GCID

Salient feature points imply that there is important structural information around them. How to extract and describe structural information make descriptors have geometric deformation and contrast invariance, which are crucial for multimodal remote image registration. To explore the solution, GCID was constructed in this section. First, we built JLFM by combing Log-Gabor filter responses over scales and orientations, which could achieve a robust description for local feature under large contrast difference and geometric deformation. Then, the regions with scale and orientation invariance based on JLFM were located. Finally, GCID were generated by using BRIEF [43] in those located regions.

Li et al. [33] proposed that phase congruency maps are not suitable for feature description because they have small value and are sensitive to noise without obvious edge features in the image. Therefore, they presented the MIM measure instead of the PC map for feature description; however, the important frequency information that can reflect structural features is discarded. To make full use of frequency information and dig out structural features in an image better, we constructed a JLFM by combing Log-Gabor filter responses over scales and orientations as follows:

$$JLFM(x,y)=\frac{1}{s}{\displaystyle \sum _s{\displaystyle \sum _o\sqrt{E_{so}{(x,y)}^2+O_{so}{(x,y)}^2},}}$$

(9)

$$\left[E_{so}(x,y),O_{so}(x,y)\right]=\left[I\left(x,y\right)\ast L{G}^e\left(x,y,s,o\right),I\left(x,y\right)\ast L{G}^o\left(x,y,s,o\right)\right].$$

(10)

where $E_{so}(x,y)$ and ${\mathrm{O}}_{so}(x,y)$ are response components produced by convolving the image $I(x,y)$ with the even-symmetric and the odd-symmetric Log-Gabor wavelets. JLFM has two advantages for feature description. First, it can integrate multi-orientation and multi-scale local frequency information that reflects structural features (Figure 5c) more comprehensively than phase congruency (Figure 5b). Second, the invariance of geometric deformation can be achieved through statistics of the response of JLFM at different scales and orientations, which is an important guarantee to obtain accurate matching of multimodal remote sensing images (Figure 5d).

3.2.1. Scale-Invariant Description Region

Scale and rotation are the main geometric deformation between multimodal remote images. Morel et al. [44] proposed that complex global transformations can be locally approximated by similarity transformations. Therefore, the scale and rotation invariance of structural feature descriptors are the main problems to be solved in this paper. First, the scale invariance of the description region is selected by the following steps:

(1)

The scale of point $(x,y)$ can be computed as follows:

$${\sigma }_p(x,y)=\underset{s\in \{1,2,\dots N_s\}}{\arg \max }\left(A_s\left(x,y\right)\right),$$

(11)

$$A_s\left(x,y\right)={\displaystyle \sum _{\mathrm{o}}\sqrt{E_{so}{(x,y)}^2+O_{so}{(x,y)}^2}},$$

(12)

where ${\sigma }_p(x,y)$ is the assigned scale of point $(x,y)$ in the neighborhood of the feature point, which means that point $(x,y)$ has the maximum response at ${\sigma }_p(x,y)$ over $N_s$ scales. $A_s(x,y)$ is the sum of filter responses of point $(x,y)$ in all orientations at scale $s$.

(2)

In the $n\times n$ neighborhood of the feature point, count the number of points that have the same ${\sigma }_p$ and use the scale ${\sigma }_F$ with the largest number of points as the scale of the feature point, which can be formulated as:

$${\sigma }_F=\underset{(x,y)\in N_F,{\sigma }_p\in \{1,2,\dots N_s\}}{\arg \max }\left(Number({\sigma }_p\left(x,y\right))\right),$$

(13)

where $N_F$ is the $n\times n$ neighborhood of the feature point.

(3)

The central frequency of Log-Gabor filter controls their scales; therefore, reciprocal of that is adopted to determine the description region radius of a feature point as follows:

$$R_F=R_0\cdot {\lambda }_{\min }\cdot {\eta }^{\sigma -1},$$

(14)

where ${\lambda }_{\min }$ is the minimum wavelength, $\eta$ is the scale factor between successive filters, and $R_0$ is the initial radius. The above adaptive determination of the description regions can realize the scale invariance of the feature.

3.2.2. Rotation-Invariant Description Region

After the description regions are selected, the rotation invariance of the feature can be achieved by the statistics of histogram on the orientation map (ORM), which is defined as follows:

$$ORM(x,y)=\arctan \left(\frac{{\displaystyle \sum _o\cos \left(\theta \right){\displaystyle \sum _s{O}_{so}(x,y)}}}{{\displaystyle \sum _o\sin \left(\theta \right){\displaystyle \sum _s{O}_{so}(x,y)}}}\right),$$

(15)

where $\theta$ is the orientation of Log-Gabor filter and $O_{so}(x,y)$ is the response component generated by the odd-symmetric Log-Gabor wavelet.

Similar to SIFT [45], the dominant orientation of a feature point is determined by the histogram statistics in its description region on ORM, which can be described as follows:

(1)

To improve the stability of description on image contrast, the range of the histogram is (0~180°); therefore, the histogram contains 36 bins, and every 5° is counted as one bin. Each bin can be calculated by ORM and Gaussian weighted JLFM as follows:

$$hist(s)=hist(s)+JLFM(i,j)\ast W(i,j),s\in [0,255],(i,j)\in N_{R_F},$$

(16)

$$W(i,j)=\exp (-(i^2+j^2)/2\ast (1.5\ast {\sigma }_F))),$$

(17)

where $N_{R_F}$ is the description region of a feature point.

(2)

Smoothing of the histogram is performed; the highest peak of the histogram is taken as the dominant orientation; the second highest peak that exceeds 80% of the highest peak is regarded as the auxiliary orientation.

Figure 6 shows the calculation process of scale and orientation invariant features from visible and NIR images. Description regions of two corresponding points are obtained by using Formulas (10)–(13), which are shown in the middle part of Figure 6a. Orientation histograms of description regions based on ORM and JLFM are given in Figure 6b. From the results of histogram statistics, dominant orientations of their description regions were assigned as 27.5° and 105°, respectively. Rotated description regions of the two corresponding points on JLFMs are shown in Figure 6d. From those results, we can see that the content of the corresponding feature point description regions on JLFMs had good similarity and consistency. After adjusting the scale and orientation of the description region of feature points, BRIEF was used for the construction of GCID as follows:

In the description region of a feature point on JLFM, a pair of points were selected randomly and compared by the formula as follows:

$$\tau (p;x,y)=\{\begin{array}{cc}1,& i{f}_{}p(x)<p(y)\\ 0,& otherwise\end{array}$$

(18)

where $p(x)$ and $p(y)$ are the response values of random points $x=\left(u_1,v_1\right)$ and $y=\left(u_1,v_1\right)$, respectively. According to the above criteria, N pairs of random points were selected in the description region, and the binary assignment was performed to form a binary code, that is, BRIEF descriptor. The selection of random points in this paper obeyed the anisotropic Gaussian distribution; the Hamming distances were adopted to match the BRIEF descriptors.

3.3. Registration Framework by Using GCID

After the feature detection and GCID construction proposed in Section 3.1 and Section 3.2, the entire multimodal remote sensing image registration framework is given in this Section.

The workflow of the proposed multimodal remote sensing image registration framework is shown in Figure 7, and its steps can be described as follows:

(1)

OCFMs are first computed from $I_{Inf}$ and $I_{Sen}$ by Formulas (1)–(8), respectively and then feature points are detected by FAST and non-maximum suppression on OCFMs.

(2)

JLFMs are obtained from $I_{Inf}$ and $I_{Sen}$ by combing Log-Gabor filter responses over scales and orientations by Formulas (9) and (10).

(3)

GCIDs from $I_{Inf}$ and $I_{Sen}$ are generated by using JLFMs and feature points obtained by steps (1) and (2) according to Formulas (11)–(18), respectively.

(4)

Matching results of GCIDs from $I_{Inf}$ and $I_{Sen}$ are computed by their distance similarity; the outliers are removed by random sample consensus (RANSAC).

(5)

Transformation is estimated according to the matching results; the registration of $I_{Inf}$ and $I_{Sen}$ are achieved.

The image registration workflow shown above was fully automatic, where parameter optimization could be done based on 2–3 image pairs from the same modal. The registration framework proposed in Figure 7 mainly included three parts: feature detection (orange), feature description (aqua green), feature matching and registration (green). Feature points localization implies that there are significant structural features around them, while the construction of GCID realizes the invariant description of structural features about the feature points. The main innovative work of this paper was the detection of salient features and the description of invariant structural features.

4. Experiment Results and Analysis

To demonstrate the effectiveness of the proposed method, four multimodal datasets were employed in comparative and evaluative experiments in this section. Those different sets of images are introduced in Section 4.1. To show the effect of parameter optimization, experiments with parameter optimization in image matching are implemented in Section 4.2. Then, robustness of the proposed method to geometric deformation is tested under different degrees of rotation and scale changes in Section 4.3. Finally, the comparative experimental results of the proposed method and state-of-the-art methods (Root-SIFT [25], RIFT [33], CFOG [38], SuperGlue [30]) are analyzed in Section 4.4. In addition, more visual matching results obtained by the proposed method are given in this section.

To quantitatively evaluate the matching performance, performance measures such as precision, root mean square error (RMSE), and median error (MEE) were adopted, which can be expressed as:

$$\mathrm{Pr}ecision=\frac{NCM}{NTM},$$

(19)

where $NCM$ and $NTM$ is the number of correct matched and total matched point pairs, respectively.

$$RMSE=\sqrt{\frac{1}{L}{\displaystyle \sum _{i=1}^L{(p_{Inf}^i-T(p_{Sen}^i))}^2}},$$

(20)

$$MEE=median{\left\{\sqrt{{(p_{Inf}^i-T(p_{Sen}^i))}^2}\right\}}_{i=1}^L,$$

(21)

where $\left\{p_{Inf}^i,p_{Sen}^i\right\}$ indicates the corresponding matched point pair, $p_{Inf}^i$ and $p_{Sen}^i$ are the feature point coordinates extracted from inference image and sensed image, respectively. $T$ is the ground-truth transformation between the two images, which can be obtained by manually selecting the corresponding point pairs. $L$ represents the number of matched point pairs, and $median\{\cdot \}$ returns the median value of a set.

4.1. Multimodal Remote Sensing Datasets

Multimodal image datasets (Figure 8) employed in comparative and evaluative experiments consist of the following four parts:

(1)

Remote sensing dataset [14]: the dataset contains 78 image pairs, which can be divided into 7 modal types, such as UAV cross-season images, visible day-night, LiDAR depth-optical, infrared-optical, map-optical, optical cross-temporal, and SAR-optical images. These images have different resolutions ranging from $500\times 500$ to $713\times 417$, while the corresponding images have the same resolution; therefore, differences in contrast and inconsistencies in detail are the main changes between them.

(2)

Computer vision dataset [14]: this dataset contains 54 image pairs, which includes 4 modal types, such as visible-infrared images, visible cross-season, day-night, and RGB-NIR images. These images have different resolutions ranging from $256\times 256$ to $2133\times 1600$; the corresponding images have the same resolution. Contrast difference and geometric deformation are the main changes between the image pairs.

(3)

UAV dataset [13]: those visible and infrared images were captured at the same time from EOP on UAV, which consisted of 160 image pairs with discontinuous focus length change from 25 to 300 mm for the infrared camera and from 6.5 to 130.2 mm for the visible camera. The infrared images were captured by a mid-wavelength infrared camera operated in the 3–5 μm waveband with a size of 640 × 512. The visible images were captured by a lightweight CCD camera with a size of 1024 × 768. Therefore, large geometric deformation and contrast differences occurred between those image pairs.

(4)

NIR-VIS dataset: this dataset was captured by Gaofen-14 satellite, which contained 40 image pairs. Near infrared (NIR) images were taken by a medium wavelength infrared camera and the visible images were taken by a visible camera; contrast difference is the main change between those image pairs.

Examples of the above four image datasets are shown in Figure 8. The dataset has a total of 9 modalities and 332 image pairs. We can see that geometric deformation and contrast changes between multimodal images are obvious. Their ground truth transformations were determined in advance by manual registration, and the match was considered correct if the RMSE was less than 3; otherwise, it was false.

4.2. Parameter Optimization Experiments

According to the parameter optimization description in Section 3.1, this section gives matching performance experiments before and after optimization to demonstrate the effectiveness of the method. The parameters to be optimized and their ranges are given in

Table 2. The range of optimization parameters.

Parameters	Optimization Range
$\sigma$	[0.1, 0.15,…, 0.95, 1] (Interval 0.05)
$\eta$	[1.3, 1.6, 2.1, 3]

; from the table, we can see that the determination of a pair of optimal parameters $[\sigma ,\eta ]$ required the calculation of 76 combinations. In this experiment, CFMs were first extracted by calculating the sum of phase congruency in multiple orientations, and then the optimal parameters could be determined according to the maximum of the cosine similarity of CFMs obtained by different parameter combinations.

Three types of multimodal remote sensing images (infrared-optical, SAR-optical, and depth-optical) were selected to test the matching performance in the case of optimized parameters adopted in the process of feature extraction. The matching results before and after parameter optimization are shown in Figure 9; from the results, we can see that the number of correct matching point pairs was significantly increased, and that the distribution was more uniform after optimization. In addition, four types of multimodal remote sensing images were employed to measure the performance of optimized parameters in terms of average NCM. The comparative results are shown in

Table 3. Average NCM in different modal datasets by using optimized and default parameters.

Datasets	Optimized Parameters $[{\mathit{\eta }}_{\mathit{I}\mathit{n}\mathit{f}},{\mathit{\sigma }}_{\mathit{I}\mathit{n}\mathit{f}},{\mathit{\eta }}_{\mathit{S}\mathit{e}\mathit{n}},{\mathit{\sigma }}_{\mathit{S}\mathit{e}\mathit{n}}]$	Average NCM	Default Parameters $[{\mathit{\eta }}_{\mathit{I}\mathit{n}\mathit{f}},{\mathit{\sigma }}_{\mathit{I}\mathit{n}\mathit{f}},{\mathit{\eta }}_{\mathit{S}\mathit{e}\mathit{n}},{\mathit{\sigma }}_{\mathit{S}\mathit{e}\mathit{n}}]$	Average NCM
Infrared-Optical	[1.6, 0.50, 1.6, 0.60]	86	[1.6, 0.70, 1.6, 0.70]	63
Depth-Optical	[1.3, 0.50, 1.6, 0.45]	204	[1.6, 0.70, 1.6, 0.70]	101
SAR-Optical	[1.6, 0.45, 1.6, 0.45]	96	[1.6, 0.70, 1.6, 0.70]	58
Map-Optical	[1.6, 0.45, 2.1, 0.65]	147	[1.6, 0.70, 1.6, 0.70]	61

, which shows a double times improvement in average NCM after parameter optimization over the pre-optimization.

4.3. Geometric Deformation Resistance Experiments

To verify the robustness of the proposed method to geometric deformations, visible and infrared image dataset captured from UAV and remote sensing image dataset were used as experimental data, which had large geometric deformations and contrast differences. Since rotation and scale transformation are the most important geometric deformations, this section mainly tested the robustness of our method to these two deformations.

4.3.1. Rotation Robustness Test

In this part, five images pairs were randomly selected from remote sensing and UAV datasets first, and then, for each pair of images, the sensed image was rotated from 30° to 180° at intervals of 30° to generate 6 rotated images. Therefore, 30 pairs of new images were obtained and used to test the robustness of the proposed method to rotational changes. The number of correct matches (NCM) of those images is shown in Figure 10; from the results, we can see that rotation had little effect on the matching results. The average fluctuation of NCM with rotation did not exceed 30%, which was mainly due to the rotation-invariant design of the proposed method. The matching results of the fifth image pair at different rotation angles are shown in Figure 11. The average NCM of the six rotated image pairs was more than 50, which was enough for accurate image registration. The proposed method successfully matched the rotated image pairs, although the sensed images had scale and rotation changes, which confirmed the effectiveness of the proposed method on rotation changes.

4.3.2. Scale Robustness Test

Robustness to scale changes is key to testing the performance of registration methods; therefore, 20 visible and infrared image pairs from the UAV dataset were employed to evaluate the matching performance of the proposed method in this section.

We first chose a visible image as the inference image and then sequentially selected 20 infrared images with discontinuously varying focus lengths from 25 to 300 mm in the same scene as the sensed images. The scale difference between the input images changed from 0.8 to 3.2, examples of which are shown in Figure 12. The upper left image is the reference image, and the other images with scale changes are the sensed images, the scale of which gradually became larger. The NCM of visible image and 20 infrared image pairs selected from the UAV dataset is given in Figure 13. Although the NCM decreases as the scale becomes larger, the NCM obtained by our method was all above 20, which could meet the needs of image registration. Matching results of Figure 12 by using the proposed method are shown in Figure 14; from the results, we can see that the proposed method could handle large-scale changes due to the scale invariance processing in the descriptor GCID. The registration results in the lower right corner demonstrate the good performance of the proposed method under large scale differences.

4.4. Comparative Experiments

To comprehensively evaluate the performance of the proposed method, comparative experiments were achieved via a comparison with state-of-the-arts (Root-SIFT [25], RIFT [33], CFOG [38], SuperGlue [30]) in terms of precision, RMSE, and MEE. The four methods chosen for comparison are excellent algorithms for image matching based on local features, among which Root-SIFT is a typical method for local feature description based on gradient distribution that has excellent performance in visual place recognition and localization. RIFT and CFOG are the representative works of multimodal remote sensing image registration in the past two years. SuperGlue won IMC-2020, as well as two other competitions at CVPR 2020.

For intuitive evaluation of matching performance, six multimodal remote sensing image pairs with geometric deformation and contrast difference were selected from a remote sensing dataset, a UAV dataset, and an NIR-VIS dataset. Figure 15 shows the matching results of six image pairs by using the five different methods, respectively. The matching result consists of two groups, each with five rows. The first group of images mainly contained contrast differences, subtle structural changes, and rotational changes, while the second group of images contained geometric deformations and contrast differences. The matching results of each group were obtained from top to bottom by Root-SIFT, RIFT, CFOG, SuperGlue, and the proposed method, respectively. For images with large contrast changes from the same sensor (the first pair of the first group), all five methods could obtain good matching point pairs, among which Root-SIFT obtained the least NCM, and RIFT obtained the most NCM. However, CFOG and SuperGlue failed for rotated images from different sensors (the third pair of the first group), mainly because they could not cope with large geometric transformations. Especially for SAR-Optical images with small structural changes and large contrast differences (the second pair of the first group), only the proposed method worked well, which was mainly due to the design of JLFM and the construction of GCID in our method. Due to the lack of ability to handle geometric changes, CFOG could not cope with infrared and visible images containing scale changes (the first pair of the second set), while other methods worked well. For multimodal images with large contrast differences (the second and third pairs of the second group), only RIFT and the proposed method were effective; obviously, our method obtained more NCM. Overall, the proposed method had the best matching performance. The registration results obtained by using the proposed method for six image pairs used above are shown in Figure 16; our method could achieve good registration results for multimodal images with large contrast differences and geometric deformations.

In addition, we selected ten representative images from the four datasets to compare the matching accuracy of the above five related algorithms. The comparative results are given in Figure 17. From the matching accuracy of remote sensing image datasets with contrast differences and inconsistencies in details, we could see that the average precisions obtained by CFOG and Root-SIFT were both below 65%, and those obtained by RIFT and SuperGlue were also below 80%. While the average precision of the proposed method was higher than 85%, because CFOG could not deal with inconsistencies in details and Root-SIFT was sensitive to the contrast differences. Due to the large geometric deformation in the UAV dataset, CFOG and RIFT had much lower average precisions than Root-SIFT and SuperGlue; however, none of these four methods achieved an average precision above 85%, while the proposed method achieved average precision over 90%. On the NIR-VIS dataset with mainly contrast difference, RIFT and the proposed method obtained higher precision than the other three methods. Due to the contrast difference and geometric deformation contained in the computer vision dataset, the matching accuracy of the other four methods fluctuated greatly, and their average precisions was lower than 80%, except that the precision of our method was higher than 85%.

To evaluate the proposed method more comprehensively, we compared the matching precision of the five methods on all image pairs of the four datasets. The comparison results are given in

Table 4. Comparison of Average Precision for five different methods.

METHOD	AVERAGE PRECISION/%
	NIR-IR	UAV	Remote Sensing	Computer Vision
SUPERGLUE	76.46	40.24	44.72	63.31
CFOG	99.17	0	61.50	63.29
RIFT	99.10	12.57	77.68	39.31
ROOT-SIFT	98.15	53.69	41.14	47.98
OURS	99.14	79.31	80.13	81.48

, the performance of the proposed method is significantly better than other methods on UAV, remote sensing, and computer vison datasets. The average precision of our method on the three datasets was 99.14%, 79.31%, 80.13% and 81.48%, respectively. On the NIR-IR dataset only, the average precision of our method was slightly lower than that of CFOG, but still higher than the other three methods.

Table 5. Quantitative evaluation of proposed method on different types of images.

Metric	Cross-Season	Day-Night	Opti-Opti	Depth-Opti	Map-Opti	SAR-Opti	IR-Opti	NIR-Opti	VIS-IR
MEE	1.446	0.9649	1.1478	1.026	1.3842	0.8799	1.1324	1.089	1.50
NCM	99	135	154	204	147	96	183	668	125
RMSE	2.4648	1.9649	1.8412	1.4847	2.378	1.2567	1.3177	1.6547	2.629
Precision	0.7678	0.8508	0.7816	0.8266	0.7595	0.8716	0.9347	0.9914	0.7931

shows the quantitative evaluation results of the proposed method on different types of images. The average precision of the proposed method on different types of multimodal images was higher than 75%; the average NCM was more than 100. The average RMSE of registration on the datasets was less than 3 pixels; the MEE was less than 1.5 pixels, which met the requirements of practical application. Finally, RMSE of registration obtained by the five methods on different types of images are shown in Figure 18. The RMSEs of SuperGlue and Root-SIFT for different types of image data fluctuated greatly; their respective average values were higher than 5 pixels. The RMSEs of CFOG and RIFT were relatively stable; however, their respective average values were still higher than 4 pixels. The average RMSE of the proposed method was less than 3 on each type of images, which indicates that our method was robust to contrast differences and large geometric deformation between the multimodal remote sensing images.

The above verification and comparison experiments show that the proposed method had good adaptability to contrast difference and geometric deformation, which was mainly due to the following two points: (1) Parameter optimization in the process of local structural information extraction; (2) Construction of joint local frequency information map and design of scale and rotation invariance in the feature description process. More matching results of the proposed method on multimodal image pairs are shown in Figure 19, which demonstrates the effectiveness of our method for multimodal image matching with different variations.

In addition, the calculation of our algorithm mainly consisted of three parts: feature detection, feature descriptor construction, feature matching and image registration. The calculation of the feature detection part was mainly the generation of OCFMs. Assume that the maximum resolution of the image is $M\times N$, then computational complexity of the feature detection is $O(M\ast N\ast \log (M\ast N))$. In the feature description part, the computational complexity of GCID is $O(L\ast R^2)$, which is related to the number of extracted feature points $L$ and their description regions $R\times R$. For the last part, the complexity of feature matching is $O(L^2)$ and image registration are $O(M\ast N)$, respectively. Therefore, the total computational complexity of our method is $O(M\ast N\ast \log (M\ast N))+O(L\ast R^2)$. The running time of the proposed registration method on the [email protected] computer was less than 3 s when the resolution of the image was lower than $1000\times 1000$ and the number of feature points was less than 1000.

5. Conclusions

To improve the robustness of multimodal remote sensing image registration with large contrast differences and geometric deformations, a robust local statistical information-based registration framework was developed in this paper. Salient feature points were firstly located according to the phase congruency response map that were optimized by control parameters. Then the geometric and contrast invariant descriptors were constructed based on a joint local frequency information map that combines Log-Gabor filter responses over multiple scales and orientations. Finally, geometric and contrast invariant descriptors are used to match the multimodal remote sensing image pairs and the registration can be achieved by the matching results. Four different multimodal image datasets were used to verify the effectiveness of the proposed method; and the results show that our method was robust to contrast and geometric variations. Through the comparative experimental analysis with the current popular four methods, it is shown that the matching accuracy and registration accuracy of the proposed algorithm are better than the four current popular methods, which confirms the superiority of the algorithm.