A Registration Method for SAR Images with Large Viewing Angle Differences

Abstract

The present invention relates to a registration method for SAR images with large view angle differences. The method uses the SAR-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. After establishing descriptors and matching them, the LSS-Flow algorithm is used for fine registration to obtain three sets of matching point pairs with different properties. Then, by constructing Voronoi polygons, the A-SAR-SIFT algorithm with affine invariance is used for registration and fine registration to obtain the matching point set again. Finally, the matching point sets obtained in the first two stages are summarized and screened, and the global transformation model is calculated. The method of the present invention deeply considers the influence of speckle noise existing in SAR images and the severe geometric distortion and radiometric distortion existing between images, and can significantly improve the registration performance of SAR images under large view angle differences.

Description

Technical Field

[0001] The present invention belongs to the technical field of SAR image registration, and particularly relates to a registration method for SAR images with large view angle differences. Background Technique

[0002] Synthetic Aperture Radar (SAR), as an active microwave remote sensing imaging system, has the imaging characteristics of all-weather, all-time, and high resolution, and can obtain ground observation data of targets of interest under different polarization modes, different bands, and different view angle conditions. At present, the registration method for SAR images has become a research hotspot for many scholars.

[0003] Most of the algorithms that have been proposed, such as the SAR-SIFT algorithm, have achieved good registration accuracy and efficiency in general SAR image registration tasks, and also performed outstandingly in subsequent tasks such as data classification, 3D reconstruction, change detection, and image fusion. However, these methods only utilize single-property features in the reference image and the image to be registered, such as corner features, while ignoring other complementary features such as lines and planes. Therefore, when the feature information contained in the imaging scene is less or the image is blurred due to poor imaging quality, these methods will fail to register because they cannot find sufficient correct matching point pairs; in addition, when the number of mismatched point pairs is large, it is difficult for existing methods to screen out more correct matching point pairs and use the fine registration algorithm to converge them to a more accurate position, resulting in an inaccurate global transformation model; finally, due to the side-looking mechanism of the SAR sensor, there are often serious geometric distortions between the reference image and the image to be registered, and the existing methods have poor registration accuracy for SAR images with large view angle differences.

[0004] The paper "Registrating Oblique SAR Images Based on Complementary Integrated Filtering and Multilevel Matching" provides an image registration method based on the combination of the MSER algorithm and the Harris-Affine algorithm. This method does not deeply consider that there are still some matching point pairs between the reference image and the image to be registered that are not easily affected by perspective changes; when using the SAR-SIFT algorithm for matching, only the affine invariance of the matching point pairs obtained by the previous MSER algorithm and Harris-Affine algorithm is utilized, while ignoring the corresponding relationship and spatial position information between the matching point pairs. There is a lack of spatial position guidance during matching, increasing a lot of interference; only the fast sample consensus (FSC) algorithm is used to remove the mismatched point pairs, and the running time of the algorithm is relatively long and the efficiency is not high; the refined registration effect obtained by the proposed refined registration algorithm still has room for further improvement. Summary of the Invention

[0005] In view of the deficiencies of the existing SAR image registration methods, the present invention proposes a registration method for SAR images with large perspective differences, and uses a three-stage registration framework to register SAR images with large perspective differences. First, different algorithms are used to fully extract different types of matching point pairs between the reference image and the image to be registered and perform refined registration; then, the A-SAR-SIFT algorithm with affine invariance is used to establish descriptors and perform local matching by constructing Voronoi polygons. Next, all the obtained matching point pairs are processed through the mismatched point pair removal step and the refined registration step; finally, all the refined matching point pairs obtained in the first two stages are aggregated, and after removing the mismatched point pairs, the global transformation model is calculated. This SAR image registration method can better complete the registration task of SAR images with large perspective differences and further improve the registration performance of SAR images.

[0006] A registration method for SAR images with large perspective differences, comprising the following steps:

[0007] (1) Use the SAR-SIFT algorithm, the multi-scale Harris-Affine algorithm, and the multi-scale MSER algorithm to extract the key points in the reference image and the image to be registered respectively. Then, establish SAR-SIFT descriptors and perform matching, and finally obtain three sets of matching point pairs with different properties: the set of point pairs S that is not easily affected by perspective changes S,c , the set of corner point pairs S with affine invariance H,c and the set of point pairs S in the region with affine invariance M,c ;

[0008] (2) Using the LSS-FLOW algorithm, we can find the point pair set S S,c , S H,c and S M,c Perform fine matching to obtain a more accurate point pair set S S,p , S H,p and S M,p ;

[0009] (3) The A-SAR-SIFT algorithm is used to extract key points on the reference image and the image to be registered and to establish descriptors. Next, the nearest neighbor principle is used to find the initial matching point pairs in each corresponding small polygon in the reference image and the image to be registered. The LSOR algorithm is used to preliminarily remove the mismatched point pairs from the matching point pair set on the entire image. Then, the FSC algorithm is used for secondary screening to obtain the matching point pair set S A,c , and then use the LSS-Flow algorithm for precise matching, and finally get the matching point pair set S A,p ;

[0010] (4) Set the point pair set S S,p , S H,p , S M,p With S A,p After summarizing, the LSOR algorithm is used to remove possible mismatched point pairs, and finally the matching point pair set S is obtained. final , using S final Compute the global transformation model.

[0011] Furthermore, the point pair set S that is not easily affected by the change of viewing angle is obtained above S,c , a set of corner point pairs S with affine invariance H,c And the set of point pairs S with affine invariant regions M,c The specific steps are:

[0012] 1a) A set of point pairs S that are not easily affected by perspective changes S,c :

[0013] The SAR-SIFT algorithm is used to extract key points in the image and establish the SAR-SIFT descriptor. The NNDR algorithm is used to find the initial matching point pair set between the reference image and the image to be registered. The FSC algorithm is then used to remove the mismatched point pairs, and finally the matching point pair set S is obtained. S,c ;

[0014] 1b) A set of corner point pairs S with affine invariance H,c :

[0015] Extract key points in the image using the multi-scale Harris-Affine algorithm and obtain the affine transformation matrix at each key point. The steps of extracting key points in the image are the same as those of the SAR-SIFT algorithm. Next, use the local affine transformation matrix to transform the image at each key point, then establish the SAR-SIFT descriptor. Use the NNDR algorithm to find the set of initial matching point pairs between the reference image and the image to be registered, and then remove the mismatched point pairs through the FSC algorithm. Finally, obtain the set of matching point pairs S H,c ;

[0016] 1c) The set of point pairs S with affine invariant regions M,c :

[0017] Extract key points in the image using the multi-scale MSER algorithm and obtain the affine transformation matrix at each key point. Use the local affine transformation matrix to transform the image at each key point, then establish the SAR-SIFT descriptor. Use the NNDR algorithm to find the set of initial matching point pairs between the reference image and the image to be registered, and then remove the mismatched point pairs through the FSC algorithm. Finally, obtain the set of matching point pairs S M,c 。

[0018] Furthermore, for the above registration method of SAR images with large view angle differences, in step (3), use the A-SAR-SIFT algorithm to extract key points on the reference image and the image to be registered and establish descriptors. The specific steps are as follows: First, use the SAR-SIFT algorithm to extract key points on the reference image and the image to be registered, and then use the set S H,p and S M,p in all point pairs as objects to establish Voronoi polygons respectively, obtain the local affine transformation matrix of each SAR-SIFT key point, and use the local affine transformation matrix to transform the image to establish the SAR-SIFT descriptor.

[0019] Furthermore, for the above registration method of SAR images with large view angle differences, in steps 1a), 1b) and (3), use the SAR-SIFT algorithm to extract key points in the image. The operation steps are as follows:

[0020] 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;

[0021] 2b) From the scale sequence of α, use the ROEWA algorithm to calculate the horizontal gradient G x,α and the vertical gradient G y,α, thus obtaining the gradient magnitude Mag of the image at different scales α and the gradient direction Ori α :

[0022]

[0023]

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

[0025]

[0026] where the parameter α is the scale parameter of the exponential weighting function, is the standard deviation of the Gaussian function;

[0027] 2d) Use C(x, y, α) to calculate the SAR - Harris response value R(x, y, α) at each pixel point in the image SH (x, y, α), SH (x, y, α),

[0028] R SH (x, y, α) = det(C SH (x, y, α)) - d·tr(C SH (x, y, α)) 2

