Urban road point cloud registration method based on salient key point extraction

Abstract

An urban road point cloud registration method based on salient key point extraction comprises the steps of taking an obtained urban road point cloud as a processing object, performing ground removal operation on data, performing down-sampling by adopting a self-defined sparse processing mode, and filtering dense regions according to the number of adjacent points of the point cloud; calculating a normal vector of the point cloud by adopting an improved adaptive normal; key points are extracted from the two pieces of point cloud data through a significant corner extraction algorithm; a fast point histogram FPFH is used for carrying out feature calculation on key points, initial transformation is calculated in combination with random sample consensus RANSAC, coarse registration is completed, and finally iterative optimization processing is carried out in combination with an average matching distance MMD and a trimmed nearest point iteration TrICP algorithm. According to the method, key points are extracted by adopting a significant angular point extraction algorithm, and points with obvious normal vector changes on curbs, street lamps and billboards can be accurately identified; redundant calculation is avoided, registration time is shortened, and the problems that road point cloud registration is tedious and low in efficiency are solved.

Description

Technical Field

[0001] The present invention relates to the technical field of three-dimensional reconstruction in vehicle autonomous driving, and specifically relates to a method for registering point clouds of urban roads based on the extraction of significant key points. Background Art

[0002] In the field of autonomous driving, the biggest current problem is to match the point cloud of the current frame with a pre-constructed high-precision map in order to achieve precise vehicle positioning. On the other hand, the decision-making and planning of the autonomous driving system are also very important. Through accurate pose estimation, the vehicle can better understand its direction and position in the environment, and thus make safer and more effective driving decisions. For the safety and efficiency of autonomous driving, it is necessary to ensure real-time communication and collaborative work between the vehicle and road infrastructure. In the field of digital smart city construction, there are difficulties in processing large-scale, multi-source, and multi-temporal large-scale scene data, and these data usually come from different sensors and platforms. The core purpose is to accurately fuse these heterogeneous data to generate a complete three-dimensional model of the city, which is of great significance for urban planning, management, and the construction of smart cities.

[0003] Laser scanning technology has the characteristics of strong initiative, good real-time performance, high precision, little influence by the environment, all-weather data acquisition, and full digitalization of the spatial information obtained. It can directly obtain high-precision three-dimensional point information of the scene in the city and has become an important means for collecting three-dimensional spatial data of large-scale urban scenes. It can provide powerful data support for tasks such as model construction, virtual simulation, and scene analysis in the construction of digital twin cities. Point cloud registration can splice three-dimensional point clouds collected from different perspectives, with different characteristics, and different platforms together to obtain complete three-dimensional scene information in the same coordinate system. It is a key step in tasks such as point cloud processing, model construction, and data analysis, and can provide strong technical support for obtaining complete three-dimensional information of the city.

[0004] However, due to the particularity of large-scale scenes such as urban roads, there are spatial deviations in the point cloud data collected at different times, locations, and by different devices, making it difficult to directly fuse and analyze these discrete data. Radar point cloud data registration has become a key link in solving this problem. Specifically, dynamic objects (such as pedestrians and vehicles) in the urban environment will cause outliers in the point cloud, resulting in low registration accuracy. In addition, due to limited perspectives and occlusion during data collection, there are few overlapping areas in the point cloud and the initial poses of the point clouds are quite different. Therefore, an efficient and accurate point cloud registration method is very important.

[0005] The registered data can accurately align the point clouds collected at different times and positions, forming a continuous and high-precision three-dimensional scene model. Therefore, solving the registration problem of lidar point cloud data in large scenes has important strategic value for improving the urban road management level and the reliability of the urban traffic system. Thus, proposing an accurate and efficient point cloud registration strategy is of top priority. Summary of the Invention

[0006] Aiming at the technical problems of low accuracy and slow efficiency in the existing large-scale road point cloud registration methods, this technical solution provides an urban road point cloud registration method based on the extraction of significant key points. The significant corner point extraction algorithm is used to extract key points, which can accurately identify the points with obvious normal vector changes on the curb, street lights, and signs; it avoids redundant calculations, reduces the registration time, and solves the problems of cumbersome road point cloud registration and slow efficiency; it can effectively solve the above problems.

[0007] The present invention is realized through the following technical solutions:

[0008] An urban road point cloud registration method based on the extraction of significant key points, comprising the steps:

[0009] Step 1: Take the obtained urban road point cloud as the processing object, perform ground removal on the data, and perform downsampling using a custom sparse processing method, and filter the dense areas according to the number of neighboring points of the point cloud;

[0010] Step 2: Calculate the normal vector of the point cloud using an improved adaptive normal method; the implementation process of the improved adaptive normal calculation is as follows:

[0011] Step 2.1: Use the adaptive normal vector calculation function to dynamically adjust the neighborhood radius to calculate the normal vector of the point cloud. Set the initial search radius base_radius to 0.05. If the number of neighborhood points is less than 10, gradually increase the search radius, increasing by a set multiple each time, and at most 5 times; and add a minimum neighbor point constraint. If the number of neighborhood points is less than 5, assign the default normal vector [0, 0, 1] to prevent abnormal calculations;

[0012] Step 2.2: Calculate the covariance matrix for the points that meet the neighbor point requirements using the method combined with the principal component analysis PCA. Taking the source point cloud P S as an example, P S ={p1, p2,..., p n}, the coordinates of each point p i are (x i , y i , z i ). Calculate the centroid of the neighborhood of each point. The formula is as follows:

[0013]

[0014] where is the calculated centroid, n is the number of neighborhood points, and x i , y i , z i are the x-axis, y-axis, and z-axis coordinates corresponding to each point respectively; i is the i-th point in the neighborhood;

[0015] Calculate the covariance matrix C, the formula is as follows:

[0016]

[0017] where, σ xx represents the covariance of the x-axis coordinate of a certain neighborhood point of a certain point p i , σ xy represents the covariance of the x-axis coordinate and the y-axis coordinate of the neighborhood point of a certain point p i , σ xz represents the covariance of the x-axis coordinate and the z-axis coordinate of the neighborhood point of a certain point p i , σ yx represents the covariance of the y-axis coordinate and the x-axis coordinate of the neighborhood point of a certain point p i , σ yy represents the covariance of the y-axis coordinate and the y-axis coordinate of the neighborhood point of a certain point p i , σ yz represents the covariance of the y-axis coordinate and the z-axis coordinate of the neighborhood point of a certain point p i , σ zx represents the covariance of the z-axis coordinate and the x-axis coordinate of the neighborhood point of a certain point p i , σ zy represents the covariance of the z-axis coordinate and the y-axis coordinate of the neighborhood point of a certain point p i , σ zz represents the covariance of the z-axis coordinate and the z-axis coordinate of the neighborhood point of a certain point p i ;

