Multi-Modal Remote Sensing Image Registration Using Curvature Scale Space Contour Point Features

1. Introduction

Multi-modal remote sensing image registration has become a key area in computer vision and image processing, as single-modal images are insufficient for comprehensive applications [1]. Captured by various sensors with distinct imaging principles, multi-modal images offer unique information that enhances earth observation. The goal is to accurately align images from different times, cameras, or perspectives [2].

Multi-modal remote sensing image registration algorithms are mainly classified into region-based and feature-based methods [3]. Region-based algorithms use image intensity information and optimization to align regions by minimizing a cost function. Feature-based algorithms, which better handle geometric distortions, extract and match features, then select a transformation model based on geometric relationships [4]. Techniques like Scale-Invariant Feature Transform (SIFT) [5] address key multi-modal registration challenges.

Gao et al. [6] proposed the Multi-Scale Partial Intensity Invariant Feature Descriptor (MS-PIIFD) algorithm to address multi-source remote sensing image differences, showing good registration accuracy but being time-consuming and producing fewer matching point pairs, sometimes failing to complete registration. Gao et al. [7] developed the Multi-Scale Histogram of Local Main Orientations (MS-HLMO) algorithm, which handles intensity, scale, and rotation differences, but is also very slow. Li et al. [8] introduced the Radiation-Invariant Feature Transform (RIFT) algorithm to address Non-linear Radiometric Distortion (NRD), offering reliable feature matching but with lower accuracy and longer processing time than MS-PIIFD. Zhu et al. [9] discovered that existing image registering algorithms for infrared and visible images of power equipment suffer from low accuracy and long processing times. To address these issues, they proposed an image registering algorithm based on Large-Gap Fracture Contours (LGFC). This algorithm has shown success in the registering of power equipment images. However, its applicability to multi-modal remote sensing images is limited because remote sensing images generally lack rich LGFC information, which leads to poor registering performance in such images.

Öfverstedt et al. [10] found that intensity-based image registration relies on similarity metrics, which are crucial for robustness and accuracy. They proposed an affine registration framework combining intensity and spatial information using symmetric, non-intensity interpolation. Jiang et al. [11] addressed challenges in multi-modal power equipment images with their Contour Angle Orientation (CAO) method and Coarse-to-Fine (C2F) algorithm (CAO-C2F). Although effective, it still faces issues with accuracy and time consumption.

In response to the above analysis, this paper develops a multi-modal remote sensing image registering algorithm based on Curvature Scale Space (CSS) point features. The main contributions of this paper are as follows.

(1)

The number of feature points is crucial for image registration quality. To ensure an adequate number of feature points, we have enhanced the CSS corner detection algorithm. Since CSS extracts feature points based on contours, and the number of edges correlates with contour quantity, maintaining a sufficient number of edges is vital. To address this, we propose a multi-scale Sobel edge detection algorithm.

(2)

Given the significant differences in intensity, resolution, and viewpoint between multi-modal remote sensing image pairs, we propose a new gradient definition and a method to determine the dominant direction of feature points for rotation invariance. This gradient definition is applied to SIFT descriptors, with segmented normalization to enhance the similarity between feature point descriptors.

The following sections are structured as follows: Section 2 details the proposed registering algorithm. Section 3 conducts an extensive experimental analysis of the algorithm. Lastly, Section 4 provides the conclusion, outlining the key findings and contributions.

2. Proposed Image Registering Algorithm

The study begins by applying the proposed multi-scale Sobel edge detection algorithm to extract image edges, followed by contour tracking to outline these edges. The improved CSS algorithm detects and extracts feature points, with their dominant direction determined for rotation invariance. Enhanced SIFT descriptors are then used to characterize the feature points for image registration, ensuring alignment between the two images.

2.1. Edge Detection

The gradient components in the x and y directions can be calculated using the convolution operation as follows:

$$G_x=S\ast I(x,y),G_y=S^{\top }\ast I(x,y).$$

(1)

$$S=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right),$$

where * denotes the convolution operation, S represents the convolution kernel and $I(x,y)$ denotes the multi-modal remote sensing image.

The gradient magnitude image G can then be obtained using the formula:

$$G=\sqrt{G_x^2+G_y^2}.$$

(2)

We introduce a multi-scale Sobel edge detection approach. Initially, the input color image undergoes a transformation into a grayscale image. Subsequently, we apply a Gaussian filter with a scale parameter of ${\sigma }_0$ to the grayscale image, yielding the first layer of the scale space, denoted as $L_1(x,y,{\sigma }_0)$. Following this, $L_2(x,y,{\sigma }_1)$ is derived by convolving a Gaussian filter with a scale of ${\sigma }_1$ over $L_1(x,y,{\sigma }_0)$, and $L_3(x,y,{\sigma }_2)$ is generated by convolving a Gaussian filter with a scale of ${\sigma }_2$ over $L_2(x,y,{\sigma }_1)$. This iterative process persists until a set of N layers, each of identical dimensions, is established, constituting the scale space. The Gaussian filter function and the relationship between consecutive layers within the scale space are depicted as follows:

