An optical and sar heterogeneous satellite image matching method fusing spatial and phase features

Abstract

This invention discloses a method for matching optical and SAR heterogeneous satellite images by fusing spatial and phase features. Using a reference image and the image to be matched as the original images, an edge map is obtained using the maximum response map of phase consistency as a phase feature. Gaussian filtering is applied to the original images as spatial features. The two images are combined as a guide map to construct a corresponding joint bilateral filtering scale pyramid. For each scale image, a maximum moment map, a minimum moment map, and a maximum index map are calculated. Feature points are detected on the maximum and minimum moment maps, and feature descriptors are calculated on the maximum index map. Two principal directions are set for each feature descriptor, constructing a dual feature description vector. All feature points are mapped onto the corresponding original images, and their corresponding dual feature description vectors are used as input. A matching method is employed to obtain an initial matching result. The initial matching result is then used as prior information for a secondary matching to obtain the final matching result.

Description

Technical Field

[0001] This invention relates to the field of laser imaging technology, and more particularly to a method for matching optical and SAR heterogeneous satellite images by fusing spatial and phase features. Background Technology

[0002] Image matching refers to the process of detecting corresponding points of the same scene in different images and matching these points one-to-one through reliable correspondences. With the development of multi-source sensors and photogrammetry technology, optical and SAR heterogeneous satellite image matching technology has become a hot topic. Unlike single-sensor image matching, optical and SAR images have significant nonlinear radiometric differences due to their different imaging modes.

[0003] Currently, image matching research methods can be divided into three categories: region-based methods, feature-based methods, and deep learning-based methods. Region-based methods mainly establish similarity metrics and construct matching templates using regions around predefined feature points or the entire image. Image matching is then performed by calculating the similarity of the matching templates. Feature-based methods extract prominent feature points, such as corner points and edge points, from the image to be matched. After obtaining the feature points, feature descriptions are performed, and matching is performed using the constructed feature descriptors. Deep learning-based methods have been increasingly used in recent years. They typically use convolutional neural networks to calculate features and then match the features.

[0004] Most matching algorithms experience a significant decrease in the number of correctly matched points when there is a rotation difference between two images compared to when there is no rotation difference. In 2020, Jiayuan Li proposed the RIFT (Radiation-variation Insensitive Feature Transform) algorithm, which achieves rotation invariance by constructing a convolutional sequence loop. However, if the rotation angle difference between the two images is about half of the two filtering angles of the sequence loop, the matching effect will be greatly reduced. Furthermore, due to the lack of scale invariance, it cannot match images of different scales. Therefore, designing a matching method that can stably match heterogeneous remote sensing images with rotation and scale differences is crucial. Summary of the Invention

[0005] This invention provides a method for matching optical and SAR heterogeneous satellite images that integrates spatial and phase features, solving the technical problems in the prior art where matching is impossible on images with significant scale differences due to the lack of scale invariance, and where rotation invariance is unstable.

[0006] This invention can be achieved through the following technical solutions:

[0007] A method for matching optical and SAR heterogeneous satellite images that integrates spatial and phase features includes the following steps:

[0008] Step 1: Construct scale pyramids corresponding to the reference image and the image to be matched.

[0009] Using the reference image and the image to be matched as the original images, a phase consistency map is first constructed using filters, and the corresponding maximum response map is obtained. Then, edge detection is performed to obtain the edge map as the phase feature map. At the same time, Gaussian filtering is applied to the original image to obtain the spatial feature map. The spatial feature map and the phase feature map are superimposed to obtain the guide map. Finally, a scale pyramid is constructed using joint bilateral filtering.

[0010] Step 2: Detect feature points on each scale of the image in the size pyramid and calculate the corresponding feature descriptors. Each feature descriptor has two principal directions: the gradient principal direction and the centroid principal direction. Then, rotate the X-axis direction of the feature descriptor to the two principal directions and construct the corresponding descriptor vector.

