Multi-view point cloud splicing method based on local feature extraction

Abstract

The invention discloses a multi-view point cloud splicing method based on local feature extraction, and the method comprises the steps: obtaining three-dimensional point cloud data in a to-be-spliced scene, and enabling a certain overlapping region to be needed between two point clouds; performing data processing on the acquired three-dimensional point cloud data; presetting a mark point in an overlapping region of the point clouds at different visual angles; key points are extracted around a preset mark point according to a local feature extraction algorithm, and matching point pairs are searched by using a fast point feature histogram descriptor; iteratively selecting an optimal matching point pair to estimate a rotation matrix between the to-be-spliced point clouds; and carrying out point cloud fusion in combination with an iterative nearest point algorithm to realize multi-view point cloud splicing. According to the invention, point cloud splicing of a low-overlapping-rate and large-data-volume scene can be realized.

Description

Technical Field

[0001] The present invention relates to the field of computer technology, in particular to point cloud stitching technology, and specifically to a multi-view point cloud stitching technology using a local feature extraction method. Background Art

[0002] Point cloud stitching is the process of aligning multiple sets of point cloud data obtained from different perspectives or positions into the same coordinate system to form a complete three-dimensional model. Existing point cloud stitching methods are difficult to achieve ideal stitching effects when dealing with point cloud data with repetitive structures or large initial errors, and the influence of noise and outliers further increases the difficulty of stitching. Although traditional ICP algorithms and their variants have solved the problem of point cloud stitching to a certain extent, they mainly rely on the geometric relationship between point clouds and ignore the local features of point clouds. This makes them vulnerable to interference when dealing with point cloud data with complex structures and noise, resulting in inaccurate stitching results.

[0003] A registration method based on local point pair features and the iterative closest point algorithm proposed in the prior art literature realizes efficient registration by improving feature description, key point detection, feature matching and other steps. However, this technology has problems such as difficulty in extracting critical features and low matching accuracy when dealing with objects with complex textures and geometric features. Aiming at the difficult problem of human body point cloud stitching with low overlap rate, the prior art literature proposes a normal constraint point-to-plane ICP stitching method based on the composite features of point clouds, which can effectively select feature points and simplify the amount of point cloud data. However, this method has a large amount of calculation and low stitching accuracy.

[0004] Therefore, the applicant proposes a technology for multi-view point cloud stitching using a local feature extraction method to improve the accuracy and robustness of stitching, while reducing the computational complexity and meeting the requirements of real-time applications. Summary of the Invention

[0005] The object of the present invention is to solve the problems that the prior art cannot complete multi-view point cloud stitching in low overlap rate and large data volume scenarios in real time, resulting in increased computational costs, increased manual processing volume and low stitching accuracy, and to propose a method for realizing multi-view point cloud stitching by using local feature extraction.

[0006] To solve the above technical problems, the technical solution adopted by the present invention is as follows:

[0007] A multi-view point cloud stitching method based on local feature extraction, comprising the following steps:

[0008] Step 1: Obtain the 3D point cloud data of the scene to be stitched, ensure that there is an overlapping area between the point cloud data obtained from two adjacent perspectives, record one of the two adjacent perspectives as scene to be stitched 1, and the 3D point cloud data of this perspective is recorded as the source point cloud, and the other perspective is recorded as scene to be stitched 2, and the 3D point cloud data of this perspective is recorded as the target point cloud;

[0009] Step 2: Process the 3D point cloud data from different perspectives obtained in step 1 to reduce the impact of noise points and abnormal points on subsequent feature extraction and point cloud stitching;

[0010] Step 3: In step 1, a marker point is preset in the overlapping area of ​​the three-dimensional point cloud data of two adjacent viewpoints, and there are distinct colors and geometric features around the marker point;

[0011] Step 4: extract key feature points around the marker points preset in step 3 according to the local feature extraction algorithm, and find a set of matching point pairs of two adjacent perspective three-dimensional point cloud data;

[0012] Step 5: Perform iterative calculation based on the set of matching point pairs obtained in step 4 to obtain the best matching point pair, and then use the total least squares algorithm to estimate the rotation matrix and translation matrix between the point clouds to be spliced;

[0013] Step 6: According to the rotation matrix and translation matrix estimated in step 5, the 3D point cloud data of scene 1 to be stitched is converted to scene 2 to achieve preliminary stitching. The converted 3D point cloud data and the 3D point cloud data in scene 2 to be stitched are fused using the iterative closest point algorithm to achieve multi-view point cloud stitching.

[0014] In step 2, the following steps are specifically included:

[0015] Step 2.1: Obtain 360° point cloud data through a 3D laser scanner. In addition to the main scene that needs to be stitched, there are also a large number of ground information points and outliers, which are removed using a pass-through filter;

[0016] Step 2.2: For the 3D point cloud data processed by the straight-through filter, the noise points generated by measurement errors, environmental interference or scanner failure have obvious distribution characteristics different from those of the surrounding points, and median filtering is used to effectively filter and remove them;

[0017] In step 2.2, the following steps are specifically included:

[0018] Step 2.2.1: For the input source point cloud P1, the template is a sphere D, and the neighborhood of the sampling point M(x0,y0,z0) in the source point cloud is D(x0,y0,z0)|(x-x0) 2 +(y-y0) 2 +(z-z0) 2<d 2 where d is the radius of the sphere, and (x, y, z) is the sorting of all points in the neighborhood of P1 in each channel:

[0019]

[0020] Step 2.2.2: Take the data at the median value as the substitute value at the sampling point M(x0, y0, z0). Median filtering can remove the redundant points in the point cloud while preserving the edge and detail information of the point cloud.