[0029] where d is a parameter of any size, generally between 0.04 and 0.06;

[0030] 2e) In the same layer of the scale - space image, compare the SAR - Harris response value of each point with the response values of the points in its 8 - neighborhood and a global threshold d SH respectively. If the response value of the central point in the neighborhood is the largest and greater than the global threshold d SH , then this point is the detected key point, where the global threshold d SH is generally taken as 0.85.

[0031] Furthermore, for the above - mentioned registration method for SAR images with large view - angle differences, in steps 1a), 1b), 1c) and (3), establish SAR - SIFT descriptors at the found key points, and the operation steps are as follows:

[0032] 3a) Use the gradient magnitude Mag α and the gradient direction Ori α, taking a circular neighborhood with a radius of 6α centered at the key point. Within the circular neighborhood, divide 0 to 360° into 36 equal parts, each part representing a range of 10°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient amplitude. Traverse all pixel points within this neighborhood. For each pixel point within the neighborhood, first find the histogram amplitude column corresponding to the gradient direction angle of this pixel point, and then accumulate the gradient amplitude of this pixel point on the histogram amplitude column of the direction angle where this point is located, thereby obtaining the main direction histogram of this key point, and perform smoothing processing on the histogram;

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

[0034] 3c) Taking a circular neighborhood with a radius of r max = 12α centered at the key point, rotate the circular neighborhood with the main direction angle as the reference to make the feature descriptor rotation-invariant, and then establish a gradient direction histogram. Divide this circular neighborhood into 17 sub-regions, and the concentric circle radii from the inside to the outside are 0.25·r max , 0.75·r max and r max , and divide 0 to 360° into 8 equal parts, calculate the gradient direction histogram of each sub-region, and splice and normalize the histogram vectors of the 17 sub-regions to generate a 136-dimensional SAR-SIFT descriptor.

[0035] Furthermore, for the above registration method for SAR images with large viewing angle differences, in step 1b), the affine transformation matrix at the key point is calculated using the multi-scale Harris-Affine algorithm, and its operation steps are as follows:

[0036] 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;

[0037] 6b) In the k-th iteration process, use the shape-adaptive matrix U (k ) to transform the reference image I1 and the key point (x i , y i ) on it to obtain the transformed image and the key point coordinates:

[0038]

[0039]

[0040] where T is using the shape-adaptive matrix U (k ) A mapping operation for transforming images or coordinates;

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

[0042] 6d) Using the ROEWA algorithm to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the square region W at the scale factor α, so as to obtain the gradient magnitude Mag α :

[0043]

[0044] Calculating the SAR-Harris matrix C SH (x, y, α) at each pixel point in W according to the calculated gradient;

[0045]

[0046] Using C SH to calculate the SAR-Harris response value R SH (x, y, α)

[0047] R SH (x, y, α) = det(C SH (x, y, α)) - d · tr(C SH (x, y, α)) 2

[0048] where d is a parameter of any size, generally between 0.04 and 0.06;

[0049] Taking the point with the maximum SAR-Harris response value within the square region W as the new key point coordinates

[0050] 6e) Updating the coordinates of the key point on the original reference image I1:

[0051]

[0052] 6f) Taking as the center, a new square region W_new with a side length of 4α is reselected on the image , and calculating the SAR-Harris matrix of the central pixel point within the new region

[0053] 6g) Updating the shape adaptive matrix:

[0054]

[0055] U (k+1) =(μ (k) ) -1 U (k)

[0056] And normalize the shape adaptive matrix so that its maximum eigenvalue is equal to 1;

[0057] 6h) Calculate the convergence rate at the k-th iteration:

[0058]

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

[0060] 6i) When ratio < 0.1, exit the loop 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 the convergence rate ratio at the key point (x i ,y i ) using the updated shape adaptive matrix; if the condition ratio < 0.1 cannot be satisfied after K iterations, discard this key point;

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

[0062] Furthermore, for the registration method of SAR images with large viewing angle differences mentioned above, in step 1c), the multi-scale MSER algorithm is used to extract key points in the image and calculate the affine transformation matrix at each key point, and its operation steps are as follows:

[0063] 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 , n = 4, which is the scale of the Gaussian kernel. Select the window length w = 4α and establish the scale space using the Gaussian blur kernel;

[0064] 8b) Sort each layer of the image according to the gray value. If it is a color image, it needs to be converted into a gray image. Assign a node to each pixel point in the image in advance, and the node index number is the gray 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 the pixel points are placed into the component tree, all the extreme value areas corresponding to the image are obtained. Among them, the extreme value area is defined as follows: If the gray values of all pixels in a certain area are greater than the gray values of its boundary pixels, then this area is defined as the maximum extreme value area; if the gray values of all pixels in the area are less than the gray values of its boundary pixels, it is defined as the minimum extreme value area.

[0065] 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 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 the area of the region, the subscript i ∈ [0,255] represents the gray level, and Δ represents a small gray change amount.

[0066] 8d) Approximately fit the irregular maximum stable value region into an elliptical region.

[0067] First, take the centroid of the maximum stable value region as the center of the ellipse and calculate the center of the ellipse:

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

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

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

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

[0072] where m 00 、m 01 and m 10 are the geometric zero-order moment and geometric first-order moment of the maximum stable value region respectively, and Ie (x, y) represents the maximally stable extremal region, from which the center coordinates of the ellipse can be obtained, i.e., the key point coordinates (x c , y c ):

[0073]

[0074]

[0075] Calculate the geometric second moment of the maximally stable extremal region:

[0076]

[0077] where

[0078] μ 20 = ∑(x - x c ) 2 I e (x, y)

[0079] μ 02 = ∑(y - y c ) 2 I e (x, y)

[0080] μ 11 = ∑(x - x c )(y - y c )I e (x, y)

[0081] Calculate the two eigenvalues of the geometric second moment:

[0082]

[0083]

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

[0085]

[0086]

[0087]

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

[0089]

[0090] 8f) Repeat the above operations by taking the negative value of the grayscale of the input image;

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

[0092] Furthermore, for the above registration method for SAR images with large viewing angle differences, in steps (2) and (3), the LSS-FLOW algorithm is used to perform fine registration on the point pair set to obtain a more accurate point pair set, and its operation steps are as follows:

[0093] 9a) Calculate the transformation model θ from the to-be-registered image to the reference image using the matching point pairs in the reference point pair set . For steps (2) and (3), the reference point pair set is the sum of S S,c , S H,c and S M,c . The reference image is I1, and the to-be-registered image is I2. Deform the to-be-registered image I2:

[0094]

[0095] where represents the image after transforming the to-be-registered image using the transformation model θ;

[0096] 9b) For the points of the to-be-finely-registered point pair set S coarse on the to-be-registered image I2, calculate their corresponding positions on the transformed image I2 T . Among them, for step (2), the to-be-finely-registered point pair set is S S,c , S H,c and S M,c . For step (3), the to-be-finely-registered point pair set is S A,c :

[0097] (x i T , y i T ) = T((x i , y i ), θ);

[0098] 9c) Use the ROEWA algorithm to calculate the horizontal gradient G T and the vertical gradient G x of the images I1 and I2 y at the scale factor α = 2, so as to obtain the gradient magnitude Mag and the gradient direction Ori of the images:

[0099]

[0100]

[0101] For a pixel point on the image I1 and taking a circular neighborhood with a radius of r = 24 centered at this pixel point, and dividing this neighborhood into 17 sub-regions. The concentric circle radii are 0.25·r, 0.75·r, and r from the inside to the outside. Within 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 pixel points within this sub-region. For each pixel point within the region, 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 column in the direction angle where this point is located, thus obtaining the gradient direction histogram of this point, and converting it into a one-dimensional histogram vector. The histogram vectors of the 17 sub-regions are concatenated and normalized to generate a 136-dimensional SAR-SIFT descriptor;

[0102] Calculate the SAR-SIFT descriptors for each pixel point on the image I1 and and concatenate them, thus forming the SAR-SIFT descriptor map I1_desc of the image I1 and and

[0103] 9d) For a matching point (x c o arse , y i , i ) on the reference image I1 in the set S lf of the set of points to be accurately matched, taking it as the center, extracting a square region descriptor image I1_squ with a side length of 2·r i T + 1 on the SAR-SIFT descriptor map I1_desc, and performing the same operation on its corresponding point (x i T ) on the image to be registered, and the obtained square region descriptor image is where r lf takes 61;

