CN112017223A - Heterologous image registration method based on improved SIFT-Delaunay - Google Patents

Heterologous image registration method based on improved SIFT-Delaunay Download PDF

Info

Publication number
CN112017223A
CN112017223A CN202010952773.5A CN202010952773A CN112017223A CN 112017223 A CN112017223 A CN 112017223A CN 202010952773 A CN202010952773 A CN 202010952773A CN 112017223 A CN112017223 A CN 112017223A
Authority
CN
China
Prior art keywords
point
image
gradient
points
feature
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010952773.5A
Other languages
Chinese (zh)
Other versions
CN112017223B (en
Inventor
梁毅
刘倩
李聪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN202010952773.5A priority Critical patent/CN112017223B/en
Publication of CN112017223A publication Critical patent/CN112017223A/en
Application granted granted Critical
Publication of CN112017223B publication Critical patent/CN112017223B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

The invention relates to the technical field of image registration, and discloses a heterogeneous image registration method based on improved SIFT-Delaunay. And then, when extracting the feature points, the method is the same as that of the traditional SIFT algorithm, and the edge strength information with higher precision is obtained by adopting an improved exponential weighted average ratio operator for calculating the gradient of the feature points. And finally, performing coarse registration by using Euclidean distance in the registration process of the feature points, and performing fine registration of the feature points by using a method of a Delaunay triangular grid, wherein the Delaunay triangular grid has uniqueness and has good invariance to translation, rotation and geometric distortion of an image.

Description

Heterologous image registration method based on improved SIFT-Delaunay
Technical Field
The invention relates to the technical field of heterologous image registration, in particular to a heterologous image registration method based on improved SIFT-Delaunay.
Background
The image registration technology is a basic part of image processing, and the heterogeneous image registration technology can make up for the defects of a single sensor image in the image registration technology and is a key technology researched in the field of current image processing. The performance of a registration algorithm can be seriously influenced by the nonlinear gray difference between different source images, and the traditional SIFT algorithm has scale invariance, rotation invariance and illumination invariance and is widely applied to practical engineering application. Due to multiplicative noise of the SAR image, when the SAR image has a high-reflectivity region and obtains a gradient with a large amplitude, a certain deviation occurs in the traditional SIFT (scale invariant feature transform) during gradient extraction, and subsequent feature point matching is influenced. In the feature description process, the SIFT algorithm constructs a descriptor by using local information in the field range of point features, and ignores the association among feature points in an image, so that the SIFT algorithm is unstable in performance in the registration of a heterogeneous image.
In order to obtain a heterogeneous image registration result with strong robustness and high accuracy, stable feature point gradient information needs to be extracted. The main gradient operators at present are the following: roberts gradient operator, Prewitt gradient operator, and Sobel gradient operator. The Roberts gradient operator is sensitive to speckle noise generated by the SAR image, the Prewitt gradient operator and the Sobel gradient operator have certain improvement on noise resistance, but due to the smoothing principle of the Prewitt gradient operator and the Sobel gradient operator, edge information of the image is lost, so that the Roberts gradient operator and the Sobel gradient operator are not suitable for a registration algorithm of the visible light image and the SAR image.
Disclosure of Invention
Aiming at the defects of the traditional SIFT algorithm in the registration of the heterogeneous images, the invention aims to provide a method for registering the heterogeneous images based on the improved SIFT-Delaunay. The method adopts an improved ROEWA (exponential weighted average ratio) operator to calculate the gradient of the characteristic point, and edge strength information with higher precision is obtained; and performing coarse registration by adopting an Euclidean distance, performing fine registration of the feature points by adopting a Delaunay triangular mesh method, enhancing the correlation between the feature points, and enabling the image to have better translation, rotation and geometric distortion invariance.
In order to achieve the purpose, the invention is realized by adopting the following technical scheme.
The heterogeneous image registration method based on the improved SIFT-Delaunay comprises the following steps:
step 1, selecting an optical image as a reference image and an SAR image as an image to be registered; PPB filtering processing is carried out on the SAR image, speckle noise in the SAR image is removed, and the SAR image after noise reduction is obtained;
step 2, extracting candidate feature points of the reference image and the SAR image subjected to noise reduction respectively by adopting an SIFT algorithm; refining the candidate characteristic points of each image, and removing low-contrast points and edge response points to obtain two corresponding groups of characteristic points;
step 3, calculating the gradient amplitude and gradient direction of each group of feature points through an improved ROEWA operator, and determining the main direction of the feature points; rotating the corresponding image according to the main direction of each group of feature points to correspondingly obtain a rotated image; extracting a feature descriptor from each feature point in a neighborhood corresponding to the rotated image;
step 4, measuring a feature descriptor corresponding to each feature point in the reference image and the image to be registered by adopting the Euclidean distance to obtain a corresponding rough matching point pair;
and 5, constructing a Delaunay triangular grid on the coarse registration point pair to eliminate a mis-registration point pair to obtain a fine registration point pair, and completing registration of the two heterogeneous images.
Compared with the prior art, the invention has the beneficial effects that:
the method firstly utilizes an improved exponential weighted mean ratio operator to calculate the gradient amplitude and the direction corresponding to the characteristic point. Because valuable characteristic points of the image are generally distributed at the edge of the image, the operator can effectively inhibit speckle noise contained in the SAR image and ensure that the gradient module value at the edge is obviously stronger than that of points in a uniform homogeneous region; secondly, by utilizing the uniqueness of the Delaunay triangular mesh, the relative position relation between the feature points in the target area can be reflected, the association between the feature points is enhanced, and the misregistration caused by the nonlinear gray difference between the feature points of the heterogeneous image is eliminated.
Drawings
The invention is described in further detail below with reference to the figures and specific embodiments.
FIG. 1 is a schematic diagram of a Gaussian pyramid scale space according to an embodiment of the present invention;
FIG. 2 is a flow chart of an implementation of the method of the present invention;
FIG. 3 is a diagram illustrating feature point detection according to an embodiment of the present invention;
fig. 4 is a heterogeneous image pair of simulation experiment data 1. Wherein (a) is an SAR image and (b) is a visible light image;
fig. 5 is a graph of simulation results connecting lines of simulation experimental data 1 based on the conventional SIFT algorithm;
FIG. 6 is a graph of simulation results for simulation experiment data 1 based on the method of the present invention;
fig. 7 is a heterogeneous image pair of simulation experiment data 2, in which (a) is an SAR image and (b) is a visible light image;
fig. 8 is a graph of simulation results connecting lines of simulation experimental data 2 based on the conventional SIFT algorithm;
fig. 9 is a line graph of simulation results of simulation experiment data 2 based on the method of the present invention.
Detailed Description
Embodiments of the present invention will be described in detail below with reference to examples, but it will be understood by those skilled in the art that the following examples are only illustrative of the present invention and should not be construed as limiting the scope of the present invention.
In the embodiment of the invention, the construction of the Gaussian pyramid scale space based on different scale space factors is shown in FIG. 1. The purpose of constructing the scale space is to simulate the multi-scale features of the image in vision, and detect the visual features which have invariance to the scale change of the image, so that the finally extracted features have the scale scaling resistance. An image can be represented in scale space as a convolution of the image and a variable Gaussian kernel function, which is represented by using a Gaussian pyramid (LOG) operator as follows:
Figure BDA0002677577470000041
in the formula, symbol
Figure BDA0002677577470000042
Representing convolution operations between functions; i (x, y) represents an input image; g (x, y, σ) is a two-dimensional gaussian function, as shown below, with σ representing a scale factor.
Figure BDA0002677577470000043
The formation of (x, y, σ) on the image scale space is performed by convolving the image I (x, y) with a gaussian filter function G (x, y, σ) with a variable kernel, which can be represented by an image gaussian pyramid.
In the conventional SIFT operator, the main idea is to extract the stable features of the image by using an approximation function of the scale normalized laplacian function, i.e., a difference gaussian function D (x, y, σ), as shown in the following formula. It is obtained by performing difference processing on convolution results of adjacent layers, and the expression is as follows:
d (x, y, σ) ═ G (x, y, k σ) -G (x, y, σ)) × I (x, y) ═ up (x, y, k σ) -, up (x, y, σ)
Wherein k is a constant, and k is 21/sDenotes the spacing between two adjacent dimensions; s is a preset value; the function G is a gaussian function, which is functionally a scale space corresponding to the image.
The difference pyramid DoG is obtained by subtracting adjacent layers of the gaussian pyramid LoG. The Gaussian pyramid LoG has O groups, and each group has S +3 layers of images; the differential pyramid has O groups, and each group has S +2 layers of images. Multiplying the image by one time by interpolation before the pyramid is constructed, obtaining a 0 th group of images, wherein M represents the row number of the 0 th group of images, N represents the column number of the 0 th group of images, and the image group number O of the pyramid is { log }2(min (M, N)) }, and { } is an integer operation.
Referring to fig. 2, the invention provides a heterologous image registration method based on improved SIFT-Delaunay, comprising the following steps:
step 1, selecting an optical image as a reference image and an SAR image as an image to be registered; PPB filtering processing (preprocessing) is carried out on the SAR image, speckle noise in the SAR image is removed, and the SAR image after noise reduction is obtained;
specifically, the statistical patch-based weights, namely the PPB filtering, is a spatial filtering method, and the algorithm obtains a weighted average formula through maximum likelihood estimation, defines a weight calculation formula of an image block, iteratively modifies prior information step by step, and finally stops when an optimal denoising result is converged.
For a SAR image, assuming no correlation between each pixel in the image, maIs the pixel value at pixel M, defining the noise distribution present in the image as
Figure BDA0002677577470000051
Model in which naIs a pixel N in the real image1The pixel value of (b), then the filtering is equivalent to looking for the pair N1Best estimated value of
Figure BDA0002677577470000052
The similarity measurement between two pixels by the PPB is shown as the following formula:
Figure BDA0002677577470000053
wherein S is 1/(2D-1), and D is the vision of the SAR image; n is a radical ofkThe size of the similarity window is generally 7 × 7; mNKFor a similarity window pixel point set, a search window is generally 21 × 21;
the PPB filtering gives the optimal image block similarity measurement under any noise distribution through probability theory derivation, so that the edge characteristics can be effectively maintained while additive noise and multiplicative noise are removed.
Step 2, extracting candidate feature points of the reference image and the SAR image subjected to noise reduction respectively by adopting an SIFT algorithm; refining the candidate characteristic points of each image, and removing low-contrast points and edge response points to obtain two corresponding groups of characteristic points;
specifically, in the scale space, the feature point is a local extreme point in the DoG scale space, and the local extreme value includes two meanings: one is an image space extreme value, namely, the extreme point is a local extreme point of 9 points in a 3 × 3 neighborhood of the same layer of the extreme point; the second is a scale space extreme value, namely a local extreme value point in 27 points in the 3 x 3 neighborhood of the point and the corresponding point in two adjacent layers. After the extreme point detection is carried out, the positions and the sizes of the feature points are preliminarily determined. As shown in fig. 3, the middle detection point is compared with 26 points, which are 8 adjacent points of the same scale and 9 × 2 points of the upper and lower adjacent scales, to ensure that extreme points are detected in both scale space and two-dimensional image space.
The specific process of the refining treatment comprises the following steps:
specifically, since the DoG operator generates a strong edge response, further refinement of the extreme points is required to generate a stable SIFT feature descriptor. The refining treatment of the extreme points comprises two parts, namely the suppression of low-contrast points and the elimination of edge response points.
Suppression of low contrast spots is performed first:
for an extreme point, translate the scale space function D (x, y, σ) to the current extreme point, perform taylor series (up to quadratic term) expansion around the current extreme point:
Figure BDA0002677577470000061
where X ═ (X, y, σ)TFor feature point positions and scale information vectors, X0The feature points are obtained for the first time,
Figure BDA0002677577470000062
is the first differential of the function D (X, y, sigma) at point X,
Figure BDA0002677577470000063
is the second order differential of the function D (X, y, σ) at point X.
The position of the extreme point after correction is obtained by taking the derivative of the above formula D and making the derivative zero:
Figure BDA0002677577470000064
judging whether the position of the extreme point after correction meets the requirement
Figure BDA0002677577470000065
If yes, the current candidate feature point and the corrected extreme point are processed
Figure BDA0002677577470000071
Interpolation is carried out between the two points to generate a new extreme point, and the current extreme point is updated to enter the next extreme point correction process; otherwise, calculating the value of the extreme point translation scale space function D after correction:
Figure BDA0002677577470000072
for all extreme points
Figure BDA0002677577470000073
For each extreme point
Figure BDA0002677577470000074
And if so, judging that the extreme point is an unstable extreme point with low contrast and removing the unstable extreme point, otherwise, keeping the extreme point.
Then, eliminating edge response points:
calculating the principal curvature of each feature point:
define a 2 × 2 Hessian matrix H:
Figure BDA0002677577470000075
the partial derivative in the formula can be approximated by the neighborhood difference of the current sampling point.
The principal curvature of each feature point is proportional to the eigenvalue of H. In the characteristic values of H, the characteristic vector corresponding to the characteristic value with the largest numerical value represents the direction with the largest curvature; and the eigenvector corresponding to the eigenvalue with the smallest value represents the place with the smallest curvature. Let α represent the maximum eigenvalue of H and β represent the minimum eigenvalue of H, resulting from the nature of the determinant:
Trac(H)＝Dxx+Dyy＝α+β
Det(H)＝DxxDyy-(Dxy)2＝αβ
trace of H
Figure BDA0002677577470000076
For the laplacian operator, let r satisfy α ═ r β, then:
Figure BDA0002677577470000077
from Trace (H)2The relationship between/Det (H) and the proportionality coefficient r is known as (r +1)2The value of/r is minimum when the two characteristic values are equal and increases with increasing r, so to detect whether the principal curvature is below a certain threshold r, it is only necessary to detect whether:
Figure BDA0002677577470000081
in this embodiment, r is 10, and feature points that do not satisfy the above formula are removed.
In the invention, the suppression of the low-contrast point and the elimination of the edge response point in the refining process have no sequence, and any step can be carried out firstly.
Step 3, calculating the gradient amplitude and gradient direction of each group of feature points through an improved ROEWA operator, and determining the main direction of the feature points; rotating the corresponding image according to the main direction of each group of feature points to correspondingly obtain a rotated image; extracting a feature descriptor from each feature point in a neighborhood corresponding to the rotated image;
specifically, the exponential weighted mean ratio operator is an improved method of the mean Ratio Operator (ROA), and the weighted mean value of the pixel gray values in the local neighborhood of the feature point is calculated by using an exponential function, so that the closer the pixel point is to the feature point, the higher the weight is, and the condition that all the pixel points in the mean ratio operator have the same weight is avoided. Given a certain pixel point (α, b), the exponentially weighted average of the gray values in two vertically opposite regions (in the domain) at this point is calculated as follows:
Figure BDA0002677577470000082
Figure BDA0002677577470000083
where α represents a parameter of the exponential function and R represents the neighborhood region.
The improved ROEWA operator in the invention has the following ratio of the exponential weighted mean values in the opposite area in the vertical direction and the normalization result:
Figure BDA0002677577470000091
Figure BDA0002677577470000092
the ratio R in the horizontal direction can be obtained by the same method1，αAnd normalizing the ratio T1，α. The edge amplitude mag (gradient modulus) and direction ori (gradient direction) of the image can be obtained according to the following formula:
Figure BDA0002677577470000093
Figure BDA0002677577470000094
the improved ROEWA operator will be the original one due to the non-linear gray scale difference between the different source images, the smaller the angular range of the gradient is, the more accurate the result is
Figure BDA0002677577470000095
Instead, it is changed into
Figure BDA0002677577470000096
In the calculation, a neighborhood window taking a feature point as a center is sampled, the gradient of a neighborhood pixel is counted by using a histogram, the horizontal axis of the histogram represents the gradient direction, the range of the gradient direction is from 0 degree to 360 degrees, the gradient direction is divided into 36 sections, and each 10 degrees is one section; the vertical axis represents the weighted gradient modulus m1. Selecting a Gaussian function
Figure BDA0002677577470000097
Is a weighting function w1Where σ is the scale value of the extreme point.
Figure BDA0002677577470000098
Where m (x, y) represents the gradient modulus at coordinate (x, y) within the neighborhood.
And finding the direction corresponding to the maximum module value in the gradient histogram as the main direction of the characteristic point. When there are other peaks corresponding to 80% of the main peak energy, these directions are taken as the secondary directions of the feature points. It is possible to assign a plurality of principal directions to one feature point, which can improve the stability of feature matching.
The purpose of calculating the main direction is to realize the rotation invariance of the descriptor, and the adopted method is to subtract the main direction value from the gradient direction values of all pixel points in the neighborhood to be used as new gradient direction values of the pixel points. If the new gradient direction value is a negative value, adding 2 pi; and if the sum is larger than 2 pi, subtracting 2 pi. Therefore, the new gradient direction can be converted into the range of 0 to 2 pi, and the rotation invariance of the descriptor is realized.
The process of constructing the feature descriptor comprises the following steps:
specifically, first, the two-dimensional coordinate axis corresponding to the scale space is rotated to the main direction of the feature point, and rotation invariance is ensured. Then, the gradient module value and gradient direction of all pixel points in a 16 × 16 window (where the row and column where the feature point is located are not taken) with the feature point as the center are calculated.
Each cell represents a pixel in the scale space where the feature point neighborhood is located, and Gaussian weighting is adopted (the closer the pixel is to the feature point, the greater the gradient direction information contribution is). Then dividing the neighborhood into 8 × 8 sub-regions, calculating gradient direction histograms of each 8 × 8 image patch in four directions of 0 degree, 90 degrees, 180 degrees and 270 degrees, calculating a gradient accumulated value in each gradient direction, and drawing the gradient histograms, wherein the horizontal axis of the histograms represents the gradient direction, the range of the gradient direction is from 0 to 2 pi, the histograms are divided into 4 segments in total, and each 90 degree is one segment; the vertical axis represents the weighted gradient modulus m2. At this time, a Gaussian function is selected
Figure BDA0002677577470000101
Is a weighting function w2Where σ is the scale value of the feature point.
Figure BDA0002677577470000102
And finally, combining the gradient information of each direction of all the sub-regions into a 256-dimensional vector to form a feature descriptor of the feature point.
Step 4, measuring a feature descriptor corresponding to each feature point in the reference image and the image to be registered by adopting the Euclidean distance to obtain a corresponding rough matching point pair;
and measuring the characteristic descriptors corresponding to each characteristic point of the reference image and the image to be registered by adopting Euclidean distance. And calculating the ratio of the nearest Euclidean distance to the next nearest Euclidean distance, and if the ratio is smaller than a certain threshold value, considering the nearest Euclidean distance and the next nearest Euclidean distance as a matching point pair.
Specifically, assuming that the feature point set of the reference image is p and the feature point set of the image to be registered is q, the distance between each point in the feature point set p and each point in the feature point set q is calculated one by one to obtain a distance set D between the feature points. Sequencing the elements in the distance set D to obtain the nearest distance DminAnd a next nearest neighbor distance dn-min。
The SIFT algorithm distinguishes correct matching pairs from incorrect matching pairs by judging the ratio of the nearest neighbor distance to the next nearest neighbor distance:
Figure BDA0002677577470000111
thus, for a correct matching pair, its nearest neighbor distance dminIs much smaller than the next nearest neighbor distance dn-minI.e. Dis tan ceRatio < 1; the false matching pair has a nearest neighbor distance d due to the high dimension of the feature spaceminAnd a next nearest neighbor distance dn-minThe difference is not large, i.e. Dis tan ceRatio ≈ 1. A distance ratio threshold value Tresh e (0, 1) can be taken to distinguish between correctly matched pairs and incorrectly matched pairs. When the pair of feature points is a correct matching point pair, the distance ratio is concentrated in a region having a smaller value, and when the pair is a wrong matching point pair, the distance ratio is concentrated in a region having a larger value. When rejecting all matching points with distance ratios greater than 0.8, nearly 90% of the false matching points can be rejected while only losing less than 5% of the correct matching points. A threshold value of 0.8 may be chosen. Namely:
dis tan percentage > Tresh, rejecting matched pair
Dis tan ceRatio is less than or equal to Tresh, and the received matching pair is a candidate matching pair
And 5, constructing a Delaunay triangular grid on the coarse registration point pair to eliminate a mis-registration point pair to obtain a fine registration point pair, and completing registration of the two heterogeneous images.
Specifically, firstly, constructing a Delaunay triangular grid for a coarse matching point pair; and constructing a Delaunay triangulation network through a triangulation network construction function carried by Matlab, and selecting homonymous point pairs of the reference image and the SAR image to be registered, namely pairing characteristic points in the two images. And then carrying out similarity measurement on the triangulation network constructed in the reference image and the image to be registered. Assuming that Δ ABC and Δ a ' B ' C ' are two triangles which need to be determined to be similar, where point a and point a ', point B and point B ', and point C ' are corresponding point pairs respectively, the similarity of two angles ═ a (assuming that its value is a) and ═ a ' (assuming that its value is x) is:
Figure BDA0002677577470000121
wherein the content of the first and second substances,
Figure BDA0002677577470000122
ρ is α P/3, and P is generally 1/2. For a pair of triangles, the similarity of three corresponding internal angles of the pair of triangles is obtained through the formula, and then the average value of the similarity of the three internal angles of the pair of triangles is obtained as the similarity of the pair of triangles by the following formula:
I＝(Ia+Ib+Ic)/3
and finally, searching a triangular pair with the similarity larger than 0.75 in the calculated similarity matrix, and listing the triangular pair as a matching triangular pair, wherein the corresponding point pair with the same name is a registration point pair.
Simulation experiment
To demonstrate the effectiveness of the present invention, the following simulation and comparative experiments were used for further illustration.
And carrying out registration simulation analysis by utilizing the measured data of the two groups of SAR images and the visible light images, and comparing the registration simulation analysis with the simulation result of the traditional SIFT algorithm. Table 1 shows registration data specific information.
TABLE 1 registration data specific information
Figure BDA0002677577470000123
Fig. 4 shows an SAR image and a visible light image in the experimental data 1, and it can be seen that the two images have good quality and the SAR image is less interfered by noise. The results of the registration by using the conventional SIFT algorithm and the method provided by the present invention are shown in fig. 5 and fig. 6, and comparing fig. 5 and fig. 6, it can be seen that the registration result obtained by using the conventional SIFT algorithm in fig. 5 has significant errors, and the correct registration result can be obtained by using the method provided by the present invention in fig. 6.
Fig. 7 shows the SAR image and the visible light image in the experimental data 2, and it can be seen from the graph that the SAR image is greatly affected by speckle noise. The results of the registration by using the conventional SIFT algorithm and the method provided by the present invention are shown in fig. 8 and fig. 9, and it can be seen from comparing fig. 8 and fig. 9 that the registration result obtained by using the conventional SIFT algorithm in fig. 8 has significant errors, and the correct registration result can be obtained by using the method provided by the present invention in fig. 9.
The above result is because the descriptor robustness is enhanced by calculating the magnitude and direction of the gradient using the improved ROEWA operator. The extracted points are precisely registered by using the Delaunay triangular mesh, and the gray information of the characteristic points is utilized and simultaneously the geometrical structure information of point distribution in the neighborhood of the characteristic points is combined, so that the finally obtained registered point pair is more accurate.
Although the present invention has been described in detail in this specification with reference to specific embodiments and illustrative embodiments, it will be apparent to those skilled in the art that modifications and improvements can be made thereto based on the present invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.

