PDE-Constrained Scale Optimization Selection for Feature Detection in Remote Sensing Image Matching

Abstract

Feature detection and matching is the key technique for remote sensing image processing and related applications. In this paper, a PDE-constrained optimization model is proposed to determine the scale levels advantageous for feature detection. A variance estimation technique is introduced to treat the observation optical images polluted by additive zero-mean Gaussian noise and determine the parameter of a nonlinear scale space governed by the partial differential equation. Additive Operator Splitting is applied to efficiently solve the PDE constraint, and an iterative algorithm is proposed to approximate the optimal subset of the original scale level set. The selected levels are distributed more uniformly in the total variation sense and helpful for generating more accurate and robust feature points. The experimental results show that the proposed method can achieve about a 30% improvement in the number of correct matches with only a small increase in time cost.

1. Introduction

Remote sensing refers to the detection of targets through various sensors mounted on different aircraft, and the obtained remote sensing image data often carry a large amount of information, requiring the use of different efficient processing techniques for practical purposes [1,2]. Image matching is an important technical foundation in remote sensing image processing, and the effectiveness of registration often directly affects the subsequent image processing results [3,4,5,6].

The work related to image matching can be traced back to P.E. Anuta’s work on cross correlation image detection [3] and Brown’s summary of previous work [4], which has since produced numerous influential works. For image pairs with significant deformation, Flusher proposed an adaptive mapping method that segments the image into multiple sub-blocks, calculates the similarity between the sub-blocks, matches them using the similarity of the sub-blocks, and registers the original image using the coordinate relationship between the sub-blocks [7]. Studholme proposed a normalized mutual information method that overcomes the sensitivity of traditional mutual information methods to the size of overlapping parts between images and achieves good matching results [8]. David Lowe proposed the classic Scale Invariant Feature Transform (SIFT) by combining the scale space pyramid and Gaussian kernel to improve the features invariance and detection rate of feature points [9,10]. The SIFT feature maintains the same changes in image grayscale, angle, scale, etc., and performs well in natural image matching. Some techniques based on SIFT have been developed over the past years. Ke and Sukthankar introduced principal component analysis to replace the weighted histograms and reduce the dimension of features to 20 [11]. Speeded-up robust features (SURFs) were presented based on a Hessian matrix and image convolutions to faster the computing [12]. Dellinger et al. proposed an improved SIFT registration algorithm for Synthetic Aperture Radar (SAR) images called SAR-SIFT [13]. This algorithm uses a new gradient calculation method to effectively suppress the noise impact of SAR images. Alcantarilla et al. proposed a novel multiscale 2D feature detection and description algorithm in nonlinear scale spaces, named KAZE [14,15]. The results revealed a step forward in performance both in detection and description compared to some traditional methods. Gong et al. utilized two matching strategies, feature matching and regional matching, and then proposed a rough-to-accurate registration method, which first uses the SIFT algorithm to coarsely register the image and then uses mutual information methods to evaluate the initial matching points and find accurate matching point pairs [16]. Wu et al. improved the Random Sample Consistency (RANSAC) algorithm and proposed using a fast iterative strategy to filter feature point pairs to obtain more accurate transformation parameters in order to achieve the accurate matching of remote sensing images [17]. In order to effectively solve the problem of significant changes in image grayscale leading to the failure of traditional registration algorithms, Ma et al. proposed a Position Scale and Orientation Scale Invariant Feature Transform (PSO-SIFT) algorithm based on the SIFT algorithm, which uses a new gradient definition to calculate feature point descriptors [6]. Accurate matching between remote sensing image pairs with significant grayscale changes has been achieved. Wu et al. proposed the Particle Swarm Optimization Sample Consistency (PSOSAC) algorithm, which uses particle swarm optimization to optimize the model parameters solved using RANSAC based on the RANSAC algorithm to obtain more accurate matching results [18]. There are some other techniques that have been also introduced to achieve better performances [19,20,21,22,23,24].

