A Registration Method for Optical and SAR Remote Sensing Images with Large View Angle Differences

Abstract

The present invention discloses a registration method for optical and SAR remote sensing images with large view angle differences. The method uses the OS-SIFT algorithm, the multi-scale Harris-Affine algorithm, and the multi-scale MSER algorithm to extract key points between the reference image and the image to be registered respectively, establish descriptors and match them. Next, the global transformation model is calculated using all the matching point pairs obtained in the previous stage. After transforming the image to be registered, the RISFM algorithm is used to match on the original reference image and the transformed image to be registered. Finally, all the matching point pairs obtained in the previous two stages are aggregated, and after fine registration, the mismatched point pairs are removed, and a new global transformation model is calculated. The method of the present invention deeply considers the influence of speckle noise existing in SAR images and the influence brought by geometric distortion and non-linear radiation distortion existing between heterologous remote sensing images, and can significantly improve the registration performance of optical and SAR remote sensing images under large view angle differences.

Description

Technical Field

[0001] The present invention belongs to the technical field of multi-source remote sensing image registration, and particularly relates to a registration method for optical and SAR remote sensing images with large viewing angle differences. Background Art

[0002] With the development of space technology and multi-source observations, the amount of image data from multiple platforms, multiple bands, and multiple spectra is continuously increasing. Remote imaging sensors include optical sensors, synthetic aperture radar (SAR), and infrared radar, each with its own advantages. As an active microwave remote sensing system, SAR is not affected by weather and time and can provide high-resolution images with multiple viewing angles, multiple bands, and multiple polarizations. Optical sensors are passive sensors that can obtain rich gray-scale and texture information under good ground conditions. Due to the complementary advantages of different sensors, multi-source remote sensing image registration has become a hot topic for many scholars.

[0003] Most of the algorithms that have been proposed, such as the OS-SIFT (Optical-to-SAR Scale-Invariant Feature Transform, OS-SIFT) algorithm, although considering the different types of noise effects in optical and SAR images and using different methods to calculate image gradients, still cannot extract enough corresponding points. Moreover, when establishing descriptors, they still cannot handle the influence of non-linear radiation distortion well, resulting in an unsatisfactory final registration result. The RIFT (Radiation-Variation Insensitive Feature Transform, RIFT) algorithm, although using phase consistency instead of image intensity for feature point detection and describing features on the maximum index map, the established descriptors have good robustness to non-linear radiation distortion, but it ignores the influence of speckle noise in SAR images, resulting in many key points detected in SAR images being distributed in smooth areas and a low repetition rate of key points. The ASS (Adjacent Self-Similarity, ASS) algorithm, although designing a weight function to suppress speckle noise with reference to the Lee filter and obtaining good feature point detection results, the established descriptors are not robust enough, and the number of correct matching point pairs obtained finally is small. In addition, most of the algorithms that have been proposed do not consider the registration situation under large viewing angle differences. When there are serious geometric distortions between the reference image and the image to be registered, the registration accuracy is often poor.

[0004] Patent CN202310420988.6 ("A Registration Method for SAR Images with Large View Angle Differences") provides a registration method for SAR images with large view angle differences. Although the three-stage registration framework proposed by this method can obtain a large number of correct matching point pairs and achieve high registration accuracy on SAR images, it performs poorly on heterologous optical and SAR remote sensing image datasets. Although this method can solve the geometric distortion problem caused by large view angle differences in the registration of optical and SAR images, it does not consider that the noise types in heterologous images are different and different methods are needed to suppress the noise, and there are also serious non-linear radiation distortion effects between heterologous images. Summary of the Invention

[0005] To solve the problem of poor registration accuracy when there is serious geometric distortion between the reference image and the image to be registered in related technologies, the present invention provides a registration method for optical and SAR remote sensing images with large view angle differences. The technical problems to be solved by the present invention are realized through the following technical solutions:

[0006] The present invention provides a registration method for optical and SAR remote sensing images with large view angle differences, as Figure 1 , the method includes:

[0007] (1) Use the OS-SIFT algorithm, multi-scale Harris-Affine algorithm, and multi-scale MSER (Maximally Stable Extremal Regions, MSER) algorithm to extract key points in the reference image I1 and the image to be registered I2 respectively. Then, establish OS-SIFT descriptors and perform matching. Finally, three sets of matching point pairs with different properties are obtained: the set of point pairs S O that is not affected by the view angle change, the set of corner point pairs S H with affine invariance, and the set of region point pairs S M with affine invariance;

[0008] (2) According to all the matching point pairs obtained in step (1), calculate the global transformation model γ, and use the global transformation model γ to transform the image I2 to obtain the image Perform matching on the reference image I1 and the image to be registered using the RISFM (Radiation-Insensitive Structural Feature matching, RISFM) algorithm to obtain the set of matching point pairs S R ; among them, the specific steps of the RISFM algorithm are: first, establish the reference image I1 and the image to be registered The smallest self-similar graph, detecting key points on the smallest self-similar graph through maximum value detection and non-maximum suppression; then, based on the two-dimensional log-Gabor wavelet transform, calculating the log-Gabor responses of the reference image I1 and the image to be registered in each direction, finding the direction with the maximum response, and taking out the index value of the direction to establish the maximum index map MIM; then, establishing descriptors based on the detected key points and the maximum index map MIM; finally, using the B-NNDR (Bidirectional nearest neighbor distance ratio, B-NNDR) algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered and then, after removing the mismatched point pairs through the FSC (fast sample consensus, FSC) algorithm, obtaining the set S of matching point pairs R ;

[0009] (3) Aggregating the point pair sets S O 、S H 、S M and S R , performing fine registration using the LOS-Flow (Local Optical-to-SAR Flow, LOS-Flow) algorithm, and then using the LSOR (length-and slope-based outlier removal method, LSOR) algorithm to remove the possible mismatched point pairs therein, finally obtaining the set S of matching point pairs final , and recalculating the global transformation model using S final .

[0010] In some embodiments, step (1) includes:

[0011] 1a) The set S of point pairs not affected by the perspective change O : Extracting key points in the image using the OS-SIFT algorithm and establishing OS-SIFT descriptors, using the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered I2, and then, after removing the mismatched point pairs through the FSC algorithm, finally obtaining the set S of matching point pairs O ;

[0012] 1b) The set S of corner point pairs with affine invariance H: Use the multi-scale Harris-Affine algorithm to extract key points in the image and obtain the local affine transformation matrix at each key point. Among them, the steps of extracting key points in the image are the same as those of the OS-SIFT algorithm; then, at each key point, use the local affine transformation matrix to transform the image and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs H ;

[0013] 1c) Set S of regional point pairs with affine invariance M : Use the multi-scale MSER algorithm to extract key points in the image and obtain the local affine transformation matrix at each key point. At each key point, use the local affine transformation matrix to transform the image and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs M 。

[0014] In some embodiments, the specific steps of using the OS-SIFT algorithm to extract key points in steps 1a) and 1b) include:

[0015] 2a) Select a set of exponentially weighted scale factors [α0, α1,..., α n-1 , where the initial value α0 = 2, α i = α0 * k i , (i ∈ [1, n - 1]), k = 2 1 / 3 , n = 8;

[0016] 2b) According to the scale sequence [α0, α1,..., α n-1 of α, use the Sobel operator and the ROEWA (Ratio of Exponentially Weighted Averages, ROEWA) operator to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the reference image I1 and the image I2 to be registered at different scale factors respectively, so as to obtain the gradient magnitude Mag α and the gradient direction Ori α of the image at different scales:

[0017]

[0018]

[0019] 2c) Calculate the Harris matrix C at each pixel point in the image based on the calculated gradient H (x, y, α):

[0020]

[0021] where, represents a Gaussian kernel with a standard deviation of , and * represents convolution;

[0022] 2d) Use C H (x, y, α) to calculate the Harris response value R H (x, y, α):

[0023] R H (x, y, α) = det(C H (x, y, α)) - d · tr(C H (x, y, α)) 2 ;

[0024] where, d is a value between 0.04 and 0.06, det represents the value of the determinant, and tr represents the trace of the determinant;

[0025] 2e) In the same layer of the scale space image, compare the Harris response value of each pixel point with the response values of the pixel points in its 8-neighborhood around it, and a global threshold d H respectively. If the response value of the central point in the neighborhood is the largest and greater than the global threshold d H , then this point is the detected key point. Among them, the global threshold d H takes 0.85.

[0026] In some embodiments, the steps of establishing the OS-SIFT descriptor in steps 1a), 1b) and 1c) specifically include:

