Point cloud registration method based on directional description and multi-feature constraint

Abstract

The invention discloses a point cloud registration method based on directional description and multi-feature constraint, and aims to solve the problem of insufficient registration precision of point cloud data due to the problems of noise, sparsity, data missing and the like. In the prior art, point cloud registration mainly depends on an end-to-end transformation matrix estimation method or a local feature-based method, and the local feature-based method still has an optimization space in key point feature description and feature matching stages. In order to solve the problem of rotation invariance in key point feature description, a directional description structure network DDN is provided on the basis of a USIP key point extraction network, and high-robustness description of key point features is achieved. Besides, in order to solve the problems of error correspondence and low efficiency in feature matching, a multi-feature constraint matching module MFC is designed, and registration precision and robustness are further improved by introducing a local coherence matching structure and an optimized closest point iterative algorithm ICP. Experimental results show that the method is superior to a traditional algorithm in registration time, rotation matrix and translation vector precision.

Description

Technical Field

[0001] The present invention relates to the technical field of computer vision and three-dimensional point cloud processing, and in particular to a point cloud registration method based on directional description and multi-feature constraints. Background Art

[0002] Point cloud registration aims to align multi-frame point cloud data to a unified coordinate system through rotation and translation, thereby restoring the geometric structure of the scene. However, sensor data is often affected by factors such as noise, sparsity, data loss, and perspective changes, making traditional registration methods challenging in terms of accuracy and robustness. In the feature-based point cloud registration process, it is usually necessary to first extract key points, build local blocks and local coordinate systems around the key points to extract feature descriptors, then match key points based on the descriptors, and finally calculate and optimize the transformation matrix between the two sets of point clouds. However, in practical applications, the point density of the key point neighborhood in point cloud data collected from different perspectives may be different, which will lead to the instability of the local reference frame, thereby increasing the registration error. Summary of the invention

[0003] In view of this, the purpose of the present invention is to provide a point cloud registration method based on directional description and multi-feature constraints to solve the problems of rotation invariance and mismatching in point cloud registration, which is widely applicable to scenarios such as SLAM and 3D reconstruction.

[0004] To achieve the above object, the present invention adopts the following technical solution: a point cloud registration method based on directional description and multi-feature constraints, comprising the following steps:

[0005] Step 1: Extract point cloud key points through USIP deep learning network;

[0006] Step 2: Use the directional description structure designed by the present invention to build a DDNet network and perform training;

[0007] Step 3: Optimize the matching process of key points of source point cloud and target point cloud using local coherence matching method;

[0008] Step 4: Use ICP optimization based on normal vector constraints to further eliminate false matches;

[0009] Step 5: Iteratively solve the transformation matrix from the source point cloud to the target point cloud;

[0010] Step 6: Experimentally prove the superiority of the point cloud registration algorithm proposed in this invention.

[0011] In a preferred embodiment, the step 1 specifically includes:

[0012] Step 11: Point Cloud Perform random transformation T to obtain point cloud Both are input into the FPN (Feature Proposal Network) layer;

[0013] Step 12: The FPN layer first extracts the input point cloud In the example, M marker points are selected according to the farthest point sampling FPS algorithm, denoted as Then each point in the input point cloud X is grouped according to the distance from Node into M groups, and the result is Where K m It represents the number of points in the original point cloud divided into m groups, and then the points divided into each Node are decentralized, that is, get Then all the points are passed to the first fully connected network in groups;

[0014] Step 13: After passing through the first fully connected network, the features of each point are obtained, and then the features of M nodes can be obtained through the maximum pooling operation Then perform KNN on M Nodes and get Then the coordinate information is decentralized, that is get Then, M nodes correspond to a feature of (C+3, M, K), which is then passed to the second fully connected network;

[0015] Step 14: After passing through the second fully connected network, we get (C 2 ,M) features and compare them with the features before KNN (C 1 ,M) performs Concat operation on the feature dimension to obtain (C 1 +C 2 ,M), then the feature vector is input into MLP to obtain the output of (4,M), the first three dimensions represent the relative position information of the key points, and the fourth dimension is transformed by the Softmax function to characterize the uncertainty level of the key point significance, where a lower uncertainty value indicates a more reliable key feature location;

[0016] Step 15: Input Point Cloud X and After the FPN layer, the final output is Q 1 and σ 1 and Q 2 and σ 2 , representing the points from point cloud X and The coordinates of the key points extracted and the uncertainty corresponding to each key point;

[0017] USIP designs two loss functions, loss function L pBy quantifying the matching degree of key points before and after rotation; the loss function L p as follows:

[0018]

[0019] Loss function L c Based on probability, the correct key point position is identified, and an uncertainty value is provided for each identified key point; the loss function L c as follows:

[0020]

[0021] In a preferred embodiment, the step 2 specifically includes: imagining the key point and its neighboring points as being contained in a regular octahedron, and then observing these points from the normal directions of the eight faces of the regular octahedron;

[0022] Step 21: First select the bottom edge of the face as the x-axis of the new coordinate system, and then use the midpoint of the bottom edge to the vertex as the y-axis of the new coordinate system; then determine the z-axis according to the right-hand rule. First, set the coordinates of the octahedron vertices as follows: The coordinates of vertex a are The coordinates of vertex b are (i, i, 0), the coordinates of vertex c are (i, -i, 0), the coordinates of vertex d are (-i, -i, 0), the coordinates of vertex e are (-i, i, 0), and the coordinates of vertex f are Compute a new coordinate system perpendicular to face ABC where Then the z-axis vector is:

[0023]

[0024] After calculating the vectors of the x, y, and z axes, each vector is normalized, so the new coordinate system of the ABC surface is:

[0025] Step 22: Rotate the key point and its surrounding neighboring points according to the new coordinate system of the octant to obtain eight new point clouds, and then send the new point clouds to the new reasoning network, namely the directional description network.

[0026] In a preferred embodiment, in step 2, the directional description network is trained using a triplet loss function; during the training process, the triplet loss function makes the network more robust by grouping anchor samples Anchor, similar key points Positive and non-similar key points Negative and training them simultaneously; the mathematical formula of the triplet loss function is as follows:

[0027]

