A point cloud registration method based on local curvature and its neighbor features

Abstract

The present invention discloses a point cloud registration method based on local curvature and its neighboring features. The present invention introduces the local curvature of the point cloud and its related feature information into the point cloud feature extraction process, thereby improving the accuracy and speed of point cloud registration. Firstly, the point cloud downsampling algorithm is completed based on the point density distribution to reduce the point cloud density, and the radius filtering algorithm is used to remove the point cloud noise points; secondly, by extracting the curvature value of the detection point, the distance between the detection point and the pseudo-plane formed by its neighboring points, and the standard deviation of the distance between the detection point and the neighboring point fitting plane, the three feature parameters are combined to complete the parameter construction of the point cloud feature, and the fast point feature histogram FPFH is used as the feature descriptor; finally, the SAC‑IA sampling consistency initial registration algorithm and the ICP iterative closest point algorithm are combined to solve the optimal rotation matrix R and translation matrix t to complete the entire point cloud registration process.

Description

Technical Field

[0001] The present invention belongs to the technical field of computer vision, and specifically designs a point cloud registration method based on local curvature and its neighbor features. Background Art

[0002] Point cloud is increasingly used in 3D vision technology. Point cloud data collection tools are developing rapidly, data collection costs are reduced, and efficiency is improved. This has led to the application of point cloud technology in machine vision, cultural relic restoration, autonomous driving, drones, and other fields. The main task of point cloud registration is to rotate and translate point clouds in different coordinate systems, and finally integrate them into a complete point cloud in the same coordinate system. According to the properties of the target point cloud and the source point cloud, it is generally divided into non-rigid registration and rigid registration. Non-rigid registration refers to the existence of non-rigid transformations such as deformation and stretching between the target point cloud and the source point cloud. In actual application scenarios, some registration problems require non-rigid registration, such as medical diagnosis and face recognition. There is still a lot of room for research on non-rigid registration; rigid registration refers to the absence of deformation between the source point cloud and the target point cloud, and only rigid transformations such as rotation and translation occur. The actual application range is wider, but the existing point cloud registration technology has the disadvantages of slow registration speed and easy to fall into local optimal solutions. How to complete point cloud registration quickly and with high quality is still an issue of great research value. Summary of the invention

[0003] In order to overcome the problems existing in the above-mentioned prior art, the purpose of the present invention is to provide a point cloud registration method based on local curvature and its neighboring features. The present invention uses three parameters, namely the curvature value of the detection point of the three-dimensional point cloud, the distance from the detection point to the pseudo-plane formed by its neighboring points, and the standard deviation fitting of the distance between the detection point and the fitting plane of the neighboring points, to jointly construct feature parameters. The feature parameters have the advantages of representativeness, universality, and fast extraction speed, which can greatly reduce the time consumption of point cloud registration.

[0004] In order to achieve the above object, the present invention adopts the following technical solution:

[0005] A point cloud registration method based on local curvature and its neighbor features includes the following steps:

[0006] Step 1: preprocessing the point cloud data to ensure that the original features of the point cloud are not destroyed, reduce the point cloud density, and remove noise points in the point cloud;

[0007] Step 2: constructing point cloud feature parameters, combining the average curvature of the detection point, the Euclidean distance from the detection point to the neighboring point to the fitted plane, and the standard deviation of the distance between the detection point and the neighboring point to the fitted plane, to jointly construct feature parameters and extract feature points;

[0008] Step 3: Point cloud registration. The FPFH fast point feature histogram algorithm is used to calculate the feature descriptor and combined with the SAC-IA sampling consistency initial registration algorithm to complete the initial registration of the point cloud. The ICP iterative closest point algorithm is used for fine registration of the point cloud.

[0009] The point cloud in the point cloud registration is divided into a target point cloud and a source point cloud. The target point cloud and the source point cloud use complete paired point clouds or multi-view point clouds collected by a laser radar, which is suitable for small point cloud data or large industrial point cloud data.

[0010] The step 1 is specifically as follows:

[0011] Step 1.1: Perform 3D voxel grid processing on the point cloud, scatter the point cloud data in the voxel grid, and use the point closest to the center of gravity of each voxel grid to replace the remaining points in the voxel grid. This method reduces the number of point clouds while ensuring that the microscopic morphology of the point cloud remains unchanged;

[0012] Step 1.2: Set different search radii r according to the density distribution of neighboring point clouds, calculate the number of points inside the sphere with radius r, and if it is less than the threshold, it will be regarded as noise points and removed.

[0013] The step 2 is specifically as follows:

[0014] Step 2.1: Calculate the average curvature H of the detection point, as described in formula (1):

[0015]

[0016] in:

[0017] k1——maximum principal curvature;

[0018] k2——minimum principal curvature;

[0019] Step 2.2: Calculate the fitting plane formed by the k nearest neighbor points of the detection point according to the least square method. The plane is as shown in equation (2):

[0020] Ax+By+Cz+D=0 (2)

[0021] in:

[0022] A, B, C, D——are known constants, and A, B, C are not zero at the same time;

[0023] Calculate the Euclidean distance d from the detection point to the fitting plane, as shown in formula (3):

[0024]

[0025] in:

[0026] d——Euclidean distance from the detection point to the fitting plane of the neighboring point;

[0027] x1, y1, z1——the coordinates of the detection point;

[0028] Step 2.3: Calculate the standard deviation σ of the projection distance between the detection point and its k nearest neighbor points to the fitting plane, as described in formula (4):

[0029]

[0030] in:

[0031] r i ——The distance between the detection point and each neighboring point to the fitting plane;

[0032] μ——the average distance between the detection point and each neighboring point to the fitting plane;

[0033] k——the sum of the number of points of the detection point and its neighboring points;

[0034] Step 2.4: Construct the characteristic parameter β based on the three parameters in steps 2.1, 2.2 and 2.3, as described in formula (5):

[0035] β=σ·H·d,Δd>0 (5)

[0036] in:

[0037] σ——Standard deviation of the projection distance between the detection point and its k nearest neighbor points to the fitting plane;

[0038] H——mean curvature;

[0039] Δd——the difference between d and μ;

[0040] Step 2.5: According to the characteristic parameter β proposed in step 2.4, when Δd is greater than 0, the detection point is the outermost point among the neighboring points in the area and is characteristic; when Δd is less than 0, the detection point and its neighboring points are intertwined and cannot be defined as a characteristic point, so the points with Δd less than 0 are first eliminated during calculation; in addition, the points with β less than the average value of β are eliminated, and the remaining points are used as characteristic points.

[0041] The step 3 is specifically as follows:

[0042] Step 3.1: Calculate the fast point feature histogram (FPFH) based on the spatial difference between the extracted feature point and its neighboring feature points to accurately describe the spatial geometric properties within the point neighborhood. First, define a three-dimensional coordinate system uvw and calculate the deviation between each point in the k-neighborhood of point P and the estimated normal. The coordinate system uvw is described as formula (6):

[0043]

[0044] in:

[0045] n s ——P s The surface normal of a point;

[0046] P t - detection point;

[0047] P s ——Any neighboring point of the detection point;

[0048] In the uvw coordinate system, calculate n s 、n t The angle deviations α, β, θ between them are called P s The three characteristic elements of the point are put into the histogram to form P s The simplified point feature histogram (SPFH) of a point is described as formula (7):

[0049]

[0050] in:

[0051] q——P s With P t Euclidean distance between two points;

[0052] n t ——P t The surface normal of a point;

[0053] α——n t The angle with the v-axis;

[0054] β——P t To P s The vector of the point and n s The angle of

[0055] θ——n t The angle between the projection vector on the u and w planes and the u axis;

[0056] Calculate the detection point P t The characteristic triple elements of each neighboring point in the k-neighborhood are calculated, and the corresponding simplified point feature histogram SPFH is calculated. The weighted statistics of each simplified point feature histogram SPFH in the neighborhood are performed to obtain the fast point feature histogram FPFH, which is described as formula (8):

