Modality-Free Feature Detector and Descriptor for Multimodal Remote Sensing Image Registration

Abstract

The nonlinear radiation distortions (NRD) among multimodal remote sensing images bring enormous challenges to image registration. The traditional feature-based registration methods commonly use the image intensity or gradient information to detect and describe the features that are sensitive to NRD. However, the nonlinear mapping of the corresponding features of the multimodal images often results in failure of the feature matching, as well as the image registration. In this paper, a modality-free multimodal remote sensing image registration method (SRIFT) is proposed for the registration of multimodal remote sensing images, which is invariant to scale, radiation, and rotation. In SRIFT, the nonlinear diffusion scale (NDS) space is first established to construct a multi-scale space. A local orientation and scale phase congruency (LOSPC) algorithm are then used so that the features of the images with NRD are mapped to establish a one-to-one correspondence, to obtain sufficiently stable key points. In the feature description stage, a rotation-invariant coordinate (RIC) system is adopted to build a descriptor, without requiring estimation of the main direction. The experiments undertaken in this study included one set of simulated data experiments and nine groups of experiments with different types of real multimodal remote sensing images with rotation and scale differences (including synthetic aperture radar (SAR)/optical, digital surface model (DSM)/optical, light detection and ranging (LiDAR) intensity/optical, near-infrared (NIR)/optical, short-wave infrared (SWIR)/optical, classification/optical, and map/optical image pairs), to test the proposed algorithm from both quantitative and qualitative aspects. The experimental results showed that the proposed method has strong robustness to NRD, being invariant to scale, radiation, and rotation, and the achieved registration precision was better than that of the state-of-the-art methods.

1. Introduction

Image registration is an essential and fundamental task for remote sensing interpretation. It is aimed at registering images obtained from different sensors, different perspectives, different times or different imaging conditions [1], and is an essential preliminary task for image fusion [2], 3D modeling [3], and change detection [4]. With the rapid advance of remote sensing systems, more and more data sources can now be acquired. As a result, the complementary information between multimodal remote sensing images can significantly improve the capacity and effectiveness of the interpretation. However, the efficiency and accuracy of the image registration result greatly affects the performance of the subsequent processing [1]. Nevertheless, for multimodal remote sensing images, a large number of nonlinear radiation distortions (NRD) will be present, as a result of the different physical imaging mechanisms, which brings significant challenges to the image registration.

Generally speaking, image registration methods can be divided into two categories according to the factors (areas and features) on which they are based. The area-based methods adopt the intensity value of the image itself, while the transformation model for the registration is calculated by optimizing a similarity measure between the image to be registered and the reference image. Correlation [5,6,7], mutual information [8,9,10,11,12,13], or frequency-domain information [14,15] can be used as metrics to measure whether the images are registered. However, these area-based methods often reach a locally optimal solution for the optimization of the model transformation, especially when there is NRD among the images. At the same time, the optimization process for the image conversion parameters is of high computational complexity. The feature-based methods have high robustness to geometric distortion and NRD, so that they are commonly adopted in the registration of multimodal remote sensing images. The feature-based methods achieve the matching goal by identifying the reliable characteristic correspondence between the images, and they are not directly based on the image intensity [16]. The features considered by the feature-based methods include point features [17,18,19], line characteristics [20], and structural features [21,22]. Scale-invariant feature transform (SIFT) [23] is a classic feature point matching method. A number of improved versions of SIFT have since been developed, including affine SIFT (ASIFT) [24], speeded-up robust features (SURF) [25], the SIFT-like algorithm for synthetic aperture radar (SAR) imaging (SAR-SIFT), and principal component analysis SIFT (PCA-SIFT). When SIFT or one of its improved algorithms is used to describe the features, the estimated reference direction must be specified, to make the descriptions more unique and robust to the rotation. Image registration based on features is confronted with two problems: (1) in the feature description stage, the estimation of the principal orientation is often prone to error when it is based on the local features of the image, and a lot of corresponding points will be removed due to incorrect principal orientation estimation; and (2) in the feature matching stage, due to the significant NRD, a feature detection result for an image can often not be found, so that there are many abnormalities in the matching result.

In recent years, research into deep learning has exploded in the computer vision field [26,27]. In the study of medical image registration, deep learning has been used to model the relationship between different modal images [28,29]. However, studies of remote sensing image registration based on deep learning are relatively rare, especially for multimodal remote sensing images [30,31,32,33]. The main reasons for this are as follows: (1) Compared with natural or medical images, remote sensing images have an extensive range, complex distortion, and weak uniqueness of targets. It is also common in remote sensing images that the same ground object type presents different forms, or the same form corresponds to different ground object types. (2) Compared with the large number of natural image datasets that can be used for training, manually annotated remote sensing image datasets are very rare. The labeling of remote sensing images needs considerable expert knowledge and manpower. Furthermore, the application of a model trained on natural images to remote sensing registration tasks is impractical. (3) Deep learning is essentially a method of supervised learning, but it is almost impossible to use a model trained on one modal image to register another modal image. For example, the performance of the LiDAR and optical image registration task can be unsatisfactory when using the SAR and optical image training model. Therefore, there is a need to develop a universal handcrafted descriptor for the multimodal remote sensing image registration task, from the perspective of the physical radiation mechanism and the imaging geometric model, in the case of limited training data.

In this paper, we propose a scale-radiation-rotation-invariant feature transform (SRIFT) algorithm for the registration of multimodal remote sensing images. The contributions of this paper can be summarized as follows:

(1) A modality-free multimodal remote sensing image registration method is proposed, which can handle the scale, radiation, and rotation distortions at the same time. The structural characteristics are captured, in which the same kind of ground object will present a similar structural distribution. Thus, the corresponding features of the images are mapped into a unified space to establish a one-to-one relationship.

