Weak feature point cloud rapid splicing method based on image guidance

Abstract

The invention discloses a weak feature point cloud quick splicing method based on image guidance, which comprises the following steps of: 1) acquiring camera parameters required by subsequent calculation by adopting a Zhang's calibration method, and defining a reference plane as a space reference for subsequent pose calculation; 2) synchronously acquiring point cloud and RGB images, performing image matching through an AKAZE algorithm for 2D images, and then screening matching pairs by adopting an angle consistency + RANSAC algorithm to obtain a better 2D matching relationship; 3) calculating an initial 3D matching relation of adjacent point clouds by combining the results of the steps 1 and 2 to realize coarse registration; and 4) taking the point clouds after applying the result of the step 3) as input, and carrying out iterative optimization by adopting a Point-to-Plane-ICP algorithm to obtain an accurate 3D matching relationship between adjacent point clouds so as to complete splicing. According to the method, the problems of high splicing difficulty and low efficiency and precision of weak 3D feature adjacent point clouds are effectively solved, and the point cloud splicing efficiency is improved while the point cloud splicing precision is considered.

Description

Technical Field

[0001] This invention relates to the field of point cloud data processing technology, and in particular to an image-guided, high-precision, and fast stitching method for point clouds with weak 3D features. Background Technology

[0002] 3D point cloud stitching technology is increasingly used in modern industrial systems centered on intelligent manufacturing, especially in precision machining fields such as automotive parts manufacturing. It is a key technology for achieving high-precision quality inspection, dimensional tolerance measurement, and automated assembly guidance. Among these technologies, monocular structured light 3D imaging systems have become the mainstream solution for acquiring 3D information about object surfaces due to their non-contact, high-efficiency, and high-data-density characteristics.

[0003] Currently, due to the inherent field of view of imaging systems and the geometric occlusion of objects, a single scan can only obtain local and incomplete point cloud data of the target. To construct a complete digital model of the target, it is necessary to scan from multiple perspectives and precisely align the various local point clouds in a unified coordinate system using point cloud stitching technology. The core of this process is to accurately solve the relative pose transformations (i.e., rotation matrices and translation vectors) between the point clouds.

[0004] Because fine registration algorithms require a good initial pose (usually provided by the coarse registration step), current point cloud stitching methods typically follow a "coarse registration-fine registration" strategy. However, existing registration techniques suffer from significant bottlenecks when handling specific industrial objects. One type of method relies on the geometric features of the point cloud itself, such as Fast Point Feature Histograms, for coarse registration. However, for "weak feature" workpieces with large flat or smooth curved surfaces, such as automotive body panels, the discriminative power of the 3D feature descriptors drops sharply, leading to a high matching failure rate, difficulty in providing a reliable initial pose, and significant time consumption when processing large-scale point clouds. Meanwhile, the ICP algorithm and its variants, as standard algorithms for fine registration, highly depend on a good initial pose for convergence and accuracy; otherwise, they are prone to getting trapped in local optima, leading to stitching failure or substandard accuracy. Another type of method uses external equipment such as high-precision turntables and robotic arms to obtain the initial pose. While this ensures accuracy, the solution is costly, inflexible, and difficult to adapt to the changing demands of industrial environments. Therefore, how to provide fast, robust and high-precision initial poses for weak feature point clouds at low cost has become a key challenge restricting the development of industrial 3D measurement technology. Summary of the Invention

[0005] To address the shortcomings of existing technologies, namely the difficulty in initial pose estimation, low success rate and efficiency, and high cost and inflexibility of relying on expensive external equipment when stitching point clouds with weak 3D features, the present invention aims to provide an image-guided fast stitching method for weak feature point clouds, which balances stitching efficiency, robustness and accuracy.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] A fast stitching method for weak feature point clouds based on image guidance, the steps of which are as follows:

[0008] 1) Camera calibration and reference plane establishment: The camera is calibrated to obtain the camera's intrinsic and extrinsic parameters required for subsequent calculations, and a reference plane is defined as the spatial reference for subsequent pose calculations.

[0009] 2) Two-dimensional image matching: Take the first of two adjacent two-dimensional images captured by the camera as the target image and the second as the source image. Perform feature point detection and matching on the source image and the target image to obtain two-dimensional matching point pairs and two-dimensional stitching matrix between the source image and the target image.

[0010] 3) Initial pose calculation: Using the camera intrinsic and extrinsic parameters and reference plane obtained in step 1), the three-dimensional spatial coordinates of the two-dimensional matching point pairs obtained in step 2) are calculated by back projection, resulting in the three-dimensional matching point pairs between the source point cloud corresponding to the source image and the target point cloud corresponding to the target image. Based on the three-dimensional matching point pairs and the two-dimensional stitching matrix obtained in step 2), the initial pose transformation relationship between the source point cloud and the target point cloud is calculated. The source point cloud is transformed using the initial pose transformation relationship to obtain the coarsely registered source point cloud.

[0011] 4) Fine Registration: Using the source point cloud obtained in step 3) after coarse registration as the initial value for iterative optimization, the Point-to-Plane Iterative Closest Point (Point-to-Plane-ICP) algorithm is used to fine register the source and target point clouds, resulting in the rotation matrix R between the coarsely registered source and target point clouds. fine Translation vector t fine It is then applied to the source point cloud for coarse registration, thus completing the point cloud stitching.

[0012] Furthermore, the camera calibration and reference plane establishment in step 1) are obtained as follows:

[0013] 1.1) Establish the coordinate correspondence of the calibration plate images: Take multiple images of the calibration plate printed with a checkerboard pattern, and define the plane where the calibration plate is located as the Z-axis of the three-dimensional spatial coordinate system. w=0 plane, the three-dimensional coordinates of each corner point P w =[X w ,Y w ,0] T Given an image of a calibration board, the two-dimensional pixel coordinates p = [u, v] of each corner point in the image are obtained using a corner detection algorithm. T ;

[0014] 1.2) Using the camera projection model shown in equation (1), establish the three-dimensional spatial coordinates of the calibration plate corner point P. w Mapping to two-dimensional pixel coordinates p; where K is the camera intrinsic parameter matrix, [R calib |t calib [ ] is a 3×4 camera extrinsic parameter matrix, and s is a scale factor representing the depth variable in the camera coordinate system:

[0015]

[0016] 1.3) Solving for camera parameters and determining the reference plane: Using the calibration board image as input, through a nonlinear optimization process, the camera intrinsic parameter matrix K in equation (1) and the extrinsic parameter R corresponding to each calibration board image are solved together. calib,j ,t calib,j From the multiple sets of extrinsic parameters obtained, select one set of extrinsic parameters {R}. ref ,t ref} as a benchmark, where R ref The third column is the normal vector n of the plane. ref n ref and t ref They jointly define a three-dimensional reference plane.

[0017] Furthermore, the two-dimensional image matching described in step 2) is performed as follows:

[0018] 2.1) Constructing the Nonlinear Scale Space: The scale space of the image is constructed by solving the nonlinear diffusion equation. This process is performed through fast explicit diffusion, and the equation is as follows:

[0019]

[0020] Where L is the image brightness, t is the scale parameter, div is the divergence operator, and g is the conductivity function, the value of which varies with the image gradient. It adapts to changes in noise to maintain sharp edges while smoothing out noise;

[0021] 2.2) Feature point detection: In the constructed nonlinear scale space, feature points are detected by calculating the determinant response of the scale-normalized Hessian matrix, as shown in equation (3):

[0022]

[0023] Where L xx L yy L xy Let σ be the second partial derivative of the image at point (x, y) with respect to the corresponding coordinates, and let σ be the scale of the feature point. Search for local maxima of equation (3) in the scale space as candidate feature points.

[0024] 2.3) Binary descriptor generation: Assign a principal direction to each candidate feature point to achieve rotation invariance; then, use local binary mode to generate feature descriptors, which generate a binary descriptor that is robust to rotation, scale and illumination changes by performing grid sampling and brightness or gradient comparison in the neighborhood of the feature point.

[0025] 2.4) Initial matching pair calculation: For a candidate feature point D in the source image... src Find the binary descriptor D with the closest Hamming distance in the target image. tgt The two constitute a matching pair; thus, a set of matching pairs is obtained. Where p i =(u i ,v i ) represents the coordinates of a feature point in the source image, p' i =(u' i ,v' i ) are the coordinates of the corresponding feature points in the target image; the matching pair set S is the two-dimensional matching point pair.

[0026] Furthermore, following step 2.4), the following steps are also included:

[0027] 2.5) First, perform angle consistency screening. For each matching pair (p... i ,p' i ), calculate its displacement vector on the image plane. and displacement angle θ i :

[0028]

[0029] When matching the displacement angle θ i With the principal direction θ median The deviation is less than the threshold T angle Points are considered to be interior points with consistent angles, and the filtered matching subset is S. angle As shown in equation (5-1):

[0030] S angle ={(p i ,p′ i )∈S||θi -θ median | <T angle} (5-1)

[0031] 2.6) The transformation model shown in equation (6) is used to describe the geometric mapping relationship between the source image and the target image:

[0032]

[0033] Where p = (u, v) is a point on the source image, and p' = (u', v') is its corresponding point on the target image; the rotation angle α and the translation vector t = [t] are solved by equation (6). x ,t y ] T These parameters together constitute a two-dimensional stitching matrix that describes the stitching relationship between the two images.

[0034] The Random Sample Consensus (RANSAC) algorithm iteratively selects matching points S from an angle-filtered set of points. angle In this process, two non-collinear matching pairs are randomly selected without replacement, and a system of linear equations is established and solved to obtain a candidate transformation matrix M. candidate ; this M candidate Matrix applied to S angle For all matching point pairs, calculate their reprojection error, which is defined as the reprojection error of the true observation point p in the target image. i 'Compared to point p in the source image i via M candidate The predicted position p obtained after transformation i The Euclidean distance between them is shown in equation (7):

[0035]

[0036] If err i Less than the reprojection error threshold T reproj Matching pairs (p) i ,p' i ) is considered to be related to the current M candidate Consistent interior points form an interior point set; throughout all iterations, record the M that yields the largest interior point set. candidate As the transformation matrix, the corresponding maximum set of interior points is the final two-dimensional matching point pair;

[0037] Based on the final two-dimensional matching point pairs, the centroids of the feature points in the source image are calculated according to equation (8). centroid of feature points in the target image

[0038]

[0039] Where p src,i p is the i-th matching point on the source image. tgt,i Let i be the i-th matching point on the target image.

