Point cloud and image registration method based on improved RIFT

Abstract

The invention discloses a point cloud and image registration method based on improved RIFT. Firstly, a point cloud is converted into a depth image, and background removal is carried out on the depth image and a rock debris image; and then, calculating the phase consistency of the RGB image and the depth image by using a phase consistency model, and carrying out feature point detection according to a calculation result. Thirdly, evaluating the preliminary registration feature points through the ratio of nearest neighbor to secondary neighbor and regional consistency; and finally, optimizing feature point matching by using an improved RANSAC (Random Sample Consensus) algorithm, and calculating a homography matrix to determine a mapping relation between the RGB image and the depth image, so as to realize registration of different-source data. Compared with the original RIFT registration, the improved method has the advantages that the residual error of the projection point and the reference point after matching is better, and the registration result is more accurate. The method can be used for rock debris particle analysis and identification.

Description

TECHNICAL FIELD

[0001] The application relates to a point cloud and image registration method based on RIFT, in particular to an improved point cloud and image registration method based on RIFT. BACKGROUND

[0002] In oil and gas exploration, cuttings logging is a key means of obtaining formation information, which has the advantages of rapid operation, simplicity and low cost. Through cuttings logging technology, the status of oil and gas resources in the formation can be effectively reflected, and through systematic observation and analysis, geologists can understand the formation structure and its oil content. The cuttings particles collected during the cuttings logging process have important significance in geological research. With the development of digital imaging and image processing technology, the analysis and identification of cuttings logging are gradually moving towards digitization. Compared with manual detection, digital technology provides higher efficiency and accuracy in cuttings lithology identification and particle size analysis, reducing human interference. Since the cuttings particles are densely arranged in the sample, it is difficult to meet the engineering requirements by relying only on image data for analysis and identification, so it is necessary to introduce point cloud depth information to improve the analysis effect of cuttings particles, and the registration of point cloud and image is the basis for the identification and analysis of cuttings particles.

[0003] Since the cuttings particles are usually small and densely packed during actual collection, it is difficult to extract and classify these particles. In order to solve the problem of insufficient precision in cuttings particle edge detection and target extraction, this study combines two-dimensional image and point cloud depth information technology. A laser profile scanner is used to obtain point cloud data of cuttings particles, and at the same time, a camera fixed beside the scanner is used to collect cuttings images in the same field of view. The registration of these two heterogeneous data is the key to combining their advantages. Laser radar technology provides rich spatial depth and geometric attribute data due to its active detection capability and low sensitivity to the environment. Therefore, this study proposes a point cloud and image registration method based on improved RIFT algorithm to realize the effective fusion of the depth information of point cloud and the color texture information of two-dimensional image. SUMMARY

[0004] The application provides a point cloud and image registration method based on improved RIFT. The improvement strategy first respectively pre-processes the point cloud and image data, including generating a depth map by index mapping of the point cloud data, background filtering and denoising operations of the two kinds of data; then, phase consistency models are used to respectively calculate the phase consistency of the RGB image and the depth image, and key point features are detected according to the calculation results; subsequently, MIM maximum index maps are used for feature description, and the nearest neighbor ratio and region consistency method is used to preliminarily realize the matching of the feature points; finally, an improved RANSAC algorithm is used to further optimize the matching results, and a homography matrix H is calculated to obtain the mapping relationship between the RGB image and the depth image, thereby realizing the registration of the heterogeneous images. BRIEF DESCRIPTION OF DRAWINGS

[0005] Figure 1 The point cloud and image registration flowchart of the application. DETAILED DESCRIPTION

[0006] The application will be further described below, including the following steps: Figure 1 The application will be further described below, including the following steps:

[0007] Firstly, the point cloud and the image are pre-processed, the point cloud is converted into a depth image, and then the background of the depth image and the debris RGB image is removed.

[0008] For point cloud preprocessing, the size of the debris point cloud is consistent, and the specific specifications are as follows: the scanned debris target area is a square dense debris particle area of about 6cm*6cm, the size of the scanning movement direction y is set to n=4000 scanning lines, and the vertical direction of the movement direction x is fixed at 3840 points. The scanning device moves horizontally in the y-axis direction and moves relative to the target object in the horizontal direction. Finally, a collected point cloud contains 15360000 (x, y, z) coordinates; the points in the point cloud are arranged according to the index i. Wherein i=1,…,15360000, the pixel point coordinates (x i ,y i ) on the XOY plane satisfy the following equation relationship:

[0009] i=3840(y i +1)+(x i +1){x=0,...,3839、y=0,...,3998}

[0010] Based on the mapping of the point index and the two-dimensional pixel coordinates, the depth value of each pixel can be calculated to complete the generation of the point cloud to the depth map. In order to standardize the Z value of the point cloud, the minimum elevation value Z min and the maximum elevation value Z max of all points in the point cloud need to be determined, and then the elevation value Zi The value can be scaled to a range of 0 to 255 using the following formula:

[0011]

[0012] Where D i Represents the pixel coordinates (x) in the depth image i ,y i The grayscale value at point ) can be used to generate a depth map of a uniform and ordered point cloud using the above method.

[0013] Background removal is performed on the generated point cloud depth map based on the target region and depth information. First, the target region of the rock debris is determined using the depth value of the horizontal plane where the rock debris particles are located, and the parts outside the target region are marked as background. At the same time, void regions are identified and segmented using depth information as background. Finally, a median filtering method is used to remove discrete points in the depth image.

[0014] In the image preprocessing stage, background filtering of rock debris images is performed using techniques such as threshold segmentation and connected component analysis.

[0015] The second step involves applying a phase consistency model to calculate the phase consistency of the RGB and depth images, and then performing feature point detection based on the calculation results.

[0016] For the calculation of the phase consistency model, it is necessary to build N s Each scale and N d A two-dimensional logarithmic Gabor filter in each direction is defined as follows:

[0017]

[0018] Where (ρ, θ) represents logarithmic polar coordinates, m and d refer to the scale and orientation parameters of the two-dimensional Log-Gabor filter, respectively, and σ ρ and σ θ ρ and θ represent the bandwidths of the polar coordinates ρ and θ, respectively, and (ρ,θ,m,d) represents the center frequency of the two-dimensional Log-Gabor filter.

[0019] Assuming I(x,y) represents a two-dimensional image, the two-dimensional logarithmic Gabor filter is defined in the spatial domain as follows:

[0020] G(x,y,m,d)=G e (x,y,m,d)+iG o (x,y,m,d)

[0021] The imaginary part G o It is an odd-symmetric Log-Gabor wavelet with real part G. e It is an even-symmetric Log-Gabor wavelet.

[0022] Convolve I(x,y) with both odd and even symmetric wavelets respectively to obtain the odd symmetric response O. md (x,y) and even-symmetric response E md (x,y):

[0023] O md (x,y)=I(x,y)*G o (x,y,m,d)

[0024] E md (x,y)=I(x,y)*G e (x,y,m,d)

[0025] Where * denotes convolution operation. Then, the amplitude and phase components of I(x,y) at scale m and direction d are expressed as follows:

[0026]

[0027] When calculating the phase consistency model, it is necessary to integrate the analysis results from various directions and scales, and add noise compensation N. The calculation formula for the phase consistency model is as follows:

[0028]

[0029] Where w d (x,y) are the weighted function components. It is a small value. The >·> operator is to prevent the inner enclosing variable from being negative. If the inner enclosing variable takes a negative value, the enclosing variable is set to 0. N is noise compensation, and ΔΦ is a phase deviation function, the specific definition of which is as follows:

[0030]

[0031] Based on the phase consistency model described above, the torque of the image is calculated. The axis with the smallest torque value is the principal axis; conversely, the axis perpendicular to the principal axis is where the maximum torque is located. Assuming the principal axis is τ, the specific calculations for the principal axis, minimum torque, and maximum torque are as follows:

[0032] τ = 0.5arctan(b / (ac))

[0033]

[0034] Where a and c represent the second-order central moments in the directions parallel to the x and y axes, respectively, and b represents the mixed second-order central moment between the x and y axes. The variables t1, t2, and t3 are calculated as follows:

[0035]

[0036] The above analysis shows that the minimum torque m τ Points with larger values ​​are most likely angular features, with the maximum torque M. τ Regions with larger values ​​are the edge regions of the image. In corner detection, the maximum pixel value in the minimum moment map is first determined, denoted as m1. Then, all points in the map with pixel values ​​greater than 0.01×m1 are selected as candidate corner points, i.e., local maxima, and non-maximum suppression is used to remove redundant candidate points. For edge feature points, the AGAST detector is used to extract edge feature points from the maximum moment map. Finally, the extracted corner points and edge points are combined for subsequent feature matching.