[0027] 3a) Using the gradient magnitude Mag α and the gradient direction Ori α , with the key point as the center, take a circular neighborhood with a radius of 6α. In the circular neighborhood, divide 0 - 360° into 18 equal parts, each part representing a range of 20°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. Traverse all the pixel points in the neighborhood. For each pixel point in the neighborhood, first find the histogram amplitude column corresponding to the gradient direction angle of this pixel point, and then accumulate the gradient magnitude of this pixel point on the histogram amplitude column of the direction angle where this point is located, so as to obtain the main direction histogram of this key point, and smooth the histogram; α represents the scale of the image layer where this key point is located;

[0028] 3b) Take the direction angle corresponding to the peak value of the smoothed histogram as the main direction of the key point, and the direction angles corresponding to the columns with peak energy greater than 80% as the secondary directions;

[0029] 3c) With the key point as the center, circular neighborhoods with radii of r max = 8α, 12α, and 16α are taken respectively. Each circular neighborhood is rotated with the main direction angle as the reference to make the feature descriptor rotation invariant. Then, each circular neighborhood is divided from the inside to the outside into a circle and two rings. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max respectively. Then, the two rings are further divided into 8 sub-regions at intervals of 45° within the range of 0 - 360°. Together with the central circle, each neighborhood is divided into a total of 17 sub-regions. The three circular neighborhoods have a total of 51 sub-regions. In each sub-region, 0 - 360° is evenly divided into 8 parts to establish a gradient direction histogram. The histogram vectors of the 51 sub-regions of the three circular neighborhoods are concatenated and normalized to generate a 408-dimensional OS-SIFT descriptor.

[0030] In some embodiments, the steps of using the B-NNDR algorithm to find the sampling set and the consistent set between the reference image and the image to be registered in steps 1a), 1b), 1c), and (2) specifically include:

[0031] 4a) In the key point space, the Euclidean distance between feature vectors is used to measure the proximity of two key points. The closer they are, the more similar they are. For a key point on the reference image, by calculating the Euclidean distance between the descriptors of the key points, the key point on the image to be registered that is the closest and the second closest to this key point are found. The closest distance and the second closest distance are represented by d min and d nd respectively; if then this key point and the key point closest to it are a pair of correct matching points; traverse each key point on the reference image. When the threshold distRatio is taken as 0.9, the sampling set is obtained. When the threshold distRatio is taken as 0.999, the consistent set

[0032] 4b) Traverse each key point on the image to be registered, and find the corresponding point on the reference image that is a pair of correct matching points with this key point. When the threshold distRatio is taken as 0.9, the sampling set is obtained. When the threshold distRatio is taken as 0.999, the consistent set

[0033] 4c) The sampling set and the consistent set between the reference image and the image to be registered are Ch and C l :

[0034]

[0035] ∪ represents taking the union set.

[0036] In some embodiments, the steps of using the FSC algorithm to remove mismatched point pairs in steps 1a), 1b), 1c) and (2) specifically include:

[0037] 5a) Set the number of iterations to N. During the t-th iteration, randomly select three pairs of matching points from the point pair set C h : t is an integer from 1 to N;

[0038] 5b) Use these three pairs of points to calculate the transformation model θ of the image t ;

[0039] 5c) Use the obtained transformation model θ t to calculate the transformation error e(c l ) of the matching point pair c in the point pair set C i ; i , θ t );

[0040]

[0041]

[0042] where, (x i , y i ) represents the key point coordinates of the matching point pair c i on the reference image, represents the corresponding point coordinates of the matching point pair c i on the image to be registered; T((x i , y i ), θ t ) represents the corresponding position on the image to be registered after transforming (x t , y i , y i ) using the transformation model θ i , θ t ); e(c t , θ i ) represents the transformation error of the matching point pair c under the transformation model θ

[0043] 5d) Traverse each pair of matching points in the point pair set C l , and summarize all matching point pairs with a transformation error e < 3 to obtain the point pair set C t ;

[0044] 5e) When the algorithm ends after N iterations, at this time, C1, C2,..., C N are obtained, a total of N point pair sets. Take out the point pair set with the largest number of point pairs among them as the matching point pair set finally obtained by the FSC algorithm.

[0045] In some embodiments, in step (2), establishing the minimum self-similarity graph of the reference image I1 and the image to be registered , and the steps of detecting key points on the minimum self-similarity graph through maximum value detection and non-maximum suppression specifically include:

[0046] 9a) Crop the reference image I1 to construct sub-images: Taking the central pixel point of the image I1 as the center, establish a search box with a size of L sub ×W sub , where L sub = L - 10, W sub = W - 10, and L and W are the length and width of the reference image I1; When the search box does not move, crop the reference image I1 according to the position of the search box to obtain the central sub-image block SubI c , and then after shifting the search box by one pixel in the directions of 0°, 45°, 90°, 135°, and 180° respectively, crop the reference image I1 according to the position of the current search box to obtain the offset sub-image blocks SubI1, SubI2, SubI3, SubI4, and SubI5;

[0047] 9b) Obtain the offset sub-image block SubI after shifting by one pixel in each direction through the following formula θ :

[0048]

[0049] θ = 180(o - 1) / N o , o = 1, 2,..., No; where θ represents the offset direction, and N o represents dividing θ ∈ [0°, 180°) into N o parts;

[0050] 9c) Calculate the weight value v(x, y) of each pixel point (x, y) on the central sub-image block SubI c ;

[0051]

[0052]

[0053]

[0054] σl (x, y) = η l (x, y) - (μ l (x, y)) 2 ;

[0055] Among them, μ l , η l and σ l respectively represent the local mean, mean square value and variance value; (x, y) represents a pixel point on the central sub-image block SubI c ; w l (x, y, r l ) represents a circular neighborhood centered at (x, y) with a radius of r l . When l takes 1, r1 is equal to 2, and when l takes 2, r2 is equal to 4; n l represents the number of pixel points in the circular neighborhood w l ;

[0056] 9d) Calculate the self-similarity map S in the θ direction using a mean filter θ :

[0057]

[0058] Among them, SubI c represents the central sub-image block, SubI θ represents the offset sub-image block shifted one pixel in the θ direction, v represents the weight value of each pixel point on the central sub-image block SubI c , w1 represents a circular neighborhood with a radius of 2, represents performing mean filtering on each pixel point with a circular neighborhood of radius 2;

[0059] 9e) Find the minimum self-similarity map S min :

[0060]

[0061] 9f) Perform maximum value detection and non-maximum suppression on the minimum self-similarity map S min to find the key points of the image I1;

[0062] 9g) Perform the same operations as above on the image to be registered to find the key points of the image .

[0063] In some embodiments, in step (2), based on the two-dimensional log-Gabor wavelet transform, calculate the reference image I1 and the image to be registered For the log-Gabor responses in all directions, the steps of finding the direction with the maximum response and extracting the index value of the direction to establish the maximum index map MIM specifically include:

[0064] 10a) Select 4 different scales s = 1, 2, 3, 4 and 6 different directions o = 0°, 30°, 60°, 90°, 120°, 150°, corresponding to index numbers 1, 2, 3, 4, 5, and 6 respectively, and establish a two-dimensional even-symmetric log-Gabor filter L even (x, y, s, o) and a two-dimensional odd-symmetric log-Gabor filter L odd (x, y, s, o);

[0065] 10b) Calculate the response amplitude A so (x, y) at the pixel point (x, y) of the image in the o direction and s scale:

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

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

[0068] where * represents convolution;

[0069]

[0070] 10c) Accumulate the response amplitudes at all different scales in the o direction to obtain the response value of the pixel point (x, y) in the o direction:

[0071]

[0072] 10d) Extract the index value corresponding to the maximum response as the value at the point (x, y), traverse each pixel point on the image, and establish the maximum index map MIM.

[0073] In some embodiments, the steps of establishing a descriptor based on the detected key points and the maximum index map MIM in step (2) specifically include:

[0074] 11a) For each key point on the maximum index map MIM, centered on it, crop an image patch of size 96×96, and then divide the image patch into 6×6 sub-grids. In each grid, establish a histogram where the abscissa represents the index number of the direction, with a value range of 1 to 6, and the ordinate represents the number of occurrences of the index number of this direction. In this way, each sub-grid can be transformed into a six-dimensional histogram vector. Finally, concatenate and normalize the six-dimensional histogram vectors of 36 sub-regions to generate a 216-dimensional descriptor for this key point.

[0075] In some embodiments, in step (3), the set of point pairs S O 、S H 、S M and S R After summarization, the steps of performing fine registration using the LOS-FLOW algorithm specifically include:

[0076] 12a) Using the set of point pairs S O 、S H 、S M Calculate the transformation model γ from the image to be registered I2 to the reference image I1, and use the transformation model γ to transform the image to be registered I2 to obtain the transformed image

