CN104318548A - Rapid image registration implementation method based on space sparsity and SIFT feature extraction - Google Patents

Abstract

The invention discloses a fast image registration implementation method based on spatial sparsity and SIFT feature extraction, which mainly overcomes feature mismatching caused by feature points that are often extracted in smooth areas and texture areas in the classic SIFT feature extraction algorithm The problem. The steps include: 1) Extracting the sparse regions of the reference image and the image to be matched respectively. 2) Extract SIFT feature points for the sparse regions of the reference image and the image to be matched. 3) Roughly match the SIFT feature point sets extracted from the reference image and the image to be matched. 4) Use a random consensus estimation algorithm to filter out false matches in the rough matching results. 5) Using the final matching point pairs of the reference image and the image to be matched, the affine transformation parameters are obtained through affine transformation to realize the registration of the two SAR images. The invention can improve the registration efficiency under the premise of ensuring the accuracy of the classical SIFT feature registration algorithm, and can be used for registration processing of SAR images.

Description

A Fast Image Registration Method Based on Spatial Sparsity and SIFT Feature Extraction

technical field

The invention belongs to the technical field of image processing and relates to image registration, in particular to a rapid image registration implementation method based on spatial sparsity and SIFT feature extraction, which can be used for early registration of SAR image change detection, fusion, splicing, etc. Work.

Background technique

The image formed by synthetic aperture radar (SAR) has the characteristics of all-weather, all-time, high resolution and strong penetrating ability. It is one of the most representative means of earth observation today. Before splicing, fusion, change detection and other operations on SAR images, it is necessary to spatially register images from the same region, at different viewpoints at the same time or at different times of the same sensor to eliminate the time, angle, and The difference in environment and sensor imaging mechanism causes problems such as translation, rotation, stretching and local deformation between images.

As a basic task in remote sensing image processing, the goal of registration is the process of geometrically aligning two or more images obtained at different times, different viewing angles or different sensors. For SAR images, the speckle noise formed in the coherent imaging process affects the gray distribution of the image. In addition, factors such as differences in imaging conditions, shooting time and seasonal changes, and changes in the scene and objects in the shooting scene cause the image to have a complex gray distribution. Due to the stable characteristics, the algorithms and techniques that can be successfully applied in optical remote sensing images are repeatedly frustrated when applied to SAR images.

Currently, a large number of registration methods have been proposed. According to different feature spaces, registration algorithms can generally be divided into two categories: grayscale-based registration methods and feature-based registration methods. Grayscale-based methods generally use the grayscale information of the entire image directly, and measure the matching degree of the surface reflection attributes of the overlapping parts of the two images by establishing a similarity measure between pixels, and then find the optimal matching. Transform parameters such as translation, rotation, and scaling. The image registration method based on grayscale, its similarity measurement function includes cross-correlation function, Fourier function and mutual information function, etc., by using the transformation relationship corresponding to the maximum similarity, the final realization of the image registration is achieved. Registration. When registering SAR images using the grayscale-based registration method, it is easily affected by the grayscale difference and noise caused by the SAR image imaging process. The feature-based image registration method selects salient and prominent objects from the image for detection, extracts the features of these feature points, establishes the transformation relationship between images through the matching of these features, and finally realizes image registration. The feature-based image registration method can solve the registration problem between images with large geometric distortion and gray information difference. Commonly used image features include: feature points (including corner points, high curvature points, inflection points, etc.), straight line segments, edges, contours, closed areas, and feature structures. The feature-based registration method first needs to extract feature points, edges, etc. from the image to be registered, and then search and match corresponding feature pairs or feature sets to obtain image transformation parameters. The registration method based on edge features needs to use some edge detection operators to extract edge features and then match the extracted edges, but often because the extracted edges are too trivial to reflect the structural features in the image well, the subsequent matching Accuracy and speed are affected. Compared with features such as lines and surfaces, point features have many advantages such as easy extraction, fast speed, high precision, and good stability, and have wider applications.

Feature-based image registration has stronger anti-interference ability to noise and deformation, so it has higher robustness. One of the most commonly used feature-based registration methods is the registration method based on SIFT feature extraction. Lowe summarized the main steps of SIFT features as: (1) scale space extreme point extraction; (2) extreme point positioning; (3) main direction assignment; (4) key point description. SIFT has strong rotation and scale invariance at small rotation angles, and is also robust to illumination changes after normalization. However, since SIFT uses a 128-dimensional vector to describe the feature points, the calculation cost and time complexity are increased when there are many feature points, and a better distinction between features is required when looking for the best match. When using the SIFT feature extraction algorithm to process SAR images, not only can more SIFT feature points be extracted in areas with sparse structures such as straight lines, broken lines, and intersections, but also many features can be extracted in smooth areas and texture areas where the features are not obvious. points, but there are many similar points around the feature points of these smooth and textured areas, so it takes a lot of time and the features of these points are not obvious, and it is easy to form a mismatch.

Contents of the invention

