Multispectral image registration method

Abstract

The invention discloses a multispectral image registration method. The method comprises the following steps: S10, preprocessing a visible light image and an infrared image; s20, carrying out canny edge detection and phase consistency detection, and detecting feature points and main directions of the infrared image and the visible light image by adopting a feature point detection method; s30, solving a gradient image, and calculating a feature descriptor of the infrared image and a multi-scale feature descriptor of the visible light image according to the gradient image; s40, performing bilateral matching by adopting cosine similarity, and screening registration points by applying an improved RANSAC (Random Sample Consensus) algorithm; and S50, mapping the matched registration point coordinates to the original image coordinates, and performing affine transformation on the visible light image to enable the visible light image to be aligned with the infrared image in a pixel level. According to the method, on the premise that the infrared image and the visible light image have the same feature points, the screening precision of the registration points is improved through the feature description and the improved RANSAC registration point screening method.

Description

Technical Field

[0001] The invention belongs to the field of image registration, and in particular relates to a multispectral image registration method. Background Art

[0002] Generally, different sensors of a multispectral camera have different fields of view and resolution. In order to effectively combine multispectral image information and realize the power transmission equipment detection function based on multispectral images, the infrared and visible light images taken by the multispectral camera need to be further registered to ensure that the infrared and visible light images correspond one to one in structure. Multispectral images are different from single-spectral images, and the feature attributes near the same key points in multispectral images are basically different, which makes the algorithm used for single-spectral image registration generally not applicable to multispectral image registration.

[0003] Infrared images use infrared radiation sources, infrared sensors and optical systems to capture infrared radiation emitted by objects, convert it into electrical signals, and generate an image after signal processing. This image displays the temperature distribution and thermal characteristics of the object in grayscale or pseudo-color. Visible light images use the reflection or transmission of light in the visible light band, and focus these lights onto light-sensitive sensors through appropriate optical lenses or reflective devices to generate image signals. Due to the different imaging principles of infrared and visible light images, the two images have great differences in color, gradient, etc. Summary of the invention

[0004] In view of this, the present invention proposes a multispectral image registration method, which uses the contour extreme points on the two images as feature points, and uses the phase image instead of the original image to detect feature points based on the similarity between the phase features of the infrared image and the visible light image; finally, combined with the edge extreme points, it makes it possible to extract the same feature points on the two images, avoiding the situation where the traditional feature point detection algorithm detects different feature points on the infrared and visible light images.

[0005] The technical solution includes the following steps:

[0006] S10, preprocessing the visible light image and the infrared image;

[0007] S20, performing canny edge detection and phase consistency detection on the visible light image and the infrared image, and using a feature point detection method to detect feature points of the infrared image and the visible light image and their main directions;

[0008] S30, obtaining a gradient image of the infrared image and the visible light image, and calculating a feature descriptor of the infrared image and a multi-scale feature descriptor of the visible light image according to the gradient image;

[0009] S40, uses cosine similarity for bilateral matching and applies the improved RANSAC algorithm to screen the registration points;

[0010] S50, mapping the matched registration point coordinates to the original image coordinates, and performing affine transformation on the visible light image to align it with the infrared image at the pixel level.

[0011] Preferably, the preprocessing in S10 includes scaling the resolution of the visible light image to the same size as the infrared image and performing Gaussian filtering on the scaled visible light image.

[0012] Preferably, the S20 includes S21, edge detection and edge extreme point detection, which specifically includes the following steps:

[0013] S211, grayscale the original image and then perform Gaussian filtering to smooth the edges and remove noise. The convolution kernel of the Gaussian filter is 7×7 in size;

[0014] S212, performing Sobel filtering on the Gaussian filtered image to obtain gradient information of the image;

[0015] S213, applying adaptive canny edge extraction to obtain an image edge binary map, using 30% and 5% of the maximum gradient as the high threshold and low threshold of the canny operator respectively;

[0016] S214, using the cv2.findContours function in OpenCV to extract continuous contours in the edge binary image and ignoring contours whose length is less than a preset value;

[0017] S215, calculating the curvature and angle of each point on the contour, and extracting the points whose curvature and angle meet the conditions as extreme points.

[0018] Preferably, the canny edge detection in S20 specifically includes:

[0019] For a two-dimensional curve, x=f(t), y=g(t), the curvature calculation formula for any point is:

[0020]

[0021] Since the image contour is a series of discrete points, , the first-order derivative and second-order derivative in formula (1) are simplified as follows:

[0022]

[0023] Where i represents the index of the current point in the contour, the first and second derivatives of y with respect to t are the same, and the two endpoints of the contour are ignored to ensure that x (i+1) and xi-1 meaningful;

[0024] The cv2.findContours function in OpenCV is used for contour extraction. The mode=cv2.RETR_LIST in this function means detecting all contours, and the method=cv2.CHAIN_APPROX_NONE means storing all contour points. This function setting ensures that the extracted contours are continuous, and that dx / dt and dy / dt are not 0 at the same time, thereby ensuring that the K value will not have NAN.

[0025] Since the points of the contour are continuous, and the coordinate difference between two adjacent points is Δx = x i+1 -x i ∈{0,1},Δy=y i+1 -y i ∈{0,1}, Δx and Δy are not 0 at the same time, so the curvature K of any point on the contour calculated according to formulas (1) and (2) is i ∈[0,1], and take points whose curvature is greater than a certain threshold as candidate extreme points.

