Point cloud registration method based on improved SHOT descriptor

Abstract

The invention discloses a point cloud registration method based on an improved SHOT descriptor, and aims to solve the defects of a traditional point cloud SHOT descriptor in the aspects of sparse point cloud or high-density point cloud feature extraction, descriptor and matching optimization. The method comprises the following steps: preprocessing a source point cloud set and a target point cloud set, calculating an average point spacing and extracting feature points; the radius of a descriptor support region is adjusted by introducing interpolation, so that the robustness of key point features is improved; then, a weighted covariance matrix is calculated, eigenvalue decomposition is carried out, and an accurate local coordinate system and normal information are obtained; sHOT descriptor matching point pairs based on the radius of the dynamic support area introduces spectrum matching to be combined with ICP, and error point pairs are eliminated. And the registration precision is further improved. Compared with the prior art, the method has the advantages that the descriptive property of the descriptor structure and the optimization of the matching process are enhanced through the dynamic adjustment of the support area, the precision, robustness and real-time property of point cloud registration are improved, and the method is suitable for the fields of computer vision, three-dimensional point cloud registration and the like.

Description

Technical Field

[0001] It belongs to the field of computer vision and 3D point cloud processing technology, specifically focusing on the improvement of point cloud descriptor algorithms. Background Art

[0002] Point cloud registration is a key technology in computer vision and 3D data processing. SHOT descriptors, which describe local features using directional histograms, offer good robustness. However, they suffer from practical limitations: a fixed support region radius does not adapt to changes in point cloud density, and global optimization matching is computationally complex and exhibits low accuracy. To address these issues, this paper proposes an improved SHOT descriptor point cloud registration method that improves registration accuracy and robustness by dynamically adjusting the support region, weighting the local coordinate system, and introducing spectral matching optimization. Summary of the Invention

[0003] The purpose of this invention is to provide a point cloud registration method based on an improved SHOT descriptor to overcome the shortcomings of the traditional SHOT method in terms of support region adaptability, descriptor discrimination, and matching efficiency, thereby improving the accuracy, robustness, and real-time performance of point cloud registration. To achieve the above objectives, the present invention adopts the following technical solutions:

[0004] A point cloud registration method based on an improved SHOT descriptor, characterized by comprising the following steps:

[0005] Step 1: Preprocessing

[0006] The two point cloud datasets to be registered are the source point cloud set Q and the target point cloud set P. The two point clouds are preprocessed and the average point distance d is calculated. avg , extract the feature point set Q={q1,q2,...,q m},P={p1,p2,...,p n};

[0007] Step 2:

[0008] Descriptor support region definition

[0009] For each feature point q i ∈Q and p i ∈P, define the support area, the support area is the feature point q i and p i is the center of the sphere, and the radius is r i and R i The support area is defined by dynamically adjusting the radius of the sphere.

[0010] 1. Calculate the initial radius:

[0011] Local density d of feature points i

[0012]

[0013] in:

[0014] N is point q i As the center, the radius is the initial value r i The number of neighboring points within r i The default value is the average point distance d of the point cloud. avg A multiple of .

[0015] q j It's Q i Neighborhood points.

[0016] ||q i -q j ||It's point q i and dot q j The Euclidean distance between them.

[0017] Support area radius r i Dynamic adjustment based on the distribution of feature points, the formula is:

[0018]

[0019] in:

[0020] r min =k1d avg ;

[0021] r max =k2d avg ;

[0022] d min =k3d avg ;

[0023] d max =k4d avg .

[0024] (k1, k2, k3, k4 are preset coefficients).

[0025] Define the support point set:

[0026] S i =||S k -q i ||≤r i ,(S k ∈Q)

[0027] For the target point cloud P, R i The calculation method of r i Same as above, and redefine the support point set:

[0028] S j =||S l -p j ||≤R i ,(S l ∈P)

[0029] S i and S j are the support point sets in the source point cloud and target point cloud respectively.

[0030] q i is a feature point in the source point cloud,

[0031] p j is a feature point in the target point cloud.

[0032] ||S k -q i || is the feature point q in the source point cloud i To support point S k The Euclidean distance between them.

