Multimodal remote sensing image registration method based on consistent gradient features

Abstract

The present invention discloses a multi-mode remote sensing image registration method based on consistent gradient features, comprising: extracting a first gradient map of a first image and a second gradient map of a second image; determining a first feature histogram and a first feature point according to the first gradient map; determining a corresponding second feature point in the second image according to multiple pairs of preprocessing control points between the first image and the second image, and the first feature point; determining a first feature descriptor according to the first feature histogram; determining a second feature histogram according to the second gradient map, and determining a second feature descriptor according to the second feature histogram; determining a minimum position error through template matching and fast Fourier transform according to the first feature descriptor, the second feature descriptor, target feature points in the first feature points and the second image; obtaining third feature points according to the second feature points and the minimum position error, each third feature point forming a pair of matching points with a first feature point; and eliminating mismatching points from the matching points.

Description

Technical Field

[0001] The invention belongs to the technical field of image processing, and in particular relates to a multi-mode remote sensing image registration method based on consistent gradient features. Background Art

[0002] Image registration technology is an important part of the image processing field and has been widely used in military and civilian fields such as aircraft navigation systems, missile positioning, urban development monitoring, terrain and mountain mapping, biological diagnosis, climatology, etc. Heterogeneous image registration technology can make up for the shortcomings of single sensor images in image registration technology and achieve the same geometric consistency between the reference image and the image to be registered. It is the basis for processing multi-source remote sensing sensor images and is also a research hotspot in the current field of image registration technology.

[0003] At present, multimodal remote sensing image registration methods can be divided into traditional registration algorithms and deep learning registration algorithms. Traditional multimodal remote sensing image registration algorithms are mainly represented by image registration methods based on grayscale information and image registration methods based on feature information. Registration algorithms based on deep learning are divided into supervised learning registration algorithms and unsupervised learning registration algorithms.

[0004] However, due to the nonlinear grayscale differences and radiation intensity differences (illumination differences) between multimodal remote sensing images, image registration methods based on grayscale information are mostly aimed at homologous images, which are sensitive to illumination differences and grayscale changes. They are difficult to process multimodal remote sensing images, have low matching accuracy, severely reduced reliability, and high time complexity. Although image registration methods based on feature information can stably process multi-source remote sensing images, due to the nonlinear radiation distortion that may exist in remote sensing images, using intensity information or gradient information for feature detection and description will cause these algorithms to fail, making it difficult to achieve stable registration results; and using phase information for feature detection and description matching has a high time complexity, and the accuracy of extracted feature points is low. Image registration algorithms based on deep learning require a large number of sample data sets, have poor transferability, and unstable registration results.

[0005] In other words, the current multi-modal remote sensing image registration method is difficult to meet the requirements of high precision and low latency for multi-modal remote sensing image registration. Summary of the invention

[0006] In order to solve the above problems existing in the related art, the present invention provides a multi-modal remote sensing image registration method based on consistent gradient features. The technical problem to be solved by the present invention is achieved through the following technical solutions:

[0007] The present invention provides a multi-mode remote sensing image registration method based on consistent gradient features, comprising:

[0008] Extracting a first gradient map of the first image and a second gradient map of the second image respectively;

[0009] Determine a first gradient consistency feature histogram and a plurality of first feature points according to the first gradient map;

[0010] Determine second position information of a corresponding second feature point in the second image according to a plurality of pairs of preprocessed control points between the first image and the second image and first position information of each first feature point;

[0011] Determining a first feature descriptor according to the first gradient consistency feature histogram;

[0012] Determine a second gradient consistency feature histogram according to the second gradient map, and determine a second feature descriptor according to the second gradient consistency feature histogram;

[0013] Determine a minimum position error by template matching and fast Fourier transform according to the first feature descriptor, the second feature descriptor, the target feature point among the plurality of first feature points, and the second image;

[0014] According to the second position information and the minimum position error, third position information is obtained; each third position information corresponds to a third feature point, and a third feature point and a one-to-one corresponding first feature point constitute a pair of matching points;

[0015] Each pair of mismatched points is removed from multiple pairs of matching points to obtain the registration result.

[0016] The present invention has the following beneficial technical effects:

[0017] By using the gradient information (gradient map) of the image to extract the feature point information of the image, it is not affected by the grayscale change and illumination difference between the heterogeneous images (the first image and the second image), so that the robust and stable feature points can be obtained; by extracting the gradient map of each image, and generating the gradient consistency feature histogram of the image according to the gradient map of the image, the influence of SAR speckle noise interference on the registration result can be reduced; by performing fast similarity measurement processing through template matching and fast Fourier transform to obtain the minimum position error, the operation efficiency can be improved, so that multiple pairs of accurately positioned matching points can be quickly obtained by quickly obtaining the position error; and by further proposing mismatching points from the obtained multiple pairs of accurately positioned matching points, high-precision matching points can be obtained. Therefore, the registration method provided by the present invention can make the registration of multi-source images more accurate and efficient.