[0077] For a point (x O 、S H 、S M on the image to be registered I2 in the set of point pairs, use the transformation model γ to calculate the corresponding position coordinates of the point (x s , y s ) on the transformed image s , y s )

[0078] Perform the above same transformation on all points in the set of point pairs S O 、S H 、S M to obtain a new set of point pairs

[0079] 12b) For image I1 and image respectively use the Sobel operator and the ROEWA operator to calculate the horizontal gradient G x and the vertical gradient G y at the scale factor α = 2, so as to obtain the gradient magnitude Mag and the gradient direction Ori of the image:

[0080] ​

[0081]

[0082] For a pixel point on the image I1 and taking a circular neighborhood with a radius of r centered at this pixel point max = 24, and dividing this circular neighborhood into a circle and two circular rings from the inside out. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max respectively from the inside out. Then, the two circular rings are divided into 8 sub-regions at intervals of 45° within the range of 0 to 360°. Together with the central circle, this neighborhood is divided into a total of 17 sub-regions. In each sub-region, 0 to 360° is evenly divided into 8 parts, each part representing a range of 45°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient amplitude. Traverse all the pixel points in this sub-region. For each pixel point in the region, first find the histogram amplitude column corresponding to the gradient direction angle of this pixel point. Then, accumulate the gradient amplitude of this pixel point on the histogram column in the direction angle where this point is located, thereby obtaining the gradient direction histogram of this point and converting it into an eight-dimensional histogram vector. After the histogram vectors of the 17 sub-regions are concatenated and normalized, the 136-dimensional OS-SIFT descriptor of this pixel point is generated;

[0083] By calculating the 136-dimensional OS-SIFT descriptors of each pixel point on the images I1 and respectively, the OS-SIFT descriptor map I1_desc of the image I1 and the OS-SIFT descriptor map of the image are formed

[0084] 12c) Set of candidate corresponding points is and the union of S R . For a point (x , y i ) on the image I1 in the set of candidate corresponding points i ), taking a local square region descriptor image I1_squ with a side length of 2·r i +1 centered at the point (x i , y lf ) on the OS-SIFT descriptor map I1_desc. Perform the same operation on the corresponding point i , y i ) of the point (x ) on the image to obtain the local square region descriptor image as where r lf takes 61;

[0085] 12d) Substitute the obtained I1_squ and into the loss function E(w) to calculate the optical flow vector w. The loss function E(w) is as follows:

[0086]

[0087] where p = (x, y) represents a certain pixel point on the local square region descriptor image, w(p) = (u(p), v(p)) represents the optical flow vector at point p, where u(p) represents the offset of point p in the horizontal direction, and v(p) represents the offset of point p in the vertical direction; ε represents the region centered at point p plus its adjacent 8 points, q represents a certain point in this region that does not include the center point p; the parameters η and α are taken as 0.001 and 0.03 respectively, and the parameters t and d are taken as 0.1 and 0.6 respectively; the above formula (1) is the data term, which constrains that the OS-SIFT descriptors along the optical flow vector w(p) should match each other; the formula (2) is the small displacement term, which constrains that in the absence of other available information, the optical flow vector should be as small as possible; the formula (3) is the smoothing term, which constrains that the optical flow vectors of adjacent pixels should be similar;

[0088] 12e) Calculate the new coordinates (x i T , y i T ) of the point (x) on the image as (x i T _ _new, y i T _new):

[0089] x i T _new = x i T + u(r lf + 1, r lf + 1);

[0090] y i T _new = y i T + v(r lf + 1, r lf + 1);

[0091] 12f) After traversing each pair of matching points in the set of points to be precisely matched , a more accurate set of point pairs

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

[0093] 1) The present invention deeply analyzes the different noise effects contained in optical and SAR remote sensing images. In the first stage, different methods are adopted to calculate the gradients of the optical map and the SAR map, so that the multi-scale Harris-Affine algorithm and the multi-scale MSER algorithm obtain more matching point pairs with affine invariance. By calculating the global transformation matrix to transform the image to be registered, the influence brought by the geometric distortion between images is preliminarily eliminated.

[0094] 2) Aiming at the influence brought by the non-linear radiation distortion in the registration of heterogeneous remote sensing images, the present invention proposes the RISFM algorithm. The algorithm uses the minimum self-similarity graph to detect key points, which greatly suppresses the influence of speckle noise in SAR images, and obtains a large number of key points with high repetition rate and stable properties in optical images and SAR images; the algorithm is based on the two-dimensional log-Gabor wavelet transform to obtain the responses of the image in multiple directions. By taking out the direction index value of the maximum response, a maximum index map is established. Descriptors are established and matched on the maximum index map, and by using the structural information, the influence brought by the non-linear radiation distortion between images is avoided, and finally a large number of correct matching point pairs are obtained.

[0095] 3) Aiming at the registration problem of optical and SAR remote sensing images with large view angle differences, the present invention proposes a three-stage registration framework, which effectively improves the accuracy of image registration.

[0096] The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. Description of the Drawings

[0097] Figure 1 It is a flowchart of the registration method for optical and SAR remote sensing images with large view angle differences provided by the embodiment of the present invention;

[0098] Figure 2 It is a schematic diagram of the stage for registering a reference image and an image to be registered provided by the embodiment of the present invention;

[0099] Figure 3 It is an image of the OS-SIFT descriptor provided by the embodiment of the present invention;

[0100] Figure 4 It is a flowchart of the implementation of the RISFM algorithm when the reference image is an optical image and the image to be registered is a SAR image provided by the embodiment of the present invention;

[0101] Figure 5 It is the reference image input for Test1 provided by the embodiment of the present invention;

[0102] Figure 6 It is the image to be registered input for Test1 provided by the embodiment of the present invention;

[0103] Figure 7 Reference image for the exemplary Test2 input provided by the embodiments of the present invention;

[0104] Figure 8 Image to be registered for the exemplary Test2 input provided by the embodiments of the present invention;

[0105] Figure 9 Registered checkerboard image for the exemplary Test1 provided by the embodiments of the present invention;

[0106] Figure 10 Correct matching point pair diagram for the exemplary Test1 provided by the embodiments of the present invention;

[0107] Figure 11 Registered checkerboard diagram for the exemplary Test2 provided by the embodiments of the present invention;

[0108] Figure 12 Correct matching point pair diagram for the exemplary Test2 provided by the embodiments of the present invention. Detailed implementation manners

[0109] The present invention will be further described in detail below in conjunction with specific embodiments, but the implementation manners of the present invention are not limited thereto.

[0110] In the description of the present invention, the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the present invention, "a plurality of" means two or more unless otherwise specifically defined.

[0111] In the description of this specification, the descriptions with reference to terms such as "an embodiment", "some embodiments", "example", "specific example", or "some examples" etc. mean that the specific features, structures, materials or characteristics described in connection 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 can be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art can combine and combine the different embodiments or examples described in this specification.

[0112] Although the present invention has been described in connection with the embodiments, however, in the process of implementing the claimed invention, those skilled in the art can understand and achieve other variations of the disclosed embodiments by viewing the accompanying drawings, the disclosure, and the appended claims. In the claims, the word "comprising" does not exclude other components or steps, and "a" or "an" does not exclude a plurality. A single processor or other unit can implement several functions recited in the claims. Certain measures are recited in mutually different dependent claims, but this does not mean that these measures cannot be combined to produce good results.

[0113] Figure 1 is a flowchart of a registration method for large-viewpoint-difference optical and SAR remote sensing images provided by an embodiment of the present invention; Figure 2 is a schematic diagram of the stage for registering a reference image and an image to be registered provided by an embodiment of the present invention; as Figure 1 and Figure 2 shown, when the present invention registers the reference image I1 and the corresponding image to be registered I2, it is divided into three stages:

[0114] (1) Use the OS-SIFT algorithm, the multi-scale Harris-Affine algorithm, and the multi-scale MSER algorithm to extract the key points in the reference image I1 and the image to be registered I2 respectively. Then, establish the OS-SIFT descriptor and perform matching. Finally, three sets of matching point pairs with different properties are obtained: the set of point pairs S O that is not affected by the perspective change, the set of corner point pairs S H with affine invariance, and the set of region point pairs S M with affine invariance;

[0115] 1a) The set of point pairs S O that is not affected by the perspective change: Use the OS-SIFT algorithm to extract the key points in the image and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered I2. After removing the mismatched point pairs through the FSC algorithm, the set of matching point pairs S O is finally obtained;

[0116] 1b) The set of corner point pairs S H: Extract key points in the image using the multi-scale Harris-Affine algorithm and obtain the local affine transformation matrix at each key point. Among them, the steps of extracting key points in the image are the same as those of the OS-SIFT algorithm. Then, at each key point, use the local affine transformation matrix to transform the image and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs H ;

