High-robustness image stitching method

Abstract

A high-robustness image stitching method disclosed by the present invention includes: dividing a preprocessed sequence of image sets into a reference image and an image to be stitched; constructing a non-linear scale space for feature detection; using a multi-scale depth hybrid feature descriptor HMD to describe feature points; obtaining an initial matching set and a candidate matching set; using a guided maximum likelihood random consensus algorithm to purify the candidate matching set and obtaining optimal model parameters through iterative calculation; eliminating the matching set that does not conform to the best parameter model; aligning the images through perspective transformation and finding the position of the best stitching line; determining the best stitching line and obtaining a seamless stitched panoramic image based on multi-band weighted fusion; outputting the final result after cropping the edges of the panoramic image. The present invention integrates the multi-scale and depth information of the target scene, effectively improves the stitching quality of panoramic images and the accuracy of geometric correction, suppresses the artifact defects caused by repeated textures and occlusions, and has higher stitching accuracy and robustness.

Description

Technical Field

[0001] The present invention relates to the technical field of pattern recognition and image processing, and more specifically, to a high-robustness image stitching method. Background Art

[0002] With the development of digital image processing technology, image stitching algorithms have been widely applied in various fields. Image stitching aims to seamlessly synthesize multiple images into a panoramic image, and is widely used in fields such as photography, map making, virtual reality, medical imaging, and surveillance systems. However, there are many challenges in the stitching process, including illumination changes, perspective distortion, feature matching in the overlapping area, and geometric distortion between different images.

[0003] In current panoramic image stitching technology, image registration and image fusion are two key technologies. Image registration refers to transforming two or more images with overlapping areas into the same coordinate system through unified transformation operations to ensure the spatial consistency between the images. Image fusion aims to achieve a smooth transition of the images to be stitched in the overlapping area, eliminate redundant pixel information in the overlapping part of the images, and thus generate a high-resolution and high-quality panoramic image. With the increasingly complex application scenarios of image stitching technology, there is an urgent need for a panoramic image stitching method with a wide adaptation range, high stitching accuracy, and strong robustness.

[0004] In the document with the application number "202410623962.6", the title is "A Panoramic Image Stitching Method and System for a Panoramic Camera". This method first uses the ORB algorithm for feature point extraction, then performs image registration based on the bidirectional KNN algorithm. After obtaining the accurate matching pairs, it calculates the homography matrix of the two views for stitching and fusion, and performs smoothing processing on the stitched image to generate the final panoramic image. However, this method has several problems: First, in the feature detection stage, the ORB algorithm used relies on the feature point direction information to achieve rotation invariance. However, at large rotation angles, the stability of the feature descriptor may decrease, resulting in a reduction in the matching accuracy. In addition, the ORB feature descriptor lacks the depth information of the target scene, thereby reducing the geometric correction accuracy of the panoramic image stitching. The descriptor is highly sensitive to noise and illumination changes and may perform poorly under complex backgrounds or low illumination conditions. Second, the bidirectional KNN algorithm mainly performs matching based on local information and does not consider the global structure and relationships of the entire image. Even if the local matching results seem good, the overall registration accuracy may still be insufficient. Relying solely on the bidirectional KNN algorithm for accurate matching without using other data purification methods to optimize the matching results may introduce false matches, which will in turn affect the subsequent registration process and ultimately lead to obvious distortion, misalignment, or poor stitching of the stitched image. Finally, the progressive fade-in and fade-out fusion method based on linear weights generates pixels using linear weights in the overlapping area. This method may cause obvious gradual edges visually, especially in the case of large illumination and color changes, and may result in an unnatural transition in the fusion area. In addition, this method may also cause blurring at the edge details, especially in areas with rich edge information, making the edge details of the image smoothed and thus losing important structural information. Summary of the Invention

[0005] To solve at least one of the above technical problems, the present invention proposes a highly robust image stitching method to achieve

[0006] The first aspect of the present invention provides a highly robust image stitching method, which is characterized by including the following steps:

[0007] S1, dividing the preprocessed sequence image set into a reference image and an image to be stitched;

[0008] S2, using the AKAZE-CHARBONNIER algorithm to construct a non-linear scale space for feature detection;

[0009] S3, using the multi-scale depth hybrid feature descriptor HMD to describe the feature points;

[0010] S4. Use BF matching to obtain an initial matching set, and use the KNN algorithm to eliminate the matching sets where the Euclidean distances of the nearest neighbor and the second nearest neighbor are greater than γ to obtain a candidate matching set;

[0011] S5. Based on the guiding mechanism, use the guided maximum likelihood random sample consensus algorithm to purify the candidate matching set, and obtain the optimal model parameters through iterative calculation;

[0012] S6. Eliminate the matching sets that do not conform to the best parameter model to obtain the best matching set;

[0013] S7. Align the images through perspective transformation and find the position of the best stitching line;

[0014] S8. Determine the best stitching line, and obtain a seamless stitched panoramic image based on multi-band weight fusion;

[0015] S9. Output the final result after cropping the edges of the panoramic image.

[0016] In a preferred embodiment of the present invention, in step S2, the AKAZE-CHARBONNIER algorithm is used to construct a non-linear scale space for feature detection, which specifically includes:

[0017] The non-linear partial differential equation for image non-linear diffusion is described as:

[0018]

[0019] Among them, I(x, y, t) represents the value of the image at the position (x, y) and time t, ▽I represents the gradient of the image, and C(||▽I||) represents a diffusion coefficient related to the gradient magnitude. The formula is as follows:

[0020]

[0021]