[0018] The calculation method of each element is as follows:

[0019]

[0020] where a, b represent the coordinate axes (x, y, z), is the centroid coordinate of the corresponding axis, and n is the number of neighborhood points;

[0021] Calculate the eigenvalues λ1, λ2, λ3 and the corresponding eigenvectors through eigenvalue decomposition, where the eigenvector corresponding to the smallest eigenvalue is the normal vector, and the formula is as follows:

[0022] Cv = λv (4);

[0023] Where C is the covariance matrix solved in the previous step, v is the eigenvector corresponding to each λ and is also the required normal vector, and λ is the eigenvalue;

[0024] Step 2.3: Adjust the search radius again in combination with the curvature information to optimize the normal vector calculation. The curvature calculation formula is as follows:

[0025] curvature = λ min / (λ sum + ∈) (5);

[0026] Where curvature is the curvature, λ min is the minimum eigenvalue obtained from the covariance matrix, λ sum is the sum of the three eigenvalues, and ∈ is a very small number to prevent division by zero;

[0027] The radius adjustment formula is:

[0028] radius = base_radius × (1 + curvature × b) (6);

[0029] Where radius is the required radius, base_radius is the initially set radius, curvature is the curvature, and b in the above formula is the amplification factor, set to 10; when the detected curvature is large, the radius increases to facilitate capturing more neighborhood points to make the normal vector calculation more accurate;

[0030] Step 2.4: Unify the direction of the calculated normal vectors. Calculate the vector from each point to the centroid. If the angle between the normal vector of the point and the vector from the point to the centroid is greater than 90°, then flip the normal vector. This step ensures the consistency of the normal vector direction;

[0031] Step Three: Use a significant corner extraction algorithm to extract key points from the two point cloud data respectively;

[0032] Step Four: Use the method of Fast Point Feature Histogram (FPFH) to calculate the features of the key point neighborhood, that is, the feature descriptor;

[0033] Step Five: Perform rough point cloud registration, and use Random Sample Consensus (RANSAC) to estimate the initial transformation matrix between the source point cloud and the target point cloud;

[0034] Step Six: Perform precise point cloud registration, use TrICP for the final optimization work, remove incorrect point pairs, and use the Mean Matching Distance (MMD) as a metric to calculate the deviation between the two point clouds for iterative optimization.

[0035] Further, the specific operation method of the first step is as follows: Since the density of trees is high, while the point cloud density of road signs, road markers or other signs is relatively uniform, the K-dimensional tree is used for neighborhood search to calculate the number of neighbors of each point. Calculate the number of neighbors of each point within a radius of about 0.35. If the number of neighbor points is less than the threshold, the point is regarded as a sparse point; otherwise, it is regarded as a dense point. Eliminate the dense points and perform cyclic filtering until the number of dense points is low, and then stop the filtering.

[0036] Further, the specific operation method of the significant corner point extraction algorithm in step three for extracting key points is as follows:

[0037] Step 3.1: Input two adjacent point clouds and the corresponding normal vectors calculated in step two respectively, and set the neighborhood search radius, non-maximum suppression, and threshold detection parameters;

[0038] Step 3.2: Obtain the eigenvalues of the covariance matrix of the points calculated at different scales s through calculation formulas (2), (3), and (4), which are used as part of the key point response formula, and normalize the obtained eigenvalue responses. The formulas are as follows:

[0039]

[0040] where R eig is the required eigenvalue response, λ1, λ2, λ3 are the eigenvalues obtained from the covariance matrix, and ∈ is a very small number to prevent division by zero errors;

[0041] Step 3.3: Calculate the change amount of the included angle between the normal vectors of the neighborhood points in the point cloud at scale s, which is used as another part of the key point response formula. The formula is as follows:

[0042]

[0043] In the above formula, for a certain point p i with the normal vector n i , for the neighborhood point p i of point p j , calculate the normal vector n j of the neighborhood point p j , θ ij is the included angle between a certain point p i and the neighborhood point p j , n i and n j are the normal vectors corresponding to the two points, n i ·n j is the dot product, and ||·|| is the norm; and normalize the above formula to obtain the normal vector response: 4444

[0044]

[0045] In the above formula, R angle is the required normal vector response, Ν(p i ) is the total number of neighborhood points of a certain point p i , and the neighborhood points of point i are p j , 1 - cosθ ij reflects the amplitude of the normal vector change;

[0046] Step 3.4: Add a scale - adaptive weight. First, calculate the average distance of neighborhood points:

[0047]

[0048] In the above formula, d avg,s is the average distance of neighborhood points at a certain scale s, Ν s (p i ) is the total number of neighborhood points of a certain point p i at a certain scale s, the neighborhood points are p j , D is the Euclidean distance between a certain point p i and its neighborhood point p j ;

[0049] Then calculate the scale - adaptive weighting factor:

[0050]

[0051] where, w s is the scale factor, γ is a hyper - parameter controlling scale adaptability, d avg,s is the average distance of neighborhood points at a certain scale s, and e is a constant;

[0052] Step 3.5: The final key - point response formula combines the results of Step 3.2, Step 3.3, and Step 3.4, calculates the response values at multiple scales, and takes the maximum value. The formula is as follows:

[0053]

[0054] where, R i is the final response formula, s represents a certain scale, N(s) represents all scales, w s is the scale weight of a certain scale, R eig,s is the eigenvalue response at a certain scale, R angle,s is the normal vector response at a certain scale, α is the weight controlling the eigenvalue response, β is the weight controlling the normal vector change response, calculate the response at different radii s and take the maximum value to improve scale invariance, and introduce w sWeighting is performed to make the response more robust to changes in point cloud density; the key point response formula is calculated for the two point clouds respectively, and if the result is greater than the threshold, it is classified as a key point.

[0055] Furthermore, the calculation of FPFH in step 4 is based on the point feature histogram PFH, which calculates the features of each point in the local neighborhood; the specific operation method in step 4 is as follows:

[0056] Step 4.1: For each key point p i , use its normal vector n i to define the local coordinate system:

[0057] The axis u that is the same as the normal vector:

[0058] u = n i (13);

[0059] where, n i is the normal vector of a certain point p i ;

[0060] Calculate the vector d from the neighborhood point p j to p i and orthogonalize it, that is, the axis v;

[0061] d = p j - p i (14);

[0062] where, d is the vector from the neighborhood point p j to the point p i ;

[0063]

[0064] where, d is the vector from the neighborhood point p j to the point p i , n i is the normal vector of a certain point p i , and here it is represented by the u axis;

[0065] Calculate the axis w:

[0066] w = u × v (16);

[0067] where, the u axis represents the normal vector n i of a certain point p i , v is the vector from the orthogonalized neighborhood point p j to p i ;

[0068] For each neighborhood point p j , calculate the following 3 rotation-invariant features:

[0069] Normal vector angle f1:

[0070]

[0071] where w j is obtained by the cross product of the normal vector of a certain point p i and the vector from the orthonormalized neighboring point p j to p i , as shown in formula (16), and n j is the normal vector of the neighboring point p j , and v j is calculated by formula (15);

[0072] Normal vector angle f2:

[0073] f2 = arccos(n i ·n j ) (18);

[0074] where n i , n j are the normal vectors of a certain point p i and the neighboring point p j ;

[0075] Angle f3 between the normal vector and the connecting vector:

[0076]

[0077] where n j is the normal vector of the neighboring point p j , v j is the vector from the orthonormalized neighboring point p j to p i , d is the vector from the neighboring point p j to p i ;

[0078] For each key point p i , calculate the distributions of (f1, f2, f3) for all p j within its neighborhood to obtain the SPFH histogram;

[0079] Step 4.2: Improve the local description by weighting the SPFH of neighboring points:

[0080]

[0081] where Ν(p i ) is the neighborhood of point p i , and w ij is the distance weight between the point pair (p i , p j ), and the formula is as follows:

[0082] w ij = ||D|| + ε (21);

[0083] Where D is the Euclidean distance between a certain point i and its neighboring point j, and ε is a very small number to prevent a zero result. The final result is the combination of the SPFH of each point itself and the weighted SPFH of the neighboring points, which improves the robustness and discrimination ability.

[0084] Furthermore, the specific operation method of the fifth step is as follows: Based on the key points extracted in the above steps and the feature descriptions of the key points, perform corresponding point matching between the source point cloud and the target point cloud. Since the overlap rate of the two point clouds is low, use the Random Sample Consensus (RANSAC) method to initially calculate the transformation matrix to align the source point cloud with the target point cloud as much as possible, achieving rough registration of the two point clouds; calculate the point cloud features, establish the feature matching relationship between the source point cloud and the target point cloud, randomly sample 3 feature points from the source point cloud, find the corresponding target points, and calculate the rigid transformation matrix of rotation R and translation t to align the sampled 3 points with the target point cloud as much as possible; then calculate the number of inliers, that is, calculate whether the transformed points fall within the reasonable range of the target point cloud, and the threshold of the target point cloud can be customized, and count the number of inliers; repeat the above operation steps 50,000 times, and select the transformation with the largest number of inliers as the final transformation result.

[0085] Furthermore, the rigid transformation matrix needs to be calculated based on the rigid transformation model, and the rigid transformation model is:

[0086] p i ' = Rp i + t (22);

[0087] Where p i is the point before transformation in the source point cloud, p' i is the point after transformation in the source point cloud, R ∈ R 3×3 is the rotation matrix, and t ∈ R 3×1 is the translation vector, provided that R and t are calculated;

[0088] Let the point p i in the source point cloud and the point q i in the target point cloud, calculate the covariance matrix:

[0089]

[0090] Where p i is a certain point in the source point cloud, q i is a certain point in the target point cloud, and respectively represent the centroids, and M is the total number of points in the source point cloud;

[0091] Calculate the Singular Value Decomposition (SVD):

[0092] UΣV T = SVD(H) (24);

[0093] where U, Σ, and V are all 3×3 matrices, and the elements on the diagonal of the obtained matrix Σ are the singular values of matrix H, and T is the transpose of the matrix;

[0094] Calculate the rotation matrix and translation respectively:

[0095] R = VU T (25);

[0096] where R is the required rotation matrix, U and V are the 3×3 matrices obtained above, and T is the transpose of the matrix;

[0097]

[0098] where R is the required rotation matrix, U and V are the 3×3 matrices obtained above, t is the translation vector, and respectively represent the centroids of the source point p i and the target point q i ;

[0099] The final transformation matrix is:

[0100]

[0101] where T r is the required transformation matrix.

[0102] Furthermore, the operation method for precise point cloud registration described in Step 6 includes the steps:

[0103] Use the Mean Matching Distance (MMD) as the basis for optimization convergence, iteratively optimize the TrICP (Trimmed Iterative Closest Point algorithm) for fine registration, achieve optimized alignment based on the result of coarse registration, and solve the transformation matrix by minimizing the error:

[0104]

[0105] where R ∈ R 3×3 is the rotation matrix, t ∈ R 3×1 is the translation vector, p' i is the point obtained after coarse registration of the source point cloud, q i is a point in the target point cloud, and M is the number of points in the source point cloud; abnormal matching points are removed in each iteration, and only 50% of the best matching points are retained to prevent the influence of incorrect matching on the result; the Mean Matching Distance (MMD) is calculated after each iteration to monitor whether the error converges:

[0106]

[0107] where M is the number of points in the source point cloud, p' i is a point after the source point cloud is transformed, and q i is the point in the target point cloud closest to p' i and ||·|| represents the Euclidean distance;

[0108] |MMD prev - MMD current | < threshold (30);

[0109] where MMD current is the calculated distance for the current iteration, MMD prev is the result of the previous iteration, and threshold is a user-defined threshold. If MMD converges, the iteration stops.

[0110] Beneficial effects

[0111] A method for registering urban road point clouds based on the extraction of significant key points proposed by the present invention has the following beneficial effects compared with the prior art:

[0112] (1) In this technical solution, the obtained urban road point cloud is used as the processing object. To improve the efficiency of processing data, the ground is first removed from the data, and then a user-defined sparsification process is used for downsampling, and the dense areas are filtered according to the number of neighboring points. Since there are too many similar objects in the point cloud data, such as trees, a significant corner extraction algorithm is used to extract key points. This method can accurately identify points with obvious changes in the normal vector on road signs such as curbs, street lights, and billboards. After the point cloud is transformed, the repetition rate of the key points is 70%. The features of the key points are also calculated using the Fast Point Feature Histogram (FPFH), and the initial transformation is calculated in combination with the Random Sample Consensus (RANSAC). Coarse registration is completed, and finally, iterative optimization processing is performed in combination with the Mean Matching Distance (MMD) and Trimmed Iterative Closest Point (TrICP) algorithms; redundant calculations are avoided, the registration time is reduced, and the problem of cumbersome and slow road point cloud registration is solved.

[0113] (2) The present invention proposes a radius-based sparse and dense downsampling strategy. This method removes the interference of ground data while retaining the characteristics of different object densities and shapes in the data, laying a good foundation for subsequent processing work, making the operation simple, and enabling fast downsampling. The points in the sparse area are distinguished by a specified threshold radius, and circular downsampling is performed until a small number of dense points are obtained. Since the amount of urban point cloud data is very large, data preprocessing work must be carried out, and the radius-based sparse and dense downsampling strategy can reduce the amount of data while retaining the characteristics of each area of the point cloud, avoiding filtering sparse but characteristic points such as poles like street lights.

[0114] (3) The adaptive normal calculation in the present invention fuses and dynamically adjusts the search neighborhood radius, uses the principal component analysis (PCA) to calculate the normal vector, and also combines the curvature information to optimize the result again, avoiding the operation of manually setting the radius. The subsequent key point extraction depends on efficient and accurate normal vector calculation.

[0115] (4) The method for extracting 3D key points of point cloud in the present invention aims at the problems of many similar objects and small data overlap in urban point cloud. With the idea of only grasping the key points, it restricts the key point response formula, so as to only focus on the obvious corner points and inflection points in the point cloud, avoiding the detection of unreliable and useless points and improving the matching performance of key points.

[0116] (5) The present invention optimizes the Trimmed Iterative Closest Point algorithm (TrICP), uses the Mean Matching Distance (MMD) as the error metric to calculate the registration deviation between point clouds for iterative optimization until the iteration converges. BRIEF DESCRIPTION OF THE DRAWINGS

[0117] Figure 1 is a schematic diagram of the overall process of the present invention.

[0118] Figure 2 is a preliminary effect diagram after downsampling of the present invention.

[0119] Figure 3 is a final effect diagram after downsampling of the present invention.

[0120] Figure 4 is a schematic diagram of key point extraction in the present invention.

[0121] Figure 5 is a point cloud registration result diagram of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0122] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. The described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Without departing from the design concept of the present invention, various modifications and improvements made by those of ordinary skill in the art to the technical solutions of the present invention shall fall within the protection scope of the present invention.

[0123] Embodiment 1:

[0124] As Figure 1 shown, a method for registering urban road point cloud based on significant key point extraction, the specific operation method includes the steps:

[0125] Step 1: Take the acquired urban road point cloud as the processing object, perform ground removal on the data, downsample it using a custom sparse processing method, and filter dense areas based on the number of neighboring points of the point cloud; input the acquired point cloud data, i.e., the source point cloud P S and the target point cloud P T . The specific operation method is as follows:

[0126] Since the number of points in the scanned point cloud data is large, it is necessary to perform preprocessing operations on the data - downsampling. Due to the high density of trees, while the point cloud densities of other objects such as road signs and street lights are relatively uniform, and to avoid most of the subsequent key point extractions being concentrated on the point cloud in the leaf area, it is necessary to perform special processing on the point cloud according to different densities - a custom sparse filtering method. Use the K-dimensional tree (KdTree) for neighborhood search to calculate the number of neighbors of each point, calculate the number of neighbors within a certain radius for each point (in this embodiment, about 0.35 is used). If the number of neighboring points is less than the threshold, the point is regarded as a sparse point; otherwise, it is regarded as a dense point, and the dense points are removed; and the loop filtering is performed until the number of dense points is relatively low and then the filtering stops. Figure 2 The figure shows the effect diagram after one filtering. The blue part is the retained part, and the green points are removed. It can be seen that most of the points slightly higher than the ground are removed, Figure 3 which is the final effect diagram, and there are very few green points.

[0127] Step 2: Calculate the normal vector of the point cloud using an improved adaptive normal method; the implementation process of the improved adaptive normal calculation is as follows:

[0128] Step 2.1: Use the adaptive normal vector calculation function to dynamically adjust the neighborhood radius to calculate the normal vector of the point cloud. Set the initial search radius base_radius to 0.05. If the number of neighborhood points is less than 10, gradually increase the search radius, each time increasing by a set multiple, and at most 5 times; and add a minimum neighbor point constraint. If the number of neighborhood points is less than 5, assign the default normal vector [0, 0, 1] to prevent abnormal calculations.

[0129] Step 2.2: For points that meet the neighbor point requirements, use the method combining principal component analysis PCA to calculate the covariance matrix. Taking the source point cloud P S as an example, P S = {p1, p2, …, p n}, the coordinates of each point p i are (x i , y i , z i ). Calculate the centroid of the neighborhood of each point. The formula is as follows:

[0130]

[0131] Among them, is the calculated centroid, n is the number of neighborhood points, and x i , y i , z i are the x-axis, y-axis, and z-axis coordinates corresponding to each point respectively; i is the i-th point in the neighborhood.

[0132] Calculate the covariance matrix C, and the formula is as follows:

[0133]

[0134] Among them, σ xx represents the covariance of the x-axis coordinate of a certain neighborhood point of a certain point p i , σ xy represents the covariance of the x-axis coordinate and the y-axis coordinate of the neighborhood points of a certain point p i , σ xz represents the covariance of the x-axis coordinate and the z-axis coordinate of the neighborhood points of a certain point p i , σ yx represents the covariance of the y-axis coordinate and the x-axis coordinate of the neighborhood points of a certain point p i , σ yy represents the covariance of the y-axis coordinate and the y-axis coordinate of the neighborhood points of a certain point p i , σ yz represents the covariance of the y-axis coordinate and the z-axis coordinate of the neighborhood points of a certain point p i , σ zx represents the covariance of the z-axis coordinate and the x-axis coordinate of the neighborhood points of a certain point p i , σ zy represents the covariance of the z-axis coordinate and the y-axis coordinate of the neighborhood points of a certain point p i , σ zz represents the covariance of the z-axis coordinate and the z-axis coordinate of the neighborhood points of a certain point p i .