(2) The nonlinear diffusion scale (NDS) space is constructed using a nonlinear diffusion filter, instead of a Gaussian filter, to preserve more structure and detail information, which is of great importance for feature extraction for the registration. The structural characteristics are captured by computing the local orientation and scale phase congruency (LOSPC) value in the NDS space of the image. The minimum and maximum moment maps of LOSPC are then used to detect the remaining points. Rotation invariance depends on a rotation-invariant coordinate (RIC) system in the feature description stage. The points in the neighborhood are statistically calculated through a continually changing local coordinate system, which itself realizes rotation invariance, without the need for the estimated orientation to be assigned.

The rest of this paper is organized as follows. The related work is introduced in Section 2. Section 3 details the process of image registration using the SRIFT algorithm. In Section 4, the experimental results obtained using both simulated and real multimodal remote sensing images are provided for the experimental verification and analysis. A summary and our conclusion are presented in Section 5.

2. Related Work

For the registration of homologous images, the pixel basis of the above methods is the intensity or the gradient of the image; however, for multimodal images, the structural information presented in the image by the same point on the ground will be entirely different, due to the NRD caused by the different imaging mechanisms. As a result, the above methods will match many pseudo-corresponding features.

With regard to multimodal remote sensing image registration, scholars have put forward various descriptors on the basis of structural information [34,35]. Compared with gradient information, structural information is less sensitive to nonlinear intensity changes [36]. For example, the edge-oriented histogram (EOH) descriptor designates the shape and contour information of the local image centered on each keypoint, instead of the gradient [37]. The partial intensity invariant feature descriptor (PIIFD) was proposed to solve the problem of the relative gradient direction of the corresponding points [38]. A descriptor for the distribution of internal geometric structures was introduced for images captured on a log-polar grid, which is known as the local self-similarity (LSS) descriptor [39]. A dense LSS (DLSS) method was used to process the registration of optical and SAR images in [40]. The histogram of oriented phase congruency (HOPC) [41] algorithm was extended on the basis of the phase congruency algorithm, and its description process adds phase direction statistics to increase the robustness. Radiation-invariant feature transform (RIFT) [42] also considers the phase congruency, and presents a maximum index map, instead of the gradient, in the feature description. The phase congruency-based structural descriptor (PCSD) [43] has also been successful in optical and SAR image registration without rotational distortion.

In this section, the HOPC and RIFT algorithms are introduced as typical registration methods for multimodal remote sensing images.

The HOPC method successfully introduces an oriented phase congruency algorithm into the automatic registration of remote sensing images, while innovatively mapping the two images acquired under different physical mechanisms into a unified space. Therefore, the corresponding features that cannot be mapped one by one can be mapped in this space. However, there are three significant deficiencies in HOPC: (1) it needs relatively accurate geographic information (or rough geometric correction) for the images in the execution process. However, many multimodal remote sensing images do not contain accurate geographic information. (2) As a template matching algorithm, it is sensitive to geometric distortion, such as rotation and scale. (3) It uses the Harris corner detector, which is very sensitive to NRD when extracting feature keypoints.

The RIFT algorithm was developed on the basis of the limitations of the HOPC algorithm. It also uses a phase congruency algorithm for reference in the stage of mapping the corresponding features, but it adopts a novel descriptor for the feature description. The RIFT algorithm is a method of describing the features per-pixel, which also takes rotation invariance into account. However, its rotation invariance and robustness are not outstanding, because of the possibility of losing the spatial information, and the time-consuming nature of its calculation process. Furthermore, it is not scale-invariant. The main disadvantages of the RIFT algorithm are as follows: (1) the RIFT algorithm adopts a “convolution sequence ring” to deal with the rotation distortion of the images, which may lose some spatial information, resulting in insufficient unique feature vectors being generated and unfavorable feature matching; and (2) the RIFT algorithm does not consider the scale invariance, so it is susceptible to scale changes.

Although image registration has been the subject of numerous studies over the last few decades, there is still no unified registration framework that can automatically register multiple multimodal images while considering the scale and rotation distortions.

3. Image Registration Based on SRIFT

In this section, we describe the proposed SRIFT method in detail. Figure 1 presents a registration flowchart based on SRIFT. The highlighted feature extraction and matching parts in the second column of the figure represent the main innovations of this paper. The first two steps of the algorithm involve solving the problem of feature extraction, and the last step involves the feature description. Initially, the NDS space is constructed using a nonlinear diffusion filter, instead of a Gaussian filter, to preserve more structural and detail information, in order to achieve scale invariance. The structural characteristics are then captured by computing the LOSPC values, in which the same kind of ground object will present similar structural distributions. Finally, the RIC system is applied, which itself realizes rotation invariance, without the need for the estimated reference orientation.

The fundamental reason that restricts the accurate registration of multiple source images is the nonlinear mapping of the corresponding features caused by the existence of NRD. Therefore, as long as the corresponding features can correspond one-to-one, the registration problem for the different source images can be transformed into a registration problem for the same source images. Hence, in the third column of Figure 1, we can use a large number of methods that have been well studied by predecessors in the field of homologous image registration for the image transformation and resampling.

3.1. Scale Invariance Through the Nonlinear Diffusion Scale (NDS) Space

In the construction process of the scale space, the SIFT algorithm uses a Gaussian space, which is generated by convolving the original image with a Gaussian filter of different scales. However, some structure and contour information in the image will be lost in the filtering process because the Gaussian filtering is a kind of smoothing operator, and the loss of information will adversely affect the feature extraction of multimodal remote sensing images with NRD.

The detailed structural information needs to be included as much as possible when the scale space is established. Inspired by anisotropic diffusion [44], a nonlinear diffusion function, instead of a Gaussian function, is adopted to generate the scale space of the image in the proposed algorithm:

$$\frac{\partial f\left(x,y\right)}{\partial t}=f_t=div\left(c\left(x,y,t\right)\nabla f\right)=c\left(x,y,t\right)\Delta f+\nabla c•\nabla f$$