[0022] Among them, the function G(▽|I(x, y, t)|) represents the diffusion coefficient related to the image gradient ▽I, which is a non-linear function, depending on the gradient magnitude of the image at the position (x, y) and time t, and controls the diffusion behavior of the image in different regions. The diffusion in the regions with larger gradients at the image edges is smaller, while in the flat regions of the image with smaller gradients, the diffusion is larger, so as to achieve effective denoising and edge preservation. The parameter K is a contrast coefficient that determines the diffusion level. Introduce the CHARBONNIER non-linear diffusion model parameters of the formula into the diffusion coefficient to approximately solve the non-linear partial differential equation;

[0023] Among them, AKAZE-CHARBONNIER generates images of different scales through a non-linear scale space method. By continuously performing convolution operations on the image, features of different scales are obtained. The scale space formula adopted by AKAZE-CHARBONNIER is as follows:

[0024]

[0025] In the formula, sx and sy respectively represent the scale changes in the x and y directions.

[0026] In a preferred embodiment of the present invention, in step S3, a multi-scale depth hybrid feature descriptor HMD is used to describe the feature points, and the formula is as follows:

[0027]

[0028] Among them, D M-LDB (P) represents the M-LDB local binary feature descriptor, D(p) is the depth value corresponding to the feature point p(x, y), μ D , D max , D min respectively represent the statistical features of the mean, variance, maximum value, and minimum value within a local window of size k*k.

[0029] In a preferred embodiment of the present invention, in step S4, BF matching is used to obtain an initial matching set, and the KNN algorithm is used to eliminate the matches where the Euclidean distance between the nearest neighbor and the second nearest neighbor is greater than γ to obtain a candidate matching set. Ratio represents the Euclidean distance between any feature descriptor in the image to be stitched and the nearest neighbor and the second nearest neighbor feature descriptors in the reference image, and the calculation formula is as follows:

[0030]

[0031] Among them, D1 and D2 respectively represent the nearest neighbor and the second nearest neighbor feature descriptors in the reference image, and Di is used to describe the components of any feature descriptor in the image to be stitched.

[0032] In a preferred embodiment of the present invention, in step S5, a guiding mechanism is designed to optimize the estimation of the model parameters. The guiding confidence of the matching points is defined through the guiding function g(di), and the calculation formula is as follows:

[0033]

[0034] Among them, g(di) is the guiding confidence at point d i , and its value range is between 0 and 1. The parameter α controls the steepness of the guiding function and is used to adjust the change of the confidence. μ represents the center or mean of this point.

[0035] The guided maximum likelihood random sample consensus algorithm (G-MLESAC) is used to purify the candidate matching set. A threshold T is set to eliminate the matching points with low confidence. When g(di)>T, this point is included in the model estimation. The matching point set D` after being screened by the guiding mechanism is as follows:

[0036]

[0037] G-MLESAC iteratively calculates the optimal parameter model. In each iteration, the likelihood function is recalculated and the parameter model is updated until the maximum number of iterations is reached. The optimal model parameter θ` after the (k + 1)-th iteration is as follows:

[0038]

[0039] Among them, θ is the currently given model parameter, and L(θ|D`) represents the likelihood function of the matching point set D` screened by the guiding mechanism.

[0040] In a preferred embodiment of the present invention, in step S6, the candidate matching points that do not conform to the optimal model parameter θ` are regarded as outliers and eliminated, and the candidate matching points that conform to the optimal model parameter θ` are regarded as inliers and retained to obtain the best matching set for subsequent image stitching.

[0041] In a preferred embodiment of the present invention, in step S7, four groups of corresponding points are randomly selected from the optimal matching point set to construct a linear equation system. The matrix is decomposed by singular value decomposition (SVD), and the smallest eigenvalue is selected to determine the perspective matrix H. The formula is as follows:

[0042] H = arg min||Ax - b|| where A is the matrix constructed after decomposing the linear equation system formed by randomly selecting four pairs of points from the matching points by singular value decomposition (SVD), and b is the value on the right side;

[0043] The perspective matrix H that contains the largest number of inliers is regarded as the optimal perspective transformation matrix H best :

[0044]

[0045] Perform perspective transformation on the image to be stitched so that it is aligned and stitched with the reference image:

[0046] Warped Image = H best ·Source Image

[0047] Among them, Source Image is the image to be stitched, and Warped Image is the reference image;

[0048] According to the suture path strength difference value function E(x, y), the optimal suture line position is found using the dynamic programming method:

[0049]

[0050] Among them, E geometry (x, y) is the structural difference intensity value between images, and E color (x, y)) is the color difference intensity value.

[0051] In a preferred embodiment of the present invention, for the image structural difference intensity value function E geometry (x, y), the calculation formula is as follows:

[0052]

[0053] In the formula: G Source x (x, y), G Source y (x, y) respectively represent the gradient values in the x and y directions obtained by the Sobel algorithm for the image to be spliced, and G warped x (x, y), G warped y (x, y) respectively represent the gradient values of the reference image in the x and y directions.

[0054] For the image color difference intensity value E color (x, y), the calculation formula is as follows:

[0055] E color (x, y) = I Source (x, y) - I Warped (X, y)

[0056] In the formula, I Source (x, y), I waeped (x, y) respectively represent the pixel value differences between the image to be spliced and the reference image.

[0057] In a preferred embodiment of the present invention, in step S8, to determine the optimal suture line and obtain a seamless spliced panoramic image based on multi-band weight fusion, it specifically includes:

[0058] Using the Gaussian weight function w(x, y), a weight is assigned to each pixel point (x, y) in the splicing area. The calculation formula of the Gaussian weight function w(x, y) is as follows:

[0059]