[0135] Among them, the calculation method of each element is:

[0136]

[0137] Among them, a and b represent the coordinate axes (x, y, z), is the centroid coordinate of the corresponding axis, and n is the number of neighborhood points.

[0138] Calculate the eigenvalues λ1, λ2, λ3 and the corresponding eigenvectors through eigenvalue decomposition. Among them, the eigenvector corresponding to the smallest eigenvalue is the normal vector, and the formula is as follows:

[0139] Cv = λv (4);

[0140] Among them, C is the covariance matrix solved in the previous step, c is the eigenvector corresponding to each λ and is also the required normal vector, and λ is the eigenvalue.

[0141] Step 2.3: Combine the curvature information to adjust the search radius again and optimize the normal vector calculation. The curvature calculation formula is as follows:

[0142] curvature = λ min / (λ sum + ∈) (5);

[0143] Among them, curvature is the curvature, λ min is the minimum eigenvalue obtained from the covariance matrix, λ sum is the sum of the three eigenvalues, and ∈ is a very small number to prevent division by zero operation.

[0144] The radius adjustment formula is:

[0145] radius = base_radius × (1 + curvature × b) (6);

[0146] Among them, radius is the required radius, base_radius is the initially set radius, curvature is the curvature, and b in the above formula is the amplification factor, which is set to 10; when the detected curvature is large, the radius increases to facilitate capturing more neighborhood points to make the normal vector calculation more accurate.

[0147] Step 2.4: Make the calculated normal vector directions consistent. Calculate the vector from each point to the centroid. If the angle between the normal vector of this point and the vector from the point to the centroid is greater than 90°, then flip the normal vector. This step ensures the consistency of the normal vector directions.

[0148] Step Three: Use the salient corner extraction algorithm to extract key points from the two point cloud data respectively. The specific operation method is as follows:

[0149] Step 3.1: Input two adjacent point clouds respectively, as well as the corresponding normals calculated in Step Two, and set the neighborhood search radius, non-maximum suppression, and threshold detection parameters.

[0150] Step 3.2: Through calculation formulas (2)(3)(4), obtain the eigenvalues of the covariance matrix of the points calculated at different scales s, which are used as part of the key point response formula, and normalize the obtained eigenvalue responses. The formula is as follows:

[0151]

[0152] Among them, R eigis the required eigenvalue response, λ1, λ2, λ3 are the eigenvalues obtained from the covariance matrix, and ∈ is a very small number to prevent division by zero errors.

[0153] Step 3.3: Calculate the change in the angle between the normal vectors of the neighborhood points in the point cloud at scale s, which is another part of the key point response formula. The formula is as follows:

[0154]

[0155] In the above formula, for a certain point p i with the normal vector n i , for the neighborhood point p i of point p j , calculate the normal vector n j of the neighborhood point p j , θ ij is the angle between a certain point p i and its neighborhood point p j , n i and n j are the corresponding normal vectors of the two points, n i ·n j is the dot product, and ||·|| is the norm. Normalize the above formula (8) to obtain the normal vector response:

[0156]

[0157] In the above formula, R angle is the required normal vector response, Ν(p i ) is the total number of neighborhood points of a certain point p i , the neighborhood points of point p i are p j , 1 - cosθ ij reflects the magnitude of the normal vector change.

[0158] Step 3.4: Add a scale adaptive weight. First, calculate the average distance of the neighborhood points:

[0159]

[0160] In the above formula, d avg,s is the average distance of the neighborhood points at a certain scale s, Ν s (p i ) is the total number of neighborhood points of a certain point p i at a certain scale s, the neighborhood points are p j , D is the Euclidean distance between a certain point p i and its neighborhood point p j .

[0161] Then calculate the scale adaptive weighting factor:

[0162]

[0163] Among them, w s is the scale factor, γ is the hyperparameter controlling scale adaptability, d avg,s is the average distance of neighborhood points at a certain scale s, and e is a constant.

[0164] Step 3.5: The final key point response formula combines the results of Step 3.2, Step 3.3, and Step 3.4, calculates the response values at multiple scales, and takes the maximum value. The formula is as follows:

[0165]

[0166] Among them, R i is the final response formula, s represents a certain scale, N(s) represents all scales, w s is the scale weight at a certain scale, R eig,s is the eigenvalue response at a certain scale, R angle,s is the normal vector response at a certain scale, α is the weight controlling the eigenvalue response, β is the weight controlling the normal vector change response, calculate the response at different radii s and take the maximum value to improve scale invariance, introduce w s for weighting, making the response more robust to changes in point cloud density; calculate the key point response formula for the two point clouds respectively, and if the result is greater than the threshold, it is classified as a key point. Figure 4 The key point extraction effect of the source point cloud shown, the yellow points are the key points, and the corner points and inflection points with obvious neighborhood feature changes can be located.

[0167] Step Four: Use the method of Fast Point Feature Histogram (FPFH) to calculate the features of the key point neighborhood, that is, the feature descriptor; having the key points cannot yet connect the two point clouds. The key points are the highly representative points found among a large number of points, and what actual features this point has requires calculating the information around the point. In this embodiment, the method of Fast Point Feature Histogram (FPFH) is used to calculate the features of the key point neighborhood, that is, the feature descriptor. The calculation of FPFH is based on the Point Feature Histogram (PFH), which calculates the features of each point in the local neighborhood. The specific operation method is as follows:

[0168] Step 4.1: For each key point p i , use its normal vector n i to define the local coordinate system:

[0169] The axis u with the same direction as the normal vector:

[0170] u = n i (13);

[0171] Among them, n i is a certain point pi Normal vector;

[0172] Calculate the neighborhood point p j to p i vector d and orthogonalize it, i.e., axis v;

[0173] d = p j - p i (14);

[0174] where d is the vector from the neighborhood point p j to the key point p i ;

[0175]

[0176] where d is the vector from the neighborhood point p j to the key point p i , n i is the normal vector of a certain point p i here the u - axis represents;

[0177] Calculate axis w:

[0178] w = u × v (16);

[0179] where the u - axis is the normal vector n i of a certain point p i , v is the vector from the orthonormalized neighborhood point p j to p i .