[0018] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 An optional flow chart of a multi-modal remote sensing image registration method based on consistent gradient features provided in an embodiment of the present invention;

[0020] Figure 2A An exemplary registration result connection diagram obtained according to the registration result of the visible light image and the SAR image provided in an embodiment of the present invention;

[0021] Figure 2B An exemplary fusion image obtained according to the registration result of a visible light image and a SAR image provided in an embodiment of the present invention;

[0022] Figure 2C An exemplary chessboard mosaic obtained according to the registration result of a visible light image and a SAR image provided in an embodiment of the present invention;

[0023] Figure 3A An exemplary registration result connection diagram obtained according to the registration result of the visible light image and the infrared image provided in the embodiment of the present invention;

[0024] Figure 3B An exemplary fusion image obtained according to the registration result of a visible light image and an infrared image provided in an embodiment of the present invention;

[0025] Figure 3C An exemplary chessboard mosaic obtained according to the registration result of the visible light image and the infrared image provided in the embodiment of the present invention;

[0026] Figure 4A An exemplary registration result connection diagram obtained according to the registration results of two different SAR images provided in an embodiment of the present invention;

[0027] Figure 4B An exemplary fusion image obtained according to the registration results of two different SAR images provided in an embodiment of the present invention;

[0028] Figure 4C An exemplary chessboard mosaic image obtained according to the registration results of two different SAR images provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0029] The present invention is further described in detail below with reference to specific embodiments, but the embodiments of the present invention are not limited thereto.

[0030] In the description of the present invention, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0031] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine different embodiments or examples described in this specification.

[0032] Although the present invention is described herein in conjunction with various embodiments, in the process of implementing the claimed invention, those skilled in the art may understand and implement other variations of the disclosed embodiments by viewing the drawings, the disclosure, and the appended claims. In the claims, the word "comprising" does not exclude other components or steps, and "one" or "an" does not exclude multiple situations. A single processor or other unit may implement several functions listed in the claims. Certain measures are recorded in different dependent claims, but this does not mean that these measures cannot be combined to produce good results.

[0033] Figure 1 is a flow chart of a multi-modal remote sensing image registration method based on consistent gradient features provided by an embodiment of the present invention. Figure 1 As shown, the method comprises the following steps:

[0034] S101 , extracting a first gradient map of a first image and a second gradient map of a second image respectively.

[0035] Here, the first image and the second image may be acquired first, and then a first gradient map of the first image and a second gradient map of the second image may be extracted respectively.

[0036] Here, the first image and the second image may be heterogeneous images of the same area or object. For example, the first image may be a visible light image, and the second image may be a radar image.

[0037] Here, the first gradient map of the first image may be a gradient amplitude map (i.e., a consistent gradient structure map, hereinafter referred to as the first gradient amplitude map) and a gradient amplitude angle map (i.e., a consistent gradient structure map, hereinafter referred to as the first gradient amplitude angle map) of the first image. The second gradient map of the second image may be a gradient amplitude map (hereinafter referred to as the second gradient amplitude map) and a gradient amplitude angle map (hereinafter referred to as the second gradient amplitude angle map) of the second image.

[0038] Here, when there is no rotation difference and scaling difference between the first image and the second image, S101 can be directly started, and when there is a rotation difference or scaling difference between the first image and the second image, the first image and the second image can be pre-processed with geometric transformation, so as to achieve correction of the first image and the second image. When correcting the first image and the second image, a radiation invariant feature transformation matching algorithm can be used to directly perform a rough geometric transformation correction, or first pre-processing control points manually marked on the first image and corresponding second pre-processing control points marked on the second image can be obtained first, and a rough geometric transformation correction can be performed according to the obtained pre-processing control points. After obtaining the corrected first image and the corrected second image, S101 is started to perform registration between the corrected first image and the corrected second image.

[0039] In some embodiments, a first edge detection operator may be used to extract a first gradient magnitude map and a first gradient argument map of a first image, respectively; a second edge detection operator may be used to extract a second gradient magnitude map and a second gradient argument map of a second image, respectively; wherein the first edge detection operator and the second edge detection operator may be the same or different. For example, the first edge detection operator or the second edge detection operator may be a ROEWA operator or a Scharr operator.

[0040] In some embodiments, the first edge detection operator is a Sobel operator corresponding to four directional templates of 0°, 45°, 90° and 135°. Based on this, the first gradient amplitude map and the first gradient amplitude angle map can be obtained by executing the following S1011 to S1013.