[0060] Among them, the parameter σ is the standard deviation of the Gaussian distribution, which represents the degree of dispersion of the values, affects the attenuation speed and influence range of the weights, and determines the degree of image smoothing. The larger the value of σ, the stronger the smoothing effect, and the image noise is strongly removed, but it will also be accompanied by the loss of some details of the image; on the contrary, the smaller the value of σ, the weaker the smoothing effect, the more details of the image are retained, but the denoising effect is weakened.

[0061] The weight function w(x, y) determines the contribution degree of each pixel. The final pixel value Ifinal(x, y) at the position (x, y) can be calculated by weighted average:

[0062]

[0063] Among them, I i (x, y) represents the pixel value of the i-th image at the position (x, y), and W i (x, y) represents the weight of the i-th image at the position (x, y), which is determined by the weight function w(x, y).

[0064] In a preferred embodiment of the present invention, in step S9, after the edge of the panorama is cropped, the final result is output, which specifically includes:

[0065] By performing edge analysis on the complete panorama image Ipanorama generated after fusion, the upper, lower, left, and right boundaries of the cropping area are obtained, and the calculation formulas are as follows:

[0066]

[0067] Among them, when performing panorama edge cropping, a buffer area (padding) is added to prevent the image content from being wrongly cropped. The calculation formulas for the upper, lower, left, and right boundaries of the buffer area are as follows:

[0068]

[0069] Among them, Top, Bottom, Left, and Right respectively represent the upper, lower, left, and right boundaries of the cropping area excluding the buffer area. According to the area defined by the cropping boundaries, the panorama is cropped. The final stitched panorama image Icropped output after edge cropping is expressed as:

[0070]

[0071] Among them, Ipanorama represents the complete panorama image generated after fusion in step S8, and Top_crop, Bottom_crop, Left_crop, and Right_crop respectively represent the upper, lower, left, and right boundaries of the cropping area including the buffer area.

[0072] Compared with the prior art, the beneficial technical effects obtained by the present invention are as follows:

[0073] (1) When the present invention performs feature detection in the non-linear scale space constructed by the AKZAZE algorithm, the CHARBONNIER non-linear diffusion model with complete functional convexity of the energy functional is introduced to accelerate the solution of the non-linear partial differential equation, ensuring the existence of a unique global optimal solution for the equation and the stability of the diffusion process. This makes the feature detection algorithm more robust when dealing with noise or artifacts, effectively suppressing noise and avoiding oscillations and instability in the solution space. In addition, in the diffusion equation for image denoising, the complete functional convexity enhances the smoothness of the solution, thereby effectively reducing the generation of artifacts in the image. Secondly, the present invention extracts the depth values of specific points, combines the depth information of the surrounding areas to generate a depth feature descriptor, and then combines it with the local binary feature descriptor to generate a multi-scale depth hybrid feature descriptor HMD to describe the feature points. The combination of such multi-modal features enhances the adaptability of the descriptor to complex scenes, improves the subsequent stitching accuracy. The M-LDB descriptor can effectively capture local texture information and has local scale invariance and rotation invariance, while the depth feature provides rich geometric structure information. By fusing the multi-scale and depth information of the target scene, the hybrid feature descriptor effectively reduces the false matching rate, accelerates the feature matching process, and improves the geometric consistency and visual quality of the stitching result.

[0074] (2) In the data purification stage of the candidate matching set, the present invention uses the guided maximum likelihood random sample consensus algorithm for data purification. By designing a guiding mechanism to optimize the estimation of model parameters, using a guiding function to assign confidence levels to the matching points, and eliminating the low-confidence matching points, and then performing iterative calculations to obtain the optimal model parameters. This guiding mechanism can effectively screen the candidate matching set, comprehensively consider the previous matching results and feature information, thereby reducing unnecessary matching pairs, reducing the risk of false matching, significantly reducing the amount of data to be processed in subsequent calculations, and accelerating the solution process of the model parameters. In addition, the guiding function can accurately quantify the credibility of the model parameters, select the matching pairs that conform to the data distribution, and further improve the robustness of the model. Through iterative calculations, the algorithm can gradually improve the model parameters, repeatedly verify and update the matching results, thereby further enhancing the robustness of the stitching algorithm.

[0075] (3) In the image stitching stage of the present invention, a stitching path intensity difference value function is adopted to find the optimal stitching line position, and a Gaussian weight function is introduced to assign weights to the stitching area. Through multi-band weight fusion, smooth transition at the image stitching part is achieved. Compared with traditional stitching methods, the optimal stitching line method used in the present invention can better align the edges of the overlapping area, reduce artifacts and geometric distortion, so that the stitched image is more consistent in details. In addition, multi-band weight fusion can dynamically adjust the weights according to the local features and lighting conditions of each image, realize weighted averaging of the overlapping area, and further improve the consistency of lighting and color. By optimizing the seam position and weight assignment, the present invention effectively suppresses the artifact defects caused by lighting changes or texture differences. By first locating the optimal stitching position, the complex situations to be considered in subsequent processing are reduced, thus ensuring both the stitching quality and the processing speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 is the overall flowchart of a highly robust image stitching method of the present invention;

[0077] Figure 2 is the panoramic image stitching result diagram of each method from the side view angle of the embodiment of the present invention;

[0078] Figure 3 is the panoramic image stitching result diagram of each method from the top view angle of the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0079] In order to more clearly understand the above objects, features and advantages of the present invention, the present invention will be further described in detail below with reference to the drawings and specific embodiments. It should be noted that, without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other.

[0080] In the following description, many specific details are set forth in order to fully understand the present invention. However, the present invention can also be implemented in other ways different from those described herein. Therefore, the protection scope of the present invention is not limited by the specific embodiments disclosed below.

