Local point cloud registration method based on geometric-spectral-pixel feature descriptor

Abstract

The application relates to the technical field of point cloud data processing, and discloses a local point cloud registration method based on a geometry-spectrum-pixel feature descriptor, which uses the following steps: firstly, taking the contour feature points on the surface of a workpiece as key points for registration, and obtaining contour key points with moderate quantity and stable position through thinning and downsampling; secondly, constructing a feature descriptor from two aspects: the relative distribution of the key point set and the local distribution of the original point cloud in the neighborhood of the key points, and constructing a stable local reference coordinate system by using improved PCA; finally, combining the above features to form a GSP feature descriptor, solving a spatial transformation matrix, and verifying the registration. The application can reduce the redundancy and misregistration during sampling point matching, improve the distinguishing ability for similar local structures, and improve the registration accuracy and reliability under the conditions of different initial poses and incomplete local point cloud data.

Description

Technical Field

[0001] This invention relates to the field of point cloud data processing technology, and in particular to a method for registering local scanned point clouds and model point clouds of free-form surface workpieces or geometrically complex workpieces, specifically a local point cloud registration method based on geometric-spectrum-pixel feature descriptors. Background Technology

[0002] With the development of laser scanning technology, 3D reconstruction technology, and point cloud processing methods, 3D point cloud data has been widely used in industrial measurement, complex part inspection, reverse modeling, robot localization, and digital assembly. For free-form surface workpieces or workpieces with complex geometry, the scanning perspective, occlusion, workpiece placement, and equipment workspace limitations often make it difficult to obtain complete surface data in a single scan. To form a complete and uniform 3D model, it is often necessary to accurately register local scan point clouds obtained from different poses with ideal CAD models, complete model point clouds, or other local scan point clouds.

[0003] The goal of point cloud registration is to solve the spatial rigid body transformation between the source and target point clouds, placing them in a unified coordinate system. Existing point cloud registration methods can be broadly categorized into iterative optimization-based methods, global sampling-based methods, deep learning-based methods, and feature descriptor-based methods. Iterative nearest-point algorithms are highly efficient when the initial pose is good and the overlap region is sufficient, but they are prone to getting trapped in local optima when the initial pose deviation is large or the local overlap rate is low. Methods such as 4PCS and SAC-IA can handle large pose differences to a certain extent, but they may still produce incorrect correspondences when dealing with regular surfaces, similar features, or locally missing data in complex workpieces. Deep learning methods have strong data representation capabilities, but generally require training data, computational resources, and model deployment costs.

[0004] Registration methods based on feature descriptors typically extract key points from the point cloud, construct local descriptors around these key points, and then find matching point pairs through descriptor similarity to solve for spatial transformation. 3D Shape Context, FPFH, and SHOT descriptors are commonly used methods for describing local point cloud features. However, traditional descriptors have certain shortcomings when applied to contour key points: Firstly, complex workpiece surfaces contain numerous surface intersections, holes, steps, boundaries, and transition regions, often failing to satisfy the assumption of locally smooth surfaces at contour key points. Traditional methods based on fitting planes or normal histograms are easily affected by normal instability. Secondly, locally scanned point clouds only cover a portion of the model; local surface similarity, inconsistent key point distribution, and scan gaps can lead to mismatches in the descriptors.

[0005] Therefore, how to leverage the saliency and discriminative characteristics of key points on the workpiece contour, while avoiding misregistration caused by instability in the contour normal and similarity in local point distribution, is a key issue in improving the accuracy and stability of local point cloud registration for complex workpieces. To address this problem, it is necessary to propose a novel feature description and matching filtering method for key points on the contour. This method should comprehensively describe the relative geometric distribution between key points and the spatial distribution of the original point cloud in the neighborhood of each key point, achieving high-precision point cloud registration under conditions of non-overlapping prior art, poor initial pose, and incomplete local data. Summary of the Invention

[0006] To overcome the problems of unstable normals at contour key points, easy mismatch of local similar features, and insufficient registration accuracy when the initial pose of the local scanned point cloud and the model point cloud differs greatly in the existing technology, the purpose of this invention is to provide a local point cloud registration method based on geometry-spectrum-pixel feature descriptors. This method uses contour feature points as registration key points, which can utilize the significant structural information of the workpiece, reduce redundancy and mismatch during sampling point matching, improve the ability to distinguish similar local structures, reduce the instability of traditional normal estimation at contour points, make spherical mesh statistics have better directional consistency, maintain the rotation invariance of the descriptor to a certain extent, and improve the registration accuracy and reliability under conditions of poor initial pose and incomplete local point cloud data.

[0007] This invention employs the following technical solution: a local point cloud registration method based on geometry-spectrum-pixel feature descriptors. First, the contour feature points of the workpiece surface are used as key points for registration, and a moderate number of stable contour key points are obtained through thinning and downsampling. Second, feature descriptions are constructed from two levels: the relative distribution of the key point set and the local distribution of the original point cloud in the neighborhood of the key points. For the relative distribution of the key point set, three statistical features—angle, distance, and shape—are used for description. For the local distribution of the original point cloud in the neighborhood of the key points, a stable local reference coordinate system is constructed using improved PCA. A spherical grid is then divided within this coordinate system, and Fourier spectrum features and pixel image features are further generated. Finally, the above features are combined to form a GSP feature descriptor (geometry-spectrum-pixel feature descriptor). Candidate matching point pairs are determined through multi-index similarity screening, curvature filtering, and an optimal matching mechanism. The spatial transformation matrix is ​​then solved, and registration verification is performed.