[0021] Step 2.3: For the three-dimensional point cloud data processed by median filtering, there is still a large amount of point cloud data. Using octree voxel downsampling for processing can reduce the data volume and improve the calculation speed.

[0022] In Step 2.3, the following steps are specifically adopted:

[0023] Step 2.3.1: Regard the three-dimensional point cloud space as a large cube, which is the 1st-level node at this time;

[0024] Step 2.3.2: Subdivide the 1st-level node into eight 2nd-level nodes to form an octree structure. If the number of points in the 2nd-level node is less than the set threshold or reaches the preset subdivision depth, stop the subdivision;

[0025] Step 2.3.3: Select the centroid P c =(x c , y c , z c ) in each final sub-node to approximately represent all points in the voxel:

[0026]

[0027] Step 4 specifically includes the following steps:

[0028] Step 4.1: Convert the point cloud RGB model to the HSV model and extract the color feature points around the preset marker points;

[0029] Step 4.2: Extract the ISS feature points according to the geometric information around the preset marker points;

[0030] Step 4.3: Use the fast point feature histogram algorithm to calculate the histogram statistical features of the color feature points and geometric feature points respectively, and fuse the two into a joint feature;

[0031] Step 4.4: Calculate the joint features of the source point cloud and the target point cloud respectively, and select the set of matching point pairs between the two point clouds;

[0032] Step 4.4 specifically includes:

[0033] Step 4.4.1: Combine the obtained color key-point histogram and geometric key-point histogram to obtain a multi-dimensional point cloud feature histogram;

[0034] Step 4.4.2: Compare the histogram features of the key points around each preset marker point to find matching point pairs.

[0035] In step 4.1, it specifically includes the following steps:

[0036] Step 4.1.1: Convert the <r, g, b> attributes to <h, s, v> attributes;

[0037] Step 4.1.2: Set any point P around the marker point i (x i , y i , z i ) with a search radius r, and calculate the weights of all the neighborhood points within the radius with P i as the center and that point;

[0038]

[0039] Step 4.1.3: According to the three-dimensional HSV color vector at P i and the corresponding color vectors at its neighborhood points, calculate the color covariance matrix cov(P i c );

[0040] Step 4.1.4: Calculate the eigenvalues and their corresponding eigenvectors of each color covariance matrix, sort the eigenvectors in ascending order, and record the product of the eigenvalues as the generalized variance value at that point;

[0041] Step 4.1.5: Calculate the generalized variance values of all points around the preset marker point, extract some points as feature points, repeat the above steps until all points around the marker point are traversed, and record the points that meet the conditions as the color feature point set P c .

[0042] In step 4.2, it specifically includes the following steps:

[0043] Step 4.2.1: Calculate the weights of the neighborhood points according to equation (4);

[0044] Step 4.2.2: Calculate the eigenvalue i and eigenvector of each point P sort them in ascending order, and construct a local coordinate system based on the eigenvectors;

[0045] Step 4.2.3: Apply non-maximum suppression to the feature points to reduce the interference of the diffused points along the main direction to the detection. Set the thresholds ε1 and ε2, and the key points need to satisfy

[0046]

[0047] Step 4.2.4: Repeat the above steps until all points are traversed. Finally, all the points that meet the conditions are recorded as the ISS feature point set P ISS 。

[0048] In Step 4.3, it specifically includes the following steps:

[0049] Step 4.3.1: Respectively extract the color feature point set and the geometric feature point set in the point cloud to be stitched, and obtain M key point sets;

[0050] Step 4.3.2: Extract the q nearest neighbor points within the radius r of any key point set in the set, and use the principal component analysis method to extract the best fitting plane including this point and the nearest neighbor points;

[0051] Step 4.3.3: Estimate the minimum eigenvector through singular value decomposition to calculate the surface normal n at M points i ;

[0052] Step 4.3.4: Construct the local coordinate system uvw of the center point and its neighborhood points, and calculate the triple based on the angle change between two points;

[0053]

[0054] Step 4.3.5: Divide the value range of the 3 angle eigenvalues into 11 intervals equally, create an 11-dimensional histogram for the eigenvalues of each interval, and count the number of eigenvalues falling into each sub-interval. This multi-dimensional histogram is denoted as SPFH;

[0055] Step 4.3.6: Calculate the triples of the neighborhood points within their neighborhoods respectively, and assign weights based on distance. Finally, combine and assign the triples to the center point, calculate the histogram, and obtain the statistical features.

[0056]

[0057] Step 5.1: Select any N pairs from the set of matching point pairs, calculate the rotation matrix and translation matrix between the point clouds to be stitched, and transform the source point cloud into the target point cloud coordinate system for rough registration;

[0058] Step 5.2: Calculate the rotation error and translation error between the source point cloud and the target point cloud after rough registration;

[0059] Step 5.3: Select another N pairs of feature points from the matching point pairs, calculate the rotation error and translation error between the two point clouds this time. If the error is higher than the previous result, select the matching point pairs again; if the error is lower than the previous result, replace them. Stop the iteration when the rotation error and translation error are lower than the preset value or the maximum number of iterations is reached.

[0060] Step 5.1 specifically includes:

[0061] Step 5.1.1: Establish a mathematical model for point cloud registration and construct a parameter matrix for the desired rotation matrix;