[0104] 9e) The loss function E(w) is:

[0105]

[0106] Among them, \(p=(x,y)\) represents a certain pixel point on the image, and \(w(p)=(u(p),v(p))\) represents the optical flow vector of point \(p\) on the image, 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; \(\varepsilon\) represents the region centered at point \(p\) plus its adjacent 8 points, and \(q\) represents a certain point in this region; the parameters \(\eta\) and \(\alpha\) are taken as 0.001 and 0.01 respectively, and the parameters \(t\) and \(d\) are both taken as 1;

[0107] Among them, equation (1) is the data item, which constrains that the SAR - SIFT descriptors along the optical flow vector \(w(p)\) should match each other; equation (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; equation (3) is the smoothing term, which constrains that the optical flow vectors of adjacent pixels should be similar;

[0108] Use the belief propagation algorithm to solve for \(w\) when the loss function \(E(w)\) takes the minimum value;

[0109] 9f) Calculate the new coordinates of the point \((x\) i T ,y\) i T ) on the image:

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

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

[0112] 9g) Use the transformation model \(\theta\) of the image to calculate the coordinate position of the point \((x\) i T _new,y i T _new) on the original image to be registered \(I2\):

[0113] (x i _new,y i _new) = T((x i T _new,y i T _new),\theta -1 )

[0114] 9h) Traverse the set S of pairs of points to be refined coarse After traversing each pair of matching points in it, a more accurate set S of point pairs is obtained precise , where, for step (2), the more accurate set of point pairs is S S,p , S H,p and S M,p , and for step (3), the more accurate set of point pairs is S A,p .

[0115] Furthermore, for the above registration method for SAR images with large view angle differences, in step (3), taking S H,p and S M,p as objects to establish Voronoi polygons, and obtaining the local affine transformation matrix of each SAR-SIFT key point, the operation steps are as follows:

[0116] 10a) Utilize all the discrete points with affine invariance on the reference image I1 in the sets S H,p and S M,p to automatically construct a Delaunay triangulation network with these discrete points as objects; number these discrete points and the formed triangles, and record which three discrete points each triangle is composed of;

[0117] 10b) Find out and record the numbers of all the triangles adjacent to each discrete point;

[0118] 10c) Sort and record the triangles adjacent to each discrete point in the clockwise direction for the next step of connecting to generate Voronoi polygons;

[0119] 10d) Calculate the circumcenter of each triangle and record it;

[0120] 10e) According to the adjacent triangles of each discrete point, connect the circumcenters of these adjacent triangles, and then Voronoi polygons can be obtained;

[0121] 10f) Assign the same affine transformation matrix as the discrete points with affine invariance in the polygon to the SAR-SIFT key points located in each small Voronoi polygon, then the local affine transformation matrix of each SAR-SIFT key point is obtained;

[0122] 10g) Utilize all the discrete point pairs with affine invariance on the image I2 to be registered in the sets S H,p and S M,p to perform the above same operations to obtain the local affine transformation matrix of each SAR-SIFT key point.

[0123] Further, for the above registration method for SAR images with large viewing angle differences, in steps (3) and (4), the LSOR algorithm is used to preliminarily remove mismatched point pairs from the matching point pairs on the entire image, and the operation steps are as follows:

[0124] 11a) Using the reference point pair set Calculate the transformation model θ' from the image to be registered to the reference image, where, for step (3), the reference point pair set is the sum of S S,p , S H,p and S M,p ; for step (4), the reference point pair set is the sum of S S,p , S H,p , S M,p and S A,p .

[0125] The reference image is I1, and the image to be registered is I2. The image to be registered I2 is deformed using the transformation model θ' to obtain the image

[0126]

[0127] 11b) For a certain point (x , y s , y s ) on the image to be registered I2 in the reference point pair set , calculate its corresponding position coordinates (x s T , y s T ) on the transformed image

[0128] (x s T , y s T ) = T((x s , y s ), θ')

[0129] Perform the above same transformation on all point pairs in the reference point pair set to obtain a new point pair set S ref T ;

[0130] 11c) Calculate the length and slope of the line connecting each pair of matching points in the point pair set S ref T , and after summing them, calculate the average length dist_ave and the average slope slope_ave;

[0131] 11d) For the set of point pairs to be screened S filAt a point (x i ,y i ), and use the transformation model θ' to calculate its transformed image The corresponding position coordinates (x i T ,y i T ). For step (3), the set of point pairs to be screened is the set of initial matching point pairs on the entire image found by the nearest neighbor principle in the previous step; for step (4), the set of point pairs to be screened is S S,p , S H,p , S M,p and S A,p The sum of:

[0132] (x i T ,y i T )=T((x i ,y i ),θ')

[0133] 11e) Calculate the point pair (x i ,y i ) and (x i T ,y i T )The length of the connection dist i and slope i ;

[0134] 11f) Set the threshold to filter the point pair set S fil Filter all point pairs in and keep the matching point pairs that meet the conditions:

[0135] |dist i -dist_ave|<Th d

[0136] |slope i -slope_ave|<Th s

[0137] Among them Th d Take 0.1, Th s Take 5 degrees.

[0138] Beneficial effects of the present invention:

[0139] 1. The present invention analyzes key points of different natures between a reference image and an image to be registered, and proposes a method of combining a multi-scale MSER algorithm, a multi-scale Harris-Affine algorithm, and a SAR-SIFT algorithm to jointly extract key points in the image, establish descriptors for matching, and effectively overcome the difficulty of registration failure due to insufficient correct matching point pairs when there is less feature information in the imaging scene or the image is blurred due to poor imaging quality;

[0140] 2. The present invention proposes to establish Voronoi polygons around the obtained multi-scale MSER matching point pairs and multi-scale Harris-Affine matching point pairs with affine invariance on the reference image and the image to be registered, endow the SAR-SIFT key points falling within each small Voronoi polygon with the same affine transformation property as the matching points with affine invariance within the polygon, and then perform registration by establishing descriptors. The proposed A-SAR-SIFT algorithm is different from the traditional SAR-SIFT algorithm, enabling each key point to have affine invariance, and more correct matching point pairs can be found for the registration problem of SAR images with large viewing angle differences;

[0141] 3. The present invention proposes an LSOR algorithm for the problem that it is difficult to screen out correct matching point pairs when there are many mismatching point pairs. By setting thresholds based on length and slope conditions, the mismatching point pairs in the set of point pairs are effectively removed; for the problem that the positions between correct matching point pairs are not accurate enough and the error is large, an LSS-Flow fine registration algorithm is proposed in combination with the idea of the optical flow algorithm. The application of the mismatching point pair removal method and the fine registration method makes the finally obtained global transformation model more accurate;

[0142] 4. The present invention proposes a three-stage registration framework for the registration problem of SAR images with large viewing angle differences, effectively improving the accuracy of image registration. Description of the Drawings

[0143] Figure 1 is a schematic diagram of the implementation process of this embodiment;

[0144] Figure 2 is an image for establishing SAR-SIFT descriptors;

[0145] Figure 3 is an image for calculating gradients using the ROEWA algorithm;

[0146] Figure 4 is the reference image input for Test1;

[0147] Figure 5 is the image to be registered input for Test1;

[0148] Figure 6 It is a reference diagram of the input of Test2;

[0149] Figure 7 It is the diagram to be registered of the input of Test2;

[0150] Figure 8 It is the checkerboard diagram after registration of Test1;

[0151] Figure 9 It is the correct matching point pair diagram of Test1;

[0152] Figure 10 It is the checkerboard diagram after registration of Test2;

[0153] Figure 11 It is the correct matching point pair diagram of Test2;

[0154] Figure 12 It is the implementation flow chart of this embodiment;

[0155] Figure 13 It is the implementation flow chart of the A-SAR-SIFT algorithm;

[0156] Figure 14 It is the input description diagram of the LSOR algorithm. Detailed implementation manners

[0157] The present invention will be further described in detail below with reference to the accompanying drawings.

[0158] Referring to Figure 1 , Figure 2 , Figure 3 and Figure 12 , this embodiment provides a registration method for SAR images with large viewing angle differences. The registration method includes the following steps:

[0159] (1) Use the SAR-SIFT (SAR scale-invariant feature transform), multi-scale Harris-Affine algorithm, and multi-scale Maximally Stable Extremal Regions (MSER) algorithm to extract key points between the reference diagram and the diagram to be registered respectively. Then, by establishing SAR-SIFT descriptors and matching, three sets of matching point pairs with different properties are finally obtained: the set of point pairs S that is not easily affected by viewing angle changes S,c , the set of corner point pairs S with affine invariance H,c and the set of point pairs S with regions of affine invariance M,c .

[0160] 1a) A set of point pairs S that are not easily affected by perspective changes S,c :

[0161] The SAR-SIFT algorithm is used to extract key points in the image and establish the SAR-SIFT descriptor. The nearest neighbor distance ratio (NNDR) algorithm is used to find the initial matching point pair set between the reference image and the image to be registered. After the FSC algorithm is used to remove the mismatched point pairs, the matching point pair set S is finally obtained. S,c .

[0162] 1b) A set of corner point pairs S with affine invariance H,c :

[0163] The multi-scale Harris-Affine algorithm is used to extract the key points in the image and obtain the affine transformation matrix at each key point. The steps of extracting the key points in the image are the same as those of the SAR-SIFT algorithm. Next, the image is transformed at each key point using the local affine transformation matrix to establish the SAR-SIFT descriptor. The NNDR algorithm is used to find the initial matching point pair set between the reference image and the image to be registered. The FSC algorithm is then used to remove the mismatched point pairs, and finally the matching point pair set S is obtained. H,c .

[0164] 1c) A set of point pairs S with affine invariant regions M,c :

[0165] The multi-scale MSER algorithm is used to extract the key points in the image and obtain the affine transformation matrix at each key point. The image is transformed using the local affine transformation matrix at each key point to establish the SAR-SIFT descriptor. The NNDR algorithm is used to find the initial matching point pair set between the reference image and the image to be registered. The FSC algorithm is then used to remove the mismatched point pairs, and finally the matching point pair set S is obtained. M,c .

[0166] (2) Using the Local-SAR-SIFT Flow (LSS-Flow) algorithm, we can obtain the point pair set S S,c , S H,c and S M,c Perform fine matching to obtain a more accurate point pair set S S,p , S H,p and S M,p .

[0167] (3) The A-SAR-SIFT algorithm is used to extract the key points on the reference image and the image to be registered and to establish the descriptor. The specific steps are as follows: first, the SAR-SIFT algorithm is used to extract the key points on the reference image and the image to be registered, and then the set SH,p and S M,p For all points on the reference image and the image to be registered, Voronoi polygons are established as objects, and the local affine transformation matrix of each SAR-SIFT key point is obtained. After transforming the image using the local affine transformation matrix, SAR-SIFT descriptors are established. Next, initial matching point pairs are found using the nearest neighbor principle within each corresponding small polygon in the reference image and the image to be registered. The outlier removal based on length and slope (LSOR) algorithm is used to preliminarily remove mis-matched point pairs from the set of matching point pairs on the entire image, and then a secondary screening is performed through the FSC algorithm to obtain the set of matching point pairs S A,c , and then the LSS-Flow algorithm is used for fine registration, and finally the set of matching point pairs S A,p .

[0168] (4) The set of point pairs S S,p , S H,p , S M,p and S A,p are summarized and the LSOR algorithm is used to remove the possible mis-matched point pairs among them, and finally the set of matching point pairs S final is obtained. Using S final to calculate the global transformation model.

[0169] This embodiment makes up for the deficiencies of existing SAR image registration methods and proposes to use a three-stage registration framework to register SAR images with large viewing angle differences. In the first-stage task, the multi-scale MSER algorithm, the multi-scale Harris-Affine algorithm, and the SAR-SIFT algorithm are used to preliminarily find the matching point pairs between the reference image and the image to be registered, and then a fine registration algorithm is used to perform fine registration on all the obtained matching point pairs. In the second-stage task, by constructing Voronoi polygons, the A-SAR-SIFT algorithm with affine invariance is used for local matching. After obtaining a large number of matching point pairs, the mis-matched point pair removal algorithm and the FSC algorithm are used to fully retain the correct matching point pairs, and finally a fine registration algorithm is used for fine registration. In the third-stage task, all the fine-registered matching point pairs obtained in the first stage and the second stage are screened and the global transformation model is calculated. This method deeply considers the influence of speckle noise existing in SAR images and the severe geometric distortion and radiometric distortion existing between images, and can significantly improve the registration performance of SAR images with large viewing angle differences.

[0170] In the SAR image registration method of this embodiment, in steps 1a), 1b) and (3), the SAR-SIFT algorithm is used to extract the key points in the image, and the operation steps are as follows:

[0171] 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.

[0172] 2b) From the scale sequence of α, use the Ratio of Exponent Weight Average (ROEWA) algorithm to calculate the horizontal gradient G x,α and vertical gradient G y,α of the image at different scale factors, so as to obtain the gradient magnitude Mag α and gradient direction Ori α of the image at different scales:

[0173]

[0174]

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

[0176]

[0177] Among them, the parameter α is the scale parameter of the exponentially weighted function, is the standard deviation of the Gaussian function.

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

[0179] R SH (x, y, α) = det(C SH (x, y, α)) - d · tr(C SH (x, y, α)) 2

[0180] Among them, d is a parameter of any size, generally between 0.04 and 0.06.

[0181] 2e) In the same layer of the scale - space image, compare the SAR - Harris response value of each point with the response values of the points in its 8 - neighborhood and a global threshold d SH respectively. If the response value of the central point in the neighborhood is the largest and greater than the global threshold d SH, then this point is the detected key point, where the global threshold d SH is generally taken as 0.85.

[0182] In the SAR image registration method of this embodiment, in steps 1a), 1b), 1c) and (3), SAR-SIFT descriptors are established at the found key points, and the operation steps are as follows:

[0183] 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. In the circular neighborhood, 0 to 360° is evenly divided into 36 parts, each part representing a range of 10°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. All pixel points in this neighborhood are traversed. 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 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.

[0184] 3b) The direction angle corresponding to the peak value of the smoothed histogram is used as the main direction of this key point, and the direction angle corresponding to the column greater than 80% of the peak energy is used as the secondary direction.

[0185] 3c) With the key point as the center, a circular neighborhood with a radius of r max = 12α is taken, and the circular neighborhood is rotated with the main direction angle as the reference to make the feature descriptor have rotational invariance, and then a gradient direction histogram is established. As shown in the appendix Figure 2 , this circular neighborhood is divided into 17 sub-regions, and the concentric circle radii from the inside to the outside are 0.25·r max , 0.75·r max and r max , and 0 to 360° is evenly divided into 8 parts. Calculate the gradient direction histogram of each sub-region, and splice and normalize the histogram vectors of the 17 sub-regions to generate a 136-dimensional SAR-SIFT descriptor.

[0186] In the SAR image registration method of this embodiment, in steps 1a), 1b) and 1c), the NNDR algorithm is used to find the set of initial matching point pairs between the reference image and the image to be registered, and the operation steps are as follows:

[0187] 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. Find the key point closest and the second closest to the key point. The closest distance and the second closest distance are represented by d min and d nd respectively.

[0188] 4b) If then this key point and the key point closest to it are a pair of correctly matched points. When the threshold distRatio is set to 0.8, a smaller set of point pairs C is obtained h , and when the threshold distRatio is set to 0.9, a larger set of point pairs C is obtained l .

[0189] In the SAR image registration method of this embodiment, in steps 1a), 1b), 1c) and (3), the FSC algorithm is used to remove mismatched point pairs, and the operation steps are as follows:

[0190] 5a) Set the number of iterations of the algorithm to N. In the t-th iteration process, randomly select three pairs of matching points from the set of point pairs C h :

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

[0192] 5c) Use the obtained transformation model to calculate the transformation error of the point pair c l in the set of point pairs C i :

[0193] c i ={(x i , y i ), (x i , y i )}

[0194] e(c i , θ t ) = ||(x i , y i ) - T((x i , y i ), θ t )||

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

[0196] 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 .

[0197] 5e) When the algorithm iterates N times and ends, at this time, C1, C2,..., C N are obtained. There are a total of N sets of point pairs. Take out the set with the largest number of point pairs as the set of matching point pairs finally obtained by the algorithm.

[0198] In the SAR image registration method of this embodiment, in step 1b), the affine transformation matrix at the key point is calculated by using the multi-scale Harris-Affine algorithm, and its operation steps are as follows:

[0199] 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.