[0081] Please refer to Figure 1As shown in the figure, the present invention first divides the preprocessed sequence images into a reference image and an image to be stitched. The AKAZE-CHARBONNIER algorithm is used to construct a non-linear scale space to obtain the hybrid feature descriptors in the two images for subsequent matching. Then, the BF matching and the KNN algorithm are combined to screen out the candidate matching set. A guiding mechanism is designed, and the guided maximum likelihood random sample consensus algorithm is used to purify the data of the candidate matching set to obtain the optimal matching set. Next, the perspective transformation between the optimal matches is calculated, and the optimal perspective matrix is solved for the image stitching process. The stitching path intensity difference value function is used to determine the position of the optimal stitching line, and multi-band weight fusion is used to obtain a seamless stitched panoramic image. Finally, after determining the edge cropping area, the final output of the panoramic image stitching result is output.

[0082] See Figure 1 , the preprocessed sequence image set is divided into a reference image and an image to be stitched. The present invention does not limit the specific category of the original image set, but it is necessary to ensure that the features of the reference image and the image to be stitched have sufficient discriminability; the AKAZE-CHARBONNIER algorithm is used for feature extraction. The proportional coefficient K value in the CHARBONNIER non-linear diffusion model can determine how much edge information in the image region can be retained. The larger the K value, the greater the diffusion intensity, and the less edge information is retained.

[0083] The multi-scale depth hybrid feature descriptor HMD is used to describe the feature points. The M-LDB feature descriptor combines the regional intensity and the gradient mean information, and uses grids of various sizes to describe the feature points.

[0084] Each size of grid can capture feature information at different scales. Larger grids are used to eliminate high-frequency noise, and smaller grids are used to capture detailed local features.

[0085] The depth feature descriptor extracts the depth value of the feature point. The combined HMD hybrid feature descriptor will include scale invariance, rotation invariance, and the depth information of the image. Calculate the Euclidean distance between the key points of the reference image and the nearest and second-nearest neighbor feature points in the image to be stitched, and filter out the matching pairs with a ratio greater than γ to obtain a candidate matching set; allocate confidence to the candidate matching points through a guiding function, remove the matching points with a confidence lower than T, and then perform iterative calculations to obtain the optimal model parameters. Consider the points that conform to the best model as inliers and those that do not as outliers. Calculate the perspective matrix with the largest number of inliers for image stitching, and find the best stitching line position by calculating the stitching path strength difference value function; use the Gaussian weight function to allocate weights to each pixel point in the stitching area and then perform image fusion to obtain a seamless stitched panoramic image. Among them, the weights allocated to the area close to the best stitching line are lower, and the weights allocated to the area far from the best stitching area are higher; first, find the minimum bounding rectangle of all edge pixels considered to be black edges or irrelevant areas, and also set a buffer area to prevent the image content from being incorrectly cropped. After edge cropping, output the final panoramic image stitching result.

[0086] A highly robust image stitching method provided by the present invention specifically includes the following steps:

[0087] S1. Divide the preprocessed sequence image set into a reference image and an image to be stitched.

[0088] S2. Use the AKAZE-CHARBONNIER algorithm to construct a non-linear scale space for feature detection.

[0089] The AKAZE-CHARBONNIER algorithm generates multi-scale images through non-linear partial differential equations, and adjusts the intensity of diffusion and the retention degree of image edge information by setting the diffusion coefficient C and the contrast coefficient K. The mathematical expression of the non-linear partial differential equation used to generate multi-scale images is shown in formula (1).

[0090]

[0091] Among them, I(x, y, t) represents the value of the image at the position (x, y) and time t, represents the gradient of the image.

[0092] represents a diffusion coefficient related to the magnitude of the image gradient, and its mathematical expression is shown in formula (2).

[0093]

[0094]

[0095] Among them, the function represents the diffusion coefficient related to the image gradient and controls the diffusion behavior of the image in different regions. The parameter K is a contrast coefficient that determines the diffusion level. By introducing the parameters of the CHARBONNIER nonlinear diffusion model in formula (3) into the diffusion coefficient of formula (2), the nonlinear partial differential equation (1) can be approximately solved, and the CHARBONNIER nonlinear diffusion model can ensure the uniqueness of the solution of the nonlinear partial differential equation (1) because its energy functional has complete functional convexity, and the complete functional convexity can be proved by finding the energy functional of the CHARBONNIER nonlinear diffusion model and its second derivative.

[0096] The energy functional of the CHARBONNIER nonlinear diffusion model and its second derivative are solved as follows:

[0097]

[0098]

[0099] As can be seen from formula (5), if the proportionality coefficient K ≥ 0, the second derivative means that the energy functional of this nonlinear diffusion model has complete functional convexity, and this property can accelerate the solution of the nonlinear partial differential equation (1) and ensure the stability of the nonlinear diffusion process.

[0100] AKAZE-CHARBONNIER generates images at different scales through the nonlinear scale space method. By continuously performing convolution operations on the image, features at different scales are obtained. The scale space formula adopted by AKAZE-CHARBONNIER is:

[0101]

[0102] where S x and S y represent the scale changes in the x and y directions respectively.

[0103] S3. The multi-scale depth hybrid feature descriptor HMD is used to describe the feature points.

[0104] First, AKAZE-CHARBONNIER locates the key points by calculating the extreme values of the Hessian matrix:

[0105]

[0106] where I xx and I yy are the results of the second derivatives of the image in the x and y directions respectively, and Ixy is the result of the mixed derivative of the image.

[0107] Centered around the key point, calculate the difference in grayscale values within its area, and perform binary encoding based on the grayscale values to form a fixed descriptor vector D, and normalize this descriptor to obtain the local binary descriptor D M-LDB (P):

[0108]

[0109]

[0110] Among them, D(i) is the calculation formula for the difference in grayscale values, I(p j ) is the grayscale value of the j-th point, and I(p0) is the grayscale value of the feature point. Through differential comparison, the M-LDB descriptor can effectively resist the influence of illumination changes and noise.