[0028] d i(a,p)=||f(a)-f(p)||,d i (a,n)=||f(a)-f(n)||

[0029] Where a represents the anchor point, the input sample anchor, p represents the positive sample, n represents the negative sample, d i (a,p) represents the distance difference between f(a) and f(p) after a and p are described by network features. Similarly, d i (a,n) represents the distance difference between f(a) and f(n) after a and n are described by network features, N represents the number of sample groups, τ is a hyperparameter, which is set to 1 as the cutoff value, and L is the loss function;

[0030] During the training process, anchors, positive samples and negative samples are grouped into triplets; the triplets are simultaneously fed into the directional description network DDNet; for each group of samples, DDNet can extract three (1,32)-dimensional feature vectors.

[0031] In a preferred embodiment, the step 3 specifically includes the following steps:

[0032] Step 31: Traverse each key point in the source point cloud and find the four other key points closest to the key point; then, select two from these four neighboring key points to form three points p with the key point 0 ,p 1 ,p 2 , find the corresponding key point q in the target point cloud according to the feature description 0 ,q 1 ,q 2 , i.e. a triple;

[0033] Step 32: Determine the three key points p in the source point cloud 0 ,p 1 ,p 2 The triangle formed corresponds to the key point q of the target point cloud 0 ,q 1 ,q 2 Whether the shapes of the triangles are roughly the same, the formula for judging whether the triangle shapes meet the conditions is as follows:

[0034]

[0035] Step 33: The qualified triples will calculate the minimum circumscribed circle according to the three key points of the source point cloud and the target point cloud, and construct two spheres with the center of each circumscribed circle as the center of the sphere and the radius of the circumscribed circle as the radius of the sphere. The calculation method for constructing the sphere is as follows:

[0036] First, constrain the three points to be at the same distance from the center of the sphere:

[0037]

[0038] Then solve jointly and eliminate the variable R:

[0039]

[0040] Simplify the formula to obtain the linear equation system about the coordinates of the sphere center:

[0041]

[0042] Solve the linear equations to obtain the coordinates of the center of the sphere and the radius:

[0043]

[0044] The two spheres obtained will cover the points in a certain area of ​​the source point cloud and the target point cloud respectively, and the points in each sphere will replace the source point cloud and the target point cloud and be sent to ICP for solution.

[0045] In a preferred embodiment, the step 4 is specifically as follows: finding the rotation matrix R and the translation vector t, and the mathematical model formula of the point cloud registration is as follows:

[0046]

[0047] Among them, q i is a point in the target point cloud, p i for q i The corresponding point in the source point cloud;

[0048] Define the centroids of the source and target point clouds:

[0049]

[0050] Simplify the optimization terms of the mathematical model:

[0051]

[0052] Optimize the project i -qR(p i -p) is zero after summing, then the function simplifies to:

[0053]

[0054] In a preferred embodiment, in step 4, the centroid coordinates of each point in the source point cloud and the target point cloud are first calculated:

[0055] q i * =q i -q,p i *=p i -p

[0056] Then optimize the first error term to calculate R:

[0057]

[0058] Expand the error term of the rotation matrix:

[0059]

[0060] Where R T R = I and the first error is independent of the rotation matrix, then the actual optimization function is:

[0061]

[0062] Define the matrix Matrix:

[0063]

[0064] Perform SVD decomposition on the matrix Matrix:

[0065] Matrix=UΣV

[0066] When the original matrix Matrix is ​​full rank, it means that the original matrix has N non-zero singular values, so the diagonal matrix Σ will have N non-zero elements; in this case, U and V represent the standard orthogonal basis of the original data space and the target data space, respectively, and Σ represents the mapping relationship between the two spaces; then the required R is:

[0067] R=UV T

[0068] Finally, substitute the obtained R into the second term of the error and calculate t:

[0069] t=q-Rp

[0070] In a preferred embodiment, in step 4, three criteria are used to eliminate false matches in the ICP process:

[0071] (1) Point distance exceeds threshold: When the distance between a point in the source point cloud and a point in the target point cloud exceeds a certain set threshold, these point pairs will be eliminated;

[0072] (2) Curvature approximation: If the curvature difference between two points is too large, these point pairs will be eliminated;

[0073] (3) Normal vector direction approximation: If the normal vector directions of two matching points are significantly different, such point pairs will also be considered as incorrect matches and removed;

[0074] According to the above formulas such as Matrix = UΣV, three eigenvalues ​​and corresponding eigenvectors are decomposed, where the eigenvalue is λ 0 ,λ 1 ,λ 2 ,λ 0 <λ 1 <λ 2 , the normal vector is λ 0 The corresponding eigenvector; the calculation formula of curvature δ is as follows:

[0075]

[0076] The ICP mathematical model based on normal vector constraints is derived as follows. First, the source point cloud is defined as The target point cloud is Then expand the midpoints of each point cloud, that is, the normal vector:

[0077]

[0078] After the source point cloud is transformed by the matrix T(R, t):

[0079]

[0080] The source point cloud is transformed into coordinates, and then the error between the point in the source point cloud and the corresponding point is defined as e ij , the calculation formula is as follows:

[0081]

[0082] The final mathematical model is defined as:

[0083]

[0084] At this point, the ICP based on the normal vector constraint becomes a nonlinear least squares problem like the traditional method, and the LM method is used to solve it as follows:

[0085] (H+λI)ΔT=b

[0086] H=∑J i T J i

[0087]

[0088] In a preferred embodiment, in step 5, first, the point cloud data of the source point cloud is set as: The key point is After the directional descriptor, the feature corresponding to the key point is described as Set the point cloud data corresponding to the target point cloud as: The key point is The feature description corresponding to the key point is Set the number of inliers inline to 0 and the inlier rate ε threshold to 0.75 to determine whether the source point cloud meets the output conditions after the transformation matrix, set the optimal inlier rate ε' to 0, the optimal transformation matrix R', and the optimal translation vector t';

[0089] 1) Traverse the key points of the source point cloud, recorded as The four key points closest to this key point are The corresponding feature description is