[0057]

[0058] in:

[0059] FPFH(P t )——P t Fast point feature histogram of points;

[0060] SPFH(Pt )——P t Simplified point feature histogram of points;

[0061] m——P t The number of neighboring points of a point;

[0062] P m ——P t The neighboring points of a point;

[0063] w m ——P t With neighboring point P m distance;

[0064] SPFH(P m )——P m Simplified point feature histogram of points;

[0065] Step 3.2: The SAC-IA sampling consistency initial registration algorithm uses the FPFH feature descriptor as input to reduce the number of iterations. It searches for similar features between the source point cloud and the target point cloud as corresponding points through the KD tree, traverses and calculates the minimum error of the registration result of each SAC-IA algorithm. When the transformation quality calculation obtains the minimum error, the rotation and translation matrix at this time is the optimal solution for the initial registration. The transformation quality D is described as formula (9):

[0066]

[0067] in:

[0068] D(e i )——transformation quality of the i-th group;

[0069] e i ——represents the error value between the i-th group of corresponding points;

[0070] t e - a predetermined value;

[0071] Step 3.3: Use the ICP iterative closest point algorithm to find the closest point between the two sets of point clouds, and iteratively solve the optimal rigid transformation matrix;

[0072] Step 3.4: According to the RMSE root mean square error between the corresponding point pairs after fine registration, the final completion degree of the point cloud registration is measured. The larger the RMSE root mean square error, the larger the error, and vice versa, the smaller the error, as described in formula (10);

[0073]

[0074] in:

[0075] RMSE(P,Q)——The overall root mean square error of the two groups of point clouds P and Q;

[0076] n——The total number of corresponding point pairs in the two point clouds P and Q;

[0077] X i ——is the Euclidean distance between each set of corresponding points;

[0078] ——is the true Euclidean distance between each set of corresponding points.

[0079] Compared with the prior art, the present invention has the following advantages:

[0080] (1) More reliable

[0081] In terms of point cloud preprocessing, the present invention sets different search radii according to the point cloud density distribution to perform point cloud downsampling operations. In this way, the defect of easily generating holes after downsampling due to uneven density distribution is avoided.

[0082] (2) More accurate

[0083] The feature point combined with the FPFH feature descriptor for SAC-IA initial alignment algorithm in the present invention can reduce the root mean square error by about 70% compared with directly using the FPFH feature descriptor combined with the SAC-IA algorithm; at the same time, it has higher accuracy than ISS, Harris and other feature points combined with the FPFH feature descriptor for SAC-IA initial alignment.

[0084] (3) Faster

[0085] The feature point extraction algorithm of the present invention is superior to the ISS algorithm, Harris algorithm and SIFT algorithm in extraction speed, with the speed increased by about 84%, 55% and 93% respectively. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1a , Figure 1b and Figure 1c The following are the effects of adding Gaussian noise to the face point cloud, using a fixed radius filtering algorithm to reduce the noise on the noisy face point cloud, and performing adaptive radius filtering on the noisy face point cloud based on the density distribution.

[0087] Figure 2 Characterization calculation flow chart.

[0088] Figure 3a , Figure 3b Figure 3c and Figure 3d These are rabbit feature extraction effect diagrams using the feature extraction algorithm of the present invention, the ISS feature extraction algorithm, the SIFT feature extraction algorithm, and the Harris feature extraction algorithm, respectively.

[0089] Figure 4Time consumption graph of feature extraction algorithm.

[0090] Figure 5a , Figure 5b Figure 5c and Figure 5d The Bunny point cloud data registration effect diagrams use the original point cloud data, the feature point set extracted by the present invention, the feature point set extracted by the ISS algorithm, and the feature point set extracted by the Harris algorithm as the registration input sets respectively. DETAILED DESCRIPTION

[0091] The present invention is described in detail below with reference to the accompanying drawings and examples.