[0033] r i and R i is the radius of the support region, which represents the range of the support point set.

[0034] Step 3:

[0035] Calculate the normal for each set of support points

[0036]

[0037] Among them, p i is a neighborhood point, p - is the neighborhood centroid, k is the number of neighboring points. The normal vector is the eigenvector corresponding to the smallest eigenvalue of C.

[0038] Step 4:

[0039] Local coordinate system and descriptor construction

[0040] In the dynamic radius r i Collect neighborhood points and calculate the centroid:

[0041]

[0042] For each support point set S i , calculate the weighted covariance matrix M i :covariance matrix M i Perform eigenvalue decomposition to obtain eigenvalues ​​λ1, λ2, λ3 (satisfying λ1≥λ2≥λ3), and obtain the corresponding eigenvectors v1, v2, v3. Through eigenvalue decomposition, the principal axis direction of the local coordinate system is obtained. If d k ≥ri , the weight is 0 to ensure the physical meaning of the weighted covariance matrix

[0043]

[0044] Among them, the weight w k =max(0,r i -||S k -c i ||),S k are the coordinates of the neighborhood points.

[0045] Step 5: Compare the SHOT descriptors of the key points of the source point cloud and the target point cloud to find the corresponding relationship:

[0046] For each source keypoint descriptor, calculate the similarity with all target keypoint descriptors:

[0047]

[0048] Among them, desc q ,desc p is a 352-dimensional descriptor vector.

[0049] Step 6: Spectral matching refinement

[0050] Filter geometrically inconsistent matches through spectral matching. Analyze the geometric constraints between matching pairs, construct a consistency matrix, and use spectral analysis to select reliable matches, improving the robustness of subsequent transformation estimation.

[0051] Construct a consistency matrix A of size n×n (n is the number of matching pairs). For each pair of matching (i, j):

[0052] dist Q =||q i -q j ||

[0053] Similarly, calculate the distance between two points in the target point cloud:

[0054] dist P =||p i -p j ||

[0055] Calculate the distance difference:

[0056] δ=|dist Q -dist P ∣

[0057] Consistency value:

[0058]

[0059] Eigen-decompose A and select the eigenvector v corresponding to the largest eigenvalue. Set threshold = 0.5 x max(v), and keep the matches with v(i) > threshold to form the refined match list rematches.

[0060] Step 7:

[0061] SVD estimates initial transformation

[0062] Based on the refined matches, compute the initial rigid body transformation using Singular Value Decomposition (SVD):

[0063] Construct point set matrices, Q = {q1, q2,..., q m} ∈ R 3*m , P = {p1, p2,..., p n} ∈ R 3*m

[0064] Compute centroids:

[0065]

[0066] Center the point sets:

[0067]

[0068] Compute the covariance matrix:

[0069] H = Q' · P' T

[0070] SVD decompose H:

[0071] H = U∑V T

[0072] Compute the rotation matrix R and translation vector t:

[0073]

[0074] Final rigid body transformation matrix:

[0075]

[0076] Optimize using Iterative Closest Point (ICP) algorithm, based on the initial point pairs, optimize the transformation matrix T between the point clouds, by minimizing the following objective function to obtain the final transformation:

[0077] E = ∑ i w i ||q i - T(p i ) || 2

[0078] Compared with the prior art, the advantages of the present invention are:

[0079] By dynamically adjusting the radius of the support area, we can adapt to changes in point cloud density and improve the robustness of feature extraction.

[0080] Use weighted covariance matrix to improve local coordinate system accuracy;

[0081] Spectral matching is introduced and combined with ICP and SE-ICP algorithms to eliminate erroneous key points, thereby improving matching accuracy and efficiency.

[0082] The present invention is applicable to fields such as computer vision and point cloud registration. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 Schematic diagram of the process of the improved registration method based on SHOT descriptor of the present invention.

[0084] Figure 2 Schematic diagram of dynamic adjustment of the key point support area of ​​the SHOT descriptor of the present invention.

[0085] Figure 3 Schematic diagram of the similarity matching process of the present invention. DETAILED DESCRIPTION

[0086] Step 1: Point cloud data acquisition and feature point extraction