Claims (10)

1. The method for registering the heterogeneous images based on the improved SIFT-Delaunay is characterized by comprising the following steps of:
step 1, selecting an optical image as a reference image and an SAR image as an image to be registered; PPB filtering processing is carried out on the SAR image, speckle noise in the SAR image is removed, and the SAR image after noise reduction is obtained;
step 2, extracting candidate feature points of the reference image and the SAR image subjected to noise reduction respectively by adopting an SIFT algorithm; refining the candidate characteristic points of each image, and removing low-contrast points and edge response points to obtain two corresponding groups of characteristic points;
step 3, calculating the gradient amplitude and the gradient direction of each group of feature points through an improved ROEWA operator, determining the main direction of the feature points, rotating the corresponding image according to the main direction of each group of feature points, and correspondingly obtaining the rotated image; extracting a feature descriptor from each feature point in a neighborhood corresponding to the rotated image;
step 4, measuring a feature descriptor corresponding to each feature point in the reference image and the image to be registered by adopting the Euclidean distance to obtain a corresponding rough matching point pair;
and 5, constructing a Delaunay triangular grid on the coarse registration point pair to eliminate a mis-registration point pair to obtain a fine registration point pair, and completing registration of the two heterogeneous images.
2. The method as claimed in claim 1, wherein the SIFT algorithm is used to extract candidate feature points of the reference image and the denoised SAR image, and the following operations are performed for using the reference image and the denoised SAR image as the images to be processed respectively:
2.1, constructing an image scale space, representing a reference image or a denoised SAR image in the scale space as the convolution of the image and a variable Gaussian kernel function, and adopting a Gaussian pyramid operator to represent as follows:
Figure FDA0002677577460000021
in the formula (I), the compound is shown in the specification,
Figure FDA0002677577460000022
representing a convolution operation; i (x, y) represents an image to be processed, namely a reference image or a denoised SAR image; g (x, y, sigma) is a two-dimensional Gaussian function, and sigma represents a scale factor;
Figure FDA0002677577460000023
2.2, dividing each layer of scale space of the image to be processed into a plurality of 3 multiplied by 3 neighborhoods, determining a maximum value point of each field, searching a local maximum value point in the neighborhood of each maximum value point and a corresponding point in two adjacent layers of scale spaces, traversing the whole image to be processed to obtain all local maximum value points, namely candidate feature points of the image to be processed.
3. The method for registering the heterogeneous images based on the improved SIFT-Delaunay as claimed in claim 1, wherein the refining processing is performed on the candidate feature points of each image, specifically:
first, the suppression of low contrast points is performed for each image:
for a certain candidate feature point, translating the scale space function D (x, y, sigma) to the current candidate feature point, and performing Taylor series expansion around the current candidate feature point:
Figure FDA0002677577460000024
where X ═ (X, y, σ)TFor feature point positions and scale information vectors, X0Are the candidate feature points and are the candidate feature points,
Figure FDA0002677577460000025
is the first differential of the function D (X, y, sigma) at point X,
Figure FDA0002677577460000026
a second order differential at point X for function D (X, y, σ);
obtaining the position of the extreme point after correction by deriving the above formula D and making the derivative zero:
Figure FDA0002677577460000027
judging whether the position of the extreme point after correction meets the requirement
Figure FDA0002677577460000031
If yes, the current candidate feature point and the corrected extreme point are processed
Figure FDA0002677577460000032
Interpolate between them to produceUpdating the current extreme point and entering the next extreme point correction process; otherwise, calculating the value of the extreme point translation scale space function D after correction:
Figure FDA0002677577460000033
for all extreme points
Figure FDA0002677577460000034
For each extreme point
Figure FDA0002677577460000035
If the contrast is smaller than the set threshold, judging that the extreme point is an unstable extreme point with low contrast and removing the unstable extreme point, otherwise, keeping the extreme point;
then, eliminating edge response points of each image:
calculating the principal curvature of each feature point:
define a 2 × 2 Hessian matrix H:
Figure FDA0002677577460000036
the principal curvature of each feature point is proportional to the eigenvalue of H; let α represent the maximum eigenvalue of H and β represent the minimum eigenvalue of H, resulting from the nature of the determinant:
Trac(H)＝Dxx+Dyy＝α+β
Der(H)＝DxxDyy-(Dxy)2＝αβ
wherein, trace of H
Figure FDA0002677577460000037
For the laplacian operator, let r satisfy α ═ r β, then:
Figure FDA0002677577460000038
judgment Trace (H)2/Det (H) and (r +1)2Whether/r satisfies the following relationship:
Figure FDA0002677577460000039
if so, retaining the current candidate feature points, otherwise, rejecting the current candidate feature points.
4. The method for registering the heterogeneous images based on the improved SIFT-Delaunay as claimed in claim 1, wherein the refining processing is performed on the candidate feature points of each image, specifically:
firstly, removing edge response points of each image:
calculating the principal curvature of each feature point:
define a 2 × 2 Hessian matrix H:
Figure FDA0002677577460000041
the principal curvature of each feature point is proportional to the eigenvalue of H; let α represent the maximum eigenvalue of H and β represent the minimum eigenvalue of H, resulting from the nature of the determinant:
Trac(H)＝Dxx+Dyy＝α+β
Det(H)＝DxxDyy-(Dxy)2＝αβ
wherein, trace of H
Figure FDA0002677577460000042
For the laplacian operator, let r satisfy α ═ r β, then:
Figure FDA0002677577460000043
judgment Trace (H)2/Det (H) and (r +1)2Whether/r satisfies the following relationship:
Figure FDA0002677577460000044
if so, retaining the current candidate feature points, otherwise, rejecting the current candidate feature points;
then, the suppression of low contrast points is performed for each image:
for a certain candidate feature point, translating the scale space function D (x, y, sigma) to the current candidate feature point, and performing Taylor series expansion around the current candidate feature point:
Figure FDA0002677577460000045
where X ═ (X, y, σ)TFor feature point positions and scale information vectors, X0Are the candidate feature points and are the candidate feature points,
Figure FDA0002677577460000051
is the first differential of the function D (X, y, sigma) at point X,
Figure FDA0002677577460000052
a second order differential at point X for function D (X, y, σ);
obtaining the position of the extreme point after correction by deriving the above formula D and making the derivative zero:
Figure FDA0002677577460000053
judging whether the position of the extreme point after correction meets the requirement
Figure FDA0002677577460000054
If yes, the current candidate feature point and the corrected extreme point are processed
Figure FDA0002677577460000055
Interpolation is carried out between the two points to generate a new extreme point, and the current extreme point is updated to enter the next extreme point correction process; otherwise, calculating the value of the extreme point translation scale space function D after correction:
Figure FDA0002677577460000056
for all extreme points
Figure FDA0002677577460000057
For each extreme point
Figure FDA0002677577460000058
And if so, judging that the extreme point is an unstable extreme point with low contrast and removing the unstable extreme point, otherwise, keeping the extreme point.
5. The method for the registration of the heterogeneous images based on the improved SIFT-Delaunay as claimed in claim 1, wherein the gradient magnitude and gradient direction of each group of feature points are calculated by an improved ROEWA operator, specifically:
3.1, given a certain pixel point (a, b), the exponentially weighted average values of the gray values in two vertically opposite regions at the point are respectively:
Figure FDA0002677577460000059
Figure FDA0002677577460000061
wherein, alpha represents the parameter of the exponential function, and R represents the neighborhood region;
3.2, using a modified ROEWARatio R of exponential weighted mean values of operators in opposite regions in vertical direction3，αAnd normalizing the result T3，αComprises the following steps:
Figure FDA0002677577460000062
Figure FDA0002677577460000063
the ratio R in the horizontal direction can be obtained by the same method1，αAnd normalizing the ratio T1，αAnd then calculating the gradient amplitude mag and the gradient direction ori of the pixel point:
Figure FDA0002677577460000064
Figure FDA0002677577460000065
wherein the angular range of the gradient is
Figure FDA0002677577460000066
Calculating all pixel points according to the steps to obtain gradient amplitudes and gradient directions of all the pixel points;
3.3, counting the gradient of pixel points in the neighborhood by using a histogram for a neighborhood window taking each feature point as a center, wherein the horizontal axis of the histogram represents the gradient direction, the range of the horizontal axis is from 0 degree to 360 degrees, the horizontal axis is divided into 36 sections, and each 10 degrees is one section; the vertical axis represents the weighted gradient modulus m1(ii) a Selecting a Gaussian function
Figure FDA0002677577460000067
Is a weighting function w1And sigma is the scale value of the characteristic point, and the weighted gradient module value is as follows:
Figure FDA0002677577460000068
where m (x, y) represents the gradient modulus at coordinate (x, y) within the neighborhood.
6. The method according to claim 5, wherein the principal directions of the feature points are determined, the corresponding images are rotated according to the principal directions of each group of feature points, and the rotated images are obtained correspondingly, specifically:
first, the principal direction of the feature points is determined: finding the maximum value of the weighted gradient module value in the gradient histogram, namely, taking the gradient direction corresponding to the main peak value of the histogram as the main direction of the characteristic point; when an auxiliary peak value which is equivalent to 80% of the energy of the main peak value exists, taking the gradient direction corresponding to the auxiliary peak value as the auxiliary direction of the characteristic point, and determining a plurality of main directions of the characteristic point;
then, the image is rotated based on the main direction: taking the direction value obtained by subtracting the main direction value from the gradient direction value of each pixel point in the neighborhood as a new gradient direction value of each pixel point; if the new gradient direction value is a negative value, adding 2 pi; if the new gradient direction value is larger than 2 pi, subtracting 2 pi to obtain a new gradient direction in the range of 0 to 2 pi, and finishing image rotation.
7. The improved SIFT-Delaunay-based heterogeneous image registration method according to claim 1, wherein each feature point extracts a feature descriptor in a neighborhood corresponding to the rotated image, and the specific steps are as follows:
firstly, rotating a two-dimensional coordinate axis corresponding to a scale space to a main direction of a characteristic point;
then, calculating the gradient module values and gradient directions of all pixel points in a 16 × 16 neighborhood with the feature point as the center;
dividing the neighborhood into 8-by-8 sub-regions, and calculating the gradient accumulation of each 8-by-8 image patch in four directions of 0 degree, 90 degrees, 180 degrees and 270 degreesAdding a value, drawing a gradient histogram, wherein the horizontal axis of the gradient histogram represents the gradient direction, the range of the gradient direction is from 0 to 2 pi, the gradient direction is divided into 4 sections in total, and each section is divided into 90 degrees; the vertical axis represents the weighted gradient modulus m2Selecting a Gaussian function
Figure FDA0002677577460000071
Is a weighting function w2σ is the scale value of the feature point;
Figure FDA0002677577460000081
wherein m (x, y) represents a gradient modulus at coordinate (x, y) within the neighborhood;
finally, the weighted gradient module values of all the subregions in each direction are combined into a 256-dimensional vector to form a feature descriptor of the feature point.
8. The improved SIFT-Delaunay-based heterogeneous image registration method according to claim 1, wherein the measuring of the feature descriptors corresponding to each feature point in the reference image and the image to be registered by Euclidean distance comprises the following specific steps: and calculating Euclidean distances between each feature point in the reference image and each feature point in the image to be registered, sequencing all the Euclidean distances, calculating a ratio of the nearest Euclidean distance to the next nearest Euclidean distance, judging the feature point pair corresponding to the nearest Euclidean distance as a candidate matching point pair if the ratio is smaller than a set threshold value Tresh (0, 1), and rejecting the matching point pair if the ratio is not smaller than the set threshold value Tresh (0, 1).
9. The improved SIFT-Delaunay-based heterogeneous image registration method according to claim 1, wherein the Delaunay triangular mesh is constructed on the coarse registration point pairs to remove misregistration point pairs, specifically:
firstly, constructing a Delaunay triangulation network on the coarse registration point pair through a triangulation network construction function, and pairing feature points in two images;
secondly, carrying out similarity measurement on the Delaunay triangulation network constructed in the reference image and the image to be registered;
and finally, taking the triangle pair with the similarity larger than 0.75 as a matching triangle pair, and taking the feature point pair corresponding to the matching triangle pair as a fine registration point pair.
10. The method for improving the SIFT-Delaunay-based heterogeneous image registration according to claim 9, wherein the similarity measure is specifically: setting delta ABC and delta A ' B ' C ' to respectively correspond to two triangles to be judged in the reference image and the image to be registered, wherein the point A and the point A ', the point B and the point B ', and the point C ' are respectively corresponding point pairs, and the similarity of the two angles & lt A & gt and & lt A & gt ' is as follows:
Figure FDA0002677577460000091
wherein the content of the first and second substances,
Figure FDA0002677577460000092
p is aP/3, P is 1/2, a is an angle value of & lt A, and x is an angle value of & lt A';
for a pair of triangles, the similarity of three corresponding internal angles of the pair of triangles is obtained through the formula, and then the average value of the similarity of the three internal angles of the pair of triangles is calculated to be used as the similarity of the pair of triangles:
I＝(Ia+Ib+Ic)/3
wherein, IbIs the similarity of corresponding internal angle B in the triangle pair, IcIs the similarity of the corresponding internal angle C in the triangle pair.
CN202010952773.5A 2020-09-11 2020-09-11 Heterologous image registration method based on improved SIFT-Delaunay Active CN112017223B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010952773.5A CN112017223B (en) 2020-09-11 2020-09-11 Heterologous image registration method based on improved SIFT-Delaunay

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010952773.5A CN112017223B (en) 2020-09-11 2020-09-11 Heterologous image registration method based on improved SIFT-Delaunay

