Rotational scale invariant multi-modal image matching method and system based on directional filtering

Abstract

The application discloses a kind of rotation scale invariant multimodal image matching method and system based on directional filtering, by constructing a new type of directional filter, and based on the direction and amplitude information of the filter, a descriptor is constructed which is insensitive to rotation between multimodal images, by performing brute force matching between each layer of pyramid image to eliminate the influence of the difference in resolution between images on matching, and then realizing the matching of multi-modal, multi-angle, multi-scale remote sensing images.The method disclosed by the application is suitable for matching between various modal images such as optical-optical, optical-SAR, optical-depth map, optical-DEM, optical-rasterized vector map, optical-infrared, optical-luminous and the like.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a rotation scale-invariant multimodal image matching method and system based on directional filtering. Background Technology

[0002] With the continuous development of optical, infrared, and microwave sensors, and the emergence of various image products such as DEM, DOM, and vector images, the fusion and complementary advantages of multi-source data have become a current trend. For example, in the field of remote sensing, optical images possess rich spectral information that can provide useful texture information for ground feature identification. However, because optical images utilize visible light imaging, they are susceptible to interference from factors such as clouds and fog, often resulting in obscured ground features, which is detrimental to ground monitoring. Compared to optical images, microwave sensors are active sensors with strong penetration capabilities through clouds and fog, enabling all-weather Earth observation. However, because they are based on distance imaging, they lack fine texture information, often making them difficult to use for ground feature interpretation. Therefore, if optical and microwave images can be effectively fused, they can provide fine texture information while achieving all-weather Earth observation. In addition, there is the elevation information from DEM data, the high-precision planar information from DOM, and the surface temperature information from infrared images. If these information can be effectively fused to form data with the advantages of multiple sources, it will significantly improve humanity's ability to acquire information about the ground. Of course, multimodal image fusion not only plays a huge role in the field of remote sensing, but also plays an extremely important role in many other fields such as medicine, computer vision, and autonomous driving.

[0003] Multimodal matching technology is crucial for achieving multi-source data fusion. However, existing research primarily focuses on matching between visible light images, which is highly sensitive to differences in geometry, radiometry, and noise between multimodal images. Therefore, accurate matching between multimodal images remains unresolved, especially when significant angular and scale differences exist, making matching even more challenging. Thus, there is an urgent need to develop a multimodal image matching method suitable for images with angular and scale differences. Summary of the Invention

[0004] This invention addresses the problem that traditional image matching algorithms cannot accurately match multimodal images with differences in angle and scale. It proposes a rotation-scale-invariant multimodal image matching method and system based on directional filters.

[0005] This invention provides a rotation-scale-invariant multimodal image matching method based on directional filtering, comprising the following steps:

[0006] Step 1: Construct pyramids for the original image and the reference image respectively;

[0007] Step 2: Calculate the maximum phase consistency moment map for each pyramid layer of the original image and the reference image;

[0008] Step 3: Use the FAST feature point extraction algorithm to extract feature points from the maximum phase consistency matrix of each pyramid layer in the original image and the reference image.

[0009] Step 4: Construct an MDMS filter that simultaneously possesses scale and orientation information;

[0010] Step 5: Apply the MDMS filter constructed in Step 4 to the original image and the reference image to obtain the response values ​​of each pixel in the original image and the reference image at different scales and in different directions.

[0011] Step 6: Accumulate the response values ​​of each pixel in the original image / reference image at different scales in the same direction to eliminate the influence of scale. The accumulation result is used as the amplitude at the position of each pixel in the original image / reference image to obtain the orientation-amplitude map of the original image / reference image.

[0012] Step 7: Based on the amplitude fitting calculation of each pixel position in the original image / reference image obtained in Step 6, calculate the precise direction of each pixel in the original image / reference image to obtain the precise direction-amplitude map of the original image / reference image.

[0013] Step 8: Based on the accurate direction-magnitude map, calculate the principal direction of each feature point using a linear interpolation method;

[0014] Step 9: Construct feature descriptors for feature points;

[0015] Step 10: Perform traversal matching between different layers of the pyramid images of the original image and the reference image;

[0016] Step 11: Use the RANSAC algorithm to remove mismatched point pairs from the matching results obtained in Step 10.

[0017] Furthermore, the method for calculating the maximum moment of phase consistency of the pyramid image in step 2 is as follows:

[0018]

[0019]

[0020]

[0021]