(1)

where $t$ is a scale parameter, $div$ is a bifurcation operator, $\nabla$ is the gradient operator,$\Delta$ is the Laplace operator, and $c\left(x,y,t\right)$ is the diffusion coefficient. In the particular case where the nonlinear diffusion is assumed to be isotropic, $c\left(x,y,t\right)$ is a constant, and the above formula is the same as a Gaussian function. Since there is no analytical solution to the nonlinear diffusion equation, a numerical method is needed to approximate the solution. Equation (1) generates the following relationship after being applied as an additive operator splitting strategy:

$$f^{k+1}={\left(I-\tau {\displaystyle \sum _{l=1}^m{A}_l}(f^k)\right)}^{-1}{f}^k$$

(2)

where $I$ is an identity matrix, $\tau$ signifies the time step, and $l$ represents the direction. Along the l-th coordinate axis, matrix $A_l$ is established accordingly.

The same tactic is applied as is used in SIFT, where the scale space is discretized into a series of $O$ octaves and $S$ sublevels. By using the original image as an initial condition, the multi-scale space for the multimodal image is generated as a series of smoothed images. The scale values are equal to:

$$s={\sigma }^2/2$$

(3)

where $\sigma$ values are calculated from the following expression:

$${\sigma }_i\left(o,s\right)={\sigma }_o{2}^{o+\frac{s}{S}},o\in \left[0,\dots ,O-1\right],s\in \left[0,\dots S+2\right],i\in \left[0,\dots W-1\right]$$

(4)

where ${\sigma }_0$ is the base scale level; $o$ and $s$ are the indices of octave $O$ and sublevel $S$, respectively; and $W$ is the total number of smoothed images. It is notable that the image is downsampled when the last sublevel is reached in each octave, and the downsampled image is used as an initial image for the next octave, as shown in Figure 2.

Through the construction of the NDS space, the attained multi-scale images can preserve the structural and detail information. Therefore, it is anticipated that the proposed SRIFT method will be able to detect many more keypoints.

3.2. Radiation Invariance Through Local Orientation and Scale Phase Congruency (LOSPC)

It should be noted that the orientation and scale referred to in this section are different from the orientation referred to in Section 3.1 and the scale referred to in Section 3.3. The orientation and scale referred to in this section are aimed at giving the spatial mapping more structural information, and the orientation referred to in Section 3.1 and the scale referred to in Section 3.3 are taken into consideration to make the feature description more stable.

3.2.1. Frequency Domain Spatial Mapping via Phase Congruency

The feature extraction can be carried out after the establishment of the multi-scale space. To solve the NRD problem, the LOSPC algorithm is proposed to construct a unified and describable space, which is a prerequisite for the subsequent extraction of the corresponding features.

Feature extraction in the spatial domain of multimodal images often fails due to the distortion of the grayscale and gradient information. Instead, in the frequency domain, an image is decomposed into amplitude and phase components, where the same kind of ground object will present similar structural distribution features in the multimodal images. Phase congruency can be used to measure the degree of local phase information consistency at various angles [45]. Instead of considering the locations with the maximum intensity gradient as being the edges, the phase congruency model regards the edge points as being where the Fourier components are maximally in phase [46]. The phase information is the measurement describing the structural distribution features of the image in the frequency domain.

A 2-D phase congruency operator is developed for the calculation of the phase congruency of any point in the plane, which is a theory that can also be applied to image processing [47]. We define an image as $I(x,y)$, and then the odd-symmetric part $O_{so}(x,y)$ and even-symmetric part $E_{so}(x,y)$ can be obtained by convolving the image $I(x,y)$ with the log-Gabor wavelet transform:

$$\begin{array}{l}\left[E_{so}(x,y),O_{so}(x,y)\right]=\\ \left[I(x,y)\otimes L^{even}(x,y,s,o),I(x,y)\otimes L^{odd}(x,y,s,o)\right]\end{array}$$

(5)

where, in scale $o$ and orientation $s$, $L^{even}(x,y,s,o)$ and $L^{odd}(x,y,s,o)$ stand for the even-symmetric and the odd-symmetric log-Gabor wavelets, respectively. The amplitude and the phase parts of image $I(x,y)$ can be expressed as:

$$A_{so}(x,y)=\sqrt{E_{so}{(x,y)}^2+O_{so}{(x,y)}^2}$$

(6)

$${\varphi }_{so}(x,y)=\arctan (O_{so}(x,y)/E_{so}(x,y))$$

(7)

When all the scales $o$ and orientations $s$ are considered, the results of the two-dimensional phase congruency are calculated as follows:

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

(8)

where $w_o(x,y)$ is a weight function; $\xi$ is a constant with a minimal number; the $\lfloor \rfloor$ action is to prevent negative values, which means that when the value is negative, its result is 0; and $A_{so}(x,y)\Delta {\mathsf{\Phi }}_{so}(x,y)$ is a phase deviation function, which is defined as:

$$\begin{array}{l}{A}_{so}(x,y)\Delta {\mathsf{\Phi }}_{so}(x,y)=(E_{so}(x,y){\overline{\varphi }}_E(x,y)\\ +O_{so}(x,y){\overline{\varphi }}_O(x,y))-\left|E_{so}(x,y){\overline{\varphi }}_O(x,y)+O_{so}(x,y){\overline{\varphi }}_E(x,y)\right|\end{array}$$

(9)

where

$${\overline{\varphi }}_E(x,y)={\displaystyle \sum _s{\displaystyle \sum _o{E}_{so}(x,y)/C(x,y)}}$$

(10)

$${\overline{\varphi }}_O(x,y)={\displaystyle \sum _s{\displaystyle \sum _o{O}_{so}(x,y)/C(x,y)}}$$

(11)