[0090] 2) The four adjacent key points and their corresponding feature descriptions are input into the local coherence matching structure for judgment. If the shape requirements are met, the point cloud of the smallest circumscribed sphere area of ​​the source point cloud divided by the three key points is output one by one as the new source point cloud. And output the point cloud of the minimum circumscribed sphere area divided by the corresponding three key points as the new target point cloud

[0091] If none of them meet the requirements, return to the first step and continue to traverse the key points;

[0092] 3) Input P' and Q' into the ICP mathematical model based on normal vector constraints and construct a matrix, then calculate the rotation matrix R and translation vector t from the key points of the source point cloud to the key points of the target point cloud; then transform the points in the input source point cloud P and rotate them to the coordinate system of the target point cloud Q;

[0093] 4) For all transformed points in P, find the closest point in the target point cloud and calculate the Euclidean distance d between them. a , if d a <0.05, then inline+=1, and after the judgment, calculate the inline rate ε=inline / N 1 Then, the internal point rate ε>0.4 is judged. If it is true, the optimal solution R and t are output. Otherwise, ε'<ε is judged. If it is true, ε', R', t' are updated. If it is not true, the key point is judged. Is the sequence number i greater than the number of key points n of the source point cloud? If so, output R', t', otherwise continue to traverse the key points of the source point cloud.

[0094] In a preferred embodiment, in step 6, the test set in the ModelNet40 standard data set is tested; certain preprocessing is performed on the experimental test set to simulate the overlap rate and sensor noise in actual working conditions, as follows:

[0095] Step 61: Conduct experiments on data with an overlap rate of 50%; the source point cloud refers to 25% of the original data in the test set randomly cropped, and the target point cloud is another 25% randomly cropped from the original data of the test set, and ensure that this part of the data does not overlap with the cropped part of the source point cloud to achieve an overlap rate of about 50%, and then randomly rotate the target point cloud x, y, and z axes in the range of [-90°, 90°], and randomly translate in the range of [-1m, 1m]; at the same time, save the rotation matrix R and translation vector t of each set of source point cloud and the corresponding target point cloud as the real transformation parameters; this group of experimental data sets is collectively referred to as data set C below;

[0096] Step 62: Experiments are conducted on data with an overlap rate of about 50% and containing noise. In this set of experimental data sets, both the source point cloud and the target point cloud are based on the point cloud data in data set C, and Gaussian noise with a mean of 0 and a variance of 0.01 is added to fit the actual sensor noise; this set of experimental data sets is collectively referred to as data set D below;

[0097] The indicators for experimental evaluation of point cloud registration effects include: mean square error RMSE R and RMSE t , root mean square error RMSE R and RMSE t , mean absolute error MAE R and MAE t And the average time Time required to complete the registration; Assume that the calculated value and the true value of the Euler angle are θ pre ,θ gt , the calculated value and true value of the translation vector are t pre ,t gt , then the specific formula for calculating the evaluation index is as follows:

[0098]

[0099]

[0100] Compared with the prior art, the present invention has the following beneficial effects: In view of the problem of poor anti-interference ability of the traditional key point local reference frame, the present invention proposes a directional description structure to ensure that the description of the key point is rotationally invariant. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] Figure 1 Schematic diagram of the USIP algorithm detection network structure of the preferred embodiment of the present invention, wherein (a) is the USIP training process, and (b) is the FPN (Feature Proposal Network) network structure;

[0102] Figure 2 A regular octahedron structure diagram of a preferred embodiment of the present invention;

[0103] Figure 3 It is a schematic diagram of the directional description structure of a preferred embodiment of the present invention;

[0104] Figure 4 A schematic diagram of a training network for describing the characteristics of a preferred embodiment of the present invention;

[0105] Figure 5 A schematic diagram of a triplet of a preferred embodiment of the present invention;

[0106] Figure 6 Schematic diagram of a spherical region of key points of a preferred embodiment of the present invention, wherein (a) is a source point cloud and key points, and (b) is a target point cloud and key points;

[0107] Figure 7 It is a schematic diagram of a local coherence matching structure of a preferred embodiment of the present invention;

[0108] Figure 8 A schematic diagram of a matching structure based on multi-feature constraints according to a preferred embodiment of the present invention;

[0109] Fig. 9 This is an example diagram of a comparative test of a preferred embodiment of the present invention. DETAILED DESCRIPTION

[0110] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0111] It should be noted that the following detailed descriptions are illustrative and are intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present application belongs.

[0112] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application; as used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or their combinations.

[0113] A point cloud registration method based on directional description and multi-feature constraints, reference Figure 1-9 , including the following steps:

[0114] Step 1: Extract point cloud key points through USIP deep learning network;

[0115] The key point extraction method used in this invention is a learning-based algorithm, using the detect module in the USIP network. The original author's open source model is used for key point extraction tasks, i.e., transfer learning, and no additional training is performed on the key point extraction module. The network structure of the detection module can be found in the attached Figure 1 , the specific key point detection network principle is as follows

[0116] (1) First, the point cloud Perform random transformation T to obtain point cloud Both are input to the FPN layer together.

[0117] (2) The FPN layer first obtains the input point cloud In the example, M nodes are selected according to the farthest point sampling (FPS) algorithm, denoted as Then each point in the input point cloud X is grouped according to the distance from Node into M groups, and the result is Where K m It represents the number of points in the original point cloud divided into m groups, and then the points divided into each Node are decentralized, that is, get Then all the points are passed to the first fully connected network in groups.

[0118] (3) After passing through the first fully connected network, the features of each point can be obtained, and then the features of M nodes can be obtained through the maximum pooling operation The paper then performs KNN (nearest neighbor classification) on M nodes and obtains Then the coordinate information is decentralized, that is get Then, M nodes correspond to a feature of (C+3, M, K), which is then passed to the second fully connected network.

[0119] (4) After passing through the second fully connected network, we get (C 2 ,M) features and compare them with the features before KNN (C 1 ,M) performs Concat operation on the feature dimension to obtain (C 1 +C 2 ,M), and then the feature vector is input into MLP to obtain the output of (4,M). The first three dimensions represent the relative position information of the key points, and the fourth dimension is transformed by the Softmax function to characterize the uncertainty level of the key point significance, where a lower uncertainty value indicates a more reliable key feature location.