[0041] S1011, using four directional templates of 0°, 45°, 90° and 135° to respectively calculate the first-order derivative of each first pixel point of the first image, and correspondingly obtain a first derivative value, a second derivative value, a third derivative value and a fourth derivative value.

[0042] Exemplarily, the four direction templates are as follows: Wherein, h1 and h2 represent the horizontal (0°) and vertical (90°) templates, respectively, and h3 and h4 represent the templates in two diagonal directions (45° and 135°), respectively. Correspondingly, the first derivative value, the third derivative value, the second derivative value, and the fourth derivative value can be respectively expressed as the following formulas (1) to (4):

[0043]

[0044]

[0045]

[0046]

[0047] Among them, D1(x1,y1), D2(x1,y1), D3(x1,y1), and D4(x1,y1) represent the first derivative value, the third derivative value, the second derivative value, and the fourth derivative value, respectively, and f(x1,y1) represents the pixel value of each first pixel point (x1,y1) in the first image.

[0048] Here, the traditional Sobel operator estimates the approximate value of the image gradient, using templates in the horizontal and vertical directions, but in practice the gradient direction is unknown, and the main body of the image and the background are not strictly distinguished. The traditional Sobel operator does not take into account the influence of the diagonal direction on the image, and the calculated result has a certain error. Therefore, on the basis of the original two-direction template, templates in the two diagonal directions of 45° and 135° can be added to obtain the gradient image, and edges in multiple directions can be detected. Therefore, the problem of limitations in only the horizontal and vertical directions can be solved to a certain extent, making the calculation result more accurate. Templates in four directions are defined here. Although continuing to increase the number of templates can further improve the accuracy of the calculation, considering the limitation of calculation efficiency and edge range, the present invention does not need to increase the number of templates here.

[0049] S1012. Calculate the gradient amplitude of each first pixel point according to the absolute values ​​of the first derivative value, the second derivative value, the third derivative value, and the fourth derivative value to obtain a first gradient amplitude map.

[0050] Exemplarily, the gradient magnitude of each first pixel point may be calculated using formula (5):

[0051]

[0052] Among them, |D1(x1,y1)|, |D2(x1,y1)|, |D3(x1,y1)|, and |D4(x1,y1)| represent the absolute values ​​of the first derivative value, the third derivative value, the second derivative value, and the fourth derivative value, respectively.

[0053] Here, the calculation formula of the gradient amplitude is defined as follows:

[0054]

[0055] However, the computational complexity of the square root sum is relatively high, and the method proposed in the present invention of calculating the gradient amplitude by summing up the absolute values ​​(Formula (5)) can make the number of processes relatively small, thereby effectively improving the computational efficiency.

[0056] S1013. Calculate the gradient amplitude of each first pixel point according to the first derivative value and the third derivative value to obtain a first gradient amplitude map.

[0057] Exemplarily, the gradient amplitude of each first pixel point may be calculated using the following formula (7):

[0058]

[0059] Wherein, φ represents the gradient magnitude angle of each first pixel.

[0060] Here, the gradient amplitude of all pixels of an image can be called the gradient amplitude map of the image.

[0061] In some embodiments, the second edge detection operator is a filter composed of a causal filter and a non-causal filter; based on this, the second gradient amplitude map and the second gradient amplitude angle map can be obtained by executing S1014 to S1017.

[0062] S1014. Calculate the causal exponential weighted mean and the non-causal exponential weighted mean in the horizontal direction and the causal exponential weighted mean and the non-causal exponential weighted mean in the vertical direction of each second pixel point of the second image according to the filter, the causal filter and the non-causal filter.

[0063] Here, the filter composed of causal filters and non-causal filters can also be called an improved ROEWA operator. The ROEWA operator is an exponentially weighted mean ratio operator that estimates the average value within a neighborhood window, and the change in the average value is positively correlated with the variance of the source image. Compared with the traditional gradient difference operator, the ROEWA operator can effectively overcome the influence of the difference between the strong scattering texture area and the low reflective texture area of ​​the SAR image on the result. In one-dimensional case, the filter expression is the formula: f(x) = Ae -αx (8); where A is a normalization constant, α is a filter coefficient, and x represents the pixel value of each pixel of the image. Under discrete conditions, the filter f can be obtained by combining the causal filter f1 and the non-causal filter f2. At this time, the filter can be discretized as follows: Where f1(n) = a·bn h(n) is a causal filter, f2(n) = a·b -n h(-n) is a non-causal filter, a and α are preset coefficients, n represents the pixel value of each pixel, 0<b=e -α <1, h(n) is the unit step function. Under discrete conditions, h(n)=1 when n≥0, otherwise h(n)=0.