Publications (2)

Publication Number Publication Date
CN112017223A true CN112017223A (en) 2020-12-01
CN112017223B CN112017223B (en) 2024-01-30

Family

ID=73522819

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010952773.5A Active CN112017223B (en) 2020-09-11 2020-09-11 Heterologous image registration method based on improved SIFT-Delaunay

Country Status (1)

Country Link
CN (1) CN112017223B (en)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112734816A (en) * 2021-01-13 2021-04-30 西安电子科技大学 Heterogeneous image registration method based on CSS-Delaunay
CN113421257A (en) * 2021-07-22 2021-09-21 凌云光技术股份有限公司 Dot matrix font text line rotation correction method and device
CN114897675A (en) * 2022-07-14 2022-08-12 南京航空航天大学 Exponential windowing method for relevance weighting in digital image correlation
CN115035326A (en) * 2022-06-09 2022-09-09 电子科技大学 Method for accurately matching radar image and optical image
CN115115678A (en) * 2022-03-24 2022-09-27 郭俊平 High-precision registration method for interferometric synthetic aperture radar (InSAR) complex images
CN116224332A (en) * 2023-01-17 2023-06-06 中国矿业大学 Radar interference phase quality estimation method for coordinating multiple indexes

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101714254A (en) * 2009-11-16 2010-05-26 哈尔滨工业大学 Registering control point extracting method combining multi-scale SIFT and area invariant moment features
CN108304883A (en) * 2018-02-12 2018-07-20 西安电子科技大学 Based on the SAR image matching process for improving SIFT
CN109409292A (en) * 2018-10-26 2019-03-01 西安电子科技大学 The heterologous image matching method extracted based on fining characteristic optimization
WO2019134327A1 (en) * 2018-01-03 2019-07-11 东北大学 Facial expression recognition feature extraction method employing edge detection and sift
CN111145228A (en) * 2019-12-23 2020-05-12 西安电子科技大学 Heterogeneous image registration method based on local contour point and shape feature fusion

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101714254A (en) * 2009-11-16 2010-05-26 哈尔滨工业大学 Registering control point extracting method combining multi-scale SIFT and area invariant moment features
WO2019134327A1 (en) * 2018-01-03 2019-07-11 东北大学 Facial expression recognition feature extraction method employing edge detection and sift
CN108304883A (en) * 2018-02-12 2018-07-20 西安电子科技大学 Based on the SAR image matching process for improving SIFT
CN109409292A (en) * 2018-10-26 2019-03-01 西安电子科技大学 The heterologous image matching method extracted based on fining characteristic optimization
CN111145228A (en) * 2019-12-23 2020-05-12 西安电子科技大学 Heterogeneous image registration method based on local contour point and shape feature fusion

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
杨雪梅;龚俊斌;王鹏;田金文;: "基于改进SIFT的SAR图像与可见光图像配准", 航天控制, no. 06 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112734816A (en) * 2021-01-13 2021-04-30 西安电子科技大学 Heterogeneous image registration method based on CSS-Delaunay
CN112734816B (en) * 2021-01-13 2023-09-05 西安电子科技大学 Heterologous image registration method based on CSS-Delaunay
CN113421257A (en) * 2021-07-22 2021-09-21 凌云光技术股份有限公司 Dot matrix font text line rotation correction method and device
CN113421257B (en) * 2021-07-22 2024-05-31 凌云光技术股份有限公司 Method and device for correcting rotation of text lines of dot matrix fonts
CN115115678A (en) * 2022-03-24 2022-09-27 郭俊平 High-precision registration method for interferometric synthetic aperture radar (InSAR) complex images
CN115035326A (en) * 2022-06-09 2022-09-09 电子科技大学 Method for accurately matching radar image and optical image
CN115035326B (en) * 2022-06-09 2024-04-12 电子科技大学 Radar image and optical image accurate matching method
CN114897675A (en) * 2022-07-14 2022-08-12 南京航空航天大学 Exponential windowing method for relevance weighting in digital image correlation
CN116224332A (en) * 2023-01-17 2023-06-06 中国矿业大学 Radar interference phase quality estimation method for coordinating multiple indexes
CN116224332B (en) * 2023-01-17 2023-10-20 中国矿业大学 Radar interference phase quality estimation method for coordinating multiple indexes