[0011] Step 3: Map all feature points onto the corresponding original image. Using the two descriptor vectors corresponding to the feature points as input, perform two matching operations on the reference image and the image to be matched to obtain the final matching result.

[0012] Furthermore, in step two, the phase consistency map corresponding to each scale image is first calculated, and then the corresponding maximum moment image and minimum moment image are calculated according to the moment analysis method. Then, the feature detector is used to detect feature points on the maximum moment image and minimum moment image. Finally, the feature descriptor is calculated using the neighboring region of the feature points as input and the feature descriptor calculation method.

[0013] Furthermore, when calculating the principal direction of the gradient,

[0014] (1) Map the feature points on the minimum moment image to the maximum moment image, and denote the corresponding points as the external feature points;

[0015] (2) Calculate the gradient magnitude and gradient direction of all pixels in the maximum moment image;

[0016] (3) Divide the gradient direction from 0 to 360 degrees into 36 columns, each column being 10 degrees, and calculate the gradient histogram of the feature point and the neighboring region of the added feature point.

[0017] (4) Take the peak direction of the gradient histogram as the main gradient direction of the feature descriptor, and rotate the x-axis of the feature descriptor to the main gradient direction.

[0018] Furthermore, when calculating the principal directions of the centroid,

[0019] (1) Find the pixel that appears most often in the neighborhood of each feature point and record it as the repeated pixel.

[0020] (2) Calculate the centroid of the region enclosed by all the repeating pixels using the gray-scale centroid method;

[0021] (3) The direction connecting the feature point to the corresponding centroid point is taken as the main centroid direction, and the x-axis of the feature descriptor is rotated to the main centroid direction.

[0022] Furthermore, when calculating the feature descriptor, based on obtaining the phase consistency map corresponding to each scale image, the corresponding maximum response map is calculated, and the corresponding maximum index map is constructed. Then, the feature descriptor is calculated using a method similar to SIFT feature descriptor calculation, with the neighboring region of the corresponding point of the feature point on the maximum index map as input.

[0023] Furthermore, in step three, during the initial matching process, the nearest neighbor matching method is first used for initial matching, then the FSC algorithm is used to remove erroneous matching points, and then the homography matrix is ​​used to perform homography transformation on the image to be matched, obtaining a new image to be matched that has no rotation or size change compared to the reference image.

[0024] Finally, the final matching is performed. The feature points and rotation-independent feature descriptors of the new image to be matched are recalculated, as well as the rotation-independent feature descriptors corresponding to the feature points retained on the reference image. Then, the pixel point corresponding to the feature point of the reference image is found in the new image to be matched. Matching search is performed in the circular region centered on the pixel point. The nearest neighbor matching method and the FSC algorithm are used to eliminate mismatches and obtain the correct matching point.

[0025] Furthermore, in step one, a two-dimensional Log-Gabor filter is first used to generate phase consistency maps corresponding to the reference image and the image to be matched. At the same time, the two-dimensional Log-Gabor filter is used to calculate the maximum response map corresponding to the phase consistency map. Then, the Canny edge detection algorithm is used to obtain the edge map corresponding to the maximum response map as the phase feature map. The original image is Gaussian filtered as the spatial feature map. The obtained phase feature map is combined with the spatial feature map as the guide map for joint bilateral filtering.

[0026] The beneficial technical effects of this invention are as follows:

[0027] 1. By utilizing the characteristic that the Log-Gabor filter has a large response value at the edge of the phase consistency map, the Canny edge detection algorithm is used to extract image edges. Combined with Gaussian filtering, a guide map of joint bilateral filtering is formed, thereby constructing a joint bilateral filtering scale pyramid. Feature points and feature descriptors are extracted at each layer of the scale pyramid, thereby achieving scale invariance.