[0117] 1c) The set S of regional point pairs with affine invariance M : Extract key points in the image using the multi-scale MSER algorithm and obtain the local affine transformation matrix at each key point. At each key point, use the local affine transformation matrix to transform the image and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs M 。

[0118] (2) According to all the matching point pairs obtained in step (1), calculate the global transformation model γ, and use the global transformation model γ to transform the image I2 to obtain the image Perform matching on the reference image I1 and the image to be registered using the RISFM algorithm to obtain the set S of matching point pairs R ; Among them, the specific steps of the RISFM algorithm are as follows: First, establish the minimum self-similarity graph of the reference image I1 and the image to be registered , and detect key points on the minimum self-similarity graph through maximum value detection and non-maximum suppression. Then, based on the two-dimensional log-Gabor wavelet transform, calculate the log-Gabor responses of the reference image I1 and the image to be registered in each direction, find the direction with the maximum response, and extract the index value of the direction to establish the maximum index map MIM. Then, based on the detected key points and the maximum index map MIM, establish the descriptor. Finally, use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered , and after removing the mismatched point pairs through the FSC algorithm, obtain the set S of matching point pairs R ;

[0119] (3) The set of point pairs S O 、S H 、S M and S RAfter summarization, the LOS-Flow algorithm is used for fine matching, and then the LSOR algorithm is used to remove the possibly included mismatched point pairs, and finally the set S of matching point pairs is obtained. final , and using S final to recalculate the global transformation model.

[0120] In the present invention, the operation steps of extracting key points in the image by using the OS-SIFT algorithm in steps 1a) and 1b) are as follows:

[0121] 2a) Select a set of exponentially weighted scale factors [α0, α1,..., α n-1 , where the initial value α0 = 2, α i = α0 * k i , (i ∈ [1, n - 1]), k = 2 1 / 3 , n = 8;

[0122] 2b) According to the scale sequence [α0, α1,..., α n-1 of α, use the Sobel operator and the ROEWA operator to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the reference image I1 and the image I2 to be registered at different scale factors, so as to obtain the gradient magnitude Mag α and the gradient direction Ori α of the image at different scales:

[0123]

[0124]

[0125] 2c) Calculate the Harris matrix C H (x, y, α) at each pixel point in the image according to the calculated gradient:

[0126]

[0127] Among them, represents a Gaussian kernel with a standard deviation of , and * represents convolution;

[0128] 2d) Use C H (x, y, α) to calculate the Harris response value R H (x, y, α) at each pixel point in the image:

[0129] R H (x, y, α) = det(C H (x, y, α)) - d · tr(C H (x, y, α)) 2 ;

[0130] Among them, d is a parameter of any size, generally between 0.04 and 0.06, det represents the value of the determinant, and tr represents the trace of the determinant.

[0131] 2e) In the same layer of the scale-space image, compare the Harris response value of each pixel point with the response values of the pixel points in its 8-neighborhood and a global threshold d H respectively. If the response value of the center point in the neighborhood is the largest and greater than the global threshold d H , then this point is the detected key point. Among them, the global threshold d H takes 0.85.

[0132] In the present invention, the operation steps of establishing the OS-SIFT descriptor at the found key points in steps 1a), 1b) and 1c) are as follows:

[0133] 3a) Using the gradient magnitude Mag α and the gradient direction Ori α , with the key point as the center, take a circular neighborhood with a radius of 6α. In the circular neighborhood, divide 0 to 360° into 18 equal parts, each part representing a range of 20°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. Traverse all the pixel points in this neighborhood. For each pixel point in the neighborhood, first find the histogram magnitude column corresponding to the gradient direction angle of this pixel point, and then accumulate the gradient magnitude of this pixel point on the histogram magnitude column of the direction angle where this point is located, so as to obtain the main direction histogram of this key point, and smooth the histogram;

[0134] 3b) Take the direction angle corresponding to the peak value of the smoothed histogram as the main direction of this key point, and the direction angles corresponding to the columns with peak energy greater than 80% as the secondary directions.

[0135] 3c) With the key point as the center, take circular neighborhoods with radii of r max = 8α, 12α, and 16α respectively. Rotate each circular neighborhood with the main direction angle as the reference to make the feature descriptor rotation-invariant. Then, divide each circular neighborhood into a circle and two rings from the inside to the outside. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max, then the two rings are respectively divided into 8 sub-regions at intervals of 45° within the range of 0 to 360°. Together with the central circle, each neighborhood is divided into a total of 17 sub-regions. The three circular neighborhoods have a total of 51 sub-regions. In each sub-region, 0 to 360° is evenly divided into 8 parts to establish a gradient direction histogram. The histogram vectors of the 51 sub-regions of the three circular neighborhoods are concatenated and normalized to generate a 408-dimensional OS-SIFT descriptor. Exemplarily, Figure 3 Is the image for generating the descriptor of a single circular domain.

[0136] In the present invention, the operation steps of using the B-NNDR algorithm to find the sampling set and the consistent set between the reference image and the image to be registered in steps 1a), 1b), 1c) and (2) are as follows:

[0137] 4a) In the key point space, the Euclidean distance between feature vectors is used to measure the distance between two key points. The closer they are, the more similar they are. For a key point on the reference image, by calculating the Euclidean distance between the descriptors of the key points, the key point closest and the second closest to this key point are found on the image to be registered. The closest distance and the second closest distance are represented by d min And d nd respectively; if then this key point and the key point closest to it are a pair of correct matching points; traverse each key point on the reference image. When the threshold distRatio is taken as 0.9, the sampling set is obtained. When the threshold distRatio is taken as 0.999, the consistent set

[0138] 4b) Traverse each key point on the image to be registered, and find the corresponding point that is a pair of correct matching points with this key point on the reference image. When the threshold distRatio is taken as 0.9, the sampling set is obtained. When the threshold distRatio is taken as 0.999, the consistent set

[0139] 4c) The sampling set and the consistent set between the reference image and the image to be registered are C h and C l :

[0140]

[0141] ∪ represents taking the union.

[0142] In the present invention, the operation steps of using the FSC algorithm to remove incorrect matching point pairs in steps 1a), 1b), 1c) and (2) are as follows:

[0143] 5a) Set the number of iterations to N. During the t-th iteration, randomly select three pairs of matching points from the set of point pairs C h where t is an integer from 1 to N. t is an integer from 1 to N.

[0144] 5b) Use these three pairs of points to calculate the transformation model θ of the image t .

[0145] 5c) Use the obtained transformation model θ t to calculate the transformation error e(c l ) of the matching point pair c in the set of point pairs C i ,θ i ,θ t ).

[0146]

[0147]

[0148] where (x i ,y i ) represents the key point coordinates of the matching point pair c on the reference image, i represents the corresponding point coordinates of the matching point pair c on the image to be registered; T((x ,y i ),θ i ,y i ),θ t ) represents the corresponding position on the image to be registered after transforming (x t ,y i ,y i ) using the transformation model θ i ,θ t ; e(c t ,θ i ) represents the transformation error of the matching point pair c under the transformation model θ

[0149] 5d) Traverse each pair of matching points in the set of point pairs C l , and summarize all the matching point pairs with a transformation error e < 3 to obtain a set of point pairs.

[0150] 5e) When the algorithm has iterated N times, it ends. At this time, C1, C2,..., C N are obtained, a total of N sets of point pairs. Take out the set of point pairs with the largest number of point pairs as the set of matching point pairs finally obtained by the FSC algorithm.

[0151] In the present invention, the operation steps of using the multi-scale Harris-Affine algorithm to calculate the affine transformation matrix at each key point in 1b) are as follows:

[0152] 6a) For a key point (x i , y i ) on the reference image I1, set the number of iterations K = 15 and initialize the shape adaptive matrix U (1) as the identity matrix E.

[0153] 6b) During the k-th iteration (here, k is an integer from 1 to K), use the shape adaptive matrix U (k) to transform the reference image I1 and the key point (x i , y i ) thereon, to obtain the transformed image and the key point coordinates

[0154]

[0155]

[0156] where T is a mapping operation that uses the shape adaptive matrix to transform the image or coordinates.

[0157] 6c) Taking as the center, take a square region W with a side length of 4α on the image , where α represents the scale of the image layer where the key point is located.

[0158] 6d) Use the Sobel operator and the ROEWA operator to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the square region W at the scale factor α for the reference image I1 and the image I2 to be registered respectively, so as to obtain the gradient magnitude Mag α of W:

[0159]

[0160] Calculate the Harris matrix C H (x, y, α) at each pixel point in W according to the already obtained horizontal and vertical gradients:

[0161]

[0162] where, represents a Gaussian kernel with a standard deviation of ;