[0026] Preferably, the phase consistency detection in S20 specifically includes:

[0027] The phase consistency calculation of the image is first performed using the Gaussian expansion of the Log-Gabor filter to obtain the 2D-Log-Gabor filter, and its formula is:

[0028]

[0029] In the formula, s represents the scale, o represents the angle, (θ,θ) is the logarithmic polar coordinate, (ρ s ,θ so ) is the center frequency, σ ρ and σ θ represents the frequency bandwidth of ρ, θ;

[0030] Perform inverse Fourier transform on equation (14) to convert it from frequency domain to spatial domain, convolve the real and imaginary parts with the image to obtain the response values ​​in scale and angle directions, and then calculate the corresponding amplitude A so and Phase After introducing the noise compensation T, the feature consistency can be calculated:

[0031]

[0032] In the formula, w0 is the weight factor, ε is the minimum value;

[0033] Calculate multiple angles θ for multiple directions o oPhase consistency PC(θ o ), and then PC(θ o ), normalize the maximum moment M and multiply it by 256 to convert it into a maximum moment map, then perform FAST feature point detection on the maximum moment map, and find the pixel center of gravity in the local area of ​​each point, and use the direction vector between the feature point and the pixel center of gravity as the main direction of the feature point;

[0034] Finally, the edge extreme points and the maximum moment graph feature points are merged as feature points.

[0035] Preferably, the gradient image of the infrared image and the visible light image is obtained in S30, specifically comprising: firstly, performing Gaussian filtering on the original image to reduce the influence of image noise, and the filter kernel adopts a size of 5×5; then, respectively obtaining the horizontal gradient and vertical gradient of the image, and calculating the overall gradient; finally, normalizing the overall gradient to 255 to obtain the gradient image.

[0036] Preferably, in S30, it is assumed that the effective radius s of each feature point (x, y) is 36 pixels, θ is the main direction of the feature point, the effective range of the feature point is divided into NBP×NBP=4×4 regions, and the direction histogram of NBO=8 directions is calculated for each region, and finally a 4×4×8=128-dimensional vector composed of the gradient direction histogram is used as a feature descriptor of the feature point;

[0037] The Gaussian weighting coefficient is used to weight the region, and the formula is as follows:

[0038]

[0039] In the formula, (g x ,g y ) is the relative coordinate of the region, and NBP is the region size;

[0040] The gradient of the gradient image within the range of the feature point is graded and normalized into six levels. The grading method is as follows: first, the regional gradient PatchMagnitudes, PM, is normalized to 0-1, and then it is sorted from small to large according to the gradient size to obtain SortedMagnitudes, SM. SM is a one-dimensional vector with a length of (2s+1)×(2s+1)=73×73. SM is evenly divided into 5 segments, that is, in SM, 4 segmentation points are taken, and the coordinates of these 4 segmentation points on the SM axis are, i=1,2,3,4;

[0041] Record segment point I i The gradient value at , using SM i Indicates that PM(x,y) is divided into 5 levels as follows:

[0042]

[0043] Calculate the relative coordinates of the pixel points within the range of the feature points relative to the feature points, and transform the relative coordinates into the main direction coordinate system. The formula is:

[0044]

[0045] Where:

[0046] d x and d y : The displacement of each pixel relative to the feature point in the x and y directions, the value is -s~s; s is the effective radius of the feature point; θ is the main direction of the feature point;

[0047] The gradient angle within the range of the feature point is also converted to the main direction coordinate system and normalized. The formula is:

[0048]

[0049] Where:

[0050] PatchAngles: gradient angle within the range of feature points; fmod: modulus operation;

[0051] The gradient of the pixel points within the range of the feature point, that is, the gradient after hierarchical normalization is:

[0052] M(x+dx,y+dy)=PM(x+dx,y+dy)#(21)

[0053] n x 、n y 、n t Decompose into integer grid cell index bins x , bin y , bin t and the offsets rbin relative to these units x , rbin y , rbin t , providing a basis for weight allocation in subsequent descriptor calculations;

[0054]

[0055]

[0056]

[0057] rbin x =n x -(bin x +0.5)#(25)

[0058] rbin y =n y -(bin y +0.5)#(26)

[0059] rbin t =(n t -bin t )#(27)

[0060] The descriptor weight is calculated using the linear interpolation method, and then the descriptor D is obtained by summing up:

[0061]

[0062] Where:

[0063]

[0064] In descriptor calculation, and Represents the offset of the grid cell relative to the current pixel in the feature point field. They represent the offset in the x direction, y direction and direction histogram respectively, and are used to traverse the adjacent grid cells around each feature point. The descriptor of each feature point is calculated by trilinear interpolation based on the contribution of the grid where it is located and its adjacent grid cells.

[0065] Preferably, the S40 includes: using the infrared image feature descriptor vector to make a dot product with each scale visible light feature descriptor vector, respectively, to find the cosine value between each vector, to form an array, and then sorting the array from small to large, comparing the ratio of the minimum cosine value to the second smallest cosine value, when the ratio is less than a certain threshold, it can be considered that the two feature points that obtain the minimum cosine value are aligned; when a feature point in the infrared image is aligned with different feature points of the visible light image at multiple scales, compare the cosine values ​​between them, find the two smallest points, and use them as the final alignment points.