[0092] The face, cat, cow point cloud data, Stanford University Bunny, Dragon point cloud data are used as input for research, and the ISS, SIFT, Harris feature extraction algorithms are used to compare with the positive feature extraction algorithm proposed in this invention. The FPFH feature descriptor is used for feature description, the SAC-IA algorithm is used for initial registration experiments, and the classic ICP algorithm is used for fine registration experiments. The entire point cloud registration process is implemented, including the following steps:

[0093] Step 1: Preprocess the point cloud of the experimental data.

[0094] Since the density distribution of face data is relatively uneven, it is easy to have a hole area when using a fixed radius value for radius filtering. In this invention, different radii r are used for radius filtering processing for different density areas in combination with the point cloud density distribution. Figure 1a , Figure 1b and Figure 1c The preprocessing effects are shown using noisy face point cloud, fixed radius filtering and density distribution radius filtering respectively.

[0095] Step 2: Feature extraction and feature description.

[0096] First, in the feature extraction stage, the curvature of the detection point, the Euclidean distance from the detection point to the neighboring point to the fitted plane, and the standard deviation of the distance between the detection point and the neighboring point to the fitted plane are extracted. The three parameters are used to jointly construct the feature parameters and then the feature points are automatically selected to obtain the feature point sets of the target point cloud and the source point cloud, respectively. t , P s .

[0097] Secondly, the fast point feature histogram FPFH is used to calculate the feature histogram of the spatial difference between the extracted feature points and their neighboring feature points, accurately describing the spatial geometric properties within the point neighborhood. Figure 2 A flowchart showing steps 1 and 2.

[0098] Figure 3a , Figure 3b , Figure 3c and Figure 3d The effect diagrams of using the feature extraction algorithm of the present invention, the ISS algorithm, the SIFT algorithm and the Harris algorithm on the Bunny point cloud data are respectively shown. It can be seen from the diagrams that: the feature points extracted by the present invention are concentrated in the ears, neck, feet, tail and other parts of the Bunny point cloud data with significant features; the feature points extracted by the ISS algorithm are scattered in the entire point cloud space; the feature points extracted by the SIFT algorithm are mostly concentrated in the edge area of ​​the point cloud; the feature points extracted by the Harris algorithm are scattered in the entire point cloud space and have no obvious features.

[0099] Figure 4 The algorithm time consumption of feature extraction using the algorithm of the present invention, ISS algorithm, Harris algorithm and SIFT algorithm for four groups of point cloud data of Bunny, Dragon, Cat and Cow respectively. It can be seen from the figure that the computational time consumption of the feature extraction algorithm of the present invention is better than that of the other three algorithms.

[0100] Step 3: initial point cloud registration and fine point cloud registration.

[0101] After the initial registration using the SAC-IA sampling consistency initial registration algorithm, the source point cloud and the target point cloud are basically overlapped, but there is still a certain offset between points. In order to achieve accurate registration and minimize the distance error between points, the ICP iterative closest point algorithm is used. Suppose the feature sets of the source point cloud and the target point cloud after initial registration are P and i * With Q i * , the general process of ICP iterative closest point algorithm is as follows:

[0102] (1) By constructing a KD tree, calculate P i With Q i For corresponding points of two point sets, using KD tree can greatly improve the search efficiency;

[0103] (2) Calculate the rotation and translation matrix values ​​based on the corresponding point pairs obtained in the first step;

[0104] (3) According to the rotation and translation matrix obtained in the second step, the source point cloud feature set P i * Perform rigid body transformation to obtain the rotated point set P i * ;

[0105] (4) Calculate P i * With Q i * The average distance of corresponding point pairs is less than a given threshold, and the iteration is stopped. Otherwise, P is used. i*’ Continue iterating with new input and stop when the average distance meets the threshold.

[0106] Figure 5 shows the registration results of Bunny when using different point sets as registration input objects.

[0107] Through analysis, it can be seen that the feature point extraction algorithm designed in the present invention reduces the amount of calculation of huge point clouds, significantly improves the speed of point cloud registration, and the feature extraction speed is better than ISS, Harris, SIFT and other methods. At the same time, it has excellent registration accuracy, and has good effects on both small-scale point cloud data and large-scale industrial point cloud data. It is efficient, universal and robust.