The purpose of the present invention is to overcome the deficiencies of the prior art, and proposes a SAR image registration method that combines spatial sparse structure and SIFT feature extraction, to realize the centralized extraction of SIFT feature points in the sparse structure area of the SAR image, and reduce algorithm operations On the basis of cost and time efficiency, the matching error rate is reduced, and the efficiency of image registration is improved.

The technical thinking of realizing the object of the present invention is: the input reference image and the image to be matched obtain the sparse regions of the two images respectively, extract the SIFT feature points for the sparse regions of the reference image and the image to be matched, and utilize the nearest short-distance ratio criterion ( NNDR) as a similarity measure to obtain a rough matching pair, and use the random consistency estimation (RANSAC) method to remove the mismatch, and finally use the affine transformation to obtain the affine transformation parameters, so as to realize the registration of the reference image and the image to be matched .

The technical scheme of the present invention is: a kind of fast image registration realization method based on spatial sparsity and SIFT feature extraction, comprises the following steps:

(1) Input two multi-temporal SAR images of the same area acquired at different times, which are respectively recorded as the reference image and the image to be matched, calculate the sparsity of each pixel in the reference image and the image to be matched, and select the one with the highest sparsity The area where the pixel is located is regarded as a sparse area, which is denoted as Mask _R and Mask _S respectively;

(2) Extract SIFT feature points from the sparse area Mask _R of the reference image and the sparse area Mask _S of the image to be matched, and obtain the SIFT feature point sets of the reference image and the image to be matched as {R _m } and {S _n } respectively;

(3) Perform feature matching on the SIFT feature point sets {R _m } and {S _n } extracted from the reference image and the image to be matched to obtain a rough matching feature pair {(pq) _i }, the reference image and the image to be matched The coarse matching point sets of are respectively {p _i } and {q _i }, i represents the index number of the matching pair;

(4) Use the random consistency estimation algorithm to remove the mismatched pairs from the rough matching results. After removing the mismatched pairs, the set of feature point pairs retained in the reference image and the image to be matched is denoted as {(PQ) _i }. At this time, the reference image and The remaining feature points in the graph to be matched are respectively {P _i } and {Q _i }, where i represents the index number of the matching pair;

(5) Using the matching point sets {P _i } and {Q _i } of the reference image and the image to be matched, affine transformation parameters are obtained through affine transformation, and finally the two images are registered according to the affine transformation parameters.

Above-mentioned step (1) comprises the steps:

1) Carry out two-layer wavelet decomposition on the reference image R and the image to be matched respectively, and take db1 as the wavelet base to obtain the detail component images of the reference image and the image to be matched, which are denoted as R _db1 and S _db1 respectively;

2) Calculate the variance coefficient for each pixel of the detail component map of the reference image and the detail component map of the image to be matched, and then take the current pixel point as the center and take a neighborhood of size m*m to calculate the variance, and the current pixel can be obtained The variance coefficient CV of the point; set a threshold T _cv , which has a linear relationship with the variance coefficient of the SAR image T _cv = β·CV, where β is a linear coefficient; the binarized variance coefficient matrix can be obtained from the threshold T _cv CV_mask:

$\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; mask</span>}<span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; \_</span><span\; class="notranslate">\; <</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\end{array}\right.$

The variance coefficient matrix of the binarization of the reference image and the image to be matched can be calculated, which are respectively denoted as CV_mask _R and CV_mask _S ;

3) Thinning the variance coefficient matrix of the reference image to obtain a thinning matrix T _R , and similarly refining the variance coefficient matrix of the matching image to obtain a thinning matrix T _S ;

4) Calculate the sparse matrix of the reference image and the image to be matched;

5) Select the top 20% of the pixels with the largest sparsity value in the sparse matrix as the sparsity point, and expand the sparse point to obtain the sparse area.

Above-mentioned step 4) comprises the steps:

(a) Calculate the size and direction of the gradient of each pixel in the refinement matrix:

(b) Analyze the direction of the gradient of each pixel point, count the direction characteristics of the current pixel point, and obtain the direction characteristic matrix of the reference image and the image to be matched, which are recorded as _DR and _DS respectively;

(c) From the refinement matrix and the direction feature matrix, the sparsity of the reference image and the image to be matched can be calculated, and the sparsity is defined as:

$\mathrm{<span\; class="notranslate">\; S</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}}\underset{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; D.</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&Center\; Dot;</span>{{\text{<span class="notranslate"> T</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Among them, (i, j) is the index of the current pixel point, Ω is the r×r neighborhood centered on the current pixel point, N is the total number of pixels in the neighborhood Ω, D is the direction feature matrix, and T is defined as fine matrix; thus obtain the sparse matrix S _R and S _S of the reference image and the image to be matched.

Above-mentioned step (b) comprises the steps:

a) Construct a 16-direction template with a size of (2r+1)*(2r+1), and the angle ranges corresponding to each sub-direction of the direction template are respectively i∈[1,16], r is the radius of the direction template;

b) Analyzing the direction characteristics of each thinning edge point in the thinning matrix with 16 direction templates;

c) Normalize the directional feature matrices of the reference image and the image to be matched to obtain the directional feature matrices D _R and D _S of the reference image and the image to be matched.

Beneficial effects of the present invention: the present invention uses the variance coefficient of the SAR image to determine the smooth area, texture area, and sparse structure area of the image, and excludes smooth areas and texture areas that are not obvious and easily cause mismatches, and only for those fold line inflection points Extract SIFT feature points for feature matching in the sparse area where pixels with high sparsity such as intersections and intersections are located, and image registration is realized by calculating affine transformation parameters. Compared with the prior art, the present invention has the following advantages:

1. The present invention adopts the classic SIFT operator. When the Gaussian scale space is established in the process of extracting SIFT feature points, the Gaussian function is used to blur the image while reducing the speckle noise of the SAR image, thereby improving the pixel point in its The distinguishability of the neighborhood reduces the occurrence of false matching to a certain extent.

2. The present invention designs a method for defining the degree of sparsity, which can accurately analyze sparse areas with obvious features such as broken lines, inflection points, and intersections, and exclude non-sparse areas with low sparsity and inconspicuous features such as smooth areas and textured areas .

3. The present invention only extracts SIFT feature points for the sparse areas of the reference image and the image to be matched, and excludes non-sparse areas with low sparsity and inconspicuous features such as smooth areas and texture areas, thereby reducing the probability of false matching and making the algorithm more efficient. fast.

4. The present invention is simple to implement and has a clear structure, and requires a small overlapping area between the reference image and the image to be matched, and the SIFT operator still has high registration accuracy for situations with rotation and scale transformation or situations including changing regions.

Description of drawings

Fig. 1 is a flow chart of the present invention;

Fig. 2 is the sub-flow chart of calculating sparse area among the present invention;

Fig. 3 is two pieces of SAR images that the experiment of the present invention inputs, and Fig. 3 (a) is reference figure, and Fig. 3 (b) is figure to be matched;

Fig. 4 is the sparse area Mask _A and Mask _B of the reference image and the image to be matched;

Figure 5 is the distribution of feature points of the reference image and the image to be matched extracted by the classic SIFT algorithm;

Fig. 6 is to extract the SIFT feature point distribution of the sparse area Mask _A and Mask _B of the reference image and the image to be matched;

Fig. 7 is the registration result of the reference image and the image to be matched obtained according to the affine transformation parameters, displayed by the chessboard method.

detailed description

With reference to Fig. 1, the concrete realization steps of the present invention are as follows:

Step 1. Input two multi-temporal SAR images R and S of the same area acquired at different times. For the convenience of description, we call the image R the reference image and the image S the image to be matched. Calculate the sparse areas of the two images respectively Mask _R and Mask _S. Refer to Figure 2 for this step, and the specific steps are as follows:

(1a) Carry out two-layer wavelet decomposition on the reference image R and the image to be matched respectively, and take db1 as the wavelet base, and obtain the detail components of the reference image and the image to be matched as R _db1 and S _db1 respectively;

(1b) Calculate the variance coefficient CV for each pixel of the reference image detail component R _db1 and the to-be-matched image detail component S _db1 . CV is a measure of the dispersion of a probability distribution, which is defined as the ratio of the standard deviation to the mean:

$\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}}{{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}}$

Among them, σ _x represents the standard deviation of the r×r neighborhood centered on the current pixel point, and _μx represents the mean value of the r×r neighborhood centered on the current pixel point. Then taking the current pixel as the center and taking a neighborhood with a size of m*m to calculate the variance, the variance coefficient CV of the current pixel can be obtained. In the example of the present invention, r=3, m=5. Using the above method, the variance coefficient matrix CV _R of the reference image detail component R _db1 can be obtained; similarly, there is a variance coefficient matrix CV _S for the image detail component S _db1 to be matched. Set a threshold T _cv , which is related to the variance coefficient CV of the SAR image (for the intensity SAR image For magnitude SAR images Where L is the visual number of the SAR image), it is generally regarded empirically as having a linear relationship with the variance coefficient, that is, T _cv =β·CV, and the value of the constant β of the coefficient in the present invention is taken as 0.55. The binarized variance coefficient matrix CV_mask can be obtained:

$\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; mask</span>}<span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; \_</span><span\; class="notranslate">\; <</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\end{array}\right.$

The variance coefficient matrices for the binarization of the reference image and the image to be matched can be calculated, which are denoted as CV_mask _R and CV_mask _S respectively.

(1c) Thinning the reference image CV_mask _R matrix to obtain a thinning matrix T _R , and refining the matching image CV_mask _S matrix to obtain a thinning matrix T _S .

(1d) Sparse matrices for the reference image and the image to be matched, respectively. The specific steps to calculate the sparse matrix are as follows:

1d1) Calculate the gradient magnitude and direction of each pixel in the refinement matrix:

${\mathrm{<span\; class="notranslate">\; m</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\sqrt{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}$

${\mathrm{<span\; class="notranslate">\; k</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; the\; tan</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}}{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; )</span>$

Where (x, y) represents the index of the current pixel in the refinement matrix, d is the gray value of the pixel in the refinement matrix, m _{(x, y)} is the gradient size, and k _{(x, y)} is the gradient direction .

1d2) Construct a 16-direction template with a size of 15×15, and the angle ranges corresponding to each sub-direction of the direction template are respectively i∈[1,16], the total number of pixels in each sub-direction on the direction template is N, and the radius of the direction template is r.

1d3) Using the 16-direction template to analyze the direction characteristics of each thinning edge point in the thinning matrix. In the refinement matrix, take the neighborhood with the current refinement edge point as the center and the size is 15×15, analyze the gradient direction of each pixel in each sub-direction, and count the number of pixels whose gradient direction is consistent with the sub-direction of the template , denoted as s _i , statistically obtain the normalized directional feature vector of the current thinning edge point $\mathrm{hist}=[\frac{{\mathrm{the\; s}}_1}{\mathrm{\&Sigma;}{\mathrm{the\; s}}_i}...\frac{{\mathrm{the\; s}}_i}{\mathrm{\&Sigma;}{\mathrm{the\; s}}_i}...\frac{{\mathrm{the\; s}}_{16}}{\mathrm{\&Sigma;}{\mathrm{the\; s}}_i}].$

1d4) Calculating the directional feature matrix of the reference image and the image to be matched. set a threshold N is the total number of pixels in each sub-direction on the direction template, r is the radius of the direction template, and the value of the threshold coefficient β in the present invention is 1.8. Analyze and refine the direction characteristics of edge points. If h _i ≥ T _hist , this direction is one of the main directions. The number of main directions of each edge point is used as its direction feature, and the directions of the reference image and the image to be matched can be obtained Characteristic matrices D _R , D _S .

1d5) Calculate the sparsity of the reference image and the image to be matched, define a sparse matrix,

$\mathrm{<span\; class="notranslate">\; S</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}}\underset{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; D.</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&Center\; Dot;</span>{{\text{<span class="notranslate"> T</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Among them, (i, j) is the index of the current pixel point, Ω is the r×r neighborhood centered on the current pixel point, N is the total number of pixels in the neighborhood Ω, D is the direction feature matrix, and T is defined as fine matrix. The sparse matrices S _R and S _S of the reference image and the image to be matched can be obtained.

(1e) Select the top 20% of the pixels with the largest sparsity value in the sparse matrix as the sparsity point, and expand the sparse point to obtain the sparse area.

Step 2. Extract SIFT feature points from the sparse regions of the reference image R and the image to be matched respectively. The specific steps of feature point extraction are as follows:

(2a) Construct a Gaussian scale space. The scale space of a two-dimensional image at different scales can be obtained by convolving the original image with the Gaussian kernel function:

L(x,y,σ)=G(x,y,σ)*I(x,y)

Among them: L(x, y, σ) represents the scale space of the image, I(x, y) represents the original image, (x, y) represents the index of the current pixel of the input image, * represents the convolution operation, G(x, y,σ) is a variable scale Gaussian kernel function, σ is a scale space factor, and the kernel function is defined as:

$\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; e</span>}}^{<span\; class="notranslate">\; -</span>\frac{{\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}}$

The scale space is discretized and sampled, that is, the Gaussian kernel under different multiples of σ is convolved with the downsampled image, and a series of Gaussian images generated are called Gaussian pyramid (LOG). Gaussian Pyramid (LOG) is divided into O groups, and each group includes S+3 layers of scale space, and O=4 and S=2 are taken in the present invention.

First find each image in the first set of Gaussian scale spaces. The input image I ₀ is convoluted by the kernel functions of Gaussian kernels of different sizes (the Gaussian kernels are σ _i =k ⁱ σ respectively) to obtain the output image I _i , i=0, 12, 3, 4, which is the first Gaussian scale space A set of images.

The image size of each layer of the second group of Gaussian pyramids is 1/4 of the size of the first group of images, and its input image is obtained by downsampling the S-th layer image in the first group of LoG with a sampling rate of 2. For example, when S=2, the input image I′ ₀ of the second group of Gaussian pyramids is obtained by downsampling the first group of Gaussian pyramids I ₁ , and then undergoes kernel functions of Gaussian kernels of different sizes (Gaussian kernels are σ _i = k ⁱ σ) The convolution obtains the output I′ _i , i=0, 12, 3, 4, that is, the second group of images in the Gaussian scale space is obtained. Similarly, the third and fourth groups of images in the Gaussian scale space can be obtained. This creates a Gaussian scale space.

(2b) Construct a difference-of-Gaussian scale space (DOG). Subtract the scale space functions of two adjacent layers in the Gaussian pyramid (LOG) to obtain the difference of Gaussian pyramid (DOG), represented by D(x,y,σ):

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ)×I(x,y))=L(x,y,kσ)-L(x,y ,σ)

Among them, (x, y) represents the current pixel in the input image, and I(x, y) represents the image to be processed, that is, the difference Gaussian scale space D(x, y, σ) consists of two phases with a difference of a constant factor k It is obtained by subtracting adjacent Gaussian scale space functions.

(2c) Extreme point detection. Compare the size of each pixel in the differential Gaussian scale space (DOG) with 8 pixels in the same layer, 9 pixels in the upper layer, and 9 pixels in the lower layer, a total of 26 pixels. If the current pixel is extremely If the value is large or small, the current pixel point is set as the key point.

(2d) Precise localization of key points. Since the edge and noise have a great influence on the differential Gaussian scale space (DOG), it is necessary to further screen the key points detected in the previous step to accurately locate the feature points with scale invariance, and it is necessary to filter out low-contrast points and Edge response points to improve the robustness and stability of feature matching. In order to accurately determine the scale and position of feature points, three-dimensional quadratic function fitting must be performed on the detected local extremum points. At the local extremum point (x ₀ ,y ₀ ,σ), the Taylor expansion of the scale space function D(x,y,σ) is:

$\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}\frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