[0028] 2. The maximum response map obtained using phase consistency has particularly large response values ​​at the edges. In this case, Canny edge detection can obtain a relatively stable edge image, i.e., a phase feature map. Superimposing the original image after Gaussian filtering creates a spatial feature map, which can sharpen details. This is used as a guide map for joint bilateral filtering. Compared with directly using Gaussian filtering, it will preserve edges. Compared with directly using the original image as a guide map, joint bilateral filtering greatly increases the image denoising ability, thereby increasing the robustness of feature point selection. At the same time, feature description is not easily interfered with, providing a good data foundation for subsequent matching.

[0029] 3. The main direction of the feature descriptor is constructed by calculating the main direction of the gradient of the maximum moment image and the main direction of the centroid of the neighboring region of the feature point. This enables the construction of a dual feature descriptor vector for the same feature point, thus ensuring robustness in matching rotated images.

[0030] 4. After matching is completed, a homography matrix between the two images is constructed using the obtained correct matching points. The homography matrix is ​​used to perform homography transformation on the image to be matched. Then, an algorithm that does not consider rotation is used to match the reference image and the new homography transformed image to be matched to obtain new correct matching points. Finally, these matching points are transformed to the coordinate system of the original image to be matched through inverse homography transformation, thereby greatly increasing the number of matching points. Attached Figure Description

[0031] Figure 1 This is a schematic diagram of the overall process of the present invention;

[0032] Figure 2 This is a schematic diagram illustrating the construction guide of the present invention;

[0033] Figure 3 This is a schematic diagram of the scale pyramid generated using joint bilateral filtering according to the present invention;

[0034] Figure 4 This is a schematic diagram of the maximum index map of the present invention, wherein (a) represents the original image; and (b) represents the maximum index map.

[0035] Figure 5 This invention describes the process of constructing a rotation-invariant dual description vector;

[0036] Figure 6 This is a schematic diagram illustrating the process of obtaining a new image to be matched without rotation or size influence after the initial matching in this invention.

[0037] Figure 7 This is a schematic diagram of the matching search range used in the final matching process of this invention;

[0038] Figure 8 These are some of the experimental data from the embodiments of the present invention;

[0039] Figure 9 This is a schematic diagram of the matching results using the existing RIFT algorithm;

[0040] Figure 10 This is a schematic diagram showing the matching result using the matching method of the present invention;

[0041] Figure 11 This is a schematic diagram showing the results of matching images rotated from 0° to 330° in an embodiment of the present invention.

[0042] Figure 12 This is a schematic diagram illustrating the relationship between the number of correctly matched points and the relative rotation angle of the image in an embodiment of the present invention. Detailed Implementation

[0043] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings and preferred embodiments.

[0044] like Figure 1 As shown, this invention provides a method for matching optical and SAR heterogeneous satellite images by fusing spatial and phase features. First, a scale pyramid corresponding to a reference image and an image to be matched is constructed. Then, feature points on each scale image in the scale pyramid are detected and corresponding feature descriptors are calculated. Each feature descriptor has two principal directions: the gradient principal direction and the centroid principal direction. The feature descriptors are then rotated to the two principal directions and corresponding descriptor vectors are constructed. Finally, all feature points are mapped onto the corresponding original images. Using the two descriptor vectors corresponding to the feature points as input, a matching method is used to perform two matching operations on the reference image and the image to be matched to obtain the final matching result.

[0045] Specifically as follows:

[0046] S1. Constructing the Scale Pyramid

[0047] Using the reference image and the image to be matched as the original images, the phase consistency map of each image is calculated using a Log-Gabor filter. At the same time, the maximum response map of the image is obtained. Canny edge detection is applied to the maximum response map, and a guide map is generated by combining Gaussian filtering. Finally, a scale pyramid is constructed using joint bilateral filtering.

[0048] Phase consistency models have been shown to be robust for multimodal image matching. This invention uses a two-dimensional Log-Gabor filter to generate phase consistency maps.

[0049] The two-dimensional Log-Gabor function is defined as follows:

[0050]

[0051] In the formula, (ρ,θ) represents logarithmic polar coordinates; s and o are the scale and direction of the two-dimensional Log-Gabor coordinates, respectively; (ρ s ,θ so ) is the center frequency of the two-dimensional Log-Gabor; s ρ and s θ These are the bandwidths in ρ and θ, respectively.

[0052] A two-dimensional Log-Gabor filter is a frequency domain filter whose spatial domain definition is:

[0053] L(x,y,s,o)=L even (x,y,s,o)+iL odd (x,y,s,o)

[0054] Where (x,y) represents the pixel coordinates, L even (x,y,s,o) denotes the real part, and L odd (x,y,s,o) represents the imaginary part, which represents the even-symmetric Log-Gabor wavelet and the odd-symmetric Log-Gabor wavelet, respectively.

[0055] The response component E can be obtained by convolving the image with even-symmetric and odd-symmetric Log-Gabor wavelets respectively. so (x,y),O so (x,y) is defined as follows:

[0056] [E so (x,y),O so (x,y)]=[I(x,y)*L even (x,y,s,o),I(x,y)*L odd (x,y,s,o)]

[0057] The amplitude component A of the pixel at different scales and in different directions can be calculated using these two response components. so (x,y) and phase component φ so (x,y):

[0058]

[0059] φ so (x,y)=arctan(E so (x,y) / O so (x,y))

[0060] By combining noise compensation terms across all directions and scales, the final phase coherence expression is obtained:

[0061]

[0062] Where, ω o (x, y) is a weighted function, where x is a local minimum, and ΔΦ so (x,y) is the phase deviation function, defined as:

[0063]

[0064]

[0065]

[0066]

[0067] After obtaining the phase consistency map, the maximum response map of the image needs to be calculated by utilizing the characteristic of the Log-Gabor filter that has a large response value at the image edges, so as to obtain a more stable image with prominent edges.

[0068] The maximum response plot is calculated as follows:

[0069] MRP(x,y)=max(A o (x,y))

[0070]

[0071] Where MRP(x,y) represents the maximum sum of amplitude components at all scales in different directions of the maximum response map, and max(·) is used to locate the maximum value in different directions in the image sequence. o (x,y) represents the sum of amplitude components of all scales in the same direction, and N represents the number of scales.

[0072] Then, the Canny edge detection algorithm is used to obtain the edge map, i.e., the phase feature map, of the maximum response map. After obtaining the phase feature map, Gaussian filtering is applied to the original image to obtain the spatial feature map. The spatial feature map and the phase feature map are then directly superimposed to obtain the final guide map required for joint bilateral filtering. Figure 2 As shown.

[0073] Guide(x,y)=Canny(MRP(x,y))+G I (x,y)

[0074] Where Canny(·) represents performing Canny edge detection on the image, G I (x,y) represents the original image after Gaussian filtering.

[0075] After obtaining the guide map, a scale pyramid can be constructed using joint bilateral filtering. The joint bilateral filtering formula is as follows:

[0076]

[0077] Among them, J p I q These are the output and input pixel values ​​of the image, respectively, where p and q are the pixel indices f(p,q) and g(Guideline) in the image, respectively. p Guide q The weights for the output pixels are represented by ), all of which are Gaussian functions. A joint bilateral filtering scale pyramid is constructed layer by layer with a ratio of 1.4. Figure 3 As shown.

[0078] S2, Feature Point Detection

[0079] Calculate the phase consistency map of each scale image in the scale pyramid of the reference image and the image to be matched. Generate the maximum moment image and the minimum moment image using the moment analysis method. Use the Fast feature detector to detect feature points on the maximum moment image and the minimum moment image respectively.

[0080] According to the moment analysis method, the axis corresponding to the minimum moment is called the principal axis. The principal axis usually represents the direction information of the feature point. The axis corresponding to the maximum moment is perpendicular to the principal axis. The magnitude of the maximum moment usually reflects the salience of the feature point. Therefore, the maximum moment and the minimum moment can better describe the edge points and corner features of the image.