[0037] When constructing the maximum index map, a background board with the same size as the original image is first created. Then, for each pixel at (x,y) on the background board, the amplitude vector values ​​in N directions are calculated, and the magnitudes of these amplitude vectors are compared. The value of the maximum amplitude vector is assigned as the corresponding value of the pixel in the maximum index map, thereby generating the final maximum index description.

[0038] After constructing the maximum index map, the feature vector is then represented using a distribution histogram. For each feature point, an n×n local pixel block is selected, and a Gaussian function with a standard deviation of n / 2 is used to assign weights to each pixel. The local region is then divided into N×N sub-grids, and a histogram distribution with N components is calculated for each sub-grid. Finally, these histograms are concatenated to generate the feature vector.

[0039] Multiple sequences are created by adjusting the order of different directions within the same scale. Then, multiple different maximum index maps are generated based on these sequences. A rotation-invariant maximum index map is constructed by changing the index order. The feature vector is then described using distribution histogram technology. Specifically, for each detected keypoint, a small neighborhood is extracted around it. A 96×96 window is selected as a fixed neighborhood centered on each keypoint. The neighborhood is then rotated to the main direction of the keypoint. Since there are eight different directions, this neighborhood is divided into 8×8 blocks. Histogram statistics are performed on the distribution of the eight directions in each block to finally achieve feature description.

[0040] The third step is to perform preliminary registration of feature points using the nearest neighbor to second nearest neighbor ratio and regional consistency assessment methods.

[0041] By calculating the distance ratio between the nearest and second nearest neighbors and setting a threshold for comparison, a smaller distance ratio indicates a larger gap between the nearest and second nearest neighbors, making the matching result more likely to be correct. For incorrect matches, the nearest and second nearest neighbors are closer together, and the ratio is closer to 1.

[0042] For each point pair obtained from the above matching, two neighborhood windows are first set: a large window of 200×200 and a small window of 32×32. Under the large window, it is only necessary to determine whether the large window regions of the two feature points intersect. If they intersect, the point pair is considered to meet spatial consistency; otherwise, the matched point pair is directly discarded. Based on satisfying spatial consistency, the consistency of region features is further evaluated under the small window using the Normalized Cross-Correlation (NCC) method. This method first calculates the average pixel value of the two small window regions. The specific calculation process is as follows:

[0043]

[0044] Where N = 32, A and B represent the average pixel values ​​of the two small window regions, and A(i,j) and B(i,j) represent the pixel values ​​at coordinates (i,j) in the two window regions, respectively. Then, normalized cross-correlation is calculated based on the average pixel values ​​of the regions:

[0045]

[0046] Where σ A and σ B It is the standard deviation of pixel values ​​in regions A and B, calculated as follows:

[0047]

[0048] The NCC value ranges from -1 to 1, where 1 indicates a perfect match between two regions, -1 indicates complete opposites, and 0 indicates no correlation. Higher NCC values ​​indicate stronger similarity. Experiments determined the optimal threshold to be 0.7. Provided that the large window region satisfies spatial consistency, if the calculated NCC value for the small window region is greater than or equal to 0.7, the point pair is considered to meet the regional feature consistency evaluation criteria and is retained; otherwise, the matching point pair is discarded.

[0049] The fourth step involves using an improved RANSAC algorithm to optimize the feature point matching results and calculating the homography matrix to establish a mapping relationship between the RGB image and the depth image, thereby achieving registration of heterogeneous data.

[0050] Even in the initial matching of point pairs, a small number of mismatches may still exist. By improving the RANSAC algorithm, these inaccurate point pairs can be further eliminated, thus calculating the mapping relationship between images more accurately. The RANSAC algorithm iteratively optimizes observation data containing outliers, estimates model parameters, and uses them to calculate the transformation matrix between two images in image registration. Its specific steps are as follows:

[0051] (1) First, sort all matching points in descending order according to NCC value, and apply the selection strategy in this order. The higher the NCC value, the greater the probability that the point pair is a correct match. Select the minimum number of points to estimate the transformation model. The homography transformation requires 4 point pairs to calculate the projection transformation matrix H.

[0052] (2) Transform all matching point pairs using the estimated transformation model and calculate the transformation error;

[0053] (3) Set a threshold, mark points smaller than the threshold as interior points, and consider these points as correct matches; otherwise, mark them as exterior points;