[0064] Exemplarily, in the horizontal direction, the image intensity of the second image can be first filtered column by column using filter f to obtain the filtered second image, and then the filtered second image can be filtered row by row using causal filter f1 and non-causal filter f2, respectively, to obtain the causal exponential weighted mean and non-causal exponential weighted mean in the horizontal direction. Similarly, in the vertical direction, the image intensity of the second image can be first filtered row by row using filter f to obtain the filtered second image, and then the filtered second image can be filtered column by column using causal filter f1 and non-causal filter f2, respectively, to obtain the causal exponential weighted mean and non-causal exponential weighted mean in the vertical direction. For example, the causal exponential weighted mean and non-causal exponential weighted mean in the horizontal direction can be expressed by formulas (10) to (11):

[0065] μ x1 (x2,y2)=f1(x2,y2)*(f(y2)*I(x2,y2)) (10);

[0066] μ x2 (x2,y2)=f2(x2,y2)*(f(y2)*I(x2,y2)) (11);

[0067] Where I(x2,y2) represents the pixel value of each second pixel point (x2,y2) of the second image, μ x1 (x2, y2) represents the causal exponential weighted mean of the second pixel (x2, y2) in the horizontal direction, μ x2 (x2, y2) represents the non-causal exponentially weighted mean of the second pixel (x2, y2) in the horizontal direction, f is a filter, f1 is a causal filter, f2 is a non-causal filter, and * represents a convolution operation.

[0068] S1015. Calculate a gradient estimation value in the horizontal direction for each second pixel point according to the causal exponential weighted mean and the non-causal exponential weighted mean in the horizontal direction.

[0069] S1016. Calculate the estimated gradient value in the vertical direction of each second pixel point according to the causal exponential weighted mean and the non-causal exponential weighted mean in the vertical direction.

[0070] Exemplarily, the following formula (12) may be used to calculate the estimated gradient value in the horizontal direction of each second pixel, and the following formula (13) may be used to calculate the estimated gradient value in the vertical direction of each second pixel:

[0071]

[0072]

[0073] Among them, r x (x2, y2) represents the estimated value of the horizontal gradient of the second pixel (x2, y2), r y (x2, y2) represents the estimated gradient value in the vertical direction of the second pixel (x2, y2), (x2-1, y2) represents the neighboring pixel on the left side of the second pixel (x2, y2), (x2+1, y2) represents the neighboring pixel on the right side of the second pixel (x2, y2), μ y1 (x2-1, y2) represents the causal exponential weighted mean in the vertical direction of the second pixel point (x2-1, y2), μ y2 (x2+1, y2) represents the non-causal exponentially weighted mean in the vertical direction of the second pixel point (x2+1, y2).

[0074] S1017. Calculate the gradient amplitude and gradient angle of each second pixel point according to the gradient estimation value in the horizontal direction, the gradient estimation value in the vertical direction and the preset gradient difference coefficient, and obtain a second gradient amplitude map and a second gradient angle map.

[0075] Here, there are gradient differences between heterogeneous images. A preset gradient difference coefficient can be used to perform linear processing on the obtained causal exponential weighted mean and non-causal exponential weighted mean to obtain the processed causal exponential weighted mean and non-causal exponential weighted mean. The gradient estimation value in the horizontal direction and the gradient estimation value in the vertical direction are calculated according to the processed causal exponential weighted mean and non-causal exponential weighted mean to obtain the processed horizontal gradient estimation value and the processed vertical gradient estimation value. The processed horizontal gradient estimation value and the processed vertical gradient estimation value are used to calculate the gradient amplitude and gradient amplitude angle of each second pixel point respectively.

[0076] Exemplarily, the estimated gradient value in the horizontal direction after processing of each second pixel point can be expressed by the following formula (14), and the estimated gradient value in the vertical direction after processing can be expressed by the following formula (15):

[0077]

[0078]

[0079] Among them, q is the preset gradient difference coefficient, r x '(x2, y2) is the estimated horizontal gradient value of each second pixel (x2, y2) after processing, r y '(x2, y2) is the estimated gradient value in the vertical direction of each second pixel point (x2, y2) after processing.

[0080] Exemplarily, the second gradient magnitude of each second pixel point (x2, y2) can be expressed by formula (16), and the second gradient amplitude angle can be expressed by formula (17):

[0081]

[0082]

[0083] Among them, r(x2,y2) is the second gradient amplitude of each second pixel point (x2,y2), is the second gradient amplitude of each second pixel point (x2, y2).

[0084] S102: Determine a first gradient consistency feature histogram and a plurality of first feature points according to the first gradient map.

