A multi-level feature fusion adaptive point cloud registration method

Abstract

The present invention discloses a multi-level feature fusion adaptive point cloud registration method, which relates to the field of point cloud processing technology and industrial visual inspection technology, comprising: preprocessing a point cloud pair; performing multi-scale feature extraction on the preprocessed point cloud pair, constructing a feature pyramid through the extracted multi-scale features, and adopting a top-down guidance mechanism between each level of the feature pyramid; performing preliminary feature fusion on the feature points of the point cloud pair based on the extracted multi-scale features; performing multi-level registration on the point cloud pair; the multi-level includes first-level registration, second-level registration and third-level registration, and after each level of registration is completed, a transformation matrix of the layer is obtained, and the point cloud to be registered is registered through the transformation matrix of each layer; the present invention adopts a multi-scale feature extraction strategy, and selects feature descriptors of three different scales, namely ESF, SHOT and FPFH feature descriptors, to address the degradation problem in point cloud registration, so as to realize feature extraction and parameter adaptive adjustment from macro to micro.

Description

Technical Field

[0001] The present invention relates to point cloud processing technology and industrial visual inspection technology, and in particular to a multi-level feature fusion adaptive point cloud registration method. Background Art

[0002] Point cloud registration technology is a key step in 3D data processing and has important applications in industrial silo monitoring and volume measurement. For example, the loading and unloading of silos is a continuous process. Measuring the difference between two point clouds in a silo is difficult due to the influence of different external environments (such as the effects of temperature, humidity, and dust on radar scanning, and the effect of physical vibrations on the point cloud angle caused by slight offsets in hardware equipment). It is difficult to maintain consistency between the current scan state and the previous scan state. Therefore, it is necessary to select appropriate features to align the point clouds, thereby combining the difference between the historical point cloud and the current point cloud to calculate the difference.

[0003] With the rapid development of LiDAR technology, continuous monitoring of material pile morphology is now possible using fixed scanning devices mounted on top of material stacks. However, in real industrial environments, acquiring and processing point cloud data still faces numerous challenges. In practical applications, point cloud data quality is often degraded by a variety of factors. First, non-rigid transformations caused by device jitter and temperature drift can lead to deviations between consecutively acquired point clouds. Second, noise interference from factors such as dust and moisture in the industrial environment. Furthermore, due to the irregularities of the material pile surface and the dust generated by falling materials, point cloud data often exhibits quality issues such as uneven density and localized point loss. Existing point cloud registration methods have significant shortcomings when processing industrial material stack scenarios. While widely used, traditional ICP algorithms are prone to falling into local optima when handling large-scale deformations. Methods based on features such as FPFH are unstable in the presence of significant noise and deformation. Most existing methods adopt fixed parameter strategies and lack adaptive mechanisms for point cloud quality, making them difficult to address the complex and changing data quality issues in industrial environments. These problems seriously restrict the application of point cloud registration technology in industrial material tower monitoring. Therefore, it is of great engineering significance to develop a registration method that can effectively handle point cloud degradation, adapt to data quality, and have high accuracy. Summary of the Invention

[0004] To solve the technical problems raised in the background technology, the present invention proposes a multi-level feature fusion adaptive point cloud registration method, comprising:

[0005] S1. Preprocess the point cloud pair;

[0006] S2. Perform multi-scale feature extraction on the pre-processed point cloud pairs, and construct a feature pyramid based on the extracted multi-scale features. A top-down guidance mechanism is adopted between the levels of the feature pyramid.

[0007] S3. Perform preliminary feature fusion on the feature points of the point cloud pair based on the extracted multi-scale features;

[0008] S4. Perform multi-level registration on the point cloud pairs; the multi-level registration includes the first-level registration, the second-level registration, and the third-level registration. After each level of registration is completed, the transformation matrix of the level is obtained, and the point cloud to be registered is registered using the transformation matrix of each level.

[0009] Specifically, the S1 includes the following steps:

[0010] S1.1 performs noise reduction on the point cloud and identifies and removes abnormal points by analyzing the neighborhood distribution characteristics of the points;

[0011] S1.2 uses voxel grid filtering algorithm to downsample the point cloud;

[0012] S1.3 calculates the average point spacing d of the downsampled point cloud and uses it as a benchmark parameter for feature extraction;

[0013] S1.4 builds an octree-based spatial index structure for neighborhood search.

[0014] Specifically, the S2 includes the following steps:

[0015] S2.1. Use ESF features to extract first-scale features and obtain the ESF feature descriptor ESF(P). Using the first-scale features as the top layer of the pyramid, perform tower structure recognition and region segmentation on the preprocessed point cloud. Cylindrical segments, conical segments, and transition segments are obtained, and the rationality of the region segmentation is verified using ESF features.

[0016] S2.2. Use the SHOT feature to extract second-scale features, obtaining the SHOT feature descriptor SHOT(p). The second-scale feature serves as the middle layer of the pyramid, used to extract local geometric structure and transmit information from the middle layer of the feature pyramid to the bottom layer, providing guidance for the bottom layer.

[0017] S2.3. Use the FPFH feature to extract third-scale features and obtain the FPFH feature descriptor FPFH(p). The third-scale feature serves as the lower layer of the pyramid to extract fine geometric structures.

[0018] For each point cloud, its feature point set is its FPFH feature point set, and the feature points in subsequent S3 and S4 are all feature points of this FPFH feature point set.

[0019] Specifically, the S2.1 includes the following steps:

[0020] S2.1.1 Perform cylindrical segment detection on the point cloud, extract the cylinder radius r and axial direction vector v, and use the axial direction as the z-axis direction;

[0021] S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud, and perform reliability verification on the cylindrical segment through the ESF feature. If the verification passes, record the cylindrical parameters, cylindrical radius r, and axial direction vector v, and mark the points in the point cloud that meet the cylindrical model features as cylindrical segments, and enter S2.1.4; if the verification fails, enter S2.1.3; the reliability verification evaluation criteria include two aspects: (1) whether the deviation between the D2 distribution of the ESF feature and the theoretical cylinder established based on the cylindrical radius r extracted by RANSAC is within the range of ε1, and (2) whether the roundness deviation between the A3 component of the ESF feature and the theoretical cylindrical cross section is within the range of ε2. Both aspects must be met at the same time for the verification to pass.

[0022] S2.1.3 Crop the point cloud along the z-axis, removing δh × H points near the highest and lowest points, where H is the total height of the point cloud and δh is the cropping ratio. After cropping, return to step S2.1.1 and accumulate the number of cropping attempts. If verification fails after reaching the maximum number of cropping attempts, the point cloud quality is insufficient for reliable cylindrical feature extraction, and registration is aborted.

[0023] S2.1.4. Based on the cylindrical segments obtained by segmentation, further divide the point cloud into cylindrical segments, transition segments, and cone segments. The transition segments are the area within the cylindrical segment boundary (-tr, tr), and the cone segments are the remaining area.

[0024] Specifically, the S2.2 includes the following steps:

[0025] S2.2.1. Define the structure indicator function I(p),

[0026]

[0027] Among them, p represents the point in the point cloud that has been divided, β is the sampling correction coefficient of the cone segment, and μ is the sampling density coefficient of the transition segment, β∈(1,2), μ∈(1,2);

[0028] S2.2.2. Calculate the basic sampling interval,

[0029] ds_base(p) = d / ρ(p)

[0030] ρ(p) = α·I(p)

[0031] Where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the cone segment sampling correction coefficient, and ρ(p) is the sampling density function;

[0032] S2.2.3. Obtain a local quality score Q(p),

[0033] Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3

[0034] Q_d(p)= min(1, |N(p)| / N_exp)

[0035] Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π)

[0036] Q_r(p) = exp(-σ / d)