[0163] Use C H (x, y, α) to calculate the Harris response value R H (x, y, α) at each pixel point in W:

[0164] R H(x, y, α) = det(C H (x, y, α)) - d·tr(C H (x, y, α)) 2 ; where d is a parameter of any size, generally between 0.04 and 0.06.

[0165] Take the point with the maximum Harris response value within the square region W as the new key point coordinates

[0166] 6e) Update the coordinates of the key point on the reference image I1:

[0167]

[0168] 6f) Centered at in the image reselect a square region W_new with a side length of 4α, and calculate the Harris matrix of the central pixel point within the new region W_new

[0169] 6g) Update the shape adaptive matrix:

[0170]

[0171] U (k+1) =(μ (k) ) -1 U (k) ;

[0172] After that, normalize the updated shape adaptive matrix U (k+1) , so that its maximum eigenvalue is equal to 1.

[0173] 6h) Calculate the convergence rate at the kth iteration:

[0174]

[0175] where λ min (μ (k) ) and λ max (μ (k) ) represent the minimum eigenvalue and the maximum eigenvalue of the matrix μ (k) respectively.

[0176] 6i) Exit the loop when ratio < 0.1, and obtain the affine transformation matrix Matrix i , y i ) at the key point (x i = U (k) , otherwise, increment the iteration count by 1 and continue the loop, and recalculate at the key point (x using the updated shape adaptive matrix U (k+1) ​i , y i ) the convergence rate ratio at; if the condition ratio < 0.1 cannot be satisfied after K iterations, then discard this key point.

[0177] 6j) Perform the same above operations on each key point on the reference image I1 and the image I2 to be registered.

[0178] In the present invention, the operation steps of transforming the image (for example, taking the reference image I1 as an example) using the affine transformation matrix at each key point in steps 1b) and 1c) are as follows:

[0179] 7a) For a certain key point (x i , y i ) on the reference image I1, use the affine transformation matrix Matrix i to transform the gradient magnitude map Mag α and the gradient direction map Ori α of the reference image I1, as well as the key point coordinates (x i , y i ):

[0180] Mag α T = T(Mag α , Matrix i );

[0181] Ori α T = T(Ori α , Matrix i );

[0182] (x i T , y i T ) = T((x i , y i ), Matrix i ).

[0183] In the present invention, the operation steps of extracting key points in the image using the multi-scale MSER algorithm and calculating the affine transformation matrix at each key point in step 1c) are as follows:

[0184] 8a) Select a set of scale space factors [α0, α1,..., α n-1 , where the initial value α0 = 2, α i = α0 * k i , (i ∈ [1, n - 1]), k = 2 1 / 3, where n = 4, which is the scale of the Gaussian kernel, select the window length w = 4α, and establish a scale space using the Gaussian blur kernel;

[0185] 8b) For each layer of the image, if the image of this layer is a color image, it is necessary to first convert the color image into a grayscale image. After that, for each layer of the image, sort according to the grayscale value; allocate a node for each pixel point in this layer of the image in advance, and the node index number is the grayscale value corresponding to the pixel point. According to the sorting result of the pixel points, place them into the component tree one by one, and the placement order is the node index number corresponding to each pixel point; during the placement process, first place the pixel point, and then check the four-neighbor positions of the pixel point. If there are nodes, find their respective root nodes and merge the two node areas. After all pixel points are placed into the component tree, all the extreme value regions corresponding to this layer of the image are obtained; among them, the extreme value region is defined as: if the grayscale values of all pixels in a certain region are greater than the grayscale values of its boundary pixels, then this region is defined as the maximum extreme value region; if the grayscale values of all pixels in the region are less than the grayscale values of its boundary pixels, then it is defined as the minimum extreme value region.

[0186] 8c) Use the maximum stability determination condition to obtain the MSER region: If Q1,...,Q i-1 ,Q i ,... are a series of mutually inclusive extreme value regions, that is If the extreme value region Q i* is the maximum stable extreme value region, if and only if the region change rate q(i) = |Q i+Δ -Q i-Δ | / |Q i | obtains a local minimum at i * , where, - represents taking the non-common part of the two regions, |·| represents the number of pixel points in the region, the subscript i ∈ [0, 255] represents the gray level, and Δ represents a small gray level change;

[0187] 8d) Approximately fit the irregular maximum stable value region into an elliptical region: First, take the centroid of the maximum stable value region as the center of the ellipse, and calculate the center of the ellipse:

[0188] Calculate the geometric zero-order moment and geometric first-order moment of the maximum stable value region:

[0189] m 00 = ∑I e (x,y);

[0190] m 01 = ∑yI e (x,y);

[0191] m 10 = ∑xI e(x, y);

[0192] where m 00 , m 01 and m 10 are respectively the geometric zero - order moment and the geometric first - order moment of the maximally stable extremal region, and I e (x, y) represents the maximally stable extremal region, so that the center coordinates of the ellipse can be obtained, that is, the key - point coordinates (x c , y c ) detected by the MSER algorithm:

[0193]

[0194]

[0195] Calculate the geometric second - order moment of the maximally stable extremal region:

[0196] where, μ 20 = ∑(x - x c ) 2 I e (x, y), μ 02 = ∑(y - y c ) 2 I e (x, y), μ 11 = ∑(x - x c )(y - y c )I e (x, y);

[0197] Calculate the two eigenvalues of the geometric second - order moment:

[0198]

[0199]

[0200] Calculate the major semi - axis w, minor semi - axis l and the major - axis direction of the ellipse :

[0201]

[0202]

[0203]

[0204] 8e) Use the major semi - axis w, minor semi - axis l and the major - axis direction of the ellipse to calculate the affine transformation matrix at the key - point (x c , y c ):

[0205]

[0206] 8f) Invert the grayscale of the original image and repeat the above operations;

[0207] 8g) Perform the same above operations on each layer of the image in the scale space.

[0208] In the present invention, in step (2), the minimum self-similarity graph of the reference image I1 and the image to be registered is obtained. The operation steps for detecting key points on the minimum self-similarity graph through maximum value detection and non-maximum suppression are as follows:

[0209] 9a) Crop the reference image I1 to construct sub-images: Taking the central pixel point of the image I1 as the center, establish a search box with a size of L sub ×W sub , where L sub = L - 10, W sub = W - 10, and L and W are the length and width of the reference image I1; when the search box does not move, crop the reference image I1 according to the position of the search box to obtain the central sub-image block SubI c . Then, after shifting the search box by one pixel in the directions of 0°, 45°, 90°, 135°, and 180° respectively, crop the reference image I1 according to the position of the current search box to obtain the offset sub-image blocks SubI1, SubI2, SubI3, SubI4, and SubI5;

[0210] 9b) Obtain the offset sub-image block SubI θ shifted by one pixel in each direction through the following formula:

[0211]

[0212] θ = 180(o - 1) / N o , o = 1, 2,..., N o

[0213] where θ represents the direction of offset, and N o represents dividing θ ∈ [0°, 180°) into N o parts, and N o can be set according to actual needs.

[0214] 9c) Calculate the weight value v(x, y) of each pixel point (x, y) on the central sub-image block SubI c ;

[0215]

[0216]

[0217]

[0218] σ l (x,y) = η l (x,y) - (μ l (x,y)) 2 ;

[0219] where μ l , η l and σ l represent the local mean, mean square value, and variance value respectively; (x,y) represents a pixel point on the central sub-image block SubI c ; w l (x,y,r l ) represents a circular neighborhood centered at (x,y) with a radius of r l ; when l takes 1, r1 is equal to 2, and when l takes 2, r2 is equal to 4; n l represents the number of pixel points in the circular neighborhood w l ;

[0220] 9d) Calculate the self-similarity graph S in the θ direction using a mean filter θ :

[0221]

[0222] where SubI c represents the central sub-image block, SubI θ represents the offset sub-image block offset by one pixel in the θ direction, v represents the weight value of each pixel point on the central sub-image block SubI c ; w1 represents a circular neighborhood with a radius of 2, indicating mean filtering for each pixel point with a circular neighborhood of radius 2;

[0223] 9e) Find the minimum self-similarity graph S min :

[0224]

[0225] 9f) Perform maximum value detection and non-maximum suppression on the minimum self-similarity graph S min to find the key points of the image I1;

[0226] 9g) Perform the same operations as above on the image to be registered to find the key points of the image ;

[0227] In the present invention, in step (2), based on the two-dimensional log-Gabor wavelet transform, calculate the reference image I1 and the image to be registered For the log-Gabor responses in all directions, find the direction with the maximum response, and extract the index value of the direction to establish the maximum index map MIM. The operation steps are as follows:

[0228] 10a) Select four different scales s = 1, 2, 3, 4 and six different directions o = 0°, 30°, 60°, 90°, 120°, 150°, corresponding to the index numbers 1, 2, 3, 4, 5, and 6 respectively, and establish a two-dimensional even-symmetric log-Gabor filter L even (x, y, s, o) and a two-dimensional odd-symmetric log-Gabor filter L odd (x, y, s, o).

[0229] 10b) Calculate the response amplitude A so (x, y) at the pixel point (x, y) of the image in the o direction and s scale:

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

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

[0232] where * represents convolution;

[0233]

[0234] 10c) Accumulate the response amplitudes at all different scales in the o direction to obtain the response value of the pixel point (x, y) in the o direction:

[0235]

[0236] 10d) Extract the index value corresponding to the maximum response as the value at the point (x, y), and traverse each pixel point on the image to establish the maximum index map MIM.

[0237] In the present invention, the operation steps for establishing a descriptor based on the detected key points and the maximum index map MIM in step (2) are as follows:

[0238] 11a) For each key point on the maximum index map MIM, a 96×96 image patch is cropped with it as the center, and then the image patch is divided into 6×6 sub-grids. In each grid, a histogram is established. The abscissa represents the index number of the direction, with a value range of 1 to 6, and the ordinate represents the number of times the index number of this direction appears. In this way, each sub-grid can be transformed into a six-dimensional histogram vector. Finally, the six-dimensional histogram vectors of 36 sub-regions are concatenated and normalized to generate a 216-dimensional descriptor for this key point.

[0239] Exemplarily, as Figure 4 shown, for the optical image and the SAR image after rough correction, on the one hand, key points are detected by establishing the minimum self-similarity map and using maximum value detection and non-maximum suppression, and on the other hand, the index map is established by calculating the log-Garbor responses in multiple directions and extracting the direction index value of the maximum response; then, descriptors are established on the index map based on the detected key points; finally, the optical image and the SAR image after rough correction are matched according to the descriptors of the key points to obtain a large number of correct matching point pairs.

[0240] In the present invention, in step (3), after aggregating the point pair sets S O , S H , S M and S R , the operation steps of fine matching using the LOS-FLOW algorithm are as follows:

[0241] 12a) Using the transformation model γ of the to-be-registered image I2 to the reference image I1 calculated from the point pair sets S O , S H , S M , the to-be-registered image I2 is transformed using the transformation model γ to obtain the transformed image

[0242] For a point (x O , y H ) on the to-be-registered image I2 in the point pair sets S M , the corresponding position coordinates of the point (x s , y s ) on the transformed image s , y s ) are calculated using the transformation model γ on the transformed image

[0243] All points in the point pair sets S O , S H , S M are subjected to the above same transformation, so as to obtain a new point pair set

[0244] 12b) For image I1 and image Calculate the horizontal gradient G and vertical gradient G at the scale factor α = 2 using the Sobel operator and ROEWA operator respectively x and y , so as to obtain the gradient magnitude Mag and gradient direction Ori of the image:

[0245]

[0246]

[0247] For image I1 and For a pixel point on, taking this pixel point as the center, take a circular neighborhood with a radius of r max = 24, and divide this circular neighborhood into a circle and two rings from the inside to the outside. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max , respectively. Then divide the two rings into 8 sub-regions at intervals of 45° within the range of 0 - 360°. Together with the central circle, this neighborhood is divided into a total of 17 sub-regions. In each sub-region, divide 0 - 360° into 8 equal parts, each part representing a range of 45°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. Traverse all the pixel points in this sub-region. For each pixel point in the region, first find the histogram magnitude column corresponding to the gradient direction angle of this pixel point. Then, accumulate the gradient magnitude of this pixel point on the histogram column in the direction angle where this point is located, so as to obtain the gradient direction histogram of this point, and convert it into an eight-dimensional histogram vector. After splicing and normalizing the histogram vectors of the 17 sub-regions, the 136-dimensional OS-SIFT descriptor of this pixel point is generated;

[0248] By calculating the 136-dimensional OS-SIFT descriptors of each pixel point on image I1 and respectively, the OS-SIFT descriptor map I1_desc of image I1 and the OS-SIFT descriptor map of image are formed

[0249] 12c) The set of candidate matching point pairs is and the union of S R . For a point (x in the set of candidate matching point pairs i , y i ) on image I1, taking the point (x i , yi ) centered, take a local square region descriptor image I1_squ with a side length of 2·r lf + 1 from the OS-SIFT descriptor graph I1_desc, for the point (x i , y i ) on the image perform the same operation on the corresponding point (x i T , y i T ), and the resulting local square region descriptor image is where r lf takes 61;

[0250] 12d) Substitute the obtained I1_squ and into the loss function E(w), and calculate the optical flow vector w. The loss function E(w) is:

[0251]

[0252] where p = (x, y) represents a certain pixel point on the local square region descriptor image, w(p) = (u(p), v(p)) represents the optical flow vector of point p, where u(p) represents the offset of point p in the horizontal direction, and v(p) represents the offset of point p in the vertical direction; ε represents the region centered on point p plus its adjacent 8 points, q represents a certain point in this region that does not include the center point p; the parameters η and α are taken as 0.001 and 0.03 respectively, and the parameters t and d are taken as 0.1 and 0.6 respectively; the above formula (1) is the data term, which constrains that the OS-SIFT descriptors along the optical flow vector w(p) should match each other; the formula (2) is the small displacement term, which constrains that in the absence of other available information, the optical flow vector should be as small as possible; the formula (3) is the smoothing term, which constrains that the optical flow vectors of adjacent pixels should be similar;

[0253] 12e) Calculate the new coordinates (x on the image for the point (x i T , y i T ) (x i T _new, y i T _new):

[0254] x i T _new = x i T + u(r lf + 1, r lf + 1);

[0255] y i T _new = y i T + v(r lf + 1, r lf + 1).

[0256] 12f) After traversing each pair of matching points in the set of pairs of points to be precisely matched, a more accurate set of pairs of points is obtained in the set.

[0257] In the present invention, the operation steps of using the LSOR algorithm to remove the possibly included mismatched pairs of points in step (3) are as follows:

[0258] 13a) Calculate the length and slope of the line connecting each pair of matching points in the set of pairs of points and average the sum thereof to obtain the average length dist_ave and the average slope slope_ave.

[0259] 13b) For each pair of matching points in the set of pairs of points calculate the length dist i and the slope slope i , and screen by setting the thresholds Th d and Th s to retain the pairs of matching points that meet the screening conditions, obtaining The screening conditions are as follows:

[0260] |dist i - dist_ave| < Th d ;

[0261] |slope i - slope_ave| < Th s ; where Th d takes 0.1 and Th s takes 5°.

[0262] 13c) For each point in the set of pairs of points in the image use the transformation model γ -1 to calculate the corresponding point of this point on the original image I2, obtaining the set of pairs of points S final .

[0263] The technical effects of the embodiments of the present invention are further described below through simulation experiment data.

[0264] Two sets of publicly available datasets are selected (the first set of datasets is from the paper "A deep translation (GAN) based change detection network for optical and SAR remote sensing images", Xinghua Li et al., 2021; the second set of datasets is from the paper "Self-Supervised Keypoint Detection and Cross-Fusion Matching Networks for Multimodal Remote Sensing Image Registration", Liangzhi Li et al., 2022). To further verify the performance of the registration method for optical and SAR remote sensing images with large viewing angle differences proposed by the present invention, one image in each of the two sets of publicly available datasets is transformed respectively to simulate the registration environment under large viewing angle differences.

[0265] Table 1 below shows the registration results of the method of the present invention on the two datasets, and compares them with the existing registration methods, the Affine-SIFT algorithm (abbreviated as ASIFT, from the paper "ASIFT: A New Framework for Fully Affine Invariant Image Comparison", SIAM Journal on Imaging Sciences, J.M. Moreal et al., 2009) and the OS-SIFT algorithm (from the paper "OS-SIFT: A Robust SIFT-Like Algorithm for High-Resolution Optical-to-SAR Image Registration in Suburban Areas", IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, Yuming Xiang et al., 2018). Precision in Table 1 represents the precision rate, that is, the percentage of correct point pairs in the finally obtained set of point pairs; RMSE represents the root mean square error; for Test1, the input reference image and the image to be registered are Figure 5 and Figure 6 , and for Test2, the input reference image and the image to be registered are Figure 7 and Figure 8 . Figure 9 and Figure 11 are the checkerboard images after registration for Test1 and Test2 respectively; Figure 10 and Figure 12They are respectively the correct matching point pair graph of Test1 and the correct matching point pair graph of Test2.

[0266] Table 1 Comparison of the registration performance between the method of the present invention and the existing methods

[0267]