[0085] Here, histogram statistics can be performed on the first gradient amplitude map and the first gradient amplitude angle map to obtain a three-dimensional first gradient consistency feature histogram, where the three dimensions include: amplitude dimension, amplitude angle dimension and scale dimension; the scale dimension characterizes the number of first gradient amplitude angle maps; then, the first gradient amplitude map is subjected to non-maximum suppression processing to obtain multiple first position information, each first position information corresponds to a first feature point.

[0086] Here, the gradient consistency feature histogram uses the gradient amplitude as the X coordinate, the gradient amplitude as the Y coordinate, and the number of gradient amplitude maps as the Z coordinate. In the embodiment of the present invention, each image corresponds to a gradient amplitude map, so the scale can be 1. The gradient consistency feature histogram calculated in this way can overcome the influence of SAR speckle noise interference and has strong robustness.

[0087] Here, by performing non-maximum suppression processing on the first gradient amplitude map, feature points in the first image can be extracted, thereby obtaining a plurality of first feature points.

[0088] S103: Determine second position information of corresponding second feature points in the second image according to multiple pairs of preprocessed control points between the first image and the second image and the first position information of each first feature point.

[0089] In some embodiments, before executing S103, a first preprocessing control point in the first image and a second preprocessing control point in the second image can be obtained; wherein each first preprocessing control point corresponds one-to-one to a second preprocessing control point; the first preprocessing control point and the second preprocessing control point having a corresponding relationship are taken as a pair of preprocessing control points to obtain multiple pairs of preprocessing control points.

[0090] Here, multiple pairs of preprocessing control points can be extracted from the first image and the second image by using a radiation invariant feature transformation matching algorithm, or multiple pairs of preprocessing control points can be obtained based on manual annotations on the first image and the second image.

[0091] In some embodiments, the first projection transformation matrix between the first image and the second image can be determined based on multiple pairs of preprocessed control points between the first image and the second image, and then the second position information of the second feature point corresponding to the first feature point in the second image can be obtained based on the first position information of each first feature point and the first projection transformation matrix.

[0092] Here, the first projection transformation matrix is ​​an optimal projection transformation matrix determined according to multiple pairs of preprocessing control points.

[0093] S104: Determine a first feature descriptor according to the first gradient consistency feature histogram.

[0094] Here, the first gradient consistency feature histogram may be normalized to obtain a first feature descriptor (ie, gradient consistency descriptor). The gradient consistency descriptor is a three-dimensional vector.

[0095] S105 . Determine a second gradient consistency feature histogram according to the second gradient map, and determine a second feature descriptor according to the second gradient consistency feature histogram.

[0096] Similarly, the second gradient consistency feature histogram can be normalized to obtain a second feature descriptor. The second feature descriptor is also a three-dimensional vector.

[0097] S106 , determining a minimum position error through template matching and fast Fourier transform according to the first feature descriptor, the second feature descriptor, the target feature point among the plurality of first feature points, and the second image.

[0098] Here, the minimum position error can be determined based on the first feature descriptor, the second feature descriptor, the target feature point among the multiple first feature points and the second image through a template matching algorithm and a fast Fourier transform, wherein the template matching algorithm can be, for example, a sum of absolute differences (SAD) algorithm, a sum of squared differences (SSD) algorithm or a mean square difference (MSD) algorithm.

[0099] In some embodiments, the minimum position error may be determined according to the first feature descriptor, the second feature descriptor, the target feature point among the plurality of first feature points, and the second image through an MSD algorithm and a fast Fourier transform, as specifically described in steps S1 to S3:

[0100] S1. Using a window template of a preset size centered on a target first feature point, determine each search window from a region of interest in a second image; the target feature point is any one of the multiple first feature points.

[0101] Here, the size of the window template can be set according to actual needs. For example, it can be a window template of size 3*3, so that the size of each corresponding search window on the second image is also 3*3.

[0102] Here, the region of interest of the second image may be the entire second image or some preset regions in the second image. For example, when the region of interest of the second image is the entire second image, a window template may be used to start from the first pixel of the second image and perform a sliding window search on each pixel in sequence.

[0103] S2. According to the first feature descriptor and the second feature descriptor, the minimum value of the average sum of square errors between the window template and the search window is calculated by fast Fourier transform; each average sum of square errors is obtained based on a position error; each average sum of square errors corresponds to a search window.

[0104] S3. The position error corresponding to the minimum value is taken as the minimum position error.

[0105] Here, when the MSD method is used for similarity measurement, the specific formula is defined as follows:

[0106]