[0037] Among them, Q_d(p) is the local point density score of point p, N(p) is the local neighborhood point set of point p, and for any point p in the point cloud, its neighborhood N(p) is defined as the set of all points in the spherical space with p as the center and radius r.

[0038] N(p) = {q | ||q - p|| ≤ r}

[0039] Where r is the neighborhood radius, ||·|| represents the Euclidean distance, N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood points and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of the local plane fitting based on principal component analysis.

[0040] S2.2.4. Adaptive sampling spacing ds(p) is obtained based on the local quality score.

[0041] ds(p) = ds_base(p)·(1 + γ(1-Q(p)))

[0042] Where γ is the sampling adjustment coefficient, Q(p)∈[0,1]. When Q(p)=1, the sampling interval is maintained. When Q(p) decreases, the sampling interval is gradually increased by 1+γ(1-Q(p)).

[0043] S2.2.5. Construct the SHOT local reference frame;

[0044] For the cylindrical segment and the transition segment on the cylindrical structure, the identified cylindrical axis is used as the main direction of the reference frame, and the secondary direction is determined by combining the normal vector of the local neighborhood point set of each point on the cylindrical segment point cloud;

[0045] For the portion of the tapered segment and transition segment on the conical structure, the principal direction of the local surface is calculated and the principal component analysis method is used to determine the complete reference frame. Specifically, the local neighborhood point set of each point on the tapered segment point cloud is first obtained, and principal component analysis is performed on this point set. The eigenvector corresponding to the largest eigenvalue is used as the principal direction of the reference frame, the eigenvector corresponding to the second largest eigenvalue is used as the secondary direction, and the eigenvector corresponding to the smallest eigenvalue is used as the third direction.

[0046] S2.2.6. Perform SHOT feature point extraction on the point cloud. For each point p in the point cloud, sort the points in the point cloud based on its local quality score Q(p). Traverse the sorted point list and find the neighborhood N1(p) of the current point p, where N1(p) = {q | ||q- p|| ≤ ds(p)}

[0047] Where ds(p) is the adaptive sampling spacing; if there is no added SHOT feature point in the neighborhood N1(p) of the current point p, then p is added to the SHOT feature point set;

[0048] S2.2.7. Calculate the SHOT feature descriptor SHOT(p) and calculate the feature response value S(p) for each point p in the point cloud SHOT feature point set, including the following steps:

[0049] (1) Based on the local reference frame established in S2.2.5, the local spherical support region of point p is divided into 32 spatial grids. The normal vector distribution histogram of 11 bins is calculated in each grid to obtain a 352-dimensional SHOT feature vector. This SHOT feature vector is the SHOT feature descriptor SHOT(p);

[0050] (2) Reorganize the eigenvector into a 32×11 matrix M with 11 elements per row to calculate the characteristic response value S(p).

[0051] S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|)

[0052] Where M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i,j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation;

[0053] Specifically, the above 2.3 includes the following steps:

[0054] S2.3.1. Sort the points in the SHOT feature point set by their feature response values ​​S(p), traverse the sorted point list, and check the neighborhood N2(p) of the current point, which is defined as,

[0055] N2(p) = {q | ||q - p|| ≤ df(p)}

[0056] df(p) is the sampling spacing, df(p)=d, that is, the average point spacing d of the point cloud downsampling. If there is no selected FPFH feature point in the neighborhood, p will be added to the FPFH feature point set of the point cloud;

[0057] S2.3.2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius.

[0058] R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]

[0059] Among them, γ and λ are radius adjustment coefficients and satisfy λ>γ>0;

[0060] Obtain the support radius neighborhood of point p for subsequent FPFH feature descriptor calculation.

[0061] N3(p) = {q | ||q - p|| ≤ R(p)}

[0062] S2.3.3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, an adaptive feature weight mapping function is constructed.

[0063] Wf(p) = η·S(p)

[0064] Among them, η is the weight adjustment coefficient;

[0065] S2.3.4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set.

[0066] FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk))

[0067] Among them, FPFH(p) is the FPFH feature descriptor of point p, SPFH(p) is the local geometric feature histogram of point p, which is used to describe the local geometric features, N3(p) is the support radius neighborhood centered on point p, |N3(p)| represents the number of support radius neighborhood points, dk is the Euclidean distance from the neighborhood point pk to the center point p, and Wf(pk) is the feature weight of the neighborhood point pk.

[0068] Specifically, S3 includes the following steps:

[0069] S3.1. Calculate the feature distance of the point cloud pair, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance between the corresponding feature points of the local feature based on the SHOT feature dist_SHOT(p,q), and the distance between the corresponding feature points of the local feature based on the FPFH feature dist_FPFH(p,q).

[0070] dist_ESF = ||ESF(P) - ESF(Q)||2

[0071] dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)||2

[0072] dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)||2

[0073] Among them, dist_ESF is the overall point cloud distance based on the global ESF feature, ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q respectively, ||||2 is the Euclidean distance calculation, p and q are the corresponding feature points in the source point cloud and the target point cloud respectively, SHOT(p) and FPFH(p) are the SHOT feature descriptor and FPFH feature descriptor of the feature point p, dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q) respectively;

[0074] S3.2. Normalized feature distance of point cloud pairs,

[0075]

[0076] in, Indicates the maximum SHOT feature distance of all feature point pairs, It represents the maximum FPFH feature distance of all feature point pairs; dist_ESF_norm is the normalized ESF feature distance, dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the feature point pair (p,q), respectively.

[0077] Specifically, the S4 includes the following steps:

[0078] S4.1. Obtain the first-layer transformation matrix T_coarse and perform feature matching globally to achieve the first-layer registration of the point cloud.

[0079] S4.2. Obtain the second transformation matrix T_optimize to achieve the second level registration of the point cloud;

[0080] S4.3. Obtain the third-layer transformation matrix T_fine to implement the third-layer registration of the point cloud;

[0081] S4.4. Obtain the registered point cloud Qp, calculated by Qp=T_fine×T_optimize×T_coarse×Q.

[0082] Specifically, in S4, the process of each layer registration is as follows:

[0083] Determine the search space for registration of this layer;

[0084] Determine the weights of the ESF feature distance, SHOT feature distance, and FPFH feature distance of the point cloud pair in this layer of registration;

[0085] The fusion feature distance is calculated in a weighted manner in the search space based on the three feature distances and their respective weights;

[0086] Determine matching point pairs based on fusion feature distance and region-adaptive matching strategy;

[0087] Verify the geometric consistency of the matching point pair results;

[0088] According to the set of matching point pairs verified by geometric consistency, the transformation matrix is ​​calculated, and the point cloud to be registered is registered through the transformation matrix of each layer.

[0089] Specifically, the S4.1 includes the following specific steps:

[0090] S4.1.1. Determine the search space; in the global scope, the search space of any feature point p is,

[0091] Ω1(p) = {q ∈ Q | ||q - p|| ≤ Rmax}

[0092] Where Ω1(p) represents the search space for feature point p within the global range Rmax, and Rmax is the maximum search radius;

[0093] S4.1.2. Determine the feature weight vector,

[0094] w = (w1, w2, w3) T

[0095] Among them, w1, w2, and w3 correspond to the initial weights of ESF features, SHOT features, and FPFH features, respectively. 0.2 ≤ wi ≤ 0.5, w1 + w2 + w3 = 1, and w1 > w2, w1 > w3.

[0096] S4.1.3. Calculate the fusion feature distance D_final(p,q) of the point pair:

[0097]

[0098] S4.1.4 Determine matching point pairs based on the fused feature distance D_final(p,q) and the region-adaptive matching strategy; including:

[0099] (1) For each feature point p in the target point cloud P, calculate the fusion distance D_final(p,q) of all feature points in the point cloud to be registered Q within its search space Ω(p);