[0054] (4) Repeat the above steps, each time randomly select a new pair of points, count the number of interior points, and use Huber estimation to fine-tune the model in each iteration;

[0055] (5) Finally, the model with the most interior points is selected as the required transformation model.

[0056] To ensure the reliability, or confidence level, of the final estimated model, the number of iterations n must be at least:

[0057]

[0058] Where p refers to the probability that all four pairs of matching points are interior points in at least one sampling, η is the probability of mismatched points, and L is the minimum value in the estimated transformation model.

[0059] The RANSAC parameter model estimation algorithm first needs to ensure the confidence level of the estimated model. To achieve a higher confidence level in the final estimated model, the RANSAC algorithm is improved in the following two ways:

[0060] (1) Optimize the sample point selection strategy: Use the normalized cross-correlation (NCC) value of the regional features of the matching point pairs as the selection criterion. First, sort all the NCC values ​​in descending order, and then select according to the sorting results, instead of the original random selection method. Point pairs with larger NCC values ​​have a higher probability of matching correctly.

[0061] (2) By combining local optimization techniques within the RANSAC framework, the model generated in each iteration is finely adjusted through Huber estimation, thereby improving the accuracy of the final model.

[0062] Affine transformations include operations such as translation, rotation, scaling, and shearing, and are typically represented by a 2×3 matrix containing six parameters. In contrast, projection transformations not only include all affine transformations but also handle perspective effects caused by changes in viewpoint. After feature point extraction, description, and matching, the RANSAC algorithm is used to eliminate erroneous matches. The RANSAC algorithm selects four pairs of matching points each time, with each pair providing two equations and a total of eight unknown parameters. By solving these eight equations, the eight unknown parameters are calculated, thus solving for the eight independent parameters h. 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 We can determine a 3×3 homography matrix H, which determines the projection of the image under the new viewpoint. This is more suitable for handling complex scenes. The specific mathematical expression of the projection transformation is as follows:

[0063]

[0064] For each matching point pair in the improved RIFT algorithm after considering the nearest neighbor to second nearest neighbor ratio, region consistency, and RANSAC parameter estimation, a specific mathematical mapping relationship can be established through projection transformation. Assuming a pair of matching points (x1, y1) and (x2, y2) in two images, the relationship of this pair under projection transformation can be expressed by the following formula:

[0065]

[0066] It can be expressed mathematically as follows:

[0067]

[0068] To solve for H, the RANSAC algorithm requires at least 4 pairs of matching points when estimating the model. In each iteration of model estimation, the transformation is fine-tuned by introducing Huber estimation.

[0069] Huber estimation fine-tunes the estimated transformation in each iteration, primarily by updating the estimated projection transformation matrix H in each iteration. First, the Huber loss function needs to be defined. This function exhibits squared loss in small residual regions and linear loss in large residual regions, and its form is:

[0070]

[0071] Where r is the residual and δ is the preset threshold.

[0072] For each inlier pair, calculate the residual under the current estimated transformation, apply the Huber loss function to each residual to convert it into weights, and apply the weights to each data point to construct a weighted least squares method to update the transformation matrix.

[0073] The residual r for each pair of matching points is defined as the Euclidean distance between the reference point and the projected point:

[0074]

[0075] Where (x3, y3) is the point obtained by mapping from the origin (x1, y1) through the projection transformation matrix H.

[0076] Next, the Huber loss function is applied to each residual r. Small residuals contribute to the total loss in a squared form, while large residuals contribute in a linear form. Then, the Huber loss is minimized through gradient descent to optimize and adjust H. Mathematically, this optimization problem can be expressed as:

[0077]

[0078] During the process of optimizing and adjusting H through gradient descent, the gradient ▽Hr i It can be calculated by r i By obtaining the partial derivative with respect to H, for the projection point (x3, y3) of (x1, y1) after the H transformation, the partial derivative with respect to h can be obtained. 11 Partial derivatives:

[0079]

[0080] Where y3 is relative to h 11 The partial derivative of h is equal to zero because h 11 It only appears in the first row of the projection transformation matrix H, while the projection point y3 is determined by the second row of H.

[0081] Once the gradient is calculated H can then be updated using gradient descent, with the specific calculation formula as follows:

[0082]

[0083] Where α is the learning rate, which is equal to 0.001, a small positive number.

Claims