[0040] Furthermore, the threshold T angle Based on dynamic calculation using the Median Absolute Deviation (MAD), the median is first used to estimate the principal direction θ of the scene motion. median Then, after calculating the absolute deviation of all data from the median, the median deviation is calculated. The dynamic calculation strategy is: when the MAD of the data is small, a set θ is used. min The minimum threshold is used to avoid overly strict thresholding; when MAD is large, C·MAD is used to relax the threshold to retain more potential interior points, where C is the scaling factor and the threshold T is the minimum threshold. angle The specific calculation formula is shown in (5-2):

[0041]

[0042] Furthermore, the initial pose calculation described in step 3) is specifically performed as follows:

[0043] 3.1) Two-dimensional point back projection: Using the camera intrinsic parameter matrix K obtained in step 1), the two centroids output in step 2) are projected back. and By back-projecting from the pixel coordinate system onto the normalized image plane in the camera coordinate system, the corresponding three-dimensional direction vector q is obtained. src and q tgt :

[0044]

[0045] 3.2) Intersection of the ray and the reference plane: Using the reference extrinsic parameter {R} determined in step 1). ref ,t ref The three-dimensional reference plane defined by} is used to calculate the distance from the camera optical center along q. src and q tgt The intersection point P of the ray in the direction and the three-dimensional reference plane. src_c and P tgt_c For a known direction vector q, the depth λ of its intersection point with the three-dimensional reference plane is calculated as follows:

[0046]

[0047] Where n ref The normal vector of the reference plane in step 1) is P; the three-dimensional coordinates of the intersection point are P. c=λ·q; Using this method, P in the camera coordinate system can be obtained. src_c and P tgt_c Then through P w =R ref ·(P c -t ref Transform it back to world coordinates to get P src_w and P tgt_w The vector difference between the two is the estimated value of the initial translation vector in three-dimensional space.

[0048] 3.3) Construct the initial three-dimensional rotation matrix R initial Using the two-dimensional rotation angle α in the two-dimensional splicing matrix obtained in step 2), the initial three-dimensional rotation matrix R is constructed according to equation (11). initial ;

[0049]

[0050] Source point cloud via R initial , Registration is performed, which yields the coarsely registered source point cloud.

[0051] Furthermore, the Point-to-Plane-ICP algorithm described in step 4) takes the source point cloud obtained from the coarse registration in step 3) as input and aims to minimize the source point cloud P after coarse registration. S To the target point cloud P T The registration error function E(R,t) is defined as the sum of squared orthogonal distances from the transformed points in the source point cloud to the tangent planes of their corresponding target points, as shown in equation (12):

[0052]

[0053] Where, p s，i ∈P S p t，i ∈P T For p s，i The corresponding point, n t，i For point p t，i The unit normal vector at point R, and R and t are the rotation matrix and translation vector to be optimized, respectively; R fine t fine This is the final optimized result.

[0054] Compared with the prior art, the present invention has the following beneficial effects:

[0055] 1. Cross-domain guidance enhances robustness and efficiency: This invention transforms the problem of coarse registration of weak feature point clouds, which is difficult to solve directly in 3D, into the more feature-rich 2D image domain. By utilizing the inherent stable image texture of the object surface, it effectively avoids the problems of sparse and unstable 3D geometric features, thereby quickly generating high-precision initial poses and greatly improving the success rate and robustness of stitching.

[0056] 2. High-quality initial pose, ensuring global optimum: Through high-precision camera calibration and rigorous 2D feature matching, this method provides an initial pose that is very close to the true value. This high-quality initial value ensures that the subsequent Point-to-Plane-ICP algorithm can quickly converge to the global optimum, avoiding the risk of traditional methods getting stuck in local optima due to poor initial position, thus significantly improving the final stitching accuracy.

[0057] 3. Algorithm Collaboration, Balancing Efficiency and Accuracy: This invention effectively combines the robustness of the Accelerated-KAZE (AKAZE) algorithm, the directness of solving the initial pose problem of point clouds, and the accuracy of the Point-to-Plane-ICP algorithm, forming a highly efficient "coarse-to-fine" stitching process. The coarse registration stage quickly locks the approximate position, while the fine registration stage performs fine adjustments, balancing the algorithm's operational efficiency with the final stitching accuracy.

[0058] This invention can control the stitching speed of multiple groups of point clouds with a total of tens of millions of points within 20 seconds without using additional auxiliary equipment such as robotic arms. Attached Figure Description

[0059] Figure 1 This is a flowchart of the image-guided rapid stitching method for weak feature point clouds according to the present invention.

[0060] Figure 2 This is a schematic diagram of the calibration board used in the camera calibration part of the stitching method of the present invention.

[0061] Figure 3 This is a 2D feature matching result diagram of the image feature matching part in a certain embodiment of the stitching method of the present invention.

[0062] Figure 4 This is the initial pose of the point cloud obtained by solving a certain embodiment of the stitching method of the present invention, i.e., the coarse registration result image.

[0063] Figure 5 This is a stitching result image of point clouds after fine registration and global optimization in a certain embodiment of the stitching method of the present invention. Detailed Implementation

[0064] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that the following embodiments are merely preferred examples of this invention, intended to explain the core technology of this invention, and not to constitute any limitation on the scope of protection of this invention.

[0065] This invention proposes an image-guided fast stitching method for weak feature point clouds. The core technology lies in transforming the initial pose problem, which is difficult to solve in 3D space due to sparse and unstable geometric features, into a more information-rich 2D image domain through an image-guided strategy. This provides a high-quality initial transformation matrix for subsequent 3D fine registration. This method effectively solves the problems of difficult initial pose estimation, low efficiency, and low success rate caused by the lack of geometric constraints when stitching weak 3D feature point clouds.