[0100] (2) Find the two points q1 and q2 with the smallest distance to the feature point p in the point cloud Q to be registered, with the corresponding distances D1 and D2, and use different matching criteria according to the area where p is located;

[0101] For the transition feature points, if D1 / D2<λ1 and the cross-validation rule is satisfied, the matching point of this feature point p in P in Q is q1; similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, the matching point of this feature point q1 in Q in P is considered to be p; if both of the above conditions are satisfied, the matching point pair (p,q1) is retained;

[0102] For the feature points of the cylindrical segment, if D1 / D2<λ2, for the feature point q1 in Q, if its nearest neighbor in P is p, the matching point of this feature point q1 in Q in P is considered to be p;

[0103] Among them, λ1<λ2<1, both are feature point application ratio thresholds;

[0104] S4.1.5. Verify the geometric consistency of the matching point pairs; calculate the geometric consistency index of any two random matching point pairs (pi, qi) and (pj, qj).

[0105] ε_ij = |||pi - pj|| / ||qi - qj|| - 1|

[0106] If ε_ij < εd, the matching point pair passes the verification, otherwise, it fails the geometric consistency verification; where (pi, pj) is the point pair in the source point cloud, (qi, qj) is the point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold;

[0107] Define Q1 = Nc / Nt, where Nc is the number of matching point pairs that pass the geometric consistency constraint, Nt is the total number of matching point pairs, and Q1 is the proportion of matching point pairs that meet the geometric constraints;

[0108] When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following updates and return to step 4.1.3 to recalculate D_final and accumulate the number of updates k.

[0109]

[0110] Among them, σpace is the weight adjustment step size, σpace=0.05, Q1>=0.6 or the cumulative update number k>5, then no weight update is performed;

[0111] S4.1.6 Based on the selected set of matching point pairs, use the RANSAC-SVD method to calculate the transformation matrix T_coarse;

[0112] The S4.2 includes the following steps:

[0113] S4.2.1. Determine the search range based on the transformation matrix T_coarse: For any feature point p, its position p = T_coarse·p after transformation by the transformation matrix T_coarse is used as the search center, and its local search space Ω2(p) is defined as,

[0114] Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}

[0115] Rmiddle=βRmax

[0116] Among them, Rmiddle is the local search radius, β is the search radius reduction coefficient;

[0117] S4.2.2. Configure the weights based on the weights defined in S4.1.2, w1 = w3 = 0.25, w2 = 0.5;

[0118] S4.2.3. Calculate the fused feature distance D_final in the registration phase based on the local search space Ω2(p), the weights, and the formula defined in S4.1.3.

[0119] S4.2.4. Calculate new matching point pairs based on the local search space Ω2(p) and the definition in S4.1.4;

[0120] S4.2.5. Calculate Q1 based on the new matching point pairs and the definition in S4.1.5. The specific judgment criteria, weight update logic, and transformation matrix calculation results are the same as in step 4.1.5.

[0121] S4.3 includes the following steps:

[0122] S4.3.1. Determine the search range based on the optimized transformation matrix T_optimize: Rmin = γRmax, γ ≈ 0.1, where Rmin is the fine search radius. For any feature point p, its position after transformation by the transformation matrix T_optimize is p = T_optimize·p as the search center. Its local search space Ω3(p) is defined as,

[0123] Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}

[0124] S4.3.2. Based on the weights defined in step 4.1.2, configure the weights w1 = 0.2, w2 = 0.3, and w3 = 0.5.

[0125] S4.3.3. Calculate D_final for the fine registration phase based on the local search space Ω3(p), the weights, and the definitions in step 4.1.3.

[0126] S4.3.4. Calculate new matching point pairs based on the local search space Ω3(p) and the definition in step 4.1.4;

[0127] S4.3.5. Define local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.4.

[0128] s_i = cos(SHOT(pi), SHOT(qi))

[0129] Among them, SHOT(pi) and SHOT(qi) are the SHOT descriptors of matching point pi and point qi respectively;

[0130] Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, and Nt is the total number of matching point pairs;

[0131] When Q2 < 0.7, it indicates that the local similarity is insufficient. Perform the following update and return to step 4.3.3 to recalculate D_final and accumulate the number of updates k.

[0132]

[0133] If Q2 >= 0.7 or the cumulative number of updates k>5, no weight update is performed. Based on the screened set of matching point pairs, the RANSAC+SVD method is used to calculate the refined optimization transformation matrix T_fine.

[0134] After adopting the above scheme, the beneficial effects of the present invention are as follows:

[0135] A weight adjustment strategy based on geometric consistency and local shape similarity achieves adaptive optimization of the feature fusion process by dynamically balancing the contributions of global and local features. When geometric consistency is poor, the weight of the ESF feature is increased to improve global structural matching; when local similarity is insufficient, the weight of the SHOT and FPFH features is increased to improve local matching accuracy. This dual verification mechanism provides a more comprehensive and reliable quality assessment criterion for feature fusion. This verification mechanism can identify matching problems: poor geometric consistency indicates a deviation in global structural matching; low local similarity indicates inaccurate local feature matching. This problem-localization capability provides clear guidance for subsequent parameter adjustments. BRIEF DESCRIPTION OF THE DRAWINGS

[0136] Figure 1 It is the main flow chart of the present invention. DETAILED DESCRIPTION

[0137] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts shall fall within the scope of protection of the present invention.

[0138] In the description of the present invention, the terms "first" and "second" are used for descriptive purposes only and should not be understood to indicate or imply relative importance or implicitly specify the number of the technical features indicated. Therefore, a feature specified as "first" or "second" may explicitly or implicitly include one or more of the specified features. In the description of the present invention, "plurality" means two or more, unless otherwise specifically defined.

[0139] In the description of the present invention, the term "for example" is used to mean "used as an example, illustration or illustration". Any embodiment of the present invention described as "for example" is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is given to enable any person skilled in the art to implement and use the present invention. In the following description, details are listed for the purpose of explanation. It should be understood that a person of ordinary skill in the art can recognize that the present invention can be implemented without using these specific details. In other examples, well-known structures and processes are not elaborated in detail to avoid obscuring the description of the present invention with unnecessary details. Therefore, the present invention is not intended to be limited to the embodiments shown, but is consistent with the widest scope consistent with the principles and features disclosed herein.

[0140] This specific implementation will provide a clearer and more complete explanation of the method of the present invention. Figure 1 As shown, the present invention provides a multi-level feature fusion adaptive point cloud registration method, comprising the following steps:

[0141] S1. Preprocessing the point cloud: Before feature extraction, the original point cloud data needs to be preprocessed to improve the reliability of subsequent processing. S1 includes the following steps:

[0142] S1.1 uses a statistical outlier removal algorithm to reduce noise on the point cloud and identifies and removes abnormal points by analyzing the neighborhood distribution characteristics of the points.

[0143] S1.2 uses a voxel grid filtering algorithm to downsample the point cloud, achieving uniformity in point cloud density while reducing data volume and improving processing efficiency.

[0144] S1.3 calculates the average point spacing d of the downsampled point cloud and uses it as the benchmark parameter for feature extraction. Specifically, this parameter will be used to determine the basic sampling spacing ds_base(p) of SHOT feature extraction and the sampling spacing df(p) of FPFH features and their support radius R(p).

[0145] S1.4 constructs an octree-based spatial index structure for neighborhood search; this structure can improve the efficiency of neighborhood search in the subsequent feature calculation process. The neighborhoods mentioned below are all based on the adjacent relationships and distance restrictions of the points determined by the octree, and are sets of adjacent points within a certain distance.

[0146] In the above steps, the point cloud pair refers to the source point cloud and the target point cloud, or the point cloud to be registered one and the point cloud to be registered two. In the subsequent steps, unless otherwise specified, this interpretation shall apply. Preprocessing lays the foundation for the subsequent feature extraction and registration process.

