A multi-modal image registration method, apparatus and storage device

Abstract

The application provides a multi-modal image registration method, device and storage device, relates to the technical field of remote sensing image registration, and uses a Gaussian filter to perform double-component decomposition on an original image, decomposes the original image into a structure layer and an energy layer, constructs a scale space with different scale factors to realize scale invariance, uses a FAST algorithm to extract feature points and uses a HOG descriptor to describe the feature points, uses a gray centroid method to calculate the directions of the feature points to realize direction invariance, uses a non-maximum suppression method to improve the distribution of the feature points, and uses an optimal transport algorithm instead of nearest neighbor matching, so that more and more accurate matching point pairs are obtained from coarse matching to fine matching. The application further provides a multi-modal image registration device and a storage device for implementing the multi-modal image registration method.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image registration technology, and in particular to a multimodal image registration method, device, and storage device. Background Technology

[0002] Currently, multimodal image registration algorithms can be broadly categorized into three types: region-based, feature-based, and deep learning-based. Region-based methods typically rely on image intensity information, such as NCC and MI. NCC struggles to handle nonlinear radiometric differences, severely impacting its performance. While MI is robust to nonlinear radiometric differences, it is prone to getting trapped in local optima due to ignoring the influence of surrounding pixels. Furthermore, region-based methods are sensitive to rotational and scale changes, requiring preprocessing to roughly eliminate geometric distortions before registration, limiting their versatility in many aspects. Deep learning-based methods have seen significant development recently, achieving excellent experimental results on some datasets. However, they remain difficult to apply in practical applications because, firstly, training relies on large-scale general-purpose data; secondly, deep neural networks require substantial computational resources and are inefficient, limiting the development of deep learning methods. Moreover, deep learning-based methods are only applicable to single modal problems, such as optical and SAR satellite imagery.

[0003] Feature-based methods rely on image features such as edges and corners. One representative method is the scale-invariant feature transform (SIFT), which is robust to rotation and scale changes and is therefore widely used in image registration. However, it is sensitive to nonlinear radiometric differences and struggles with multimodal image registration. The PSO-SIFT algorithm optimizes the calculation of image gradients, increasing robustness to intensity differences. The absolute phase-consistent gradient histogram (HAPCG) algorithm and the radiometrically insensitive feature transform (RIFT) algorithm use frequency domain information instead of spatial domain information to avoid the influence of nonlinear radiometric differences, but require significant computational and memory overhead. The locally normalized image feature transform (LNIFT) algorithm uses a locally normalized filter to filter the original image, greatly reducing nonlinear radiometric differences between original images and exhibiting robustness to translation and rotation changes, but it is not suitable for image pairs with complex geometric distortions. Summary of the Invention

[0004] To address the issue of large simultaneous differences in radiometric and structural information in multimodal images, and to improve the registration accuracy of multimodal images with significant nonlinear radiometric differences, this invention provides a multimodal image registration method, device, and storage device. The multimodal image registration method mainly includes:

[0005] S1: Use a Gaussian filter to perform bicomponent decomposition on the original image to obtain the structural layer and the energy layer;

[0006] S2: Construct the scale space of the structural layer and energy layer using different scale factors;

[0007] S3: Perform LPB low-pass filtering on the structural layer scale space, use the directional FAST operator to extract feature points, use the non-maximum suppression algorithm to improve the feature point distribution, and calculate the feature point orientation according to the gray centroid method to achieve rotation invariance.

[0008] S4: Describe the feature points processed in step S3 at the energy layer, calculate the feature point descriptors, use the optimal transfer method to perform feature point matching, obtain the optimal transfer coarse matching, use FSC to remove outliers, obtain the fitted transformation matrix H', add the distance metric based on the matrix H' calculation results to obtain a new cost function, and use the new cost function to find the optimal transfer scheme for the source distribution and target distribution to obtain the optimal transfer fine matching, that is, obtain the final multimodal image registration result.

[0009] A storage device that stores instructions and data for implementing a multimodal image registration method.

[0010] A multimodal image registration device includes: a processor and a storage device; the processor loads and executes instructions and data in the storage device to implement a multimodal image registration method.