[0120] (5) Input point cloud X and After the FPN layer, the corresponding Q 1 and σ 1 and Q 2 and σ 2 . Represents the point cloud X and The key points extracted from the image and the uncertainty corresponding to each key point.

[0121] To ensure that the network can perform its tasks effectively, USIP designed two loss functions to supervise the training process of the key point extraction network. One of the loss functions focuses on the point-to-point correspondence, with the aim of ensuring that the extracted key points are highly consistent after the corresponding inverse coordinate transformation before and after the point cloud data undergoes a rotation transformation. This means that even if the input point cloud is rotated, the network can identify the same or similar key points, thereby maintaining the consistency of the key points. This loss function achieves this by quantifying the degree of matching of the key points before and after the rotation, so that the key points learned by the network are rotation invariant. The loss function is as follows:

[0122]

[0123] Another loss function focuses on the uncertainty of keypoint extraction. This loss function is probability-based and aims to optimize the uncertainty score given by the network when extracting each keypoint. With this design, the network not only needs to identify the correct keypoint location, but also provides an uncertainty value for each identified keypoint. This approach encourages the network to output high-quality keypoint predictions and ensures that these predictions have reliable confidence scores, thereby improving the overall keypoint detection performance. The loss function is as follows:

[0124]

[0125] In summary, the USIP key point detection module shows significant advantages over other traditional algorithms. The module is specially designed to identify key points that remain unchanged under rotation transformation. In addition, it effectively prevents excessive aggregation between key points and avoids the problem of overly dense distribution of key points. Through such a strategy, the USIP key point detection module can also maintain the stability of the number of key point extraction, and can output a consistent number of key points even under different input conditions, which provides a reliable basis for subsequent processing and analysis.

[0126] Step 2: Use the directional description structure designed by the present invention to build a DDNet network and perform training;

[0127] The design principle of the directional description structure is to ensure that the feature description of the key points and their neighboring points is invariant to rotation transformation. This means that even if these points are rotated in space, their feature descriptions, which usually include information such as direction, size, and shape, should remain unchanged or change very little. To achieve this goal, this paper proposes a method: imagine the key points and their neighboring points as being contained in a regular octahedron, and then observe these points from the normal direction of the eight faces of the regular octahedron. See the attached diagram for a schematic diagram. Figure 2 .

[0128] A regular octahedron is a polyhedron with eight equilateral triangular faces, and each vertex is a common vertex of four faces. For each face of a regular octahedron, you can first select the base of the face as the x-axis of the new coordinate system, and then use the midpoint of the base to the vertex as the y-axis of the new coordinate system. Then determine the z-axis according to the right-hand rule, first set the coordinates of the octahedron vertices as follows:

[0129]

[0130] Compute a new coordinate system perpendicular to face ABC where Then the z-axis vector is:

[0131]

[0132] After calculating the vectors of the x, y, and z axes, each vector is normalized, so the new coordinate system of the ABC surface is: Similarly, according to the above calculation method, the new coordinate system of each face of the regular octahedron is as follows:

[0133]

[0134]

[0135] After completing this step, the key point and its surrounding neighboring points are rotated according to the new coordinate system of the octant to obtain eight new point clouds. Next, these point clouds are fed into an inference network, which analyzes the eight point clouds and outputs eight 32-dimensional feature vectors. In order to obtain a compact feature description, these feature vectors generated by the network are processed through a maximum pooling process to provide a comprehensive feature representation for each key point. The Directed Description Network (DDNet) proposed in this invention is shown in the attached figure. Figure 3 .

[0136] In order to train the network, the present invention adopts a triplet loss function. During the training process, the triplet loss function makes the network more robust by grouping anchor samples (Anchor), similar key points (Positive) and non-similar key points (Negative) and training them simultaneously. Its core goal is to ensure that the network generates similar feature vectors for similar instances, and generates feature vectors with significant differences for different instances, thereby improving the ability to distinguish features. See the attached training network for details. Figure 4 The mathematical formula of the triplet loss function is as follows:

[0137]

[0138] d i (a,p)=||f(a)-f(p)||,d i (a,n)=||f(a)-f(n)||

[0139] Where a represents anchor (i.e., anchor point, input sample), p represents positive sample (sample of the same type as anchor), n represents negative sample (sample of different key points), d i (a,p) represents the distance (Euclidean distance) difference between f(a) and f(p) after a and p are described by network features. Similarly, d i (a,n) represents the distance difference between f(a) and f(n) after a and n are described by network features, N represents the number of sample groups, τ is a hyperparameter, which is set to 1 as the dividing value, and L is the loss function.

[0140] During the training process, anchors, positive samples (samples similar to anchors), and negative samples (samples completely different from anchors) are grouped into triplets. These combinations are fed into the Directed Description Network (DDNet) at the same time. For each set of samples, DDNet can extract three (1,32)-dimensional feature vectors. The goal of training is to optimize the network weights through the triplet loss function, minimize the feature distance of samples of the same type, and maximize the feature distance of samples of different types. In this way, the feature description learned by DDNet has a compact representation for key points of similar shapes, while there are significant differences for key points of different shapes. This not only enhances the discriminability of the feature description, but also improves its generalization ability, enabling the network to better handle unseen input data and achieve higher accuracy and robustness in key point detection and matching tasks.

[0141] Step 3: Optimize the matching process of key points of source point cloud and target point cloud using local coherence matching method;

[0142] In actual work, the overlap between the source point cloud and the target point cloud is usually less than 50% and contains noise. At this time, matching only by describing key points will inevitably establish wrong point pair matching relationships, thus affecting the accuracy of point cloud registration.

[0143] Considering that the correctly corresponding key points are locally correlated and the incorrectly corresponding key points tend to be randomly distributed, the present invention proposes a matching structure based on local coherence. The specific idea is:

[0144] (1) Traverse each key point in the source point cloud and find the four other key points closest to the key point. Then, select two of the four neighboring key points to form three points p with the key point. 0 ,p 1 ,p 2 , find the corresponding key point q in the target point cloud according to the feature description 0 ,q 1 ,q 2 , that is, a triple. See the attached diagram Figure 5 .