[0147] S2. Perform multi-scale feature extraction on the preprocessed point cloud pairs, and construct a feature pyramid based on the extracted multi-scale features; a top-down guidance mechanism is adopted between each level of the feature pyramid.

[0148] This solution adopts a multi-scale feature extraction strategy. To address the degradation problem in point cloud registration, three mature feature descriptors at different scales are selected: the macro-scale ESF feature descriptor is used to capture global structural information, the meso-scale SHOT feature descriptor describes local geometric structure, and the micro-scale FPFH feature descriptor provides fine geometric description. A feature pyramid is constructed using the ESF, SHOT, and FPFH feature descriptors to achieve feature extraction and parameter adaptive adjustment from macro to micro. The hierarchical structure of the feature pyramid is as follows from top to bottom:

[0149] 1) ESF features (top layer): capture global structural information and provide a basis for structural classification;

[0150] 2) SHOT features (middle layer): describe the local geometry and are guided by ESF;

[0151] 3) FPFH features (bottom layer): provide fine geometric description and accept SHOT guidance;

[0152] In this invention, macro, meso, and micro are relative concepts. Their purpose is to distinguish the scales of the three feature descriptors and do not represent a specific size. In the following description, these three concepts are written as the first scale, the second scale, and the third scale, respectively.

[0153] S2.1. ESF features are used to extract first-scale features, obtaining the ESF feature descriptor ESF(P). The first-scale features serve as the top layer of the pyramid, and the preprocessed point cloud is subjected to tower structure recognition and region segmentation. Cylindrical, conical, and transitional segments are obtained, and the rationality of the region segmentation is verified using ESF features. ESF features, through their global structural description capabilities, provide structural-level guidance for SHOT features. This top-down guidance mechanism ensures that mid-level features can maintain local description capabilities while gaining support from global structural information. Here, structure recognition and region segmentation specifically segment the point cloud into cylindrical, conical, and transitional segments. LiDAR operates on the principle of light signal reflection, so transition and conical segments are significantly affected by inclination angles. In actual operation, point cloud defects (i.e., expansion, contraction, irregular deformation, missing points, etc.) are more likely to occur. To address this issue, the structure is divided into three segments, which will be processed separately in subsequent steps based on the point cloud distribution.

[0154] S2.1.1 Perform cylindrical segment detection on the point cloud using RANSAC technology, extract the cylinder radius r and axial direction vector v, and use the axial direction as the z-axis direction.

[0155] S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud and perform reliability verification on the cylindrical segment using the ESF feature. If the verification passes, record the cylindrical parameters (cylinder radius r) and axial direction vector v, and mark the points in the point cloud that meet the cylindrical model features as cylindrical segments, and proceed to S2.1.4. If the verification fails, proceed to S2.1.3. The reliability verification evaluation criteria include two aspects: (1) whether the deviation between the D2 distribution of the ESF feature and the theoretical cylinder established based on the cylinder radius r extracted by RANSAC is within the range of ε1, and (2) whether the roundness deviation between the A3 component of the ESF feature and the theoretical cylinder cross section is within the range of ε2. Both of these aspects must be met for the verification to pass. In this specific implementation, ε1 is set to 0.1 and ε2 is set to 0.05, which means that the D2 distribution deviation is allowed to be within 10% and the roundness deviation is allowed to be within 5%.

[0156] S2.1.3 Crop the point cloud in the z-axis direction, removing the point cloud data of δh×H near the highest point and the lowest point, where H is the total height of the point cloud and δh is the cropping ratio; after cropping, return to step S2.1.1 and accumulate the number of cropping attempts; when the number of cropping attempts reaches the maximum number of attempts N, if it still fails to pass the verification, the point cloud quality is not sufficient to extract reliable cylindrical features, and the alignment is exited.

[0157] S2.1.4. Based on the cylindrical segments obtained by segmentation, the point cloud is further divided into cylindrical segments, transition segments and cone segments, where the transition segment is the area within the cylindrical segment boundary (-tr, tr), and the cone segment is the remaining area. In this specific implementation, the value of tr is 0.2m.

[0158] S2.2. SHOT features are used to extract second-scale features to obtain a SHOT feature descriptor, SHOT(p). The second-scale features serve as the middle layer of the pyramid to extract local geometric structure. Information from the middle layer of the feature pyramid is transferred to the bottom layer to provide guidance. The SHOT features utilize their stable description capability of local geometric structure to provide local geometric saliency guidance for the FPFH features, thereby enabling the bottom layer features to obtain more reliable local structure support while maintaining their fine description capability. S2.2 includes:

[0159] S2.2.1. Define the structure indicator function I(p),

[0160]

[0161] Among them, p represents the point in the point cloud that has been divided, β is the sampling correction coefficient of the cone segment, and μ is the sampling density coefficient of the transition segment, β∈(1,2), μ∈(1,2).

[0162] S2.2.2. Calculate the basic sampling interval:

[0163] ds_base(p) = d / ρ(p)

[0164] ρ(p) = α·I(p)

[0165] Where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the cone segment sampling correction coefficient, and ρ(p) is the sampling density function.

[0166] S2.2.3. To improve the reliability of feature extraction, a local quality evaluation is designed. The local quality score Q(p) is obtained by comprehensively evaluating the three key features.

[0167] Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3

[0168] Q_d(p)= min(1, |N(p)| / N_exp)

[0169] Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π)

[0170] Q_r(p) = exp(-σ / d)

[0171] Among them, Q_d(p) is the local point density score of point p, N(p) is the local neighborhood point set of point p, and for any point p in the point cloud, its neighborhood N(p) is defined as the set of all points in the spherical space with p as the center and radius r.

[0172] N(p) = {q | ||q - p|| ≤ r}

[0173] Where r is the neighborhood radius, ||·|| represents the Euclidean distance, and subsequent neighborhood definitions are based on this. N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood points and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of the local plane fitting based on principal component analysis.

[0174] S2.2.4. Based on the local quality score, an adaptive mapping is constructed to achieve dynamic adjustment of the sampling spacing, and the adaptive sampling spacing ds(p) is obtained:

[0175] ds(p) = ds_base(p)·(1 + γ(1-Q(p)))

[0176] Where γ is the sampling adjustment coefficient, Q(p)∈[0,1]. When Q(p)=1, the sampling interval is maintained. When Q(p) decreases, the sampling interval is gradually increased by 1+γ(1-Q(p)) to cope with unstable areas.

[0177] This adaptive sampling strategy realizes the dynamic optimization of sampling parameters through the organic combination of structural feature classification and local quality assessment, which not only ensures the reliability of feature extraction but also improves computational efficiency.

[0178] S2.2.5. Construct the SHOT local reference frame; the construction of the local reference frame also adopts different strategies according to the area type. The local reference frame refers to the local coordinate system established for each point when calculating the SHOT feature. This local coordinate system is used to transform the local area of ​​the point cloud into a standardized reference space, making the feature description rotationally invariant, ensuring that the same geometric structure obtains consistent feature descriptions under different perspectives, and improving the accuracy of feature matching.

[0179] Specifically, for the cylindrical segment and the transition segment on the cylindrical structure, the identified cylindrical axis is used as the main direction (z-axis) of the reference frame, and the secondary directions (x-axis and y-axis) are determined by combining the normal vectors of the local neighborhood point set of each point on the cylindrical segment point cloud, thereby establishing a stable local coordinate system.

[0180] For the conical segment and transition segment on the conical structure, the main direction of the local surface is calculated, and the principal component analysis method is used to determine the complete reference frame (x, y, and z directions); specifically, the local neighborhood point set of each point on the conical segment point cloud is first obtained (this article specifically uses N(p) and local neighborhood point set, and the two concepts are completely equivalent. When point p is not emphasized, the local neighborhood point set is used, and if point p is emphasized, N(p) is used.), and the point set is subjected to principal component analysis. The eigenvector corresponding to the maximum eigenvalue is used as the main direction of the reference frame (z-axis), the eigenvector corresponding to the second largest eigenvalue is used as the secondary direction (x-axis or y-axis), and the eigenvector corresponding to the minimum eigenvalue is used as the third direction (y-axis or x-axis), thereby establishing a complete reference frame; this method can effectively capture the main change directions of the local geometric structure.