[0008] Specifically, the local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention includes the following steps: Step S1: Fix the workpiece position, use a 3D scanner to scan and measure, obtain the surface information of the workpiece and store it as 3D point cloud data; The aforementioned 3D scanner is a non-contact optical scanner that emits structured light for measurement. During operation, the scanner laser emits multiple laser lines to acquire 3D local point cloud data of the surface of the workpiece from different perspectives.

[0009] Step S2: Using point cloud curvature as an evaluation index, extract the three-dimensional point set located in the linear contour region of the workpiece from the three-dimensional point cloud data obtained from the actual scanning. The specific calculation process for the point cloud curvature includes, for any point in the point cloud p Searching its neighborhood point set based on KD-tree Using the least squares method to fit this m The quadratic surface at point After solving for the coefficients, the principal curvatures are obtained based on the differential geometric properties of spatial surfaces. k 1, k 2. Then, the mean curvature at that point is calculated. Set a curvature threshold. If the average curvature calculation result is greater than the threshold, it indicates that the point cloud changes significantly at that point and is identified as a data point located in the linear contour area of ​​the scanned workpiece.

[0010] Step S3: Extract single-layer fine contour feature points from the 3D point set of the linear contour region extracted in Step 2 by point cloud thinning. The specific process of point cloud slimming includes any point in the point cloud. p Using a KD-tree search, all nearest neighbor points within a three-dimensional spherical region defined by the search radius are obtained, and calculations are performed. x The two points with the largest differences in distance from the axis to the search point are then compared. A distance threshold is set; if the difference is less than this threshold, the judgment condition is met. y Axial direction or z The same applies to the axial directions; if the criteria are met in at least two directions, then the search point... p These are fine contour feature points.

[0011] Step S4: Downsample the fine contour feature points, and use the downsampled point set containing contour feature information as the key point set for matching the scanned point cloud and the model point cloud; The downsampling calculation process includes creating a three-dimensional voxel mesh based on a fine contour point set, replacing all points within a voxel with the point closest to the voxel center, and using the simplified contour points as the key point set for matching the scanned point cloud with the model point cloud.

[0012] Step S5: Based on the relative positional distribution between key points, statistically analyze the angle, distance, and shape indices between adjacent key points to describe the geometric features between key points; The specific index calculation process for the geometric features between the key points includes, using any key point Centered on a target, a search radius is set to determine all nearest keypoints within the neighborhood. Vectors from the search keypoint to its nearest keypoints are searched, and the angle information between each pair of normalized vectors is calculated to create an angle distribution histogram. The Euclidean distance between the search keypoint and its nearest keypoints is calculated, and a distance distribution histogram is created. The nearest keypoints in the outermost circle within the search radius are extracted, and the geometric shape differences are calculated to describe the spatial distribution of the nearest keypoints.