[0066] Preferably, in S40, a structural similarity index is introduced into RANSAC to replace the confidence level for determining the registration points, and three pairs of registration points with a confidence level greater than a certain value are selected. There are several groups, and in order to reduce the amount of calculation, only one group with the same confidence level is selected to participate in the calculation.

[0067] Preferably, in S50, an affine transformation is performed on each group with the same confidence selected in S40 to obtain the affine transformation parameters, and then an affine transformation is performed on the visible light edge binary image to obtain a transformed visible light edge binary image, and the edge binary image of the infrared image is obtained at the same time; then, the structural similarity index value of the transformed visible light edge binary image and the infrared edge binary image is obtained, and finally, the group of points with the largest structural similarity index value replaces the original group with the largest confidence as the benchmark group for screening.

[0068] Compared with the prior art, the method of the present invention uses three-point sampling derivatives to estimate the curvature of the contour, and then combines the contour angle to determine the extreme point of the contour curve and uses phase consistency to find the maximum moment. The maximum moment image is subjected to FAST feature point detection, and the two feature points are merged as the feature points of the image, ensuring the premise that the infrared image and the visible light image have the same feature points. At the same time, a feature description suitable for infrared and visible light images is designed, and an improved RANSAC registration point screening method is designed to improve the screening accuracy of the registration points. The feasibility of the algorithm in this paper is verified by quantitative and visual means, especially in rotation and similar structure scenes, the infrared and visible light image registration accuracy is significantly improved compared with algorithms such as CAO-C2F and improved RIFT. The comprehensive registration accuracy for infrared and visible light image registration is better than that of image registration algorithms such as CAO-C2F, and can be used for the registration of infrared images and visible light images of power transmission equipment. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] In order to make the purpose, technical solution and beneficial effects of the present invention clearer, the present invention provides the following drawings for illustration:

[0070] Figure 1 A flowchart of the steps of the multispectral image registration method according to an embodiment of the present invention;

[0071] Figure 2 It is a flow chart of feature point detection of a multispectral image registration method according to an embodiment of the present invention;

[0072] Figure 3 It is a feature descriptor calculation and registration flow chart of the multispectral image registration method according to an embodiment of the present invention;

[0073] Figure 4 This is a flow chart of edge extreme point and main direction calculation of the multispectral image registration method according to an embodiment of the present invention;

[0074] Figure 5 A diagram of true extreme points and false extreme points of the multispectral image registration method according to an embodiment of the present invention;

[0075] Figure 6 It is a contour direction map of candidate extreme points of the multispectral image registration method according to an embodiment of the present invention. DETAILED DESCRIPTION

[0076] The preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0077] See also Figure 1 The multispectral image registration method of the present invention comprises the following steps:

[0078] S10, preprocessing the visible light image and the infrared image;

[0079] S20, performing canny edge detection and phase consistency detection on the visible light image and the infrared image, and using a feature point detection method to detect feature points of the infrared image and the visible light image and their main directions;

[0080] S30, obtaining a gradient image of the infrared image and the visible light image, and calculating a feature descriptor of the infrared image and a multi-scale feature descriptor of the visible light image according to the gradient image;

[0081] S40, uses cosine similarity for bilateral matching and applies the improved RANSAC algorithm to screen the registration points;

[0082] S50, mapping the matched registration point coordinates to the original image coordinates, and performing affine transformation on the visible light image to align it with the infrared image at the pixel level.

[0083] See also Figure 2 , Figure 4 ,The preprocessing in S10 includes scaling the resolution of the visible light image to the same size as the infrared image and performing Gaussian filtering on the scaled visible light image.

[0084] S20 includes S21, edge detection and edge extreme point detection;

[0085] S22, phase consistency detection and feature point detection.

[0086] Wherein S21, edge detection and edge extreme point detection, specifically includes the following steps:

[0087] S211, grayscale the original image and then perform Gaussian filtering to smooth the edges and remove noise. The convolution kernel of the Gaussian filter is 7×7 in size;

[0088] S212, performing Sobel filtering on the Gaussian filtered image to obtain gradient information of the image;

[0089] S213, applying adaptive canny edge extraction to obtain an image edge binary map, using 30% and 5% of the maximum gradient as the high threshold and low threshold of the canny operator respectively;

[0090] S214, using the cv2.findContours function in OpenCV to extract continuous contours in the edge binary image and ignoring contours whose length is less than a preset value;

[0091] S215, calculating the curvature and angle of each point on the contour, and extracting the points whose curvature and angle meet the conditions as extreme points.

[0092] The canny edge detection in S20 specifically includes:

[0093] For a two-dimensional curve, x=f(t), y=g(t), the curvature calculation formula for any point is:

[0094]

[0095] Since the image contour is a series of discrete points and the curvature does not need to be accurately calculated in the curvature calculation of the present invention, the first-order derivative and the second-order derivative in formula (1) are simplified as follows:

[0096]

[0097] Where i represents the index of the current point in the contour, the first and second derivatives of y with respect to t are the same, and the two endpoints of the contour are ignored to ensure that x (i+1) and x i-1 meaningful;