[0200] 6b) In the k-th iteration process, use the shape adaptive matrix U (k) to transform the reference image I1 and the key point (x i , y i ) on it to obtain the transformed image and the key point coordinates:

[0201]

[0202]

[0203] where T is a mapping operation that uses the shape adaptive matrix U( k ) to transform the image or coordinates;

[0204] 6c) Centered on , take a square area W with a side length of 4α on the image , where α represents the scale of the image layer where the key point is located.

[0205] 6d) Use the ROEWA algorithm to calculate the horizontal gradient G x,α and the vertical gradient G y,α of the square area W at the scale factor α, so as to obtain the gradient magnitude Mag α of W:

[0206]

[0207] Calculate the SAR-Harris matrix C at each pixel point in W according to the calculated gradient SH (x, y, α):

[0208]

[0209] Among them, the parameter α is the parameter of the exponential weighting function, is the standard deviation of the Gaussian function;

[0210] Use C SH (x, y, α) to calculate the SAR-Harris response value R at each pixel point in W SH (x, y, α)

[0211] R SH (x, y, α) = det(C SH (x, y, α)) - d · tr(C SH (x, y, α)) 2

[0212] Among them, d is a parameter of any size, generally between 0.04 and 0.06.

[0213] Take the point with the largest SAR-Harris response value in the square region W as the new key point coordinates

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

[0215]

[0216] 6f) Take as the center, and reselect a square region W_new with a side length of 4α on the image and calculate the SAR-Harris matrix of the central pixel point in the new region

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

[0218]

[0219] U (k+1) = (μ (k) ) -1 U (k )

[0220] And normalize the shape adaptive matrix so that its maximum eigenvalue is equal to 1.

[0221] 6h) Calculate the convergence rate at the k-th iteration:

[0222]

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

[0224] 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 the convergence rate ratio at the key point (x i , y i ) using the updated shape adaptive matrix; if the condition ratio < 0.1 cannot be satisfied after K iterations, discard this key point;

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

[0226] In the SAR image registration method of this embodiment, in steps 1b), 1c) and (3), the image is transformed using the affine transformation matrix at each key point. For a 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 α , the gradient direction map Ori α and the key point coordinates:

[0227] Mag α T = T(Mag α , Matrix i )

[0228] Ori α T = T(Ori α , Matrix i )

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

[0230] In the SAR image registration method of this embodiment, in step 1c), the multi-scale MSER algorithm is used to extract key points in the image and calculate the affine transformation matrix at each key point. The operation steps are as follows:

[0231] 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 = 21 3 , n = 4, which is the scale of the Gaussian kernel. Select the window length w = 4α, and use the Gaussian blur kernel to establish the scale space;

[0232] 8b) Sort each layer of the image according to the gray value. If it is a color image, the color image needs to be converted into a gray image; assign a node to each pixel point in the image in advance, and the node index number is the gray 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 regions. After all pixel points are placed into the component tree, all the extreme value regions corresponding to the image are obtained. Among them, the extreme value region is defined as: if the gray values of all pixels in a certain region are greater than the gray values of its boundary pixels, then this region is defined as the maximum extreme value region; if the gray values of all pixels in the region are less than the gray values of its boundary pixels, it is defined as the minimum extreme value region;

[0233] 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 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 the area of the region, the subscript i ∈ [0, 255] represents the gray level, and Δ represents a small gray change amount;

[0234] 8d) Approximately fit the irregular maximum stable value region into an elliptical region.

[0235] First, take the centroid of the maximum stable value region as the center of the ellipse and calculate the center of the ellipse:

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

[0237] m 00 = ΣI e (x,y)

[0238] m 01 = ΣyI e (x,y)

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

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

[0241]

[0242]

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

[0244]

[0245] where

[0246] μ 20 = ∑(x - x c ) 2 I e (x,y)

[0247] μ 02 = ∑(y - y c ) 2 I e (x,y)

[0248] μ 11 = ∑(x - x c )(y - y c )I e (x,y)

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

[0250]

[0251]

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

[0253]

[0254]

[0255]

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

[0257]

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

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

[0260] In the SAR image registration method of this embodiment, in steps (2) and (3), the LSS - FLOW algorithm is used to perform fine registration on the point - pair set to obtain a more accurate point - pair set. The operation steps are as follows:

[0261] 9a) Calculate the transformation model θ from the to - be - registered image to the reference image using the matching point - pairs in the reference point - pair set . For steps (2) and (3), the reference point - pair set is the sum of S S,c , S H,c and S M,c . The reference image is I1, and the to - be - registered image is I2. Deform the to - be - registered image I2:

[0262]

[0263] where represents the image obtained by transforming the to - be - registered image using the transformation model θ.

[0264] 9b) For the points of the to - be - finely registered point - pair set S coarse on the to - be - registered image I2, calculate their corresponding positions on the transformed image . Among them, for step (2), the to - be - finely registered point - pair set is S S,c , S H,c and S M,c . For step (3), the to - be - finely registered point - pair set is S A,c :

[0265] (x i T , y i T)=T((x i ,y i ),θ).

[0266] 9c) Use the ROEWA algorithm to calculate the image I1 and The horizontal gradient G at the scale factor α = 2 x With vertical gradient G y , thus obtaining the gradient amplitude Mag and gradient direction Ori of the image:

[0267]

[0268]

[0269] For image I1 and Take a pixel point on the circle, take the pixel point as the center, take a circular neighborhood with a radius of r = 24, and divide this neighborhood into 17 sub-areas, as shown in the figure. Figure 2 As shown in the figure, the radii of the concentric circles are 0.25·r, 0.75·r and r from the inside to the outside respectively. In each sub-region, 0~360° is evenly divided into 8 parts, each representing a range of 45°. The horizontal axis represents the gradient direction angle, and the vertical axis represents the gradient amplitude. All the pixels in the sub-region are traversed. For each pixel in the region, the histogram amplitude column corresponding to the gradient direction angle of the pixel point is first found, and then the gradient amplitude of the pixel point is accumulated on the histogram column of the direction angle of the point, so as to obtain the gradient direction histogram of the point, and convert it into a one-dimensional histogram vector. The histogram vectors of 17 sub-regions are spliced ​​and normalized to generate a 136-dimensional SAR-SIFT descriptor.

[0270] Calculate image I1 and The SAR-SIFT descriptor of each pixel on the image is concatenated to form the image I1 and The SAR-SIFT description sub-image I1_desc and

[0271] 9d) For the set S of point pairs to be precisely matched c o arse A matching point (x i ,y i ), taking it as the center, take out a region with a side length of 2·r on the SAR-SIFT description subgraph I1_desc lf +1 square area describes the sub-image I1_squ, and its corresponding point (x i T ,y i T ) and the obtained square region description sub-image is where r lf take 61;

[0272] 9e) The loss function E(w) is as follows:

[0273]

[0274] where p = (x, y) represents a certain pixel point on the image, w(p) = (u(p), v(p)) represents the optical flow vector at point p on the image, 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; the parameters η and α are taken as 0.001 and 0.01 respectively, and the parameters t and d are both taken as 1.

[0275] Among them, equation (1) is the data term, which constrains that the SAR - SIFT descriptors along the optical flow vector w(p) should match each other; equation (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; equation (3) is the smooth term, which constrains that the optical flow vectors of adjacent pixels should be similar.

[0276] 9f) Calculate the new coordinates of the point (x on the image i T , y i T ):

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

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

[0279] 9g) Use the transformation model θ of the image to calculate the coordinate position of the point (x i T _new, y i T _new) on the original image to be registered I2:

[0280] (x i _new, y i _new) = T((x i T _new, yi T _new), θ -1 )

[0281] 9h) After traversing each pair of matching points in the set S of points to be precisely matched, a more accurate set of point pairs S is obtained coarse Among them, for step (2), the more accurate set of point pairs is S precise , where, for step (2), the more accurate set of point pairs is S S,p , S H,p and S M,p , for step (3), the more accurate set of point pairs is S A,p .

[0282] In the SAR image registration method of this embodiment, in step (3), taking S H,p and S M,p as objects to establish Voronoi polygons, and obtaining the local affine transformation matrix of each SAR-SIFT key point. The operation steps are as follows:

[0283] 10a) Using all the discrete points with affine invariance on the reference image I1 in the set of point pairs S H,p and S M,p , automatically construct a Delaunay triangulation network with these discrete points as objects; number these discrete points and the formed triangles, and record which three discrete points each triangle is composed of.