[0180] For each neighborhood point p j , calculate the following 3 rotation - invariant features:

[0181] Normal vector angle f1:

[0182]

[0183] where w j is obtained by the cross - product of the normal vector of a certain point p i and the vector from the orthonormalized neighborhood point p j to p i , as shown in formula (16), n j is the normal vector of the neighborhood point p j , v j is calculated by formula (15);

[0184] Normal vector angle f2:

[0185] f2 = arccos(n i ·n j ) (18);

[0187] Among them, n i , n j is the normal vector of a certain point p i and its neighborhood point p j ;

[0188] The included angle f3 between the normal vector and the connection vector:

[0189]

[0190] Among them, n j is the normal vector of the neighborhood point p j , v j is the vector from the orthonormalized neighborhood point p j to p i , d is the vector from the neighborhood point p j to p i .

[0191] For each key point p i , calculate the (f1, f2, f3) distribution of all p j within its neighborhood to obtain the SPFH histogram.

[0192] Step 4.2: Improve the local description by weighting the SPFH of neighborhood points:

[0193]

[0194] Among them, Ν(p i ) is the neighborhood of point p i , w ij is the distance weight between the point pair (p i , p j ), and the formula is as follows:

[0195] w ij = ||D|| + ε (21);

[0196] Among them, D is the Euclidean distance between a certain point p i and its neighborhood point j, ε is a very small number to prevent a zero result, and the final result is the combination of the SPFH of each point itself and the weighted SPFH of neighborhood points, which improves the robustness and discrimination ability.

[0197] Step Five: Perform rough registration of the point cloud, and use Random Sample Consensus (RANSAC) to estimate the initial transformation matrix between the source point cloud and the target point cloud.

[0198] Based on the key points extracted from the above steps and the feature descriptions of the key points, corresponding point matching between the source point cloud and the target point cloud is performed. Since the overlap rate of the two point clouds is low, the random sample consensus (RANSAC) method is used to initially calculate the transformation matrix to align the source point cloud with the target point cloud as much as possible, achieving rough registration of the two point clouds; calculate the point cloud features, establish the feature matching relationship between the source point cloud and the target point cloud, randomly sample 3 feature points from the source point cloud, find the corresponding target points, and calculate the rigid transformation matrix of rotation R and translation t to align the 3 sampled points with the target point cloud as much as possible; then calculate the number of inliers, that is, calculate whether the transformed points fall within the reasonable range of the target point cloud, and the threshold of the target point cloud can be customized, and count the number of inliers; repeat the above operation steps 50,000 times, and select the transformation with the largest number of inliers as the final transformation result.

[0199] The rigid transformation matrix needs to be calculated based on the rigid transformation model, and the rigid transformation model is:

[0200] p i ' = Rp i + t (22);

[0201] where p i is the point before transformation in the source point cloud, p' i is the point after transformation in the source point cloud, R ∈ R 3×3 rotation matrix, t ∈ R 3×1 translation vector, provided that R and t are calculated.

[0202] Let the point p i in the source point cloud and the point q i in the target point cloud, calculate the covariance matrix:

[0203]

[0204] where p i is a certain point in the source point cloud, q i is a certain point in the target point cloud, and respectively represent the centroid, and M is the number of points in the source point cloud.,

[0205] Calculate the singular value decomposition (SVD):

[0206] UΣV T = SVD(H) (24);

[0207] where U, Σ, V are all 3×3 matrices, and the elements on the diagonal of the obtained matrix Σ are the singular values of the matrix H, and T is the transpose of the matrix.

[0208] Calculate the rotation matrix and translation respectively:

[0209] R = VU T (25);

[0210] Wherein, R is the required rotation matrix, U and V are the 3×3 matrices obtained above, and T represents the transpose of the matrix.

[0211]

[0212] Wherein, R is the required rotation matrix, U and V are the 3×3 matrices obtained above, t is the translation vector, and respectively represent the centroids of the source point p i and the target point q i of.

[0213] The final transformation matrix is:

[0214] Wherein, T r is the required transformation matrix.

[0215] Step 6: Perform precise point cloud registration, use TrICP for the final optimization work, remove incorrect point pairs, and use the mean matching distance MMD as a metric to calculate the deviation between two point clouds for iterative optimization; the operation method includes the steps:

[0216] Use the mean matching distance MMD as the basis for optimization convergence, iteratively optimize TrICP (Pruned Iterative Closest Point algorithm) for fine registration, achieve optimized alignment based on the result of rough registration, and solve the transformation matrix by minimizing the error:

[0217]

[0218] Wherein, R ∈ R 3×3 is the rotation matrix, t ∈ R 3×1 is the translation vector, p' i is the point obtained after rough registration of the source point cloud, q i is the point of the target point cloud, M is the number of points in the source point cloud; abnormal matching points are removed in each iteration, and only 50% of the best matching points are retained to prevent incorrect matching from affecting the result; the mean matching distance MMD is calculated after each iteration to monitor whether the error converges:

[0219]

[0220] Wherein, M is the number of points in the source point cloud, p' i is the point after transformation of the source point cloud, q i is the point in the target point cloud closest to p' i , ||·|| represents the Euclidean distance;

[0221] |MMD prev -MMD current |<threshold (30);

[0222] where MMD current is the calculated distance of the current iteration, and MMD prev is the result of the previous iteration. threshold is a user-defined threshold. If MMD converges, the iteration stops.

[0223] Figure 5 Shown are the source point cloud in blue, the target point cloud in green, the extracted key points in yellow, and the source point cloud after registration movement in red. It can be seen that the two point clouds overlap.

Claims

1. A method for urban road point cloud registration based on salient key point extraction, characterized by: Includes steps: Step 1: Take the acquired urban road point cloud as the processing object, remove the ground operation on the data, use a custom sparse processing method to downsample, and filter the dense area according to the number of neighboring points of the point cloud; Step 2: Calculate the normal vector of the point cloud using the improved adaptive normal. The implementation process of the improved adaptive normal calculation is as follows: Step 2.1: Use the adaptive normal vector calculation function to dynamically adjust the neighborhood radius to calculate the normal vector of the point cloud. Set the initial search radius base_radius to 0.