$$\left\{\begin{array}{c}{g}_k\left(x,y,{\sigma }_{k-1}\right)={\displaystyle \frac{1}{2\pi {\sigma }_{k-1}^2}}{e}^{-{\displaystyle \frac{x^2+y^2}{2{\sigma }_{k-1}^2}}},\hfill \\ L_k\left(x,y,{\sigma }_{k-1}\right)=g_k\left(x,y,{\sigma }_{k-1}\right)\ast I(x,y),k=1,2,\cdots ,N,\hfill \end{array}\right.$$

(3)

where $L_k(x,y,{\sigma }_{k-1})$ signifies the image obtained through Gaussian filtering.

The selection of an appropriate value for ${\sigma }_{k-1}$ is crucial in this process. Increasing ${\sigma }_{k-1}$ leads to more blurring, resulting in the loss of fine details and nuances. On the other hand, using a smaller ${\sigma }_{k-1}$ value produces a clearer image with more details.

To ensure the required relationship ${\sigma }_{k-1}=f\left(k\right)$ with $f^{\prime }\left(k\right)<0$, we can assume a function $f\left(k\right)={\sigma }_0{k}^u$. Calculating the first derivative of $f\left(k\right)$ yields $f^{\prime }\left(k\right)=u{\sigma }_0{k}^{u-1}<0$, indicating that $u<0$. Consequently, the following relationship between ${\sigma }_{k-1}$ and ${\sigma }_0$ is established:

$${\sigma }_{k-1}={\sigma }_0{k}^u,u<0,k=1,2,3,\cdots ,N.$$

(4)

After constructing the scale space $L_k(x,y,{\sigma }_{k-1})$, $k=1,2,\cdots ,N.$ We apply the single-scale Sobel edge detection method to each scale space image $L_k(x,y,{\sigma }_{k-1})$, resulting in N edge images $E_k(x,y,{\sigma }_{k-1})$. Then, we perform an overlay operation on these edge images $E_k(x,y,{\sigma }_{k-1})$ to obtain an overlaid edge image $E(x,y)$. Last but not least, the overlaid edge image $E(x,y)$ is normalized using Equation (5):

$$\tilde{E}(x,y)={\displaystyle \frac{E(x,y)-a}{b-a}},$$

(5)

where $\tilde{E}(x,y)$ represents the normalized image. The minimum value in matrix $E(x,y)$ is denoted as a, and the maximum value is denoted as b. $E(x,y)-a$ refers to subtracting a from each element in $E(x,y)$, and then dividing each element in the resulting matrix $E(x,y)-a$ by $b-a$.

Considering the binary nature of the edge image, the following technique can be employed:

$$\dot{E}(x,y)=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\tilde{E}(x,y)>T_1,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(6)

where $\dot{E}(x,y)$ represents the pixel value at coordinate $(x,y)$, and $T_1$ is a threshold value. By performing the aforementioned operations, we obtain the resulting edge image $\dot{E}(x,y)$.

2.2. Contour Extraction

Contours offer more continuity than edges, providing geometric details like area and perimeter. Before corner detection, using a contour tracking method is recommended to extract contour sets from multi-modal remote sensing images.

$$L_c=\left\{l^i\mid l^i=\left\{{\left(x_1,y_1\right)}^i,{\left(x_2,y_2\right)}^i,\cdots ,{\left(x_n,y_n\right)}^i\right\}\right\},$$

(7)

where $i=1,2,\cdots ,M,$$L_c$ denotes the set of extracted contours, l represents the elements within each contour, i signifies the i-th contour, and M is the total number of contours in $L_c$.

2.3. Feature Point Detection

The classical corner detection algorithm, CSS [12], computes the absolute curvature of identified corners along the contour at a low scale. It uses the local maximum of the absolute curvature as initial candidate corner points and employs an adaptive algorithm to eliminate rounded corners, thus reducing false corner detections.

Each contour extracted from a multi-modal remote sensing image can be treated as a curve, represented by an arc length parameter v:

$$\beta \left(v\right)=\left[x\right(v),y(v\left)\right],$$

(8)

where $x\left(v\right)$ and $y\left(v\right)$ denote the sequence of horizontal and vertical coordinates constituting each curve.

The curve $\beta \left(v\right)$ is smoothed using a Gaussian function with a scale parameter $\alpha$ to derive the smooth filter curve ${\beta }_{\alpha }\left(v\right)$.

$$\left\{\begin{array}{c}{\beta }_{\alpha }\left(v\right)=[x\left(v\right)\ast g(v,\alpha ),y\left(v\right)\ast g(v,\alpha )]=[X(v,\alpha ),Y(v,\alpha )],\hfill \\ g(v,\alpha )={\displaystyle \frac{1}{\sqrt{2\pi }\alpha }}{e}^{-{\displaystyle \frac{v^2}{2{\alpha }^2}}},\hfill \end{array}\right.$$

(9)

where the convolution operator symbol * is employed in the context, and $g(v,\alpha )$ represents a one-dimensional Gaussian function.

The absolute curvature of each point on the curve can be determined by:

$$K(u,\alpha )={\displaystyle \frac{\left|X_v(v,\alpha )Y_{vv}(v,\alpha )-X_{vv}(v,\alpha )Y_v(v,\alpha )\right|}{{\left[X_v{(v,\alpha )}^2+Y_v{(v,\alpha )}^2\right]}^{\frac{3}{2}}}}.$$

(10)

where $X_v(v,\alpha )=x\left(v\right)\ast g_v(v,\alpha )$, $Y_v(v,\alpha )=y\left(v\right)\ast g_v(v,\alpha )$, $X_{vv}(v,\alpha )=x\left(v\right)\ast g_{vv}(v,\alpha )$, and $Y_{vv}(v,\alpha )=y\left(v\right)\ast g_{vv}(v,\alpha )$.

By utilizing Equations (8)-(10), we can calculate the absolute curvature for each point along the contour. Subsequently, the CSS [12] selects initial corner candidates at a lower scale based on the maximum value of the absolute curvature $K(v,\alpha )$.

Since the CSS is based on contours, it cannot detect a sufficient number of feature points if the image lacks rich edge detail information. Therefore, we have improved CSS to increase the detection of corner information. The possible relationships between the line segments connecting adjacent feature points on each contour $l^i(i=1,2,\cdots ,M)$ are illustrated in Figure 1.

As shown in Figure 1(a), for the case where the line segment connecting feature points $P_L^i(x_L^i,y_L^i)$ and $P_R^i(x_R^i,y_R^i)$ is a straight line, we use the midpoint $P_M^i(x_M^i,y_M^i)$ of the line segment $P_L^i{P}_R^i$ as a feature point. The midpoint $P_M^i$ can be calculated using the following formula:

$$\begin{array}{cc}\hfill P_M^i& =(x_M^i,y_M^i)={\displaystyle \frac{1}{2}}(P_L^i+P_R^i)={\displaystyle \frac{1}{2}}(x_L^i+x_R^i,y_L^i+y_R^i).\hfill \end{array}$$

(11)

As shown in Figure 1(b) and Figure 1(c), the connecting curve between feature points $P_L^i$ and $P_R^i$ is a curved segment, exhibiting concave and convex shapes respectively. This connecting line segment is continuous between feature points $P_L^i$ and $P_R^i$. According to Lagrange’s theorem, as long as the function is continuous on the closed interval and differentiable on the open interval, there must be at least one point $P_M^i$ where Equation (12) holds.

$$F\left(x_M^i\right)={\displaystyle \frac{y_L^i-y_R^i}{x_L^i-x_R^i}},$$

(12)

where $F(·)$ is the derivative of the function $f(·)$ that represents the curved segment between feature points $P_L^i$ and $P_R^i$.

From Figure 1(b) and Figure 1(c), it can be observed that the straight line segment $P_L^i{P}_R^i$ connecting $P_L^i$ and $P_R^i$ is parallel to the tangent $L_t^i$ at point $P_M^i$, and there is only one such tangent. Therefore, the tangent point $P_M^i$ is used as a feature point, and $P_M^i$ is calculated as shown in Equation (13).

$$\begin{array}{cc}\hfill P_M^i& =(x_M^i,y_M^i)=(F^{-1}\left(x_m^i\right),f\left(F^{-1}\left(x_m^i\right)\right)),\hfill \end{array}$$

(13)

where $F^{-1}(·)$ is the inverse function of $F(·)$.

After the above processing, the number of feature points extracted by the refined CSS corner detection algorithm is ${\sum }_{i=1}^M\left(2N^i-1\right)$. Here, $N^i$ represents the number of feature points detected by the original CSS algorithm in the i-th contour, and M denotes the total number of contours extracted from the image.

2.4. Dominant Direction

Assigning a dominant direction linked to the local gradient for each feature point is essential for achieving rotation, translation, and scaling invariance, which helps in generating effective descriptors. We propose a novel method to determine the dominant direction, specifically optimized to handle significant intensity variations in multi-modal remote sensing image pairs. For the image $I(x,y)$, the gradient magnitude image G is computed using Equation (2). Subsequently, Equation (1) is utilized to calculate the first horizontal and vertical steps, $G_{1x}$ and $G_{1y}$, of the gradient magnitude image G, and the gradient magnitude $G_1$ is determined as:

$$G_1=\sqrt{G_{1x}^2+G_{1y}^2},$$

(14)

where $G_1$ denotes the gradient magnitude image of G.

Next, proceed to normalize the gradient magnitude of the gradient magnitude image $G_1$, as defined by the operation:

$$\overline{G}=U\left(G_1\right),$$

(15)

where the symbol U denotes the normalization operation, which involves dividing all values in $G_1$ by its maximum value. The resulting $\overline{G}$ represents the normalized gradient magnitude image.

Subsequently, calculate the horizontal and vertical partial derivatives based on the normalized gradient image using Equation (1). This can be expressed as:

$$\begin{array}{c}\left(\begin{array}{c}{\overline{G}}_x\hfill \\ {\overline{G}}_y\hfill \end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\partial }{\partial x}}\left(\overline{G}\right)\\ {\displaystyle \frac{\partial }{\partial y}}\left(\overline{G}\right)\end{array}\right),\hfill \end{array}$$