[0145] (2) Determine the three key points p in the source point cloud 0 ,p 1 ,p 2 The triangle formed corresponds to the key point q of the target point cloud 0 ,q 1 ,q 2 Whether the shapes of the triangles are roughly the same, the formula for judging whether the triangle shapes meet the conditions is as follows:

[0146]

[0147] (3) The qualified triples will calculate the minimum circumscribed circle based on the three key points of the source point cloud and the target point cloud, and construct two spheres with the center of each circumscribed circle as the center of the sphere and the radius of the circumscribed circle as the radius of the sphere. The calculation method for constructing the sphere is as follows:

[0148] First, constrain the three points to be at the same distance from the center of the sphere:

[0149]

[0150] Then solve jointly and eliminate the variable R:

[0151]

[0152] Simplify the formula to obtain the linear equation system about the coordinates of the sphere center:

[0153]

[0154] Solve the linear equations to obtain the coordinates of the center of the sphere and the radius:

[0155]

[0156] The two spheres obtained will cover the points in a certain area of ​​the source point cloud and the target point cloud respectively, and the points in each sphere will replace the source point cloud and the target point cloud and be sent to ICP for solution. Figure 6 .

[0157] Assuming that the number of key points in the source point cloud is N, this method not only reduces N(N-1)(N-2) / 6 possibilities to 6N, but also greatly reduces the number of feature transformation solutions by judging the similarity of each triangle. This structure can find the best triplet faster, thereby improving the accuracy of point cloud registration and reducing the number of iterations. See the attached figure for the overall structure flow chart of local coherence matching. Figure 7 .

[0158] Step 4: Use ICP optimization based on normal vector constraints to further eliminate false matches;

[0159] The ultimate goal of point cloud registration is to find the rotation matrix R and translation vector t. The mathematical model formula of point cloud registration is as follows:

[0160]

[0161] Among them, q i is a point in the target point cloud, p i Q i At the corresponding point in the source point cloud, the meaning of this formula can be expressed as solving a set of rotation matrices R and translation vectors t so that the error between the key points in the source point cloud and the corresponding points in the target point cloud after transformation is minimized.

[0162] Define the centroids of the source and target point clouds:

[0163]

[0164] Simplify the optimization items of the mathematical model:

[0165]

[0166] Optimize the project i -qR(p i -p) is zero after summing, then the function simplifies to:

[0167]

[0168] The above optimization function has two error terms. The first error variable is only the rotation matrix R, and the second error variable is the rotation matrix R and the translation vector t. However, the second term is only related to the centroid of each point cloud. Therefore, after obtaining R in the solution process, the second error term is set to zero to get t. The detailed calculation process is as follows:

[0169] First, calculate the centroid coordinates of each point in the source point cloud and the target point cloud:

[0170] q i * =q i -q,p i * =p i -p

[0171] Then optimize the first error term and calculate R:

[0172]

[0173] Expand the error term of the rotation matrix:

[0174]

[0175] Where R T R = I and the first error is independent of the rotation matrix, then the actual optimization function is:

[0176]

[0177] Define the matrix Matrix:

[0178]

[0179] Perform SVD decomposition on the matrix Matrix:

[0180] Matrix=UΣV

[0181] When the original matrix Matrix is ​​full rank, it means that the original matrix has N non-zero singular values, so the diagonal matrix Σ will have N non-zero elements. In this case, U and V represent the standard orthogonal basis of the original data space and the target data space, respectively, and Σ represents the mapping relationship between the two spaces. Then the required R is:

[0182] R=UV T

[0183] Finally, substitute the obtained R into the second term of the error and calculate t:

[0184] t=q-Rp

[0185] The above is the calculation principle of the iterative solution of the iterative closest point algorithm (ICP). In each iteration, ICP finds the nearest matching point in the target point cloud for each point in the source point cloud, and then uses the least squares method to solve the rigid body transformation matrix and update the position of the source point cloud to gradually reduce the spatial distance error between the two sets of point clouds.

[0186] In order to further eliminate false matches and improve the registration accuracy, the present invention introduces an optimization method based on normal vector constraints, taking into full consideration the problem that the points in the source point cloud and the target point cloud are difficult to correspond one to one in real situations. Therefore, the present invention adopts three criteria to eliminate false matches in the ICP process:

[0187] (1) Point distance exceeds threshold: Considering the distribution characteristics of points in the point cloud, when the distance between the points in the source point cloud and the points in the target point cloud exceeds a certain set threshold, these point pairs will be eliminated because they are likely not to belong to the same object surface or are not under the same rigid body transformation.

[0188] (2) Curvature approximation: Curvature approximation is also an important criterion. Since points on the same object surface usually have similar curvature properties, if the curvature difference between two points is too large, the match between them may be wrong, so such point pairs will also be excluded.

[0189] (3) Normal vector direction approximation: The normal vector of a point on the surface can well reflect the local geometric structure of the point. If the normal vector directions of two matching points are significantly different, they are likely to come from different surfaces or in different orientations. Such point pairs will also be considered as mismatches and removed.

[0190] According to the above formulas such as Matrix = UΣV, three eigenvalues ​​and corresponding eigenvectors are decomposed, where the eigenvalue is λ 0 ,λ 1 ,λ 2 ,λ 0 <λ 1 <λ 2 , the normal vector is λ 0 The corresponding eigenvector. The calculation formula of curvature δ is as follows:

[0191]

[0192] The ICP mathematical model based on normal vector constraints is derived as follows. First, the source point cloud is defined as The target point cloud is Then expand the midpoints of each point cloud, that is, the normal vector:

[0193]

[0194] After the source point cloud is transformed by the matrix T(R, t):

[0195]

[0196] The source point cloud is transformed into coordinates, and then the error between the point in the source point cloud and the corresponding point is defined as e ij , the calculation formula is as follows:

[0197]

[0198] The final mathematical model is defined as:

[0199]

[0200] At this time, the ICP based on the normal vector constraint becomes a nonlinear least squares problem as in the traditional method. The present invention adopts the LM method to solve it, as follows:

[0201] (H+λI)ΔT=b

[0202] H=∑J i T J i

[0203]

[0204] Step 5: Iteratively solve the transformation matrix from the source point cloud to the target point cloud;