[0087] The workpiece point cloud data is obtained by a 3D scanner, and two point cloud sets to be registered are obtained, namely the source point cloud set Q and the target point cloud set P. The intrinsic shape signature (ISS) method is used to extract the feature point set Q = {q1,q2,...,q m},P={p1,p2,...,p n}; Calculate the average point spacing d avg ,

[0088] Step 2:

[0089] Descriptor support region definition

[0090] For each feature point q i ∈Q and p i ∈P, define the support area, the support area is the feature point q i and p i is the center of the sphere, and the radius is r i and R i The support area is defined by dynamically adjusting the radius of the sphere.

[0091] 1. Calculate the initial radius:

[0092] Local density d of feature points i

[0093]

[0094] in:

[0095] N is point q i As the center, the radius is the initial value r i The number of neighboring points within r i The default value is the average point distance d of the point cloud. avg A multiple of .

[0096] q j It's Q i 's neighborhood points.

[0097] ||q i -q j ||It's point q i and dot q j The Euclidean distance between .

[0098] Support area radius r i Dynamic adjustment based on the distribution of feature points, the formula is:

[0099]

[0100] in:

[0101] r min =k1d avg ;

[0102] r max =k2d avg ;

[0103] d min =k3d avg ;

[0104] d max =k4d avg .

[0105] (k1, k2, k3, k4 are preset coefficients). (Initial values ​​k1 = 5, k2 = 20, k3 = 0.5, k4 = 1.5)

[0106] Define the support point set:

[0107] S i =||S k -q i ||≤r i ,(S k ∈Q)

[0108] For the target point cloud P, R i The calculation method of r i Same as above, and redefine the support point set:

[0109] S j =||S l -p j ||≤R i ,(S l ∈P)

[0110] S i and S j are the support point sets in the source point cloud and target point cloud respectively.

[0111] q i is a feature point in the source point cloud,

[0112] p j is a feature point in the target point cloud.

[0113] ||S k -q i || is the feature point q in the source point cloud i To support point S k The Euclidean distance between .

[0114] r i and R i is the radius of the support region, which represents the range of the support point set.

[0115] Step 3:

[0116] Calculate the normal for each set of support points

[0117]

[0118] Among them, p i is a neighborhood point, p - is the neighborhood centroid, k is the number of neighboring points. The normal vector is the eigenvector corresponding to the smallest eigenvalue of C.

[0119] Step 4:

[0120] Local coordinate system and descriptor construction

[0121] In the dynamic radius r i Collect neighborhood points and calculate the centroid:

[0122]

[0123] For each support point set S i , calculate the weighted covariance matrix M i :covariance matrix M iPerform eigenvalue decomposition to obtain eigenvalues ​​λ1, λ2, λ3 (satisfying λ1≥λ2≥λ3), and obtain the corresponding eigenvectors v1, v2, v3. Through eigenvalue decomposition, the principal axis direction of the local coordinate system is obtained. If d k ≥r i , the weight is 0 to ensure the physical meaning of the weighted covariance matrix

[0124]

[0125] Among them, the weight w k =max(0,r i -||S k -c i ||),S k are the coordinates of the neighborhood points.

[0126] Covariance matrix M i Perform eigenvalue decomposition to obtain its eigenvalues ​​and eigenvectors. Eigenvalues ​​λ1, λ2, λ3, and corresponding eigenvectors v1, v2, v3, through eigenvalue decomposition, the covariance matrix M i The eigenvectors represent the main directions of the local point cloud distribution, and the eigenvalues ​​reflect the variance of these directions. To further construct the local coordinate system, the eigenvectors v1, v2, and v3 corresponding to the eigenvalues ​​are selected, and the local coordinate system is constructed according to the order of their eigenvalues ​​(λ1 corresponds to the maximum principal axis, λ2 and λ3 correspond to the secondary principal axis and the minimum principal axis, respectively). Since the direction of the eigenvector may be affected by the positive and negative signs, resulting in uncertainty in the direction, a voting mechanism needs to be introduced to eliminate the ambiguity of the coordinate axis direction. The core of the voting mechanism is to use the normal information of the support point set to check the direction consistency. The specific method is to calculate the dot product of the normal of each support point and the eigenvector. If the dot product is positive, it indicates that the direction of the eigenvector is correct; if the dot product is negative, the inconsistency is eliminated by reversing the direction of the eigenvector.