[0022] In the formula, M ψ The maximum moment of phase consistency for any image point on the pyramid image at the current scale. This represents the sum of the squares of the phase coherence components in the x-direction from different directions. This represents the sum of the products of the x-direction and y-direction components of phase coherence in different directions. PC2(θ) represents the sum of the squares of the phase consistency components in the y-direction across different directions. It refers to the phase consistency of the current-scale 2D pyramid image at direction θ, calculated using the following formula:

[0023]

[0024] In the formula, W(θ) is a weighting factor, which assigns larger weights to multiple frequencies with a high degree of consistency; A n ξ is the amplitude at the nth scale, which represents the effective scale of the filter; ξ is a minimum value used to prevent division by zero. For noise compensation; The symbol indicates that the value is set to zero when it is less than 0, and is not processed when it is greater than 0; ΔΦ n (θ) represents the phase deviation function with scale n and direction θ, and is calculated as follows:

[0025]

[0026] In the formula, Φ n (θ) represents the phase component with scale n and direction θ, where θ represents the direction used by the directional filter. It represents the average value of phase components in different directions with a scale of n, and |·| represents taking the absolute value.

[0027] The maximum phase consistency moment is calculated for each pixel in each layer of the pyramid image of the original image and the reference image, respectively, to obtain the maximum phase consistency moment map of each layer of the pyramid image of the original image and the reference image.

[0028] Furthermore, the MDMS filter is constructed in step 4 as follows:

[0029] Lg = Lg s ·Lg o (7)

[0030]

[0031] Lg o =(cos(β-a)+1) / 2 (9)

[0032] In the formula, Lg represents the MDMS filter; Lg sFor a filter containing scale information; e is the natural base; r represents the radius from each pixel to the image center; w1 represents the wavelength of the minimum scale filter; μ represents the scale factor; s represents the current scale; σ represents the Gaussian standard deviation; Lg o The diffusion function is represented by β = 360 / N, where N is the number of directions specified by the user; and a is the polar angle of each image point in the polar coordinate system.

[0033] Furthermore, in step 5, N1 directions and N2 scales of the MDMS filter are selected to form N1×N2 filters. These filters are then multiplied by the original image / reference image after Fourier transform, and then an inverse Fourier transform is performed to calculate N1×N2 response values ​​corresponding to different directions and scales for each pixel of the original image / reference image.

[0034] Furthermore, in step 7, it is assumed that the direction of the maximum amplitude is θ, corresponding to an amplitude of A, and that its adjacent directions are θ1 and θ2, corresponding to amplitudes A1 and A2, respectively. Then:

[0035] d θ =0.5(A2-A1) / (A1+A2-2A) (10)

[0036]

[0037] θ max =θ1+d θ (12)

[0038] In the formula, N1 represents the number of directions of the filter; d θ This indicates that the direction θ1 is related to the direction of maximum fitted amplitude θ. max The offset; θ max This indicates the direction of the maximum amplitude after fitting.

[0039] θ max This refers to the precise direction of the fitted image points. A precise direction-amplitude map is constructed using the precise direction and maximum amplitude of each fitted image point.

[0040] Furthermore, in step 8, all feature points extracted in step 3 are traversed, and a square region with a side length of δ1 is extracted on the precise direction-amplitude map centered on each feature point. The amplitudes in the same direction within the square region are accumulated, and the direction with the highest accumulated value is taken as the principal direction of the corresponding feature point. To further improve the estimation accuracy, linear interpolation is performed using local statistical values. The specific calculation method is as follows: Assume the principal direction is θ. * The corresponding amplitude is A * And its adjacent directions are θ3 and θ4, with corresponding amplitudes A3 and A4, respectively. Then we have:

[0041]

[0042]

[0043]

[0044] In the formula, N3 represents the number of directions used by the filter. Indicates the direction θ3 and the principal direction after linear interpolation. The offset, This indicates the principal direction after linear interpolation.

[0045] Furthermore, in step 9, a square region with a side length of δ2 is extracted from the precise direction-amplitude map centered on the feature point and rotated by the angle of the main direction. The rotated square region is further subdivided into N4×N4 small grids. The amplitudes in the same direction in each small grid are accumulated. The accumulated amplitudes in each direction are merged and represented in vector form. This vector is the descriptor corresponding to each small grid. By merging the descriptors of the N4×N4 small grids, the rotation-invariant feature descriptor corresponding to the feature point for matching is obtained.