[0098] The cv2.findContours function in OpenCV is used for contour extraction. The mode=cv2.RETR_LIST in this function means detecting all contours, and the method=cv2.CHAIN_APPROX_NONE means storing all contour points. This function setting ensures that the extracted contours are continuous, and that dx / dt and dy / dt are not 0 at the same time, thereby ensuring that the K value will not have NAN.

[0099] Since the points of the contour are continuous, and the coordinate difference between two adjacent points is Δx = x i+1 -x i ∈{0,1},Δy=y i+1 -y i ∈{0,1}, Δx and Δy are not 0 at the same time, so the curvature K of any point on the contour calculated according to formulas (1) and (2) is i ∈[0,1], and take points whose curvature is greater than a certain threshold as candidate extreme points.

[0100] Many of the candidate extreme points extracted only by curvature are false extreme points due to edge noise or image blur, so the candidate extreme points need to be further screened. This paper screens the edge angles at both ends of the candidate extreme points. The angles of the edges at both ends of the true extreme points in the local neighborhood should be small, while the angles of the edges at both ends of the false extreme points in the local neighborhood should be large. Figure 5 The shapes of the true extreme point A0 and the false extreme point B0 within one pixel are the same, so the curvatures of the two points calculated by formula (1) and formula (2) are the same. Taking point A0 as an example, the curvature result is as follows:

[0101]

[0102] The main idea of ​​calculating the contour angle is to select the contour of the neighborhood near the candidate extreme point, take the candidate extreme point as the center point, divide the contour into two parts (assuming left and right), and then use two straight lines passing through the candidate extreme point to fit the left and right contours. The fitting method can use the least squares method to fit the straight line or use the weighted average method to find the contour average point and determine the straight line from two points. Since the farther the point is from the candidate point, the greater the influence on the direction of the fitted straight line, if the weighted average is used, the farther the distance from the candidate point, the greater the weight should be, giving the farther points a greater influence. Common weighted average methods include: simple weighted average, exponential weighted average, Gaussian weighted average, etc.

[0103] Simple weighted average multiplies each data point by its corresponding weight, then adds all the products together and finally divides by the sum of the weights to calculate the average. The mathematical expression is as follows:

[0104]

[0105] Exponentially weighted average: uses an exponential function to distribute weights, usually using an exponential smoothing coefficient to control the weights of different points.

[0106] Gaussian weighted average: Using Gaussian distribution (normal distribution) weight distribution, more recent data points still have a certain weight, but the weight decreases more smoothly.

[0107] Simple weighted average requires manual weighting and is not easy to control; exponential weighted average and Gaussian weighted average have similar effects and can automatically assign weights, but exponential weighted average focuses on capturing the latest trends more sensitively, while Gaussian weighted average weight assignment is smoother and more robust to outliers, focusing on being more robust to outliers. Therefore, if the weighted average method is used to find the average point of the contour, and then the average point and the extreme point are used to find the vector (straight line), the Gaussian weighted method is used in the specific embodiment. The Gaussian function is:

[0108]

[0109] According to formula (4), the weight of the contour points in the local neighborhood of the candidate extreme point is:

[0110]

[0111] Where e is the base of the natural logarithm, L is the number of left (right) contour points in the local neighborhood of the candidate extreme point (including the candidate extreme point), and i is the index distance of the left (right) contour point from the candidate extreme point, which is 0, 1, ..., L-1.

[0112] Then, the coordinates of the weighted average point can be obtained according to formula (3).

[0113] If a straight line fit is used, the ideal straight line equation should be (taking A0 as an example):

[0114]

[0115] The error can be regarded as the sum of the differences between the y values ​​of the real point and the predicted point or the sum of the squares of the differences, that is:

[0116]

[0117] After that, the least square method or gradient descent method is used to find the slope a that minimizes error1, or the minimize function in "scipy.optimize" is used to optimize error2, which can also get the slope a that meets the conditions. After finding the fitting line, the unit direction vector of the line is found, and the direction is from point A0 to the contour.

[0118] Whether the linear fitting method or the Gaussian weighted average method is used, the purpose of finding the contour angle can be achieved. Figure 6 Gaussian weighted averaging is used to obtain the weighted average points A-(B-) and A+(B+) of the contours on both sides of the candidate extreme point. Therefore, the angle between the contours on both sides of the candidate extreme point can be obtained according to formulas (9)(10)(11), and then a threshold is set for the angle φ. Only when φ is less than the threshold will it be considered as a true extreme point. The threshold set in the specific embodiment is 135°.

[0119] Taking point A0 as an example, the calculation formula of angle φ1 is:

[0120]

[0121]

[0122] When the true extreme point is extracted, the opening direction of the contour, i.e., the main direction of the extreme point, can be obtained according to the contour vectors on both sides obtained by formula (9), and the corresponding and Unit vectors with the same direction and Then the opening direction, i.e. the main direction, is The direction is as follows:

[0123]

[0124] Directions:

[0125]

[0126] After judging the contour angle in the neighborhood of the candidate extreme point, the false extreme points can be greatly reduced. At this point, the image contour extreme points and main directions are calculated. The feature point extraction and main direction calculation process proposed in this paper are summarized in Algorithms 1.1 and 1.2.

[0127]

[0128]

[0129]

[0130]