$$C(x,y)=\sqrt{{({\displaystyle \sum _s{\displaystyle \sum _o{E}_{so}(x,y)}})}^2,{({\displaystyle \sum _s{\displaystyle \sum _o{O}_{so}(x,y)}})}^2}$$

(12)

Through the above formulas, each pixel in the image can acquire a statistical value for the phase congruency, which contains orientation and scale information and is based on the structural distribution. In different multimodal remote sensing images, the same ground object will have the same structural distribution. Although it has undergone different physical radiation mechanisms, the value of the phase congruency will be the same.

3.2.2. Feature Point Extraction

In the feature extraction stage, it is necessary to design a feature extraction method based on the statistical index of the phase congruency. Equation (8) is used to calculate the phase congruency of the image pixel by pixel, so an edge graph can be obtained, which is robust to the various multimodal remote sensing images. However, when the image is calculated with log-Gabor filters in different directions, the phase congruency should change with the direction of the filter. Unfortunately, this information is not recorded. Therefore, in order to prevent the phase congruency information changing with the direction, it is necessary to calculate the moment of the phase congruency in each direction, and to record the values with the change of direction. In moment statistics, the principal axis indicates that the moment at this axis is the smallest [46]. The moment perpendicular to the principal axis is the maximum moment. If the maximum moment is large, it is likely to be an edge point in the image, while if the minimum moment is large, it is likely to be a corner point in the image.

According to the moment analysis algorithm, the following three statistics can be obtained from the typical moment calculation formula:

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

(13)

$$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)}$$

(14)

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

(15)

The angle of the principal axis $\psi$ can then be calculated by the following formula:

$$\psi =\frac{1}{2}\arctan \left(\frac{b}{a-c}\right)$$

(16)

After obtaining the principal axis, the minimum moment $m_{\psi }$ and maximum moment $M_{\psi }$ can are calculated as follows:

$$m_{\psi }=\frac{1}{2}\left(c+a-\sqrt{b^2{\left(a-c\right)}^2}\right)$$

(17)

$$M_{\psi }=\frac{1}{2}\left(c+a+\sqrt{b^2+{\left(a-c\right)}^2}\right)$$

(18)

By detecting the extreme values of the maximum moment and the minimum moment on the image, a group of feature points can be obtained, which are called keypoints. The feature extraction step has now been completed. The detection of feature points is undertaken in the frequency domain, instead of the traditional approach of the feature points being detected directly from the image intensity or gradient value, so that the proposed method can better deal with the NRD in multimodal remote sensing images. These keypoints are then used in the subsequent feature matching.

3.3. Rotation Invariance Through a Rotation-Invariant Coordinate (RIC) System

After extracting a large number of stable feature keypoints, these points then need to be described. The process of description should consider the change of features with the transformation of the various influencing factors, and should highlight the uniqueness of the features, to ensure the uniqueness of the subsequent feature matching. In order to achieve rotation invariance in the feature description stage, inspired by the work of “aggregating gradient distributions into intensity orders” [48], a RIC system is proposed. In the feature description, the sample points neighboring a keypoint are statistically calculated through a constantly changing local coordinate system, which itself realizes rotation invariance, without the need for the dominant direction. Firstly, we select several candidate support regions, according to a certain proportion. The difference between local orientation and scale phase congruency (DLOSPC) histogram is then calculated in the local RIC system at each sub-region. Finally, by connecting each DLOSPC vector in the image neighborhood of the multiple support regions, the descriptor is constructed. A flowchart of the construction of the descriptor is shown in Figure 3.

3.3.1. The Local Rotation-Invariant Coordinate (RIC) System

In order to realize the rotation invariance of the proposed algorithm, the descriptor of each support region is calculated by the local RIC system. Specifically, a RIC system is built around each keypoint. $p_i$ is a point in one support region of the keypoint $p$, where the line connecting $p$ and $p_i$ is set as the y-axis, and the direction of the vector $\stackrel{\rightharpoonup }{p{p}_i}$ is the y-axis direction. We then construct a local Cartesian (x-y) coordinate system. For the sample points $p_i$, the pixels in the field participate in the calculation in a rotation-invariant manner; that is to say, the local structure in the field of the sample points is retained. Therefore, the feature description in the locally invariant coordinate system is rotation-invariant. A local RIC system for keypoint $p$ is then constructed. $p_i$ is set as the origin, i.e., the first pixel along the direction of the y-axis is set to $p_{i1}$, and then the pixels in the eight fields of $p_i$ are marked as $p_{i2}$, $p_{i3}$,... $p_{i8}$, as shown in Figure 3b.

3.3.2. The Difference between Local Orientation and Scale Phase Congruency (DLOSPC) in the RIC System

For each sample point $p_i$, its difference of local orientation scale phase congruency can be computed in the local RIC system. The calculation formula is as follows:

$$Dx\left(p_i\right)=I\left(p_{i3}^{}\right)-I\left(p_{i7}^{}\right),$$

(19)

$$Dy\left(p_i\right)=I\left(p_{i1}^{}\right)-I\left(p_{i5}^{}\right),$$

(20)

where $p_{ij}$ are point $p_i$’s neighboring points along the x-axis and y-axis in the local x-y coordinate system, and $I\left(p_{ij}\right)$ stands for the intensity at $p_{ij}$ on the LOSPC map. The difference $D\left(p_i\right)$ and orientation $\theta \left(p_i\right)$ can then be computed as:

$$D\left(p_i\right)=\sqrt{Dx{\left(p_i\right)}^2+Dy{\left(p_i\right)}^2},$$

(21)

$$\theta \left(p_i\right)={\tan }^{-1}\left(Dy\left(p_i\right)/Dx\left(p_i\right)\right).$$

(22)