[0107] Where C1 and C2 represent the feature descriptor of the window template i and the feature descriptor of each search window respectively, p represents the coordinates of each first pixel point in the first image, T i() is a mask function, when p is the first pixel in the window template i, T i (p) = 1, otherwise T i (p) = 0, v represents the window template and position error. i (v) represents the position error, M i (v) is the MSD value between window template i and each search window. Since there are several search windows in a region of interest, M i (v) is in the form of a one-dimensional vector. By minimizing M i (v) The position of the registration window corresponding to the template window i can be found. Therefore, the matching function can be defined as follows:

[0108]

[0109] Substituting the matching function into formula (18), we can obtain the following formula:

[0110]

[0111] The first term in formula (19) is a constant, and the second and third terms are variables. Therefore, when the second and third terms are minimum, M i (v) is the smallest, then the corresponding v is the smallest, and the fast Fourier transform FFT can be used to accelerate the processing to obtain the smallest v. The specific formula is as follows:

[0112]

[0113] Among them, F(·), F -1 (·) and F * (·) represent the Fourier transform, inverse Fourier transform and complex conjugate of Fourier transform, respectively, where v i is the minimum position error.

[0114] Here, since it takes a long time to perform template matching using a sliding window operation, the MSD algorithm is used for similarity measurement, and the MSD algorithm is accelerated by using fast Fourier transform, which can greatly shorten the matching time and thus improve the efficiency of obtaining the minimum position error.

[0115] S107. Obtain third position information according to the second position information and the minimum position error; each piece of third position information corresponds to a third feature point, and a third feature point and a one-to-one corresponding first feature point constitute a pair of matching points.

[0116] Here, each second position information is actually a two-dimensional coordinate, for example, (x, y), and correspondingly, each third position information is also a two-dimensional coordinate, and the minimum position error is a two-dimensional position error, for example, it can be Δx and Δy. The sum of the x coordinate of each second position information and Δx can be used as the x coordinate in the third position information, and the sum of the y coordinate of each second position information and Δy can be used as the y coordinate in the third position information, so as to obtain a third position information corresponding to each second position information one by one, and the pixel point represented by each third position information is the third feature point in the second image, and the third feature point is the feature point with accurate position in the second image, and each third feature point and a first feature point corresponding to each other in the first image constitute a pair of matching points, so that multiple pairs of matching points can be obtained.

[0117] S108, eliminating each pair of mismatched points from the multiple pairs of matching points to obtain a registration result.

[0118] In some embodiments, an improved Fast Sample Consensus (FSC) algorithm or a Random Sampling Consensus (RANSAC) algorithm may be directly used to remove each pair of mismatched points from multiple pairs of matching points, thereby using the remaining matching points as the registration results.

[0119] Exemplarily, the principle of using the improved FSC algorithm to eliminate each pair of mismatched points from multiple pairs of matching points is as follows: (1) Randomly select m pairs of matching points from multiple pairs of matching points n between the first image and the second image as m sample data, and calculate the transformation matrix H; (2) Use the position information of all first feature points of the first image and the transformation matrix H to estimate the position information of the corresponding feature points in the second image. If the error between the feature point position of the second image and the estimated feature point position information is less than a threshold, add the feature point position to the inlier set I; (3) If the number of elements s in the current inlier set is greater than the optimal number of inlier sets I_best, update I_best=I and update the number of iterations k at the same time; (4) If the number of iterations is greater than k, exit; otherwise, increase the number of iterations by 1 and repeat the above steps.

[0120] In some embodiments, the second projection transformation matrix between the first image and the second image can be first determined based on multiple pairs of matching points; the estimated feature point corresponding to the first feature point can be determined based on the second projection transformation matrix and the first feature point in each pair of matching points; when the difference between the second position information of the third feature point corresponding to the first feature point and the position information of the estimated feature point is greater than or equal to a preset error threshold, the first feature point and the third feature point are eliminated to obtain multiple pairs of matching points after error detection; thereafter, each pair of mismatched points in the multiple pairs of matching points after error detection is continued to be eliminated to obtain the alignment result.

[0121] Here, after obtaining multiple pairs of matching points after error detection, the improved FSC algorithm or the RANSAC algorithm can also be used to remove each pair of mismatched points in the multiple pairs of matching points after error detection. The accuracy of the obtained registration result can be further improved by error detection estimation and removal of mismatched points.

[0122] In the embodiment of the present invention, by using the gradient information (gradient map) of the image to extract the feature point information of the image, it can be unaffected by the grayscale change and illumination difference between heterogeneous images (first image and second image), so that robust and stable feature points can be obtained; by extracting the gradient map of each image, and generating the gradient consistency feature histogram of the image according to the gradient map of the image, the influence of SAR speckle noise interference on the registration result can be reduced; by performing fast similarity measurement processing through template matching and fast Fourier transform to obtain the minimum position error, the operation efficiency can be improved, so that multiple pairs of accurately positioned matching points can be quickly obtained by quickly obtaining the position error; and by further proposing mismatching points from the obtained multiple pairs of accurately positioned matching points, high-precision matching points can be obtained. Therefore, the registration method provided by the present invention can make the registration of multi-source images more accurate and efficient.