Claims

1. A point cloud registration method based on local curvature and its neighbor features, characterized by: The steps include: Step 1: preprocessing the point cloud data to ensure that the original features of the point cloud are not destroyed, reduce the point cloud density, and remove noise points in the point cloud; Step 2: constructing point cloud feature parameters, combining the average curvature of the detection point, the Euclidean distance from the detection point to the neighboring point to the fitted plane, and the standard deviation of the distance between the detection point and the neighboring point to the fitted plane, to jointly construct feature parameters and extract feature points; Step 3: Point cloud registration. The FPFH fast point feature histogram algorithm is used to calculate the feature descriptor and combined with the SAC-IA sampling consistency initial registration algorithm to complete the initial registration of the point cloud. The ICP iterative closest point algorithm is used for fine registration of the point cloud. The point cloud in the point cloud registration is divided into a target point cloud and a source point cloud. The target point cloud and the source point cloud use complete paired point clouds or multi-view point clouds collected by a laser radar, which is suitable for small point cloud data or large industrial point cloud data. The step 3 is specifically as follows: Step 3.1: Calculate the fast point feature histogram FPFH based on the spatial difference between the extracted feature point and its neighboring feature points to accurately describe the spatial geometric properties within the point neighborhood. First, define a three-dimensional coordinate system uvw and calculate the deviation between each point in the k-neighborhood of point P and the estimated normal. In the uvw coordinate system, calculate P s The surface normal of the point n s , P t The surface normal of the point n t The angle deviations α, β, θ between them are called P s The three characteristic elements of the point are put into the histogram to form P s Simplified point feature histogram SPFH of the point; calculate the detection point P t The characteristic triple elements of each neighboring point in the k-neighborhood are calculated, and the corresponding simplified point feature histogram SPFH is calculated. The weighted statistics of each simplified point feature histogram SPFH in the neighborhood are performed to obtain the fast point feature histogram FPFH. Step 3.2: The SAC-IA sampling consistency initial registration algorithm uses the FPFH feature descriptor as input to reduce the number of iterations. It searches for similar features between the source point cloud and the target point cloud as corresponding points through the KD tree, traverses and calculates the minimum error of the registration result of each SAC-IA algorithm. When the transformation quality calculation obtains the minimum error, the rotation and translation matrix at this time is the optimal solution for the initial registration; Step 3.3: Use the ICP iterative closest point algorithm to find the closest point between the two sets of point clouds, and iteratively solve the optimal rigid transformation matrix; Step 3.4: Measure the final completion degree of point cloud registration based on the RMSE root mean square error between corresponding point pairs after fine registration. The larger the RMSE root mean square error, the larger the error, and vice versa.

2. The point cloud registration method based on local curvature and its neighboring features according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1: Perform 3D voxel grid processing on the point cloud, scatter the point cloud data in the voxel grid, and use the point closest to the center of gravity of each voxel grid to replace the remaining points in the voxel grid. This method reduces the number of point clouds while ensuring that the microscopic morphology of the point cloud remains unchanged; Step 1.2: Set different search radii r according to the density distribution of neighboring point clouds, calculate the number of points inside the sphere with radius r, and if it is less than the threshold, it will be regarded as noise points and removed.