[0131] The phase consistency detection in S20 focuses on the phase information of the image rather than the amplitude information. It is usually relatively insensitive to changes in the brightness and contrast of the image. When there are factors such as local light and dark, the phase consistency feature is relatively stable. At the same time, the phase consistency feature will consider the local structure and phase information of the image when extracting features, so the edge information of the image can be retained to a certain extent. Especially in the case of low edge contrast, the phase consistency feature may more prominently show the position and direction of the edge, which is helpful for detecting and identifying edge features in the image. Therefore, phase consistency can compensate for the influence of factors such as illumination on canny edge detection to a certain extent.

[0132] Specifically include:

[0133] The phase consistency calculation of the image is first performed using the Gaussian expansion of the Log-Gabor filter to obtain the 2D-Log-Gabor filter, and its formula is:

[0134]

[0135] In the formula, s represents the scale, o represents the angle, (ρ,θ) is the logarithmic polar coordinate, (ρ s ,θso ) is the center frequency, σ ρ and σ θ represents the frequency bandwidth of ρ, θ;

[0136] Perform inverse Fourier transform on equation (14) to convert it from frequency domain to spatial domain, convolve the real and imaginary parts with the image to obtain the response values ​​in scale and angle directions, and then calculate the corresponding amplitude A so and Phase After introducing the noise compensation T, the feature consistency can be calculated:

[0137]

[0138] In the formula, w0 is the weight factor, ε is the minimum value;

[0139] Calculate multiple angles θ for multiple directions o o Phase consistency PC(θ o ), and then PC(θ o ), normalize the maximum moment M and multiply it by 256 to convert it into a maximum moment map, then perform FAST feature point detection on the maximum moment map, and find the pixel center of gravity in the local area of ​​each point, and use the direction vector between the feature point and the pixel center of gravity as the main direction of the feature point;

[0140] Finally, the edge extreme points and the maximum moment graph feature points are merged as feature points.

[0141] In S30, feature descriptors are generally calculated on the original image. This calculation method is effective for the feature description of single-spectrum images, but for multi-spectral images, such as infrared and visible light images, due to the large difference between the two spectral images, this method of feature description will lead to excessive differences in descriptors at the same position, thus reducing the registration rate. In order to solve this problem, this paper proposes a feature descriptor calculation method based on gradient images, and proves the effectiveness and practicality of this calculation method through experiments. In the registration of infrared and visible light images, using gradient images instead of original images to calculate descriptors for SURF can enhance its registration ability.

[0142] In S30, the gradient image of the infrared image and the visible light image is obtained, which specifically includes: firstly, Gaussian filtering is performed on the original image to reduce the influence of image noise, and the filter kernel size is 5×5; then, the horizontal gradient and vertical gradient of the image are respectively obtained, and the overall gradient is calculated; finally, the overall gradient is normalized to 255 to obtain the gradient image.

[0143] Assume that the effective radius s of each feature point (x, y) is 36 pixels, θ is the main direction of the feature point, divide the effective range of the feature point into NBP×NBP=4×4 regions, calculate the direction histogram of NBO=8 directions for each region, and finally use the 4×4×8=128-dimensional vector composed of the gradient direction histogram as the feature descriptor of the feature point;

[0144] The Gaussian weighting coefficient is used to weight the region, and the formula is as follows:

[0145]

[0146] In the formula, (g x ,g y ) is the relative coordinate of the region, and NBP is the region size;

[0147] The gradient of the gradient image within the range of the feature point is graded and normalized into six levels. The grading method is as follows: first, the regional gradient PatchMagnitudes, PM, is normalized to 0-1, and then it is sorted from small to large according to the gradient size to obtain SortedMagnitudes, SM. SM is a one-dimensional vector with a length of (2s+1)×(2s+1)=73×73. SM is evenly divided into 5 segments, that is, in SM, 4 segmentation points are taken, and the coordinates of these 4 segmentation points on the SM axis are, i=1,2,3,4;

[0148] Record segment point I i The gradient value at , using SM i Indicates that PM(x,y) is divided into 5 levels as follows:

[0149]

[0150] Calculate the relative coordinates of the pixel points within the range of the feature points relative to the feature points, and transform the relative coordinates into the main direction coordinate system. The formula is:

[0151]

[0152] Where:

[0153] d x and d y : The displacement of each pixel relative to the feature point in the x and y directions, the value is -s~s; s is the effective radius of the feature point; θ is the main direction of the feature point;

[0154] The gradient angle within the range of the feature point is also converted to the main direction coordinate system and normalized. The formula is:

[0155]

[0156] Where:

[0157] PatchAngles: gradient angle within the range of feature points; fmod: modulus operation;

[0158] The gradient of the pixel points within the range of the feature point, that is, the gradient after hierarchical normalization is:

[0159] M(x+dx,y+dy)=PM(x+dx,y+dy)#(21)

[0160] n x 、n y 、n t Decompose into integer grid cell index bins x , bin y , bin t and the offsets rbin relative to these units x , rbin y , rbin t , providing a basis for weight allocation in subsequent descriptor calculations;

[0161]

[0162]

[0163]

[0164] rbin x =n x -(bin x +0.5)#(25)

[0165] rbin y =n y -(bin y +0.5)#(26)

[0166] rbin t =(n t -bin t )#(27)