[0013] Step S6: Optimize the orientation of the normal vector at the key points based on the improved principal component analysis method, construct a local coordinate system and divide a spherical mesh in the coordinate system; The specific process for constructing a local coordinate system based on the improved principal component analysis method includes: Step S321: Set the search radius, based on key points Set the search radius using a KD-tree centered on the target area. r Determine the local neighborhood of the sphere , n This represents the number of points within the neighborhood. Step S322: Update the coordinate system with the key points as the origin. Relative coordinates, eliminating the effects of translation, using vectors It indicates that there are: ; Step S323: Calculate the covariance matrix and perform eigenvalue decomposition, as shown in the following formula: (1) In the formula, l 0, l 1, λ2 are eigenvalues. v 0, v 1, v 2 represents the 3×1 eigenvectors corresponding to the eigenvalues, and the eigenvectors are mutually orthogonal; Step S324: Sort the feature values, assuming ,but v 0 is the direction of the first principal axis ( ); v 1 represents the direction of the second principal axis ( ); v 2 is the direction of the third principal axis ( For the data at the intersection of two planes, geometrically speaking, the X-axis is along the direction of the intersection line, and the Z-axis points to the direction of the angle bisector of the "V"-shaped opening. Step S325: To ensure the Z-axis points are unique and consistent, perform statistical projection and force the Z-axis to point to the side with more points, as shown in the following formula: (2) In the formula, This indicates the final determined Z-axis direction. N + Indicates that the projection satisfies Similarly, the total number of points, N - This indicates the number of points whose projection is less than 0. Step S326: Determined by sorting, Determined by disambiguation, This ensures that the three coordinates are orthogonal and satisfy the right-hand rule, thus obtaining a local coordinate system that can be referenced.

[0014] The spherical mesh segmentation process includes: using key points as centers, based on a pre-constructed local coordinate system and referencing Earth's latitude and longitude, performing azimuth and polar angle segmentation; simultaneously generating multiple concentric spheres at different radial distances; and setting the maximum radius of the radial concentric spheres, which is the search radius of the spherical mesh. R The total number of concentric spheres constructed is L Then the radii of the concentric spheres from the inside out are respectively The azimuth angle is divided into equal parts from 0° to 360°, starting from the positive X-axis direction. J A sector-shaped region; starting from the North Pole (positive Z-axis direction), the range of 0° to 180° is divided into equal parts. K Each region; the final segmentation yields... L × J × K Each grid.

[0015] Step S7: Use Fourier spectrum to describe the azimuth angle distribution change of the spherical grid corresponding to the distribution of adjacent points, and at the same time represent the point distribution in the form of pixel image. Each pixel value corresponds to the number of points in the azimuth-polar angle combined grid, thereby describing the local neighborhood features of the original data at the key point. The specific calculation process for describing local neighborhood features using the Fourier spectrum includes, assuming that the 360° azimuth angle is equally divided into 12 intervals ( J = 12), construct a periodic azimuth signal ,in c 0 indicates the number of points within the 0th azimuth angle interval. The rotated point distribution is as follows: At this point, the Fourier transform spectrum of the original azimuth signal is: The rotated signal spectrum , i 0 represents the phase angle of the translation, but the amplitude spectrum remains unchanged, i.e. Assuming radial stratification I= 3, simultaneously calculate radial layer 0, merge radial layers 0 and 1, merge the data distribution of all polar angle layers of radial layers 0, 1 and 2 at the same azimuth angle, and perform Fourier spectrum transformation to solve the spectrum diagram according to the above steps.

[0016] The specific calculation process for describing local neighborhood features in a pixel image includes sorting the grid in a radial order from inside to outside, an azimuth angle clockwise around the Z-axis, and a polar angle from north to south, denoted as [radial]. r [Azimuth] a [Polar Angle] e According to the sorted spherical grid, for each fixed azimuth-polar angle ( a , e The image features are formed by combining pixels with azimuth as the horizontal axis and polar angle as the vertical axis. Each pixel value corresponds to the number of points within the combined grid. This is assumed to create radially layered structures. L = 3. In order to more comprehensively describe the changes in point distribution, the point distribution image features of radial layers 0, radial layers 0 and 1, and radial layers 0, 1 and 2 are calculated simultaneously.

[0017] Step S8: Combine the features of the model point cloud and the scanned point cloud described by the two dimensions. The geometric histogram features are compared by the Bach coefficient, the Fourier spectrum features are compared by the Pearson correlation coefficient, and the pixel image features are compared by the structural similarity index. Based on the comprehensive similarity, candidate matching key point pairs are determined. The geometric histogram feature comparison process includes, based on the fact that a histogram is a probability distribution, comparing the similarity of two probability distributions using the Bach coefficient (BC), as shown in the following formula: (3) In the formula, k This indicates the number of intervals in the histogram. If the range of 0 to 180° is divided into 12 equal intervals, then... k = 12; H 1( i ) indicates the first i The ratio of the number of occurrences of an angle in a given interval to the total number of angles is the probability distribution. H 2( i Similarly, this refers to comparing the probability information at key points; if the geometric feature histograms at two key points are perfectly similar (the difference is 0), then... .

[0018] The Fourier spectral feature comparison process includes, based on the fact that spectral features are numerical vectors, using the Pearson correlation coefficient. The correlation of the azimuth distribution spectrum maps of different radial layers is compared using the following formula: (4) In the formula, nRepresents the number of azimuth intervals. J ; x i Represents the spectrum F 1 in the i Amplitude values ​​at each frequency point; y i Indicates the spectrum to be compared. F 2 in the i Amplitude values ​​at each frequency point; It is the spectrum F 1. The average amplitude at all frequency points; It is the spectrum F 2. The average amplitude at all frequency points. ;when A result of 1 indicates that the two spectra have exactly the same shape, and the closer the result is to 1, the more similar the two spectra are.

[0019] The pixel image feature comparison process includes, based on the fact that pixel image features are essentially two-dimensional images, using the structural similarity index (SSIM) to measure similarity, as shown in the following formula: (5) In the formula, l ( A , B ) is the brightness contrast function; C ( A , B () is the contrast function; S ( A , B ) is a structure comparison function; m A , m B It is an image A , B Average pixel value; stability constant , W It refers to the dynamic range of pixel values, which is often taken as an empirical value. K 1 = 0.01; s A , s B It is an image A , B Pixel standard deviation; stability constant Experience often takes K 2 = 0.03; s AB It is an image A , BThe pixel covariance; the SSIM index ranges from [-1, 1], and SSIM(A, B) = 1 if and only if image A equals image B. That is, the closer the SSIM index is to 1, the more similar the local neighborhood point distribution images at the two key points are.

[0020] Step S9: Filter candidate matching point pairs by curvature magnitude, determine the correct match, and calculate the coordinate transformation matrix based on the correct match to complete the point cloud registration.

[0021] The correct matching calculation process includes randomly selecting 3 pairs from all candidate matching keypoint pairs, solving for the rotation matrix R and the translation matrix T, and denoting the 3 pairs of candidate matching keypoint pairs as follows: The rotation matrix R can be obtained by following the process of centroid calculation, decentralization, covariance matrix construction, and SVD decomposition. The translation matrix T is obtained by the difference between the centroid rotation transformations of the points. For other key points not involved in the matrix calculation, coordinate transformation is performed according to the solution results. For the transformed key points, the nearest key point in the model point cloud is searched, and the distance is calculated. If the RMSE of the transformed registration error is greater than the threshold, it is determined that there is an incorrect matching point cloud among the three candidate matching point pairs selected at this time, and the candidate point pair library is returned for re-selection and judgment. The rotation and translation matrix determined at this time is considered as the target matrix if and only if the RMSE of the neighborhood transformation of the three key points involved in the calculation of the scanned point cloud and the model point cloud is less than the threshold. The formula for calculating the registration error RMSE is as follows: (6) In the formula, l This indicates the number of key points in the scanned point cloud that were not included in the calculation. d i This indicates the first point cloud after coordinate transformation. i The distance between each neighboring point and the nearest point found in the model point cloud using a KD-tree; default. .

[0022] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention uses contour feature points as registration key points, which can utilize significant structural information such as workpiece holes, steps, boundaries, and surface intersections to reduce redundancy and mismatch when using ordinary surface points or random sampling points for matching.

[0023] 2. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention proposes a GSP feature descriptor that simultaneously includes key point geometric features, Fourier spectrum features and pixel image features, which can describe the morphological features at key points from multiple dimensions and improve the ability to distinguish similar local structures.

[0024] 3. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention uses improved PCA to construct a local reference coordinate system, and reduces the instability of traditional normal estimation at contour points by principal axis sorting and Z-axis direction disambiguation, so that the spherical grid statistics have better directional consistency.

[0025] 4. The local point cloud registration method based on geometry-spectrum-pixel feature descriptor of the present invention describes the azimuth distribution with Fourier amplitude spectrum, so that rotation around the local Z-axis can be equivalent to azimuth signal translation, thereby maintaining the rotation invariance of the descriptor to a certain extent.

[0026] 5. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention further reduces mismatched point pairs through curvature information filtering and optimal matching mechanism, and combines RMSE threshold verification of transformation matrix, which can improve the registration accuracy and reliability under the conditions of poor initial pose and incomplete local point cloud data. Attached Figure Description

[0027] Figure 1 This is an overall flowchart of the local point cloud registration method based on geometry-spectrum-pixel feature descriptors of the present invention; Figure 2 This is a schematic diagram illustrating the workpiece three-dimensional point cloud data scanning experiment conducted in this invention; Figure 3 These are schematic diagrams of three local point cloud models obtained by scanning the workpiece in this invention; Figure 4 A visualization of candidate matching generated by the local point cloud registration method of the present invention and traditional feature descriptors (3DCS, FPFH, SHOT); Figure 5 This is a visualization of the results of registering three local point clouds (green) obtained from workpiece scanning with the model point cloud (blue) using the local point cloud registration method of the present invention and various traditional registration methods. Detailed Implementation

[0028] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0029] The purpose of this invention is to address the shortcomings of existing technologies by providing a local point cloud registration method based on geometric-spectrum-pixel feature descriptors.

[0030] Example 1 The local point cloud registration method based on geometry-spectrum-pixel feature descriptors in this embodiment refers to... Figure 1 As shown, it includes the following steps: (i) Acquiring local point cloud data of the workpiece surface using a 3D scanner The workpiece used in this embodiment is a calibration test piece. The surface of the workpiece contains multiple round holes, countersunk holes, groove structures and planar features, with rich geometric features, making it suitable as a test object for laser scanning point cloud acquisition and model point cloud registration experiments.

[0031] Fix the workpiece position, refer to Figure 2 As shown, a 3D scanner is used to scan and measure the workpiece, obtaining its surface information and storing it as 3D point cloud data. (Refer to...) Figure 3 As shown, in this embodiment, three local point cloud data I, II, and III of the workpiece were measured; The 3D scanner is a non-contact optical scanner that emits structured light for measurement. During operation, the scanner laser emits multiple laser lines to acquire 3D local point cloud data of the surface of the workpiece from different perspectives.

[0032] (II) Extraction of linear contour feature points and construction of key point sets First, using point cloud curvature as an evaluation index, a three-dimensional point set located in the linear contour region of the workpiece is extracted from the data obtained from the actual scanning. The specific calculation process for point cloud curvature includes, for any point in the point cloud... p Searching its neighborhood point set based on KD-tree Using the least squares method to fit this m The quadratic surface at point After solving for the coefficients, the principal curvatures are obtained based on the differential geometric properties of spatial surfaces. k 1, k 2. Then, the mean curvature at that point is calculated. Set a curvature threshold. If the curvature calculation result is greater than the threshold, it means that the point cloud changes significantly at that point and is identified as a data point located in the linear contour area of ​​the scanned workpiece.

[0033] Secondly, using the center point as the linear contour feature point of the transition area at the intersection of the surfaces, single-layer fine contour feature points are extracted by point cloud thinning. The specific process of point cloud slimming includes any point in the point cloud... p Using a KD-tree search, all nearest neighbor points within a three-dimensional spherical region defined by the search radius are obtained, and calculations are performed. xThe two values ​​with the largest differences between the axis direction and the search point are then compared to determine if the difference between these two values ​​is less than a threshold. y Axial direction or z The same principle applies to axis direction; select the decision axis based on the point cloud data, and if the decision criteria are met, the search point is considered valid. p Center point.

[0034] Finally, the fine contour feature points are downsampled, and the downsampled point set containing contour feature information is used as the key point set for matching the scanned point cloud and the model point cloud. The downsampling calculation process includes creating a three-dimensional voxel mesh based on a fine contour point set, replacing all points within a voxel with the point closest to the voxel center, and using the simplified contour points as the key point set for matching the scanned point cloud with the model point cloud.

[0035] (III) Construction of GSP Feature Descriptors (Geometric-Spectrum-Pixel Feature Descriptors) After obtaining the set of contour key points, GSP feature descriptors are further constructed to characterize the contour key points from two dimensions: the relative geometric distribution between key points and the local distribution of the original point cloud in the neighborhood of the key points.

[0036] In terms of geometric feature description, based on the relative positional distribution between key points, the angles, distances, and shape indices between adjacent key points are statistically analyzed to describe the geometric features between key points.

[0037] The calculation process for the geometric features between key points specifically includes, using any key point... Centered on a search radius, all adjacent key points within the neighborhood are determined. Vectors from the search key point to its nearest neighbor are calculated, and the angle information between each pair of normalized vectors is plotted as an angle distribution histogram. The Euclidean distance between the search key point and its nearest neighbor is calculated, and the distance distribution histogram is plotted. The nearest neighbor key points in the outermost circle within the search radius are extracted, and the geometric shape differences are calculated to describe the spatial distribution of the nearest neighbor key points.

[0038] In terms of local neighborhood feature description, the orientation of the normal vector at the key point is first optimized based on the improved principal component analysis method, a local coordinate system is constructed, and a spherical mesh is divided in the coordinate system.

[0039] The specific process for constructing a local coordinate system based on the improved principal component analysis method includes: Step S321: Set the search radius based on the contour key points. Set the search radius using a KD-tree centered on the target area. r Determine the local neighborhood of the sphere , n This represents the number of points within the neighborhood. Step S322: Update the coordinate system with the key points as the origin. Relative coordinates, eliminating the effects of translation, using vectors It indicates that there are: ; Step S323: Calculate the covariance matrix and perform eigenvalue decomposition, as shown in the following formula: (1) In the formula, l 0, l 1, λ2 are eigenvalues. v 0, v 1, v 2 represents the 3×1 eigenvectors corresponding to the eigenvalues, and the eigenvectors are mutually orthogonal; Step S324: Sort the feature values, assuming ,but v 0 is the direction of the first principal axis ( ); v 1 represents the direction of the second principal axis ( ); v 2 is the direction of the third principal axis ( For the data at the intersection of two planes, geometrically speaking, the X-axis is along the direction of the intersection line, and the Z-axis points to the direction of the angle bisector of the "V"-shaped opening. Step S325: To ensure the Z-axis points are unique and consistent, perform statistical projection and force the Z-axis to point to the side with more points, as shown in the following formula: (2) In the formula, This indicates the final determined Z-axis direction. N + Indicates that the projection satisfies Similarly, the total number of points, N - This indicates the number of points whose projection is less than 0. Step S326: Determined by sorting, Determined by disambiguation, This ensures that the three coordinates are orthogonal and satisfy the right-hand rule, thus obtaining a local coordinate system that can be referenced.

[0040] The process of performing spherical mesh segmentation in a pre-constructed local coordinate system includes: using key points as centers, based on the pre-constructed local coordinate system and referencing Earth's latitude and longitude, performing azimuth and polar angle segmentation; simultaneously generating multiple concentric spheres at different radial distances; and setting the maximum radius of the radial concentric spheres, which is the search radius of the spherical mesh. R The total number of concentric spheres constructed is L Then the radii of the concentric spheres from the inside out are respectively The azimuth angle is divided into equal parts, starting from the positive X-axis direction, and (0~360)°. J A sector-shaped region; starting from the North Pole (positive Z-axis direction), (0~180)° is divided into equal parts. K Each region; the final segmentation yields... L × J × K Each grid.

[0041] Based on this, Fourier spectrum is used to describe the azimuth distribution changes of the spherical grid corresponding to the distribution of adjacent points, while the point distribution is represented in the form of pixel images. The specific calculation process of Fourier spectrum describing local neighborhood features includes assuming that the 360° azimuth angle is equally divided into 12 intervals ( J = 12), construct a periodic azimuth signal ,in c 0 indicates the number of points within the 0th azimuth angle interval. The rotated point distribution is as follows: At this point, the Fourier transform spectrum of the original azimuth signal is: The rotated signal spectrum , i 0 represents the phase angle of the translation, but the amplitude spectrum remains unchanged, i.e. Assuming radial stratification I = 3, we simultaneously calculate radial layer 0, merge radial layers 0 and 1, merge the data distribution of all polar angle layers of radial layers 0, 1 and 2 at the same azimuth angle, and perform Fourier spectrum transformation to solve the spectrum diagram according to the above steps.

[0042] The specific calculation process for describing local neighborhood features in a pixel image includes sorting the grid in a radial order from inside to outside, an azimuth angle clockwise around the Z-axis, and a polar angle from north to south, denoted as [radial]. r [Azimuth] a [Polar Angle] e According to the sorted spherical grid, for each fixed azimuth-polar angle ( a , e The image features are formed by combining pixels with azimuth as the horizontal axis and polar angle as the vertical axis. Each pixel value corresponds to the number of points within the combined grid. This is assumed to create radially layered structures. L = 3. In order to more comprehensively describe the changes in point distribution, the point distribution image features of radial layers 0, radial layers 0 and 1, and radial layers 0, 1 and 2 are calculated simultaneously.

[0043] The GSP feature descriptor is constructed through the above process. In this embodiment, based on the density of the point cloud, the search radius for calculating curvature is set to 3 times the average point distance, the search radius for geometric features is set to 5 times the average point distance, the search radius for the local coordinate system is set to 10 times the average point distance, and the azimuth angle is equally divided. J=12, polar angles are equally divided K =6, set the search radius for data I. R =5.0, number of radial layers I =4; Data II settings R =3.5, I =3; Data III Settings R =5.5, I =4.

[0044] (iv) Feature matching and spatial transformation solution based on GSP feature descriptors After constructing the GSP feature descriptor, the similarity between the model point cloud and the scanned point cloud features described by the two dimensions is compared to determine the candidate matching key point pairs. Based on this, the spatial transformation matrix is ​​solved to complete the point cloud registration.

[0045] For feature similarity comparison, different indices are used to compare the three types of features. The geometric histogram feature comparison process includes, based on the fact that a histogram is essentially a probability distribution, comparing the similarity of two probability distributions using the Bach coefficient (BC), as shown in the following formula: (3) In the formula, k This indicates the number of intervals in the histogram. If the range of 0 to 180° is divided into 12 equal intervals, then... k = 12; H 1( i ) indicates the first i The ratio of the number of occurrences of an angle in a given interval to the total number of angles is the probability distribution. H 2( i Similarly, this refers to comparing the probability information at key points; if the geometric feature histograms at two key points are perfectly similar (the difference is 0), then... .

[0046] The Fourier spectral feature comparison process includes, based on the fact that spectral features are numerical vectors, using the Pearson correlation coefficient. The correlation of the azimuth distribution spectrum maps of different radial layers is compared using the following formula: (4) In the formula, n Represents the number of azimuth intervals. J ; x i Represents the spectrum F 1 in the i Amplitude values ​​at each frequency point; y i Indicates the spectrum to be compared. F 2 in the i Amplitude values ​​at each frequency point; It is the spectrum F 1. The average amplitude at all frequency points; It is the spectrum F 2. The average amplitude at all frequency points. ;when A result of 1 indicates that the two spectra have exactly the same shape, and the closer the result is to 1, the more similar the two spectra are.

[0047] The pixel image feature comparison process includes, based on the fact that pixel image features are the essence of two-dimensional images, using the structural similarity index (SSIM) to measure similarity, as shown in the following formula: (5) In the formula, l ( A , B ) is the brightness contrast function; C ( A , B () is the contrast function; S ( A , B ) is a structure comparison function; m A , m B It is an image A , B Average pixel value; stability constant , W It refers to the dynamic range of pixel values, which is often taken as an empirical value. K 1 = 0.01; s A , s B It is an image A , B Pixel standard deviation; stability constant Experience often takes K 2 = 0.03; s AB It is an image A , B The pixel covariance; the SSIM index ranges from [-1, 1], and SSIM(A, B) = 1 if and only if image A equals image B. That is, the closer the SSIM index is to 1, the more similar the local neighborhood point distribution images at the two key points are.

[0048] The overall similarity is determined based on the similarity calculation results of the above three types of features, and candidate matching key point pairs are selected accordingly.

[0049] After candidate matching point pairs are determined, they are filtered by curvature magnitude. The curvature difference at key points is used to determine if it is less than a threshold, eliminating incorrect matches with excessively large curvature differences. Finally, the coordinate transformation matrix is ​​calculated based on the correct matches, thus completing point cloud registration.

[0050] The correct matching calculation process includes randomly selecting 3 pairs from all candidate keypoint matching pairs, solving for the rotation matrix R and the translation matrix T, and denoting the 3 pairs of candidate matching keypoint pairs as follows: The rotation matrix R can be obtained by following the process of centroid calculation, decentralization, covariance matrix construction, and SVD decomposition. The translation matrix T is obtained by the difference between the centroid rotation transformations of the points. For other key points not involved in the matrix calculation, coordinate transformation is performed according to the solution results. For the transformed key points, the nearest key point in the model point cloud is searched, and the distance is calculated. If the RMSE of the transformed registration error is greater than the threshold, it is determined that there is an incorrect matching point cloud among the three candidate matching point pairs selected at this time, and the candidate point pair library is returned for re-selection and judgment. The rotation and translation matrix determined at this time is considered as the target matrix if and only if the RMSE of the neighborhood transformation of the three key points involved in the calculation of the scanned point cloud and the model point cloud is less than the threshold. The formula for calculating the registration error RMSE is as follows: (6) In the formula, l This indicates the number of key points in the scanned point cloud that were not included in the calculation. d i This indicates the first point cloud after coordinate transformation. i The distance between each neighboring point and the nearest point found in the model point cloud using a KD-tree; default. .

[0051] The above process completes the registration of the local scanned point cloud with the model point cloud. In this embodiment, the RMSE threshold is set to 0.5 mm.

[0052] (V) Comparison of traditional methods and the local point cloud registration method of the present invention First, taking Data I as an example, the selection results of matching point pairs are compared between traditional feature descriptors and the GSP feature descriptors of this invention, as shown in Table 1: Table 1. Comparison of matching point pair selection results between traditional feature descriptors and the GSP feature descriptors of this invention.

[0053] In Table 1, Recall is the percentage of point pairs that correctly match TP out of the total number of point pairs that match KP (which is consistent with the number of key points in the scanned point cloud).

[0054] Analysis of the results in Table 1 shows that, in terms of recall, the GSP feature descriptor of this invention has the highest recall rate, followed by the FPFH descriptor, while the SHOT descriptor has the worst matching effect.

[0055] SHOT, based on local features, encodes the statistical histogram of normal features for each point in the search space and then combines them. This approach is not well-suited for local features where contour points are used as keypoints. This is because the local surface variations found in the contour feature region fluctuate significantly, making it difficult to determine an optimal fitting plane. This results in unstable calculations of the normal information for contour points using traditional methods, making it difficult to achieve unified orientation. Therefore, SHOT's judgment based on local neighborhood normal information performs poorly.

[0056] 3DCS descriptors also segment the spherical coordinate system to describe the point distribution in local areas, which can more accurately describe the local neighborhood of contour key points. However, since some local features of local data I have similar features in mismatched areas, it is impossible to completely determine the matching relationship globally based solely on the local neighborhood.

[0057] Compared to SHOT and 3DCS, the FPFH descriptor does not use a segmented spherical coordinate system. It mainly calculates the geometric features of the local neighborhood of the contour key points and combines the weighted sum of adjacent key points, thus it has stronger localization and resolution capabilities and a decent recall rate.

[0058] In contrast, the GSP feature descriptor of this invention not only considers the geometric features of the distribution of contour key points, but also establishes a more stable spherical coordinate system and comprehensively describes the point distribution when calculating the local features of key points, resulting in the highest recall rate and better performance.

[0059] Regarding runtime, the GSP feature descriptor of this invention calculates more comprehensive features, and its runtime is slightly longer than existing descriptors, but still within a reasonable range, with minimal difference. While the SHOT descriptor has the lowest runtime, it also has the worst recall rate. (See reference...) Figure 4 As shown, by observing the overall candidate matching visualization results generated by the descriptors, it can be found that the more unique ellipse and small step area key points are always identified first, and the parts with the most correct matches are the most.

[0060] Secondly, the registration results of traditional point cloud registration methods and the local point cloud registration method based on geometry-spectrum-pixel feature descriptors of this invention are compared. A visualization of the registration results is shown below. Figure 5 As shown in Table 2, the specific data is as follows: α , β , c These correspond to the rotation angles along the X, Y, and Z axes, respectively. X t , Yt , Z t These correspond to the translation amounts in the X, Y, and Z axis directions, respectively.

[0061] Table 2 Comparison of registration results between traditional point cloud registration methods and the method of the present invention

[0062] Analysis of Table 2 reveals that the 4PCS, CPD, and ICP algorithms all have shorter registration runtimes (faster registration speeds). This is because these three algorithms do not require calculating a complete feature descriptor before determining matching point pairs, resulting in generally faster iteration speeds compared to feature descriptor-based point cloud registration algorithms. However, observation shows that the method of this invention is also fast for data II, relatively slower for data I, and slowest for data III. In fact, the main factor affecting the registration speed of the method of this invention is the search radius R. In this embodiment, the search radius R is set to 5.0 for data I, 3.5 for data II, and 5.5 for data III. Therefore, the runtimes for data I and III are relatively longer. Simultaneously, data III requires a larger scanned point cloud with more key points, which slightly affects the algorithm's runtime, but remains within a reasonable range. Comparing the rotation angle, translation amount, and registration accuracy evaluation index RMSE of the coordinate transformation reveals that the method of this invention has the best registration effect and the highest accuracy.

[0063] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A local point cloud registration method based on geometry-spectrum-pixel feature descriptors, characterized in that: include Step 1: Fix the position of the workpiece, use a 3D scanner to scan and measure, obtain the surface information of the workpiece and store it as 3D point cloud data; Step 2: Extract linear contour feature points and construct a contour key point set; Step 3: Construct GSP feature descriptors: After obtaining the set of contour key points, construct GSP feature descriptors to characterize the contour key points from two dimensions: the relative geometric distribution between key points and the local distribution of the original point cloud in the neighborhood of the key points. Step 4: Perform feature matching and spatial transformation solution based on GSP feature descriptors: After completing the construction of GSP feature descriptors, compare the similarity between the model point cloud and the scanned point cloud features described by the two dimensions to determine candidate matching key point pairs, and solve the spatial transformation matrix to complete point cloud registration. Step 3 includes the following steps: Step 31: Based on the relative positional distribution between key points, statistically analyze the angle, distance, and shape indices between adjacent key points to describe the geometric features between key points; Step 32: Optimize the orientation of the normal vector at the key points based on the improved principal component analysis method, construct a local coordinate system and divide a spherical mesh in the coordinate system; Step 33: Use Fourier spectrum to describe the azimuth angle distribution change of the spherical grid corresponding to the distribution of adjacent points, and at the same time, use pixel image to represent the point distribution. Each pixel value corresponds to the number of points in the azimuth-polar angle combination grid, thereby describing the local neighborhood features of the original data at the key points. In step 33, the specific calculation process for describing the local neighborhood features using the Fourier spectrum includes: dividing the 360° azimuth angle into 12 equal intervals to construct a periodic azimuth angle signal and obtain the rotated point distribution; calculating the Fourier transform spectrum of the original azimuth angle signal and the spectrum of the rotated signal; simultaneously calculating the radial layer, merging the data distributions of all polar angle layers of the radial layer at the same azimuth angle, and performing the Fourier spectrum transform to solve the spectrum diagram according to the above steps; The specific calculation process for describing local neighborhood features in a pixel image includes: sorting the grid according to the following order: radially from inside to outside, azimuth angle clockwise around the Z-axis, and polar angle from north to south, denoted as [radial]. r [Azimuth] a [Polar Angle] e Based on the sorted spherical grid, for each fixed azimuth-polar angle combination, with the azimuth angle as the horizontal axis and the polar angle as the vertical axis, each pixel value corresponds to the number of points within the corresponding grid combination, forming image features; radial layering. L = 3, and simultaneously calculate the point distribution image features of radial layers 0, radial layers 0 and 1, and radial layers 0, 1 and 2.

2. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 1, characterized in that, Step 2 includes the following steps: Step 21: Using point cloud curvature as an evaluation index, extract the three-dimensional point set located in the linear contour region of the workpiece from the three-dimensional point cloud data obtained from the actual scanning. Step 22: Extract single-layer fine contour feature points from the 3D point set of the linear contour region extracted in Step 21 by point cloud thinning. Step 23: Downsample the fine contour feature points and use the downsampled point set containing contour feature information as the key point set for matching the scanned point cloud and the model point cloud.

3. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 2, characterized in that, Step 4 includes the following steps: Step 41: Combine the features of the model point cloud and the scanned point cloud described by the two dimensions, compare the geometric histogram features with the Bach coefficient, compare the Fourier spectrum features with the Pearson correlation coefficient, and compare the pixel image features with the structural similarity index. Determine the candidate matching key point pairs based on the comprehensive similarity. Step 42: Filter candidate matching point pairs by curvature magnitude, determine the correct match, and calculate the coordinate transformation matrix based on the correct match to complete the point cloud registration.

4. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 3, characterized in that: In step 21, the calculation process for the curvature of the point cloud includes, for any point in the point cloud... p Searching its neighborhood point set based on KD-tree Using the least squares method to fit this m The quadratic surface at point A After solving for the coefficients, the principal curvatures are obtained based on the differential geometric properties of spatial surfaces. k 1, k 2. Then, the mean curvature at that point is calculated. ; A curvature threshold is set. If the average curvature calculation result is greater than the threshold, it indicates that the point cloud changes significantly at that point and is identified as a data point located in the linear contour area of ​​the scanned workpiece. In step 22, the point cloud slimming process includes any point in the point cloud... p Using a KD-tree search, all nearest neighbor points within a three-dimensional spherical region defined by the search radius are obtained, and calculations are performed. x The two points with the largest differences in distance from the axis to the search point are then compared. A distance threshold is set; if the difference is less than this threshold, the judgment condition is met. y Axial direction or z The same applies to the axial direction; If the decision condition is met in at least two directions, then the search point... p For fine contour feature points; In step 23, the downsampling calculation process includes creating a three-dimensional voxel mesh based on a fine contour point set, replacing all points within a voxel with the point closest to the voxel center, and using the simplified contour points as the key point set for scanning the point cloud and matching the model point cloud.

5. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 4, characterized in that: In step 31, the specific index calculation process for the geometric features between key points includes setting a search radius with any key point as the center and determining all neighboring key points in the neighborhood. Search for the vectors from keypoints to nearest keypoints, calculate the angle information between all pairs of normalized vectors, and draw an angle distribution histogram. Calculate the Euclidean distance between the search keypoint and its nearest neighbor keypoint, and draw a distance distribution histogram; extract the nearest neighbor keypoints in the outermost circle within the search radius, and calculate the geometric shape differences to describe the spatial distribution of the nearest neighbor keypoints.

6. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 5, characterized in that: In step 32, the specific process for constructing a local coordinate system based on the improved principal component analysis method includes: Step S321: Set the search radius using a KD-tree, centered on the key point. r Determine the local neighborhood points of the sphere; Step S322: Update the relative coordinates with the key point as the origin to eliminate the effect of translation; Step S323: Calculate the covariance matrix and perform eigenvalue decomposition; Step S324: Sort the feature values; Step S325: To ensure that the Z-axis points are unique and consistent, perform statistical projection and force the Z-axis to point to the side with more points distributed. Step S326: Obtain a local coordinate system that can be referenced; The spherical mesh segmentation process includes: using key points as the center, based on the constructed local coordinate system and referencing Earth's latitude and longitude, performing azimuth and polar angle segmentation; simultaneously generating multiple concentric spheres at different radial distances; and setting the maximum radius of the radial concentric spheres, which is the search radius of the spherical mesh. R The total number of concentric spheres constructed is L Then the radii of the concentric spheres from the inside out are respectively The azimuth angle is divided into equal parts from the positive X-axis direction, with 0° to 360°. J A sector-shaped region; the polar angle starts from the North Pole, dividing 0° to 180° into equal parts. K Each region; the final segmentation yields... L × J × K Each grid.

7. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 6, characterized in that: In step 41, the geometric histogram feature comparison process includes: based on the fact that a histogram is a probability distribution, comparing the similarity of two probability distributions through the Bach coefficient BC; The Fourier spectral feature comparison process includes: based on the fact that spectral features are numerical vectors, using the Pearson correlation coefficient (PCC) to compare the correlation of azimuth distribution spectral maps of different radial layers; The pixel image feature comparison process includes: based on the fact that pixel image features are the essence of two-dimensional images, the structural similarity index (SSIM) is used to measure similarity.

8. The local point cloud registration method based on geometry-spectrum-pixel feature descriptors according to claim 7, characterized in that: In step 42, the correct matching calculation process includes: randomly selecting from all candidate keypoint matching pairs. i Solve for the rotation matrix R and translation matrix T using a set of methods. i Group candidate matching key point pairs The rotation matrix R can be obtained by following the process of centroid calculation, decentralization, covariance matrix construction, and SVD decomposition; the translation matrix T is obtained by the difference between points after centroid rotation transformation; for other key points not involved in matrix calculation, coordinate transformation is performed according to the solution results; for the transformed key points, the nearest key point in the model point cloud is searched, and the distance is calculated; if the registration error RMSE after transformation is greater than a threshold, then the selected key point is determined to be the nearest key point. i If an incorrectly matched point cloud exists in the candidate point pair group, return to the candidate point pair library for re-selection and judgment; this applies if and only if the scanned point cloud is involved in the calculation. i If the RMSE of the neighborhood transformation at each key point and the model point cloud is less than the threshold, then the rotation and translation matrix determined at this time is considered the target matrix.

Images