[0181] S2.2.6. Perform SHOT feature point extraction on the point cloud. For each point p in the point cloud, sort the points in the point cloud based on its local quality score Q(p). Traverse the sorted point list and find the neighborhood N1(p) of the current point p, where N1(p) = {q | ||q- p|| ≤ ds(p)}

[0182] Where ds(p) is the adaptive sampling spacing; if there is no added SHOT feature point in the neighborhood N1(p) of the current point p, then p will be added to the SHOT feature point set.

[0183] S2.2.7. Calculate the SHOT feature descriptor SHOT(p) and calculate the feature response value S(p) for each point p in the point cloud SHOT feature point set, including the following steps:

[0184] (1) Based on the local reference frame established in S2.2.5, the local spherical support area of ​​point p is divided into 32 spatial grids. The normal vector distribution histogram of 11 bins is calculated in each grid to obtain a 352-dimensional SHOT feature vector. This SHOT feature vector is the SHOT feature descriptor SHOT(p).

[0185] (2) Reorganize the eigenvector into a 32×11 matrix M with 11 elements per row to calculate the characteristic response value S(p).

[0186] S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|)

[0187] Where M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i,j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation. The response value S(p) reflects the significance of the local geometric features at point p: when S(p) is large, it indicates that the point is located in an area with obvious geometric features, such as an edge or corner; when S(p) is small, it indicates that the area where the point is located has relatively gentle geometric changes, such as a flat surface. This method of calculating feature response values ​​provides an important basis for the subsequent adaptive adjustment of feature extraction parameters.

[0188] S2.3. Use FPFH features to extract third-scale features and obtain the FPFH feature descriptor FPFH(p). The third-scale features are used as the lower layer of the pyramid to achieve fine geometric structure extraction.

[0189] The above 2.3 specifically includes:

[0190] S2.3.1. Take dense features, and use a larger spacing value νd in places where the features are not obvious to improve computational efficiency;

[0191] Sort the points in the SHOT feature point set by the feature response value S(p), traverse the sorted point list, and check the neighborhood N2(p) of the current point, which is defined as,

[0192] N2(p) = {q | ||q - p|| ≤ df(p)}

[0193] df(p) is the sampling spacing, df(p)=d, that is, the average point spacing d of the point cloud downsampling. If there is no selected FPFH feature point in the neighborhood, p will be added to the FPFH feature point set of the point cloud.

[0194] S2.3.2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius.

[0195] R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]

[0196] Where γ and λ are radius adjustment coefficients and satisfy λ>γ>0; this allows the use of a smaller support radius γd to preserve details in regions with significant geometric features (S(p)≥θ), and a larger support radius λd to improve stability in regions with less significant features (S(p)<θ);

[0197] Obtain the support radius neighborhood of point p for subsequent FPFH feature descriptor calculation.

[0198] N3(p) = {q | ||q - p|| ≤ R(p)}

[0199] S2.3.3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, an adaptive feature weight mapping function is constructed.

[0200] Wf(p) = η·S(p)

[0201] Where η is the weight adjustment coefficient; the weight is proportional to the local geometric significance, ensuring that the geometrically significant region has a greater contribution in the feature calculation;

[0202] S2.3.4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set.

[0203] FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk))

[0204] Where FPFH(p) is the FPFH feature descriptor for point p, SPFH(p) is the local geometric feature histogram of point p, used to describe local geometric features, N3(p) is the support radius neighborhood centered on point p, |N3(p)| represents the number of points in the support radius neighborhood, dk is the Euclidean distance from neighboring point pk to the center point p, and Wf(pk) is the feature weight of neighboring point pk. The specific calculation process is based on existing mature theory and will not be repeated here.

[0205] Finally, for each point cloud, its feature point set is its FPFH feature point set. The feature points in subsequent S3 and S4 are all feature points of this FPFH feature point set. These points are only used to select the registration points. The specific registration features also need to use the ESF features and SHOT features calculated previously.

[0206] S3. Perform preliminary feature fusion on the feature points of the point cloud pair based on the extracted multi-scale features; specifically, the following steps are included:

[0207] S3.1. Calculate the feature distance of the point cloud pair, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance between the corresponding feature points of the local feature based on the SHOT feature dist_SHOT(p,q), and the distance between the corresponding feature points of the local feature based on the FPFH feature dist_FPFH(p,q).

[0208] dist_ESF = ||ESF(P) - ESF(Q)||2

[0209] dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)||2

[0210] dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)||2

[0211] Among them, dist_ESF is the overall point cloud distance based on the global ESF feature, which is used to measure the global structural similarity. ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q, respectively. ||||2 is the Euclidean distance calculation, p and q are the corresponding feature points in the source point cloud and the target point cloud, respectively. SHOT(p) and FPFH(p) are the SHOT feature descriptor and FPFH feature descriptor of the feature point p, and dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q), respectively.

[0212] S3.2. Normalized feature distance of point cloud pairs,

[0213]

[0214] in, Indicates the maximum SHOT feature distance of all feature point pairs, Represents the maximum FPFH feature distance of all feature point pairs; dist_ESF_norm is the normalized ESF feature distance, dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the feature point pair (p,q), respectively. By dividing by their respective maximum values, all feature distances are mapped to a unified [0,1] interval, providing a comparable metric for subsequent feature fusion.

[0215] S4. Perform multi-level registration on the point cloud pairs. The multi-level registration proposed in this invention realizes a progressive registration process from coarse to fine through the dynamic adjustment of feature dominant weights and the information transmission mechanism between levels. The multi-level registration includes the first level registration, the second level registration, and the third level registration. After each level registration is completed, the transformation matrix of that level is obtained, and the registration point cloud to be registered is performed using the transformation matrix of each level. These three levels of registration can be regarded as coarse registration, medium registration, and fine registration. The specific process of each level registration is as follows:

[0216] (1) Determine the search space for registration of this layer.

[0217] (2) Determine the weights of the ESF feature distance, SHOT feature distance, and FPFH feature distance of the point cloud pair in this layer of registration.

[0218] (3) The fusion feature distance is calculated in a weighted manner in the search space based on the three feature distances and their respective weights.

[0219] (4) Determine the matching point pairs based on the fusion feature distance and region-adaptive matching strategy.

[0220] (5) Verify the geometric consistency of the matching point pair results.

[0221] (6) Based on the set of matching point pairs verified by geometric consistency, the transformation matrix is ​​calculated, and the point cloud to be registered is registered through the transformation matrix of each layer.

[0222] The steps include:

[0223] S4.1. The first step in multi-level registration is coarse registration. The main goal of this stage is to obtain the initial transformation matrix, namely the first-level transformation matrix T_coarse, and perform feature matching on a global scale to achieve the first-level registration of the point cloud. In terms of feature configuration, the ESF global feature is dominant. The registration process includes the following steps:

[0224] S4.1.1. Determine the search space; in the global scope, the search space of any feature point p is,

[0225] Ω1(p) = {q ∈ Q | ||q - p|| ≤ Rmax}

[0226] Where Ω1(p) represents the search space for feature point p within the global range Rmax, that is, the set of points in the point cloud Q to be registered that may match p. Rmax is the maximum search radius, which is the threshold used to limit the search range.

[0227] S4.1.2. Determine the feature weight vector,

[0228] w = (w1, w2, w3) T