1. A point cloud and image registration method based on improved RIFT, characterized by the following steps: (1) converting the point cloud into a depth image, and performing background removal processing on the depth image and the debris RGB image; (2) applying a phase consistency model to perform phase consistency calculation on the RGB image and the depth image respectively, and detecting feature points based on the calculation results; (3) using the nearest neighbor and the second nearest neighbor ratio and the regional consistency evaluation method to preliminarily register the feature points; (4) using an improved RANSAC algorithm to optimize the feature point matching results, and calculating a homography matrix to establish a mapping relationship between the RGB image and the depth image, thereby realizing the registration of heterogeneous data.

2. The method of claim 1, wherein In step (1), the point cloud is converted into a depth image, and the depth image and the debris RGB image are subjected to background removal processing. The specific operation method is as follows: For point cloud preprocessing, the size of the debris point cloud remains consistent, assuming that each point in the point cloud is arranged by index i, the search range is i = 1, …, 15360000, x and y are the scanning moving directions, and then the pixel point coordinates (x i ,y i ) on the XOY plane satisfy the following relationship: i = 3840(y i +1)+(x i +1){x = 0,..., 3839, y = 0,..., 3998} Based on the mapping of the point index and the two-dimensional pixel coordinate, the depth value Z of each pixel is calculated; in order to standardize the Z value of the point cloud, the minimum elevation value Z of all points in the point cloud needs to be determined min and the maximum elevation value Z max Then, the elevation value Z i of the point index i can be scaled to the range of 0 to 255 by the following formula: where D i represents the gray value at the pixel coordinate (x i ,y i ) in the depth image, and a uniform and ordered point cloud depth map can be generated by this method; then, according to the target region and the depth information, the background is removed from the generated point cloud depth map by using median filtering; For image preprocessing, the background removal of the debris image can be achieved through image processing techniques such as threshold segmentation and connected region analysis.

3. The method of claim 1, wherein In step (2), the phase consistency model is applied to calculate the phase consistency of the RGB image and the depth image respectively, and the feature points are detected based on the calculation results. The operation details are as follows: For the calculation of the phase coherence model, N s two-dimensional log-Gabor filters in N d directions are built, defined as follows: where (p, q) denotes log-polar coordinates, m and d refer to the scale and orientation parameters of the two-dimensional Log-Gabor filter, and ρ and s θ represent the bandwidths of the polar coordinates p and q, respectively, and (p, q, m, d) denotes the center frequency of the two-dimensional Log-Gabor filter; Assuming that I(x, y) represents a two-dimensional image, for a two-dimensional log Gabor filter, its definition in the spatial domain is as follows: G(x, y, m, d) = G e (x, y, m, d) + iG o (x, y, m, d) where the imaginary part G o is an odd-symmetric Log-Gabor wavelet and the real part G e is an even-symmetric Log-Gabor wavelet; The odd-symmetric response O(x,y) is obtained by convolving I(x,y) with the odd-symmetric wavelet md (x,y) and the even-symmetric response E(x,y) is obtained by convolving I(x,y) with the even-symmetric wavelet md (x,y) : O md (x,y) = I(x,y) * G o (x,y,m,d) E md (x,y) = I(x,y) * G e (x,y,m,d) where * represents convolution operation, and then the amplitude and phase components of I(x, y) at scale m and direction d are represented as: The phase consistency model is calculated by integrating the analysis results from each direction and scale and adding noise compensation N. The specific calculation formula is as follows: where w d (x,y) is a weight function component, is a small value, the <·> operator is to prevent the inner enclosed variable from taking negative values, if the inner enclosed variable takes negative values, the enclosed variable is set to 0, N is a noise compensation, ΔΦ is a phase deviation function, ΔΦ is defined as follows: According to the above phase consistency model, the moment of the image is calculated, the axis with the minimum moment value is the principal axis, and the axis perpendicular to the principal axis is the maximum moment. Assuming that the principal axis is τ, the specific calculation of the principal axis, the minimum moment and the maximum moment is as follows: τ = 0.5arctan(b / (a-c)) where a and c represent the second-order central moments in the x and y axis parallel directions, respectively, and b represents the mixed second-order central moment between the x and y axes. In addition, the variables t1, t2 and t3 are calculated as follows: In the corner point detection, first, the maximum value of the pixels in the minimum moment mapping is determined, denoted as m1. Then, all points with pixel values greater than 0.01×m1 in the mapping are selected as candidate corner points, i.e. local maximum values, and non-maximum suppression is used to remove redundant candidate points. For edge feature points, an AGAST detector is used to extract edge feature points from the maximum moment mapping. Finally, the extracted corner points and edge points are combined for subsequent feature matching. In constructing the maximum index map, first, a background board with the same size as the original image is created. Then, for each (x, y) pixel point on the background board, the amplitude vector values in N directions are calculated, and the maximum amplitude vector value is designated as the corresponding value of the pixel in the maximum index map, thereby generating the final maximum index description. The feature vector is characterized by using a distribution histogram. For each feature point, a local pixel block of n x n size is selected, and a Gaussian function with a standard deviation equal to n / 2 is used to assign a weight to each pixel. Then, the local region is divided into a sub-grid of N x N, and a histogram distribution containing N components is counted in each sub-grid. Finally, the histograms are connected to generate the final feature vector. By adjusting the order of different directions under the same scale, a plurality of sequences are generated. Based on the sequences, a plurality of different maximum index map descriptions are constructed. By changing the index order, a rotation-invariant maximum index map is obtained. Then, a distribution histogram technique is used to describe the feature vector.