In the above formula, X = (x, y, σ) ^T represents the position and scale information of the feature point, X ₀ represents the position and scale information of the feature point to be selected, and the first-order and second-order partial derivatives are calculated as follows,

$\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}\\ \frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}\\ \frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}\end{array}\right]<span\; class="notranslate">\; ,</span>\frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}\frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; xy</span>}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x\&sigma;</span>}}\\ \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; yx</span>}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; y\&sigma;</span>}}\\ \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; \&sigma;x</span>}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; \&sigma;y</span>}}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}\end{array}\right]$

Then take the derivative of this fitting function and set the first derivative to zero, then the extremum x of x can be obtained

$\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; D.</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}}{{<span\; class="notranslate">\; \&PartialD;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}$

Will Substituting the scale space function D(x,y,σ), we can get

$\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\frac{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; D.</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}$

like It means that the contrast of this point is low, and it will be removed at this time.

An important step to obtain stable and accurate feature points is to filter out edge response feature points. When filtering edge response feature points, refer to the following criteria: at the edge intersection and peak of the Gaussian difference space function D(x, y, σ), a larger principal curvature value exists on the feature point of the image edge, but in The value of the principal curvature in the vertical direction is smaller. Among them, the main curvature value can be obtained through the 2×2 Hessian matrix H, Hessian matrix H,

$\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; xx</span>}}& {\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}\\ {\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; yx</span>}}& {\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}\end{array}\right]$

Since the eigenvalues of the H matrix are proportional to the principal curvature of D, the method of not specifically calculating the eigenvalues but only the ratio ratio is adopted. α and β are the largest and smallest eigenvalues respectively, and γ=α/β, then have:

$\mathrm{<span\; class="notranslate">\; ratio</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; Tr</span>}{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; Det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&alpha;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&beta;</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; \&alpha;\&beta;</span>}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}$

The calculation formula of the sum Tr(H) of the main diagonal elements of the matrix H is: Tr(H)=D _xx +D _yy =α+β; the calculation formula of the value Det(H) of the matrix H determinant is: Det( H)=D _xx D _yy −(D _xy ) ² =αβ.

Therefore we only need to detect

$\frac{\mathrm{<span\; class="notranslate">\; Tr</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; Det</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; <</span>\frac{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}$

where γ=10.

(2e) Generate feature descriptors. In order to ensure the rotation invariance of feature points, it is necessary to use the gradient distribution characteristics of the pixel points around the feature point to calculate the main direction of the feature point. The gradient magnitude and direction of each pixel are calculated as follows:

$\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\sqrt{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}$

$\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>$

Among them, the size of the gradient of the feature point at (x, y) is m(x, y), the direction is θ(x, y), and the scale of the feature point (x, y) in the pyramid is L(x, y) To represent. In the neighborhood centered on the current feature point, the gradient direction of the pixels in the neighborhood is counted using the gradient direction histogram. The gradient direction range is [0,360°), where every 10° is a histogram cylinder. There are 36 cylinders in total. The peak of the gradient direction histogram represents the main direction of the neighborhood gradient at the feature point, which is taken as the main direction of the feature point. If there is a peak of more than 80% of the main peak capacity in the gradient direction histogram, the direction corresponding to the peak is used as the auxiliary direction, and each feature point has a main direction and several auxiliary directions.

After determining the main direction of the feature points, the feature descriptor can be generated. Take a 16×16 window centered on the feature point, calculate the gradient direction histogram of 8 directions on the 4×4 image block, and draw the cumulative value of each gradient direction. Gradient information is calculated for each 4×4 sub-block, and 4×4 sub-blocks are generated in a 16×16 neighborhood, so the feature descriptor has a total of 4×4×8=128 dimensions. The feature descriptor consists of the gradient orientation histograms of all sub-blocks, and finally forms a 128-dimensional SIFT feature description vector.

Step 3, perform feature matching on the SIFT feature points of the reference image and the image to be matched, and then use the random consistency estimation (RANSAC) method to filter out the mismatched pairs in the matching results. Specific steps are as follows:

(3a) Perform feature point matching on the SIFT feature points extracted from the sparse regions of the reference image and the image to be matched, and use the nearest neighbor distance to the next nearest neighbor distance ratio (NNDR) criterion to achieve. Suppose A is a certain feature point in the reference image, and D _A is its corresponding feature vector. In the image to be registered, the nearest neighbor (the smallest Euclidean distance) of the feature point A is the feature point B and the second nearest neighbor (the second smallest Euclidean distance) is the feature point _{C, and their corresponding feature vectors are DB and D C ,} and the nearest neighbor The distance to next nearest neighbor distance ratio (NNDR) criterion is defined as:

$\frac{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; B</span>}}<span\; class="notranslate">\; |</span>}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; C</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; dis</span>}}$