[0284] 10b) Find out and record the numbers of all the triangles adjacent to each discrete point.

[0285] 10c) Sort the triangles adjacent to each discrete point in clockwise order and record them for the next step of connecting to generate Voronoi polygons.

[0286] 10d) Calculate the circumcenter of the circumcircle of each triangle and record it.

[0287] 10e) According to the adjacent triangles of each discrete point, connect the circumcenters of these adjacent triangles, and then the Voronoi polygons can be obtained.

[0288] 10f) Assign the same affine transformation matrix of the discrete points with affine invariance in the polygon to the SAR-SIFT key points located in each small Voronoi polygon, and then the local affine transformation matrix of each SAR-SIFT key point is obtained.

[0289] 10g) Using the set of point pairs S H,p and S M,pPerform the same operation as above on all discrete point pairs with affine invariance on the image I2 to be registered, and obtain the local affine transformation matrix of each SAR-SIFT key point.

[0290] In the SAR image registration method of this embodiment, in steps (3) and (4), the LSOR algorithm is used to preliminarily remove the mismatched point pairs among the matching point pairs on the entire image, and the operation steps are as follows:

[0291] 11a) Use the reference point pair set to calculate the transformation model θ' from the image to be registered to the reference image, where, for step (3), the reference point pair set is the sum of S S,p 、S H,p and S M,p ; for step (4), the reference point pair set is the sum of S S,p 、S H,p 、S M,p and S A,p ;

[0292] The reference image is I1, and the image to be registered is I2. Deform the image to be registered I2 using the transformation model θ' to obtain the image

[0293]

[0294] where represents the image after transforming the image to be registered using the transformation model θ';

[0295] 11b) For a certain point (x ,y s ) on the image to be registered I2 in the reference point pair set s , calculate its corresponding position coordinates (x s T ,y s T ) on the transformed image ) using the transformation model θ':

[0296] (x s T ,y s T ) = T((x s ,y s ), θ')

[0297] Perform the same transformation on all point pairs in the reference point pair set to obtain a new point pair set S ref T ;

[0298] 11c) Calculate the point pair set Sref T For each pair of matching points, calculate the length and slope of the connecting line, and then sum them up and average to obtain the average length dist_ave and the average slope slope_ave;

[0299] 11d) For the set S of point pairs to be screened fil For a certain point (x i , y i ) on the image I2 to be registered, use the transformation model θ' to calculate its corresponding position coordinates (x on the transformed image i T , y i T ). Among them, for step (3), the set of point pairs to be screened is the set of initial matching point pairs on the entire image found by the nearest neighbor principle in the previous step; for step (4), the set of point pairs to be screened is the sum of S S,p , S H,p , S M,p and S A,p :

[0300] (x i T , y i T ) = T((x i , y i ), θ')

[0301] 11e) Calculate the length dist i and slope slope i of the connecting line between the point pair (x i T , y i T ) and (x i ; i ;

[0302] 11f) Use the threshold to screen all point pairs in the set S fil of point pairs to be screened, and retain the matching point pairs that meet the conditions:

[0303] |dist i - dist_ave| < Th d

[0304] |slope i - slope_ave| < Th s

[0305] where Th d is taken as 0.1, and Th s is taken as 5 degrees.

[0306] In the SAR image registration method of this embodiment, the horizontal gradient G and the vertical gradient G of the image at different scale factors are calculated using the ROEWA algorithm in steps 2b), 6d), and 9c). The operation steps are as follows: x,α and the vertical gradient G y,α , and the operation steps are as follows:

[0307] 12a) Calculate the horizontal gradient G x,α : For any pixel point (i, j), first calculate the exponentially weighted mean of the gray values of the pixel points within the range of (4α + 1)×2α on the left and right sides of this pixel point and Then divide by and take the logarithm to obtain the horizontal gradient G x,α , where α is the exponential weighting factor. The calculation formula of the horizontal gradient G x,α is as follows:

[0308]

[0309]

[0310]

[0311] where I(·) represents the gray value of the pixel point in the SAR image.

[0312] 12b) Calculate the vertical gradient G y,α : For any pixel point (i, j), first calculate the exponentially weighted mean of the gray values of the pixel points within the range of 2α×(4α + 1) above and below this pixel point and Then divide by and take the logarithm to obtain the vertical gradient G y,α , where α is the exponential weighting factor. The calculation formula of the vertical gradient G y,α is as follows:

[0313]

[0314]

[0315]

[0316] where I(·) represents the gray value of the pixel point in the SAR image.

[0317] To address the deficiencies of existing SAR image registration methods, the present invention proposes a three-stage registration framework for registering SAR images with large view angle differences. In the first-stage task, the multi-scale MSER algorithm, the multi-scale Harris-Affine algorithm, and the SAR-SIFT algorithm are used to preliminarily find the matching point pairs between the reference image and the image to be registered, and then a fine registration algorithm is used to perform fine registration on all the obtained matching point pairs. In the second-stage task, the A-SAR-SIFT algorithm with affine invariance is used to perform local matching by constructing Voronoi polygons. After obtaining a large number of matching point pairs, the false matching point pair removal algorithm and the FSC algorithm are used to fully retain the correct matching point pairs, and finally a fine registration algorithm is used for fine registration. In the third-stage task, all the matching point pairs obtained in the first stage and the second stage are screened, and then the global transformation model is calculated. This method takes into account the speckle noise in SAR images and the impacts brought by severe geometric distortion and radiometric distortion between images, and can significantly improve the registration performance of SAR images with large view angle differences.

[0318] Two publicly available datasets (from the paper "KAZE-SAR: SAR Image Registration Using KAZE Detector and Modified SURF Descriptor for Tackling Speckle Noise", Mohammadreza Pourfard et al., 2021) are selected. To further verify the performance of the proposed registration method for SAR images with large view angle differences in the present invention, the images to be registered in the two publicly available datasets are deformed to simulate the registration environment under large view angle differences.

[0319] Figures 4 - 7 These are the two sets of registration images used in the experiments of the present invention. Among them, Figure 4 represents the input reference image of the first group of experiments, Figure 5 represents the input image to be registered in the first group of experiments; Figure 6 represents the input reference image of the second group of experiments, Figure 7 represents the input image to be registered in the second group of experiments.

[0320] Figures 8 - 11 These are the registration results of the present invention on the two sets of registration images. Among them, Figure 8 represents the checkerboard image after registration in the first group of experiments, Figure 9 represents the correct matching point pairs found in the two images in the first group of experiments, Figure 10 represents the checkerboard image after registration in the second group of experiments, Figure 11 represents the correct matching point pairs found in the two images in the second group of experiments.

[0321] Table 1 below presents the registration results of the method of the present invention on two datasets and compares them with the existing registration method, 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). In Table 1, Precision represents the precision rate, that is, the percentage of correct point pairs among the finally obtained point pairs; RMSE represents the root mean square error. The reference image and the image to be registered input for Test1 are Figure 4 and Figure 5 , respectively. The reference image and the image to be registered input for Test2 are Figure 6 and Figure 7 , respectively. It can be seen from the results that the method proposed by the present invention has better registration performance.

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

[0323]

[0324] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. All contents that do not depart from the technical solution of the present invention should be included within the protection scope of the present invention.

Claims

1. A registration method for SAR images with large view angle differences, characterized in that, It includes the following steps: (1) The key points between the reference image and the image to be registered are extracted by using the SAR-SIFT algorithm, the multi-scale Harris-Affine algorithm, and the multi-scale MSER algorithm respectively. Descriptors are established and matched, and finally three sets of matching point pairs with different properties are obtained: the set of point pairs S that is not easily affected by perspective changes S,c , the set of corner point pairs S with affine invariance H,c and the set of point pairs S with regions of affine invariance M,c ; (2) Use the LSS-Flow algorithm to perform fine matching on the point pair sets S S,c , S H,c and S M,c to obtain a more accurate point pair set S S,p , S H,p and S M,p ; (3) Extract key points on the reference image and the image to be registered using the A-SAR-SIFT algorithm and establish descriptors. Next, use the nearest neighbor principle to find initial matching point pairs within each corresponding small polygon in the reference image and the image to be registered. Use the LSOR algorithm to preliminarily remove false matching point pairs from the set of matching point pairs on the entire image, and then perform a secondary screening through the FSC algorithm to obtain the set of matching point pairs S A,c , and then use the LSS-Flow algorithm for fine registration to finally obtain the set of matching point pairs S A,p ; (4) Aggregate the point pair sets S S,p , S H,p , S M,p and S A,p , and then use the LSOR algorithm to remove the potentially included mismatched point pairs, finally obtaining the set of matching point pairs S final . Use S final to calculate the global transformation model.