[0229] Among them, w1, w2, and w3 correspond to the initial weights of ESF features, SHOT features, and FPFH features, respectively. 0.2 ≤ wi ≤ 0.5 ensures the basic contribution of each feature. w1 + w2 + w3 = 1 ensures weight normalization, and w1 > w2, w1 > w3.

[0230] In the coarse registration stage of this specific implementation, the weight configuration satisfies: w1 = 0.5, w2 = w3 = 0.25, and the reliability of the overall pose estimation is ensured by increasing the weight of the global features.

[0231] S4.1.3. Calculate the fusion feature distance D_final(p,q) of the point pair,

[0232]

[0233] Through this weighted fusion method, an adaptive combination of multi-feature distances is achieved, which not only ensures that the contribution of various features will not deviate excessively, but also maintains a certain adjustment space.

[0234] S4.1.4 Determine matching point pairs based on the fusion feature distance D_final(p,q) and the region-adaptive matching strategy; including:

[0235] (1) For each feature point p in the target point cloud P, calculate the fusion distance D_final(p,q) of all feature points in the point cloud Q to be registered within its search space Ω(p).

[0236] (2) Find the two points q1 and q2 with the smallest distance to the feature point p in the point cloud Q to be registered, with corresponding distances D1 and D2, and use different matching criteria according to the area where p is located.

[0237] For the feature points of the transition section, if D1 / D2<λ1 and the cross-validation rule is required to be met (the nearest neighbor of q1 in P is p), then the matching point of this feature point p in P in Q is q1; similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p; if both of the above are met, the matching point pair (p,q1) is retained.

[0238] For cylindrical segment feature points, if D1 / D2<λ2, the matching point of feature point p in P is considered to be q1 in Q; for feature point q1 in Q, if its nearest neighbor in P is p, the matching point of feature point q1 in Q is considered to be p in P. If both of the above conditions are met, the matching point pair (p,q1) is retained.

[0239] Among them, λ1<λ2<1, both are ratio thresholds for feature point application. The smaller the values ​​of λ1 and λ2, the stricter the ratio test, that is, the higher the degree of matching of the point pairs that pass the test. However, too small a value will result in too few matching points, and it needs to be adjusted based on actual experience.

[0240] S4.1.5. Verify the geometric consistency of the matching point pairs and select them; use them to examine the distance ratio relationship between the matching point pairs; calculate the geometric consistency index of any two random matching point pairs (pi, qi) and (pj, qj).

[0241] ε_ij = |||pi - pj|| / ||qi - qj|| - 1|

[0242] If ε_ij <εd, the matching point pair passes the verification, otherwise, it fails the geometric consistency verification; where (pi, pj) is the point pair in the source point cloud, (qi, qj) is the point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold.

[0243] Define Q1 = Nc / Nt, where Nc is the number of matching point pairs that pass the geometric consistency constraint, Nt is the total number of matching point pairs, and Q1 is the proportion of matching point pairs that meet the geometric constraints;

[0244] When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following updates and return to step 4.1.3 to recalculate D_final and accumulate the number of updates k.

[0245]

[0246] Among them, σpace is the weight adjustment step size, σpace=0.05, Q1>=0.6 or the cumulative update number k>5, then no weight update is performed.

[0247] S4.1.6 Based on the screened set of matching point pairs, the RANSAC-SVD method is used to calculate the transformation matrix T_coarse (the specific implementation method is the classic transformation matrix calculation method, which will not be repeated here).

[0248] The second step of multi-level registration is mid-registration. The main goal of this stage is to obtain the optimized transformation matrix, namely the second transformation matrix T_optimize, to achieve the second-level registration of the point cloud and improve the alignment accuracy of the structural features. In terms of feature configuration, SHOT features are the main focus. The registration process includes the following steps:

[0249] S4.2.1. Determine the search range based on the transformation matrix T_coarse: For any feature point p, its position p = T_coarse·p after transformation by the transformation matrix T_coarse is used as the search center, and its local search space Ω2(p) is defined as,

[0250] Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}

[0251] Rmiddle=βRmax

[0252] Wherein, Rmiddle is the local search radius, β is the search radius reduction coefficient, which is 0.3 in this specific implementation, indicating that the search radius is reduced to 30% of the global search radius by the coefficient β.

[0253] S4.2.2. Based on the weights defined in step 4.1.2, configure the weights as w1 = w3 = 0.25 and w2 = 0.5 to highlight the dominant role of structural features.

[0254] S4.2.3. Based on the above local search space Ω2(p), weights and the definition in step 4.1.3, calculate the fusion feature distance D_final in the registration stage.

[0255] S4.2.4. Calculate new matching point pairs based on the above local search space Ω2(p) and the definition in step 4.1.4. The new matching point pairs are based on a smaller search range and weights that are more concerned with structural features. Compared with the matching point pairs in the coarse alignment stage of step 4.1.4, more accurate feature point matching results are ensured at the structural feature level.

[0256] S4.2.5. Calculate Q1 based on the new matching point pairs and the definition in step 4.1.5. That is, geometric consistency must also be ensured in the mid-registration stage. The specific judgment criteria, weight update logic, and transformation matrix calculation results are the same as those in step 4.1.5.

[0257] S4.3. The third step of multi-level registration is fine registration. The main goal of this stage is to obtain a finely optimized transformation matrix, namely the third-level transformation matrix T_fine, to achieve the third-level registration of the point cloud, namely fine alignment. In terms of feature configuration, the FPFH feature is the main feature, and the registration process includes the following steps:

[0258] S4.3.1. Determine the search range based on the optimization transformation matrix T_optimize:

[0259] For any feature point p, its position p=T_optimize·p after transformation by the transformation matrix T_optimize is used as the search center, and its local search space Ω3(p) is defined as,

[0260] Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}

[0261] Rmin = γRmax

[0262] Where Rmin is the fine search radius, γ ≈ 0.1, which is taken as 0.1 in this implementation. The coefficient γ is used to further narrow the search range to 10% of the global range.

[0263] In this fine registration stage, the weights defined in step 4.1.2 are configured as w1 = 0.2, w2 = 0.3, and w3 = 0.5, respectively; this emphasizes the role of local fine features.

[0264] S4.3.3. Based on the above local search space Ω3(p), weights and the definition in step 4.1.3, calculate D_final of the fine alignment stage.

[0265] S4.3.4. Calculate new matching point pairs based on the above local search space Ω3(p) and the definition in step 4.1.4. The new matching point pairs are based on a smaller search range and are more concerned with weights of local fine features. Compared with the matching point pairs in the mid-registration stage in step 4.2.4, they ensure more accurate feature point matching results at the local fine feature level.

[0266] S4.3.5. Define local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.4.

[0267] s_i = cos(SHOT(pi), SHOT(qi))

[0268] Among them, SHOT(pi) and SHOT(qi) are the SHOT descriptors of matching point pi and point qi respectively;

[0269] Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, which is obtained by calculating their cosine similarity and averaging all matching point pairs, where Nt is the total number of matching point pairs;

[0270] When Q2 < 0.7, it indicates that the local similarity is insufficient. Perform the following update and return to step 4.3.3 to recalculate D_final and accumulate the number of updates k.

[0271]

[0272] If Q2 >= 0.7 or the cumulative number of updates k>5, no weight update is performed. Based on the screened set of matching point pairs, the RANSAC+SVD method is used to calculate the refined optimization transformation matrix T_fine.

[0273] S4.4. Obtain the registered point cloud Qp, calculated by Qp=T_fine×T_optimize×T_coarse×Q.

[0274] Although the preferred embodiments of the present invention have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present invention.

[0275] Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if such changes and modifications fall within the scope of the claims and their equivalents, the present invention is intended to include such changes and modifications.

Claims