4. The method of claim 1, wherein In step (3), the nearest neighbor and the second nearest neighbor ratio and the region consistency evaluation method are used to preliminarily register the feature points. The improvement strategy is as follows: By calculating the distance ratio of the nearest neighbor and the second nearest neighbor, and setting a threshold for comparison, when the distance ratio is small, it indicates that the difference between the nearest neighbor and the second nearest neighbor is large, and in this case, the matching result is more likely to be correct. For false matching, the distance between the nearest neighbor and the second nearest neighbor is closer, and the ratio is closer to 1. For each point pair obtained from the above matching, first set two neighborhood windows: a large window of 200 x 200 and a small window of 32 x 32. In the large window, only need to judge whether the large window region of the two feature points has intersection, if there is intersection, it is considered that the point pair meets the spatial consistency, otherwise, directly eliminate the matching point pair. On the basis of meeting the spatial consistency, further evaluate the consistency of the region features in the small window, and use the normalized cross correlation (NCC) method for measurement. The method first calculates the pixel average values of the two small window regions, and the specific calculation process is as follows: Where N = 32, A and B represent the pixel average values of the two small window regions, A(i,j) and B(i,j) represent the pixel value size at position coordinates (i,j) in the two window regions, and then the normalized cross correlation is calculated according to the region pixel average value: where σ A and σ B are the standard deviations of the pixel values in the A and B regions, respectively, and are calculated as follows: The NCC value ranges from -1 to 1, 1 indicates that the two regions are completely matched, -1 indicates that the two regions are completely opposite, and 0 indicates that there is no correlation between the two. A higher NCC value indicates a stronger similarity. The best threshold is determined by experiment to be 0.

7. On the premise of ensuring that the large window region meets the spatial consistency, if the NCC value calculated in the small window region is greater than or equal to 0.7, it is considered that the point pair meets the region feature consistency evaluation condition and is retained. Otherwise, the matching point pair is eliminated.

5. The method of claim 1, wherein In step (4), the improved RANSAC algorithm is used to optimize the feature point matching result, and the specific improvement strategy is as follows: (1) First, sort all matching points in descending order according to the NCC value, and apply the selection strategy according to the order. The higher the NCC value, the more likely the point pair is a correct match. The least number of points is selected to estimate the transformation model. Four point pairs are needed to calculate the projection transformation matrix H for homography transformation; (2) Use the estimated transformation model to transform all matching point pairs, and calculate the transformation error; (3) Set a threshold, mark the points less than the threshold as inliers, consider these points as correct matching; otherwise mark them as outliers; (4) Repeat the above steps, each time randomly select a new point pair, count the number of inliers, and use Huber estimation to fine-tune the model in each iteration; (5) Finally, select the model with the most inliers as the required transformation model; In order to ensure the reliability of the final estimated model, the number of iterations n is at least: Where p refers to the probability that all 4 matching point pairs are inliers in at least one sampling, η is the probability of mis-matching point pairs, and L is the minimum value in the estimated transformation model; Finally, for each matching point pair, a specific mathematical mapping relationship is established through the projection transformation. Assuming a matching point pair (x1, y1) and (x2, y2) in two images, the relationship of the matching point pair under the projection transformation can be expressed as follows: Where H is a 3x3 homography matrix, defined as follows: Express it in the form of mathematical calculation: In each iteration of the estimated model, fine-tune the transformation by introducing Huber estimation.

Images