[0081] Maximum moment M max and minimum moment M mim The calculation method is as follows:

[0082]

[0083]

[0084]

[0085]

[0086]

[0087] After obtaining the moment map, the Fast feature detector is used to detect feature points on the maximum moment image and the minimum moment image, which yields a large number of edge points and corner points.

[0088] S3, Calculate feature descriptors

[0089] After obtaining the feature points, a feature descriptor needs to be created for each feature point. Based on obtaining the phase consistency map corresponding to each scale image, the corresponding maximum response map is calculated using the methods described above. Since max(A) has already been obtained when calculating the maximum response map... o Given (x, y), we only need to find the index of the corresponding maximum value, i.e., which 'o' makes A = 0. o The maximum index (x, y) is found at the coordinate (x, y), and this value is used as the pixel value at that coordinate. The maximum index (MaxIndex(x, y)) can be obtained by calculating the index of the corresponding pixel across the entire image. Figure 4 As shown, this graph can be directly used for feature description calculations.

[0090] MaxIndex(x,y)=index(max(A o (x,y)))

[0091] The index(·) function is used to obtain the index of the maximum value in the image sequence.

[0092] Feature descriptors are calculated based on the maximum index map and detected feature points. Two principal directions for each feature point are obtained by calculating the principal gradient direction of the moment map at the feature point and the principal centroid direction of the feature descriptor. The feature descriptors are then rotated to these two directions to obtain two rotated feature descriptors, forming a dual feature description vector. Figure 5 As shown, the details are as follows:

[0093] After obtaining the maximum index map, a feature description method similar to SIFT is used to divide the neighboring region centered on the feature point, such as a 96×96-sized image block, into 8×8 sub-grids at equal intervals. A distribution histogram without partitions is built for each sub-grid, and the number of different index values ​​is counted. Since the value range of the maximum index map is between 0 and N-1, the dimension of the feature descriptor is 8×8×N.

[0094] To avoid the effects of lighting changes, the feature descriptors are finally normalized. However, the feature descriptors at this point do not possess rotation invariance because the maximum index map cannot, in essence, be like the gradient. Figure 1 The principal direction of a feature descriptor is determined by using its gradient histogram. Therefore, the principal direction of the feature descriptor is determined by calculating the gradient histogram of the maximum moment image at the corresponding feature point. The main steps are as follows:

[0095] (1) Map the feature points on the minimum moment image to the maximum moment image, and denote the corresponding points as the external feature points;

[0096] (2) Calculate the gradient magnitude and gradient direction of all pixels in the maximum moment image;

[0097] (3) Divide the gradient direction from 0 to 360 degrees into 36 columns, each column being 10 degrees, and calculate the gradient histogram of the feature point and the neighboring region of the added feature point.

[0098] (4) Take the peak direction of the gradient histogram as the main gradient direction of the feature descriptor, and rotate the x-axis of the feature descriptor to the main gradient direction.

[0099] Meanwhile, since the maximum moment image is relatively sparse, in order to enhance the effect of rotation invariance, the centroid direction of the dominant region of the feature descriptor is increased. The main steps are as follows:

[0100] (1) Find the pixel that appears most often in the neighborhood of each feature point and record it as the repeated pixel.

[0101] (2) Calculate the centroid (x) of the region enclosed by all repeating pixels using the gray-scale centroid method. c ,y c The calculation formula is as follows;

[0102]

[0103]

[0104] (3) The direction connecting the feature point to the corresponding centroid point is taken as the main centroid direction, and the x-axis of the feature descriptor is rotated to the main centroid direction.

[0105] S4, Two-way matching

[0106] S4.1 Initial Matching