Where T _dis is the distance ratio threshold, D _A represents the 128-dimensional feature descriptor of the current feature point A in the reference image, DB represents the 128-dimensional feature descriptor of the feature point B most similar to A in the image to be matched, and _{D C represents} For the 128-dimensional feature descriptor of the feature point C similar to A in the image to be matched, when the ratio between the nearest neighbor Euclidean distance and the second nearest neighbor Euclidean distance is greater than the threshold T _dis , it is considered that the nearest neighbor Euclidean distance corresponds to a A pair of feature points is a pair of matching feature points. The larger the threshold value, the looser the matching requirements, and the higher the probability of wrong matching points; the smaller the threshold value, the stricter the matching requirements, the lower the probability of wrong matching points, and some correct matches will be removed. . Therefore, the threshold value is set according to the actual situation such as the type of image to be registered. After the matching is completed, the matching point sets in the reference image and the image to be matched are respectively recorded as {p _i } and {q _i }, that is, the rough matching pair {(pq) _i } is obtained, and i represents the index number of the matching pair.

(3b) Due to the influence of cluttered background and noise, there are often many wrong matching point pairs in the SIFT matching pairs obtained in the previous step, which will cause errors in the geometric correction of the image. We use the random consistency estimation (RANSAC) method to remove the rough matching results The wrong matching pair in , after removing the wrong matching, the point sets in the matching pair in the reference image and the image to be matched are denoted as {P( _xi ,y _j )} and {Q(x′ _i ,y′ _i )}, respectively, i represents the index number of the matched pair.

Step 4, use the matching point sets {P( _xi ,y _i )} and {Q(x′ _i ,y′ _i )} of the reference image and the image to be matched to obtain accurate affine transformation parameters through the affine transformation function , the radiation transformation model is as follows:

$\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; x</span>}\\ \mathrm{<span\; class="notranslate">\; the\; y</span>}\end{array}\right]<span\; class="notranslate">\; =</span>\left[\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]<span\; class="notranslate">\; \&Center\; Dot;</span>\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; \&prime;</span>}\\ {\mathrm{<span\; class="notranslate">\; the\; y</span>}}^{<span\; class="notranslate">\; \&prime;</span>}\end{array}\right]<span\; class="notranslate">\; +</span>\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]$

The transformation matrix M is:

$\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]$

Wherein, a ₁ , a ₂ , b ₁ , and b ₂ are rotation variables, c ₁ is a horizontal offset, and c ₂ is a vertical offset. Finally, the images are registered according to the solved affine transformation matrix.

Effect of the present invention can further confirm by following experiment:

1. Experimental conditions and content

The hardware platform is: Intel(R) Pentium(R) 1CPU2.4GHz;

The software platform is: WindowXP Professional, MATLABR2010;

The input images used in the experiment are shown in Figure 3, where Figure 3(a) and Figure 3(b) are SAR images.

Experiment content:

Simulation content 1, under the above-mentioned experimental conditions, use the classic SIFT algorithm and the fast registration algorithm in the present invention to extract SIFT feature points respectively, and then according to the nearest short-distance ratio criterion (NNDR), the conditions of different distance ratio thresholds T _dis Next, match the feature points extracted by the two methods.

Simulation content 2, on the basis of the rough matching results obtained, the random consistency estimation algorithm (RANCAC) is used to remove mismatching. Use the last remaining matching pairs to calculate the affine transformation parameters, and finally realize the image registration.

Considering that if the distance is too small than the threshold _Tdis , the number of correct matching pairs in the rough matching result will be too small to meet the requirements of affine transformation; if the distance is too large than the threshold _Tdis , the correct matching in the rough matching result If the proportion of pairs in all matching pairs is too low, the robustness of the random consensus estimation (RANCAC) algorithm is too poor. Therefore, the following values are selected for the distance ratio threshold T _dis in the experiment: 0.70, 075, 0.80, and 0.85.

2. Experimental results

The sparse area obtained by the present invention is shown in Figure 4, Figure 4(a) is the sparse area of the reference image; Figure 4(b) is the sparse area of the image to be matched. It can be seen from the figure that most of the sparse areas are concentrated in areas such as inflection points and intersections of polylines, and these areas have high sparsity and obvious characteristics.