Note that $\theta \left(p_i\right)$ is converted into the range of $\left[0,2\pi \right)$, along with the values of $Dx\left(p_i\right)$ and $Dy\left(p_i\right)$. The gradient of $p_i$ is then constructed as a d-dimensional vector represented as $F_G\left(p_i\right)=\left(f_{G1},f_{G2},\dots ,f_{Gd}\right)$. To do this, $\left[0,2\pi \right)$ is divided into $d$ equivalent boxes as $di{r}_i=\left(2\pi /d\right)\times \left(i-1\right),i=1,2,\dots ,d$, and then $\theta \left(p_i\right)$ is assigned to the different boxes by linear distance weighting $D\left(p_i\right)$:

$$f_j^G=\{\begin{array}{l}D\left(p_i\right)\frac{\left(2\pi /d-\alpha \left(\theta \left(p_i\right),di{r}_j\right)\right)}{\frac{2\pi }{d}},\alpha \left(\theta \left(p_i\right),di{r}_j\right)<\frac{2\pi }{d}\\ 0,\mathrm{otherwise}\end{array}$$

(23)

where $\alpha \left(\theta \left(p_i\right),di{r}_j\right)$ is the angle between $\theta \left(p_i\right)$ and $di{r}_j$.

3.3.3. Construction of the Keypoint Descriptors Through Multiple Supported Regions

As shown in Figure 3, it is not sufficient to distinguish correct matches from a large number of wrong matches with a single support region. Furthermore, two non-corresponding keypoints may exhibit similarity in some support regions. However, the two corresponding keypoints should have a similar appearance in all the support regions of different sizes, although there may be some small differences due to positioning errors of the keypoints and area detection. That is, when multiple support regions are used, mismatches can be better handled than when only a single support region is used.

N nested support regions are selected with the different radius $r_i$, centered on a keypoint, as shown in Figure 3a. The minimum support region is defined as $A\in ƛ{\Re }^{2\times 2}$, and then the other support regions can be expressed as $A_i=(1/r_i)A$, where $r_i$ represents the size of the $i$-th support region. $r_i$ is defined as $r_i=1+0.5\times \left(i-1\right)$ in this paper, so that the radius increments of the support regions are equal.

The cumulative vectors are combined to form a vector in each support region, and then the cumulative vectors of the different support regions are connected together to describe the features of the keypoint. $F\left(R_i\right)$ is used to represent the cumulative vector for a support region:

$$F\left(R_i\right)={\displaystyle \sum _{p\in R_i}{F}_G\left(p\right)},$$

(24)

Finally, all the vectors calculated in the N support regions are connected together to form the final descriptor $\left\{F_1,F_2\dots F_N\right\}$.

4. Experiments and Analyses

The performance of the proposed SRIFT method was tested in both simulated and real-data experiments. In the simulated experiments, the robustness of the algorithm to geometric distortion was tested by artificially adding various scale and rotation distortions. In the real-data experiments, the registration performance obtained by the SRIFT method was compared to that of eight state-of-the-art image registration methods (SIFT [23], ASIFT [24], SAR-SIFT [49], PSO-SIFT [50], DLSS [40], HOPC [41], PCSD [43], and RIFT [51]), with different modal image pairs.

4.1. Experimental Description

The aim of the simulated experiments was to analyze the ability of SRIFT to resist different image geometric distortions, including scale, rotation, and combined distortions.

4.1.1. Simulated Dataset Construction and Evaluation

The simulated dataset experiments involved two sets of images, in which the ground truth had been geometrically corrected and manually checked, so that the positioning accuracy was better than one pixel, as shown in Figure 4 and Figure 5. Figure 4a,b show the first case of a SAR and optical image pair from Shanghai, China, which were acquired by the GF-3 and GF-2 remote sensing satellites, respectively. Figure 5a,b show the second case of a SAR and optical image pair from Leshan, Sichuan province, China, which were acquired by the Sentinel-1 and Sentinel-2 remote sensing satellites, respectively.

The experimental parameters were set as follows. The simulated datasets were generated by scale and rotation transforms with regard to the ground-truth images, where the images were resampled by 0.5, 0.3, and 0.25 times for the scale transforms and rotated by 10, 30, and 90 degrees, respectively. Therefore, the geometric distortion parameters of the images were known and could be used to calculate the accuracy of the registration results.

If the sensed image $I_{}$ can be regarded as an initial condition, and transform matrices of the simulated dataset can be denoted as $T$, then image $I_{}$ can be transformed into an image $I_{}^{T_{}}$. When the image $I_{}$ has $m$ keypoints $\left\{p_1^{},p_2^{},\dots ,p_m^{}\right\}$, the corresponding keypoints in the image $I_{}^{T_{}}$ are $\left\{p_1^{T_{}},p_2^{T_{}},\dots ,p_m^{T_{}}\right\}$, and are regarded as the ground truth.

The transform matrices calculated after registration are appropriately denoted as $\tilde{T}$, and the root-mean-square error (RMSE) is used to evaluate the registration accuracy of the simulated data images. The higher the RMSE, the worse the accuracy. The RMSE is defined as:

$$RMSE=\sqrt{\frac{1}{\mathrm{m}}{\displaystyle {\sum }_{i=1}^{\mathrm{m}}{(p_i^{\tilde{T}}-p_i^T)}^2}}$$

(25)

where $(p_i^{\tilde{T}}-p_i^T)$ is the residual error, which is calculated between the transformation parameters $\tilde{T}$ and the transformation parameters $T_{}$. $\tilde{T}$ is solved by the corresponding points after registration, while $T_{}$ is given by the ground truth. The RMSE is measured in pixels. The keypoints in which the residual error is less than two pixels are regarded as correctly matched. The number of correctly matched ($NCM$) corresponding keypoints is an important evaluation metric for image matching.

4.1.2. The Overall Performance Comparison

Qualitative evaluation: the overall performances for the corresponding points obtained by the proposed SRIFT method on the simulated datasets are shown in Figure 6 and Figure 7, in which the best results of the state-of-the-art comparison methods for each distortion case are selected for display. As can be seen, the corresponding points obtained by SRIFT are abundant, evenly distributed, and accurately located, and can thus be used to calculate the image transformation model, and then register the image with the calculated model to obtain the registration result.

From the experimental results, it can be seen that, in the process of keypoint extraction, the SRIFT algorithm can extract points characterized by a uniform distribution and sufficient quantity, but after the feature matching and error elimination, the regions with rich shape and texture, and high structural uniqueness, retain more keypoints. In contrast, flat regions with less texture or regions with less structural uniqueness and more repeated textures have fewer corresponding keypoints. Compared with the registration methods based on intensity or gradient, the SRIFT algorithm is essentially a kind of structural descriptor. The description vector of each keypoint is the geometric statistics in an image patch of a specific size centered on this keypoint. In other words, a SRIFT vector describes all the structural information of an image block centered on the keypoint, with a specific sized region (the region sizes are described in Section 3.3.3) as the radius. When the structural information of two image blocks is highly similar, their center points are the registered keypoints.

Quantitative evaluation:

Table 1. Root-mean-square error (RMSE) comparison for simulated dataset 1 (GF-3 SAR image and GF-2 optical image).

Transform	RMSE (Pixels)
	SIFT	ASIFT	SAR-SIFT	PSO-SIFT	DLSS	HOPC	PCSD	RIFT	SRIFT
Scale
0.5	-	-	3.89	2.16	-	-	1.93	-	1.36
0.3	-	-	-	-	-	-	2.85	-	1.77
0.25	-	-	-	-	-	-	2.22	-	1.69
Rotation
10	-	-	2.98	2.88	1.91	1.87	2.06	1.82	1.62
30	-	-	3.33	3.10	-	-	2.96	1.97	1.82
90	-	-	3.12	2.79	-	-	-	1.92	1.49

and

Table 2. RMSE comparison for simulated dataset 2 (Sentinel-1 SAR image and Sentinel-2 optical image).

Transform	RMSE (Pixels)
	SIFT	ASIFT	SAR-SIFT	PSO-SIFT	DLSS	HOPC	PCSD	RIFT	SRIFT
Scale
0.5	-	-	3.65	1.99	-	-	1.89	-	1.21
0.3	-	-	-	-	-	-	2.23	-	1.43
0.25	-	-	-	-	-	-	1.98	-	1.36
Rotation
10	-	-	2.22	2.14	1.98	1.99	1.96	1.95	1.47
30	-	-	3.84	3.27	-	-	2.66	1.97	1.72
90	-	-	3.01	2.52	-	-	-	1.99	1.36

list the RMSEs of the simulated data image registration results of the SRIFT method, as well as those of the other state-of-the-art image registration methods.

As can be seen from Table 1 and Table 2, the registration accuracy of the SRIFT algorithm is consistently the highest. Compared with SIFT and ASIFT, the SRIFT algorithm has a better feature extraction and description ability for multimodal images. SAR-SIFT, PSO-SIFT, DLSS, HOPC, PCSD, and RIFT are specially designed for multimodal image registration. Compared with the template matching algorithms such as DLSS and HOPC, SRIFT can resist the various scale and rotation distortions. Compared with feature matching algorithms such as SAR-SIFT, PSO-SIFT, PCSD, and RIFT, SRIFT can overcome more complex image geometric distortions. For a more detailed analysis of the results of the state-of-the-art algorithms, see Section 4.2.3.

4.1.3. The Ability of the Algorithm to Resist the Scale and Rotation Distortion

To test the robustness of SRIFT to image scale and rotation distortion, simulated images were generated by resizing the image using a scale change factor from 1 to 10 with different intervals and rotating the image using a rotation change factor from 0 to 360 with an interval of 30. The matching performance of the proposed method with scale and rotation distortion is shown in Figure 8.

As shown in Figure 8a, when the image scale difference is between one time and four times, the precision of the proposed method does not decrease significantly when the scale factor increases. The proposed method maintains a good performance when the scale factor is less than four times. We can, therefore, conclude that the proposed method is robust for scale difference. However, when the scale difference increases to more than five times, the correct matching point number plummets. From Figure 8b, for rotation distortion at any angle, the proposed SRIFT algorithm extracts relatively abundant corresponding points.

4.2. Experiments with Real Images

In the real-data experiments, the registration results obtained by the proposed method were compared to those obtained by eight state-of-the-art methods: SIFT, ASIFT, SAR-SIFT, PSO-SIFT, DLSS, HOPC, PCSD, and RIFT.

4.2.1. Real Datasets

In the real-data experiments, nine sets of multimodal images were selected: SAR/optical, DSM/optical, LiDAR/optical, NIR/optical, SWIR/optical, classification/optical, and map/optical images, which are shown in Figure 9 and described in

Table 3. The real datasets.

No.	Data Source	Image Type	Size	GSD	Date	Image Content	Description
1	GF-3	SAR	4865 × 3504	4 m		Urban	SAR/optical
	GF-2	Visible	3979 × 3619	4 m
2	GF-3	SAR	6000 × 6000	8 m		Mountains and water	SAR/optical
	GF-1	MSS	6000 × 6000	8 m
3	Sentinel-1	SAR	2000 × 2000	10 m	04/11/2017	Mountains and water	SAR/optical
	Sentinel-2	Visible	2000 × 2000	10 m	06/11/2017
4	Manual	DSM	1200 × 1200	1 m	15/05/2017	Urban	DSM/optical
	UAV	Visible	1200 × 1200	1 m	15/05/2017
5	LiDAR	Intensity	349 × 349	2.5 m	22/06/2012	Urban	LiDAR/optical
	Airborne	Hyperspectral	349 × 349	2.5 m	23/06/2012
6	GF-2	NIR	1202 × 1011	3.2 m	11/04/2016	Farmland and water	NIR/optical
	Google Earth	Visible	1014 × 950	4 m
7	Sentinel-2	SWIR	1000 × 1000	20 m		Farmland and water	SWIR/optical
	Sentinel-2	Visible	1000 × 1000	10 m
8	Manual	Classification	640 × 400	4 m	05/03/2013	Urban	Classification/optical
	GF-1	Visible	640 × 400	4 m	05/03/2013
9	Google Earth	Map	1867 × 1018	4 m		Urban	Map/optical
	Google Earth	Visible	1867 × 1018	4 m

.

Figure 9a shows the first SAR/optical image pair covering an urban area. These images were acquired by the GF-3 and GF-2 remote sensing satellites, respectively. The resolution of the SAR image was set to 1 m, referring to the 4-m resolution of the optical image through panchromatic and multispectral fusion. The image sizes are 4865 × 3504 and 3979 × 3619 pixels, respectively.

Figure 9b shows the second SAR/optical image pair covering a mountainous and water area. These images were acquired by the GF-3 and GF-1 remote sensing satellites, respectively. The resolution of the SAR image was set to 8 m, referring to the original 8-m resolution of the multispectral optical image. The image sizes are both 6000 × 6000 pixels.

Figure 9c shows the third SAR/optical image pair, again covering a mountainous and water area. These images were acquired by the Sentinel-1 and Sentinel-2 remote sensing satellites, respectively, in November 2017. The resolution of the SAR image was set to 10 m, referring to the original 10-m resolution of the multispectral optical image. The image sizes are both 2000 × 2000 pixels.

Figure 9d shows the DSM/optical image pair covering an urban area. These images were acquired by manual production and an unmanned aerial vehicle (UAV), respectively, in May 2017. The resolution of the DSM was set to 1 m, referring to the original 1-m resolution of the optical image. The image sizes are both 1200 × 1200 pixels.

Figure 9e shows the LiDAR/optical image pair covering an urban area. These images were acquired by a LiDAR system and aerial photography, respectively, in June 2017. The image resolutions are 2.5 m, and the image sizes are 349 × 349 pixels.

Figure 9f shows the NIR/optical image pair covering a farmland and water area. The NIR image was acquired by the GF-2 remote sensing satellite in April 2016, and the optical image was downloaded from Google Earth. The image resolutions are 3.2 m and 4 m, and the sizes are 1202×1011 and 1014 × 950 pixels, respectively.

Figure 9g shows the SWIR/optical image pair covering a farmland and water area. These images were both acquired by the Sentinel-2 remote sensing satellite. The image resolutions are 20 m and 1 m, respectively, and the image sizes are 1000 × 1000 pixels.

Figure 9h shows the classification/optical image pair covering an urban area. These images were respectively acquired by manual production and the GF-1 remote sensing satellite in March 2013. The resolution of the classification image was set to 4 m, referring to the original 4-m resolution of the multispectral optical image. The image sizes are both 640 × 400 pixels.

Figure 9i shows the map/optical image pair covering an urban area. These images were downloaded from Google Earth. The image resolutions are 4 m, and the image sizes are 1867 × 1018 pixels.

There is severe distortion between all these images, especially radiation distortion. They therefore pose a significant challenge for the image registration algorithms.

4.2.2. Ground-Truth Setting and Evaluation Metrics

Ground truth is essential for the quantitative evaluation of registration. However, due to the different sensors and/or different perspectives in multimodal images, it is often impossible to achieve a one-to-one correspondence between pixels. Therefore, precise geometric correction is required for every set of image pairs in real datasets, and then the geometric correction transformation parameters are approximated as the ground truth. In detail, a certain number of precise corresponding points are artificially selected in each pair of images, and the image transformation parameters are then solved using these points. The images are registered using these parameters, and the registration result is manually checked. If the result is not accurate, the artificial selection is repeated until the two image pairs are accurately registered. When evaluating the registration accuracy, the image transformation parameters estimated by the registration algorithm are used to calculate the residuals of the artificially selected corresponding points.

For the real-data experiments, three evaluation metrics were selected to evaluate the registration performance. The precision is expressed as $Precision=NCM/NM$, where $NM$ is the total number of corresponding keypoints. The definitions of $NCM$ and RMSE have been given in Section 4.1

4.2.3. Registration Performance Comparison

Qualitative comparison: Figure 9 intuitively shows the corresponding point-line diagrams obtained by the SRIFT algorithm in the real-image experiments, where the number and distribution of the corresponding points reflect the robustness and applicability of the algorithm. Figure 10 uses checkerboard mosaicked images of the nine groups for the qualitative evaluation, where the continuousness of the sub-region edges of the images directly reflect the accuracy of the registration.

Through the analysis of the above experimental results, the following conclusions can be drawn. The proposed SRIFT registration algorithm was able to obtain satisfactory results on all nine datasets of various multimodal remote sensing images. When dealing with the multimodal image registration task, the SRIFT algorithm fully considers the influence of NRD in the feature extraction, as well as the feature description, so that it can extract a large number of stable and evenly distributed feature points. The SRIFT algorithm can also resist image scale and rotation distortion, and the registration image has a high coincidence degree with the reference image.

Quantitative comparison: Figure 11 quantitatively reflects the registration effect and accuracy of the eight algorithms on the nine sets of data with the three measurement indices introduced in Section 4.2. Figure 11a is the line chart of $NCM$, where the higher the value of $NCM$, the more keypoints are correctly matched, which reflects the ability of the different algorithms in the feature matching stage. Figure 11b is the line chart of the precision, where the higher the value, the higher the proportion of correct matching points in all the matching points, reflecting the ability of the different algorithms in the feature description stage. Figure 11c is the line chart of the RMSE, where the lower the value, the higher the registration accuracy, reflecting the overall registration ability of the different algorithms.

As can be seen in Figure 11, SRIFT achieves the best precision. RIFT ranks second and PCSD ranks third. The basic idea of the HOPC and DLSS algorithms is to divide the images into blocks by counting the local information of each image block, and to then integrate the blocks into the overall information. Their performance generally lies in the middle level among the eight compared methods. The SIFT and ASIFT algorithms are not designed for multimodal data, so that they perform the worst of all. A detailed analysis of each method is presented in the following.

Due to the SIFT algorithm detecting feature points directly based on the intensity, and using gradient information for the feature description, which is sensitive to NRD, SIFT only obtains good registration results for the SWIR/optical case. In the SWIR/optical case, the sensors are on the same satellite, and the two sources have little difference in radiation mechanism, so that the image registration is relatively easy.

The results of the ASIFT algorithm are slightly better than those of SIFT. The registration results cannot be obtained for most of the images, but the registration results are better than those of SIFT in the complex set of geometric transforms in the NIR/optical case, because ASIFT is specially designed for affine transformation. ASIFT simulates the scale and the camera direction and normalizes the rotation and the translation.

The SAR-SIFT method was specially designed for SAR imagery, and it relies on a new gradient computation method adapted to SAR images. Therefore, the image registration results for the first three SAR datasets are satisfactory. However, the redefined gradient probability has difficulty in dealing with complex radiation distortion, and the multi-scale Harris detector has insufficient resistance to NRD.

PSO-SIFT achieves the best registration effect among the SIFT-related algorithms, which is due to the fact that PSO-SIFT applies multiple constraints, e.g., the feature distance, and hence results in a better registration.

The DLSS algorithm is an improved version of the LSS algorithm, which divides the template window into spatial regions called “cells.” Each cell contains n × n pixels, and has an overlapping region of half a cell width with the neighboring cell. This method of division is essentially a template matching method, rather than a feature matching method. As a result, DLSS cannot resist the complex geometric distortion, and the registration performance in datasets 1, 4, 6, 7, and 9 is poor.

HOPC uses the Harris detector to detect the feature points. However, the Harris detector is very sensitive to NRD, and it is not universally suitable for all the different types of multimodal images. Therefore, overall, its registration effect is slightly worse than that of DLSS. Moreover, HOPC is similar to DLSS in the blocking strategy, so that HOPC also has a poor effect in images that cannot be registered by the DLSS algorithm.

The RIFT algorithm does not have scale invariance, so its registration effect in the SWIR/optical case is inferior. In order to achieve rotation invariance, RIFT transforms the initial layer to reconstruct a set of convolution sequences with different initial layers, and then calculates a maximum index map (MIM) from each convolution sequence to obtain a set of MIMs. This statistical method only establishes the relationship between the neighborhood pixels of the keypoints and the keypoints themselves, and it destroys the structural relationship between the neighborhood pixels of the keypoints. Therefore, for images with both scale and rotation transformation, the registration effect is weak.

In the feature description stage, the PCSD algorithm adopts a method similar to SIFT, which requires estimation of the main direction, and error in the primary direction estimation causes the extracted corresponding points to be deleted by mistake. In the feature extraction stage of PCSD, the points with the closest cosine similarity to the keypoints in the reference image are matched. However, this method does not take advantage of a phase congruency algorithm in the feature extraction, and its stable alignment points are insufficient.

The analysis of the experimental results confirms the powerful registration ability of the SRIFT algorithm. If the structural information of the images to be registered is similar, the SRIFT algorithm can register the images without considering the distortion of the image intensity.

4.2.4. Computational Cost

An experiment addressing the computational cost of the proposed method and that of some of the comparison methods (SIFT, SAR-SIFT, PCSD, HOPC, RIFT, SRIFT) was undertaken. The average run times of the different methods on the real datasets are shown in

Table 4. Computational cost comparisons (seconds).

Method	SIFT	SAR-SIFT	PCSD	HOPC	RIFT	SRIFT
Computational cost	2.04	15.25	46.82	57.63	38.65	64.70

. The experiment is carried out on a computer with an Intel i7-6200U CPU @ 2.30GHz and 8 GB of RAM. All the methods were implemented in MATLAB.

According to Table 4, the SIFT algorithm has the highest computational efficiency, followed by the SAR-SIFT algorithm, because neither algorithm performs phase congruency calculation. Among the four methods that adopt phase congruency (PCSD, HOPC, RIFT, and SRIFT), RIFT shows the highest computational efficiency because it replaces phase congruency with maximum index graph. PCSD and HOPC are in the middle level, but the calculation time of the proposed method is relatively high. Due to the use of the multi-scale space, phase congruency, and RIC system to deal with the NRD, as well as the complex scale and rotation variation, the proposed method is relatively time-consuming. Algorithm optimization and efficiency improvement will be carried out in the future.

5. Conclusions

In this paper, we have presented a modality-free multimodal remote sensing image registration method named SRIFT, which has the advantage of scale, radiation, and rotation invariance, making it suitable for use with different multimodal remote sensing images. The building of a robust NDS space, the definition of a new concept called LOSPC, and the development of a new RIC system are the three main contributions of the proposed SRIFT method. The NDS space is constructed to resist the scale distortion of image pairs with a large difference in gradient distribution. LOSPC is computed in the NDS space of the images, in which the same kind of ground object will present similar structural distributions, and thus the features of these images are mapped into the same space. The idea of the RIC system is that the points in the neighborhood are statistically calculated through a continually changing local coordinate system, which is more suitable for the feature matching task than a global coordinate system, and realizes rotation invariance, without the need for the estimated orientation to be assigned. In the experimental analysis, two simulated datasets and nine sets of real data were used to qualitatively and quantitatively compare the registration performance of SIFT, ASIFT, SAR-SIFT, PSO-SIFT, DLSS, HOPC, PCSD, RIFT, and the proposed SRIFT method. The registration performance of the SRIFT algorithm on the multimodal images with NRD was superior to that of the other state-of-the-art image registration methods. Our future study will focus on research into a correction model and error elimination for multimodal image registration.