[0111] Secondly, define the depth map D as a two-dimensional matrix, where each element D(i,j) represents the depth value of a certain pixel in the scene. If you want to obtain the local depth feature descriptor in this scene, you need to perform corner detection and calculate the response function R(x,y) of each pixel:

[0112] R(x, y) = det(M) - k·trace(M) 2 (10)

[0113] Among them, M is the autocorrelation matrix of the local window, and k is an empirical constant with a value of 0.04 - 0.06. The autocorrelation matrix M can be defined as:

[0114]

[0115] Among them, I x and I y are the gradients of the depth map D in the x and y directions respectively.

[0116] For each detected feature point p, extract the depth value of this point and combine the depth information of the surrounding area to generate a depth feature descriptor, where the depth value D(p) corresponding to the feature point p(x,y) = D(x,y). If considering the depth information around the feature point, a local window W(p) of size k*k can be defined:

[0117]

[0118] Calculate the statistical features within the window, including: mean μ D , variance σ D 2 , maximum value D max and minimum value D min :

[0119]

[0120] Combine the above statistical features into a deep feature descriptor D depth (p):

[0121]

[0122] Finally, combine the local binary descriptor DM-LDB(P) with the depth feature descriptor Ddepth(p) to generate a multi-scale deep hybrid feature descriptor HMD (Hybrid Multi-Scale Deep Descriptor):

[0123]

[0124] The HMD hybrid feature descriptor has scale invariance, rotation invariance, and the depth information of the image, and is more suitable for the panoramic image stitching task in complex scenes.

[0125] S4. Use BF matching to obtain an initial matching set and use the KNN algorithm to eliminate the matches where the Euclidean distance between the nearest neighbor and the second nearest neighbor is greater than γ to obtain a candidate matching set.

[0126] BF matching finds the optimal match between feature points by traversing all possible matching pairs. For each feature point, BF matching calculates the distances to all other feature points and selects the feature point with the smallest distance as the matching object. It can find the global optimal match but is only applicable when the number of feature points is small. BF matching combined with the KNN algorithm can filter out most of the wrong matches. The KNN algorithm calculates the distance ratio between the nearest neighbor and the second nearest neighbor for each query feature point. If the distance of the nearest neighbor is significantly smaller than the distance of the second nearest neighbor, it indicates that the match is reliable. The calculation of eliminating the matches where the Euclidean distance between the nearest neighbor and the second nearest neighbor is greater than γ to obtain a candidate matching set in the present invention is shown in formula (15).

[0127]

[0128] Among them, D1 and D2 respectively represent the feature descriptors of the nearest neighbor and the second nearest neighbor in the reference image, and Di is used to describe the components of any feature descriptor in the image to be stitched.

[0129] S5. Design a guiding mechanism, use the guided maximum likelihood random sample consensus algorithm to purify the candidate matching set, and obtain the optimal model parameters through iterative calculation.

[0130] If the homography matrix H between the reference image and the image to be stitched is calculated, then for the given model parameter θ, define a likelihood function L(θ|D) to describe the probability of observing the data D under the model parameter θ:

[0131]

[0132] Among them, P(d i |θ) represents the probability calculated for each matching point d i , according to the model parameter θ. The likelihood function simplifies the calculation through logarithmic transformation to obtain the optimal model parameter θ`:

[0133]

[0134] The present invention designs a guiding mechanism to optimize model parameter estimation, and defines the confidence of matching points through the guiding function g(d i ):

[0135]

[0136] Among them, g(d i ) is the guiding confidence at point di, and its value range is between 0 and 1. The parameter α controls the steepness of the guiding function and is used to adjust the change of confidence, and μ represents the center or mean of this point.

[0137] By setting the threshold T, the matching points with low confidence are removed. Only when g(d i ) > T, will this point be included in the model estimation. The set of matching points D` after being screened by the guiding mechanism is:

[0138] D` = {d i | g(d i ) > T} (19)

[0139] The parameter T is used to determine which matching points are credible and include them in the model estimation. Regarding the value principle of the parameter T, when the stitching scenario is a simple scenario, that is, the feature points are obvious and the matching is clear, the value of T is usually higher (0.7 or 0.8), and only the extremely reliable matching points are retained; when the stitching scenario is a complex scenario, that is, the feature points are blurred and the matching is difficult, the value of T is usually lower (0.3 or 0.4) to ensure that more feature points are retained and participate in the subsequent iterative optimization. In short, choosing the appropriate threshold T is the key step to improve the robustness and accuracy of the model. G - MLESAC calculates the optimal parameter model through iteration. In each iteration, the likelihood function is recalculated and the parameter model is updated until the maximum number of iterations is reached. The optimal model parameter θ` after the (k + 1)-th iteration is:

[0140]