[0167] The descriptor weight is calculated using the linear interpolation method, and then the descriptor D is obtained by summing up:

[0168]

[0169] Where:

[0170]

[0171] In descriptor calculation, and Represents the offset of the grid cell in the feature point field relative to the current pixel. They represent the offset in the x direction, y direction and direction histogram respectively, and are used to traverse the adjacent grid cells around each feature point; the descriptor of each feature point is calculated by trilinear interpolation based on the contribution of its grid and its adjacent grid cells. These offsets are used to calculate the relationship between the grid cell center and the direction histogram to ensure smooth transition and more accurate feature description.

[0172] After the calculated descriptor is normalized, it can be used as the feature descriptor of the feature point. This paper maps the coordinates and direction to the main direction of the feature point for calculation, so that the descriptor has good rotation invariance. The hierarchical normalization process is used to increase the similarity between feature descriptors in multispectral image registration.

[0173] S40 includes: using the infrared image feature descriptor vector to do a dot product with each scale visible light feature descriptor vector, calculating the cosine value between each vector, forming an array, and then sorting the array from small to large, comparing the ratio of the minimum cosine value to the second smallest cosine value, when the ratio is less than a certain threshold, it can be considered that the two feature points that obtain the minimum cosine value are aligned; when a feature point in the infrared image is aligned with different feature points of the visible light image at multiple scales, compare the cosine values ​​between them, find the two smallest points, and use them as the final alignment points.

[0174] After obtaining the infrared image feature descriptor and the multi-scale visible light image feature descriptor, cosine similarity is used for rough matching. The basic idea is to do dot product with the infrared image feature descriptor vector and each scale visible light feature descriptor vector, find the cosine value between each vector, form an array, and then sort the array from small to large, compare the ratio of the minimum cosine value to the second minimum cosine value, when the ratio is less than a certain threshold, it can be considered that the two feature points with the minimum cosine value are registered. When a feature point in the infrared image is registered with different feature points in the visible light image at multiple scales, compare the cosine values ​​between them, find the smallest two points, and use them as the final registration points. When point A in the infrared image is registered with point B in the visible light image of scale 1 and point C in scale 2, compare the cosine value between the descriptor vector of point A and point B in scale 1 and the cosine value between the descriptor vector of point A and point C in scale 2, and take a group of points with small cosine values ​​as the final registration points. After screening by cosine similarity, many registration points can be obtained, among which there are a large number of wrong registration points, so further screening is needed.

[0175] For the selection of coarse registration points, the commonly used method is to use the Random Adoption Consensus Algorithm (RANSAC), the core idea of ​​which is to estimate the model by randomly selecting a subset of data points (i.e. random sampling), and then to determine which data points are "inliers" (fit the model) and which are "outliers" (outliers) by calculating the fitting error between the data points and the estimated model. The algorithm iterates multiple times to find the best model and finally outputs a robust model estimate.

[0176] In the feature point registration screening process, the algorithm flow of the present invention is written according to the core idea of ​​RANSAC: randomly select 3 pairs of non-collinear and sufficiently far-distance registration points (P1, Q1), (P2, Q2), (P3, Q3) from the rough registration point set (P, Q) (P is the infrared image I1 point set, Q is the optical image I2 point set), and calculate the affine transformation matrix of these 3 pairs of data (using cv2.estimateAffinePartial2D function in OpenCV to calculate, assuming that (Q1, Q2, Q3) is used as the reference), perform affine transformation on the points in P except the three points (P1, P2, P3), and obtain the error (square of the distance) between the transformed coordinates and the reference point (the corresponding point in Q). When the error of a certain point is within the maximum allowable error, the confidence of the combination is +1, and then the above process is repeated K times ( n is the total number of elements in Q), and several possible groups are obtained. Then their confidences are compared to find the three pairs of registration point combinations with the largest confidence. Based on the affine transformation matrix of these three pairs of points, all points in Q are radially transformed, and the error between the transformed point and the reference point (P) is calculated. When the error is less than the maximum allowable error, the point can be considered to be an accurately registered point.

[0177] Traditional RANSAC requires manual setting of the maximum allowable error. If this value is large, it will lead to more mismatches in the screening effect. If it is small, it may lead to insufficient number of final registration points. For this problem, the present invention proposes an improved RANSAC algorithm: introducing the average error and maximum error of the registration points. After the first screening using RANSAC, the average error and maximum error between the screened registration points are calculated. When the maximum error is greater than 1.5 times the average error, the points after this screening are continued to be screened by RANSAC. At this time, the average error of the previous time is used instead of the maximum allowable error given by the user. Through this setting, the influence of the maximum allowable error given by the user can be greatly eliminated, and the stability of the algorithm is increased.

[0178] In order to avoid the occurrence of screening errors, the present invention further improves the improved RANSAC algorithm. Regardless of whether it is the traditional RANSAC or the improved RANSAC algorithm, the main basis for screening the alignment points is to randomly extract three pairs of alignment points as mentioned above, calculate the confidence that these three pairs of points are precise alignment points, and select the three pairs of points with the largest confidence as the basis to determine which points are precise alignment points, which will eventually lead to incorrect alignment.