[0011] The beneficial effects of the technical solution provided by this invention are as follows: Firstly, this invention uses a Gaussian filter to perform bicomponent decomposition on the image, dividing it into an energy layer and a structure layer, and constructing a scale space to achieve scale invariance. The energy layer contains rich radiation intensity information, while the structure layer captures sufficient structural detail information. In the structure layer, directional FAST is used to extract feature points, a non-maximum suppression strategy is used to improve the distribution of feature points, and the gray-scale centroid method is used to calculate the orientation of feature points to achieve rotation invariance. In the energy layer, the extracted feature points are described using a similar HOG descriptor, and blocks are extracted at intervals centered on the feature points to expand the range of image blocks to be described. When calculating the image gradient, [0°, 360°) is normalized to [0°, 180°) to alleviate gradient problems, and feature point descriptors are calculated. In feature point matching, optimal transfer is used instead of traditional nearest neighbor matching, moving from coarse matching to fine matching to obtain more and more accurate feature point matching pairs. Ultimately, this improves the registration accuracy of multimodal images with large nonlinear radiation differences. Attached Figure Description

[0012] The present invention will be further described below with reference to the accompanying drawings and embodiments. In the accompanying drawings:

[0013] Figure 1This is a flowchart of a multimodal image registration method according to an embodiment of the present invention.

[0014] Figure 2 This is a schematic diagram of the two-component scale decomposition in an embodiment of the present invention.

[0015] Figure 3 This is a schematic diagram of the low-pass Butterworth filter before and after in an embodiment of the present invention.

[0016] Figure 4 This is a visual image of the HOG class descriptor in an embodiment of the present invention.

[0017] Figure 5 This is a schematic diagram of the hardware device working in an embodiment of the present invention. Detailed Implementation

[0018] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0019] like Figure 1 As shown, a multimodal image registration method mainly includes:

[0020] S1: Use a Gaussian filter to perform bicomponent decomposition on the original image to obtain the structural layer and the energy layer;

[0021] S2: Construct the scale space of the structural layer and energy layer using different scale factors;

[0022] S3: Perform LPB low-pass filtering on the structural layer scale space, use the Feature from Accelerated SegmentTest (FAST) operator to extract feature points, use the non-maximum suppression algorithm to improve the feature point distribution, and calculate the feature point orientation according to the gray centroid method to achieve rotation invariance.

[0023] S4: Describe the feature points processed in step S3 at the energy layer, calculate the feature point descriptors, use the optimal transfer method to perform feature point matching, obtain the optimal transfer coarse matching, use the Fast Sample Consensus (FSC) algorithm to remove outliers, obtain the fitted transformation matrix H', add a distance metric based on the result of matrix H' to obtain a new cost function, and use the new cost function to find the optimal transfer scheme for the source distribution and target distribution to obtain the optimal transfer fine matching, that is, the final multimodal image registration result.

[0024] This embodiment mainly explains the process from three parts: bicomponent decomposition scale space, feature point extraction and description, and feature point matching.

[0025] 1. Bicomponent decomposition and scale space construction

[0026] Two-component decomposition technology has been widely used in image processing. It uses averaging or low-pass filtering for basic detail separation, and the resulting base layer is typically smoothed using filtering techniques. The detail layer, which complements the base layer, contains rich high-pass information. Intensity and gradient structure information are crucial components of an image, and the Gaussian function also has wide applications in image processing.

[0027] Gaussian filtering smooths an image by applying a Gaussian function kernel. This blurs the image, removes some high-frequency information, reduces image details, and preserves global information. Therefore, a Gaussian filter is used to perform a two-component decomposition on the image, yielding energy layer I. E (x,y) and structural layer I S (x,y):

[0028] Energy Layer I E The formula for calculating (x, y) is as follows:

[0029]

[0030] Where (x,y) represents the coordinates of the feature point, σ represents the standard deviation, and (x0,y0) represents the coordinates of the center point of the image.

[0031] The high-frequency components of the image contain details from the original image I, namely subtle features such as edges and textures, thus enhancing or highlighting image details to a certain extent. (Structural layer I) S The formula for calculating (x, y) is as follows:

[0032] I S (x,y)=I(x,y)-I E (x,y)

[0033] Among them, I E (x,y) represents the energy layer, and I(x,y) represents the gray value of the coordinate point (x,y) in the original image.

[0034] A scale space is constructed similarly to a Gaussian pyramid to achieve scale invariance. The original image is downsampled using k as a scaling factor, constructing a scale space with L layers, resulting in an energy layer scale space and a structure layer scale space. The energy layer scale space contains intensity information from multiple scales, while the structure layer scale space contains detail information from multiple scales. The scale definitions for each layer are as follows:

[0035] σ n =k n-1 σ0, (n=1,2...,L)

[0036] Where, σ n Let σn represent the scale of the nth layer image in scale space, σ0 represent the scale of the first layer image in scale space, and L represent the number of layers in scale space.

[0037] Figure 2 This is a schematic diagram of the original image (left) and the structure map (right). An optical image was selected. The original image contains both intensity and gradient information. After being filtered by a Gaussian filter, the structure map obtained based on the filtered high-frequency information contains rich details. Subtle features such as edges and textures are enhanced to a certain extent, making it easier to extract the required feature point information.