1. A multi-level feature fusion adaptive point cloud registration method, characterized by: include: S1. Preprocess the point cloud pair; S2. Perform multi-scale feature extraction on the pre-processed point cloud pairs, and construct a feature pyramid based on the extracted multi-scale features. A top-down guidance mechanism is adopted between the levels of the feature pyramid. S3. Perform preliminary feature fusion on the feature points of the point cloud pair based on the extracted multi-scale features; S4. Perform multi-level registration on the point cloud pairs; the multi-level registration includes first-level registration, second-level registration, and third-level registration. After each level of registration is completed, the transformation matrix of that level is obtained, and the registration point cloud to be registered is registered using the transformation matrix of each level. The S2 comprises the following steps: S2.

1. Use ESF features to extract first-scale features and obtain the ESF feature descriptor ESF(P). Using the first-scale features as the top layer of the pyramid, perform tower structure recognition and region segmentation on the preprocessed point cloud. Cylindrical segments, conical segments, and transition segments are obtained, and the rationality of the region segmentation is verified using ESF features. S2.

2. Use the SHOT feature to extract second-scale features, obtaining the SHOT feature descriptor SHOT(p). The second-scale feature serves as the middle layer of the pyramid, used to extract local geometric structure and transmit information from the middle layer of the feature pyramid to the bottom layer, providing guidance for the bottom layer. S2.

3. Use the FPFH feature to extract third-scale features and obtain the FPFH feature descriptor FPFH(p). The third-scale feature serves as the lower layer of the pyramid to extract fine geometric structures. For each point cloud, its feature point set is its FPFH feature point set, and the feature points in subsequent S3 and S4 are all feature points of this FPFH feature point set; The S3 includes the following steps: S3.

1. Calculate the feature distance of the point cloud pair, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance between the corresponding feature points of the local feature based on the SHOT feature dist_SHOT(p,q), and the distance between the corresponding feature points of the local feature based on the FPFH feature dist_FPFH(p,q). dist_ESF = ||ESF(P) - ESF(Q)|| 2； dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)|| 2； dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)|| 2； Among them, dist_ESF is the overall point cloud distance based on the global ESF feature, ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q respectively, ||||2 is the Euclidean distance calculation, p and q are the corresponding feature points in the source point cloud and the target point cloud respectively, SHOT(p) and FPFH(p) are the SHOT feature descriptor and FPFH feature descriptor of the feature point p, dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q) respectively; S3.

2. Normalized feature distance of point cloud pairs, ； in Indicates the maximum SHOT feature distance of all feature point pairs, Indicates the maximum FPFH feature distance of all feature point pairs; dist_ESF_norm is the normalized ESF feature distance, dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the feature point pair (p,q) respectively; The S4 comprises the following steps: S4.

1. Obtain the first-layer transformation matrix T_coarse and perform feature matching globally to achieve the first-layer registration of the point cloud. S4.

2. Obtain the second transformation matrix T_optimize to achieve the second level registration of the point cloud; S4.

3. Obtain the third-layer transformation matrix T_fine to implement the third-layer registration of the point cloud; S4.

4. Obtain the registered point cloud Qp, calculated by Qp=T_fine×T_optimize×T_coarse×Q.

2. The multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that: The S1 comprises the following steps: S1.1 performs noise reduction on the point cloud and identifies and removes abnormal points by analyzing the neighborhood distribution characteristics of the points; S1.2 uses voxel grid filtering algorithm to downsample the point cloud; S1.3 calculates the average point spacing d of the downsampled point cloud and uses it as a benchmark parameter for feature extraction; S1.4 builds an octree-based spatial index structure for neighborhood search.

3. The multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that: The S2.1 includes the following steps: S2.1.1 Perform cylindrical segment detection on the point cloud, extract the cylinder radius r and axial direction vector v, and use the axial direction as the z-axis direction; S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud, and perform reliability verification on the cylindrical segment through the ESF feature. If the verification passes, record the cylindrical parameters, the cylindrical radius r, and the axial direction vector v, and mark the points in the point cloud that meet the cylindrical model features as cylindrical segments, and enter S2.1.4; if the verification fails, enter S2.1.3; the reliability verification evaluation criteria include two aspects: (1) whether the deviation between the D2 distribution of the ESF feature and the theoretical cylinder established based on the cylindrical radius r extracted by RANSAC is within the range of ε1, and (2) whether the roundness deviation between the A3 component of the ESF feature and the theoretical cylindrical cross section is within the range of ε2. Both aspects must be met at the same time for the verification to pass; S2.1.3 Crop the point cloud along the z-axis, removing δh × H points near the highest and lowest points, where H is the total height of the point cloud and δh is the cropping ratio. After cropping, return to step S2.1.1 and accumulate the number of cropping attempts. If verification fails after reaching the maximum number of cropping attempts, the point cloud quality is insufficient for reliable cylindrical feature extraction, and registration is aborted. S2.1.

4. Based on the cylindrical segments obtained by segmentation, the point cloud is further divided into cylindrical segments, transition segments and cone segments, where the transition segment is the area within the cylindrical segment boundary (-tr, tr) and the cone segment is the remaining area.

4. The multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that: The S2.2 includes the following steps: S2.2.

1. Define the structure indicator function I(p), ； Among them, p represents the point in the point cloud that has been divided, β is the sampling correction coefficient of the cone segment, and μ is the sampling density coefficient of the transition segment, β∈(1,2), μ∈(1,2); S2.2.

2. Calculate the basic sampling interval, ds_base(p) = d / ρ(p); ρ(p) = α·I(p); Where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the cone segment sampling correction coefficient, and ρ(p) is the sampling density function; S2.2.

3. Obtain a local quality score Q(p), Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3; Q_d(p)= min(1, |N(p)| / N_exp); Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π); Q_r(p) = exp(-σ / d); Among them, Q_d(p) is the local point density score of point p, N(p) is the local neighborhood point set of point p, and for any point p in the point cloud, its neighborhood N(p) is defined as the set of all points in the spherical space with p as the center and radius r. N(p) = {q | ||q - p|| ≤ r}; Where r is the neighborhood radius, ||·|| represents the Euclidean distance, N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood points and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of the local plane fitting based on principal component analysis. S2.2.

4. Adaptive sampling spacing ds(p) is obtained based on the local quality score. ds(p) = ds_base(p)·(1 + γ(1-Q(p))); Where γ is the sampling adjustment coefficient, Q(p)∈[0,1]. When Q(p)=1, the sampling interval is maintained. When Q(p) decreases, the sampling interval is gradually increased by 1+γ(1-Q(p)). S2.2.

5. Construct the SHOT local reference frame; For the cylindrical segment and the transition segment on the cylindrical structure, the identified cylindrical axis is used as the main direction of the reference frame, and the secondary direction is determined by combining the normal vector of the local neighborhood point set of each point on the cylindrical segment point cloud; For the portion of the tapered segment and transition segment on the conical structure, the principal direction of the local surface is calculated and the principal component analysis method is used to determine the complete reference frame. Specifically, the local neighborhood point set of each point on the tapered segment point cloud is first obtained, and principal component analysis is performed on this point set. The eigenvector corresponding to the largest eigenvalue is used as the principal direction of the reference frame, the eigenvector corresponding to the second largest eigenvalue is used as the secondary direction, and the eigenvector corresponding to the smallest eigenvalue is used as the third direction. S2.2.

6. Perform SHOT feature point extraction on the point cloud. For each point p in the point cloud, sort the points in the point cloud based on its local quality score Q(p). Traverse the sorted point list and find the neighborhood N1(p) of the current point p, where N1(p) = {q | ||q - p|| ≤ ds(p)} Where ds(p) is the adaptive sampling spacing; if there is no added SHOT feature point in the neighborhood N1(p) of the current point p, then p is added to the SHOT feature point set; S2.2.