[0205] This step will introduce in detail the feature matching structure based on local coherence matching in step 3 and ICP based on normal vector constraints in step 4. The specific process is as follows:

[0206] First, set the point cloud data of the source point cloud as: The key point is After the directional descriptor, the feature corresponding to the key point is described as Set the point cloud data corresponding to the target point cloud as: The key point is The feature description corresponding to the key point is The number of inliers inline is set to 0, and the inlier rate ε threshold is set to 0.75 to determine whether the source point cloud meets the output conditions after the transformation matrix. The optimal inlier rate ε' is set to 0, the optimal transformation matrix R', and the optimal translation vector t'.

[0207] (1) Traverse the key points of the source point cloud, recorded as The four key points closest to this key point are The corresponding feature description is

[0208] (2) The four adjacent key points and their corresponding feature descriptions are input into the local coherence matching structure for judgment. If the shape requirements are met, the point cloud of the smallest circumscribed sphere area of ​​the source point cloud divided by the three key points is output one by one as the new source point cloud. And output the point cloud of the minimum circumscribed sphere area divided by the corresponding three key points as the new target point cloud If none of them meet the requirements, return to the first step and continue traversing the key points.

[0209] (3) Input P' and Q' into the ICP mathematical model based on normal vector constraints and construct the matrix, then calculate the rotation matrix R and translation vector t from the key points of the source point cloud to the key points of the target point cloud. Then transform the points in the input source point cloud P and rotate them to the coordinate system of the target point cloud Q.

[0210] (4) For all transformed points in P, find the closest point in the target point cloud and calculate the Euclidean distance d between them. a , if d a <0.05, then inline+=1, and after the judgment, calculate the inline rate ε=inline / N 1 Then, the internal point rate ε>0.4 is judged. If it is true, the optimal solution R and t are output. Otherwise, ε'<ε is judged. If it is true, ε', R', t' are updated. If it is not true, the key point is judged. Whether the sequence number i is greater than the number of key points n of the source point cloud, if so, output R', t', otherwise continue to traverse the key points of the source point cloud.

[0211] The flowchart corresponding to the above process is attached Figure 8 .

[0212] Step 6: Prove the superiority of the point cloud registration algorithm proposed by the present invention through experiments;

[0213] The experiment uses the test set in the ModelNet40 standard data set for testing. Some preprocessing is performed on the experimental test set to simulate the overlap rate and sensor noise in actual working conditions, as follows:

[0214] (1) Experiments are conducted on data with an overlap rate of 50%. Here, the source point cloud refers to 25% of the original data in the test set that is randomly cropped, and the target point cloud refers to another 25% of the original data in the test set that is randomly cropped, and it is ensured that this part of the data does not overlap with the cropped part of the source point cloud to achieve an overlap rate of approximately 50%. The target point cloud x, y, and z axes are randomly rotated in the range of [-90°, 90°], and randomly translated in the range of [-1m, 1m]. At the same time, the rotation matrix R and translation vector t of each set of source point clouds and the corresponding target point clouds are saved as the actual transformation parameters. This group of experimental data sets is collectively referred to as data set C below.

[0215] (2) Experiments were conducted on data with an overlap rate of about 50% and containing noise. In this set of experimental data sets, both the source point cloud and the target point cloud are based on the point cloud data in data set C, and Gaussian noise with a mean of 0 and a variance of 0.01 is added to fit the actual sensor noise. This set of experimental data sets is collectively referred to as data set D below.

[0216] The indicators used to evaluate the effect of point cloud registration include: mean square error (MSE R and MSE t ), root mean square error (RMSE R and RMSE t ), mean absolute error (MAE R and MAE t ) and the average time required to complete the registration (Time). Assume that the calculated and true values ​​of the Euler angles are θ pre ,θ gt , the calculated value and true value of the translation vector are t pre ,t gt , then the specific formula for calculating the evaluation index is as follows:

[0217]

[0218] First, in order to verify the effectiveness of the directional description structure and multi-feature constraint MFC (local coherence matching structure, ICP based on normal vector constraint) proposed in this invention in the point cloud registration algorithm, an ablation experiment will be carried out. The four experimental objects are:

[0219] (1) Base algorithm: a point cloud registration algorithm based on the USIP key point extraction module, the local reference frame (LRF) key point feature description module, and ICP feature matching.

[0220] (2) A-Base algorithm: Based on the Base algorithm, this algorithm combines the point cloud registration algorithm with the directional description structure. Other structures are consistent with the Base algorithm.

[0221] (3) B-Base algorithm: Based on the Base algorithm, the feature matching module uses a multi-feature constraint method, and the other structures remain consistent with Base.

[0222] (4) AB-Base algorithm: Based on the Base algorithm, this algorithm combines the directional description structure with the point cloud registration algorithm with multi-feature constraints, which is the complete form of the algorithm of the present invention.

[0223] In the dataset C where the overlap rate is about 50%, the ablation experiment results are shown in the following table:

[0224] Ablation experiment (without noise)

[0225]

[0226] In the dataset D containing noise and with an overlap rate of about 50%, the ablation experiment results are shown in the following table:

[0227] Ablation experiment (including noise)

[0228]

[0229] The test results in the table show that the key point directional description structure and multi-feature constraint module proposed in the present invention show significant performance improvement compared with the Base algorithm when processing point cloud data with an overlap rate of about 50%, regardless of whether the data set contains noise. It not only improves the accuracy of registration, but also shortens the time consumption of algorithm operation.

[0230] At the same time, in order to further verify the superiority of the point cloud registration algorithm proposed in this paper, it is compared with mainstream algorithms such as ICP, Go-ICP, FGR, PointNetLK and DCP. The experiment also uses the C and D datasets, and the results are as follows:

[0231] Comparative test (excluding noise)

[0232]

[0233] Comparative experiment (including noise)

[0234]

[0235]

[0236] The experimental results show that the proposed point cloud registration method based on key point orientation description and multi-feature constraints performs better than mainstream algorithms such as ICP, Go-ICP, FGR, PointNetLK and DCP when dealing with point cloud registration tasks with an overlap rate of about 50%, regardless of whether the data contains noise. All evaluation indicators have been improved to a certain extent, verifying the superiority of the proposed method. Fig. 9 .