[0046] Furthermore, in step 10, the first layer of the original image pyramid is sequentially matched with the first, second, and third layers of the reference image pyramid. Then, the second layer of the original image pyramid is sequentially matched with the first, second, and third layers of the reference image pyramid. Finally, the third layer of the original image pyramid is sequentially matched with the first, second, and third layers of the reference image pyramid, ultimately obtaining initially matched feature point pairs. Brute-force matching refers to calculating the Euclidean distance between the descriptor of each feature point in each layer of the original image pyramid and the descriptor of each point in each layer of the reference image pyramid, and selecting the point with the smallest Euclidean distance as the matching point. Assume the feature point A of the original image pyramid has a descriptor F... A =(a1,a2,a3,…a n ), the feature point descriptor F of the pyramid image of the reference image (B-descriptor). B = (b1,b2,b3,…b n The Euclidean distance is calculated as follows:

[0047]

[0048] In the formula, D(F) A ,F B ) represents the Euclidean distance between two descriptors.

[0049] Furthermore, in step 11, assuming that the matching point set obtained in step 10 is X and Y, a random sampling consistency strategy based on an affine model is used to select the optimal point set as the correctly matched point set. The specific method is as follows:

[0050] Randomly select three or more pairs of points with the same name from point sets X and Y, and match the point pairs (x, y, y). i ,y i ) and (u i ,v i The coordinates are represented as homogeneous coordinates and expressed in matrix form as follows:

[0051]

[0052] In the formula, X and Y represent the sets of matching points in the source image and the target image, respectively;

[0053] According to the formula H = Y T X(X T X) -1 Calculate the least squares solution to obtain the transformation matrix H. Use the calculated affine transformation matrix H to map all matching points X to obtain the predicted matching point position X′. Calculate the Euclidean distance L between the predicted matching point position X′ and the actual matching point X. If L is less than the threshold η, it is considered an interior point; if L is greater than the threshold η, it is considered an exterior point. Count the total number of interior points contained in X. Repeat the above operation and take the set of data with the most interior points as the best dataset. The point pairs contained in this dataset are the corresponding points correctly matched by this invention.

[0054] The present invention also provides a rotation scale invariant multimodal image matching system based on directional filtering, for implementing the rotation scale invariant multimodal image matching method based on directional filtering as described above.

[0055] Furthermore, it includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the stored instructions in the memory to execute a rotation scale-invariant multimodal image matching method based on directional filtering as described above.

[0056] Compared with the prior art, the present invention has the following advantages:

[0057] 1) Unlike traditional gradient-based matching methods, this invention constructs a novel multi-directional filter. The response values ​​of filters with different directions are used to fit the direction of each pixel, and these response values ​​are used as amplitudes in descriptor construction. Descriptors constructed based on this filter can resist radiation, noise, and scale differences between multimodal images. 2) This invention employs a pyramid-based layer-by-layer matching strategy. Since the descriptor has some resistance to different scales, when the sampling rate between pyramid layers is small, a layer in the original image may successfully pair with feature points in two layers of the reference image, thus effectively increasing the number of correctly matched points. Attached Figure Description

[0058] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0059] Figure 1 This is a flowchart of the rotation scale-invariant multimodal image matching method based on directional filtering proposed in this invention.

[0060] Figure 2(a) shows the effect of image matching for microwave-optical image pairs with angular differences using the method proposed in this invention. Figure 2(b) shows the effect of image matching for vector-optical image pairs with angular and scale differences using the method proposed in this invention. Figure 2(c) shows the effect of image matching for optical-optical image pairs with angular differences using the method proposed in this invention. Figure 2(d) shows the effect of image matching for nighttime-daytime image pairs with scale differences using the method proposed in this invention. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0062] Example 1

[0063] like Figure 1 As shown, this embodiment of the invention provides a rotation-scale-invariant multimodal image matching method based on directional filtering, comprising the following steps:

[0064] Step 1: Construct pyramids for the original image and the reference image respectively.

[0065] Using one image from the multimodal image pair to be matched as the original image and the other image as the reference image, a pyramid is constructed by downsampling both the original and reference images. In this embodiment, a three-layer pyramid is constructed, requiring two downsampling operations. The first layer is the original image, the second layer is sampled to 3 / 4 of the original image, and the third layer is sampled to 1 / 2 of the original image. The nearest neighbor sampling strategy is used for the sampling method.

[0066] Step 2: Calculate the maximum phase consistency matrix for each pyramid layer of the original and reference images.