[0179] In order to avoid the above situation, the present invention introduces the SSIM (Structural Similarity Index (SSIM) to measure the structural similarity between images, including brightness, contrast and structure) index in RANSAC to replace the confidence to judge the registration points, and selects three pairs of registration points with confidence greater than a certain value. There are several groups. In order to reduce the amount of calculation, for the groups with the same confidence, only one group is selected to participate in the calculation. Affine transformation is performed on each selected group to obtain the affine transformation parameters, and then the visible light edge binary image is affine transformed to obtain the transformed visible light edge binary image, and the edge binary image of the infrared image is obtained at the same time. Then the SSIM values ​​of the transformed visible light edge binary image and the infrared edge binary image are calculated, and finally the group of points with the largest SSIM value replaces the original group with the largest confidence as the benchmark group for screening. In this way, the influence of similar structures can be avoided. As for why the SSIM value is calculated for the edge binary image instead of the original image, it is because the SSIM value basically cannot evaluate the relationship between the infrared image and the visible light image, while for the edge binary image, the infrared image and the visible light image have similar structures, so SSIM can be used to evaluate. The improved RANSAC screening method proposed in the present invention is also applicable to single spectrum image registration, and for a single spectrum, the SSIM can be directly calculated for the original image, and it is unnecessary to calculate the SSIM for the edge binary image. In RANSAC, a structural similarity index is introduced to replace the confidence to judge the registration points, and three pairs of registration points with confidence greater than a certain value are picked out. There are several groups. In order to reduce the amount of calculation, for the groups with the same confidence, only one group is selected to participate in the calculation.

[0180] In addition to the above embodiments, the present invention may also have other implementation modes. Any technical solution formed by equivalent replacement or equivalent transformation is within the protection scope required by the present invention.

[0181] The present invention is described in detail above, but the specific implementation of the present invention is not limited thereto. Those skilled in the art may make various modifications or alterations without departing from the spirit and scope of the claims of the present application.

Claims

1. A multispectral image registration method, characterized in that: The following steps are involved: S10, preprocessing the visible light image and the infrared image; S20, performing canny edge detection and phase consistency detection on the visible light image and the infrared image, and using a feature point detection method to detect feature points of the infrared image and the visible light image and their main directions; S30, obtaining a gradient image of the infrared image and the visible light image, and calculating a feature descriptor of the infrared image and a multi-scale feature descriptor of the visible light image according to the gradient image; S40, uses cosine similarity for bilateral matching and applies the improved RANSAC algorithm to screen the registration points; S50, mapping the matched registration point coordinates to the original image coordinates, and performing affine transformation on the visible light image to align it with the infrared image at the pixel level.

2. The multispectral image registration method according to claim 1, characterized in that: The preprocessing in S10 includes scaling the resolution of the visible light image to the same size as the infrared image and performing Gaussian filtering on the scaled visible light image.

3. The multispectral image registration method according to claim 1, characterized in that: The S20 includes S21, edge detection and edge extreme point detection, which specifically includes the following steps: S211, grayscale the original image and then perform Gaussian filtering to smooth the edges and remove noise. The convolution kernel of the Gaussian filter is 7×7 in size; S212, performing Sobel filtering on the Gaussian filtered image to obtain gradient information of the image; S213, applying adaptive canny edge extraction to obtain an image edge binary map, using 30% and 5% of the maximum gradient as the high threshold and low threshold of the canny operator respectively; S214, using the cv2.findContours function in OpenCV to extract continuous contours in the edge binary image and ignoring contours whose length is less than a preset value; S215, calculating the curvature and angle of each point on the contour, and extracting the points whose curvature and angle meet the conditions as extreme points.

4. The multispectral image registration method according to claim 1, characterized in that: The canny edge detection in S20 specifically includes: For a two-dimensional curve, x=f(t), y=g(t), the curvature calculation formula for any point is: Since the image contour is a series of discrete points, , the first-order derivative and second-order derivative in formula (1) are simplified as follows: Where i represents the index of the current point in the contour, the first and second derivatives of y with respect to t are the same, and the two endpoints of the contour are ignored to ensure that x (i+1) and x i-1 meaningful; The cv2.findContours function in OpenCV is used for contour extraction. The mode=cv2.RETR_LIST in this function means detecting all contours, and the method=cv2.CHAIN_APPROX_NONE means storing all contour points. This function setting ensures that the extracted contours are continuous, and that dx / dt and dy / dt are not 0 at the same time, thereby ensuring that the K value will not have NAN. Since the points of the contour are continuous, and the coordinate difference between two adjacent points is Δx = x i+1 -x i ∈{0,1},Δy=y i+1 -y i ∈{0,1}, Δx and Δy are not 0 at the same time, so the curvature K of any point on the contour calculated according to formulas (1) and (2) is i ∈[0,1], and take points whose curvature is greater than a certain threshold as candidate extreme points.