[0268] It can be seen from Table 1 and the above result graphs that the method proposed by the present invention has better registration performance.

[0269] The above content is a further detailed description of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is only limited to these descriptions. For those of ordinary skill in the technical field to which the present invention pertains, without departing from the concept of the present invention, several simple deductions or substitutions can be made, and all should be regarded as belonging to the protection scope of the present invention.

Claims

1. A registration method for optical and SAR remote sensing images with large viewing angle differences, characterized in that, Including: (1) Extract the key points in the reference image I1 and the image I2 to be registered using the OS-SIFT algorithm, the multi-scale Harris-Affine algorithm, and the multi-scale MSER algorithm respectively. Then, establish the OS-SIFT descriptor and perform matching. Finally, three sets of matching point pairs with different properties are obtained: the set of point pairs S that is not affected by the perspective change O , the set of corner point pairs S with affine invariance H , and the set of region point pairs S with affine invariance M ; (2) Calculate the global transformation model γ based on all the matching point pairs obtained in step (1), and use the global transformation model γ to transform the image I2 to obtain an image Perform matching on the reference image I1 and the image to be registered using the RISFM algorithm to obtain a set of matching point pairs S R ; Among them, the specific steps of the RISFM algorithm are as follows: First, establish the minimum self-similarity graph of the reference image I1 and the image to be registered , and detect key points on the minimum self-similarity graph through maximum value detection and non-maximum suppression; then, based on the two-dimensional log-Gabor wavelet transform, calculate the log-Gabor responses of the reference image I1 and the image to be registered in each direction, find the direction of the maximum response, and extract the index value of the direction to establish the maximum index map MIM; then, based on the detected key points and the maximum index map MIM, establish descriptors; finally, use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered , and then remove the mismatched point pairs through the FSC algorithm to obtain the set of matching point pairs S R ; (3) Aggregate the point pair sets S O , S H , S M and S R , then perform fine matching using the LOS-Flow algorithm, and use the LSOR algorithm to remove the potentially included mismatched point pairs. Finally, obtain the matching point pair set S final , and use S final to recalculate the global transformation model.

2. The registration method for large-viewpoint-difference optical and SAR remote sensing images according to claim 1, wherein Step (1) includes: 1a) Set S of point pairs unaffected by perspective change O : Use the OS-SIFT algorithm to extract key points in the image and establish OS-SIFT descriptors. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. Then, after removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs O ; 1b) Set S of corner point pairs with affine invariance H : Use the multi-scale Harris-Affine algorithm to extract key points in the image and obtain the local affine transformation matrix at each key point. Among them, the steps of extracting key points in the image are the same as those of the OS-SIFT algorithm; then, use the local affine transformation matrix to transform the image at each key point and establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image to be registered I2, and then remove the mismatched point pairs through the FSC algorithm to finally obtain the set S of matching point pairs H ; 1c) Set S of region point pairs with affine invariance M : Use the multi-scale MSER algorithm to extract key points in the image, and obtain the local affine transformation matrix at each key point. After transforming the image with the local affine transformation matrix at each key point, establish the OS-SIFT descriptor. Use the B-NNDR algorithm to find the sampling set and the consistent set between the reference image I1 and the image I2 to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set S of matching point pairs M .

3. The registration method for large-viewpoint-difference optical and SAR remote sensing images according to claim 2, characterized in that, The specific steps of using the OS-SIFT algorithm to extract key points in the image in steps 1a) and 1b) include: 2a) Select a set of exponentially weighted scale factors [α0, α1,..., α n-1 , where the initial value α0 = 2, α i = α0 * k i , (i ∈ [1, n - 1]), k = 2 1 / 3 , n = 8; 2b) According to the scale sequence of α [α0, α1,..., α n-1 , for the reference image I1 and the image to be registered I2, use the Sobel operator and the ROEWA operator respectively to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the image at different scale factors, so as to obtain the gradient magnitude Mag α and the gradient direction Ori α : 2c) Calculate the Harris matrix C at each pixel point in the image according to the already obtained gradient H (x, y, α): Among them, represents a Gaussian kernel with a standard deviation of , and * represents convolution; 2d) Use C H to calculate the Harris response value R at each pixel point in the image H (x, y, α): R H (x, y, α) = det(C H (x, y, α)) - d·tr(C H (x, y, α)) 2 ; Where d is a value between 0.04 and 0.06, det represents the value of the determinant, and tr represents the trace of the determinant; 2e) In the same layer of the scale space image, compare the Harris response value of each pixel point with the response values of the pixel points in its surrounding 8-neighborhood, as well as a global threshold d H respectively. If the response value of the center point in the neighborhood is the largest and greater than the global threshold d H , then this point is the detected key point, where the global threshold d H is taken as 0.

85.

4. The registration method for large-viewpoint-difference optical and SAR remote sensing images according to claim 3, wherein The specific steps of establishing the OS-SIFT descriptor in steps 1a), 1b) and 1c) specifically include: 3a) Using the gradient magnitude Mag α and the gradient direction Ori α , with the key point as the center, a circular neighborhood with a radius of 6α is taken. Within the circular neighborhood, the range from 0 to 360° is evenly divided into 18 parts, each part representing a range of 20°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. All pixel points within the neighborhood are traversed. For each pixel point within the neighborhood, first, the histogram magnitude bin corresponding to the gradient direction angle of this pixel point is found. Then, the gradient magnitude of this pixel point is accumulated on the histogram magnitude bin of the direction angle where this point is located, thereby obtaining the main direction histogram of this key point, and the histogram is smoothed; 3b) Use the direction angle corresponding to the peak of the smoothed histogram as the main direction of the key point, and the direction angles corresponding to the bins with peak energy greater than 80% as the secondary directions; 3c) Taking the key point as the center, circular neighborhoods with radii of r max = 8α, 12α, and 16α are respectively taken. Each circular neighborhood is rotated with the main direction angle as the reference to make the feature descriptor rotation-invariant. Then, each circular neighborhood is divided from the inside to the outside into a circle and two circular rings. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max , respectively. Then, the two circular rings are each divided into 8 sub-regions at intervals of 45° within the range of 0 to 360°. Together with the central circle, each neighborhood is divided into a total of 17 sub-regions. The three circular neighborhoods have a total of 51 sub-regions. In each sub-region, 0 to 360° is evenly divided into 8 parts to establish a gradient direction histogram. The histogram vectors of the 51 sub-regions of the three circular neighborhoods are concatenated and normalized to generate a 408-dimensional OS-SIFT descriptor.

5. The registration method for large-viewpoint-difference optical and SAR remote sensing images according to claim 1, characterized in that, The specific steps of using the B-NNDR algorithm to find the sampling set and the consistent set between the reference image and the image to be registered in steps 1a), 1b), 1c) and (2) specifically include: 4a) In the key point space, the Euclidean distance between feature vectors is used to measure the distance between two key points. The closer they are, the more similar they are. For a key point on the reference image, by calculating the Euclidean distance between the descriptors of the key points, the key point closest and the second closest to this key point are searched for on the image to be registered. The closest distance and the second closest distance are represented by d min and d nd respectively; if then this key point and the key point closest to it are a pair of correctly matched key points; traverse each key point on the reference image. When the threshold distRatio is taken as 0.9, the sampling set is obtained. When the threshold distRatio is taken as 0.999, the consensus set 4b) Traverse each key point on the image to be registered, and find the corresponding point on the reference image that is a correct matching point pair with this key point; when the threshold distRatio is taken as 0.9, the sampling set is obtained When the threshold distRatio is taken as 0.999, the consensus set is obtained 4c) The sampling set and the consensus set between the reference image and the image to be registered are C h and C l : The symbol ∪ represents taking the union.

6. The registration method for large-view-angle-difference optical and SAR remote sensing images according to claim 5, characterized in that, The specific steps of using the FSC algorithm to remove mismatched point pairs in steps 1a), 1b), 1c) and (2) specifically include: 5a) Set the number of iterations to N. During the t-th iteration, randomly select three pairs of matching points from the set of point pairs C h ; t is an integer from 1 to N; 5b) Calculate the transformation model θ of the image using these three pairs of points t ; 5c) Use the obtained transformation model θ t to calculate the matching point pairs c l in the point pair set C i and the transformation error e(c i , θ t ); Among them, (x i , y i ) represents the key point coordinates of the matching point pair c i on the reference image, and represents the corresponding point coordinates of the matching point pair c i on the image to be registered; T((x i , y i ), θ t ) represents the corresponding position on the image to be registered after transforming (x t , y i , y i ) using the transformation model θ i , θ t ) represents the transformation error of the matching point pair c t under the transformation model θ i ; 5d) Traverse the set of point pairs C l For each pair of matching points in it, summarize all the matching point pairs with a transformation error e < 3 to obtain the set of point pairs C t ; 5e) When the algorithm ends after iterating N times, at this time, C1, C2,..., C N , a total of N point pair sets are obtained. Take out the point pair set with the largest number of point pairs among them as the matching point pair set finally obtained by the FSC algorithm.