2. The registration method for large-view-angle-difference SAR images according to claim 1, wherein The point pair set S that is not easily affected by perspective changes S,c , the corner point pair set S with affine invariance H,c , the point pair set S of the region with affine invariance M,c The specific steps are as follows: 1a) Set S of point pairs that are not easily affected by changes in the viewing angle S,c : Extract key points in the image using the SAR-SIFT algorithm and establish SAR-SIFT descriptors. Use the NNDR algorithm to find the initial set of matching point pairs between the reference image and the image to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set of matching point pairs S S,c ; 1b) Set S of corner point pairs with affine invariance H,c : Extract key points in the image using the multi-scale Harris-Affine algorithm and obtain the affine transformation matrix at each key point. The steps of extracting key points in the image are the same as those of the SAR-SIFT algorithm. Next, use the local affine transformation matrix to transform the image at each key point, then establish the SAR-SIFT descriptor. Use the NNDR algorithm to find the set of initial matching point pairs between the reference image and the image to be registered. After removing the mismatched point pairs through the FSC algorithm, finally obtain the set of matching point pairs S H,c ; 1c) Set S of point pairs with an affine invariant region M,c : The key points in the image are extracted by using the multi-scale MSER algorithm, and the affine transformation matrix at each key point is obtained. After the image is transformed by using the local affine transformation matrix at each key point, the SAR-SIFT descriptor is established. The initial set of matching point pairs between the reference image and the image to be registered is found by using the NNDR algorithm. After removing the mismatched point pairs through the FSC algorithm, the set S of matching point pairs is finally obtained M,c 。 3. The registration method for large-view-angle-difference SAR images according to claim 2, characterized in that In step (3), the A-SAR-SIFT algorithm is used to extract key points on the reference image and the image to be registered and establish descriptors. The specific steps are as follows: First, the SAR-SIFT algorithm is used to extract key points on the reference image and the image to be registered, and then, taking all the point pairs in set S H,p and S M,p as objects, Voronoi polygons are established respectively to obtain the local affine transformation matrix of each SAR-SIFT key point. After transforming the image using the local affine transformation matrix, the SAR-SIFT descriptor is established.

4. The registration method for large-viewpoint-difference SAR images according to claim 3, characterized in that In steps 1a), 1b) and (3), the SAR-SIFT algorithm is used to extract key points in the image, and its operation steps are as follows: 2a) Select a set of exponentially weighted scaling 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) From the scale sequence of α, use the ROEWA algorithm to calculate the horizontal gradient G of the image at different scale factors x,α and the vertical gradient G y,α , thereby obtaining the gradient magnitude Mag α and the gradient direction Ori α : 2c) Calculate the SAR-Harris matrix C at each pixel point in the image according to the already obtained gradient SH (x, y, α): Among them, the parameter α is the scale parameter of the exponential weighting function, which is the standard deviation of the Gaussian function; 2d) Use C SH (x, y, α) to calculate the SAR-Harris response value R at each pixel point in the image SH (x, y, α): R SH (x, y, α) = det(C SH (x, y, α)) - d·tr(C SH (x, y, α)) 2 where d is a parameter of any size, generally between 0.04 and 0.06; 2e) In the same layer of the scale-space image, compare the SAR-Harris response value of each point with the response values of the points in its 8-neighborhood and a global threshold d SH respectively. If the response value of the central point in the neighborhood is the largest and greater than the global threshold d SH , then this point is the detected key point, where the global threshold d SH is generally taken as 0.

85.

5. The registration method for large-view-angle-difference SAR images according to claim 3, characterized in that In steps 1a), 1b), 1c) and (3), a SAR-SIFT descriptor is established at the found key points, and its operation steps are as follows: 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 36 parts, each part representing a range of 10°. 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 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 column in the direction angle where this point is located, so as to obtain the main direction histogram of this key point, and smooth the histogram; 3b) The direction angle corresponding to the peak value of the smoothed histogram gradient is used as the main direction of the key point, and the direction angles corresponding to the bins with more than 80% of the peak energy are used as the secondary directions; 3c) Take a circular neighborhood with a radius of r max = 12α centered at the key point, and rotate the circular neighborhood with the main direction angle as the reference to make the feature descriptor rotation invariant. Then, establish a gradient direction histogram, divide this circular neighborhood into 17 sub-regions, and the concentric circle radii are 0.25·r max , 0.75·r max and r max respectively. Divide 0 - 360° into 8 equal parts on average, calculate the gradient direction histogram of each sub-region, and splice and normalize the histogram vectors of the 17 sub-regions to generate a 136-dimensional SAR-SIFT descriptor.

6. The registration method for large-viewpoint-difference SAR images according to claim 2, wherein In step 1b), the multi-scale Harris-Affine algorithm is used to calculate the affine transformation matrix at the key points, and its operation steps are as follows: 6a) For a key point (x i , y i ) on the reference diagram I1, set the number of iterations K = 15 and initialize the shape adaptive matrix U (1) as the identity matrix E; 6b) During the k-th iteration, use the shape-adaptive matrix U (k) to transform the reference image I1 and the key points (x i , y i ) thereon, obtaining the transformed image and the key point coordinates: where T is a mapping operation that transforms an image or coordinates using the shape-adaptive matrix U (k) for image or coordinate transformation 6c) Taking as the center, a square region W with a side length of 4α is taken on the image , where α represents the scale of the image layer where the key point is located; 6d) Calculate the horizontal gradient G of the square region W at the scale factor α using the ROEWA algorithm x,α and the vertical gradient G y,α to obtain the gradient magnitude Mag of W α : Calculate the SAR-Harris matrix C at each pixel in W according to the calculated gradient SH (x, y, α): Among them, the parameter α is the parameter of the exponentially weighted function, is the standard deviation of the Gaussian function; Using C SH Calculate the SAR-Harris response value R at each pixel point in W using (x, y, α) SH (x, y, α) R SH (x,y,α) = det(C SH (x,y,α)) - d·tr(C SH (x,y,α)) 2 where d is a parameter of any size, generally between 0.04 and 0.06; Take the point with the maximum SAR-Harris response value within the square region W as the new key point coordinates 6e) Update the coordinates of the key points on the original reference image I1: 6f) Centered on in the image, reselect a square region W_new with a side length of 4α, and calculate the SAR-Harris matrix of the central pixel point within the new region 6g) Update the shape adaptive matrix: U (k+1) = (μ (k) ) -1 U (k) And normalize the shape adaptive matrix so that its maximum eigenvalue is equal to 1; 6h) Calculate the convergence rate at the k-th iteration: where λ min (μ (k) ) and λ max (μ (k) ) represent the minimum eigenvalue and the maximum eigenvalue of matrix μ (k) respectively; 6i) Exit the loop when ratio < 0.1, and obtain the affine transformation matrix Matrix at the key point (x i , y i ), where Matrix i = U (k) . Otherwise, increment the iteration count by 1 and continue the loop. Recalculate the convergence rate ratio at the key point (x i , y i ) using the updated shape adaptation matrix. If the condition ratio < 0.1 cannot be satisfied after K iterations, discard this key point; 6j) Perform the same above operations on each key point on the reference image I1 and the image I2 to be registered.