Most of these researches are based on a Gaussian Scale Space and have shown success in many problems. However, it is still to be further explored in more complex scenes, such as heterogeneous image matching, multi-modal image matching, and medical image matching [25,26,27]. A remote sensing image is much different from a nature image, not only in acquisition and transport but also in information presentation and vision characteristics [1,28,29]. It is necessary to develop more suitable and efficient feature describing methods for remote sensing image matching. Scale selection techniques can be found in various image processing tasks. Olson presented a method where the scale of the filtering and feature detection is varied locally according to the distance to the scene pixel estimated through stereoscopy [30]. Energy models inspired from Markov random fields and segmentation quality measures can be introduced to determine scales for different purposes [31,32]. PDE-constrained optimization has generated a great deal of attention in many fields [33,34]. In related problems, a certain functional is aimed to be optimized with one or more PDEs serving as constraints.

This paper aimed to present a novel scale optimization selection technique for local feature describing in remote sensing image matching. A nonlinear partial differential equation is configured with the noise estimation results to generate a continuous nonlinear scale space. Then, a PDE-constrained model is proposed to determine the optimized scale levels. It can be numerically solved based on Additive Operator Splitting (AOS) and subset selection. After some calculation steps similar to those of KAZE and SIFT are performed on the scale space, the feature points can be found for remote sensing image matching.

The rest of this paper is arranged as follows. In Section 2, a linear scale space and nonlinear scale space are introduced. In Section 3, an evolution partial differential equation configured with noise estimation is applied to describe a nonlinear scale space. And then, a PDE-constrained model is presented to optimally select the scale levels. In Section 4, the several experiments implemented in this study to verify the accuracy and efficiency of proposed algorithm are discussed.

2. Fundamentals and Basis

2.1. Linear Scale Space

Scale space theory has been introduced in many models, and it can help to distinguish the features at different levels. After a scale parameter is added to the related model, a multi-scale representation sequence can be computed following a change in the parameter. Many tasks can be achieved by analyzing or processing the sequence, such as extracting features.

The scale space method can simulate the process of people observing objects from near to far, and the continuous change in scale can effectively describe the essential features of the image. The Scale Invariant Feature Transform (SIFT) uses a Gaussian pyramid to represent the linear scale space, which means to convolve a serial of Gaussian kernel functions with the original image function to obtain smoothed images at different scales [10,35]. It can be formulated as

$$\begin{array}{c}\hfill I(x,t)=I\left(x\right)\ast g_t\left(x\right)=I\left(x\right)\ast \frac{1}{4\pi t^2}exp\left(-\frac{{\left|x\right|}^2}{4t^2}\right),(x,t)\in \Omega \times \Gamma .\end{array}$$

(1)

Here, t denotes the scale parameter, and x denotes the pixel position. $\Omega \subset R^2$ denotes the image domain, and $\Gamma$ denotes the time interval. $I\left(x\right)$ is the original image, and $I(x,t)$ means the related linear scale space. The value of t represents the degree of smoothness of the original image. The larger the value, the higher the degree of smoothness, and then, different scale image features can be obtained through some specific steps, which are the following: calculating the Difference of Gaussian (DoG), searching for extreme points, determining the key points (interesting points) and the main direction, and generating feature descriptors.

The discrete difference of Gaussian space can be denoted as

$$\begin{array}{c}\hfill DoG(x,t)=I(x,t+\Delta t)-I(x,t).\end{array}$$

(2)

Here, $\Delta t$ can be preset or determined according to the image data. The first two rows in Figure 1 show different Gaussian kernels and the corresponding smoothed images. The last row shows the difference between two neighbor images. With the scale parameter t increasing, the original image will be gradually smoothed, and the feature information will be redundant. A downsampling technique is adopted in SIFT to obtain a more compact version (small-sized) for later calculations.

Based on the scale space and the difference space, the original SIFT can be summarized as two main steps (see [10,36]).