(16)

where ∂ denotes the operation of partial derivative.

Equation (16) represents the new gradient definition proposed in this paper. The weighted squared gradient in the average squared gradient method is defined as follows:

$$\left[\begin{array}{c}{G}_{w_r,s,x}\hfill \\ G_{w_r,s,y}\hfill \end{array}\right]=\left[\begin{array}{c}\sum _{w_r}{\overline{G}}_x^2-{\overline{G}}_y^2\\ 2\sum _{w_r}{\overline{G}}_x{\overline{G}}_y\end{array}\right],$$

(17)

where $w_r$ represents a Gaussian window with variance r. Subsequently, the dominant direction of the feature points can be determined by:

$$D=\left|\angle \left(G_{w_r,s,x},G_{w_r,s,y}\right)\right|.$$

(18)

In Equation (18), the absolute value operation $|·|$ is applied to confine the dominant direction range of feature points from $(-\pi ,\pi )$ to $(0,\pi )$. The variable D denotes the dominant direction of the specific feature point. The calculation of $\angle \left(G_{w_r,s,x},G_{w_r,s,y}\right)$ is determined using:

$$\angle (X,Y)=\left\{\begin{array}{c}arctan\left({\displaystyle \frac{Y}{X}}\right),X<0,Y<0,\hfill \\ arctan\left({\displaystyle \frac{Y}{X}}\right)+\pi ,X⩾0,\hfill \\ arctan\left({\displaystyle \frac{Y}{X}}\right)+2\pi ,X<0,Y\ge 0.\hfill \end{array}\right.$$

(19)

2.4.1. Feature Descriptor Construction

Due to variations in resolution, intensity, and other attributes, the SIFT descriptor struggles with multi-modal images. The traditional SIFT, based on the first-order gradient, fails to effectively address these disparities. To improve this, we propose a novel gradient definition within the SIFT framework, using Equation (16) to calculate a more robust gradient magnitude $G_2$, enhancing the similarity of feature descriptors across multi-modal remote sensing images.

$$G_2=\sqrt{{\overline{G}}_x^2+{\overline{G}}_y^2}.$$

(20)

Utilizing a methodology akin to that outlined in Equation (15), the process of normalization is executed as represented by:

$${\overline{G}}_2=U\left(G_2\right).$$

(21)

In addition, we applied piecewise normalization to the gradient magnitude of the normalized gradient image ${\overline{G}}_2$. Specifically, the gradient magnitudes were sorted in descending order and then normalized accordingly.

2.4.2. Coarse-to-Fine Feature Matching

We use the Best Bin First (BBF) [13] matching algorithm with generated feature descriptors, performing two-sided matching to establish initial correspondences between feature points.

In the context of bilateral matching of two multi-modal remote sensing images, denoted as $R(x,y)$ (reference image) and $F(x,y)$ (floating image), we identify the nearest p and second-nearest q neighbors related to a point o in $R(x,y)$, as well as the nearest f and second-nearest v neighbors corresponding to a point e in $F(x,y)$. The distance between two points is computed:

$$d_{ij}=\sqrt{{\sum }_{k=1}^{128}{\left[R_i^k(x,y)-F_j^k(x,y)\right]}^2}.$$

(22)

We enforce a distance ratio threshold $T_2$, as defined in Equation (23), where satisfying the conditions ensures successful matching between points from both images.

$$\left\{\begin{array}{c}{\displaystyle \frac{d_{op}}{d_{oq}}}\le T_2,\hfill \\ {\displaystyle \frac{d_{ef}}{d_{ev}}}\le T_2.\hfill \end{array}\right.$$

(23)

We use the line consistency theory [11] to filter out most of the incorrect matches. Finally, to achieve more robust feature point pairs, we apply the Random Sample Consensus (RANSAC) algorithm [14] for fine-tuning the matches.