[0107] The nearest neighbor matching method is used to perform initial matching on the obtained dual feature description vectors, and the FSC algorithm is used to quickly remove mismatched points. The remaining matching points are recorded as initial matching points. The homography matrix of the image to be matched to the reference image is calculated using the initial matching points. The homography matrix is ​​then used to perform homography transformation on the image to be matched, obtaining a new image to be matched that has no rotation or size change compared to the reference image. Figure 6 As shown.

[0108] S4.2 Final Match

[0109] After obtaining a new image to be matched, due to the rotation angle between the two images and the change in image size, which has been largely eliminated through homography matrix transformation, it is necessary to re-detect feature points using a Fast feature detector and construct a feature descriptor that does not consider rotation (i.e., calculation without principal direction) using a method similar to SIFT feature description. For the reference image, the feature descriptor corresponding to the retained feature points that does not consider rotation needs to be recalculated, and then the nearest neighbor matching method is used again for final matching. However, after the initial matching, ideally, pixels with the same coordinates in the two images should represent the same ground position. However, due to the low accuracy of the initial matching and the small number of matching points, there will be errors in the pixel correspondence of the images. This error will be within a certain range. Therefore, in the final matching, it is only necessary to limit the matching search range of the reference image in the new image to be matched. That is, if the feature point coordinates of the reference image are (x, y), when performing the matching search, it is only necessary to find the pixel with the same coordinates in the new matching image and search within a circular area centered on that pixel, instead of searching and matching all feature points in the entire image as in the initial matching. Figure 7 As shown, limiting the search range can further improve the robustness of the matching results and greatly increase the number of matched feature points.

[0110] Next, the FSC algorithm is used again to remove incorrect matching points. Then, a uniform homography transformation matrix is ​​calculated using the remaining correct matching points. All matching points that conform to the homography transformation matrix are added to the set of correct matching points. This operation is repeated until no more matching points are added. Finally, the iterated matching points are projected back onto the original image to be matched to obtain the final matching result.

[0111] To verify the feasibility of this invention, we conducted the following experiment:

[0112] All experiments were conducted in a Windows 11 environment with the following CPU: AMD Ryzen 9 5900X 12-Core Processor, GPU: NVIDIA GeForce RTX 3080, RAM: 64GB, programming language: C++11, and OpenCV version: 3.4.5.

[0113] In matching multiple pairs of optical and SAR images, the matching method of this invention exhibits good performance. Figure 8 The experiment data is shown in part. Each pair of data includes one optical image and one SAR image. The experimental results using the existing RIFT algorithm and the matching method of this invention are as follows: Figure 9 , 10 As shown.

[0114] from Figure 9 It is evident that the RIFT algorithm, lacking scale invariance and exhibiting significant scale differences in pairs 2, 3, and 5, loses its matching ability on these images. Furthermore, on pair 4, the RIFT algorithm achieves rotation invariance by constructing a loop in the convolution sequence, but this method's resistance to rotation is unstable, thus failing to match successfully. In contrast, the algorithm of this invention demonstrates excellent matching performance on all images. Figure 10 As shown, even with significant scale variations in the image, a large number of correct matching points can be obtained, and the overall performance is significantly better than the RIFT algorithm. The quantitative matching metrics are shown in the table below:

[0115] Quantitative analysis and comparison of experimental results

[0116]

[0117] To test the matching algorithm's performance when images are rotated at different angles, images were rotated from 0° to 330° for matching. Figure 11 It can be observed that the number of matching points changes relatively steadily. Figure 12 .

[0118] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these are merely illustrative examples. Various changes or modifications can be made to these embodiments without departing from the principles and essence of the present invention. Therefore, the scope of protection of the present invention is defined by the appended claims.

Claims