[0123] The overall idea of ​​the registration method proposed in the present technical invention is different from the traditional heterogeneous image registration method that detects and describes features through the grayscale or phase information of the image. Instead, it fully extracts the consistent gradient information of the image, obtains feature points with high stability, and then obtains feature descriptors with a high degree of consistency. Then, the MSD algorithm is accelerated by using the fast Fourier transform. Finally, the improved FSC algorithm is used to remove the same-name points with large errors, and the root mean square error of the algorithm is reduced, thereby greatly improving the registration efficiency and the accuracy of the registration results.

[0124] The following further illustrates the registration effect of the registration method provided by the embodiment of the present invention through experimental data.

[0125] Figure 2A is a connection diagram of the registration results obtained based on the registration results of the visible light image and the SAR image; Figure 2B It is the fusion image obtained based on the registration results of the visible light image and the SAR image; Figure 2C This is a checkerboard mosaic image obtained based on the registration results of the visible light image and the SAR image.

[0126] Figure 3A is a connection diagram of the registration results obtained according to the registration results of the visible light image and the infrared image; Figure 3B It is a fusion image obtained based on the registration results of the visible light image and the infrared image; Figure 3C It is a checkerboard mosaic image obtained based on the registration results of the visible light image and the infrared image.

[0127] Figure 4A It is a connection diagram of the registration results obtained based on the registration results of two different SAR images; Figure 4B It is a fusion image obtained based on the registration results of two different SAR images; Figure 4C The checkerboard mosaic image is obtained based on the registration results of two different SAR images.

[0128] The CFOG algorithm and the registration method provided by the present invention are used to simulate and compare the above three pairs of images, and their root mean square error (RMSE) and running time are shown in Table 1. Obviously, the method provided by the present invention has a smaller root mean square error, better registration effect, and shorter running time, which can verify the effectiveness of the method proposed in this article.

[0129] CFOG The present invention RMSE / pixel Time / t RMSE / pixel Time / t The first pair of images 0.7679 4.5978 0.3814 3.3357 The second pair of images 0.3417 4.1095 0.1603 3.1981 The third pair of images 0.5457 3.6967 0.4015 2.6595

[0130] Table 1

[0131] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as falling within the scope of protection of the present invention.

Claims

1. A multi-modal remote sensing image registration method based on consistent gradient features, characterized in that: include: Extracting a first gradient map of the first image and a second gradient map of the second image respectively; Determine a first gradient consistency feature histogram and a plurality of first feature points according to the first gradient map; Determine second position information of a corresponding second feature point in the second image according to a plurality of pairs of preprocessed control points between the first image and the second image and first position information of each first feature point; Determining a first feature descriptor according to the first gradient consistency feature histogram; Determine a second gradient consistency feature histogram according to the second gradient map, and determine a second feature descriptor according to the second gradient consistency feature histogram; Determine a minimum position error by template matching and fast Fourier transform according to the first feature descriptor, the second feature descriptor, the target feature point among the plurality of first feature points, and the second image; According to the second position information and the minimum position error, third position information is obtained; each third position information corresponds to a third feature point, and a third feature point and a one-to-one corresponding first feature point constitute a pair of matching points; Each pair of mismatched points is removed from multiple pairs of matching points to obtain the registration result.

2. The multimodal remote sensing image registration method based on consistent gradient features according to claim 1, characterized in that: The first gradient map includes: a first gradient magnitude map and a first gradient amplitude angle map; and determining a first gradient consistency feature histogram and a plurality of first feature points according to the first gradient map includes: Performing histogram statistics on the first gradient amplitude map and the first gradient argument map to obtain a three-dimensional first gradient consistency feature histogram, wherein the three dimensions include: an amplitude dimension, an argument dimension, and a scale dimension; the scale dimension represents the number of the first gradient argument maps; Non-maximum suppression processing is performed on the first gradient amplitude map to obtain a plurality of first position information, each of which corresponds to a first feature point.

3. The multimodal remote sensing image registration method based on consistent gradient features according to claim 1 or 2, characterized in that: The determining, according to the plurality of pairs of preprocessed control points between the first image and the second image and the first position information of each first feature point, the second position information of the corresponding second feature point in the second image comprises: determining a first projection transformation matrix between the first image and the second image according to a plurality of pairs of preprocessed control points between the first image and the second image; According to the first position information of each first feature point and the first projection transformation matrix, the second position information of the second feature point corresponding to the first feature point in the second image is obtained.