[0066] The specific implementation steps of this invention are detailed below, and please also refer to... Figure 1 The flowchart shown is as follows.

[0067] 1) Camera Calibration and Reference Plane Establishment: The camera is calibrated to obtain the intrinsic and extrinsic parameters required for subsequent calculations, and a reference plane is defined as the spatial benchmark for subsequent pose calculation. By processing the checkerboard calibration board image, an optimization problem with reprojection error E as the objective function is established and solved to accurately solve for the camera's intrinsic and extrinsic parameters and distortion coefficients containing higher-order terms, establishing a precise mathematical mapping model from two-dimensional pixel coordinates to three-dimensional spatial coordinates. This step provides the necessary mathematical model and parameter benchmark for subsequent three-dimensional pose calculation. Its specific implementation relies on Zhang's calibration method, and the specific process is as follows:

[0068] 1.1) Establish the coordinate correspondence between the calibration plate images: Take multiple images of the calibration plate printed with a checkerboard pattern (e.g., ...). Figure 2 As shown), the plane where the calibration plate is located is defined as the Z-axis of the three-dimensional spatial coordinate system. w =0 plane, the three-dimensional coordinates of each corner point P w =[X w ,Y w ,0] T Given an image of a calibration board, the two-dimensional pixel coordinates p = [u, v] of each corner point in the image are obtained using a corner detection algorithm. T ;

[0069] 1.2) Using the camera projection model shown in equation (1), establish the three-dimensional spatial coordinates of the calibration plate corner point P. w Mapping to two-dimensional pixel coordinates p; where K is the camera intrinsic parameter matrix, [R calib |t calib[ ] is a 3×4 camera extrinsic parameter matrix, and s is a scale factor representing the depth variable in the camera coordinate system:

[0070]

[0071] 1.3) Solving for camera parameters and determining the reference plane: Using the calibration board image as input, through a nonlinear optimization process, the camera intrinsic parameter matrix K in equation (1) and the extrinsic parameter R corresponding to each calibration board image are solved together. calib,j ,t calib,j From the multiple sets of extrinsic parameters obtained, select one set of extrinsic parameters {R}. ref ,t ref} as a benchmark, where R ref The third column is the normal vector n of the plane. ref n ref and t ref A three-dimensional reference plane is defined together, which serves as a fixed spatial reference in the subsequent calculation in step 3.

[0072] 2) Two-dimensional image matching: For RGB images acquired from two adjacent scans (the first image as the target image and the second as the source image), the AKAZE algorithm, which is highly robust to scale and rotation changes, is used for two-dimensional image feature matching and matching pair selection. This step aims to obtain two-dimensional matching point pairs and a two-dimensional stitching matrix (e.g., ...) between two images from adjacent viewpoints. Figure 3 (As shown). The specific implementation process is as follows:

[0073] 2.1) Constructing the Nonlinear Scale Space: The scale space of the image is constructed by solving the nonlinear diffusion equation. This process is performed through fast explicit diffusion, and the equation is as follows:

[0074]

[0075] Where L is the image brightness, t is the scale parameter, div is the divergence operator, and g is the conductivity function, the value of which varies with the image gradient. It adapts to changes in noise to maintain sharp edges while smoothing out noise;

[0076] 2.2) Feature point detection: In the constructed nonlinear scale space, feature points are detected by calculating the determinant response of the scale-normalized Hessian matrix, as shown in equation (3):

[0077]

[0078] Where L xx L yy L xyLet σ be the second partial derivative of the image at point (x, y) with respect to the corresponding coordinates, and let σ be the scale of the feature point. Search for local maxima of equation (3) in the scale space as candidate feature points.

[0079] 2.3) Binary descriptor generation: Assign a principal direction to each candidate feature point to achieve rotation invariance; then, use local binary mode to generate feature descriptors, which generate a binary descriptor that is robust to rotation, scale and illumination changes by performing grid sampling and brightness or gradient comparison in the neighborhood of the feature point.

[0080] 2.4) Initial matching pair calculation: Hamming distance d H (D1, D2) is defined as the difference between corresponding bits of two binary descriptors D1 and D2. The smaller the distance (difference), the more similar the two feature points are. For a candidate feature point in the source image, the descriptor D1... src Find the descriptor D with the closest Hamming distance in the target image. tgt These two constitute an initial matching pair;

[0081] 2.5) Two-stage screening strategy and matrix solution: After initial matching, a set of matching pairs is obtained. Where p i =(u i ,v i ) represents the coordinates of a feature point in the source image, p' i =(u' i ,v' i (p) represents the coordinates of the corresponding feature points in the target image. First, angle consistency filtering is performed; for each matching pair (p)... i ,p' i ), calculate its displacement vector on the image plane. and its displacement angle θ i :

[0082]

[0083] All matching pairs of displacement angles form a set Furthermore, to avoid the influence of extreme values, a more robust median is used to estimate the principal direction θ of scene motion. median During the screening process, the deviation of each angle from the main direction is calculated. When the deviation of its displacement angle from the main direction is less than the threshold T, the result is considered a screening result. angle The inliers are considered to be "angle-consistent", and the filtered matching subset is S. angle :

[0084]