[0127] Determine the three principal axes x, y, and z of the local coordinate system; define the positive direction of the eigenvector as x+, y+, and z+, and the negative direction as x-, y-, and z-;

[0128] Voting on the x-axis, define the set:

[0129]

[0130] The y and z axes are voted in the same way to determine their directions, and finally the local coordinate system (x, y, z) is obtained. In this way, a consistent and accurate local coordinate system is finally obtained, where v1, v2, and v3 are the three main axes of the local coordinate system, forming an orthogonal coordinate system. The points of the support point set are projected according to the local coordinate system, and their coordinates are projected into the local coordinate system through the following process:

[0131]

[0132] By projection, the coordinates of the support point skare transformed to the local coordinate system, resulting in new coordinates (xk',yk',zk'). These three coordinates correspond to the three coordinate axes in the local coordinate system. Based on the key point, its normal, neighborhood points and LRF, a 352-dimensional SHOT descriptor is computed, describing the histogram distribution of the local surface normal.

[0133] Step 5: Compare the SHOT descriptors of the key points of the source and target point clouds to find the correspondence:

[0134] For each source key point descriptor, the similarity with all target key point descriptors is computed:

[0135]

[0136] where desc q ,desc p are the 352-dimensional descriptor vectors.

[0137] The best and second-best matches are recorded in terms of similarity and corresponding indices. If the best similarity sim≥0.9 and the ratio of the similarities of the best and second-best matches is greater than 0.8, the match is accepted, forming a list of matched pairs matches. This preliminary set of matches can still contain false matches, so further optimization is needed to make the final set of matched pairs m more reliable.

[0138] Step 6: Spectral matching refinement

[0139] Geometrically inconsistent matches are filtered out by spectral matching. By analyzing the geometric constraints between the matched pairs, a consistency matrix is constructed and spectral analysis is used to select reliable matches, improving the robustness of subsequent transformation estimation.

[0140] A consistency matrix A is constructed, with size n x n (n is the number of matched pairs). For each matched pair (i,j):

[0141] dist Q = ||q i -q j ||

[0142] Similarly, the distance between two points in the target point cloud is calculated:

[0143] dist P = ||p i -p j ||

[0144] The distance difference is calculated:

[0145] δ = |dist P - dist P |Q -dist P ∣

[0146] Consistency value:

[0147]

[0148] Perform eigendecomposition on A and select the eigenvector v corresponding to the maximum eigenvalue. Set the threshold threshold = 0.5 × max(v)threshold = 0.5 × max(v) and retain the matching pairs with v(i) ≥ threshold to form the refined matching list rematches.

[0149] Step 7:

[0150] SVD estimates the initial transformation

[0151] Based on the refined matching pairs, singular value decomposition (SVD) is used to calculate the initial rigid body transformation:

[0152] Construct a point set matrix, Q={q1,q2,...,q m}∈R 3*m , P={p1,p2,...,p n}R 3*m

[0153] Calculate the centroid:

[0154]

[0155] Centralized point set:

[0156]

[0157] Compute the covariance matrix:

[0158] H=Q'·P' T

[0159] Perform SVD decomposition on H:

[0160] H=U∑V T

[0161] Calculate the rotation matrix R and translation vector t:

[0162]

[0163] Final rigid body transformation matrix:

[0164]

[0165] The Iterative Closest Point (ICP) algorithm is used for optimization. Based on the initial point pairs, the transformation matrix T between the point clouds is optimized, and the final transformation is obtained by minimizing the following objective function:

[0166] E=∑ i w i ||q i -T(p i )|| 2

[0167] Compared to traditional SHOT descriptor registration methods, this method improves the adaptability and descriptiveness of key points in point clouds of varying density by utilizing a dynamic descriptor radius. It also introduces spectral matching to improve the accuracy of key point pair matching between source and target point clouds. This method meets the requirements of high precision and robustness, achieving stable and accurate point cloud registration in various scenarios and providing a new approach to point cloud registration technology.

Claims

1. A point cloud registration method based on an improved SHOT descriptor, comprising the following steps: Step 1: Preprocessing The two point cloud datasets to be registered are the source point cloud set Q and the target point cloud set P. The two point clouds are preprocessed and the average point distance d is calculated. avg , extract the feature point set Q={q1,q2,...,q m},P={p1,p2,...,p n}; Step 2: Descriptor support region definition For each feature point q i ∈Q and p i ∈P, define the support area, the support area is the feature point q i and p i is the center of the sphere, and the radius is r i and R i The support area is defined by dynamically adjusting the radius of the sphere.

1. Calculate the initial radius: Local density d of feature points i in: N is point q i As the center, the radius is the initial value r i The number of neighboring points within r i The default value is the average point distance d of the point cloud. avg A multiple of . q j It's Q i 's neighborhood points. ||q i -q j ||It's point q i and dot q j The Euclidean distance between them. Support area radius r i Dynamic adjustment based on the distribution of feature points, the formula is: in: r min =k1d avg ; r max =k2d avg ; d min =k3d avg ; d max =k4d avg . (k1, k2, k3, k4 are preset coefficients). Define the support point set: S i =||S k -q i ||≤r i ,(S k ∈Q) For the target point cloud P, R i The calculation method of r i Same as above, and redefine the support point set: S j =||S l -p j ||≤R i ,(S l ∈P) S i and S j are the support point sets in the source point cloud and target point cloud respectively. q i is a feature point in the source point cloud, p j is a feature point in the target point cloud. ||S k -q i || is the feature point q in the source point cloud i To support point S k The Euclidean distance between them. r i and R i is the radius of the support region, which represents the range of the support point set. Step 3: Calculate the normal for each set of support points Among them, p i is a neighborhood point, p - is the neighborhood centroid, k is the number of neighboring points. The normal vector is the eigenvector corresponding to the smallest eigenvalue of C. Step 4: Local coordinate system and descriptor construction In the dynamic radius r i Collect neighborhood points and calculate the centroid: For each support point set S i , calculate the weighted covariance matrix M i :covariance matrix M i Perform eigenvalue decomposition to obtain eigenvalues ​​λ1, λ2, λ3 (satisfying λ1≥λ2≥λ3), and obtain the corresponding eigenvectors v1, v2, v3. Through eigenvalue decomposition, the principal axis direction of the local coordinate system is obtained. If d k ≥r i , the weight is 0 to ensure the physical meaning of the weighted covariance matrix Among them, the weight S k are the coordinates of the neighborhood points. Step 5: Compare the SHOT descriptors of the key points of the source point cloud and the target point cloud to find the corresponding relationship: For each source keypoint descriptor, calculate the similarity with all target keypoint descriptors: Among them, desc q ,desc p is a 352-dimensional descriptor vector. Step 6: Spectral matching refinement Filter geometrically inconsistent matches through spectral matching. Analyze the geometric constraints between matching pairs, construct a consistency matrix, and use spectral analysis to select reliable matches, improving the robustness of subsequent transformation estimation. Construct a consistency matrix A of size n×n (n is the number of matching pairs). For each pair of matching (i, j): dist Q =||q i -q j || Similarly, calculate the distance between two points in the target point cloud: dist P =||p i -p j || Calculate the distance difference: δ=∣dist Q -dist P ∣ Consistency value: Perform eigendecomposition on A and select the eigenvector v corresponding to the maximum eigenvalue. Set the threshold threshold = 0.5 × max(v)threshold = 0.5 × max(v) and retain the matching pairs with v(i) ≥ threshold to form the refined matching list rematches. Step 7: SVD estimates the initial transformation Based on the refined matching pairs, the initial rigid body transformation is calculated using singular value decomposition (SVD): the source and target keypoint coordinates are extracted from the matching pairs to form 3×m matrices Q and P respectively, where m is the number of matching pairs. The Iterative Closest Point (ICP) algorithm is used for optimization. Based on the initial point pairs, the transformation matrix T between the point clouds is optimized, and the final transformation is obtained by minimizing the following objective function: E=∑ i w i ||q i -T(p i )|| 2 。 2. The method according to claim 1, characterized in that The feature point extraction adopts the ISS algorithm.

Images