[0062] Specifically, approximate the coordinate transformation process as a linear model \((A + E)θ=(b + e)\), where \(E\) and \(e\) are the errors generated during the measurement process, the \(m*n\) - dimensional data matrix \(A\) composed of the source point cloud coordinates and the \(m*1\) - dimensional data vector \(b\) composed of the corresponding points in the target coordinate system, and the \(n*1\) - dimensional vector \(θ\) of the parameters to be estimated. Let the augmented matrix \(B = [-b,A]\), and the matrix \(A\), the parameter vector \(θ\), and the data vector \(b\) are as shown in formula (8), where \((x a ,y a ,z a ) is the coordinate in the source point cloud coordinate system, and \((x b ,y b ,z b ) is the coordinate in the Earth Centered Earth Fixed coordinate system;

[0063]

[0064] Step 5.1.2: Calculate the SVD of the augmented matrix \(B\) as shown in formula (9), and store the right singular matrix \(V\), where \(\sum p =diag(σ1,σ2,...,σ p ,0,...,0)\) is a diagonal matrix;

[0065]

[0066] Step 5.1.3: Determine the effective rank \(p\) of the augmented matrix \(B\), that is, only \(p\) of the parameter vectors are linearly independent;

[0067] Step 5.1.4: Calculate the \((p + 1)*(p + 1)\) matrix \(S (p) , where \(ν(i,j)\) is the element in the \(i\) - th row and \(j\) - th column of the right singular matrix \(V\);

[0068]

[0069] Step 5.1.5: Find the inverse matrix \(S (p) of \(S\) and find the optimal least - squares approximation solution according to formula (10), where \(S -(p) , and \(S -1(p)(i + 1, 1) is an element on the (i + 1)-th row and 1st column of matrix S -(p) The element on the (i + 1)-th row and 1st column of the matrix;

[0070] θ TLS θ(i) = S -1(p) (i + 1, 1) / S -1(p) (1, 1) where i = 1, 2,..., p(7).

[0071] Step 6 specifically includes:

[0072] Step 6.1: Use the Iterative Closest Point (ICP) algorithm to register the source point cloud to the target point cloud coordinate system;

[0073] Step 6.1.1: Select a point set p ∈ P in the target point cloud P;

[0074] Step 6.1.2: In the source point cloud Q, select the corresponding point set q ∈ Q such that the value of |q - p| is minimized;

[0075] Step 6.1.3: Calculate the rotation matrix R and translation matrix T based on the corresponding point sets to minimize the distance error calculated by formula (12);

[0076]

[0077] Step 6.1.4: Use the rotation matrix R and translation matrix T to correct the pose of p, obtaining a new corresponding point set p' = {Rp' + t, p' ∈ p};

[0078] Step 6.1.5: Calculate the average distance d between p' and the corresponding point set q using formula (8);

[0079] Step 6.1.6: If d is less than the given distance threshold or the number of iterations reaches the preset maximum value, the iteration stops, and the calculated rotation matrix and translation matrix can achieve point cloud position calibration. Otherwise, return to Step 6.1 until the stop condition is met;

[0080] Step 6.2: Fuse the registered source point cloud and target point cloud in the same coordinate system to achieve point cloud stitching.

[0081] Compared with the prior art, the present invention has the following technical effects:

[0082] 1) By means of the feature point extraction method, the present invention focuses on the neighborhood around the preset marked points and extracts matching point pairs therefrom, reducing the workload of manual processing during the stitching of large amounts of point clouds and improving the feasibility and convenience of the entire stitching process;

[0083] 2) By presetting marker points in the overlapping area, the present invention can achieve low-overlap-rate point cloud stitching, which plays a positive role in promoting subsequent three-dimensional point cloud reconstruction and broadens the applicability of point cloud stitching methods in different complex scenarios;

[0084] 3) Different from the traditional feature extraction method that needs to traverse all point cloud data to find feature points, the present invention accurately locates local feature points related to parameter estimation and more representative by limiting the range of feature extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] The present invention will be further described below in conjunction with the drawings and embodiments:

[0086] Figure 1 is a flowchart of the multi-viewpoint cloud stitching method based on local feature extraction provided by the present invention;

[0087] Figure 2 is a schematic diagram of the time required for each algorithm to achieve point cloud stitching with different densities;

[0088] Figure 3 is a schematic diagram of the error of each algorithm to achieve point cloud stitching with different densities;

[0089] Figure 4 is a schematic diagram of the time required for each algorithm to achieve point cloud stitching with different overlap rates;

[0090] Figure 5 is a schematic diagram of the error of each algorithm to achieve point cloud stitching with different overlap rates. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0091] The object of the present invention is to provide a multi-viewpoint cloud stitching method based on local feature extraction, which can apply a three-dimensional laser scanner to obtain point cloud images in real time and quickly and accurately achieve multi-viewpoint cloud stitching, and output a complete scene to be stitched.

[0092] As Figure 1 shown, the present invention provides a multi-viewpoint cloud stitching method based on local feature extraction, and the stitching method includes:

[0093] Step 1: Select the position of the three-dimensional laser scanner according to the size of the scene to be stitched, ensure that the three-dimensional point cloud data obtained from each view can cover the entire scene to be stitched, and there is an overlapping area between the point cloud data obtained from two adjacent views. Taking the position facing the camera of the three-dimensional laser scanner as the positive direction, the left view among two adjacent views is denoted as the scene to be stitched 1, the three-dimensional point cloud data of this view is denoted as the source point cloud, the right view is denoted as the scene to be stitched 2, and the three-dimensional point cloud data of this view is denoted as the target point cloud. Step 1 specifically includes:

[0094] Step 1.1: There should be no obstacles in front of the selected position for setting up the 3D laser scanner. Before acquiring the 3D point cloud data, the attitude of the 3D laser scanner needs to be adjusted to be parallel to the ground, and the attitude angle at this time is recorded using a gyroscope.

[0095] Step 1.2: Considering that the 3D laser scanner will lose the information at its own position when acquiring the 3D point cloud data, adjacent point clouds need to contain each other's position information.

[0096] In practical applications, there is at least a 10% overlap rate between the point cloud data of two adjacent viewpoints, and 5 marker points are preset in the overlapping area of each scene to be stitched.

[0097] Step 2: Process the 3D point cloud data of different viewpoints obtained in Step 1 to reduce the influence of noise points and abnormal points on subsequent feature extraction and point cloud stitching. Step 2 specifically includes:

[0098] Step 2.1: The 3D laser scanner can acquire 360° point cloud data. In addition to the main scene to be stitched, there are a large number of ground information points and outliers, which are removed using a pass-through filter.

[0099] In practical applications, first determine the filtering range of Figure 3 each dimension of the point cloud and set it in the filter, then let the filter traverse each point, decide whether to keep or discard the point according to whether the coordinates are within the filtering range, and finally output the filtered point cloud.

[0100] Step 2.2: For the 3D point cloud data processed by the pass-through filter, the distribution characteristics of the noise points generated due to measurement errors, environmental interference or scanner failures are significantly different from those of the surrounding points. Median filtering is used to effectively filter and remove them.

[0101] Step 2.2.1: For the input source point cloud P1, the template is a sphere D, and the neighborhood of the sampling point M(x0, y0, z0) in the source point cloud is D(x0, y0, z0)|(x - x0) 2 +(y - y0) 2 +(z - z0) 2 <d 2 where d is the radius of the sphere, and (x, y, z) is the sorting of all points in the neighborhood of P1 in each channel:

[0102]

[0103] Step 2.2.2: Take the data at the median as the replacement value at the sampling point M(x0, y0, z0). Median filtering removes the redundant points in the point cloud while retaining the edge and detail information of the point cloud.

[0104] In practical applications, iterative updating can be used to determine the size of the neighborhood. Set the initial neighborhood value d and the cut-off threshold ε; perform median filtering on the noisy point cloud based on this initial value; calculate the root mean square error of the original point cloud and the denoised point cloud, and the calculation formula is shown in Equation (2):

[0105]

[0106] where (x i , y i , z i ) are the points in the original point cloud, (x′ i , y′ i , z′ i ) are the corresponding points in the denoised point cloud, m is the number of points. If the value of RMSE is greater than the threshold ε, continue the iteration. When the RMSE value is less than the threshold ε or reaches the maximum number of iterations, stop the iteration.

[0107] Step 2.3: For the three-dimensional point cloud data after the median filtering process, there is still a large amount of point cloud data. Use octree voxel downsampling for processing, which can reduce the data volume and improve the calculation speed.

[0108] Step 2.3.1: Consider the three-dimensional point cloud space as a large cube, which is the 1st-level node at this time;

[0109] Step 2.3.2: Subdivide the 1st-level node into eight 2nd-level nodes to form an octree structure. If the number of points in the 2nd-level node is less than the set threshold or reaches the preset subdivision depth, stop the subdivision;

[0110] Step 2.3.3: Select the centroid P c = (x c , y c , z c ) in each final sub-node to approximately represent all the points in this voxel.

[0111]

[0112] In practical applications, the threshold is set to 10 and the subdivision depth is set to 5.

[0113] Step 3: Preset marker points in the overlapping area of the three-dimensional point cloud data of two adjacent viewpoints in Step 1. There are distinct color and geometric features around the marker points;

[0114] Step 4: Extract key feature points around the marker points preset in Step 3 according to the local feature extraction algorithm, and find the set of matching point pairs of the three-dimensional point cloud data of two adjacent viewpoints. Step 4 specifically includes:

[0115] Step 4.1: Convert the point cloud RGB model to the HSV model and extract the color feature points around the preset marker points. Step 4.1 specifically includes:

[0116] Step 4.1.1: Convert the <r, g, b> attributes to <h, s, v> attributes;

[0117] Step 4.1.2: Set any point P i (x i , y i , z i ) around the marker point, calculate the search radius r, and calculate the weights of all neighborhood points within the radius centered at P i and this point;

[0118]

[0119] Step 4.1.3: According to the three-dimensional HSV color vector at P i and the corresponding color vectors at its neighborhood points, calculate the color covariance matrix cov(P i c );

[0120] Step 4.1.4: Calculate the eigenvalues and their corresponding eigenvectors of each color covariance matrix, sort the eigenvectors in ascending order, and record the product of the eigenvalues as the generalized variance value at this point;

[0121] Step 4.1.5: Calculate the generalized variance values of all points around the preset marker points, extract some points as feature points, repeat the above steps until all points around the marker points are traversed, and record the points that meet the conditions as the color feature point set P c .

[0122] In practical applications, extract the points with a generalized variance value greater than 0.1 as candidate feature points.

[0123] Step 4.2: Extract ISS feature points according to the geometric information around the preset marker points. Step 4.2 specifically includes:

[0124] Step 4.2.1: Calculate the weights of neighborhood points according to Equation (4);

[0125] Step 4.2.2: Calculate the eigenvalues and eigenvectors of each point P i , sort them in ascending order, and construct a local coordinate system according to the eigenvectors;

[0126] Step 4.2.3: Use non-maximum suppression for feature points to reduce the interference of diffused points along the main direction on detection, set thresholds ε2 and ε3, and the key points need to meet

[0127]

[0128] Step 4.2.4: Repeat the above steps until all points are traversed. Finally, all points that meet the conditions are recorded as the ISS feature point set P ISS 。

[0129] In practical applications, set the thresholds ε2 = 1 and ε3 = 1.

[0130] Step 4.3: Use the fast point feature histogram algorithm to calculate the histogram statistical features of color feature points and geometric feature points respectively, and fuse the two into a joint feature; Step 4.3 specifically includes:

[0131] Step 4.3.1: Extract the color feature point set and geometric feature point set in the point cloud to be stitched respectively to obtain N key point sets;

[0132] Step 4.3.2: Extract q nearest neighbor points within the radius r of any key point set in the set, and use the principal component analysis method to extract the best fitting plane including this point and the nearest neighbor points;

[0133] Step 4.3.3: Estimate the minimum eigenvector through singular value decomposition to calculate the surface normal n at N points i ;

[0134] Step 4.3.4: Construct the local coordinate system uvw of the center point and its neighborhood points, and calculate the triple based on the angular change between two points;

[0135]

[0136] Step 4.3.5: Divide the value ranges of the 3 angular feature values into 11 intervals equally, create an 11-dimensional histogram for each interval feature value, and count the number of feature values falling within each sub-interval. This 33-dimensional histogram is denoted as SPFH;

[0137] Step 4.3.6: Calculate the triples of neighborhood points within their neighborhoods respectively, and assign weights based on distance. Finally, combine and assign the triples to the center point, calculate the histogram, and obtain the statistical features.

[0138]

[0139] Step 4.4: Calculate the joint features of the source point cloud and the target point cloud respectively, and select the set of matching point pairs between the two point clouds. Step 4.4 specifically includes:

[0140] Step 4.4.1: Combine the obtained color key point histogram and geometric key point histogram to obtain a 66-dimensional point cloud feature histogram;

[0141] Step 4.4.2: Compare the histogram features of the key points around each preset marker point to find matching point pairs.

[0142] Step 5: Perform iterative calculations on the set of matching point pairs obtained in Step 4 to obtain the best matching point pair, and then use the total least squares algorithm to estimate the rotation matrix and translation matrix between the point clouds to be stitched; Step 5 specifically includes:

[0143] Step 5.1: Select any 5 pairs from the set of matching point pairs, calculate the rotation matrix and translation matrix between the point clouds to be stitched, and transform the source point cloud into the target point cloud coordinate system for rough registration; Step 5.1 specifically includes:

[0144] Step 5.1.1: Establish a mathematical model for point cloud registration and construct the parameter matrix of the required rotation matrix;

[0145] In practical applications, the coordinate transformation process is approximately regarded as a linear model (A + E)θ = (b + e), where E and e are the errors generated during the measurement process, the m*n-dimensional data matrix A composed of the source point cloud coordinates and the m*1-dimensional data vector b composed of the corresponding points in the target coordinate system, the n*1-dimensional vector θ composed of the parameters to be estimated, let the augmented matrix B = [-b, A], the matrix A, the parameter vector θ and the data vector b are as shown in formula (8), where (x a , y a , z a ) is the coordinate in the source point cloud coordinate system, and (x b , y b , z b ) is the coordinate in the Earth-Centered Earth-Fixed coordinate system;

[0146]

[0147] Step 5.1.2: Calculate the SVD of the augmented matrix B, as shown in formula (9), and store the right singular matrix V, where ∑ p = diag(σ1, σ2,..., σ p , 0,..., 0) is a diagonal matrix;

[0148]

[0149] Step 5.1.3: Determine the effective rank p of the augmented matrix B, that is, only p of the parameter vectors are linearly independent;

[0150] Step 5.1.4: Calculate the (p + 1)*(p + 1) matrix S (p) , where ν(i, j) is the element in the i-th row and j-th column of the right singular matrix V;

[0151]

[0152] Step 5.1.5: Calculate S (p) 's inverse matrix S -(p) , and calculate the optimal least squares approximation solution according to formula (10), where S -1(p) (i + 1, 1) is the element in the (i + 1)-th row and 1st column of the S -(p) matrix.

[0153] θ TLS (i) = S -1(p) (i + 1, 1) / S -1(p) (1, 1) for i = 1, 2,..., p (18)

[0154] Step 5.2: Calculate the rotation error and translation error between the source point cloud and the target point cloud after rough registration;

[0155] In practical applications, calculate the normal vectors of the ground on the left and right sides of the overlapping area after stitching, record the difference between the two as the rotation error, convert the key points in the source point cloud to the target point cloud coordinate system, calculate the Euclidean distance between the converted key points and the corresponding key points of this point in the target point cloud, and record this distance as the translation error.

[0156] Step 5.3: Select another 5 pairs of feature points from the matching point pairs, calculate the rotation error and translation error between the two point clouds this time. If the error is higher than the previous result, select the matching point pairs again. If the error is lower than the previous result, replace them. Stop the iteration when the rotation error and translation error are lower than the preset value or reach the maximum number of iterations.

[0157] Step 6: Convert the 3D point cloud data of the scene 1 to be stitched into the scene 2 to be stitched according to the rotation matrix and translation matrix estimated in Step 5 to achieve preliminary stitching, and use the Iterative Closest Point (ICP) algorithm to perform point cloud fusion on the converted 3D point cloud data and the 3D point cloud data in the scene 2 to be stitched to achieve multi-viewpoint cloud stitching. Step 6 specifically includes:

[0158] Step 6.1: Use the Iterative Closest Point (ICP) algorithm to register the source point cloud to the target point cloud coordinate system;

[0159] Step 6.1.1: Take a point set p ∈ P in the target point cloud P;

[0160] Step 6.1.2: In the source point cloud Q, take out the corresponding point set q ∈ Q such that the value of |q - p| is the smallest;

[0161] Step 6.1.3: Calculate the rotation matrix R and translation matrix T according to the corresponding point set to minimize the distance error calculated by formula (12);

[0162]

[0163] Step 6.1.4: Use the rotation matrix R and the translation matrix T to correct the pose of p, and obtain a new corresponding point set p′ = {Rp′ + t, p′ ∈ p};

[0164] Step 6.1.5: Use formula (8) to calculate the average distance d between p′ and the corresponding point set q;

[0165] Step 6.1.6: If d is less than the given distance threshold or the number of iterations reaches the preset maximum value, the iteration stops, and the calculated rotation matrix and translation matrix can achieve point cloud position calibration. Otherwise, return to Step 6.1 until the stop condition is met.

[0166] Step 6.2: Fuse the registered source point cloud and target point cloud in the same coordinate system to achieve point cloud stitching.

[0167] Example:

[0168] The present invention uses the Windows 10 operating system, the CPU is i5-7500, a Leica P30 three-dimensional laser scanner is selected, the programming language is MATLAB, and the selected point cloud databases are the TUM dataset and the indoor lidar RGBD dataset.

[0169] Select three point cloud stitching methods, namely FPFH-ICP, ISS-FGR, and RANSAC-TrICP, for comparative experiments. The method proposed by the present invention is denoted as LEPE-ICP. Stitching experiments are carried out on the point cloud dataset and the collected multi-view point clouds respectively, and the evaluation is carried out from four aspects: data volume, overlap rate, stitching accuracy, and stitching efficiency. In the experiment, the search radius r for color feature extraction is 0.05, the threshold ε1 is 0.01, the thresholds ε2 and ε3 in ISS feature extraction are 0.5, 0.5. For median filtering, when the increment of the results of two consecutive iterations is less than 1×10 -6 Or when the number of iterations is greater than 1000, it is judged to converge. The minimum number of points set for octree voxel downsampling is 5, and the maximum subdivision degree is 5. The stitching metrics used are RMSE, MSE, and MAE, and the calculation methods are shown in Equation (13).

[0170]

[0171] The RMSE, MAE, and MSE of each method are compared in Tables 2 and 3. From the data in the comparison charts, for FPFH-ICP, a feature histogram is constructed by calculating the geometric relationship between a point and the points within its neighborhood, mainly focusing on the local geometric information of the point cloud. There is a certain inclination in the stitching results of both datasets, and the stitching effect is slightly worse in the semantically rich RGBD dataset; for ISS-FGR, it has rotational invariance and scale invariance, and the rotation matrix is estimated well on both datasets, with the main error coming from the translational error; for RANSAC-TrICP, it cannot guarantee finding the global optimal solution every time. Although it can obtain a good approximate solution, it may miss the optimal registration result in some cases, and the discrimination threshold and the number of iterations need to be determined through experiments; for LEPE-ICP, good stitching effects can be achieved on both datasets, with small stitching errors.

[0172] Table 1 Comparison of Different Stitching Algorithms in TMU Dataset

[0173] Method RMSE / m MSE / m MAE / m FPFH-ICP 0.3848 0.1481 0.2214 ISS-FGR 0.3361 0.1134 0.1659 RANSAC-TrICP 0.4179 0.1746 0.2814 LEPE-ICP 0.1217 0.0148 0.0962

[0174] Table 2 Comparison of Different Stitching Algorithms in Indoor LiDAR RGBD Dataset

[0175] Method RMSE / m MSE / m MAE / m FPFH-ICP 0.5847 0.3419 0.5295 ISS-FGR 0.3807 0.3088 0.1449 RANSAC-TrICP 0.5016 0.4221 0.3154 LEPE-ICP 0.1332 0.1082 0.0177

[0176] Differences in point cloud density will, to a large extent, affect the accuracy and efficiency of point cloud stitching. It is difficult to extract feature points from low-density point clouds, and it is hard to accurately identify representative feature points. While for high-density point clouds, the feature points are too dense, making the calculation difficult. When matching feature points in low-density point clouds, due to the small number of feature points, mis-matching is likely to occur. For high-density point clouds, the calculation amount of feature point matching is large, and the matching process takes a long time. Therefore, the stitching efficiency for different densities is an important indicator for evaluating the effectiveness of stitching algorithms. Figure 2 is the time required for each algorithm to achieve stitching of point clouds with different densities, Figure 3 is the error of each algorithm in achieving stitching of point clouds with different densities.

[0177] Comparison Figure 2 and Figure 3It can be seen from the data in [the figure] that the stitching time increases with the increase in the amount of point cloud data, and the stitching error decreases with the increase in the point cloud density. The FPFH-ICP algorithm has a relatively high computational complexity when dealing with large-scale or high-density point cloud data, and also has higher accuracy compared to the other two comparison methods; the ISS-FGR algorithm has a high time cost range, but the accuracy improvement is not significant; the RANSAC-TrICP algorithm has the smallest time range and a high accuracy improvement, and is suitable for point cloud stitching of large data volumes; while the LEPE-ICP algorithm can complete the pose correction of the point cloud without performing feature comparison on all points in the point cloud, greatly reducing the time cost of large-scale point cloud registration. The stitching efficiency is increased by 49.4%, the stitching error is reduced by 56.4%, and the error range is small and the stability is high.

[0178] When the overlap rate of two point clouds is relatively high, there are more corresponding points available for use in the feature extraction and matching process, providing richer geometric constraint information. However, too high an overlap rate will increase the workload and computational amount of subsequent point cloud reconstruction. The stitching accuracy at different overlap rates is an important indicator for evaluating the effectiveness of the stitching algorithm. Figure 4 is the time required for each algorithm to achieve point cloud stitching at different overlap rates, Figure 5 is the error of each algorithm in achieving point cloud stitching at different overlap rates. It can be seen from the figure that the stitching efficiency decreases with the increase in the point cloud overlap rate, and the stitching accuracy increases with the increase in the overlap rate. Compared with the other three stitching algorithms, in the case of a low overlap rate, the stitching time of LEPE-ICP is at least 12.471 s, an increase of 73.3%, and the stitching error is at least 0.2279 s, a decrease of 65.9%.

[0179] In summary, the present invention effectively improves the point cloud stitching accuracy when the overlap rate between two point clouds is too low in a large data volume scenario, facilitating subsequent three-dimensional point cloud reconstruction. By using the local feature point extraction algorithm, key points can be extracted around the manually set marker points, effectively describing the color information and geometric information here, and obtaining matching point pairs according to the fast point feature histogram algorithm. Compared with the existing feature matching algorithms, it does not need to traverse the global point cloud for feature extraction, relies on a small number of reference points to estimate the transformation matrix of the entire point cloud, and does not require complex iterative operations, improving the efficiency of point cloud stitching. Compared with the method of manually setting reference points commonly used in engineering, which has problems such as low accuracy of matching feature points and large workload, this method can automatically match feature points, and uses total least squares estimation in parameter estimation to take into account noise and measurement errors of three-dimensional laser scanners, enhancing the stitching effect.

Claims

1. A multi-view point cloud stitching method based on local feature extraction, characterized in that: The following steps are involved: Step 1: Obtain the 3D point cloud data of the scene to be stitched, ensure that there is an overlapping area between the point cloud data obtained from two adjacent perspectives, record one of the two adjacent perspectives as scene to be stitched 1, and the 3D point cloud data of this perspective is recorded as the source point cloud, and the other perspective is recorded as scene to be stitched 2, and the 3D point cloud data of this perspective is recorded as the target point cloud; Step 2: Process the 3D point cloud data from different perspectives obtained in step 1 to reduce the impact of noise points and abnormal points on subsequent feature extraction and point cloud stitching; Step 3: In step 1, a marker point is preset in the overlapping area of ​​the three-dimensional point cloud data of two adjacent viewpoints, and there are distinct colors and geometric features around the marker point; Step 4: extract key feature points around the marker points preset in step 3 according to the local feature extraction algorithm, and find a set of matching point pairs of two adjacent perspective three-dimensional point cloud data; Step 5: Perform iterative calculation based on the set of matching point pairs obtained in step 4 to obtain the best matching point pair, and then estimate the rotation matrix and translation matrix between the point clouds to be spliced; Step 6: According to the rotation matrix and translation matrix estimated in step 5, the 3D point cloud data of the scene to be stitched 1 is converted to the scene to be stitched 2 to achieve preliminary stitching, and the converted 3D point cloud data is fused with the 3D point cloud data in the scene to be stitched 2 to achieve multi-view point cloud stitching.

2. The method according to claim 1, characterized in that In step 2, the following steps are specifically included: Step 2.1: Obtain 360° point cloud data through a 3D laser scanner. In addition to the main scene that needs to be stitched, there are also a large number of ground information points and outliers, which are removed using a pass-through filter; Step 2.2: For the 3D point cloud data processed by the straight-through filter, the noise points generated by measurement errors, environmental interference or scanner failure have obvious distribution characteristics different from those of the surrounding points, and median filtering is used to effectively filter and remove them; In step 2.2, the following steps are specifically included: Step 2.2.1: For the input source point cloud P1, the template is a sphere D, and the neighborhood of the sampling point M(x0,y0,z0) in the source point cloud is D(x0,y0,z0)|(x-x0) 2 +(y-y0) 2 +(z-z0) 2 <d 2 Where d is the radius of the sphere and (x, y, z) is the order of all points in the neighborhood of P1 in each channel: Step 2.2.2: Take the data at the median as the replacement value at the sampling point M (x0, y0, z0). The median filter removes redundant points in the point cloud while retaining the edge and detail information of the point cloud; Step 2.3: For the three-dimensional point cloud data processed by median filtering, there is still a large amount of point cloud data. Using octree voxel downsampling for processing can reduce the amount of data and improve the calculation speed.

3. The method according to claim 2, characterized in that In step 2.3, the following steps are specifically adopted: Step 2.3.1: Consider the 3D point cloud space as a large cube, which is a level 1 node at this time; Step 2.3.2: Subdivide the level 1 node into eight level 2 nodes to form an octree structure. If the number of points in the level 2 node is less than the set threshold or reaches the preset subdivision depth, stop subdividing; Step 2.3.3: Select the centroid P in each final child node c =(x c ,y c ,z c ) to approximate all points within the voxel:

4. The method according to any one of claims 1 to 3, characterized in that: Step 4 specifically includes the following steps: Step 4.1: Convert the point cloud RGB model to the HSV model and extract the color feature points around the preset marker points; Step 4.2: Extract ISS feature points based on the geometric information around the preset marker points; Step 4.3: Use the fast point feature histogram algorithm to calculate the histogram statistical features of color feature points and geometric feature points respectively, and merge the two into a joint feature; Step 4.4: Calculate the joint features of the source point cloud and the target point cloud respectively, and select the set of matching point pairs between the two point clouds; Step 4.4 specifically includes: Step 4.4.1: Combine the obtained color key point histogram and geometric key point histogram to obtain a multi-dimensional point cloud feature histogram; Step 4.4.2: Compare the histogram features of the key points around each preset marker point to find matching point pairs.

5. The method according to claim 4, characterized in that In step 4.1, the following steps are specifically included: Step 4.1.1:<r,g,b> Attributes are converted to<h,s,v> property; Step 4.1.2: Set any point P around the marker point i (x i ,y i ,z i ) with a search radius r, and calculate all the i is the weight of the neighborhood point of the circle center and the point Step 4.1.3: According to P i The three-dimensional HSV color vector at the point and the corresponding color vector at its neighboring point are used to calculate the color covariance matrix cov(P i c ); Step 4.1.4: Calculate the eigenvalues ​​of each color covariance matrix and its corresponding eigenvector, sort the eigenvectors in ascending order, and record the result of multiplying the eigenvalues ​​as the generalized variance value at that point; Step 4.1.5: Calculate the generalized variance value of all points around the preset marking point, extract some points as feature points, repeat the above steps until all points around the marking point are traversed, and record the points that meet the conditions as the color feature point set P c .

6. The method according to claim 5, characterized in that In step 4.2, the following steps are specifically included: Step 4.2.1: Calculate the weight of the neighborhood points according to formula (4); Step 4.2.2: Calculate each point P i The eigenvalue of and the eigenvector Sort them in ascending order, and construct a local coordinate system based on the eigenvectors; Step 4.2.3: Use non-maximum suppression on feature points to reduce the interference of diffusion points along the main direction on detection, set thresholds ε1 and ε2, and the key points must meet Step 4.2.4: Repeat the above steps until all points are traversed, and finally all points that meet the conditions are recorded as the ISS feature point set P ISS .

7. The method according to claim 5, characterized in that In step 4.3, the following steps are specifically included: Step 4.3.1: Extract the color feature point set and the geometric feature point set in the point cloud to be spliced ​​respectively to obtain M key point sets; Step 4.3.2: Extract the q nearest neighbor points within the radius r of any key point in the set, and use the principal component analysis method to extract the best fit table including the point and the nearest neighbor points; Step 4.3.3: Calculate the surface normal n at M points by estimating the minimum eigenvector through singular value decomposition i ; Step 4.3.4: Construct the local coordinate system uvw of the center point and its neighborhood points, and calculate the triplet based on the angle change between the two points; Step 4.3.5: Divide the value range of the three angle eigenvalues ​​into 11 intervals, create an 11-dimensional histogram for each interval eigenvalue, and count the number of eigenvalues ​​falling in each sub-interval. This multidimensional histogram is recorded as SPFH; Step 4.3.6: Calculate the triplets of the neighborhood points in their neighborhoods respectively, and assign weights based on distance. Finally, assign the triplet combination to the center point, calculate the histogram, and obtain the statistical features.

8. The method according to claim 4, characterized in that Step 5.1: Select N pairs of matching points from the set, calculate the rotation matrix and translation matrix between the point clouds to be stitched, and transform the source point cloud into the target point cloud coordinate system for rough registration; Step 5.2: Calculate the rotation error and translation error between the source point cloud and the target point cloud after rough registration; Step 5.3: Select another N pairs of feature points from the matching point pairs, calculate the rotation error and translation error between the two point clouds, if the error is higher than the previous result, select the matching point pair again, if the error is lower than the previous result, replace it, when the rotation error and translation error are lower than the preset value or the maximum number of iterations is reached, stop the iteration.

9. The method according to claim 8, characterized in that Step 5.1 specifically includes: Step 5.1.1: Establish a mathematical model for point cloud registration and construct the parameter matrix of the required rotation matrix; Specifically, the coordinate transformation process is approximated as a linear model (A+E)θ=(b+e), where E and e are the errors generated during the measurement process, the m*n dimensional data matrix A consisting of the source point cloud coordinates and the m*1 dimensional data vector b consisting of the corresponding points in the target coordinate system, and the n*1 dimensional vector θ consisting of the parameters to be estimated. Let the augmented matrix B=[-b,A], the matrix A, the parameter vector θ and the data vector b are shown in formula (8), where (x a ,y a ,z a ) is the coordinate in the source point cloud coordinate system, (x b ,y b ,z b ) are the coordinates in the Earth-centered Earth-fixed coordinate system; Aθ=b b=[x b ,y b ,z b ,1] T <h2 style=";text-align:left;direction:ltr">θ = [a1,a2,a3,t1,a4,a5,a6,t2,a7,a8,a9,t3,0,0,0,1]<h2 style=";text-align:left;direction:ltr"> T Step 5.1.2: Calculate the SVD of the augmented matrix B, as shown in formula (9), and store the right singular matrix V, where ∑ p =diag(σ1,σ2,...,σ p ,0,...,0) is a diagonal matrix; Step 5.1.3: Determine the effective rank p of the augmented matrix B, that is, only p of the parameter vectors are linearly independent; Step 5.1.4: Calculate the (p+1)*(p+1) matrix S (p) , where ν(i,j) is the element in the i-th row and j-th column of the right singular matrix V; Step 5.1.5: Find S (p) The inverse matrix S -(p) , and the optimal least squares approximate solution is obtained according to formula (10), where S -1(p) (i+1,1) is S -(p) The element in the i+1th row and 1st column of the matrix; θ TLS (i)=S -1(p) (i+1,1) / S -1(p) (1,1)i=1,2,...,p(7)。 10. The method according to any one of claims 5 to 9, characterized in that Step 6 specifically includes: Step 6.1: Use the iterative closest point algorithm to register the source point cloud to the target point cloud coordinate system; Step 6.1.1: Select a point set p∈P in the target point cloud P; Step 6.1.2: In the source point cloud Q, take out the corresponding point set q∈Q so that the value of |qp| is minimized; Step 6.1.3: Calculate the rotation matrix R and the translation matrix T based on the corresponding point set so that the distance error calculated by formula (12) is minimized; Step 6.1.4: Use the rotation matrix R and translation matrix T to correct the position of p and obtain a new corresponding point set p′={Rp′+t,p′∈p}; Step 6.1.5: Use formula (8) to calculate the average distance d between p′ and the corresponding point set q; Step 6.1.6: If d is less than the given distance threshold or the number of iterations reaches the preset maximum value, the iteration stops, and the calculated rotation matrix and translation matrix can realize the point cloud position calibration, otherwise return to step 6.1 until the stop condition is met; Step 6.2: Fuse the registered source point cloud and target point cloud in the same coordinate system to achieve point cloud stitching.

Images