7. Calculate the SHOT feature descriptor SHOT(p) and calculate the feature response value S(p) for each point p in the SHOT feature point set of the point cloud, including the following steps: (1) Based on the local reference frame established in S2.2.5, the local spherical support region of point p is divided into 32 spatial grids. The normal vector distribution histogram of 11 bins is calculated in each grid to obtain a 352-dimensional SHOT feature vector. This SHOT feature vector is the SHOT feature descriptor SHOT(p); (2) Reorganize the eigenvector into a 32×11 matrix M with 11 elements per row to calculate the characteristic response value S(p). S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|); Where M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i,j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation.

5. The multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that: The 2.3 step includes the following steps: S2.3.

1. Sort the points in the SHOT feature point set by their feature response values ​​S(p), traverse the sorted point list, and check the neighborhood N2(p) of the current point, which is defined as, N2(p) = {q | ||q - p|| ≤ df(p)}; df(p) is the sampling spacing, df(p)=d, that is, the average point spacing d of the point cloud downsampling. If there is no selected FPFH feature point in the neighborhood, p will be added to the FPFH feature point set of the point cloud; S2.3.

2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius. R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]; Among them, γ and λ are radius adjustment coefficients and satisfy λ>γ>0; Obtain the support radius neighborhood of point p for subsequent FPFH feature descriptor calculation. N3(p) = {q | ||q - p|| ≤ R(p)}; S2.3.

3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, an adaptive feature weight mapping function is constructed. Wf(p) = η·S(p); Among them, η is the weight adjustment coefficient; S2.3.

4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set. FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk)); Among them, FPFH(p) is the FPFH feature descriptor of point p, SPFH(p) is the local geometric feature histogram of point p, which is used to describe the local geometric features, N3(p) is the support radius neighborhood centered on point p, |N3(p)| represents the number of support radius neighborhood points, dk is the Euclidean distance from the neighborhood point pk to the center point p, and Wf(pk) is the feature weight of the neighborhood point pk.

6. The multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that: In S4, the process of each layer registration is as follows: Determine the search space for registration of this layer; Determine the weights of the ESF feature distance, SHOT feature distance, and FPFH feature distance of the point cloud pair in this layer of registration; The fusion feature distance is calculated in a weighted manner in the search space based on the three feature distances and their respective weights; Determine matching point pairs based on fusion feature distance and region-adaptive matching strategy; Verify the geometric consistency of the matching point pair results; According to the set of matching point pairs verified by geometric consistency, the transformation matrix is ​​calculated, and the point cloud to be registered is registered through the transformation matrix of each layer.

7. The multi-level feature fusion adaptive point cloud registration method according to claim 6, characterized in that: S4.1 includes the following specific steps: S4.1.

1. Determine the search space; in the global scope, the search space of any feature point p is, Ω1(p) = {q ∈ Q | ||q - p|| ≤ Rmax}; Where Ω1(p) represents the search space for feature point p within the global range Rmax, and Rmax is the maximum search radius; S4.1.

2. Determine the feature weight vector, <h2 style=";text-align:left;direction:ltr">w = (w1, w2, w3)<h2 style=";text-align:left;direction:ltr"> T <h2 style=";text-align:left;direction:ltr"> ; Among them, w1, w2, and w3 correspond to the initial weights of ESF features, SHOT features, and FPFH features, respectively, 0.2 ≤ wi ≤ 0.5, w1 + w2 + w3 = 1, and w1 > w2, w1 > w3; S4.1.

3. Calculate the fusion feature distance D_final(p,q) of the point pair: D_final(p,q) = w1×dist_ESF_norm + w2×dist_SHOT_norm(p,q) + w3×dist_FPFH_norm(p,q); S4.1.4 Determine matching point pairs based on the fused feature distance D_final(p,q) and the region-adaptive matching strategy; including: (1) For each feature point p in the target point cloud P, calculate the fusion distance D_final(p,q) of all feature points in the point cloud to be registered Q within its search space Ω(p); (2) Find the two points q1 and q2 with the smallest distance to the feature point p in the point cloud Q to be registered, with the corresponding distances D1 and D2, and use different matching criteria according to the area where p is located; For the transition feature points, if D1 / D2<λ1 and the cross-validation rule is satisfied, the matching point of this feature point p in P in Q is q1; similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, the matching point of this feature point q1 in Q in P is considered to be p; if both of the above conditions are satisfied, the matching point pair (p,q1) is retained; For the feature points of the cylindrical segment, if D1 / D2<λ2, for the feature point q1 in Q, if its nearest neighbor in P is p, the matching point of this feature point q1 in Q in P is considered to be p; Among them, λ1<λ2<1, both are feature point application ratio thresholds; S4.1.

5. Verify the geometric consistency of the matching point pairs; calculate the geometric consistency index of any two random matching point pairs (pi, qi) and (pj, qj). ε_ij = |||pi - pj|| / ||qi - qj|| - 1|; If ε_ij < εd, the matching point pair passes the verification, otherwise, it fails the geometric consistency verification; where (pi, pj) is the point pair in the source point cloud, (qi, qj) is the point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold; Define Q1 = Nc / Nt, where Nc is the number of matching point pairs that pass the geometric consistency constraint, Nt is the total number of matching point pairs, and Q1 is the proportion of matching point pairs that meet the geometric constraints; When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following update and return to step 4.1.

3. Recalculate D_final and accumulate the number of updates k. ； Among them, σpace is the weight adjustment step size, σpace=0.05, Q1>=0.6 or the cumulative update number k>5, then no weight update is performed; S4.1.6 Based on the selected set of matching point pairs, use the RANSAC-SVD method to calculate the transformation matrix T_coarse; The S4.2 includes the following steps: S4.2.

1. Determine the search range based on the transformation matrix T_coarse: For any feature point p, its position p = T_coarse·p after transformation by the transformation matrix T_coarse is used as the search center, and its local search space Ω2(p) is defined as, Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}; Rmiddle=βRmax; Among them, Rmiddle is the local search radius, β is the search radius reduction coefficient; S4.2.

2. Configure the weights based on the weights defined in S4.1.2, w1 = w3 = 0.25, w2 = 0.5; S4.2.

3. Calculate the fused feature distance D_final in the registration phase based on the local search space Ω2(p), the weights, and the formula defined in S4.1.

3. S4.2.

4. Calculate new matching point pairs based on the local search space Ω2(p) and the definition in S4.1.4; S4.2.

5. Calculate Q1 based on the new matching point pairs and the definition in S4.1.

5. The specific judgment criteria, weight update logic, and transformation matrix calculation results are the same as in step 4.1.

5. S4.3 includes the following steps: S4.3.

1. Determine the search range based on the optimized transformation matrix T_optimize: Rmin = γRmax, γ ≈ 0.1, where Rmin is the fine search radius. For any feature point p, its position after transformation by the transformation matrix T_optimize is p = T_optimize·p as the search center. Its local search space Ω3(p) is defined as, Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}; S4.3.

2. Based on the weights defined in step 4.1.2, configure the weights: w1 = 0.2, w2 = 0.3, w3 = 0.5; S4.3.

3. Calculate D_final for the fine registration phase based on the local search space Ω3(p), the weights, and the definitions in step 4.1.

3. S4.3.

4. Calculate new matching point pairs based on the local search space Ω3(p) and the definition in step 4.1.4; S4.3.

5. Define local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.

4. s_i = cos(SHOT(pi), SHOT(qi)) Among them, SHOT(pi) and SHOT(qi) are the SHOT descriptors of matching point pi and point qi respectively; Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, and Nt is the total number of matching point pairs; When Q2 < 0.7, it indicates that the local similarity is insufficient. Perform the following update and return to step 4.3.3 to recalculate D_final and accumulate the number of updates k. ； If Q2 >= 0.7 or the cumulative number of updates k>5, no weight update is performed. Based on the screened set of matching point pairs, the RANSAC+SVD method is used to calculate the refined optimization transformation matrix T_fine.

Images