Claims

1. A point cloud registration method based on directional description and multi-feature constraints, characterized in that: The following steps are involved: Step 1: Extract point cloud key points through USIP deep learning network; Step 2: Use the directional description structure designed by the present invention to build a DDNet network and perform training; Step 3: Optimize the matching process of key points of source point cloud and target point cloud using local coherence matching method; Step 4: Use ICP optimization based on normal vector constraints to further eliminate false matches; Step 5: Iteratively solve the transformation matrix from the source point cloud to the target point cloud; Step 6: Experimentally prove the superiority of the point cloud registration algorithm proposed in this invention.

2. The point cloud registration method based on directional description and multi-feature constraints according to claim 1, characterized in that: The step 1 specifically includes: Step 11: Point Cloud Perform random transformation T to obtain point cloud Both are input to the FPN layer together; Step 12: The FPN layer first extracts the input point cloud In the example, M marker points are selected according to the farthest point sampling FPS algorithm, denoted as Then each point in the input point cloud X is grouped according to the distance from Node into M groups, and the result is Where K m It represents the number of points in the original point cloud divided into m groups, and then the points divided into each Node are decentralized, that is, get Then all the points are passed to the first fully connected network in groups; Step 13: After passing through the first fully connected network, the features of each point are obtained, and then the features of M nodes can be obtained through the maximum pooling operation Then perform KNN on M Nodes and get Then the coordinate information is decentralized, that is get Then, M nodes correspond to a feature of (C+3, M, K), which is then passed to the second fully connected network; Step 14: After passing through the second fully connected network, the feature (C2, M) is obtained, and it is concat-operated with the feature (C1, M) before KNN in the feature dimension to obtain (C1+ C2,M), then the feature vector is input into MLP to get the output of (4,M), where the first three dimensions represent the relative position information of the key points, and the fourth dimension is transformed by the Softmax function to characterize the uncertainty level of the key point significance, where lower uncertainty values ​​indicate more reliable key feature positioning; Step 15: Input Point Cloud X and After the FPN layer, the final outputs are Q1 and σ1, Q2 and σ2, which represent the points from point cloud X and The coordinates of the key points extracted and the uncertainty corresponding to each key point; USIP designs two loss functions, loss function L p By quantifying the matching degree of key points before and after rotation; the loss function L p as follows: Loss function L c Based on probability, the correct key point position is identified, and an uncertainty value is provided for each identified key point; the loss function L c as follows:

3. The point cloud registration method based on directional description and multi-feature constraints according to claim 1, characterized in that: The step 2 specifically includes: imagining the key point and its neighboring points as being contained in a regular octahedron, and then observing these points from the normal directions of the eight faces of the regular octahedron; Step 21: First select the bottom edge of the face as the x-axis of the new coordinate system, and then use the midpoint of the bottom edge to the vertex as the y-axis of the new coordinate system; then determine the z-axis according to the right-hand rule. First, set the coordinates of the octahedron vertices as follows: The coordinates of vertex a are The coordinates of vertex b are (i, i, 0), the coordinates of vertex c are (i, -i, 0), the coordinates of vertex d are (-i, -i, 0), the coordinates of vertex e are (-i, i, 0), and the coordinates of vertex f are Compute a new coordinate system perpendicular to face ABC where Then the z-axis vector is: After calculating the vectors of the x, y, and z axes, each vector is normalized, so the new coordinate system of the ABC surface is: Step 22: Rotate the key point and its surrounding neighboring points according to the new coordinate system of the octant to obtain eight new point clouds, and then send the new point clouds to the new reasoning network, namely the directional description network.

4. The point cloud registration method based on directional description and multi-feature constraints according to claim 3, characterized in that: In step 2, the directional description network is trained using a triplet loss function. During the training process, the triplet loss function makes the network more robust by grouping anchor samples, similar key points Positive, and non-similar key points Negative and training them simultaneously. The mathematical formula of the triplet loss function is as follows: d i (a,p)=f(a)-f(p),d i (a,n)=f(a)-f(n) Where a represents the anchor point, the input sample anchor, p represents the positive sample, n represents the negative sample, d i (a,p) represents the distance difference between f(a) and f(p) after a and p are described by network features. Similarly, d i (a,n) represents the distance difference between f(a) and f(n) after a and n are described by network features, N represents the number of sample groups, τ is a hyperparameter, which is set to 1 as the cutoff value, and L is the loss function; During the training process, anchors, positive samples and negative samples are grouped into triplets; the triplets are simultaneously fed into the directional description network DDNet; for each group of samples, DDNet can extract three (1,32)-dimensional feature vectors.

5. The point cloud registration method based on directional description and multi-feature constraints according to claim 1, characterized in that: The step 3 specifically comprises the following steps: Step 31: Traverse each key point in the source point cloud and find the four other key points closest to the key point; then, select two from the four neighboring key points to form three points p0, p1, p2 with the key point, and find the corresponding key points q0, q1, q2 in the target point cloud according to the feature description, that is, a triplet; Step 32: Determine whether the triangle formed by the three key points p0, p1, and p2 in the source point cloud is roughly the same as the triangle formed by the key points q0, q1, and q2 corresponding to the target point cloud. The formula for determining whether the triangle shape meets the conditions is as follows: Step 33: The qualified triples will calculate the minimum circumscribed circle according to the three key points of the source point cloud and the target point cloud, and construct two spheres with the center of each circumscribed circle as the center of the sphere and the radius of the circumscribed circle as the radius of the sphere. The calculation method for constructing the sphere is as follows: First, constrain the three points to be at the same distance from the center of the sphere: Then solve jointly and eliminate the variable R: Simplify the formula to obtain the linear equation system about the coordinates of the sphere center: Solve the linear equations to obtain the coordinates of the center of the sphere and the radius: The two spheres obtained will cover the points in a certain area of ​​the source point cloud and the target point cloud respectively, and the points in each sphere will replace the source point cloud and the target point cloud and be sent to ICP for solution.