1. A method for matching optical and SAR heterogeneous satellite images by fusing spatial and phase features, characterized in that... Includes the following steps: Step 1: Construct scale pyramids corresponding to the reference image and the image to be matched. Using the reference image and the image to be matched as the original images, a phase consistency map is first constructed using filters, and the corresponding maximum response map is obtained. Then, edge detection is performed to obtain the edge map as the phase feature map. At the same time, Gaussian filtering is applied to the original image to obtain the spatial feature map. The spatial feature map and the phase feature map are superimposed to obtain the guide map. Finally, a scale pyramid is constructed using joint bilateral filtering. Step 2: Detect feature points on each scale of the image in the size pyramid and calculate the corresponding feature descriptors. Each feature descriptor has two principal directions: the gradient principal direction and the centroid principal direction. Then, rotate the X-axis direction of the feature descriptor to the two principal directions and construct the corresponding descriptor vector. Step 3: Map all feature points onto the corresponding original image. Using the two descriptor vectors corresponding to the feature points as input, perform two matching operations on the reference image and the image to be matched to obtain the final matching result.

2. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 1, characterized in that: In step two, the phase consistency map corresponding to each scale image is first calculated, and then the corresponding maximum moment image and minimum moment image are calculated according to the moment analysis method. Then, the feature detector is used to detect feature points on the maximum moment image and minimum moment image. Finally, the feature descriptor is calculated using the neighboring region of the feature points as input.

3. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 2, characterized in that: When calculating the principal direction of the gradient (1) Map the feature points on the minimum moment image to the maximum moment image, and denote the corresponding points as the external feature points; (2) Calculate the gradient magnitude and gradient direction of all pixels in the maximum moment image; (3) Divide the gradient direction from 0 to 360 degrees into 36 columns, each column being 10 degrees, and calculate the gradient histogram of the feature point and the neighboring region of the added feature point. (4) Take the peak direction of the gradient histogram as the main gradient direction of the feature descriptor, and rotate the x-axis of the feature descriptor to the main gradient direction.

4. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 2, characterized in that: When calculating the principal directions of the centroid (1) Find the pixel that appears most often in the neighborhood of each feature point and record it as the repeated pixel. (2) Calculate the centroid of the region enclosed by all the repeating pixels using the gray-scale centroid method; (3) The direction connecting the feature point to the corresponding centroid point is taken as the main centroid direction, and the x-axis of the feature descriptor is rotated to the main centroid direction.

5. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 2, characterized in that: When calculating feature descriptors, based on obtaining the phase consistency map corresponding to each scale image, the corresponding maximum response map is calculated, and the corresponding maximum index map is constructed. Then, the feature descriptor is calculated using a method similar to SIFT feature descriptor calculation, with the neighboring region of the corresponding point on the maximum index map as input.

6. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 1, characterized in that: In step three, during the initial matching process, the nearest neighbor matching method is first used for initial matching, then the FSC algorithm is used to remove erroneous matching points, and finally, the homography matrix is ​​used to perform homography transformation on the image to be matched, obtaining a new image to be matched that has no rotation or size change compared to the reference image. Finally, the final matching is performed. The feature points and rotation-independent feature descriptors of the new image to be matched are recalculated, as well as the rotation-independent feature descriptors corresponding to the feature points retained on the reference image. Then, the pixel point corresponding to the feature point of the reference image is found in the new image to be matched. Matching search is performed in the circular region centered on the pixel point. The nearest neighbor matching method and the FSC algorithm are used to eliminate mismatches and obtain the correct matching point.

7. The optical and SAR heterogeneous satellite image matching method fusing spatial and phase features according to claim 1, characterized in that: In step one, a two-dimensional Log-Gabor filter is first used to generate phase consistency maps corresponding to the reference image and the image to be matched. At the same time, the two-dimensional Log-Gabor filter is used to calculate the maximum response map corresponding to the phase consistency map. Then, the Canny edge detection algorithm is used to obtain the edge map corresponding to the maximum response map as the phase feature map. The original image is Gaussian filtered as the spatial feature map. The obtained phase feature map is combined with the spatial feature map as the guide map for joint bilateral filtering.

Images