[0038] 2. Feature point extraction and feature point description

[0039] 2.1 Feature Point Extraction

[0040] Feature points are extracted based on the structural layer scale space. To address the sensitivity of traditional image gradients to nonlinear radiation differences, a low-pass Butterworth (LPB) filter is introduced as a filtering method to generate new gradients and reduce the impact of nonlinear radiation differences.

[0041]

[0042] LPB represents the convolution kernel of the low-pass Butterworth filter, and D(u,v) represents the distance from the feature point to the reference point (u,v). (u,v) is a binary vector specifying the size of the filter to be constructed (here set to [3,3]); cut off The cutoff frequency of the filter is represented, ranging from 0 to 0.5; n represents the filter order. Low-pass Butterworth filtering effectively reduces grayscale differences between structural layers, better highlighting detailed features. It has a significant effect on edge and corner areas of the image, extracting more and more accurate feature point information.

[0043] Based on the transformed structure diagram, feature points are extracted using the FAST algorithm. While FAST is an efficient feature point extraction algorithm, it is prone to clustering. Therefore, a non-maximum suppression algorithm is used to suppress clustering features, thus obtaining uniformly distributed feature points. The non-maximum suppression algorithm requires obtaining the intensity of the FAST feature points. Let the radius of the non-maximum suppression algorithm be r. Then, the region of non-maximum suppression is (2r+1)×(2r+1).

[0044] I'(x,y)=max(I(xr:x+r,yr:y+r))

[0045]

[0046] I(x,y) represents the FAST intensity of coordinate point (x,y) in image region B; I'(x,y) represents the result of taking the maximum value within the specified region; logical(x,y) indicates whether it is a feature point after non-maximum suppression. If it is 1, it means that coordinate point (x,y) is a feature point, and if it is 0, it means that it is not a feature point.

[0047] The principal orientation of feature points is calculated using the gray-scale centroid method. First, an image region B is obtained with the feature points as the center and radius R. Then, the moments m of image region B are calculated. pq :

[0048]

[0049] (x,y) represents the coordinates of the feature point, B represents the image region to be calculated, and the method for calculating the centroid C0 is as follows:

[0050]

[0051] Where, m 00 The zeroth-order spatial moment of the image, m 01 and m 10 Represents the moments of an image region about the x-axis and y-axis;

[0052] The principal direction of the feature point can be obtained by finding the centroid coordinates of the feature point pointing to the image region B. When describing the feature point, the image block is rotated to the principal direction to achieve direction invariance.

[0053] Figure 3 This diagram illustrates the use of the same FAST algorithm to extract feature points from both the structural image (left) and the structural image after low-pass Butterworth filtering (right). It can be observed that while the structural image contains rich details, some intensity information still remains. Using a low-pass Butterworth filter effectively reduces grayscale differences in the structural layers, better highlighting detailed features. It has a significant effect on edge and corner areas of the image, extracting more and more accurate feature point information.

[0054] 2.2 Feature Point Description

[0055] When describing feature points, use, for example Figure 4The histogram of oriented gradients (HOG) descriptor shown uses a local image patch of a certain size centered on the feature point to be described, and constructs the descriptor for the feature point using the feature information of this patch. A local image patch is divided into 10×10 grids centered on the feature point. To expand the image patch description range without increasing the descriptor dimension, the grids are extracted at intervals as valid grids. Since multimodal images often have opposite gradient directions, and the intensity usually reverses, the gradient direction is normalized from [0°, 360°) to [0°, 180°) when calculating the gradient direction, unlike traditional gradient-based descriptors. The HOG descriptor encodes each grid as an 8-bin histogram; therefore, this algorithm encodes each valid grid as a 4-bin histogram. Assuming the size of a local image patch is J×J pixels, the patch is divided into a 10×10 grid. A 4-bin histogram is constructed for each grid. The histograms of all 100 grids are concatenated and normalized to a unit vector to obtain invariance to illumination changes. Furthermore, 60 grids are extracted at intervals to construct the final descriptor describing the feature points. Therefore, a HOG-like descriptor has a length of 60×4=240 dimensions.

[0056] 3. Feature point matching

[0057] Optimal transport theory provides a method for inferring the correspondence between two distributions. In recent years, it has received considerable attention in various computer vision tasks. Optimal transport aims to solve the problem of finding the correspondence between two distributions in a given source distribution μ. s and target distribution μ t The problem of minimizing transmission cost between two images is addressed by using optimal transmission instead of the traditional nearest neighbor matching method for feature point matching. Optimal transmission describes the best way to transmit data between two different distributions, thus finding the best match between feature points of the image to be registered and the reference image can be achieved through optimal transmission.