6. The point cloud registration method based on directional description and multi-feature constraints according to claim 1, characterized in that: The step 4 is specifically as follows: finding the rotation matrix R and the translation vector t. The mathematical model formula of point cloud registration is as follows: Among them, q i is a point in the target point cloud, p i Q i The corresponding point in the source point cloud; Define the centroids of the source and target point clouds: Simplify the optimization terms of the mathematical model: Optimize the project i -qR(p i -p) is zero after summing, then the function simplifies to:

7. The point cloud registration method based on directional description and multi-feature constraints according to claim 6, characterized in that: In step 4, the centroid coordinates of each point in the source point cloud and the target point cloud are first calculated: q i * =q i -q,p i * =p i -p Then optimize the first error term to calculate R: Expand the error term of the rotation matrix: Where R T R = I and the first error is independent of the rotation matrix, then the actual optimization function is: Define the matrix Matrix: Perform SVD decomposition on the matrix Matrix: Matrix=UΣV When the original matrix Matrix is ​​full rank, it means that the original matrix has N non-zero singular values, so the diagonal matrix Σ will have N non-zero elements; in this case, U and V represent the standard orthogonal basis of the original data space and the target data space, respectively, and Σ represents the mapping relationship between the two spaces; then the required R is: R=UV T Finally, substitute the obtained R into the second term of the error and calculate t: t=q-Rp 8. The point cloud registration method based on directional description and multi-feature constraints according to claim 7, characterized in that: In step 4, three criteria are used to eliminate false matches in the ICP process: (1) Point distance exceeds threshold: When the distance between a point in the source point cloud and a point in the target point cloud exceeds a certain set threshold, these point pairs will be eliminated; (2) Curvature approximation: If the curvature difference between two points is too large, these point pairs will be eliminated; (3) Normal vector direction approximation: If the normal vector directions of two matching points are significantly different, such point pairs will also be considered as incorrect matches and removed; According to the above formulas such as Matrix = UΣV, three eigenvalues ​​and corresponding eigenvectors are decomposed, where the eigenvalues ​​are λ0, λ1, λ2, λ0<λ1<λ2, and the normal vector is the eigenvector corresponding to λ0; the calculation formula of curvature δ is as follows: The ICP mathematical model based on normal vector constraints is derived as follows. First, the source point cloud is defined as The target point cloud is Then expand the midpoints of each point cloud, that is, the normal vector: After the source point cloud is transformed by the matrix T(R, t): The source point cloud is transformed into coordinates, and then the error between the point in the source point cloud and the corresponding point is defined as e ij , the calculation formula is as follows: The final mathematical model is defined as: At this point, the ICP based on the normal vector constraint becomes a nonlinear least squares problem like the traditional method, and the LM method is used to solve it as follows: (H+λI)ΔT=b H=∑J i T J i 9. The point cloud registration method based on directional description and multi-feature constraints according to claim 1, characterized in that: In step 5, first, the point cloud data of the source point cloud is set as: The key point is After the directional descriptor, the feature corresponding to the key point is described as Set the point cloud data corresponding to the target point cloud as: The key point is The feature description corresponding to the key point is Set the number of inliers inline to 0 and the inlier rate ε threshold to 0.75 to determine whether the source point cloud meets the output conditions after the transformation matrix, set the optimal inlier rate ε' to 0, the optimal transformation matrix R', and the optimal translation vector t'; 1) Traverse the key points of the source point cloud, recorded as The four key points closest to this key point are The corresponding feature description is 2) The four adjacent key points and their corresponding feature descriptions are input into the local coherence matching structure for judgment. If the shape requirements are met, the point cloud of the smallest circumscribed sphere area of ​​the source point cloud divided by the three key points is output one by one as the new source point cloud. And output the point cloud of the minimum circumscribed sphere area divided by the corresponding three key points as the new target point cloud If none of them meet the requirements, return to the first step and continue to traverse the key points; 3) Input P' and Q' into the ICP mathematical model based on normal vector constraints and construct a matrix, then calculate the rotation matrix R and translation vector t from the key points of the source point cloud to the key points of the target point cloud; then transform the points in the input source point cloud P and rotate them to the coordinate system of the target point cloud Q; 4) For all transformed points in P, find the closest point in the target point cloud and calculate the Euclidean distance d between them. a , if d a <0.05, then inline+=1, after the judgment, calculate the inlier rate ε=inline / N1, then judge the inlier rate ε>0.4, if it is true, output the optimal solution R and t, otherwise judge ε'<ε, if it is true, update ε', R', t', if not, judge the key point Is the sequence number i greater than the number of key points n of the source point cloud? If so, output R', t', otherwise continue to traverse the key points of the source point cloud.

10. The point cloud registration method based on orientation description and multi-feature constraints according to claim 1, characterized in that: In step 6, the test set in the ModelNet40 standard data set is tested; certain preprocessing is performed on the experimental test set to simulate the overlap rate and sensor noise in actual working conditions, as follows: Step 61: Conduct experiments on data with an overlap rate of 50%; the source point cloud refers to 25% of the original data in the test set randomly cropped, and the target point cloud refers to another 25% of the original data in the test set randomly cropped, and ensure that this part of the data does not overlap with the cropped part of the source point cloud to achieve an overlap rate of about 50%, and then randomly rotate the target point cloud x, y, and z axes in the range of [-90°, 90°], and randomly translate in the range of [-1m, 1m]; at the same time, save the rotation matrix R and translation vector t of each set of source point clouds and the corresponding target point cloud as the real transformation parameters; This set of experimental data sets is collectively referred to as data set C below; Step 62: Experiments are conducted on data with an overlap rate of about 50% and containing noise. In this set of experimental data sets, both the source point cloud and the target point cloud are based on the point cloud data in data set C, and Gaussian noise with a mean of 0 and a variance of 0.01 is added to fit the actual sensor noise. This set of experimental data sets is collectively referred to as data set D below; The indicators for experimental evaluation of point cloud registration effects include: mean square error RMSE R and RMSE t , root mean square error RMSE R and RMSE t , mean absolute error MAE R and MAE t And the average time Time required to complete the registration; Assume that the calculated value and the true value of the Euler angle are θ pre ,θ gt , the calculated value and true value of the translation vector are t pre ,t gt , then the specific formula for calculating the evaluation index is as follows:

Images