[0067] The method for calculating the maximum moment of phase consistency in pyramid images is as follows:

[0068]

[0069]

[0070]

[0071]

[0072] In the formula, M ψ The maximum moment of phase consistency for any image point on the pyramid image at the current scale. This represents the sum of the squares of the phase coherence components in the x-direction from different directions. This represents the sum of the products of the x-direction and y-direction components of phase coherence in different directions. PC2(θ) represents the sum of the squares of the phase consistency components in the y-direction across different directions. It refers to the phase consistency of the current-scale 2D pyramid image at direction θ, calculated using the following formula:

[0073]

[0074] In the formula, W(θ) is a weighting factor, which assigns larger weights to multiple frequencies with a high degree of consistency; A n ξ is the amplitude at the nth scale, which represents the effective scale of the filter; ξ is a minimum value used to prevent division by zero. For noise compensation; The symbol indicates that the value is set to zero when it is less than 0, and is not processed when it is greater than 0; ΔΦ n (θ) represents the phase deviation function with scale n and direction θ, and is calculated as follows:

[0075]

[0076] In the formula, Φ n (θ) represents the phase component with scale n and direction θ, where θ represents the direction used by the directional filter. It represents the average value of phase components in different directions with a scale of n, and |·| represents taking the absolute value.

[0077] The maximum phase consistency moment is calculated for each pixel in each layer of the pyramid image of the original image and the reference image, respectively, to obtain the maximum phase consistency moment map of each layer of the pyramid image of the original image and the reference image.

[0078] Step 3: Use the FAST feature point extraction algorithm to extract feature points from the maximum phase consistency matrix of each pyramid layer in the original image and the reference image.

[0079] Step 4: Construct an MDMS filter that has both scale and orientation information.

[0080] The MDMS filter (Multi-Directional Multi-Scale Filter) is constructed as follows:

[0081] Lg = Lg s ·Lg o (7)

[0082]

[0083] Lg o =(cos(β-a)+1) / 2 (9)

[0084] In the formula, Lg s For a filter containing scale information; e is the natural base; r represents the radius from each pixel to the image center; w1 represents the wavelength of the minimum scale filter, which is taken as an empirical value of 3 in this embodiment; μ represents the scale factor, which is taken as an empirical value of 2.1 in this embodiment; s represents the current scale, which is generally taken as a natural integer; σ represents the Gaussian standard deviation, which is taken as an empirical value of 0.55 in this embodiment; Lg o β represents the diffusion function; β = 360 / N, where N is the number of directions specified by the user. In this embodiment, N is set to 8, meaning 8 directions are used; a represents the polar angle of each image point in the polar coordinate system; Lg represents the MDMS filter, which is a directional filter formed by combining a filter containing scale information and a diffusion function. The combined filter has both scale and direction information.

[0085] Step 5: Apply the MDMS filter constructed in Step 4 to the original image and the reference image to obtain the response values ​​of each pixel in the original image and the reference image at different scales and in different directions.

[0086] N1 directions and N2 scales from the MDMS filter are selected to form N1×N2 filters. These filters are then multiplied by the original image / reference image after Fourier transform, and then an inverse Fourier transform is performed to calculate the N1×N2 response values ​​corresponding to each pixel in the original image / reference image with different directions and scales.

[0087] Step 6: Accumulate the response values ​​of each pixel in the original image / reference image at different scales in the same direction to eliminate the influence of scale. The accumulation result is used as the amplitude at the position of each pixel in the original image / reference image to obtain the N1-axis-amplitude map of the original image / reference image.

[0088] Step 7: Based on the amplitude values ​​at each pixel location in the original image / reference image calculated in Step 6, the precise orientation of each pixel in the original image / reference image is obtained by fitting the amplitude values.

[0089] Assume the direction of the maximum amplitude is θ, corresponding to an amplitude of A, and its adjacent directions are θ1 and θ2, corresponding to amplitudes A1 and A2 respectively, then:

[0090] d θ =0.5(A2-A1) / (A1+A2-2A) (10)

[0091]

[0092] θ max =θ1+d θ (12)

[0093] In the formula, N1 represents the number of directions of the filter; d θ This indicates that the direction θ1 is related to the direction of maximum fitted amplitude θ. max The offset; θ max This indicates the direction of the maximum amplitude after fitting.