4. The multimodal remote sensing image registration method based on consistent gradient features according to claim 3 is characterized in that: Before determining a first projection transformation matrix between the first image and the second image according to a plurality of pairs of preprocessed control points between the first image and the second image, the method further includes: Acquire a first preprocessing control point in the first image and a second preprocessing control point in the second image; wherein each first preprocessing control point corresponds one-to-one to one second preprocessing control point; The first preprocessing control point and the second preprocessing control point having a corresponding relationship are taken as a pair of preprocessing control points to obtain the plurality of pairs of preprocessing control points.

5. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 1, characterized in that: The determining the minimum position error by template matching and fast Fourier transform according to the first feature descriptor, the second feature descriptor, the target feature point among the plurality of first feature points and the second image comprises: Using a window template of a preset size centered on a target first feature point, each search window is determined from the region of interest of the second image; the target feature point is any one of the multiple first feature points; According to the first feature descriptor and the second feature descriptor, the minimum value of the average sum of square errors between the window template and the search window is calculated by fast Fourier transform; each average sum of square errors is obtained based on a position error; each average sum of square errors corresponds to a search window; The position error corresponding to the minimum value is used as the minimum position error.

6. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 1, characterized in that: The step of removing each pair of mismatched points from the plurality of pairs of matching points to obtain a registration result includes: Determining a second projection transformation matrix between the first image and the second image according to the multiple pairs of matching points; Determine, according to the second projection transformation matrix and the first feature point in each pair of matching points, an estimated feature point corresponding to the first feature point; When the difference between the second position information of the third feature point corresponding to the first feature point and the position information of the estimated feature point is greater than or equal to a preset error threshold, the first feature point and the third feature point are eliminated to obtain multiple pairs of matching points after error detection; Each pair of mismatched points in the plurality of pairs of matching points after the error detection is eliminated to obtain the registration result.

7. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 1, characterized in that: The first gradient map includes: a first gradient magnitude map and a first gradient amplitude angle map; the second gradient map includes: a second gradient magnitude map and a second gradient amplitude angle map; the first gradient map of the first image and the second gradient map of the second image are respectively extracted, including: Using a first edge detection operator to respectively extract a first gradient magnitude map and a first gradient argument angle map of the first image; A second edge detection operator is used to respectively extract a second gradient magnitude map and a second gradient angle map of the second image; wherein the first edge detection operator is the same as or different from the second edge detection operator.

8. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 7, characterized in that: The first edge detection operator is a Sobel operator corresponding to four directional templates of 0°, 45°, 90° and 135°; the first gradient magnitude map and the first gradient angle map of the first image are respectively extracted by using the first edge detection operator, including: The four directional templates of 0°, 45°, 90° and 135° are used to respectively calculate the first-order derivative of each first pixel point of the first image, and correspondingly obtain a first derivative value, a second derivative value, a third derivative value and a fourth derivative value; Calculate the gradient amplitude of each first pixel point according to the absolute values ​​of the first derivative value, the second derivative value, the third derivative value, and the fourth derivative value to obtain the first gradient amplitude map; The gradient amplitude of each first pixel is calculated according to the first derivative value and the third derivative value to obtain the first gradient amplitude map.

9. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 7, characterized in that: The second edge detection operator is a filter composed of a causal filter and a non-causal filter; The step of respectively extracting a second gradient magnitude map and a second gradient angle map of the second image by using a second edge detection operator comprises: Calculate, according to the filter, the causal filter and the non-causal filter, a causal exponential weighted mean and a non-causal exponential weighted mean in the horizontal direction, and a causal exponential weighted mean and a non-causal exponential weighted mean in the vertical direction for each second pixel of the second image; Calculate a gradient estimate value in the horizontal direction of each second pixel point according to the causal exponential weighted mean and the non-causal exponential weighted mean in the horizontal direction; Calculate a gradient estimate in the vertical direction for each second pixel point according to the causal exponential weighted mean and the non-causal exponential weighted mean in the vertical direction; According to the gradient estimation value in the horizontal direction, the gradient estimation value in the vertical direction and the preset gradient difference coefficient, the gradient amplitude and gradient amplitude angle of each second pixel are calculated respectively to obtain the second gradient amplitude map and the second gradient amplitude angle map.

10. The multi-modal remote sensing image registration method based on consistent gradient features according to claim 1, characterized in that: The determining a first feature descriptor according to the first gradient consistency feature histogram includes: The first gradient consistency feature histogram is normalized to obtain the first feature descriptor.

Images