High-robustness image splicing method

Abstract

The invention discloses a high-robustness image splicing method. The method comprises the following steps: dividing a preprocessed sequence image set into a reference image and an image to be spliced; constructing a nonlinear scale space for feature detection; describing feature points by adopting a multi-scale depth mixed feature descriptor HMD; obtaining an initial matching set and a candidate matching set; purifying the candidate matching set by adopting a guide type maximum likelihood random consistency algorithm, and obtaining an optimal model parameter through iterative calculation; removing a matching set which does not conform to the optimal parameter model; aligning the images through perspective transformation, and searching an optimal suture line position; determining an optimal suture line, and obtaining a seamless spliced panorama based on multi-band weight fusion; according to the method, multi-scale and depth information of the target scene is fused, the splicing quality and geometric correction accuracy of the panoramic image are effectively improved, artifact defects caused by repeated textures and shielding are inhibited, and the method has higher splicing precision 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 highly robust image splicing method. Background Art

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

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

[0004] In the document with application number "202410623962.6", the title is "A panoramic image stitching method and system for panoramic cameras". This method first uses the ORB algorithm to extract feature points, and then performs image registration based on the bidirectional KNN algorithm. After obtaining the exact matching pair, the homography matrix of the two views is calculated for stitching and fusion, and the stitched image is smoothed 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 a large rotation angle, the stability of the feature descriptor may decrease, resulting in reduced matching accuracy. In addition, the ORB feature descriptor lacks depth information of the target scene, which reduces the geometric correction accuracy of the panoramic image stitching. The descriptor is highly sensitive to noise and illumination changes and may perform poorly in complex backgrounds or low-light conditions. Secondly, the bidirectional KNN algorithm mainly matches based on local information, without considering the global structure and relationship of the entire image. Even if the local matching results seem good, the accuracy of the overall registration may still be insufficient. Relying only on the bidirectional KNN algorithm for accurate matching without using other data purification methods to optimize the matching results may introduce mismatches, which in turn affects the subsequent registration process and ultimately leads to obvious distortion, misalignment or poor stitching in the stitched image. Finally, the linear weighted progressive fusion method uses linear weights to generate pixels in the overlapping area. This method may visually cause obvious gradient edges, especially in the case of large changes in lighting and color, which may cause the transition of the fused area to appear unnatural. In addition, this method may also cause blurring at the edge details, especially in areas with rich edge information, so that the edge details of the image are smoothed, thereby losing important structural information. Summary of the invention

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

[0006] A first aspect of the present invention provides a highly robust image stitching method, characterized in that it comprises the following steps:

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

[0008] S2, uses the AKAZE-CHARBONNIER algorithm to construct a nonlinear scale space for feature detection;

[0009] S3, uses a multi-scale deep hybrid feature descriptor HMD to describe feature points;

[0010] S4, using BF matching to obtain the initial matching set, and using the KNN algorithm to eliminate the matching sets whose nearest neighbor and next nearest neighbor Euclidean distance are greater than γ to obtain the candidate matching set;

[0011] S5, based on the bootstrapping mechanism, uses the bootstrap maximum likelihood random consensus algorithm to purify the candidate matching set, and obtains the optimal model parameters through iterative calculation;

[0012] S6, eliminate the matching sets that do not meet the optimal parameter model to obtain the best matching set;

[0013] S7, aligns images through perspective transformation and finds the best stitching line position;

[0014] S8, determine the best stitching line based on multi-band weighted fusion to obtain a seamless stitched panorama;

[0015] S9, the panorama is cropped at the edges and the final result is output.

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

[0017] The nonlinear partial differential equation for nonlinear diffusion of images is described as:

[0018]

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

[0020]

[0021]

[0022] The function G(▽|I(x,y,t)|) represents the diffusion coefficient related to the image gradient ▽I. It is a nonlinear function that depends on the gradient of the image at the (x,y) position and time t. It controls the diffusion behavior of the image in different areas. The diffusion is small in the area with large gradient at the edge of the image, while in the flat area with small gradient, the diffusion is large, thereby achieving effective denoising and edge preservation. The parameter K is the contrast coefficient that determines the diffusion level. The formula The CHARBONNIER nonlinear diffusion model parameters are introduced into the diffusion coefficient In , the nonlinear partial differential equation is approximately solved;