Also Published As

Publication number Publication date
CN112017223B (en) 2024-01-30

Similar Documents

Publication Publication Date Title
CN112017223A (en) Heterologous image registration method based on improved SIFT-Delaunay
CN109409292B (en) Heterogeneous image matching method based on refined feature optimization extraction
CN107301661B (en) High-resolution remote sensing image registration method based on edge point features
CN110097093B (en) Method for accurately matching heterogeneous images
CN104318548B (en) Rapid image registration implementation method based on space sparsity and SIFT feature extraction
CN110569861B (en) Image matching positioning method based on point feature and contour feature fusion
CN113256653B (en) Heterogeneous high-resolution remote sensing image registration method for high-rise ground object
CN107862708A (en) A kind of SAR and visible light image registration method
CN109978848A (en) Method based on hard exudate in multiple light courcess color constancy model inspection eye fundus image
CN110222661B (en) Feature extraction method for moving target identification and tracking
CN104899888A (en) Legemdre moment-based image subpixel edge detection method
CN111369605A (en) Infrared and visible light image registration method and system based on edge features
CN112233116A (en) Concave-convex mark visual detection method based on neighborhood decision and gray level co-occurrence matrix description
CN114331879A (en) Visible light and infrared image registration method for equalized second-order gradient histogram descriptor
Zhang et al. Saliency-driven oil tank detection based on multidimensional feature vector clustering for SAR images
CN114897705A (en) Unmanned aerial vehicle remote sensing image splicing method based on feature optimization
CN117764983A (en) Visual detection method for binocular identification of intelligent manufacturing production line
CN109767442B (en) Remote sensing image airplane target detection method based on rotation invariant features
CN112734816B (en) Heterologous image registration method based on CSS-Delaunay
CN114943744A (en) Edge detection method based on local Otsu thresholding
CN112801141B (en) Heterogeneous image matching method based on template matching and twin neural network optimization
CN113781413A (en) Electrolytic capacitor positioning method based on Hough gradient method
CN117853510A (en) Canny edge detection method based on bilateral filtering and self-adaptive threshold
CN113298725A (en) Correction method for superposition error of ship icon image
Olson Adaptive-scale filtering and feature detection using range data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Liang Yi

Inventor after: Liu Qian

Inventor after: Xiang Cong

Inventor after: Li Cong

Inventor before: Liang Yi

Inventor before: Liu Qian

Inventor before: Li Cong

GR01 Patent grant
GR01 Patent grant