[0085] The threshold T angle The selection strategy is made more adaptable by dynamically calculating MAD. MAD is a statistic used to measure the dispersion of data. The calculation method is as follows: first, calculate the absolute deviation of all data from the median, and then calculate the median deviation. The dynamic calculation strategy is: when the MAD of the data is small, a set θ is used. min The minimum threshold is used to avoid overly strict thresholding. When MAD is large, C·MAD is used to relax the threshold to retain more potential inliers, where C is a scaling factor. When the data follows a normal distribution, C·MAD (where C=3) is approximately two standard deviations, which theoretically can cover about 95% of inliers, further enhancing the robustness of the algorithm. Through this stage, mismatches with directional deviations can be effectively eliminated.

[0086] 2.6) To further filter matching pairs and accurately solve the transformation matrix between images, this scheme uses the RANSAC algorithm for the final filtering and stitching matrix solution. Considering that when the camera acquires images at adjacent positions, its relative pose change is mainly manifested as rigid body motion (translation and rotation), the similarity transformation model shown in equation (6) is selected to describe the geometric mapping relationship between the source image and the target image:

[0087]

[0088] Where p = (u, v) is a point on the source image, and p' = (u', v') is its corresponding point on the target image; the rotation angle α and the translation vector t = [t] are solved by equation (6). x ,t y ] T These parameters together constitute a two-dimensional stitching matrix that describes the stitching relationship between the two images.

[0089] The RANSAC algorithm iteratively selects matching points S from the angle-filtered set of points. angle In this process, two non-collinear matching pairs are randomly selected without replacement, and a system of linear equations is established and solved to obtain a candidate transformation matrix M. candidate . This M candidate Matrix applied to S angle For all matching point pairs, calculate their reprojection error, which is defined as the reprojection error of the true observation point p in the target image. i 'Compared to point p in the source image i via M candidate The predicted position p obtained after transformation i The Euclidean distance between them is shown in equation (7):

[0090]

[0091] If err i Less than the reprojection error threshold T reproj Matching pairs (p) i ,p' i ) is considered to be related to the current M candidate Consistent interior points form an interior point set. During all iterations, record the M that yields the largest interior point set. candidate As the transformation matrix (containing α), the corresponding maximum set of interior points contains all the finally selected matching pairs (source image p). src,i Target image p tgt,i ).

[0092] Based on the final high-precision two-dimensional matching point pairs, the centroid of feature points in the source image is calculated. centroid of feature points in the target image

[0093]

[0094] 3) Initial Pose Calculation: This step utilizes the camera parameters and reference plane from step 1) to calculate the 3D spatial coordinates of the 2D matching point pairs obtained in step 2) through back projection and other steps. This yields the 3D matching point pairs between the source point cloud corresponding to the source image and the target point cloud corresponding to the target image. Using the 3D matching point pairs and the 2D stitching matrix information, the initial pose transformation relationship between the source and target point clouds is calculated. The source point cloud is then transformed using this initial pose transformation relationship to obtain a coarsely registered source point cloud (e.g., ...). Figure 4 As shown in the image, the coarse registration process is completed. The specific implementation process is as follows:

[0095] 3.1) Two-dimensional point back projection: Since only the centroid pixel coordinates (2 known variables) are known, while s in Equation 1 and its corresponding 3D coordinates have 4 unknown variables, it is necessary to use a reference plane and back projection to solve for the 3D spatial points. Using the camera intrinsic parameter matrix K obtained in step 1), the two centroids output in step 2) are... and By back-projecting from the pixel coordinate system onto the normalized image plane (s=1) in the camera coordinate system, the corresponding three-dimensional direction vector q is obtained. src and q tgt :

[0096]

[0097] 3.2) Intersection of the ray and the reference plane: Using the reference extrinsic parameter {R} determined in step 1). ref ,t ref The defined three-dimensional reference plane is used to calculate the distance from the camera optical center (origin) along q. src and q tgtThe intersection point P of the ray in the direction and the three-dimensional reference plane. src_c and P tgt_c For a known direction vector q, the depth λ of its intersection point with the reference plane is calculated as follows:

[0098]

[0099] Where n ref The normal vector of the reference plane in step 1) is P; the three-dimensional coordinates of the intersection point are P. c =λ·q; Using this method, P in the camera coordinate system can be obtained. src_c and P tgt_c Then through P w =R ref ·(P c -t ref Transform it back to world coordinates to get P src_w and P tgt_w The vector difference between the two is the estimated value of the initial translation vector in three-dimensional space.

[0100] 3.3) Construct the initial three-dimensional rotation matrix R initial Using the two-dimensional rotation angle α in the two-dimensional stitching matrix calculated in step 2), an initial three-dimensional rotation matrix R is constructed. initial Rotation in the image plane corresponds to rotation about its principal optical axis (Z-axis) in the 3D camera coordinate system. Therefore, the initial rotation matrix is:

[0101]

[0102] Source point cloud via R initial , Registration is performed, which yields the coarsely registered source point cloud.

[0103] 4) Fine registration based on the Point-to-Plane-ICP algorithm: To further improve the stitching accuracy to meet industrial application requirements, this step performs fine registration based on the coarse registration results (e.g., ...). Figure 5 (As shown). The point cloud data after coarse registration in step 3 is used as the initial input for the Point-to-Plane-ICP algorithm. Through iterative optimization, the alignment relationship between the source point cloud and the target point cloud is further refined, and finally a high-precision rotation matrix R is obtained. fine Translation vector t fine The specific implementation process is as follows:

[0104] 4.1) Iterative optimization: The Point-to-Plane-ICP algorithm is adopted. Its core lies in minimizing the source point cloud P after coarse registration. S To the target point cloud P TThe registration error function E(R,t) is defined as the sum of squared orthogonal distances from the transformed points in the source point cloud to the tangent planes of their corresponding target points, as shown in equation (12):

[0105]

[0106] Where, p s，i ∈P S p t，i ∈P T For p s，i The corresponding point, n t，i For point p t，i The unit normal vector at point R is given, and R and t are the rotation matrix and translation vector to be optimized, respectively.

[0107] 4.2) Convergence Criterion: The algorithm converges when the increment of the transformation parameter is less than a preset threshold or when the maximum number of iterations is reached. The final transformation (R) is obtained. fine ,t fine The final optimized result is obtained by combining and applying it to the source point cloud of coarse registration in step 3, thus obtaining the final high-precision stitching result.

[0108] Repeat steps 2) to 4) to align and merge all local point cloud data into a unified global coordinate system. This will stitch all local point cloud data together to complete a complete and accurate 3D reconstruction of the weak 3D feature point cloud.

[0109] To verify the technical effectiveness of the proposed rapid point cloud stitching method for weak features, the experimental subject was a car fender, a common material in the automotive manufacturing industry, characterized by its weak feature surface. This fender is L-shaped, with a small number of marker points affixed to one side (the marker point side) to verify registration accuracy, while the other side (the unmarked side) is unmarked. A fixed camera was used, and the fender was moved to simulate an industrial assembly line scanning scenario, continuously acquiring nine sets of point cloud data to cover one complete side of the fender. Detailed data information is shown in Table 1. This invention does not require prior model training and can directly process the input point cloud data.

[0110] Table 1. Experimental Data (Point Cloud Points in Ten Thousands)

[0111]

[0112] The experimental environment consisted of a Windows 11 operating system, an Intel Core i7-10750H CPU, an NVIDIA GeForce RTX 2060 GPU, and 16GB of RAM; the software environment was Python 3.9.

[0113] To verify the effectiveness of the method of this invention, five mainstream or representative point cloud registration methods were selected for comparison. Since the selected methods primarily target the registration of two sets of point clouds, data 4 and data 5 from the "Identifier Point Side" section of Table 1 were used as input for pairwise registration tests. Each method was repeated five times, and the average value was taken as the final experimental index. A brief description of the selected comparison methods is as follows:

[0114] Point-to-Plane ICP algorithm: Iteratively optimizes by minimizing the distance between the source and target points on the plane, and is one of the benchmark methods for fine registration; Principal Component Analysis (PCA): A fast coarse registration method. It achieves initial alignment by calculating the covariance matrix of the point cloud and aligning its principal directions (eigenvectors). This method is extremely fast, but sensitive to the symmetry and overlap rate of the point cloud; Generalized Iterative Process (GICP): A further generalization of the ICP algorithm, it considers the local planar structure (covariance) of both the source and target point clouds in error calculation, and has higher robustness to noise and outliers, but the computational cost is correspondingly increased; Fast Global Registration (FGR): A global registration algorithm that can quickly find a better initial transformation matrix by constructing a correspondence on the downsampled point cloud and optimizing a global objective function; Coherent Point Shift (CPD): A non-rigid registration method based on probability theory, which aligns the target point cloud with the target point cloud by maximizing the likelihood function.

[0115] The experiment used the following three core metrics to quantitatively evaluate the performance of all methods:

[0116] 1. Average deviation of marker points (mm): After registration of the "marker point side" data, the Euclidean distance between the corresponding marker points in the two sets of point clouds is measured. This value can intuitively and accurately reflect the results; the smaller the value, the higher the accuracy.

[0117] 2. Root Mean Square Error (RMSE) Average (mm): Used to evaluate the geometric alignment between the registered source and target point clouds. RMSE calculates the root mean square distance between a point in the transformed source point cloud and its nearest neighbor in the target point cloud. The smaller the value, the better the overall geometric alignment. The calculation formula is as follows:

[0118]

[0119] Where N is the number of valid corresponding point pairs, p i It is a point in the source point cloud, q i It is its corresponding point in the target point cloud, and ||·|| represents the Euclidean distance between the two points.

[0120] 3. Average total algorithm time (s): Records the total time required from input point cloud data to output of the final registration result. This metric directly reflects the computational efficiency of the algorithm.

[0121] The five existing methods mentioned above were compared with the method of this invention, and the experimental results are shown in Table 2:

[0122] Table 2. Comparative Experiment Results (Data processed are data from marker points 4 and 5)

[0123]

[0124] Analysis of the experimental results shows that most methods simply overlap adjacent point clouds as much as possible to bring the distance between two point clouds closer, while ignoring that point cloud stitching only requires overlapping the common parts of the two point clouds. Therefore, in terms of stitching accuracy (marker point deviation), the average marker point deviation of this invention is only 1.16mm, significantly better than all the comparative methods. The second best performing CPD method has a deviation of 18.76mm, while the deviations of the other methods are all above 20mm. This indicates that this invention can achieve high-fidelity point cloud stitching and effectively solves the problem of inaccurate stitching of weak feature surfaces.