5. The multispectral image registration method according to claim 1, characterized in that: The phase consistency detection in S20 specifically includes: The phase consistency calculation of the image is first performed using the Gaussian expansion of the Log-Gabor filter to obtain the 2D-Log-Gabor filter, and its formula is: In the formula, s represents the scale, o represents the angle, (θ,θ) is the logarithmic polar coordinate, (ρ s ,θ so ) is the center frequency, σ ρ and σ θ represents the frequency bandwidth of ρ, θ; Perform inverse Fourier transform on equation (14) to convert it from frequency domain to spatial domain, convolve the real and imaginary parts with the image to obtain the response values ​​in scale and angle directions, and then calculate the corresponding amplitude A so and Phase After introducing the noise compensation T, the feature consistency can be calculated: In the formula, w0 is the weight factor, ε is the minimum value; Calculate multiple angles θ for multiple directions o o Phase consistency PC(θ o ), and then PC(θ o ), normalize the maximum moment M and multiply it by 256 to convert it into a maximum moment map, then perform FAST feature point detection on the maximum moment map, and find the pixel center of gravity in the local area of ​​each point, and use the direction vector between the feature point and the pixel center of gravity as the main direction of the feature point; Finally, the edge extreme points and the maximum moment graph feature points are merged as feature points.

6. The multispectral image registration method according to claim 1, characterized in that: The step S30 of obtaining the gradient image of the infrared image and the visible light image specifically includes: firstly, performing Gaussian filtering on the original image to reduce the influence of image noise, and the filter kernel adopts a size of 5×5; then, respectively obtaining the horizontal gradient and the vertical gradient of the image, and calculating the overall gradient; finally, normalizing the overall gradient to 255 to obtain the gradient image.

7. The multispectral image registration method according to claim 1, characterized in that: In the S30, it is assumed that the effective radius s of each feature point (x, y) is 36 pixels, θ is the main direction of the feature point, the effective range of the feature point is divided into NBP×NBP=4×4 regions, and the direction histogram of NBO=8 directions is calculated for each region. Finally, a 4×4×8=128-dimensional vector composed of the gradient direction histogram is used as the feature descriptor of the feature point; The Gaussian weighting coefficient is used to weight the region, and the formula is as follows: In the formula, (g x ,g y ) is the relative coordinate of the region, and NBP is the region size; The gradient of the gradient image within the range of the feature point is graded and normalized into six levels. The grading method is as follows: first, the regional gradient PatchMagnitudes, PM, is normalized to 0-1, and then it is sorted from small to large according to the gradient size to obtain SortedMagnitudes, SM. SM is a one-dimensional vector with a length of (2s+1)×(2s+1)=73×73. SM is evenly divided into 5 segments, that is, in SM, 4 segmentation points are taken, and the coordinates of these 4 segmentation points on the SM axis are, i=1,2,3,4; Record segment point I i The gradient value at , using SM i Indicates that PM(x,y) is divided into 5 levels as follows: Calculate the relative coordinates of the pixel points within the range of the feature points relative to the feature points, and transform the relative coordinates into the main direction coordinate system. The formula is: Where: d x and d y : The displacement of each pixel relative to the feature point in the x and y directions, the value is -s~s; s is the effective radius of the feature point; θ is the main direction of the feature point; The gradient angle within the range of the feature point is also converted to the main direction coordinate system and normalized. The formula is: Where: PatchAngles: gradient angle within the range of feature points; fmod: modulus operation; The gradient of the pixel points within the range of the feature point, that is, the gradient after hierarchical normalization is: M(x+dx,y+dy)=PM(x+dx,y+dy)#(21) n x 、n y 、n t Decompose into integer grid cell index bins x , bin y , bin t and the offsets rbin relative to these units x , rbin y , rbin t , providing a basis for weight allocation in subsequent descriptor calculations; rbin x =n x -(am x +0.5)#(25) rbin y =n y -(am y +0.5)#(26) rbin t =(n t -am t )#(27) The descriptor weight is calculated using the linear interpolation method, and then the descriptor D is obtained by summing up: Where: In descriptor calculation, and Represents the offset of the grid cell relative to the current pixel in the feature point field. They represent the offset in the x direction, y direction and direction histogram respectively, and are used to traverse the adjacent grid cells around each feature point. The descriptor of each feature point is calculated by trilinear interpolation based on the contribution of the grid where it is located and its adjacent grid cells.

8. The multispectral image registration method according to claim 1, characterized in that: The S40 includes: using the infrared image feature descriptor vector to perform dot product with each scale visible light feature descriptor vector, respectively, to obtain the cosine value between each vector, to form an array, and then sorting the array from small to large, comparing the ratio of the minimum cosine value to the second smallest cosine value, when the ratio is less than a certain threshold, it can be considered that the two feature points that obtain the minimum cosine value are aligned; when a feature point in the infrared image is aligned with different feature points of the visible light image at multiple scales, compare the cosine values ​​between them, find the two smallest points, and use them as the final alignment points.

9. The multispectral image registration method according to claim 1, characterized in that: In the S40, a structural similarity index is introduced into RANSAC to replace the confidence level for determining the registration points, and three pairs of registration points with a confidence level greater than a certain value are selected. There are several groups. To reduce the amount of calculation, only one group with the same confidence level is selected to participate in the calculation.

10. The multispectral image registration method according to claim 9, characterized in that: In the S50, an affine transformation is performed on each group with the same confidence selected in S40 to obtain the affine transformation parameters, and then an affine transformation is performed on the visible light edge binary image to obtain a transformed visible light edge binary image, and the edge binary image of the infrared image is obtained at the same time; then the structural similarity index value of the transformed visible light edge binary image and the infrared edge binary image is obtained, and finally the group of points with the largest structural similarity index value is used to replace the original group with the largest confidence as the benchmark group for screening.

Images