3. The point cloud registration method based on local curvature and its neighboring features according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2.1: Calculate the average curvature H of the detection point, as described in formula (1): in: k1——maximum principal curvature; k2——minimum principal curvature; Step 2.2: Calculate the fitting plane formed by the k nearest neighbor points of the detection point according to the least square method. The plane is as shown in equation (2): Ax+By+Cz+D=0 (2) in: A, B, C, D——are known constants, and A, B, C are not zero at the same time; Calculate the Euclidean distance d from the detection point to the fitting plane, as shown in formula (3): in: d——Euclidean distance from the detection point to the neighboring point fitting plane; x1, y1, z1——coordinates of the detection point; Step 2.3: Calculate the standard deviation σ of the projection distance between the detection point and its k nearest neighbor points to the fitting plane, as described in formula (4): in: r i ——The distance between the detection point and each neighboring point to the fitting plane; μ——the average distance between the detection point and each neighboring point to the fitting plane; k——the sum of the number of points of the detection point and its neighboring points; Step 2.4: Construct the characteristic parameter β based on the three parameters in steps 2.1, 2.2 and 2.3, as described in formula (5): β=σ·H·d,Δd>0 (5) in: σ——Standard deviation of the projection distance between the detection point and its k nearest neighbor points to the fitting plane; H——mean curvature; Δd——the difference between d and μ; Step 2.5: According to the characteristic parameter β proposed in step 2.4, when Δd is greater than 0, the detection point is the outermost point among the neighboring points in the area and is characteristic; when Δd is less than 0, the detection point and its neighboring points are intertwined and cannot be defined as a characteristic point, so the points with Δd less than 0 are first eliminated during calculation; in addition, the points with β less than the average value of β are eliminated, and the remaining points are used as characteristic points.

4. According to the point cloud registration method based on local curvature and its neighboring features according to claim 1, the step 3 is specifically: Step 3.1: Calculate the fast point feature histogram FPFH based on the spatial difference between the extracted feature point and its neighboring feature points to accurately describe the spatial geometric properties within the point neighborhood. First, define a three-dimensional coordinate system uvw and calculate the deviation between each point in the k-neighborhood of point P and the estimated normal. The coordinate system uvw is described as formula (6): in: n s ——P s The surface normal of a point; P t - detection point; P s ——Any neighboring point of the detection point; In the uvw coordinate system, calculate n s 、n t The angle deviations α, β, Call it P s The three characteristic elements of the point are put into the histogram to form P s The simplified point feature histogram SPFH of a point is described as formula (7): in: α——n t The angle with the v-axis; q——P s With P t Euclidean distance between two points; β——P t To P s The vector of points and n s The angle of n t ——P t The surface normal of a point; ——n t The angle between the projection vector on the u and w planes and the u axis; Calculate the detection point P t The characteristic triple elements of each neighboring point in the k-neighborhood are calculated, and the corresponding simplified point feature histogram SPFH is calculated. The weighted statistics of each simplified point feature histogram SPFH in the neighborhood are performed to obtain the fast point feature histogram FPFH, which is described as formula (8): in: FPFH(P t )——P t Fast point feature histogram of points; SPFH(P t )——P t Simplified point feature histogram of points; m——P t The number of neighboring points of a point; P m ——P t The neighboring points of a point; w m ——P t With neighboring point P m distance; SPFH(P m )——P m Simplified point feature histogram of points; Step 3.2: The SAC-IA sampling consistency initial registration algorithm uses the FPFH feature descriptor as input to reduce the number of iterations. It searches for similar features between the source point cloud and the target point cloud as corresponding points through the KD tree, traverses and calculates the minimum error of the registration result of each SAC-IA algorithm. When the transformation quality calculation obtains the minimum error, the rotation and translation matrix at this time is the optimal solution for the initial registration. The transformation quality D is described as formula (9): in: D(e i )——transformation quality of the i-th group; e i ——represents the error value between the i-th group of corresponding points; t e - a predetermined value; Step 3.3: Use the ICP iterative closest point algorithm to find the closest point between the two sets of point clouds, and iteratively solve the optimal rigid transformation matrix; Step 3.4: According to the RMSE root mean square error between the corresponding point pairs after fine registration, the final completion degree of the point cloud registration is measured. The larger the RMSE root mean square error, the larger the error, and vice versa, the smaller the error, as described in formula (10); in: RMSE(P,Q)——The overall root mean square error of the two groups of point clouds P and Q; n——The total number of corresponding point pairs in the two point clouds P and Q; X i ——is the Euclidean distance between each set of corresponding points; ——is the true Euclidean distance between each set of corresponding points.

Images