[0023] Among them, AKAZE-CHARBONNIER generates images of different scales through nonlinear scale space methods. By continuously performing convolution operations on images, features of different scales are obtained. AKAZE-CHARBONNIER uses the scale space formula as follows:

[0024]

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

[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 They represent the statistical features of the mean, variance, maximum and minimum values ​​within a local window of size k*k respectively.

[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 whose nearest neighbor and next nearest neighbor Euclidean distance 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 next nearest neighbor feature descriptors in the reference image, and the calculation formula is as follows:

[0030]

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

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

[0033]

[0034] Where g(di) is the point d i The bootstrap confidence at , ranging from 0 to 1, the parameter α controls the steepness of the bootstrap function and is used to adjust the change in confidence, μ represents the center or mean of the point.

[0035] The guided maximum likelihood stochastic consensus algorithm (G-MLESAC) is used to purify the candidate matching set. The threshold T is set to remove the matching points with low confidence. When g(di)>T, the point is included in the model estimation. The matching point set D` after screening by the guided mechanism is:

[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+1th iteration is:

[0038]

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

[0040] In a preferred embodiment of the present invention, in step S6, the candidate matching set that does not meet the optimal model parameter θ` is regarded as an outlier and is eliminated, and the candidate matching set that meets the optimal model parameter θ` is regarded as an inlier and retained to obtain the best matching set for subsequent image stitching.

[0041] In a preferred embodiment of the present invention, in step S7, 4 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, and the formula is as follows:

[0042] H = arg min||Ax-b||, where A is the matrix constructed by randomly selecting 4 pairs of points from the matching points to construct a linear equation system, and b is the value on the right side;

[0043] Select the perspective matrix H containing the largest number of internal points as the optimal perspective transformation matrix H best :

[0044]

[0045] Perform perspective transformation on the image to be stitched so that it is aligned 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 function E(x,y), the optimal suture line position is found using dynamic programming:

[0049]

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

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

[0052]

[0053] Where: G Source x (x,y),G Source y (x, y) represent the gradient values ​​of the image to be stitched in the x and y directions obtained by the Sobel algorithm, G warped x (x,y),G warped y (x, y) represents the gradient value of the reference image in the x and y directions respectively.

[0054] Image color difference intensity value E color The calculation formula for (x, y) 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) represent the pixel value difference between the image to be stitched and the reference image.

[0057] In a preferred embodiment of the present invention, in step S8, determining the best stitching line and obtaining a seamless stitching panorama based on multi-band weighted fusion specifically includes:

[0058] Using the Gaussian weight function w(x,y), a weight is assigned to each pixel point (x,y) in the stitching 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 indicates the discrete degree of the value, affects the decay speed and influence range of the weight, and determines the degree of image smoothing. The larger the σ value, the stronger the smoothing effect, and the image noise is strongly removed, but it will also be accompanied by the loss of some image details; conversely, the smaller the σ value, the weaker the smoothing effect, and more image details are retained, but the denoising effect is weakened;

[0061] The weight function w(x,y) determines the contribution of each pixel. The final pixel value Ifinal(x,y) at 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 position (x, y), W i (x, y) represents the weight of the i-th image at position (x, y), which is determined by the weight function w(x, y).

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

[0065] By performing edge analysis on the complete panoramic image Ipanorama generated after fusion, the upper, lower, left and right boundaries of the cropped area are obtained. The calculation formula is as follows:

[0066]

[0067] When cropping the edge of the panorama, a buffer area (padding) is added to prevent the image content from being cropped incorrectly. The calculation formula for the upper, lower, left and right boundaries of the buffer area is as follows:

[0068]

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

[0070]

[0071] Wherein, 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.

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

[0073] (1) When the AKZAZE algorithm constructs a nonlinear scale space for feature detection, the present invention introduces the CHARBONNIER nonlinear diffusion model with complete function convexity of the energy functional to accelerate the solution of the nonlinear partial differential equation, ensuring that the equation has a unique global optimal solution and the stability of the diffusion process. The feature detection algorithm is more robust when dealing with noise or artifacts, and can effectively suppress noise and avoid oscillation and instability in the solution space. In addition, in the diffusion equation of image denoising, the complete function convexity enhances the smoothness of the solution, thereby effectively reducing the generation of artifacts in the image. Secondly, the present invention extracts the depth value of a specific point and generates a deep feature descriptor in combination with the depth information of the surrounding area, and then combines it with the local binary feature descriptor to generate a multi-scale deep hybrid feature descriptor HMD to describe the feature point. This combination of multimodal features enhances the adaptability of the descriptor to complex scenes and improves the subsequent splicing accuracy. The M-LDB descriptor can effectively capture local texture information and has local scale invariance and rotation invariance, while the deep 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 mismatch rate, accelerates the feature matching process, and improves the geometric consistency and visual quality of the stitching results.

[0074] (2) In the data purification stage of the candidate matching set, the present invention adopts a guided maximum likelihood random consensus algorithm to purify the data. The estimation of the model parameters is optimized by designing a guiding mechanism, and the confidence is assigned to the matching points using a guiding function, and the matching points with low confidence are eliminated, and then iterative calculation is performed 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 mismatching, significantly reducing the amount of data required for 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 matching pairs that conform to the data distribution, and further improve the robustness of the model. Through iterative calculation, the algorithm can gradually improve the model parameters, repeatedly verify and update the matching results, thereby further enhancing the robustness of the splicing algorithm.

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

[0076] Figure 1 This is an overall flow chart of a highly robust image stitching method of the present invention;

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

[0078] Figure 3 It is a panoramic image stitching result diagram of each method under a top-down perspective according to an embodiment of the present invention. DETAILED DESCRIPTION

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

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

[0081] Please refer to Figure 1As shown, the present invention first divides the pre-processed sequence images into reference images and images to be stitched, uses the AKAZE-CHARBONNIER algorithm to construct a nonlinear scale space to obtain the mixed feature descriptors in the two images for subsequent matching, and then combines the BF matching with the KNN algorithm to screen out the candidate matching set, design a guiding mechanism, and use the guided maximum likelihood random consistency algorithm to purify the candidate matching set to obtain the best matching set. Then calculate the perspective transformation between the best matches, solve the best perspective matrix for the image stitching link, use the stitching path intensity difference value function to determine the best stitching line position, use multi-band weight fusion to obtain a seamless stitching panorama, and finally determine the edge cropping area and output the final output panoramic image stitching result.

[0082] See also Figure 1 , the preprocessed sequence image set is divided into a reference image and an image to be spliced. 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 spliced ​​have sufficient discernibility; the AKAZE-CHARBONNIER algorithm is used for feature extraction. The proportional coefficient K value in the CHARBONNIER nonlinear diffusion model can determine how much edge information can be retained in the image area. The larger the K value, the greater the diffusion intensity, and the less edge information is retained.

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

[0084] Each grid size can capture feature information at different scales. Larger grids are used to eliminate high-frequency noise, while 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 contain scale invariance, rotation invariance and depth information of the image. The Euclidean distance between the nearest neighbor and the next nearest neighbor feature point between the key point of the reference image and the image to be stitched is calculated, and the matching pairs with a ratio greater than γ are screened out to obtain a candidate matching set. The confidence of the candidate matching points is assigned by the guiding function, and the matching points with a confidence lower than T are removed and iterative calculation is performed to obtain the optimal model parameters. The points that meet the best model are regarded as inliers, and those that do not meet the model are regarded as outliers. The perspective matrix containing the largest number of inliers is calculated for image stitching, and the best stitching line position is found by calculating the stitching path intensity difference value function. The Gaussian weight function is used to assign weights to each pixel point in the stitching area, and then image fusion is performed to obtain a seamless stitching panorama, where the weight assigned to the area close to the best stitching line is lower, and the weight assigned to the area far from the best stitching area is higher. It is necessary to first find the minimum bounding rectangle boundary of all edge pixels that are considered to be black edges or irrelevant areas, and a buffer area must be set to prevent the image content from being incorrectly cropped. After edge cropping, the final panoramic image stitching result image is output.

[0086] The present invention provides a highly robust image stitching method, which specifically comprises the following steps:

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

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

[0089] The AKAZE-CHARBONNIER algorithm generates multi-scale images through nonlinear partial differential equations, and adjusts the intensity of diffusion and the degree of retention of image edge information by setting the diffusion coefficient C and the contrast coefficient K. The mathematical expression of the nonlinear 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 position (x,y) and time t, Represents the gradient of the image.

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

[0093]

[0095]

[0096] Among them, the function Representation and Image Gradient The related diffusion coefficient controls the diffusion behavior of the image in different regions. The parameter K is the contrast coefficient that determines the diffusion level. The CHARBONNIER nonlinear diffusion model parameters of formula (3) are introduced 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. The complete functional convexity can be proved by calculating the energy functional and its second-order derivative of the CHARBONNIER nonlinear diffusion model.

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

[0098]

[0099]

[0100] From formula (5), we can see that if the proportional coefficient K ≥ 0, the second-order derivative This means that the energy functional of the nonlinear diffusion model has complete functional convexity, which can accelerate the solution of the nonlinear partial differential equation (1) and ensure the stability of the nonlinear diffusion process.

[0101] AKAZE-CHARBONNIER generates images of different scales through nonlinear scale space methods. By continuously performing convolution operations on images, features of different scales are obtained. The scale space formula used by AKAZE-CHARBONNIER is:

[0102]

[0103] Among them, S x , S y Represent the scale changes in the x and y directions respectively.

[0104] S3. Use a multi-scale deep hybrid feature descriptor HMD to describe feature points.

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

[0106]

[0107] Among them, I xx ,I yy They 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.

[0108] Taking the key point as the center, the gray value difference of the area is calculated within its domain, and binary encoding is performed according to the gray value to form a fixed descriptor vector D, and the descriptor is normalized to obtain the local binary descriptor D M-LDB (P):

[0109]

[0110]

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

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

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

[0114] Where 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:

[0115]

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

[0117] For each detected feature point p, extract the depth value of the point and combine it with 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 the depth information around the feature point is considered, a local window W(p) of size k*k can be defined:

[0118]

[0119] The statistical features within the calculation window include: mean μ D , variance σ D 2 , maximum value D max and the minimum value D min :

[0120]

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

[0122]

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

[0124]

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

[0126] S4. Use BF matching to obtain the initial matching set and use the KNN algorithm to eliminate the matches whose nearest neighbor and next nearest neighbor Euclidean distance is greater than γ to obtain the candidate matching set.

[0127] BF matching searches for the best match between feature points by traversing all possible matching pairs. For each feature point, BF matching calculates the distance to all other feature points and selects the feature point with the smallest distance as the matching object. It can find the global best match but is only suitable for a small number of feature points. 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 means that the match is reliable. The present invention uses the KNN algorithm to eliminate matches whose nearest neighbor and second nearest neighbor Euclidean distance are greater than γ to obtain the calculation of the candidate matching set as shown in formula (15).

[0128]

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

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

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

[0132]

[0133] Among them, P(d i |θ) represents that for each matching point d i , calculate its probability according to the model parameters θ. The likelihood function simplifies the calculation by logarithmic transformation and obtains the optimal model parameters θ`:

[0134]

[0135] The present invention designs a guidance mechanism to optimize the model parameter estimation, through the guidance function g(d i ) to define the matching point confidence:

[0136]

[0137] Among them, g(d i ) is the guided confidence at point di, ranging from 0 to 1. The parameter α controls the steepness of the guided function and is used to adjust the change in confidence. μ represents the center or mean of the point.

[0138] By setting the threshold T to remove low confidence matching points, only when g(d i )>T, the point is included in the model estimation, and the matching point set D` after screening by the guidance mechanism is:

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

[0140] The parameter T is used to determine which matching points are credible and incorporate them into the model estimation. Regarding the value principle of parameter T, when the splicing scene is a simple scene, 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 extremely reliable matching points are retained; when the splicing scene is a complex scene, that is, the feature points are fuzzy and 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 subsequent iterative optimization. In short, choosing a suitable threshold T is a 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+1th iteration is:

[0141]

[0142] Among them, θ is the current given model parameter, and L(θ|D`) represents the likelihood function of the matching point set D` after screening by the guidance mechanism.

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

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

[0145] S7. Align images through perspective transformation and find the optimal stitching line position.

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

[0147]

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

[0149] Among them, h ij It is the matrix element of perspective change, 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.

[0150] Each iteration will recalculate the perspective matrix H, and the matching points that meet the optimal parameter model θ' are regarded as inliers, and those that do not meet the optimal parameter model θ' are regarded as outliers. The perspective matrix H containing the largest number of inliers is regarded as the optimal perspective transformation matrix H. best :

[0151]

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

[0153] Warped Image=H best Source Image (24)

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

[0155]

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

[0157]

[0158] Where: G Source x (x,y),G Source y (x, y) represent the gradient values ​​of the image to be stitched in the x and y directions obtained by the Sobel algorithm, G warped x (x,y),G warped y (x, y) represents the gradient value of the reference image in the x and y directions respectively.

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

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

[0161] Where ISource(x, y) and Iwaeped(x, y) represent the pixel value differences between the image to be stitched and the reference image respectively.

[0162] Where Egeometry(x,y) is the intensity value of the image structure difference, and Ecolor(x,y) is the intensity value of the image color difference. GSourcex(x,y) and GSourcey(x,y) represent the gradient values ​​of the image to be stitched in the x and y directions obtained by the Sobel algorithm. Similarly, Gwarpedx(x,y) and Gwarpedy(x,y) represent the gradient values ​​of the reference image in the x and y directions. ISource(x,y) and Iwaeped(x,y) represent the pixel value differences between the image to be stitched and the reference image. The gradient values ​​Gx and Gy of the Sobel operator in the x and y directions are:

[0163]

[0164] Where I is the pixel value of the image. The intensity difference values ​​of all points on each stitching path are calculated, and the path with the smallest intensity difference is selected as the best stitching line.

[0165] S8. After determining the best stitching line, multi-band weighted fusion is performed to obtain a seamless stitching panorama.

[0166] After confirming the best stitching line, it is necessary to perform multi-band weighted fusion (Multi-bandblending) on ​​the image to be stitched and the reference image, using different frequency bands to fuse the images, and assigning a weight to each pixel in the image. The area close to the best stitching line has a lower weight, while the area far from the best stitching line has a higher weight. In the process of stitching and stitching, there are common parts between the two images, which need to be weightedly 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:

[0167]

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

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

[0170]

[0171] Among them, I i (x, y) represents the pixel value of the i-th image at position (x, y), and wi(x, y) represents the weight of the i-th image at 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 that the weights are normalized, making the stitching of overlapping areas smoother and reducing the abruptness of the boundaries.

[0172] S9, output the final result after cropping the edges of the panorama.

[0173] The panorama generated after fusion often has black edges, and the edges need to be detected and identified for cropping. If the complete panoramic image generated after fusion is defined as Ipanorama, the cropping area is obtained by analyzing the edge of Ipanorama, and the minimum bounding rectangle of all edge pixels that are 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:

[0174]

[0175] A buffer area (padding) needs to be added to prevent the image content from being incorrectly cropped. The boundaries of the buffer area are:

[0176]

[0177] The panorama is cropped according to the area defined by the cropping boundary. The final stitched panoramic image Icropped output after edge cropping can be expressed as:

[0178]

[0179] The following is an illustration of a highly robust image stitching method provided by the present invention, and the stitching effect of a panoramic image under a side view and a top view: The operating conditions of the embodiment: 64-bit Windows 10 operating system, Intel (R) Core (TM) i5-9400F processor, CUDA 10.1 version of the GPU acceleration library, environment: Visual Studio 2015, OpenCV 4.4.0, OpenCV_contrib 4.4.0, cmake-1.17.2-win64-x64, etc.

[0180] The public data set used in the embodiment is the "UC Merced Land Use Dataset". Taking into account factors such as geometric deformation, lighting differences, feature matching difficulties, and background complexity that may be caused by different perspective changes, two types of images under side view and top view were selected for stitching experiments, and a comparative experiment was set up. In the comparison method, the mature SIFT algorithm was used for feature extraction, and the BF algorithm and KNN algorithm were combined to obtain the candidate matching point set. At the same time, the classic random consensus sampling algorithm (RANSAC) was used for data purification, and the optimal perspective transformation matrix was obtained based on the optimal model parameters calculated by RANSAC. In the overlapping area, the linear weighted progressive fusion method was used to generate pixels, and the final panoramic stitching image was output after edge cropping. Figure 2 and Figure 3 The panoramic stitching results of the comparison method and the method of the present invention under the side view and top view are shown. Figure 2 and Figure 3 It can be observed from the local magnified image that the method of the present invention significantly reduces artifacts, chromatic aberration and splicing misalignment under different viewing angles. This is due to the feature descriptor provided by the AKAZE-CHARBONNIER algorithm, which contains multi-scale and depth information of the current scene.

[0181] Improves the geometric accuracy of image stitching. Compared with the RANSAC method, the present invention uses G-MLESAC for data purification. G-MLESAC is based on the principle of maximum likelihood estimation, which can more effectively utilize 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 accuracy of subsequent perspective transformation matrix solutions. Among the comparative methods, although the linear weighted gradual fusion method is simple and effective in processing overlapping areas, Figure 2 and Figure 3 However, the method in the previous paper shows obvious gradient edges, which makes the transition of the fused area unnatural. The optimal stitching line and multi-band weighted fusion method adopted by the present invention significantly improves the transition effect while maintaining the details, alleviating the shortcomings of the linear weight method, making the stitching area visually flatter and the transition of the fused area more natural. These improvements make the stitching effect more excellent and enhance the overall quality and visual experience of the image.

[0182] The device embodiments described above are only illustrative. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as: 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 can be through some interfaces, and the indirect coupling or communication connection of the devices or units can be electrical, mechanical or other forms.

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

[0184] In addition, all functional units in the embodiments of the present invention may be integrated into one processing unit, or each unit may be separately used as a unit, or two or more units may be integrated into one unit; the above-mentioned integrated units may be implemented in the form of hardware or in the form of hardware plus software functional units.

[0185] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A highly robust image stitching method, characterized in that: The steps include: S1, dividing the preprocessed sequence image set into reference images and images to be stitched; S2, uses the AKAZE-CHARBONNIER algorithm to construct a nonlinear scale space for feature detection; S3, uses a multi-scale deep hybrid feature descriptor HMD to describe feature points; S4, using BF matching to obtain the initial matching set, and using the KNN algorithm to eliminate the matching sets whose nearest neighbor and next nearest neighbor Euclidean distance are greater than γ to obtain the candidate matching set; S5, based on the bootstrapping mechanism, uses the bootstrap maximum likelihood random consensus algorithm to purify the candidate matching set, and obtains the optimal model parameters through iterative calculation; S6, eliminate the matching sets that do not meet the optimal parameter model to obtain the best matching set; S7, aligns images through perspective transformation and finds the best stitching line position; S8, determine the best stitching line based on multi-band weighted fusion to obtain a seamless stitched panorama; S9, the panorama is cropped at the edges and the final result is output.

2. The highly robust image stitching method according to claim 1, characterized in that: In step S2, the AKAZE-CHARBONNIER algorithm is used to construct a nonlinear scale space for feature detection, which specifically includes: The nonlinear partial differential equation for nonlinear diffusion of images is described as: Among them, I(x,y,t) represents the value of the image at the (x,y) position and time t, ▽I represents the gradient of the image, and C(||▽I||) represents a diffusion coefficient related to the gradient size. The formula is as follows: The function G(▽|I(x,y,t)|) represents the diffusion coefficient related to the image gradient ▽I. It is a nonlinear function that depends on the gradient of the image at the (x,y) position and time t. It controls the diffusion behavior of the image in different areas. The diffusion is small in the area with large gradient at the edge of the image, while in the flat area with small gradient, the diffusion is large, thereby achieving effective denoising and edge preservation. The parameter K is the contrast coefficient that determines the diffusion level. The formula The parameters of the CHARBONNIER nonlinear diffusion model are introduced into the diffusion coefficient C(||▽I||) to approximately solve the nonlinear partial differential equation; Among them, AKAZE-CHARBONNIER generates images of different scales through nonlinear scale space methods. By continuously performing convolution operations on images, features of different scales are obtained. AKAZE-CHARBONNIER uses the scale space formula as follows: Where sx and sy represent the scale changes in the x and y directions, respectively.

3. The highly robust image stitching method according to claim 2, characterized in that: In step S3, a multi-scale depth hybrid feature descriptor HMD is used to describe the 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 They represent the statistical features of the mean, variance, maximum and minimum values ​​within a local window of size k*k respectively.

4. The highly robust image stitching method according to claim 3, characterized in that: In step S4, BF matching is used to obtain an initial matching set and the KNN algorithm is used to eliminate the matches whose nearest neighbor and next nearest neighbor Euclidean distance 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 next nearest neighbor feature descriptors in the reference image. The calculation formula is as follows: Among them, D1 and D2 represent the feature descriptors of the nearest neighbor and the second nearest neighbor in the reference image respectively, and Di is used to describe the component of any feature descriptor in the image to be spliced.

5. The highly robust image stitching method according to claim 4, characterized in that: In step S5, a guidance mechanism is designed to optimize the estimation of model parameters, and the guidance confidence of the matching point is defined by the guidance function g(di), and the calculation formula is as follows: Where g(di) is the point d i The guide confidence at , the value range is between 0 and 1, the parameter α controls the steepness of the guide function, which is used to adjust the change of confidence. μ Indicates the center or mean of the point; The guided maximum likelihood stochastic consensus algorithm (G-MLESAC) is used to purify the candidate matching set. The threshold T is set to remove the matching points with low confidence. When g(di)>T, the point is included in the model estimation. The matching point set D` after screening by the guided mechanism is: D`={d i |g(d i )>T} 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+1th iteration is ` for: Among them, θ is the current given model parameter, and L(θ|D`) represents the matching point set D after being screened by the guidance mechanism. ` The likelihood function of .

6. The highly robust image stitching method according to claim 5, characterized in that: In step S6, the candidate matches that do not meet the optimal model parameters θ are ` is regarded as an outlier and removed, and the optimal model parameter θ ` They are regarded as inliers and retained to obtain the best matching set for subsequent image stitching.

7. The highly robust image stitching method according to claim 6, characterized in that: 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, which is as follows: H=argmin||Ax-b|| Where A is the matrix constructed by randomly selecting 4 pairs of points from the matching points to construct a linear equation system and decomposing it through singular value decomposition (SVD), and b is the value on the right side; Select the perspective matrix H containing the largest number of internal points as the optimal perspective transformation matrix H best : Perform perspective transformation on the image to be stitched so that it is aligned with the reference image: Warped Image=H best ·Source Image Among them, Source Image is the image to be stitched, and Warped Image is the reference image; According to the suture path strength difference function E(x,y), the optimal suture line position is found using dynamic programming: Among them, E geometry (x, y) is the intensity value of the structural difference between images, and Ecolor(x, y)) is the intensity value of the color difference.

8. The highly robust image stitching method according to claim 7, characterized in that: Image structure difference intensity value function E geometry The calculation formula for (x, y) is as follows: Where: G Source x (x,y),G Source y (x, y) represent the gradient values ​​of the image to be stitched in the x and y directions obtained by the Sobel algorithm, G warped x (x,y),G warped y (x, y) represents the gradient value of the reference image in the x and y directions respectively; The calculation formula of the image color difference intensity value Ecolor(x, y) is as follows: E color (x,y)=I Source (x,y)-I Warped (x,y) In the formula, I Source (x, y), I waeped (x, y) represent the pixel value difference between the image to be stitched and the reference image.

9. The highly robust image stitching method according to claim 8, characterized in that: In step S8, the best stitching line is determined, and a seamless stitching panorama is obtained based on multi-band weight fusion, which specifically includes: Using the Gaussian weight function w(x,y), a weight is assigned to each pixel point (x,y) in the stitching area. The calculation formula of the Gaussian weight function w(x,y) is as follows: Among them, the parameter σ is the standard deviation of the Gaussian distribution, which indicates the discrete degree of the value, affects the decay speed and influence range of the weight, and determines the degree of image smoothing. The larger the σ value, the stronger the smoothing effect, and the image noise is strongly removed, but it will also be accompanied by the loss of some image details; conversely, the smaller the σ value, the weaker the smoothing effect, and more image details are retained, but the denoising effect is weakened; The weight function w(x,y) determines the contribution of each pixel. The final pixel value Ifinal(x,y) at 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 position (x, y), W i (x, y) represents the weight of the i-th image at position (x, y), which is determined by the weight function w(x, y).

10. The highly robust image stitching method according to claim 9, characterized in that: In step S9, the panorama image is cropped at the edges and the final result is output, which specifically includes: By performing edge analysis on the complete panoramic image Ipanorama generated after fusion, the upper, lower, left and right boundaries of the cropped area are obtained. The calculation formula is as follows: When cropping the edge of the panorama, a buffer area (padding) is added to prevent the image content from being cropped incorrectly. The calculation formula for the upper, lower, left and right boundaries of the buffer area is as follows: Among them, Top, Bottom, Left, and Right represent the upper, lower, left, and right boundaries of the cropping area excluding the buffer area, respectively. The panorama is cropped according to the area defined by the cropping boundary. The final stitched panoramic image Icropped output after edge cropping is expressed as: Wherein, 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