[0058] Due to the feature point source distribution μ s and target distribution μ t μ is a discrete distribution rather than a continuous distribution. s and μ t This can be expressed as:

[0059]

[0060] Where δ(·) is the Dirac function, n s and n t p represents the number of samples in the source distribution and the number of samples in the target distribution. i s and p i tIt represents the probability quality of the i-th source distribution sample and the i-th target distribution sample, and

[0061] Since the probability mass of each location in the feature point distribution is unknown and difficult to obtain, it is assumed that the probability mass of each location in both the source and target distributions follows a uniform distribution, i.e., p. i s and p i t The probability masses are respectively and

[0062] The source distribution μ was calculated using the Kantorovich formula. s and target distribution μ t Optimal transmission result P between * :

[0063]

[0064]

[0065] Among them, P * For the optimal transfer matrix, C ij P represents the matching cost between the i-th feature point in the image to be registered and the j-th feature point in the reference image. ij This represents the matching weight between the i-th feature point in the image to be registered and the j-th feature point in the reference image. Indicates that n s line n t A matrix consisting of columns of positive real numbers, n s and n t These are the number of samples in the source distribution and the number of samples in the target distribution, 1 n If P is an n-dimensional vector of 1s, the sum of its rows is equal to μ. s The distribution vector P has a column sum equal to μ. t Distribution vector.

[0066] In the coarse matching stage, only the cost function is calculated using a similarity metric, resulting in the coarse matching result. However, coarse matching neglects the global location of feature points. Correctly matched feature point pairs need to satisfy the same transformation matrix. Therefore, a distance metric is added as part of the cost in the fine matching stage. The matching points obtained in the coarse matching stage contain a large number of outliers that do not correspond to each other. These outliers need to be removed. Fast Sample Consensus (FSC) is used to remove small outliers, with an outlier removal threshold set to 3 pixels, resulting in a fitted transformation matrix H'. However, this H' still has some error compared to the true transformation matrix H. Therefore, in the fine matching stage, a distance metric is added, which is obtained by analyzing the feature points x in the original image. i Ideally, the coordinate transformation obtained by fitting the transformation matrix H' should correspond to the coordinates y of the feature point in the target image. j Since they overlap, the difference between them yields the distance metric D:

[0067]

[0068] The distance metric D and similarity metric M are normalized to between 0 and 1, respectively:

[0069]

[0070] Cost D Cost represents the normalized distance metric. M Min represents the normalized similarity measure. D Max represents the minimum value in the distance metric. D min represents the maximum value in the distance metric. M Max represents the minimum value in the similarity metric. M This represents the maximum value in the similarity metric.

[0071] Then, the two are weighted and summed to obtain a new cost function C:

[0072] C = α × Cost D +(1-α)×Cost M ,α∈(0,1)

[0073] Here, α represents the assigned weight, with a value between 0 and 1.

[0074] Please see Figure 5 , Figure 5 This is a schematic diagram of the hardware device in operation according to an embodiment of the present invention. The hardware device specifically includes: a multimodal image registration device 401, a processor 402, and a storage device 403.

[0075] A multimodal image registration device 401: The multimodal image registration device 401 implements the multimodal image registration method.

[0076] Processor 402: The processor 402 loads and executes the instructions and data in the storage device 403 to implement the multimodal image registration method.

[0077] Storage device 403: The storage device 403 stores instructions and data; the storage device 403 is used to implement the multimodal image registration method.