The results of extracting feature points from the reference image and the image to be matched using the classic SIFT feature extraction method are shown in Figure 5. The white dots in the figure represent the extracted feature points, and Figure 5(a) shows the distribution of SIFT feature points on the reference image; 5(b) is the distribution of SIFT feature points of the image to be matched. It can be seen from the result figure that the classic SIFT feature extraction method can detect more SIFT feature points in areas with high sparsity such as the centerline target area and intersection area in the image, and can also detect smooth and textured areas with inconspicuous features. To some feature points, there are often many similar points around the points in these smooth areas and texture areas, so it is easy to form a mismatch.

The results of extracting feature points from the reference image and the image to be matched using the fast image registration method based on spatial sparsity and SIFT feature extraction are shown in Figure 6. The white dots in the figure represent the extracted feature points, and Figure 6(a) is the reference The distribution of feature points of the graph; Figure 6(b) is the distribution of feature points of the graph to be matched. It can be seen from the result figure that the feature points extracted by the fast image registration method based on spatial sparsity and SIFT feature extraction only exist in the inflection point and intersection area of the polyline with high sparsity and obvious features, and the point sparsity in these areas is high , stable features, and not easy to form mismatches.

The feature points extracted by the two methods are roughly matched using the feature matching method based on the Euclidean distance ratio. Then the random consistency estimation method RANSAC is used to remove the wrong match from the rough matching result. The results of the two methods of rough matching and removing false matching are listed in Table 1, and the matching accuracy is used as the evaluation index of the matching effect. Definition of matching accuracy:

Table 1 Results of feature point matching of two methods

It can be found from Table 1 that under different nearest short distance ratio thresholds T _dis =0.70, T _dis =0.75, T _dis =0.80, T _dis =0.85, the proportion of correct matching pairs in all matching pairs in the present invention The ratio is higher than the classic SIFT algorithm. Moreover, the advantage of the present invention over the classic SIFT algorithm is more obvious when the mismatch rate in the rough matching result is high. It can be seen from the statistics of the matching accuracy rate that under the conditions of different Euclidean distance ratio thresholds, the matching accuracy rate of the present invention is higher than that of the classic SIFT algorithm.

The affine transformation parameters are obtained by affine transformation, and the affine transformation parameters obtained by the final matching result when the threshold is T _dis =0.75 are shown in Table 2:

Table 2: Affine transformation parameter results

From table 2, it can be found that the affine transformation parameters obtained by the method of the present invention are not much different from the results of the SIFT algorithm, but the time efficiency has been greatly improved, and the algorithm time in the present invention has been reduced from the original 81.04 seconds to 20.06 seconds, which is shortened to a quarter of the original time, greatly improving the registration efficiency. When the threshold is T _dis =0.75, the registration result of the two images is shown in Figure 7, which is displayed by the checkerboard method. Because the present invention only extracts SIFT feature points from the sparse areas with obvious features such as broken lines, inflection points, and intersection points of the SAR image, and excludes areas with low sparsity such as smooth areas and texture areas, where features are not obvious and are easy to form false matches. On the premise of reducing the false matching rate, the registration efficiency is improved.

The present invention designs a new rapid image registration method based on spatial sparsity and SIFT feature extraction. First, the image sparsity matrix is analyzed, and the SIFT feature points are extracted from the sparse areas with obvious features in the image. Accuracy is improved while improving the efficiency of the algorithm.

To sum up, the present invention uses the variance coefficient of the SAR image to judge the smooth area, texture area, and sparse structure area of the image, and excludes smooth areas and texture areas that are not obvious and likely to cause mismatches, and only for those polyline inflection points and intersections Extract SIFT feature points for feature matching in the sparse area where the pixels with higher sparsity are located, and achieve image registration by calculating affine transformation parameters. Compared with the prior art, the present invention has the following advantages:

1. The present invention adopts the classic SIFT operator. When the Gaussian scale space is established in the process of extracting SIFT feature points, the Gaussian function is used to blur the image while reducing the speckle noise of the SAR image, thereby improving the pixel point in its The distinguishability of the neighborhood reduces the occurrence of false matching to a certain extent.

2. The present invention designs a method for defining the degree of sparsity, which can accurately analyze sparse areas with obvious features such as broken lines, inflection points, and intersections, and exclude non-sparse areas with low sparsity and inconspicuous features such as smooth areas and textured areas .

3. The present invention only extracts SIFT feature points for the sparse regions of the reference image and the image to be matched, and excludes non-sparse regions with low sparsity and inconspicuous features such as smooth regions and texture regions, thereby reducing the probability of false matching and making the algorithm more efficient. fast.

4. The present invention is simple to implement and has a clear structure, and requires a small overlapping area between the reference image and the image to be matched, and the SIFT operator still has high registration accuracy for situations with rotation and scale transformation or situations including changing regions.