[0094] θ max This refers to the precise direction of the fitted image points. A precise direction-amplitude map is constructed using the precise direction and maximum amplitude of each fitted image point.

[0095] Step 8: Based on the accurate direction-amplitude map, calculate the principal direction of each feature point using a linear interpolation method.

[0096] Traverse all feature points extracted in step 3, and extract a square region with a side length of δ1 on the precise direction-amplitude map centered on each feature point. Accumulate the amplitudes of the same direction within the square region, and take the direction with the highest accumulated value as the principal direction of the corresponding feature point. In this embodiment, direction histograms from 0 to 360° are statistically analyzed at 10-degree intervals. The principal direction accuracy is 10 degrees at this point. However, to further improve the estimation accuracy, linear interpolation is performed using local statistical values ​​to achieve an angle accuracy within 1°. The calculation method is as follows: Assume the principal direction is θ. * The corresponding amplitude is A * And its adjacent directions are θ3 and θ4, with corresponding amplitudes A3 and A4, respectively. Then we have:

[0097]

[0098]

[0099]

[0100] In the formula, N3 represents the number of directions used by the filter. Since the interval is 10 degrees, N3 is 360 / 10 = 36. Indicates the direction θ3 and the principal direction after linear interpolation. The offset; This indicates the principal direction after linear interpolation.

[0101] Step 9: Construct feature descriptors for feature points.

[0102] Centered on the feature point, a square region with a side length of δ2 is extracted on the precise direction-amplitude map and rotated by the main direction (in this embodiment, the square region of each feature point is rotated clockwise by the main direction). The rotated square region is further subdivided into N4×N4 small grids. The amplitudes in the same direction in each small grid are accumulated. The accumulated amplitudes in each direction are merged and represented in vector form. This vector is the descriptor corresponding to each small grid. By merging the descriptors of the N4×N4 small grids, the rotation-invariant feature descriptor corresponding to the feature point for matching is obtained.

[0103] Step 10: Perform traversal matching between different layers of the pyramid images of the original image and the reference image.

[0104] The process involves brute-force matching of the first layer of the original pyramid image with the first, second, and third layers of the reference pyramid image, then again, and finally, the third layer of the original pyramid image with the first, second, and third layers of the reference pyramid image. This process yields initial matched feature point pairs. Brute-force matching involves calculating the Euclidean distance between the descriptor of each feature point in each layer of the original pyramid image and the descriptor of each point in each layer of the reference pyramid image, and selecting the point with the smallest Euclidean distance as the matched point.

[0105] Assume the feature point descriptor F of the original pyramid image is A. A =(a1,a2,a3,…a n ), the feature point descriptor F of the pyramid image of the reference image (B-descriptor). B = (b1,b2,b3,…b n The Euclidean distance is calculated as follows:

[0106]

[0107] In the formula, D(F) A ,F B ) represents the Euclidean distance between two descriptors.

[0108] Step 11: Use the RANSAC algorithm to remove mismatched point pairs from the matching results obtained in Step 10.

[0109] Assuming the matching point set obtained in step 10 is X and Y, the optimal point set is selected as the correctly matched point set using the Random Sampling Consensus (RANSAC) strategy based on an affine model. The specific method is as follows:

[0110] Randomly select three pairs of points with the same name from point sets X and Y, and solve for the affine transformation matrix H using the least squares method. Then, match the points (x... i ,y i ) and (u i ,v i The coordinates are represented as homogeneous coordinates and expressed in matrix form as follows:

[0111]

[0112] In the formula, X and Y represent the sets of matching points in the source image and the target image, respectively.

[0113] According to the formula H = Y T X(X T X) -1Calculate the least squares solution to obtain the transformation matrix H. Use the calculated affine transformation matrix H to map all matching points X to obtain the predicted matching point positions X′. Calculate the Euclidean distance L between the predicted matching point position X′ and the actual matching point X. If L is less than the threshold η (η is 3 in this embodiment), it is considered an interior point; if L is greater than the threshold η, it is considered an exterior point. Count the total number of interior points contained in X. Repeat the above operation, and take the set of data with the most interior points as the optimal dataset. The point pairs contained in this dataset are the correctly matched corresponding points obtained by this invention. To improve the robustness of the algorithm, this invention sets the number of iterations to 10,000.

[0114] To verify the reliability of the matching method proposed in this invention, a total of 45,668 image pairs from publicly available WHU-OPT-SAR, OS, LEVIR-CD, LLVIP, MRSIs, and multimodal image datasets were used for testing. For ease of analysis, this invention divides these data into 10 groups based on different factors such as size, category, and dataset. The specific parameters of each group are shown in Table 1.

[0115] Table 1 Parameters of each group of data

[0116] 1 Optical-SAR 2781 600×600 2 SAR-optics 1208 512×512 3 Optics - Optics 1335 1024×1024 4 Optical-Depth Map 3000 1280×1024 5 SAR-optics 29220 256×256 6 SAR-optics 7305 512×512 7 Optical Elevation Map 30 600×600 8 Optical Point Cloud 63 600×600 9 Optical - Vector Graphics 669 400×400 10 Daytime-Nighttime Light Data 57 500×500

[0117] Figure 2 shows the matching results of four combinations of multimodal image pairs. As can be seen from Figure 2, the method proposed in this invention can achieve good matching results for multimodal images of different categories. Statistical analysis of all experimental data matching results shows that the method proposed in this invention achieves a matching success rate of 94.716% for multimodal images with different rotation angles and scales, with a matching accuracy of 1.94 pixels.

[0118] Example 2

[0119] Based on the same inventive concept, the present invention also provides a rotation scale invariant multimodal image matching system based on directional filtering, including a processor and a memory. The memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the rotation scale invariant multimodal image matching method based on directional filtering as described above.

[0120] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0121] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to replace them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A rotation-scale-invariant multimodal image matching method based on directional filtering, characterized in that, Includes the following steps: Step 1: Construct pyramids for the original image and the reference image respectively; Step 2: Calculate the maximum phase consistency moment map for each pyramid layer of the original image and the reference image; Step 3: Use the FAST feature point extraction algorithm to extract feature points from the maximum phase consistency matrix of each pyramid layer in the original image and the reference image. Step 4: Construct an MDMS filter that simultaneously possesses scale and orientation information; Step 5: Apply the MDMS filter constructed in Step 4 to the original image and the reference image to obtain the response values ​​of each pixel in the original image and the reference image at different scales and in different directions. Selecting MDMS filters One direction and Combination of scales forms Each of the following filters is used to perform a dot product with the original / reference image after Fourier transform, followed by an inverse Fourier transform to calculate the pixel value of each image in the original / reference image. Each response value corresponds to a different direction and a different scale; Step 6: Accumulate the response values ​​of each pixel in the original image / reference image at different scales in the same direction to eliminate the influence of scale. The accumulation result is used as the amplitude at the position of each pixel in the original image / reference image to obtain the orientation-amplitude map of the original image / reference image. Step 7: Based on the amplitude fitting calculation of each pixel position in the original image / reference image obtained in Step 6, calculate the precise direction of each pixel in the original image / reference image to obtain the precise direction-amplitude map of the original image / reference image. Step 8: Based on the accurate direction-magnitude map, calculate the principal direction of each feature point using a linear interpolation method; Step 9: Construct feature descriptors for feature points; Step 10: Perform traversal matching between different layers of the pyramid images of the original image and the reference image; Step 11: Use the RANSAC algorithm to remove mismatched point pairs from the matching results obtained in Step 10.

2. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: The method for calculating the maximum moment of phase consistency of the pyramid image in step 2 is as follows: （1） （2） （3） （4） In the formula, The maximum moment of phase consistency for any image point on the pyramid image at the current scale. Indicating phase coherence in different directions x The sum of the squares of the directional components. Indicating phase coherence in different directions x Directional components and y The sum of the products of the directional components. Indicating phase coherence in different directions y The sum of the squares of the directional components. This refers to the orientation of the current scale two-dimensional pyramid image. The phase consistency at a given point is calculated using the following formula: （5） In the formula, As a weighting factor, its function is to assign larger weights to multiple frequencies that have a high degree of consistency; It is the first n The magnitude of a scale, which represents the effective scale of the filter; It is used to prevent extremely small values ​​where the divisor is 0; For noise compensation; The symbol indicates that the value is set to zero when it is less than 0, and is not processed when it is greater than 0; To show the scale is n , direction is The phase deviation function is calculated as follows: （6） In the formula, The scale is represented as n、 Direction is phase components, Indicates the direction used by the directional filter. The scale is represented as n The average value of phase components in different directions, Indicates taking the absolute value; The maximum phase consistency moment is calculated for each pixel in each layer of the pyramid image of the original image and the reference image, respectively, to obtain the maximum phase consistency moment map of each layer of the pyramid image of the original image and the reference image.

3. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: The MDMS filter is constructed in step 4 as follows: （7） （8） （9） In the formula, Indicates the MDMS filter; A filter containing scale information; e The base is the natural number; r This represents the radius from each image point to the center of the image; The wavelength represents the minimum-scale filter; μ Indicates the scale factor; s Indicates the current scale; Indicates the Gaussian standard deviation; Represents the diffusion function; , N The number of directions is manually specified; The polar angle of each image point in the polar coordinate system.

4. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: In step 7, it is assumed that the direction of the maximum amplitude is... The corresponding amplitude is And their adjacent directions are respectively , The corresponding amplitudes are respectively , ,but: （10） （11） （12） In the formula, Indicates the number of directions of the filter; Indicates direction The direction of the maximum amplitude after fitting The offset; Indicates the direction of maximum amplitude after fitting; This refers to the precise direction of the fitted image points. A precise direction-amplitude map is constructed using the precise direction and maximum amplitude of each fitted image point.

5. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: In step 8, all feature points extracted in step 3 are traversed, and for each feature point, a side length of [missing value] is extracted on the precise orientation-amplitude map. Within a square region, the amplitudes in the same direction are accumulated, and the direction with the highest accumulated value is taken as the principal direction of the corresponding feature point. To further improve the estimation accuracy, linear interpolation is performed using local statistical values. The specific calculation method is as follows: Assume the principal direction is... The corresponding amplitude is And their adjacent directions are respectively , The corresponding amplitudes are respectively , Then we have: （13） （14） （15） In the formula, Indicates the number of directions used by the filter. Indicates direction With the principal direction after linear interpolation The offset, This indicates the principal direction after linear interpolation.

6. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: In step 9, using the feature point as the center, extract the side length from the precise direction-amplitude map. The square region is then rotated by an angle equal to the magnitude of the principal direction, and the rotated square region is further subdivided into... The process involves dividing the data into several small cells, summing the amplitude values ​​in the same direction within each cell, merging the sums from all directions, and representing them as a vector. This vector serves as the descriptor for each small cell. By combining the descriptions of each small cell, we obtain the rotation-invariant feature descriptors for matching the feature points.

7. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: In step 10, the first layer of the original pyramid image is sequentially matched with the first, second, and third layers of the reference pyramid image. Then, the second layer of the original pyramid image is sequentially matched with the first, second, and third layers of the reference pyramid image. Finally, the third layer of the original pyramid image is sequentially matched with the first, second, and third layers of the reference pyramid image, resulting in the initial matched feature point pairs. Brute-force matching involves calculating the Euclidean distance between the descriptor of each feature point in each layer of the original pyramid image and the descriptor of each point in each layer of the reference pyramid image, and selecting the point with the smallest Euclidean distance as the matching point. Assume the feature point A descriptor of the original pyramid image... Feature point B descriptor of the pyramid image of the reference image The Euclidean distance is calculated as follows: （16） In the formula, This represents the Euclidean distance between two descriptors.

8. The rotation-scale-invariant multimodal image matching method based on directional filtering as described in claim 1, characterized in that: In step 11, assuming the matching point set obtained in step 10 is X and Y, a random sampling consistency strategy based on an affine model is used to select the optimal point set as the correctly matched point set. The specific method is as follows: Randomly select three or more pairs of points with the same name from point sets X and Y, and match the point pairs. and Represented as homogeneous coordinates, and expressed in matrix form as follows: （17） In the formula, X and Y represent the sets of matching points in the source image and the target image, respectively; According to the formula Calculate the least squares solution to obtain the transformation matrix H. Use the calculated affine transformation matrix H to apply the transformation to all matching points. By performing mapping, the predicted matching point positions can be obtained. Calculate the predicted matching point location Matching points Euclidean distance between L ,like L Less than the threshold If , then it is considered an interior point. L Greater than the threshold If it is an outlier, it is considered an outlier and statistically analyzed. The total number of interior points is counted. The above operation is repeated, and the set of data with the most interior points is taken as the best dataset. The point pairs contained in this dataset are the corresponding points that are correctly matched by this invention.

9. A rotation-scale-invariant multimodal image matching system based on directional filtering, characterized in that, It includes a processor and a memory, the memory being used to store program instructions, and the processor being used to call the program instructions in the memory to execute the rotation scale-invariant multimodal image matching method based on directional filtering as described in any one of claims 1-8.

Images