[0141] Among them, θ is the currently given model parameter, and L(θ|D`) represents the likelihood function of the set of matching points D` screened by the guiding mechanism.

[0142] By introducing a guiding mechanism, maximum likelihood estimation, and step-by-step iterative optimization, G-MLESAC can more effectively eliminate mismatched points in feature matching, improving the accuracy and robustness of model estimation. Compared with traditional RANSAC, when evaluating the rationality of samples, G-MLESAC takes into account the estimation error of model parameters and the probability distribution of data points, making it more robust in dealing with noise and outliers. G-MLESAC uses the previous fitting results to guide the sample selection in the next round, which can reduce invalid sample combinations and thus improve the efficiency of data purification of the algorithm.

[0143] S6. Eliminate the matches that do not conform to the best parameter model to obtain the best matching set. The matches in the candidate matching set that do not conform to the optimal model parameter θ` are regarded as "outliers" and eliminated, and the matches that conform to the optimal model parameter θ` are regarded as "inliers" and retained to obtain the best matching set for subsequent image stitching.

[0144] S7. Align the images through perspective transformation and find the position of the best stitching line.

[0145] First, perspective transformation is used to map one image to the coordinate system of another image to align the feature points in the images for stitching. In the present invention, 4 pairs of points are randomly selected from the matching points to construct a linear equation system, the matrix is decomposed by singular value decomposition (SVD), and the smallest eigenvalue is selected to determine the perspective matrix H:

[0146]

[0147] H = H = arg min||Ax - b|| (22)

[0148] where h ij is the matrix element of the perspective transformation, which determines the image rotation, translation, scaling, and perspective effects, (x, y) represents the point coordinates in the image to be stitched, (x`, y`) represents the point coordinates in the reference image, A is the constructed matrix, and b is the value on the right side.

[0149] In each iteration, the perspective matrix H is recalculated. The points in the matching points that conform to the optimal parameter model θ` are regarded as inliers, and those that do not conform are regarded as outliers. The perspective matrix H with the largest number of inliers is regarded as the optimal perspective transformation matrix H best :

[0150]

[0151] Use the calculated optimal perspective transformation matrix Hbest to perform perspective transformation on the image to be stitched (Source Image) to align it with the reference image (Warped Image):

[0152] Warped Image = H best ·Source Image (24)

[0153] The purpose of finding the optimal suture line is to find the suture path with the smallest strength difference value in the splicing area. The calculation formula of the suture path strength difference value function E(x, y) is:

[0154]

[0155] The calculation formula of the image structure difference strength value function Egeometry(x, y) is as follows:

[0156]

[0157] In the formula: G Source x (x, y), G Source y (x, y) represent the gradient values in the x and y directions obtained by the Sobel algorithm for the image to be spliced, respectively. G warped x (x, y), G warped y (x, y) represent the gradient values in the x and y directions of the reference image, respectively.

[0158] The calculation formula of the image color difference strength value Ecolor(x, y) is as follows:

[0159] E color (x, y) = I Source (x, y) - I Warped (x, y) (27)

[0160] In the formula, ISource(x, y) and Iwaeped(x, y) represent the pixel value differences between the image to be spliced and the reference image, respectively.

[0161] Among them, Egeometry(x, y) is the image structure difference strength value, and Ecolor(x, y) is the image color difference strength value. GSourcex(x, y) and GSourcey(x, y) represent the gradient values in the x and y directions obtained by the Sobel algorithm for the image to be spliced, respectively. Similarly, Gwarpedx(x, y) and Gwarpedy(x, y) represent the gradient values in the x and y directions of the reference image. ISource(x, y) and Iwaeped(x, y) represent the pixel value differences between the image to be spliced and the reference image. The gradient values Gx and Gy of the Soble operator in the x and y directions:

[0162]

[0163] Among them, I is the pixel value of the image. Calculate the intensity difference values of all points on each stitching path, and select the path with the smallest intensity difference as the optimal stitching line.

[0164] S8. After determining the optimal stitching line, perform multi-band weight fusion to obtain a seamless stitched panoramic image.

[0165] After confirming the optimal stitching line, it is necessary to perform multi-band weight fusion (Multi-band blending) on the images to be stitched and the reference image. Use different frequency bands to fuse the images, assign a weight to each pixel in the images, with a lower weight in the area close to the optimal stitching line and a higher weight in the area far from the optimal stitching line. During the stitching process, there is a common part between the two images, and this part needs to be weighted and fused to achieve a smooth transition. The present invention uses the Gaussian weight function w(x, y) to assign a weight to each pixel point (x, y) in the stitching area:

[0166]

[0167] Among them, σ controls the width of the Gaussian function and determines the degree of smoothness.

[0168] The final fused pixel is the result of weighted averaging. The weight function W(x, y) determines the contribution degree of each pixel. The final pixel value I final (x, y) at the position (x, y) can be calculated by weighted averaging:

[0169]

[0170] Among them, I i (x, y) represents the pixel value of the i-th image at the position (x, y), and wi(x, y) represents the weight of the i-th image at the position (x, y), which is determined by the weight function w(x, y). In formula (30), the numerator is the sum of the weighted values of all images at the current position, and the denominator is the sum of all weights, ensuring weight normalization, making the stitching in the overlapping area smoother and reducing the abruptness at the boundary.

[0171] S9. After cropping the edges of the panoramic image, output the final result.

[0172] The panoramic image generated after fusion often has a black edge part, and it is necessary to detect and identify the edge part for cropping. If the complete panoramic image generated after fusion is defined as Ipanorama, the cropping area is obtained through the edge analysis of Ipanorama, and the minimum circumscribed rectangle of all edge pixels considered to be black edges or irrelevant areas is found, that is, the upper, lower, left, and right boundaries of the cropping area. The boundaries of the cropping area can be expressed as:

[0173]

[0174] A buffer area (padding) also needs to be added to prevent the image content from being wrongly cropped. The boundaries of the buffer area are as follows for top, bottom, left, and right:

[0175]

[0176] According to the area defined by the cropping boundaries, perform cropping operations on the panoramic image. The final stitched panoramic image Icropped output after edge cropping can be expressed as:

[0177]

[0178] The following uses specific application examples to illustrate the stitching effects of a highly robust image stitching method provided by the present invention for panoramic images from side view and top view perspectives: The operating conditions of the embodiments: 64-bit Windows 10 operating system, Intel(R) Core(TM) i5-9400F processor, GPU acceleration library of CUDA 10.1 version, environment: Visual Studio 2015, OpenCV 4.4.0, OpenCV_contrib 4.4.0, cmake-1.17.2-win64-x64, etc.

[0179] The publicly available dataset used in the embodiments is the "UC Merced Land Use Dataset". Considering factors such as geometric deformation, illumination differences, feature matching difficulties, and background complexity that may be caused by different perspective changes, two types of images from side view and top view perspectives are respectively selected for stitching experiments, and a comparative experiment is set up. In the comparative method, the mature SIFT algorithm is used for feature extraction, and the BF algorithm and KNN algorithm are combined to obtain a candidate matching point set. At the same time, the classic Random Sample Consensus algorithm (RANSAC) is used for data purification, and the best perspective transformation matrix is obtained according to the optimal model parameters calculated by RANSAC. In the overlapping area, a linear weight-based fade-in and fade-out fusion method is applied to generate pixels, and the final panoramic stitched image is output after edge cropping. Figure 2 and Figure 3 show the panoramic stitching results of the comparative method and the method of the present invention from side view and top view perspectives. From the Figure 2 and Figure 3 local enlarged views, it is observed that under different perspectives, the method of the present invention significantly reduces artifacts, color differences, and stitching misalignment phenomena. This benefits from the feature descriptors provided by the AKAZE-CHARBONNIER algorithm, which contains multi-scale and depth information of the current scene.

[0180] The geometric accuracy of image stitching is improved. Compared with the RANSAC method, the present invention uses G-MLESAC for data purification. Based on the principle of maximum likelihood estimation, G-MLESAC can make more effective use of the information of all data points, is more accurate in parameter estimation, and has a stronger ability to eliminate mismatches, which will significantly improve the solution accuracy of the subsequent perspective transformation matrix. In the comparison method, although the progressive fade-in and fade-out fusion method based on linear weights is simple and effective in processing overlapping regions, Figure 2 and Figure 3 in this method, obvious gradient edges appear, making the transition of the fusion region not natural enough. However, the optimal stitching line and multi-band weight fusion method adopted by the present invention significantly improves the transition effect while maintaining details, alleviates the deficiencies brought by the linear weight method, thereby making the stitching region visually flatter and the transition of the fusion region more natural. These improvements make the stitching effect more excellent and enhance the overall quality and visual experience of the image.

[0181] The device embodiments described above are merely illustrative. For example, the division of the units is only a logical function division, and there may be other division methods in actual implementation. For example, multiple units or components can be combined, or can be integrated into another system, or some features can be ignored, or not executed. In addition, the coupling, direct coupling, or communication connection between the components shown or discussed with each other can be through some interfaces, and the indirect coupling or communication connection of the device or unit can be electrical, mechanical, or other forms.

[0182] The units described above as separate components may or may not be physically separated, and the components shown as units may or may not be physical units; they can be located in one place or distributed to multiple network units; some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0183] In addition, each functional unit in the embodiments of the present invention can be all integrated in a processing unit, or each unit can be separately used as a unit, or two or more units can be integrated in a unit; the above integrated unit can be implemented in the form of hardware, or in the form of a combination of hardware and software functional units.

[0184] The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A high-robustness image stitching method, characterized in that, It includes the following steps: S1. Divide the preprocessed sequence image set into a reference image and an image to be stitched; S2. Use the AKAZE-CHARBONNIER algorithm to construct a non-linear scale space for feature detection; S3. Use the multi-scale depth hybrid descriptor HMD to describe feature points; S4. Use BF matching to obtain an initial matching set, and use the KNN algorithm to eliminate the matching sets where the Euclidean distances between the nearest neighbor and the second nearest neighbor are greater than γ to obtain a candidate matching set; S5. Based on the guiding mechanism, use the guided maximum likelihood random sample consensus algorithm to purify the candidate matching set, and obtain the optimal model parameters through iterative calculation; S6. Eliminate the matching sets that do not conform to the optimal parameter model to obtain the best matching set; S7. Align the images through perspective transformation and find the position of the best stitching line; S8. Determine the best stitching line, and obtain a seamless stitched panoramic image based on multi-band weight fusion; S9. Output the final result after cropping the edges of the panoramic image; In step S2, the AKAZE-CHARBONNIER algorithm is used to construct a non-linear scale space for feature detection, which specifically includes: The non-linear partial differential equation for image non-linear diffusion is described as: where \(I(x, y, t)\) represents the value of the image at position \((x, y)\) and time \(t\), represents the gradient of the image, represents a diffusion coefficient related to the gradient magnitude, and the formula is as follows: Among them, the function represents the diffusion coefficient related to the image gradient . It is a non-linear function that depends on the gradient magnitude of the image at the (x, y) position and time t, and controls the diffusion behavior of the image in different regions. In regions with a large gradient at the image edges, the diffusion is small, while in flat regions of the image with a small gradient, the diffusion is large, thus achieving effective denoising and edge preservation. The parameter K is a contrast coefficient that determines the diffusion level. By introducing the CHARBONNIER non-linear diffusion model parameters of the formula into the diffusion coefficient , the non-linear partial differential equation is approximately solved; Among them, AKAZE-CHARBONNIER generates images of different scales through the non-linear scale space method by continuously performing convolution operations on the image to obtain features of different scales. The scale space formula adopted by AKAZE-CHARBONNIER is as follows: In the formula, sx and sy respectively represent the scale changes in the x and y directions.

2. The high-robustness image stitching method according to claim 1, wherein, In step S3, the multi-scale depth hybrid descriptor HMD is used to describe feature points, and the formula is as follows: Among them, D M-LDB (P) represents the M-LDB local binary feature descriptor, D(p) is the depth value corresponding to the feature point p(x, y), μ D , D max , D min respectively represent the statistical features of the mean, variance, maximum value, and minimum value defined within a local window of size k*k.

3. The high-robustness image stitching method according to claim 2, wherein, In step S4, BF matching is used to obtain an initial matching set and the KNN algorithm is used to eliminate the matches where the Euclidean distances between the nearest neighbor and the second nearest neighbor are greater than γ to obtain a candidate matching set. Ratio represents the Euclidean distance between any feature descriptor in the image to be stitched and the nearest neighbor and the second nearest neighbor feature descriptors in the reference image. The calculation formula is as follows: Among them, D1 and D2 respectively represent the nearest neighbor and the second nearest neighbor feature descriptors in the reference image, and Di is used to describe the components of any feature descriptor in the image to be stitched.

4. The high-robustness image stitching method according to claim 3, characterized in that In step S5, a guiding mechanism is designed to optimize the estimation of model parameters. The guiding confidence of the matching points is defined through the guiding function g(di), and the calculation formula is as follows: where g(di) is the guiding confidence at point d i with a value range between 0 and 1, the parameter α controls the steepness of the guiding function and is used to adjust the change in confidence, and μ represents the center or mean of this point; Use the guided maximum likelihood random sample consensus algorithm (G-MLESAC) to purify the candidate matching set, set a threshold T to eliminate the low-confidence matching points. When g(di)>T, include this point in the model estimation. The matching point set D` after being screened by the guiding mechanism is: D` = {d i | g(d i ) > T} The G-MLESAC iteratively calculates the optimal parameter model. In each iteration, the likelihood function is recalculated and the parameter model is updated until the maximum number of iterations is reached. The optimal model parameter θ` after the (k + 1)-th iteration is: Among them, θ is the currently given model parameter, and L(θ|D`) represents the likelihood function of the matching point set D` screened by the guiding mechanism.

5. The high-robustness image stitching method according to claim 4, wherein In step S6, the candidates in the candidate matching set that do not conform to the optimal model parameter θ` are regarded as outliers and removed, while those that conform to the optimal model parameter θ` are regarded as inliers and retained to obtain the best matching set for subsequent image stitching.

6. The high-robustness image stitching method according to claim 5, wherein In step S7, four groups of corresponding points are randomly selected from the optimal matching point set to construct a system of linear equations, and the system of linear equations is decomposed by singular value decomposition (SVD). The perspective matrix H is determined by selecting the smallest eigenvalue, and the formula is as follows: H = arg min||Ax - b|| where A is the matrix constructed after decomposing the system of linear equations by randomly selecting four pairs of points from the matching points through singular value decomposition (SVD), and b is the value on the right side; Select the perspective matrix H with the largest number of interior points as the optimal perspective transformation matrix H best : Perform a perspective transformation on the image to be stitched to align and stitch it with the reference image: Warped Image=H best ·Source Image where Source Image is the image to be stitched and Warped Image is the reference image; According to the stitching path intensity difference value function E(x, y), use dynamic programming to find the position of the best stitching line: Among them, E geometry (x, y) is the intensity value of the structural difference between images, and E color (x, y)) is the intensity value of the color difference.

7. The high-robustness image stitching method according to claim 6, wherein Image structure difference intensity value function E geometry (x, y) is calculated as follows: Where: G Source x (x, y), G Source y (x, y) represent the gradient values in the x and y directions obtained by the Sobel algorithm for the image to be spliced, respectively, and G warped x (x, y), G warped y (x, y) represent the gradient values in the x and y directions of the reference image, respectively; Image color difference intensity value E color (x, y) is calculated as follows: E color (x,y) = I Source (x,y) - I Warped (x,y) Where, I Source (x, y) and I waeped (x, y) respectively represent the pixel value differences between the image to be spliced and the reference image.

8. The high-robustness image stitching method according to claim 7, wherein In step S8, determine the best stitching line and use multi-band weighted fusion to obtain a seamless stitched panoramic image, specifically including: Use the Gaussian weight function w(x, y) to assign a weight to each pixel point (x, y) in the stitching area. The calculation formula of the Gaussian weight function w(x, y) is as follows: where the parameter σ is the standard deviation of the Gaussian distribution, representing the degree of dispersion of the values, affecting the attenuation speed and influence range of the weights, and determining the degree of image smoothing. The larger the value of σ, the stronger the smoothing effect, and the image noise is strongly removed, but some details of the image will also be lost; on the contrary, the smaller the value of σ, the weaker the smoothing effect, and more details of the image are retained, but the denoising effect is weakened; The weight function w(x, y) determines the contribution degree of each pixel. The final pixel value Ifinal(x, y) at the position (x, y) can be calculated by weighted average: Among them, I i (x,y) represents the pixel value of the i-th image at the position (x,y), and W i (x,y) represents the weight of the i-th image at the position (x,y), which is determined by the weight function w(x,y).

9. The high-robustness image stitching method according to claim 8, characterized in that In step S9, the final result is output after cropping the edges of the panoramic image, specifically including: By analyzing the edges of the complete panoramic image Ipanorama generated after fusion, the upper, lower, left, and right boundaries of the cropping area are obtained, and the calculation formula is as follows: where when cropping the edges of the panoramic image, a buffer area (padding) is added to prevent the image content from being wrongly cropped. The calculation formula for the upper, lower, left, and right boundaries of the buffer area is as follows: where Top, Bottom, Left, and Right respectively represent the upper, lower, left, and right boundaries of the cropping area excluding the buffer area. According to the area defined by the cropping boundaries, perform a cropping operation on the panoramic image. The finally stitched panoramic image Icropped output after edge cropping is expressed as: where Ipanorama represents the complete panoramic image generated after fusion in step S8, and Top_crop, Bottom_crop, Left_crop, and Right_crop respectively represent the upper, lower, left, and right boundaries of the cropping area including the buffer area.

Images