The parts that are not described in detail in this embodiment are commonly known and commonly used means in this industry, and will not be described here one by one. The above examples are only illustrations of the present invention, and do not constitute a limitation to the protection scope of the present invention. All designs that are the same as or similar to the present invention fall within the protection scope of the present invention.

Claims (4)

1. A fast image registration implementation method based on spatial sparsity and SIFT feature extraction, is characterized in that: comprise the steps: (1) Input two multi-temporal SAR images of the same area acquired at different times, which are respectively recorded as the reference image and the image to be matched, calculate the sparsity of each pixel in the reference image and the image to be matched, and select the one with the highest sparsity The area where the pixel is located is regarded as a sparse area, which is denoted as Mask _R and Mask _S respectively; (2) Extract SIFT feature points from the sparse area Mask _R of the reference image and the sparse area Mask _S of the image to be matched, and obtain the SIFT feature point sets of the reference image and the image to be matched as {R _m } and {S _n } respectively; (3) Perform feature matching on the SIFT feature point sets {R _m } and {S _n } extracted from the reference image and the image to be matched to obtain a rough matching feature pair {(pq) _i }, the reference image and the image to be matched The coarse matching point sets of are respectively {p _i } and {q _i }, i represents the index number of the matching pair; (4) Use the random consistency estimation algorithm to remove the mismatched pairs from the rough matching results. After removing the mismatched pairs, the set of feature point pairs retained in the reference image and the image to be matched is denoted as {(PQ) _i }. At this time, the reference image and The remaining feature points in the graph to be matched are respectively {P _i } and {Q _i }, where i represents the index number of the matching pair; (5) Using the matching point sets {P _i } and {Q _i } of the reference image and the image to be matched, affine transformation parameters are obtained through affine transformation, and finally the two images are registered according to the affine transformation parameters. 2. a kind of fast image registration realization method based on spatial sparsity and SIFT feature extraction according to claim 1, is characterized in that: described step (1), comprises the steps: 1) Carry out two-layer wavelet decomposition on the reference image R and the image to be matched respectively, and take db1 as the wavelet base to obtain the detail component images of the reference image and the image to be matched, which are denoted as R _db1 and S _db1 respectively; 2) Calculate the variance coefficient for each pixel of the detail component map of the reference image and the detail component map of the image to be matched, and then take the current pixel point as the center and take a neighborhood of size m*m to calculate the variance, and the current pixel can be obtained The variance coefficient CV of the point; set a threshold T _cv , which has a linear relationship with the variance coefficient of the SAR image T _cv = β·CV, where β is a linear coefficient; the binarized variance coefficient matrix can be obtained from the threshold T _cv CV_mask: $\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; mask</span>}<span\; class="notranslate">\; =</span>\left\{\begin{array}{c}\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; cv</span>}<span\; class="notranslate">\; <</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; cv</span>}}\end{array}\right.$ Calculate the variance coefficient matrix for the binarization of the reference image and the image to be matched, which are respectively denoted as CV_mask _R and CV_mask _S ; 3) Thinning the variance coefficient matrix of the reference image to obtain a thinning matrix T _R , and similarly refining the variance coefficient matrix of the matching image to obtain a thinning matrix T _S ; 4) Calculate the sparse matrix of the reference image and the image to be matched; 5) Select the top 20% of the pixels with the largest sparsity value in the sparse matrix as the sparsity point, and expand the sparse point to obtain the sparse area. 3. a kind of fast image registration realization method based on spatial sparsity and SIFT feature extraction according to claim 2, is characterized in that: described step 4), comprises the steps: (a) Calculate the size and direction of the gradient of each pixel in the refinement matrix: (b) Analyze the direction of the gradient of each pixel point, count the direction characteristics of the current pixel point, and obtain the direction characteristic matrix of the reference image and the image to be matched, which are recorded as _DR and _DS respectively; (c) From the refinement matrix and the direction feature matrix, the sparsity of the reference image and the image to be matched can be calculated, and the sparsity is defined as: $\mathrm{<span\; class="notranslate">\; S</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}}\underset{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; D.</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&Center\; Dot;</span>{{\mathrm{<span\; class="notranslate">\; T</span>}}_{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$ Among them, (i, j) is the index of the current pixel point, Ω is the r×r neighborhood centered on the current pixel point, N is the total number of pixels in the neighborhood Ω, D is the direction feature matrix, and T is defined as fine matrix; thus obtain the sparse matrix S _R and S _S of the reference image and the image to be matched. 4. a kind of fast image registration realization method based on spatial sparsity and SIFT feature extraction according to claim 3, is characterized in that: described step (b), comprises the steps: a) Construct a 16-direction template with a size of (2r+1)*(2r+1), and the angle ranges corresponding to each sub-direction of the direction template are respectively i∈[1,16], r is the radius of the direction template; b) Analyzing the direction characteristics of each thinning edge point in the thinning matrix with 16 direction templates; c) Normalize the directional feature matrices of the reference image and the image to be matched to obtain the directional feature matrices D _R and D _S of the reference image and the image to be matched.