7. The registration method for large-view-angle-difference optical and SAR remote sensing images according to claim 1, characterized in that In step (2), establishing the minimum self-similarity graph of the reference image I1 and the image to be registered The steps of detecting key points on the minimum self-similarity graph through maximum value detection and non-maximum suppression specifically include: 9a) Crop the reference image I1 to construct sub-images: centered on the central pixel of the image I1, establish a search box of size L sub ×W sub , where L sub = L - 10, W sub = W - 10, and L and W are the length and width of the reference image I1; when the search box does not move, crop the reference image I1 according to the position of the search box to obtain the central sub-image block SubI c , and then offset the search box by one pixel at 0°, 45°, 90°, 135°, and 180° respectively, and crop the reference image I1 according to the position of the current search box to obtain the offset sub-image blocks SubI1, SubI2, SubI3, SubI4, and SubI5; 9b) The offset sub-image block SubI after being offset by one pixel in each direction is obtained through the following formula θ :[[]]END]] θ = 180(o - 1) / N o , o = 1, 2, ..., N o ; where θ represents the direction of the offset, and N o represents dividing θ ∈ [0°, 180°) into N o parts; 9c) Calculate the weight value v(x, y) of each pixel point (x, y) on the central sub-image block SubI c ; σ l (x,y) = η l (x,y) - (μ l (x,y)) 2 ; Among them, μ l , η l and σ l represent the local mean, mean square value, and variance value respectively; (x, y) represents a pixel point on the central sub-image block SubI c ; w l (x, y, r l ) represents a circular neighborhood centered at (x, y) with a radius of r l ; when l takes 1, r1 is equal to 2, and when l takes 2, r2 is equal to 4; n l represents the number of pixel points within the circular neighborhood w l ; 9d) Calculate the self-similarity graph S in the θ direction using a mean filter θ : Among them, SubI c represents the central sub-image block, SubI θ represents the offset sub-image block offset by one pixel in the θ direction, and v represents the weight value of each pixel point on the central sub-image block SubI c ; w1 represents a circular neighborhood with a radius of 2, indicating performing mean filtering on a circular neighborhood with a radius of 2 for each pixel point; 9e) Search for the minimum self-similar graph S min : 9f) Perform maximum value detection and non-maximum suppression on the minimum self-similarity graph S min to find the key points of the image I1; 9g) Image to be registered Perform the same operations as above to find the keypoints of the image.

8. The registration method for large-view-angle-difference optical and SAR remote sensing images according to claim 1, wherein In step (2), based on the two-dimensional log-Gabor wavelet transform, calculate the log-Gabor responses of the reference image I1 and the image to be registered The steps of finding the direction with the maximum response in the log-Gabor responses in each direction and extracting the index value of the direction to establish the maximum index map MIM specifically include: 10a) Select 4 different scales s = 1, 2, 3, 4 and 6 different directions o = 0°, 30°, 60°, 90°, 120°, 150°, corresponding to index numbers 1, 2, 3, 4, 5 and 6 respectively, and establish two-dimensional even-symmetric log-Gabor filter L even (x, y, s, o) and two-dimensional odd-symmetric log-Gabor filter L odd (x, y, s, o); 10b) Calculate the response amplitude A of the pixel point (x, y) on the image in the o direction and at the s scale so (x, y): E so (x,y) = I(x,y) * L even (x,y,s,o), O so (x,y) = I(x,y) * L odd (x,y,s,o), Where * represents convolution; 10c) Accumulate the response amplitudes at all different scales in the o direction to obtain the response value of the pixel point (x, y) in the o direction: 10d) Take out the direction index value corresponding to the maximum response as the value at the point (x, y), traverse each pixel point on the image, and establish the maximum index map MIM.

9. The registration method for large-view-angle-difference optical and SAR remote sensing images according to claim 1, wherein The specific steps of establishing the descriptor based on the detected key points and the maximum index map MIM in step (2) specifically include: 11a) For each key point on the maximum index map MIM, centered on it, crop an image patch with a size of 96×96, then divide the image patch into 6×6 sub-grids. In each grid, establish a histogram, where the abscissa represents the index number of the direction, and the value range is 1 to 6, and the ordinate represents the number of times the index number of this direction appears. In this way, each sub-grid can be transformed into a six-dimensional histogram vector. Finally, concatenate and normalize the six-dimensional histogram vectors of 36 sub-regions to generate a 216-dimensional descriptor for this key point.

10. The registration method for large-view-angle-difference optical and SAR remote sensing images according to claim 1, characterized in that, In step (3), after aggregating the point pair sets S O , S H , S M and S R , the specific steps of performing fine matching using the LOS-FLOW algorithm include: 12a) Using the set of point pairs S O , S H , S M The transformation model γ from the image I2 to be registered to the reference image I1 calculated, and the image to be registered I2 is transformed using the transformation model γ to obtain the transformed image For the point-to-point set S O , S H , S M For a point (x s , y s ) on the image I2 to be registered, use the transformation model γ to calculate the corresponding position coordinates of the point (x s , y s ) on the transformed image ​ Apply the same transformation to all the points in the set of point pairs S O , S H , S M to obtain a new set of point pairs 12b) For image I1 and the image respectively use the Sobel operator and the ROEWA operator to calculate the horizontal gradient G at the scale factor α = 2 x and the vertical gradient G y , so as to obtain the gradient magnitude Mag and the gradient direction Ori of the image: For an image I1 and a pixel on it, taking this pixel as the center, a circular neighborhood with a radius of r max = 24 is taken, and this circular neighborhood is divided into a circle and two circular rings from the inside to the outside. The radii of the concentric circles are 0.25·r max , 0.75·r max and r max respectively. Then, the two circular rings are each divided into 8 sub-regions at intervals of 45° within the range of 0 to 360°. Together with the central circle, this neighborhood is divided into a total of 17 sub-regions. In each sub-region, 0 to 360° is evenly divided into 8 parts, each part representing a range of 45°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient amplitude. All pixel points within the sub-region are traversed. For each pixel point within the region, first, the histogram amplitude column corresponding to the gradient direction angle of this pixel point is found. Then, the gradient amplitude of this pixel point is accumulated on the histogram column of the direction angle where this point is located, thereby obtaining the gradient direction histogram of this point and converting it into an eight-dimensional histogram vector. After the histogram vectors of the 17 sub-regions are spliced and normalized, a 136-dimensional OS-SIFT descriptor of this pixel point is generated; By separately calculating the 136-dimensional OS-SIFT descriptors for each pixel point in the image I1 and the OS-SIFT descriptor graph I1_desc of the image I1 is formed, and for the image the OS-SIFT descriptor graph 12c) Set of pairs of points to be refined is the union of R and S. For a point (x in the set of pairs of points to be refined i , y i ) on the image I1, taking the point (x i , y i ) as the center, a local square region descriptor image I1_squ with a side length of 2·r lf +1 is taken from the OS-SIFT descriptor graph I1_desc. The same operation is performed on the corresponding point (x i , y i ) of the point (x on the image i T , y i T ), and the resulting local square region descriptor image is where r lf is taken as 61; 12d) Substitute the obtained I1_squ and into the loss function E(w) to calculate the optical flow vector w, where the loss function E(w) is: Among them, p = (x, y) represents a certain pixel point on the local square region descriptor image, and w(p) = (u(p), v(p)) represents the optical flow vector of point p, where u(p) represents the offset of point p in the horizontal direction, and v(p) represents the offset of point p in the vertical direction; ε represents the region centered at point p and composed of its adjacent 8 points, q represents a certain point in this region that does not include the center point p; the parameters η and α are taken as 0.001 and 0.03 respectively, and the parameters t and d are taken as 0.1 and 0.6 respectively; the above formula (1) is a data item, which constrains the OS along the optical flow vector w(p). - The SIFT descriptors need to match each other; formula (2) is a small displacement term, which constrains that in the absence of other available information, the optical flow vector should be as small as possible; formula (3) is a smoothing term, which constrains that the optical flow vectors of adjacent pixels should be similar; 12e) Calculate the image of the point (x i T , y i T ) on the image to get the new coordinates (x i T _new, y i T _new): x i T _new = x i T + u(r lf + 1, r lf + 1); y i T _new = y i T + v(r lf + 1, r lf + 1); 12f) Traverse the set of pairs of points to be precisely matched After traversing each pair of matching points in the set, a more accurate set of point pairs is obtained

Images