2.1.1. Calculating the DoG and Determining the Key Points

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle K=\left\{x\in P,\left|-{\left(\frac{{\partial }^{2}DoG}{\partial x^2}\right)}^{-1}\frac{\partial DoG}{\partial x}\right|<\frac{1}{2}\wedge \left|DoG\left(x\right)\right|<0.03\wedge \frac{Tr{\left(H\right)}^2}{\left|H\right|}<\frac{{(r+1)}^2}{r}\right\}.}\hfill \end{array}\end{array}$$

(3)

Here, H means the Hessian of the DoG, and $P=\left\{{argmax}_{z\in B}DoG\left(z\right)\right\}$ denotes a set of extreme points found in some basis subareas of the image domain. For example, a subarea can be set to a circle-centered pixel x with radius d, that is, $B(x,d)$. The second derivative of the DoG is actually the Hessian in a 3D space (spatial 2D and scale 1D), and it is different from the Hessian in the 2D spatial space. ∧ denotes the Boolean operator AND. $Tr$ denotes the trace operator. r is a threshold used to determine good key points.

2.1.2. Orientation Assignment and Generating Feature Descriptors

$$\begin{array}{c}\hfill Feat{s}_{SIFT}\left(x\right)={\left[\mathrm{histogram}\left(\nabla R_{\theta }I\left(B_{i,j}\right),{\theta }_d\right)\right]}_{4\times 4\times 8}\end{array}$$

(4)

$$\begin{array}{c}\hfill \theta \left(B\right)=\underset{{\theta }_d}{argmin}\mathrm{histogram}(\nabla B\left(x\right),{\theta }_d).\end{array}$$

Here, ${\theta }_d,d=1,2,\cdots ,8$ means the preset eight direction angles. $B_{i,j}$ is a subdomain of $B\left(x\right)$, which is the neighborhood centered keypoint x. $R_{\theta }$ denotes the rotate matrix with angle $\theta$. ∇ denotes the gradient operator. The image is divided into $4\times 4$ sub-regions, and eight direction statistics are calculated on each corresponding gradient field. Then, the statistics on 16 sub-regions are collected to form a feature vector with dimensions $4\times 4\times 8$.

2.2. Nonlinear Scale Space

Stable feature points at different scales can be extracted on the previously established linear scale space. However, while smoothing noise, Gaussian convolution also smooths image details and edges, affecting the accuracy and uniqueness of feature points. Some methods have been introduced to improve this defect in the past years and achieved good results in applications [11,14,19,37,38]. The KAZE features are extracted from a nonlinear scale space that is generated using a related partial differential equation. And it has been validated to have a comparable efficiency on some problems in image denoising and restoration. In general, a nonlinear scale space can be denoted as the solution of a nonlinear partial differential equation:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\partial I}{\partial t}(x,t)=\mathrm{div}\left(c\left(|\nabla I|\right)·\nabla I\right)}\hfill \end{array}\end{array}$$

(5)

Here, $c\left(|\nabla I|\right)$ means the diffusion velocity, and the equation will degenerate into a linear version if c degenerates into a constant. Unlike the linear scale space, there is no explicit formulation of the nonlinear scale space but an implicit formulation described using the above Equation (5).

It is clear that the continuous nonlinear scale space is determined using the diffusion velocity and evolution time. There are lots of different configurations for them. For the velocity, there are several useful formulations, such as the following [14]:

$$\begin{array}{c}\hfill c_1\left(s\right)=exp\left(-\frac{s^2}{b^2}\right),c_2\left(s\right)=\frac{1}{1+\frac{s^2}{b^2}},c_3\left(s\right)=\left\{\begin{array}{cc}{\displaystyle 1,}\hfill & s=0\hfill \\ {\displaystyle 1-exp\left(-\frac{3.315}{s^8/b^8}\right),}\hfill & s>0\hfill \end{array}\right.\end{array}$$

(6)

As there is no analytical solution for Equation (5), it is necessary to introduce some numerical technique such as Additive Operator Splitting (AOS) to approximate the solution at a different scale t. The discretization of the equation can be denoted as

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

(7)

and an efficient semi-implicit scheme can be obtained as

$$I^{n+1}={\left[1-\tau \sum _{l=1}^m{A}_l\left(I^n\right)\right]}^{-1}{I}^n.$$

(8)

Here, $A_l$ is a matrix that encodes the image conductivities for each dimension, and this scheme is advantageous for large time steps. After taking some steps similar to those in SIFT, the image features can be generated.

Figure 2 shows three nonlinear scale spaces based on the velocity $c_2\left(s\right)$ with different values of k. It can be found that the total variation decreased at different rates. And if such a space is divided according to a fixed mode (such as uniformly or proportionally in the scale parameter), the differences between two adjacent parts are various and difficult to be balanced. It is not advantageous to detect the image feature points and achieve a satisfactory performance in image matching.

3. Proposed Method

It can be found that the continuous nonlinear scale space $I(x,t)$ is absolutely determined through the velocity and evolution time. Different diffusion velocities are very different in smoothing or preserving the potential features, so a proper velocity should be selected for each image. Of course, the noise level is also an important factor affecting the determination. The number of scale levels is another important parameter affecting the feature point detection. It is easy to know that if more scale levels are set, a higher computation cost is required though less feature points may be lost. Even if a proper number of levels is set, it is still difficult to determine which levels should be set. Though some experience and recommend selection methods can be taken from previous works, we still tried to present an optimization model to solve these difficulties.

The surroundings in remote sensing images is a complex variable, and the signal is susceptible to being polluted. There are also many changes in the type and intensity of image features; a flexible pattern of an alternating scale space is advantageous for the detection of more characteristics, which benefits later applications.

It is important to know some information about noise and feature intensity before the parameters of the scale space are determined. Here, we introduce a filter to compute a low-density part of the image [39]:

$$\begin{array}{c}\hfill \widehat{n}=I\ast h=I\ast \left[\begin{array}{ccc}\hfill 1& \hfill -2& \hfill 1\\ \hfill -2& \hfill 4& \hfill -2\\ \hfill 1& \hfill -2& \hfill 1\end{array}\right]\end{array}$$

(9)

Then, the noise mean, noise variance, and Signal-to-Noise Ratio (SNR) can be estimated by

$$\begin{array}{c}\hfill {\displaystyle {\mu }_{est}=\mathrm{Mean}\left(\widehat{n}\right),{\sigma }_{est}^2=\mathrm{Var}\left(\widehat{n}\right),SN{R}_{est}=10{log}_{10}\frac{\sum {[I-\widehat{n}]}^2}{\sum {\widehat{n}}^2}.}\end{array}$$

(10)

Based on these, a new parameter configuration is formulated as (11), for Equation (5) to generate a continuous nonlinear scale space:

$$b=\frac{{\mu }_{est}·SN{R}_0}{SN{R}_{est}}.$$

(11)

Here, $SN{R}_0$ means a reference SNR value that can be preset using some prior knowledge of the image data.

However, it is still difficult to determine which scale levels should be selected to compute the feature descriptor for the interesting points. In this paper, it is assumed that the interesting points are uniformly distributed in the scale space, so we tried to divide the scale space into several parts in which the changed amount of some metric is equal. Then, a related PDE-constrained optimization problem can be expressed as

$$\begin{array}{c}\hfill (t_1,t_2,\cdots ,t_N)=\underset{t_1,t_2,\cdots ,t_N}{argmin}\sum _{p=0}^N{\left[\mathrm{Ene}(I(x,t_{p+1})-\mathrm{Ene}\left(I(x,t_p)\right)\right]}^2\end{array}$$

(12)

$$\begin{array}{c}\hfill s.t.\left\{\begin{array}{ccc}{\displaystyle \frac{\partial I}{\partial t}(x,t)}\hfill & =\mathrm{div}\left(\frac{\nabla I}{1+{\left(\frac{SN{R}_{est}}{{\mu }_{est}\ast SN{R}_0}\right)}^2{|\nabla I|}^2}\right),\hfill & (x,t)\in \Omega \times (0,T)\hfill \\ {\displaystyle I(x,0)}\hfill & =I_0\left(x\right),\hfill & x\in \Omega ,\hfill \\ {\displaystyle \frac{\partial I}{\partial \mathbf{n}}(x,t)}\hfill & =0,\hfill & x\in \partial \Omega .\hfill \end{array}\right.\end{array}$$

(13)

Here, $\mathrm{Ene}\left(I\right)={\left|\right|\nabla I\left|\right|}_2^2$ denotes the image energy needed to measure the variant of $I(x,t)$. The parameters $t_0=0$ and $t_{N+1}=T$ are preset. The initial condition means that the scale space should be generated from the observed image $I_0$. div denotes the divergence operator. This model tries to determine N scale levels $t_1,t_2,\cdots ,t_N$ that uniformly divide the scale space $I(x,t),t\in [t_0,t_{N+1}]$.

To find the optimal solution, the continuous scale space $I(x,t)$ must be solved from Equation (13), which requires a fair amount of calculation. To reduce the cost for solving the whole scale space, a discrete version is introduced to approximate the scale space, and then the model (12) can be rewritten as

$$\begin{array}{cc}\hfill (k_1,k_2,\cdots ,k_N)=& \underset{k_1,k_2,\cdots ,k_N}{argmin}{\sum }_{p=0}^N{\left[\mathrm{Ene}\left(I^{k_{p+1}}\right)-\mathrm{Ene}\left(I^{k_p}\right)\right]}^2\end{array}$$

(14)

$$\begin{array}{cc}\hfill s.t.\hspace{1em}\hspace{1em}& I^{n+1}={\left[1-\tau \sum _{l=1}^m{A}_l\left(I^n\right)\right]}^{-1}{I}^n\end{array}$$

(15)

where $\tau$ denotes a small enough time step, and (15) means to solve Equation (13) solved using the AOS scheme.

Attributed to the reference step size $\tau$, the original optimized scale levels $t_1,t_2,\cdots ,t_N$ are approximated by $k_1\tau ,k_2\tau ,\cdots ,k_N\tau$. The variables to be solved are changed from continuous type to integer type, and then the cost will be greatly reduced.

However, it is still not easy to properly choose N levels from the level set with about $\frac{T}{\tau }$ elements, and this is similar to a subset selection problem to some extent. To efficiently solve such a problem, an iterative approximation is introduced to the proposed algorithm, Algorithm 1.

After setting the initial image $I^0=I(x,0)$, maximal scale level T, maximal allowed error $\epsilon$, maximal iteration number $maxIter$, initial solution $k_1^0,k_2^0,\cdots ,k_N^0$ and then applying Algorithm 1, optimal scale levels $k_1,\cdots ,k_N$ will be determined, and the nonlinear scale space can be uniformly divided into several parts in terms of total variation. As shown in Figure 3, the nonlinear scale space of a remote sensing image is divided by five optimal scale levels, and the differences between two adjacent scales are almost uniform.

Algorithm 1 Optimal Scale Iterative Approximation

Require: $I^0=I(x,0)$, T, $\tau$, $\epsilon$, $maxIter$, initial solution $k_1^0,k_2^0,\cdots ,k_N^0$. Ensure: Optimal solution $k_1,k_2,\cdots ,k_N$ 1: for $i=1$ to N do 2:       $k_i=k_i^0$; 3: end for 4: $E_0=\left|\right|\nabla I^0\left|\right|,k_0=0,k_{N+1}=\left[\frac{T}{\tau }\right],\Delta E=inf,j=0$; 5: while $\Delta E<\epsilon$ and $j<maxIter$ do 6:       for $i=1$ to $N+1$ do 7:             $I^{k_i}={\left[1-\tau (k_i-k_{i-1}){\displaystyle \sum _{l=1}^m}{A}_l\left(I^{k_{i-1}}\right)\right]}^{-1}{I}^{k_{i-1}}$; 8:             $E_i=\left|\right|\nabla I^{k_i}\left|\right|$; 9:       end for 10:      $\Delta E=\sqrt{\frac{{\sum }_{i=1}^{N+1}{(E_i-E_{i-1})}^2}{N+1}}$, $j=j+1$; 11:      for $i=1$ to N do 12:            $k_i=k_i-\left[(k_{i+1}-k_{i-1})\frac{E_i-(E_0+i\Delta E)}{E_{i+1}-E_{i-1}}\right]$; 13:      end for 14: end while

4. Experiments and Results

In this study, we arranged several experiments to illustrate the performance and advantages of the proposed method. In real-world remote sensing, the same scene is often observed from different positions/viewpoint withs different surroundings. Subsequently, it may be more difficult to detect the feature points and search the matches. Here, we focused on remote sensing images collected with optical sensors. For a comprehensive comparison, the number of feature points was used to evaluate the cost. The match number and correct number are used to evaluate the efficiency. We also present the time cost of the comparison methods.

In the first experiment, we prepared one remote sensing images and generated four pairs for the test. One of each pair of images was the original image, and the other was generated by rotating and re-scaling the original image. The detected feature points and matching results are shown in Figure 4 and Figure 5. The detailed data (including time cost) can be found in

Table 1. Feature numbers under different rotations.

Method	Keypoint Number 1 (1000)				Keypoint Number 2 (1000)
	0	0	0	0	$\frac{\mathbf{25}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{55}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{85}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{115}\mathbf{\pi }}{\mathbf{180}}$
SIFT	1.146	1.146	1.146	1.146	0.641	0.645	0.624	0.638
SURF	0.170	0.170	0.170	0.170	0.089	0.099	0.091	0.084
BRISK	0.813	0.813	0.813	0.813	0.264	0.271	0.252	0.273
KAZE	0.303	0.303	0.303	0.303	0.168	0.166	0.164	0.172
Proposed	0.402	0.402	0.402	0.402	0.242	0.240	0.238	0.248

and

Table 2. Match numbers under different rotations.

Method	Matches Number				Correct Number				Total Time Cost (s)
	$\frac{\mathbf{25}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{55}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{85}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{115}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{25}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{55}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{85}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{115}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{25}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{55}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{85}\mathbf{\pi }}{\mathbf{180}}$	$\frac{\mathbf{115}\mathbf{\pi }}{\mathbf{180}}$
SIFT	19	21	22	20	19	19	18	16	1.034	1.141	1.114	1.115
SURF	22	16	39	26	20	14	37	19	0.014	0.016	0.016	0.015
BRISK	11	15	18	11	11	15	18	11	0.418	0.441	0.439	0.464
KAZE	30	25	23	33	29	24	23	31	0.063	0.062	0.064	0.067
Proposed	50	44	44	56	48	43	42	53	0.061	0.066	0.064	0.078

.

For the first column of Figure 4, the rotation angle was $\frac{25\pi }{180}$, and the scaling factor was $0.75$. For the next three columns, the scaling factor was fixed, and the rotation angle was set to be $\frac{55\pi }{180}$, $\frac{85\pi }{180}$, and $\frac{115\pi }{180}$, respectively. The rows indicate different methods: SIFT, SURF, BRISK, KAZE, and the proposed method.

It can be found that SIFT finds the most feature points but few matches. Surf finds fewer feature points but obtains more matches than SIFT ol average. BRISK finds the least matches, but there is no error. KAZE can find more correct matches than the three previously mentioned methods but less than the proposed method. In general, because of the optimized scale space, the proposed method achieves the best performance with cost and efficiency being well balanced.

In the second experiment, we tried to test the performances of the different methods under a scaling scene. The rotation angle was fixed to be $\frac{25\pi }{180}$. The original image was scaled using different factor separately to generate four image pairs. The detected feature points and matching results are shown in Figure 6 and Figure 7. The detailed data (including time cost) can be found in

Table 3. Feature numbers under different scales.

Method	Keypoint Number 1 (1000)				Keypoint Number 2 (1000)
	$\mathbf{100}\%$	$\mathbf{100}\%$	$\mathbf{100}\%$	$\mathbf{100}\%$	$\mathbf{55}\%$	$\mathbf{65}\%$	$\mathbf{75}\%$	$\mathbf{85}\%$
SIFT	1.146	1.146	1.146	1.146	0.303	0.421	0.606	0.798
SURF	0.170	0.170	0.170	0.170	0.062	0.078	0.085	0.105
BRISK	0.813	0.813	0.813	0.813	0.111	0.199	0.263	0.410
KAZE	0.303	0.303	0.303	0.303	0.110	0.129	0.173	0.194
Proposed	0.402	0.402	0.402	0.402	0.143	0.188	0.246	0.283

and

Table 4. Match numbers under different scales.

Method	Matches Number				Correct Number				Total Time Cost (s)
	$\mathbf{55}\%$	$\mathbf{65}\%$	$\mathbf{75}\%$	$\mathbf{85}\%$	$\mathbf{55}\%$	$\mathbf{65}\%$	$\mathbf{75}\%$	$\mathbf{85}\%$	$\mathbf{55}\%$	$\mathbf{65}\%$	$\mathbf{75}\%$	$\mathbf{85}\%$
SIFT	8	9	18	32	8	8	17	29	0.755	1.018	1.120	1.224
SURF	14	16	19	35	8	11	18	30	0.013	0.017	0.016	0.016
BRISK	12	7	16	29	12	7	16	29	0.466	0.449	0.455	0.406
KAZE	5	21	30	35	5	18	29	32	0.050	0.064	0.064	0.069
Proposed	12	38	51	56	11	35	49	52	0.056	0.067	0.062	0.076

.

The columns correspond to different scaling factors ($55\%$, $65\%$, $75\%$, $85\%$), and the rows correspond to different methods (SIFT, SURF, BRISK, KAZE, proposed). It can be found that SIFT finds more feature points but fewer matches. Surf finds more matches than SIFT but shows less accuracy, especially when the scaling factor is smaller. BRISK finds the second most feature points but is not advantageous in terms of the number of correct matches. KAZE and the proposed method still find the most correct matches, but KAZE holds no obvious advantage in the match number compared with the other three methods. Our method finds not as many feature points but the most correct matches.

In the third experiment, the rotation angle was fixed to be $\frac{25\pi }{180}$, and the scale factor was fixed to be $75\%$. Gaussian noise was added to the original image to prepare one image of the image pair, and then, the polluted image was rotated and scaled to generate the other image of the pair. The detected feature points and matching results are shown in Figure 8 and Figure 9. The detailed data (including time cost) can be found in

Table 5. Featur numbers under different noise levels.

Method	Keypoint Number 1 (1000)				Keypoint Number 2 (1000)
	-	-	-	-	20 db	18 db	16 db	14 db
SIFT	1.001	0.969	0.945	0.873	0.801	0.810	0.780	0.753
SURF	0.179	0.192	0.213	0.219	0.125	0.118	0.124	0.134
BRISK	0.996	1.219	1.636	2.155	0.440	0.442	0.497	0.570
KAZE	0.303	0.307	0.303	0.301	0.199	0.197	0.195	0.188
Proposed	0.412	0.420	0.414	0.398	0.292	0.285	0.295	0.293

and

Table 6. Match numbers under different noise levels.

Method	Matches Number				Correct Number				Total Time Cost (s)
	20 db	18 db	16 db	14 db	20 db	18 db	16 db	14 db	20 db	18 db	16 db	14 db
SIFT	25	22	23	22	25	21	23	22	0.957	0.967	0.948	0.922
SURF	35	33	39	37	28	30	26	32	0.017	0.015	0.015	0.016
BRISK	35	24	22	23	35	24	22	23	0.392	0.375	0.383	0.395
KAZE	28	29	29	35	26	27	26	34	0.057	0.055	0.056	0.055
Proposed	41	46	52	49	38	44	47	47	0.056	0.058	0.057	0.058

.

The columns correspond to different noise intensities (SNR = 20 db, 18 db, 16 db, 14 db), and the rows correspond to different methods (SIFT, SURF, BRISK, KAZE, proposed).

As shown in Figure 9, SIFT and BRISK find the most key points in each pair, but there is no advantage in finding correct matches. SURF finds the least key points, but there is no obvious disadvantage in finding correct matches. KAZE finds a few more key points than SURF but achieves a similar correct number of matches. The proposed method finds key points a little more than KAZE but achieves the largest correct number of matches.

In the fourth experiment, the rotation angle is fixed to be $\frac{55\pi }{180}$ and the scale factor is fixed to be $65\%$. Four optical remote sensing images are prepared to be scaled and rotated. The feature detection and matching identify are performed on four image pairs. The results are shown in Figure 10 and Figure 11,

Table 7. Features number of different remote sensing images.

Method	Keypoint Number 1 (1000)				Keypoint Number 2 (1000)
	rs1	rs2	rs3	rs4	rs1	rs2	rs3	rs4
SIFT	1.523	1.604	1.896	1.496	0.631	0.698	0.686	0.663
SURF	0.241	0.157	0.188	0.176	0.095	0.049	0.038	0.012
BRISK	0.796	0.697	0.972	0.621	0.252	0.208	0.288	0.133
KAZE	0.163	0.084	0.12	0.053	0.051	0.018	0.015	0.008
Proposed	0.413	0.255	0.349	0.194	0.128	0.07	0.077	0.026

and

Table 8. Matches number of different remote sensing images.

Method	Matches Number				Correct Number				Total Time Cost (s)
	rs1	rs2	rs3	rs4	rs1	rs2	rs3	rs4	rs1	rs2	rs3	rs4
SIFT	20	6	8	2	17	5	8	2	1.035	0.911	0.963	0.915
SURF	24	13	10	6	21	12	10	6	0.014	0.016	0.012	0.012
BRISK	10	3	2	5	10	3	2	5	0.367	0.362	0.368	0.359
KAZE	13	5	4	2	13	5	4	2	0.044	0.043	0.044	0.043
Proposed	27	20	27	8	26	20	27	8	0.051	0.046	0.049	0.045

.

Similarly, the columns are corresponding to different remote sensing image (rs1, rs2, rs3, rs4) and rows are related to different methods (SIFT, SURF, BRISK, KAZE, proposed).

As shown in Figure 11, SIFT finds the most feature points in each pair but no more matches than other methods. BRISK finds the second most feature points and achieve the least correct match number in average. SURF finds more feature points and correct matches than KAZE. The proposed method finds key points a little more than KAZE and achieves the biggest correct number.

From above experiments, the proposed method shows the significantly best comprehensive performance exceed other four.

5. Conclusions

In this paper, a constrained optimization model is introduce to determine some scale slices of a nonlinear scale space for remote sensing image feature detection and matching. With noise estimation introduced to determine the parameter values, a partial differential equation is applied to generate a continuous nonlinear scale space advantageous to represent the important features. The continuous scale space is approximated by a discrete one by Additive Operator Splitting solving the PDE. Then an discrete model is present to fast approach the optimal scale levels for feature detection and image matching. Four experiments are arranged to test the performance of proposed method from views of rotation, scaling, noise and the mixed version. The results show the accuracy and efficiency of the propose method. About 30% improvement in correct matches number with only a small increase in time cost can be found in the proposed method results against the others. In the future work, more efforts will be made to improve the proposed method and develop more efficient and robust method for complex scenes.