05. If the number of neighborhood points is less than 10, gradually increase the search radius by the set multiple each time, and up to 5 times; and add a minimum neighbor point constraint. If the number of neighborhood points is less than 5, assign the default normal vector [0,0,1] to prevent calculation anomalies; Step 2.2: For the points that meet the neighbor point requirements, the covariance matrix is ​​calculated by combining the principal component analysis PCA method, with the source point cloud P S For example, P S ={p1,p2,…,p n }, each point p i The coordinates of (x i ,y i ,z i ), calculate the centroid of each point’s neighborhood using the following formula: in is the calculated centroid, n is the number of neighborhood points, x i ,y i 、z i are the x-axis, y-axis, and z-axis coordinates of each point respectively; i is the i-th point in the neighborhood; Calculate the covariance matrix C, the formula is as follows: Among them, σ xx Represents a point p i The covariance of the x-axis coordinates of a neighborhood point, σ xy Represents a point p i The covariance of the x-axis coordinates and y-axis coordinates of the neighborhood points, σ xz Represents a point p i The covariance of the x-axis coordinates and z-axis coordinates of the neighborhood points, σ yx Represents a point p i The covariance between the y-axis coordinates and the x-axis coordinates of the neighborhood points, σ yy Represents a point p i The covariance between the y-axis coordinates of the neighborhood points and the y-axis coordinates, σ yz Represents a point p i The covariance of the y-axis coordinates and the z-axis coordinates of the neighborhood points, σ zx Represents a point p i The covariance of the z-axis coordinates and x-axis coordinates of the neighborhood points, σ zy Represents a point p i The covariance of the z-axis coordinates and y-axis coordinates of the neighborhood points, σ zz Represents a point p i The covariance between the z-axis coordinates of the neighborhood points and the z-axis coordinates; Each element is calculated as: Where a and b represent the coordinate axes (x, y, z). is the centroid coordinate of the corresponding axis, n is the number of neighborhood points; The eigenvalues ​​λ1, λ2, λ3 and the corresponding eigenvectors are calculated by eigenvalue decomposition, where the eigenvector corresponding to the minimum eigenvalue is the normal vector. The formula is as follows: Cv=λv (4); Where C is the covariance matrix solved in the previous step, v is the eigenvector corresponding to each λ and is also the required normal vector, and λ is the eigenvalue; Step 2.3: Adjust the search radius again based on the curvature information to optimize the normal vector calculation. The curvature calculation formula is as follows: curvature=λ min / (l sum +∈) (5); Where curvature is the curvature, λ min is the minimum eigenvalue obtained from the covariance matrix, λ sum is the sum of the three eigenvalues, ∈ is a very small number to prevent division by zero; The radius adjustment formula is: radius=base_radius×(1+curvature×b) (6); Where radius is the radius to be sought, base_radius is the radius initially set, curvature is the curvature, and b in the above formula is the magnification factor, which is set to 10. When the detected curvature is large, the radius is increased to capture more neighborhood points and make the normal calculation more accurate. Step 2.4: Make the calculated normal vector direction consistent, calculate the vector from each point to the centroid, and if the angle between the normal vector of the point and the vector from the point to the centroid is greater than 90°, flip the normal vector. This step ensures the consistency of the normal vector direction; Step 3: Use the salient corner point extraction algorithm to extract key points from the two point cloud data respectively; Step 4: Use the fast point feature histogram (FPFH) method to calculate the features of the key point neighborhood, i.e., the feature descriptor; Step 5: Perform rough registration of point clouds and use random sampling consistency (RANSAC) to estimate the initial transformation matrix between the source point cloud and the target point cloud; Step 6: Perform precise point cloud registration, use TrICP for final optimization, remove erroneous point pairs, and use the average matching distance MMD as a metric to calculate the deviation between the two point clouds for iterative optimization.

2. The urban road point cloud registration method based on salient key point extraction according to claim 1 is characterized by: The specific operation method of the step 1 is as follows: since the density of trees is relatively high, while the point cloud density of road signs, road signs or other signs is relatively uniform, a K-dimensional tree is used to perform a neighborhood search to calculate the number of neighbors of each point, and the number of neighbors within a radius of approximately 0.35 of each point is calculated. If the number of neighbor points is less than a threshold, the point is regarded as a sparse point, otherwise it is regarded as a dense point, and the dense points are removed and cyclic filtering is performed until the number of dense points is lower, and filtering is stopped.