[0078] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A multimodal image registration method, characterized in that: include: S1: Use a Gaussian filter to perform bicomponent decomposition on the original image to obtain the structural layer and the energy layer; S2: Construct the scale space of the structural layer and the energy layer using different scale factors; S3: Perform LPB low-pass filtering on the structural layer scale space, extract feature points using the directional FAST operator, improve the feature point distribution using the non-maximum suppression algorithm, and calculate the feature point orientation based on the gray centroid method to achieve rotation invariance. S4: Describe the feature points processed in step S3 at the energy layer, calculate the feature point descriptors, use the optimal transfer method to perform feature point matching, obtain the optimal transfer coarse matching, use FSC to remove outliers, obtain the fitted transformation matrix H', add the distance metric based on the matrix H' calculation result to obtain a new cost function, and use the new cost function to find the optimal transfer scheme for the source distribution and target distribution to obtain the optimal transfer fine matching, that is, obtain the final multimodal image registration result; In step S4, feature point descriptors are first used as similarity metrics, and the Euclidean distance between descriptors is denoted as the cost function for optimal transmission. Assuming that the extracted feature points in the image pair are uniformly distributed, the optimal transmission scheme is solved based on the source and target distributions to obtain the feature point matching results. This is the coarse matching process. In order to consider not only similarity metrics but also distance metrics between feature points in the cost function and obtain more accurate matching results, fine matching is required. Outliers are removed from the feature point matching results obtained from coarse matching using FSC, resulting in a fitted transformation matrix H'. The feature points in the source image are transformed according to matrix H', and the distances between them and the corresponding feature points in the target image are calculated. The distance and similarity metrics are normalized and then weighted and summed to obtain a new cost function. Fine matching is performed again to obtain more accurate feature point matching pairs, thereby improving the accuracy of image registration. The new cost function C is: in, Indicates the weights assigned. This represents the normalized distance metric. This represents the normalized similarity measure. Source distribution and target distribution The problem is defined as a linear optimization problem, and the source distribution is calculated using the Kantorovich formula. and target distribution Optimal transmission results between : in, For the optimal transfer matrix, This represents the matching cost between the i-th feature point in the image to be registered and the j-th feature point in the reference image. This represents the matching weight between the i-th feature point in the image to be registered and the j-th feature point in the reference image. Indicates by OK A matrix consisting of columns of positive real numbers. and These are the number of samples in the source distribution and the number of samples in the target distribution. It is an n-dimensional vector of 1. This represents the row sum of matrix P. This represents the column sum of matrix P.

2. The multimodal image registration method as described in claim 1, characterized in that: In step S1, the structural layer contains detailed information; structural layer I S The formula for calculating (x,y) is: Among them, I E (x,y) represents the energy level. Represents the coordinates of the points in the original image grayscale value; The energy layer contains image intensity information; energy layer I E The formula for calculating (x,y) is: Where (x,y) represents the coordinates of the feature point. Let (x0, y0) represent the standard deviation, and (x0, y0) represent the coordinates of the center point of the image.

3. The multimodal image registration method as described in claim 1, characterized in that: In step S2, the image is downsampled using k as a scaling factor to obtain images at different scales. Gaussian filters with different scale parameters are then used to filter the images, yielding the energy layer scale space and the structure layer scale space. The scale definitions for each layer are as follows: in, This represents the scale of the nth layer image in scale space. L represents the scale of the first layer of the scale space image, and L represents the number of layers in the scale space. The energy layer scale space contains intensity information of multi-scale images, while the structure layer scale space contains detail information of multi-scale images.

4. The multimodal image registration method as described in claim 1, characterized in that: In step S3, the formula for LPB low-pass filtering is: Where LPB represents the LPB convolution kernel, and D(u,v) represents the distance from the feature point to the reference point (u,v). This represents the cutoff frequency of the filter, and n represents the order of the filter.

5. The multimodal image registration method as described in claim 1, characterized in that: In step S3, feature points are detected using a FAST-9 detector for the structural layer, and the principal orientation of the feature points is calculated using the gray-scale centroid method. First, an image region B is obtained with the feature points as the center and radius R. The moments of image region B are then calculated using the formula shown below. : Where B represents the image region to be calculated, and x and y represent the x-coordinate and y-coordinate of a point in image region B, respectively. Represents the coordinates of point B in image region B. FAST intensity; The formula for calculating the centroid C0 is as follows: in, Represents the zeroth-order spatial moment of the image. and Represents the moments of an image region about the x-axis and y-axis; By pointing from the coordinates of the feature point to be calculated to the centroid coordinates of the image region B, the principal direction of the feature point can be obtained. When describing the feature point, the image block is rotated to the principal direction to achieve direction invariance.

6. The multimodal image registration method as described in claim 1, characterized in that: In step S4, the HOG descriptor is used to describe the feature points. The process involves constructing a local image patch centered on the coordinates of the feature point to be described. Assume the size of the local image patch is... Pixels divide a local image patch into multiple grids. In multimodal images, intensity often reverses, therefore the gradient direction is changed from... Return to One A 4-bin histogram is constructed for each grid, and the histograms of multiple grids are concatenated and normalized to a unit vector to obtain invariance to illumination changes. In order to obtain a large description range while controlling the dimension of the descriptor, the method of extracting grids at intervals is used to construct descriptors describing feature points.

7. A storage device, characterized in that: The storage device stores instructions and data for implementing the multimodal image registration method according to any one of claims 1 to 6.

8. A multimodal image registration device, characterized in that: include: A processor and a storage device; the processor loads and executes instructions and data in the storage device to implement the multimodal image registration method according to any one of claims 1 to 6.

Images