7. The registration method for large-viewpoint-difference SAR images according to claim 2, characterized in that In step 1c), the multi-scale MSER algorithm is used to extract key points in the image and calculate the affine transformation matrix at each key point, and its operation steps are as follows: 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 , n = 4, which is the scale of the Gaussian kernel. Select the window length w = 4α i , and establish a scale space using the Gaussian blur kernel; 8b) Sort each layer of the image according to the gray value. If it is a color image, it is necessary to convert the color image into a gray image; a node is assigned to each pixel point in the image in advance, and the node index number is the gray value corresponding to the pixel point; according to the sorting result of the pixel points, they are placed 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 regions; after all pixel points are placed into the component tree, all the extreme value regions corresponding to the image are obtained; among them, the extreme value region is defined as: if the gray values of all pixels in a certain region are greater than the gray values of its boundary pixels, then this region is defined as the maximum extreme value region; if the gray values of all pixels in the region are less than the gray values of its boundary pixels, then it is defined as the minimum extreme value region; 8c) Use the maximum stability determination condition to obtain MSER regions: If Q1,...,Q i-1 ,Q i ,... is a series of mutually inclusive extreme regions, that is If the extreme region is the maximum stable extreme 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 the area of the region, the subscript i ∈ [0, 255] represents the gray level, and Δ represents a small gray level change amount; 8d) Approximately fit the irregular maximally stable extremal regions into an elliptical region; First, use the centroid of the maximally stable extremal region as the center of the ellipse to calculate the center of the ellipse: Calculate the geometric zero-order moment and geometric first-order moment of the maximally stable extremal region: m 00 = ∑I e (x,y) m 01 = ∑yI e (x, y) m 10 = ∑xI e (x, y) 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 region, and I e (x, y) represents the maximally stable 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: Calculate the geometric second-order moment of the maximally stable extremal region: 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) Calculate the two eigenvalues of the geometric second-order moment: Calculate the major semi-axis w, minor semi-axis l of the ellipse and the major axis direction 8e) Using the major semi-axis w, minor semi-axis l of the ellipse and the major axis direction Calculate the affine transformation matrix at the key point (x c , y c ): 8f) Invert the gray values of the original image and repeat the above operations; 8g) Perform the same above operations on each layer of the image in the scale space.

8. The registration method for large-viewpoint-difference SAR images according to claim 3, wherein In steps (2) and (3), the LSS-FLOW algorithm is used to refine the point pair set to obtain a more accurate point pair set, and its operation steps are as follows: 9a) Calculate the transformation model θ from the to - be - registered image to the reference image using the matching point pairs in the set of reference point pairs For steps (2) and (3), the set of reference point pairs is the sum of S S,c , S H,c and S M,c . The reference image is I1 and the to - be - registered image is I2. Deform the to - be - registered image I2: wherein denotes the image obtained by transforming the image to be registered using the transformation model θ; 9b) For the points on the image to be registered I2 in the set S of pairs of points to be precisely matched coarse calculate their corresponding positions on the transformed image where, for step (2), the set of pairs of points to be precisely matched is S S,c S H,c and S M,c and for step (3), the set of pairs of points to be precisely matched is S A,c : 9c) Calculate the horizontal gradient G of image I1 and at the scale factor α = 2 x and the vertical gradient G y , thus obtaining the gradient magnitude Mag and the gradient direction Ori of the image: For the image I1 and a pixel point on it, taking this pixel point as the center, a circular neighborhood is taken and this neighborhood is divided into 17 sub-regions. The concentric circle radii are 0.25·r, 0.75·r, and r from the inside to the outside respectively. Within each sub-region, 0 to 360° is evenly divided into 8 parts, and each part represents a range of 45°. The abscissa represents the gradient direction angle, and the ordinate represents the gradient magnitude. All pixel points within this sub-region are traversed. For each pixel point within the region, 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 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 a one-dimensional histogram vector. The histogram vectors of the 17 sub-regions are spliced together and normalized to generate a 136-dimensional SAR-SIFT descriptor; Calculate the SAR-SIFT descriptors for each pixel point on the image I1 and splice them together to form the SAR-SIFT descriptor map I1_desc of the image I1 and and and the SAR-SIFT descriptor map I1_desc of 9d) For the set S of pairs of points to be precisely matched coarse For a matching point (x i , y i ) on the reference image I1 in the set S, a square region descriptor image I1_squ with a side length of 2·r lf +1 is taken from the SAR - SIFT descriptor graph I1_desc. The same operation is performed on its corresponding point on the image to be registered, and the resulting square region descriptor image is where r lf is taken as 61; 9e) The loss function E(w) is: Among them, \(p=(x,y)\) represents a certain pixel point on the image, and \(w(p)=(u(p),v(p))\) represents the optical flow vector of point \(p\) on the image, 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; \(\varepsilon\) represents the region centered at point \(p\) plus its adjacent 8 points, and \(q\) represents a certain point in this region; the parameters \(\eta\) and \(\alpha\) are taken as 0.001 and 0.01 respectively, and the parameters \(t\) and \(d\) are both taken as 1; Among them, equation (1) is the data term, which constrains the SAR - SIFT descriptors along the optical flow vector \(w(p)\) to match each other; equation (2) is the small displacement term, which constrains the optical flow vector to be as small as possible when there is no other available information; equation (3) is the smoothing term, which constrains the optical flow vectors of adjacent pixels to be similar; Use the belief propagation algorithm to solve for \(w\) when the loss function \(E(w)\) takes the minimum value; 9f) Calculate the image of the point with new coordinates: 9g) Calculate the coordinate position of the point on the original image I2 to be registered: 9h) Traverse the set S of pairs of points to be precisely matched coarse After each pair of matching points in it, a more accurate set S of pairs of points is obtained precise , where, for step (2), the more accurate set of pairs of points is S S,p 、S H,p and S M,p , for step (3), the more accurate set of pairs of points is S A,p .

9. The registration method for large-viewpoint-difference SAR images according to claim 3, characterized in that In step (3), taking S H,p and S M,p as objects to establish Voronoi polygons, and obtaining the local affine transformation matrix of each SAR-SIFT key point. The operation steps are as follows: 10a) Using the set of point pairs S H,p and S M,p All discrete points with affine invariance on the reference image I1 are used to automatically construct a Delaunay triangulation with these discrete points as objects; these discrete points and the formed triangles are numbered, and it is recorded which three discrete points each triangle is composed of; 10b) Find and record the numbers of all triangles adjacent to each discrete point; 10c) Sort and record the triangles adjacent to each discrete point in the clockwise direction for the next step of connecting to generate the Voronoi polygon; 10d) Calculate the circumcenter of each triangle and record it; 10e) According to the adjacent triangles of each discrete point, connect the circumcenters of these adjacent triangles, and the Voronoi polygon can be obtained; 10f) Assign the same affine transformation matrix of the discrete point with affine invariance within each small Voronoi polygon to the SAR - SIFT key points located within the polygon, then the local affine transformation matrix of each SAR - SIFT key point is obtained; 10g) Using the set of point pairs S H,p and S M,p Perform the same operations as above on all discrete point pairs with affine invariance on the image I2 to be registered, and obtain the local affine transformation matrix of each SAR-SIFT key point.

10. The registration method for large-viewpoint-difference SAR images according to claim 3, characterized in that In steps (3) and (4), the LSOR algorithm is used to preliminarily remove the mismatched point pairs on the entire image, and the operation steps are as follows: 11a) Set of reference point pairs Calculate the transformation model θ' from the image to be registered to the reference image. For step (3), the set of reference point pairs is S S,p , S H,p and S M,p ; for step (4), the set of reference point pairs is the sum of S S,p , S H,p , S M,p and S A,p ; The reference image is I1, and the image to be registered is I2. The image to be registered I2 is deformed using the transformation model θ' to obtain an image 11b) For a set of reference point pairs For a certain point (x s , y s ) on the image I2 to be registered in the set, use the transformation model θ' to calculate its corresponding position coordinates (x ), y s T , y s T ) on the transformed image: (x s T ,y s T ) = T((x s ,y s )), θ') The set of reference point pairs Perform the same transformation as above on all point pairs in it to obtain a new set of point pairs S ref T ; 11c) Calculate the lengths and slopes of the lines connecting each pair of matching points in the set S of point pairs, and after summing them up and averaging, obtain the average length dist_ave and the average slope slope_ave; ref T ​ 11d) For the set S of point pairs to be screened fil For a certain point (x i , y i ) on the image I2 to be registered, use the transformation model θ' to calculate its corresponding position coordinates (x , y i T , y i T ) on the transformed image. Among them, for step (3), the set of point pairs to be screened is the set of initial matching point pairs on the entire image found by the nearest neighbor principle in the previous step; for step (4), the set of point pairs to be screened is the sum of S S,p , S H,p , S M,p and S A,p : (x i T ,y i T ) = T((x i ,y i )), θ') 11e) Calculate the length dist i of the line connecting the point pairs (x i , y i T ) and (x i T , y i ), and the slope slope i ; 11f) Use a threshold to screen all the point pairs in the set S of point pairs to be screened, and retain the matching point pairs that meet the conditions: fil ​ |dist i -|dist_ave| < Th d |slope i -slope_ave|<Th s Among them, Th d Take 0.1, Th s Take 5 degrees.

Images