[0125] Regarding alignment quality (RMSE): the average RMSE value of this invention is 0.2644 mm, which is on par with the best-performing PCA method (0.2599 mm) and significantly better than methods such as ICP, GICP, FGR, and CPD. This demonstrates that this invention achieves high precision while also ensuring tightness and high quality in local point cloud alignment, without sacrificing accuracy for speed.

[0126] In terms of computational efficiency (total algorithm time): the average total time for processing the two sets of point clouds using the method of this invention is only 2.02 seconds. Compared with the fastest comparison method PCA (8.23 seconds), the efficiency of this invention is improved by about 4 times. Compared with other methods with higher accuracy such as ICP and GICP, the efficiency is improved by more than 6 to 12 times.

[0127] Therefore, it can be seen that the present invention outperforms the comparative methods in the field of weak feature point cloud stitching.

[0128] In summary, the technical solution proposed in this invention transforms the three-dimensional registration problem into a two-dimensional image domain solution, and combines high-precision camera calibration and initial pose solving to provide reliable initial conditions for subsequent fine registration. This forms a logically rigorous, step-by-step optimization process that uses images to guide coarse registration and ICP to optimize fine registration. As a result, the robustness and efficiency of stitching weak three-dimensional feature point clouds are significantly improved while ensuring high accuracy.

[0129] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for rapid stitching of weak feature point clouds based on image guidance, characterized in that: The steps are as follows: 1) Camera calibration and reference plane establishment: The camera is calibrated to obtain the camera's intrinsic and extrinsic parameters required for subsequent calculations, and a reference plane is defined as the spatial reference for subsequent pose calculations. 2) Two-dimensional image matching: Take the first of two adjacent two-dimensional images captured by the camera as the target image and the second as the source image. Perform feature point detection and matching on the source image and the target image to obtain two-dimensional matching point pairs and two-dimensional stitching matrix between the source image and the target image. 3) Initial pose calculation: Using the camera intrinsic and extrinsic parameters and reference plane obtained in step 1), the three-dimensional spatial coordinates of the two-dimensional matching point pairs obtained in step 2) are calculated by back projection, resulting in the three-dimensional matching point pairs between the source point cloud corresponding to the source image and the target point cloud corresponding to the target image. Based on the three-dimensional matching point pairs and the two-dimensional stitching matrix obtained in step 2), the initial pose transformation relationship between the source point cloud and the target point cloud is calculated. The source point cloud is transformed using the initial pose transformation relationship to obtain the coarsely registered source point cloud. 4) Fine registration: Using the source point cloud obtained in step 3) after coarse registration as the initial value for iterative optimization, the point-to-surface iterative nearest point algorithm is used to fine register the source and target point clouds, resulting in the rotation matrix R between the source and target point clouds. fine Translation vector t fine It is then applied to the source point cloud for coarse registration, thus completing the point cloud stitching.

2. The image-guided rapid stitching method for weak feature point clouds according to claim 1, characterized in that, The camera calibration and reference plane establishment in step 1) are obtained as follows: 1.1) Establish the coordinate correspondence of the calibration plate images: Take multiple images of the calibration plate printed with a checkerboard pattern, and define the plane where the calibration plate is located as the Z-axis of the three-dimensional spatial coordinate system. w =0 plane, the three-dimensional coordinates of each corner point P w =[X w ,Y w ,0] T Given an image of a calibration board, the two-dimensional pixel coordinates p = [u, v] of each corner point in the image are obtained using a corner detection algorithm. T ; 1.2) Using the camera projection model shown in equation (1), establish the three-dimensional spatial coordinates of the calibration plate corner point P. w Mapping to two-dimensional pixel coordinates p; where K is the camera intrinsic parameter matrix, [R calib |t calib [ ] is a 3×4 camera extrinsic parameter matrix, and s is a scale factor representing the depth variable in the camera coordinate system: 1.3) Solving for camera parameters and determining the reference plane: Using the calibration board image as input, through a nonlinear optimization process, the camera intrinsic parameter matrix K in equation (1) and the extrinsic parameter R corresponding to each calibration board image are solved together. calib,j ,t calib,j From the multiple sets of extrinsic parameters obtained, select one set of extrinsic parameters {R}. ref ,t ref } as a benchmark, where R ref The third column is the normal vector n of the plane. ref n ref and t ref They jointly define a three-dimensional reference plane.

3. The image-guided rapid stitching method for weak feature point clouds according to claim 1, characterized in that, The two-dimensional image matching described in step 2) is performed as follows: 2.1) Constructing the Nonlinear Scale Space: The scale space of the image is constructed by solving the nonlinear diffusion equation. This process is performed through fast explicit diffusion, and the equation is as follows: Where L is the image brightness, t is the scale parameter, div is the divergence operator, and g is the conductivity function, the value of which varies with the image gradient. It adapts to changes in noise to maintain sharp edges while smoothing out noise; 2.2) Feature point detection: In the constructed nonlinear scale space, feature points are detected by calculating the determinant response of the scale-normalized Hessian matrix, as shown in equation (3): Where L xx L yy L xy Let σ be the second partial derivative of the image at point (x, y) with respect to the corresponding coordinates, and let σ be the scale of the feature point. Search for local maxima of equation (3) in the scale space as candidate feature points. 2.3) Binary descriptor generation: Assign a principal direction to each candidate feature point to achieve rotation invariance; Subsequently, a feature descriptor is generated using a local binary mode. This feature descriptor generates a binary descriptor that is robust to changes in rotation, scale, and illumination by performing mesh sampling and brightness or gradient comparison within the neighborhood of the feature point. 2.4) Initial matching pair calculation: For a candidate feature point D in the source image... src Find the binary descriptor D with the closest Hamming distance in the target image. tgt The two constitute a matching pair; thus, a set of matching pairs is obtained. Where p i =(u i ,v i ) represents the coordinates of a feature point in the source image, p' i =(u' i ,v' i ) are the coordinates of the corresponding feature points in the target image; the matching pair set S is the two-dimensional matching point pair.