3. The urban road point cloud registration method based on salient key point extraction according to claim 1 is characterized by: The specific operation method of extracting key points by the salient corner point extraction algorithm described in step 3 is: Step 3.1: Input two adjacent point clouds and the corresponding normals calculated in step 2, and set the neighborhood search radius, non-maximum suppression, and threshold detection parameters; Step 3.2: By calculating formula (2), (3), (4), the eigenvalues ​​of the covariance matrix of the calculation points at different scales s are obtained as part of the key point response formula, and the obtained eigenvalue response is normalized. The formula is as follows: Where R eig is the required eigenvalue response, λ1,λ2,λ3 are the eigenvalues ​​obtained from the covariance matrix, and ∈ is a very small number to prevent division by zero errors; Step 3.3: Calculate the change in the angle between the normal lines of the neighboring points in the point cloud at scale s as another part of the key point response formula to obtain the normal vector of a point i: n i =(n x ,n y ,n z ), the neighboring points of point i are: p j ∈Ν(p i ), calculate the neighborhood point p j Normal vector: n h =(n' x ,n' y ,n' z ), calculate the angle change, the formula is as follows: In the above formula, a point p i The normal vector n i , for point p i The neighboring point is p j , calculate the neighborhood point p j The normal vector n j ,θ ij is a certain point p i With neighboring point p j The angle, n i and n j is the normal vector corresponding to the two points, n i ·n j is the dot product, ||·|| is the norm, and the above formula is normalized to get the normal vector response: In the above formula, R angle is the required normal vector response, Ν(p i ) is a point p i The total number of neighborhood points, point p i The neighboring point is p j , 1-cosθ ij Reflects the magnitude of the change in the normal vector; Step 3.4: Add scale-adaptive weights and first calculate the average distance of neighborhood points: In the above formula, d avg,s is the average distance of neighborhood points at a certain scale s, Ν s (p i ) is a point p at a certain scale s i The total number of neighboring points is p j , D is a point p i and its neighboring point p j The Euclidean distance of Then calculate the scale-adaptive weighting factor: Among them, w s is the scale factor, γ is a hyperparameter that controls scale adaptivity, and d avg,s is the average distance of neighborhood points at a certain scale s, and e is a constant; Step 3.5: The final key point response formula combines the results of steps 3.2, 3.3, and 3.4, calculates the response value at multiple scales, and takes the maximum value as follows: Among them, R i is the final response formula, s represents a certain scale, N(s) represents all scales, and w s is the scale weight of a certain scale, R eig,s is the eigenvalue response at a certain scale, R angle,s is the normal vector response at a certain scale, α is the weight for controlling the eigenvalue response, β is the weight for controlling the normal vector change response, the response is calculated at different radii s and the maximum value is taken to improve scale invariance, and w is introduced s Weighting is performed to make the response more robust to changes in point cloud density; the key point response formula is calculated for each of the two point clouds, and if the result is greater than a threshold, it is classified as a key point.

4. The urban road point cloud registration method based on salient key point extraction according to claim 1 is characterized by: The calculation of FPFH in step 4 is based on the point feature histogram PFH, which calculates the features of each point in the local neighborhood; the specific operation method in step 4 is: Step 4.1: For each key point p i , using its normal vector n i Define the local coordinate system: Same axis u as the normal vector: u=n i (13); Among them, n i is a certain point p i The normal vector of Calculate the neighborhood point p j to p i The vector d is orthogonalized, that is, the axis v; d=p j -p i (14); Where d is the neighborhood point p j To the key point p i A vector of Where d is the neighborhood point p j To the key point p i The vector of n i is a certain point p i The normal vector of , represented by the u axis here, has been transformed in formula (13); Calculate axis w: w = u × v (16); Among them, the u axis represents a point p i The normal vector n i , v is the orthogonalized neighborhood point p j to p i A vector of For each neighborhood point p j , calculate the following three rotation invariant features: Normal vector angle f1: Among them, w j is a point p i The normal vector of the orthogonalized neighborhood point p j to p i The cross product of the vector is obtained, as shown in formula (16), n j is the neighborhood point p j The normal vector, v j Calculated by formula (15); Normal vector angle f2: f2=arccos(n i ·n j ) (18); Among them, n i , n j is a certain point p i and the neighborhood point p j The normal vector of The angle f3 between the normal vector and the connection vector: Among them, n j is the neighborhood point p j The normal vector, v j is the orthogonalized neighborhood point p j to p i vector, d is the neighborhood point p j to p i A vector of For each key point p i , calculate all p in its neighborhood j The (f1, f2, f3) distribution is obtained, and the SPFH histogram is obtained; Step 4.2: Improve the local description by weighting the SPFH of the neighborhood points: Among them, Ν(p i ) is point p i Neighborhood, w ij is a point pair (p i ,p j ), the distance weight between them is as follows: w ij =||D||+ε (21); Where D is a point p i The Euclidean distance to its neighboring point j, ε is a very small number to prevent the result of zero. The final result is the combination of the SPFH of each point itself and the weighted SPFH accumulation of the neighboring points, which improves the robustness and discrimination ability.

5. The urban road point cloud registration method based on salient key point extraction according to claim 1 is characterized by: The specific operation method of step 5 is as follows: based on the key points extracted in the above steps and the feature description of the key points, the corresponding points of the source point cloud and the target point cloud are matched. Since the overlap rate of the two point clouds is low, the random sampling consistency RANSAC method is used to preliminarily calculate the transformation matrix so that the source point cloud is aligned with the target point cloud as much as possible, so as to achieve a rough registration of the two point clouds; Calculate the point cloud features, establish the feature matching relationship between the source point cloud and the target point cloud, randomly sample 3 feature points from the source point cloud, find the corresponding target point, rotate R and translate t to calculate the rigid transformation matrix, so that the sampled 3 points are aligned to the target point cloud as much as possible; Then calculate the number of inliers, that is, calculate whether the transformed points fall within the reasonable range of the target point cloud. The threshold of the target point cloud can be customized, and count the number of inliers. Repeat the above steps 50,000 times, and select the transformation with the largest number of inliers as the final transformation result.

6. The urban road point cloud registration method based on salient key point extraction according to claim 5 is characterized by: The rigid transformation matrix needs to be calculated based on the rigid transformation model, and the rigid transformation model is: p i '=Rp i +t (22); Among them, p i is the point in the source point cloud before transformation, p' i is the transformed point in the source point cloud, R∈R 3×3 The rotation matrix, t∈R 3×1 The translation vector is obtained after R and t are calculated. Let point p in the source point cloud be i and point q in the target point cloud i , calculate the covariance matrix: where p i is a point in the source point cloud, q i is a point in the target point cloud, and They represent the center of mass respectively; Compute the singular value decomposition (SVD): UΣV T =SVD(H) (24); Where U, Σ, V are all 3×3 matrices, the elements on the diagonal of the matrix Σ are the singular values ​​of the matrix H, and T is the transpose of the matrix; Calculate the rotation matrix and translation separately: R=VU T (25); Where R is the required rotation matrix, U, V are the 3×3 matrices obtained above, and T is the transpose of the matrix; Among them, R is the required rotation matrix, U, V are the 3×3 matrices obtained above, and t is the translation vector. and Represent the source point p i and the target point q i The centroid of The final transformation matrix is: Among them, T r This is the required transformation matrix.

7. The urban road point cloud registration method based on salient key point extraction according to claim 6 is characterized by: The operation method for performing accurate point cloud registration described in step 6 includes the following steps: Use the average matching distance MMD as the basis for optimizing convergence, iteratively optimize TrICP, prune the iterative closest point algorithm, perform fine registration, optimize alignment based on the results of coarse registration, and solve the transformation matrix by minimizing the error: Where R∈R 3×3 The rotation matrix, t∈R 3×1 The translation vector, p' i is the point obtained after the rough registration of the source point cloud, q i is a point in the target point cloud, and M is the number of points in the source point cloud. In each iteration, abnormal matching points are removed and only 50% of the best matching points are retained to prevent incorrect matching from affecting the results. After each iteration, the average matching distance MMD is calculated to monitor whether the error converges: Where M is the number of points in the source point cloud, p' i is a point in the source point cloud after transformation, q i is the target point cloud with p' i The nearest point, ||·|| represents the Euclidean distance; |MMD prev -MMD current |<threshold (30); Among them, MMD current is the calculated distance of the current iteration, MMD prev is the result of the previous iteration, threshold is a custom threshold, and the iteration stops if MMD converges.

Images