4. The image-guided rapid stitching method for weak feature point clouds according to claim 3, characterized in that, After step 2.4), also Includes the following steps, 2.5) First, perform angle consistency screening. For each matching pair (p... i ,p' i ), calculate its displacement vector on the image plane. and displacement angle θ i : When matching the displacement angle θ i With the principal direction θ median The deviation is less than the threshold T angle Points are considered to be interior points with consistent angles, and the filtered matching subset is S. angle As shown in equation (5-1): S angle ={(p i ,p′ i )∈S||θ i -θ median | <T abgle } (5-1) 2.6) The transformation model shown in equation (6) is used to describe the geometric mapping relationship between the source image and the target image: Where p = (u, v) is a point on the source image, and p' = (u', v') is its corresponding point on the target image; the rotation angle α and the translation vector t = [t] are solved by equation (6). x ,t y ] T These parameters together constitute a two-dimensional stitching matrix that describes the stitching relationship between the two images. The random sampling consensus algorithm iteratively selects from the set of matching points S that have been filtered by angles. angle In this process, two non-collinear matching pairs are randomly selected without replacement, and a system of linear equations is established and solved to obtain a candidate transformation matrix M. candidate ; this M candidate Matrix applied to S angle For all matching point pairs, calculate their reprojection error, which is defined as the reprojection error of the true observation point p in the target image. i 'Compared to point p in the source image i via M candidate The predicted position p obtained after transformation i The Euclidean distance between them is shown in equation (7): If err i Less than the reprojection error threshold T reproj Matching pairs (p) i ,p' i ) is considered to be related to the current M candidate Consistent interior points form an interior point set; throughout all iterations, record the M that yields the largest interior point set. candidate As the transformation matrix, the corresponding maximum set of interior points is the final two-dimensional matching point pair; Based on the final two-dimensional matching point pairs, the centroids of the feature points in the source image are calculated according to equation (8). centroid of feature points in the target image Where p src,i p is the i-th matching point on the source image. tgt,i Let i be the i-th matching point on the target image.

5. The image-guided rapid stitching method for weak feature point clouds according to claim 4, characterized in that, The threshold T angle Based on dynamic calculation using the median absolute deviation, the median is first used to estimate the principal direction θ of the scene motion. median Then, after calculating the absolute deviation of all data from the median, the median deviation is calculated. The dynamic calculation strategy is: when the MAD of the data is small, a set θ is used. min The minimum threshold is used to avoid overly strict thresholding; when MAD is large, C·MAD is used to relax the threshold to retain more potential interior points, where C is the scaling factor and the threshold T is the minimum threshold. angle The specific calculation formula is shown in (5-2):

6. The image-guided rapid stitching method for weak feature point clouds according to claim 4, characterized in that, Step 3) involves the initial pose calculation, and the specific steps are as follows: 3.1) Two-dimensional point back projection: Using the camera intrinsic parameter matrix K obtained in step 1), the two centroids output in step 2) are projected back. and By back-projecting from the pixel coordinate system onto the normalized image plane in the camera coordinate system, the corresponding three-dimensional direction vector q is obtained. src and q tgt : 3.2) Intersection of the ray and the reference plane: Using the reference extrinsic parameter {R} determined in step 1). ref ,t ref The three-dimensional reference plane defined by} is used to calculate the distance from the camera optical center along q. src and q tgt The intersection point P of the ray in the direction and the three-dimensional reference plane. src_c and P tgt_c For a known direction vector q, the depth λ of its intersection point with the three-dimensional reference plane is calculated as follows: Where n ref The normal vector of the reference plane in step 1) is P; the three-dimensional coordinates of the intersection point are P. c =λ·q; Using this method, P in the camera coordinate system can be obtained. src_c and P tgt_c Then through P w =R ref ·(P c -t ref Transform it back to world coordinates to get P src_w and P tgt_w The vector difference between the two is the estimated value of the initial translation vector in three-dimensional space. 3.3) Construct the initial three-dimensional rotation matrix R initial Using the two-dimensional rotation angle α in the two-dimensional splicing matrix obtained in step 2), the initial three-dimensional rotation matrix R is constructed according to equation (11). initial ; Source point cloud via R initial , Registration is performed, which yields the coarsely registered source point cloud.

7. The image-guided rapid stitching method for weak feature point clouds according to claim 1, characterized in that, The Point-to-Plane-ICP algorithm described in step 4) takes the source point cloud obtained from the coarse registration in step 3) as input and aims to minimize the source point cloud P after coarse registration. S To the target point cloud P T The registration error function E(R,t) is defined as the sum of squared orthogonal distances from the transformed points in the source point cloud to the tangent planes of their corresponding target points, as shown in equation (12): Where, p s，i ∈P S p t，i ∈P T For p s，i The corresponding point, n t，i For point p t，i The unit normal vector at point R, and R and t are the rotation matrix and translation vector to be optimized, respectively; R fine t fine This is the final optimized result.

Images