Method for extracting local surface features of point cloud based on morphology

Abstract

The application discloses a kind of based on morphological extraction point cloud local surface feature method, steps are as follows: S1 acquires point cloud data and carries out data pretreatment to point cloud data;S2 searches the neighborhood point of each point of point cloud in three-dimensional space, calculates normal information, then based on normal information, the curvature of each point is calculated, and curvature value is obtained;Then refer to the curvature value range corresponding to local surface feature and filter out target surface;S3 uses progressive morphological filter to adjust filter parameter constantly, according to the curvature change of selected local area, the size of filter window is adjusted accordingly, non-target point in point cloud data after curvature screening is filtered out, and the point of target local area is extracted;S4 uses dilatation, corrosion three-dimensional morphological operator and dilatation corrosion combination open operation and close operation to optimize target local area, and extract the local surface feature area of complete smooth point cloud data.The method is extracted to local area with obvious curvature change, high accuracy, and wide adaptability.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of point cloud data processing, and in particular to a method for extracting local surface features of point cloud based on morphology. BACKGROUND

[0002] At present, point cloud filtering algorithms mainly focus on LiDAR data research, especially for terrain regions with obvious features, such as large buildings in urban clusters, roads, forest vegetation, and mountain ridges and valleys, aiming to remove noise and non-ground information in point cloud data. However, when it is necessary to extract local regions with similar terrain features from point cloud data, the use of point cloud filtering for processing shows poor performance, and it is difficult to accurately extract and separately view each local feature.

[0003] Morphological filtering uses the opening operation of erosion followed by dilation as the basis for processing LiDAR data, and can remove discrete points and smooth boundaries without significantly changing the overall position and shape, and can filter out objects with elevation mutations within a specified window size. Morphological filtering requires repeated parameter setting during the filtering process, and it is difficult to select an appropriate filter window size. If the filter window size is too small, it is easy to retain objects with large sizes, and if the filter window size is too large, it will incorrectly filter out information similar to mountain tops, and the set parameters have poor universality for different terrain features. The introduction of a height difference threshold for point cloud judgment can retain points smaller than the height difference threshold, but a slope parameter needs to be preset, and the height difference threshold cannot automatically adapt to terrain changes.

[0004] Therefore, the present application proposes a method for extracting local surface features of point cloud based on morphology, which is suitable for extracting local regions with obvious curvature changes, has high accuracy and wide adaptability. SUMMARY

[0005] The technical problem to be solved by the present application is to provide a method for extracting local surface features of point cloud based on morphology, which uses three-dimensional morphological operations to repair information loss and noise points caused during the processing process, and optimizes the target feature region to make it complete and smooth. The method for extracting local surface features of point cloud provided by the present application is suitable for extracting local regions with obvious curvature changes, has high accuracy and wide adaptability.

[0006] To solve the above technical problems, the technical scheme adopted by the present application is as follows: the method for extracting local surface features of point cloud based on morphology, specifically comprising the following steps:

[0007] S1: collecting point cloud data and performing data preprocessing on the point cloud data;

[0008] S2: searching neighborhood points of each point of the point cloud in a three-dimensional space, calculating normal information, and then calculating curvature of each point based on the normal information to obtain a curvature value; and filtering out a target surface according to a curvature value range corresponding to a local surface feature;

[0009] S3: continuously adjusting a filtering parameter by using a progressive morphological filter, adjusting a size of a filtering window according to a curvature change of a selected local region, filtering out non-target points in the point cloud data after the curvature screening, and extracting points of the target local region;

[0010] S4: optimizing the target local region by using an inflation, erosion three-dimensional morphological operator and an open operation and a close operation of the inflation-erosion combination to extract a complete and smooth local surface feature region of the point cloud data.

[0011] By using the technical scheme, the obtained point cloud data is preprocessed to obtain clear point cloud data for subsequent processing; points of a target region of the point cloud are filtered according to a curvature range corresponding to a local surface of different terrain types, so that the local feature extraction in the next step is facilitated; the target feature region is extracted by using the progressive morphological filter, the size of the window is adjusted according to the curvature change in the local region selected by the filtering window, the filtered slope value is recorded to maintain the dynamic update of the height difference threshold, and the non-ground points are more accurately filtered out; the progressive morphological filter is applied to extract the local region with different curvature changes in the point cloud by using the curvature change characteristics of each point on the point cloud, and the three-dimensional morphological operator is used to fill in the missing part of the local feature region and filter out the point cloud noise, so that the smooth and complete local feature region is obtained; and the method for extracting the local surface feature of the point cloud provided by the application extracts the local region with obvious curvature change, has high accuracy and wide adaptability.

[0012] Preferably, the data preprocessing of the point cloud data in step S1 specifically includes the following steps:

[0013] S11: calculating distances between points in the point cloud data, screening and deleting repeated points in the point cloud data according to the minimum distance between the points; dividing points with too close distances from two points and one point into repeated points and deleting the repeated points; and preventing the repeated points from interfering with subsequent processing operations;

[0014] S12: dividing the point cloud space by using a point cloud downsampling method to obtain a plurality of spaces; and selecting a point in each space after the division to replace all points in the space;

[0015] S13: Remove sparse outliers existing in the down-sampled point cloud data by point cloud filtering, and then obtain the main points by searching the neighborhood of each point in the point cloud. The point cloud filtering removes the sparse outliers after down-sampling to avoid affecting the subsequent search and identification of points; the preferred method is to use statistical filtering to remove outliers or gross error points caused by measurement errors on the down-sampled point cloud data. This filtering method does not need to set the filtering parameters in advance, avoiding increasing the difficulty of target classification and affecting the subsequent search and identification of points.

[0016] Preferably, in the step S12, specifically, the point cloud data is uniformly sampled on the curved surface, the point cloud space is divided to obtain a plurality of spaces, and then in each space after division, a point is selected to replace all points within a sphere with the point as the center and R as the radius. The uniform sampling on the curved surface replaces all points within the sphere with the point as the center with one point, and the sampled points are uniformly distributed, but the position of the point cloud is not changed, and good accuracy is achieved. The implementation of the algorithm depends on the Open3d library, and the larger the radius of the selected sphere, the fewer the points after sampling. The sampled points are uniformly distributed, but the position and basic outline of the point cloud are not changed, and good accuracy is achieved.

[0017] Preferably, after the statistical filtering of the point cloud data in the step S13, the spatial topological relationship is established, the neighborhood point set of the target point is obtained by k-neighbor searching, and if the neighborhood points of the target point satisfy the Gaussian distribution, the average distance between the target point and the k neighborhood points is calculated, the standard range of the average distance is determined by the mean and standard deviation of the global distance, the target point with an average distance exceeding the standard range from the neighborhood points is defined as an outlier and is removed, and if the average distance between the target point and the k neighborhood points does not exceed the standard range, the target point is retained as the main point of the model.

[0018] Preferably, the specific steps of the step S2 are as follows:

[0019] S21: Obtain the normal equation based on the least square fitting plane, and obtain the normal value of the plane represented by the eigenvector through the normal equation; it is determined that the normal of a point on the plane is similar to the normal of a tangent plane of the estimated plane, so after conversion, it becomes a least square plane fitting estimation problem;

[0020] S22: Based on the obtained point cloud data and the calculated normal value of each point in the point cloud data, the curvature information of the target point is approximated using the surface variation in the k-neighborhood of the target point, and the geodesic curvature of each point in the point cloud data is calculated.

[0021] S23: Set the curvature interval by combining the curvature value range corresponding to different types of local surfaces, and screen out the target local surface feature.

[0022] Preferably, the KdTree nearest neighbor algorithm is used in the step S21, the KdTree is a data structure based on binary space partitioning, the entire k-dimensional vector space is divided into several regions, the position of each local space is determined by continuously performing branch judgment in the search process; the normal vector estimation NormalEstimation in the point cloud library (PCL library) is used to calculate the normal feature of the surface where the target point is located, specifically: the plane expression is converted from the general formula to the normal formula equation, the optimal solution of the equation is solved by introducing the Lagrange multiplier method, the partial derivatives of the parameters in the expression are solved and the obtained expressions are unified, the eigenvalues are obtained by performing eigenvalue decomposition on the covariance matrix in the expression, and the eigenvector corresponding to the minimum eigenvalue in each eigenvalue is the normal vector of the fitted plane. The KDTree is an index tree structure evolved from the binary search tree for large-scale high-dimensional data search scenarios, and is mainly used in nearest neighbor search and approximate nearest neighbor search scenarios. The position of each local space can be determined by the KdTree nearest neighbor algorithm, which can avoid searching the global space; the nearest neighbor algorithm based on KdTree starts from the root node and recursively searches the binary tree downward until the leaf node is found, and then the search path is backtracked, the KdTree is constructed by iteratively backtracking the leaf node to the root node, and then the KdTree is used to search the point cloud data, to check whether there is a closer node in the other subspace, to search the child node area and backtrack to the root node, which can not only save the calculation time, but also can ensure the matching and searching of the point cloud associated points at any time. The neighborhood points of the target point are searched based on the KdTree, a direction n is sought by PCA analysis so that the distribution of the projection points of all neighborhood points in the direction n is the most concentrated, that is, the variance of the projection of the target point in the direction is the smallest, and the normal corresponds to the least important direction in the principal component. The least square rule is used to fit the plane of the target point and its neighborhood points; therefore, in order to obtain the best fitted plane, the sum of the square distances of any point to the plane should be minimized, therefore the plane expression is converted from the general formula to the normal formula equation, the problem is converted to solving the extreme value of the normal equation, the optimal solution of the equation is solved by introducing the Lagrange multiplier method, the partial derivatives of the parameters in the expression are solved and the obtained expressions are unified, the eigenvalues are obtained by performing eigenvalue decomposition on the covariance matrix in the expression, and the eigenvector corresponding to the minimum eigenvalue in each eigenvalue is the normal vector of the fitted plane; from solving the normal to solving the eigenvector corresponding to the minimum eigenvalue, the normal value of each point of the point cloud is represented; finally, the normal direction is checked, if the position of the viewpoint of the point cloud is known, all the normals are consistent with the viewpoint direction.

[0023] Preferably, in the step S22, a point cloud curvature estimation method (Principal Curvatures Estimation) in a point cloud library (PCL library) is adopted to search the neighborhood points of each point in the point cloud based on a KdTree, and the principal curvatures of each point in the point cloud are calculated according to the following formula:

[0024]

[0025] wherein σ represents the surface variation of the target point in the k neighborhood; λ0represents the surface variation along the normal vector; λ1and λ2respectively represent the distribution of the target point on the tangent plane; the surface variation of the target point in the k neighborhood is used to approximately represent the curvature information, and the concave-convex condition of the surface in the local region is determined according to the sign of the curvature value. The curvature is a second-order differential quantity of the surface in differential geometry, which is represented by a differential form to represent the degree of deviation of the curve from a straight line. The greater the value, the higher the degree of deviation. The curvature of the object reflects the concave-convex degree of the surface of the object. The normal curvature represents the bending degree of the surface at a vertex in different directions. The principal curvatures are the maximum k1and minimum k2of the normal curvature; the mean curvature H k is the average value of the two orthogonal curvatures at a vertex, and the Gaussian curvature H is the product of the two principal curvatures, reflecting the overall bending degree of the vertex.

[0026] Preferably, the specific steps of the step S3 are as follows:

[0027] S31: using progressive morphological filtering to distinguish the ground points and non-ground points in the point cloud;

[0028] S32: using a three-dimensional morphological dilation operator d and a three-dimensional morphological erosion operator e to filter out the non-ground points marked by the filtering window on the target local surface, so as to extract the ground points of the target local surface.

[0029] Preferably, the specific steps of the step S31 are as follows:

[0030] S311: increasing the filtering window and performing iterative operation by linearly increasing the window size, and then adjusting the size of the filtering window according to the curvature change of the local region selected by the filtering window;

[0031] S312: setting the slope value in advance according to the average slope of the point cloud data, and extracting the target surface based on the slope change;

[0032] S313: adjusting the height difference threshold value by using the slope value set in advance in step S312, and adjusting the height difference threshold value adaptively according to the change of the terrain slope and the window size, and determining the points with the height difference value greater than the set height difference threshold value as the non-ground points; and recording the slope information of the filtered curved surface, and calculating the height difference threshold value based on the slope information to keep the dynamic update of the height difference threshold value. When the size of the filtering window is adjusted according to the change of the curvature of the local area selected by the filtering window, if the curvature change of each point on the local area is approximately zero, the curved surface is flat, and the size of the filtering window can be further enlarged, and if the curvature change of the local area corresponds to the convex or concave area characteristics, the increase speed of the filtering window is adjusted.

[0033] Preferably, the formula for increasing the filtering window size when performing the iterative operation in step S311 with the linearly increased window size is as follows:

[0034] w k = nk+2;

[0035] wherein w k is the window size of the kth filtering; k is the iteration number, k=0, 1, 2, …m; and n is the initial window size. n is generally 2. The window size parameter determines the detail degree of the filtering, but it is time-consuming. Therefore, the iteration number is reduced to save time by performing the iterative operation with the linearly increased window size.

[0036] Preferably, in step S313, the slope value s is assumed to be constant, and the current height difference threshold value, the window size and the terrain slope value s satisfy the following formula:

[0037]

[0038] wherein h max(t)，k represents the maximum height difference threshold value of the kth iteration at t time, w k is the window size of the kth filtering, w k-1 is the window size of the (k-1)th filtering; h k represents the height difference threshold value under the current filtering window; h0 represents the initial height difference threshold value; c is the grid unit; and h max represents the maximum height difference threshold value.

[0039] Preferably, the formula of the three-dimensional morphological dilation operator d in step S32 is as follows:

[0040]

[0041] The calculation formula of the erosion operator e is as follows:

[0042]

[0043] Wherein, w is a structure element window; z(x, y) represents the elevation value of the xth row and yth column of the generated image after calculation.

[0044] Preferably, in the step S32, if the position of the non-target point cannot be filtered out by the filter window, the rotate function in the Open3d library is used to set the rotation center and the rotation matrix to rotate the point cloud data, and the required filtered point is rotated above the target surface; If a rotation of θ degrees around the x-axis is selected, the x-coordinate remains unchanged, and a two-dimensional rotation is performed on the yOz plane composed of y and z, then the formula is:

[0045] x' = x

[0046] y' = y cosθ-z sinθ

[0047] z' = y sinθ+z cosθ.

[0048] Further, the rotation matrix is obtained;

[0049]

[0050] Preferably, the step S4 is specifically: using a three-dimensional dilation operator d to connect the discontinuous part of the target surface, using a three-dimensional erosion operator e to separate the adhered objects and shrink the region, using an opening operation to smooth the boundary of the target region, using a closing operation to connect the boundary of the object, filling the disconnected connected region, and eliminating the boundary points and noise points.

[0051] Wherein, the dilation is a morphological operation for local maximum value, and the effect is to reconnect the discontinuous part of the surface and expand the target region; the definition of A being dilated by S is: all points in the set Z are displaced by the points in the structure element S, and then the union is calculated to obtain the dilated image, and the calculation formula is:

[0052]

[0053] Wherein, A is a set of subspaces Z 2 , S is a structure element of the subspace Z 2 , S being dilated by S is a set of all displacements, is a dilation operator;

[0054] The erosion operation is a morphological operation for local minimum value, and the main function is to separate the adhered objects and shrink the boundary line inward; using S to erode A is to translate the set of points contained in A in all structure elements S by Z, and the erosion expression is:

[0055]

[0056] Wherein, is an erosion operator.

[0057] The open operation is to erode the image first and then dilate, which can separate the adhered regions and also can smooth the boundary of the image without changing the original area, and the expression is:

[0058]

[0059] Wherein, The open operator is represented by;

[0060] The close operation is to dilate the image first and then erode, which can connect the boundary between two adjacent objects and fill the disconnected connected regions, and the expression is:

[0061]

[0062] Wherein, the close operator is represented by;

[0063] By adopting the technical scheme, the point cloud data is easy to cause local missing, hollow and fine seam of the curved surface, and the originally non-adhered region appears adhesion after processing, at this time, the three-dimensional morphological operator is used to fill and smooth the local feature region of the extracted point cloud with defects, and the noise points are eliminated; the three-dimensional dilation operator is used to connect the disconnected part of the target curved surface, the three-dimensional erosion operator is used to separate the adhered objects and shrink the region, the open operation is used to smooth the boundary of the target region, the close operation is used to connect the boundary of the object, fill the disconnected connected region, and eliminate the boundary points and noise points, and the relatively complete and smooth point cloud local curved surface feature region is extracted.

[0064] Compared with the prior art, the present application has the beneficial effects that:

[0065] (1) The current point cloud filtering algorithm is mainly used for filtering the terrain regions with obvious features, such as large buildings in urban areas, roads, forest vegetation, ridges, valleys, etc., the present application uses the characteristics of the curvature change of each point on the point cloud to apply the progressive morphological filter to extract the local regions with different curvature changes in the point cloud, which is convenient for subsequent research on each local region; it can also be used to filter most of the point clouds with obvious features first, to prevent interference when processing the remaining small and complex local regions later;

[0066] (2) Since the data points on the point cloud are very large, preprocessing is a very necessary process; removing duplicate points can reduce the data amount and save calculation time for subsequent operations, uniform sampling of the curved surface filters out redundant points so that the remaining points are uniformly distributed and do not change position, statistical filtering removes obvious outliers to avoid interference with subsequent operations, and the preprocessed point cloud data is small and uniformly distributed but has no noise points;

[0067] (3) A difficulty in the processing of the progressive morphological filter is the selection of the filter window size, and a too large window is likely to filter out part of the feature information, and a too small window is likely to retain the large-sized object as a ground point; therefore, the filter window is gradually linearly increased, and the size of the filter window is adjusted according to the curvature variation of each point in the local area selected by the filter window, so that the size of the filter window can be more accurately set, and thus the non-ground point is filtered out;

[0068] (4) The obtained point cloud data is relatively rough, and there are noise points and incomplete areas of collection, and in addition, the area is missing in the filtering process, so the extracted target feature area usually needs to be optimized; the three-dimensional morphological operator can smooth the boundary, fill the missing area, separate the adhered area, dilate and shrink the local area, so that the finally extracted target feature area is relatively complete. BRIEF DESCRIPTION OF DRAWINGS

[0069] Figure 1 is a flowchart of the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0070] Figure 2 is a three-dimensional point cloud graph of a tooth original of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0071] Figure 3 is a three-dimensional tooth point cloud graph after preprocessing of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0072] Figure 4 is a rough tooth segmentation line graph after curvature screening of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0073] Figure 5 is a tooth segmentation line graph after progressive morphological filtering of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0074] Figure 6 is a tooth segmentation line graph after optimization of a three-dimensional morphological operator of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application;

[0075] Figure 7 is a three-dimensional point cloud erosion graph of a specific embodiment in the method for extracting a local surface feature of a point cloud based on morphology of the present application. DETAILED DESCRIPTION

[0076] The technical solutions in the embodiments of the present application will be described clearly and completely in combination with the drawings in the embodiment diagram of the present application.

[0077] Embodiment: As shown in the figure, the method for extracting local surface features of point cloud based on morphology specifically includes the following steps: Figure 1

[0078] S1: Collecting point cloud data and performing data preprocessing on the point cloud data;

[0079] The data preprocessing on the point cloud data in step S1 specifically includes the following steps:

[0080] S11: Calculate the distance between points in the point cloud data, and filter and delete the repeated points in the point cloud data according to the minimum distance between points; points with a distance too close to two points are divided into repeated points and deleted; to prevent repeated points from causing interference in subsequent processing operations;

[0081] S12: Divide the point cloud space by using a point cloud downsampling method to obtain a plurality of spaces; then select a point in each space after division to replace all points in the space;

[0082] In step S12, specifically: uniformly sampling the surface of the point cloud, dividing the point cloud space to obtain a plurality of spaces; then in each space after division, a point is selected to replace all points within the sphere with the point as the center and R as the radius. Uniformly sampling the surface of the point cloud replaces all points within the sphere with the point as the center with the point. The sampled points are uniformly distributed, but do not change the position of the point cloud and have good accuracy; the implementation of this algorithm depends on the Open3d library. The larger the radius of the selected sphere, the fewer the points after sampling. The sampled points are uniformly distributed, but do not change the position of the point cloud and have good accuracy;

[0083] S13: Remove sparse outliers in the point cloud data after downsampling by point cloud filtering, and then obtain the main points by searching the neighborhood of each point in the point cloud. The point cloud filtering removes sparse outliers after downsampling to avoid affecting the search and identification of subsequent points; the preferred method here is to use statistical filtering to remove sparse outliers after downsampling; statistical filtering is used to select and remove outliers or gross error points caused by measurement errors on the point cloud data after downsampling. This filtering method does not require the setting of filtering parameters in advance, avoiding the increase in the difficulty of target classification and the impact on the search and identification of subsequent points;

[0084] ​After the point cloud data is statistically filtered in the step S13, a spatial topological relation is established, a neighborhood point set of the target point is obtained through k-neighbor searching, if the neighborhood points of the target point satisfy a Gaussian distribution, an average distance between the target point and the k neighborhood points is calculated, a standard range of the average distance is determined by a mean value and a standard deviation of the global distance, the target point whose average distance with the neighborhood points exceeds the standard range is defined as an outlier and is removed, if the average distance between the target point and the k neighborhood points does not exceed the standard range, the target point is reserved as a main point of the model;

[0085] S2: searching neighborhood points of each point of the point cloud in the three-dimensional space, calculating normal information, then calculating curvature of each point based on the normal information to obtain a curvature value, and screening out a target surface according to a curvature value range corresponding to a local surface feature;

[0086] The specific steps of the step S2 are as follows:

[0087] S21: obtaining a normal equation based on a least square fitting plane, and obtaining a normal value of the plane through the normal equation; it is determined that the normal of a point of the plane is similar to a tangent plane normal of the estimated plane, so the problem is converted into a least square plane fitting estimation problem;

[0088] In the step S21, the KdTree nearest neighbor algorithm is used, the KdTree is a data structure based on binary space partitioning, the entire k-dimensional vector space is divided into several regions, the position of each local space is determined by continuously performing branch judgment in the search process; then the normal vector estimation Normal Estimation in the point cloud library (PCL library) is used to calculate the normal feature of the surface where the target point is located, specifically: the plane expression is converted from the general formula to the normal formula equation, the optimal solution of the equation is solved by introducing the Lagrange multiplier method, the partial derivatives of the parameters in the expression are calculated and the obtained expressions are unified, the covariance matrix in the expression is subjected to eigenvalue decomposition to obtain each eigenvalue, wherein the eigenvector corresponding to the minimum eigenvalue in each eigenvalue is the normal vector of the fitted plane, and the normal vector represents the normal value of the point; the KDTree is an index tree structure evolved from the binary search tree for large-scale high-dimensional data search scenarios, mainly used in nearest neighbor search and approximate nearest neighbor search scenarios, and the position of each local space can be determined by the KdTree nearest neighbor algorithm to avoid searching the global space; the nearest neighbor algorithm based on KdTree starts from the root node and recursively searches the binary tree downward, until the leaf node is found, and then the search path is backtracked, the leaf node is iterated and backtracked to the root node to complete the construction of the KdTree, and then the KdTree is used to search the point cloud data, to check whether there is a closer node in another subspace, to search the subspace and backtrack to the root node, so that the calculation time can be saved, and the matching and searching of the point cloud associated points can be ensured at any time. The neighborhood points of the target point are searched based on the KdTree, a direction n is sought by PCA analysis so that the distribution of the projection points of all neighborhood points in the direction n is the most concentrated, that is, the variance of the projection of the target point in the direction is the smallest, and the normal corresponds to the least important direction in the principal component; the least square rule is used to fit the plane of the target point and its neighborhood points; therefore, to obtain the best fitted plane, the sum of the square distances of any point to the plane should be minimized, therefore the plane expression is converted from the general formula to the normal formula equation, the problem is converted to solving the extreme value of the normal equation, the optimal solution of the equation is solved by introducing the Lagrange multiplier method, the partial derivatives of the parameters in the expression are calculated and the obtained expressions are unified, the covariance matrix in the expression is subjected to eigenvalue decomposition to obtain each eigenvalue, and the eigenvector corresponding to the minimum eigenvalue is the normal vector of the fitted plane; the normal is converted to the eigenvector corresponding to the minimum eigenvalue, which represents the normal value of each point of the point cloud; finally, the normal direction is checked, if the position of the viewpoint of the point cloud is known, all the normals are consistent with the viewpoint direction;

[0089] S22: Based on the obtained point cloud data and the normal value of each point in the calculated point cloud data, the curvature information of the target point in the k neighborhood is approximated using the surface variation, and the geodesic curvature of each point in the point cloud data is calculated;

[0090] The step S22 adopts a point cloud curvature estimation method (Principal Curvatures Estimation) in a point cloud library (PCL library) to search the neighborhood points of each point in the point cloud based on KdTree, calculate the principal curvatures of each point in the point cloud, and the calculation formula is as follows:

[0091]

[0092] wherein σ represents the surface variation of the target point in the k neighborhood; λ0 represents the surface variation along the normal vector; λ1 and λ2 respectively represent the distribution of the target point on the tangent plane; the surface variation of the target point in the k neighborhood is used to approximately represent the curvature information; and the concave-convex condition of the surface in the local region is determined according to the sign of the curvature value. The curvature is a second-order differential quantity of the surface in differential geometry, which is expressed by a differential form to represent the degree of deviation of the curve from a straight line. The greater the numerical value, the higher the degree of deviation. The curvature of the object reflects the concave-convex degree of the surface of the object. The normal curvature represents the bending degree of the surface at a vertex along different directions. The principal curvatures are the maximum value k1 and the minimum value k2 of the normal curvature; the mean curvature H k is the average value of the two orthogonal curvatures at a vertex, and the Gaussian curvature H is the product of the two principal curvatures, reflecting the overall bending degree of the vertex;

[0093] S23: Set the curvature interval in combination with the curvature value range corresponding to the local surface of different types, and screen out the target local surface feature;

[0094] S3: Adjust the filtering parameters by using the progressive morphological filtering, adjust the size of the filtering window according to the curvature change of the selected local region, filter out the non-target points in the point cloud data after the curvature screening, and extract the points of the target local region;

[0095] The specific steps of the step S3 are as follows:

[0096] S31: Distinguish the ground points and non-ground points in the point cloud by using the progressive morphological filtering;

[0097] The specific steps of the step S31 are as follows:

[0098] S311: Increase the filtering window and perform iterative operation by linearly increasing the window size, and then adjust the size of the filtering window according to the curvature change of the selected local region in the filtering window;

[0099] The formula for increasing the filtering window when performing iterative operation by linearly increasing the window size in the step S311 is as follows:

[0100] w k = nk+2;

[0101] wherein wk Wk is the window size of the kth filtering; k is the iteration number, k = 0, 1, 2, … m; n is the initial window size. n is generally 2. The parameter of window size determines the degree of filtering details, but it is time-consuming, so in order to reduce the iteration number to save time, the iteration operation is carried out by linearly increasing the window size.

[0102] In the step S313, it is assumed that the slope value s is constant, and the current height difference threshold value, the window size and the terrain slope value s have the following formula relationship:

[0103]

[0104] Wherein, h max(t)，k represents the maximum height difference threshold value of the kth iteration at t time, w k is the window size of the kth filtering, w k-1 is the window size of the k-1th filtering; h k represents the height difference threshold value under the current filtering window; h0 represents the initial height difference threshold value, c is the grid unit; h max represents the maximum height difference threshold value.

[0105] S312: The slope value is set in advance according to the average slope of the point cloud data, and the target surface is extracted based on the slope change;

[0106] S313: The height difference threshold value is adjusted by using the slope value set in advance in the step S312, and the height difference threshold value is adaptively adjusted according to the change of the terrain slope and the window size. The points with the height difference value greater than the set height difference threshold value are determined as non-ground points; and the slope information of the filtered surface is recorded, and the height difference threshold value is calculated based on the slope information, so as to keep the dynamic update of the height difference threshold value setting. When the size of the filtering window is adjusted with the aid of the curvature change of the local area selected by the filtering window, if the curvature change of each point on the local area is approximately zero, the surface is flat, and the size of the filtering window can be further enlarged, and if the curvature change of the local area corresponds to the convex or concave area characteristics, the increase speed of the filtering window is adjusted;

[0107] S32: The three-dimensional morphological dilation operator d and the three-dimensional morphological erosion operator e are used to filter out the non-ground points marked by the target local surface in the filtering window, so as to extract the ground points of the target local surface;

[0108] The formula of the three-dimensional morphological dilation operator d in the step S32 is:

[0109]

[0110] The calculation formula of the erosion operator e is:

[0111]

[0112] wherein w is a structure element window; z(x, y) represents the elevation value of the xth row and yth column of the generated image after calculation;

[0113] If the position of the non-target point cannot be filtered out using the filter window in the step S32, the rotate function in the Open3d library (Open3D is prior art, source: A Modern Library for 3D Data Processing-Open3D master (4c5ded8) documentation) is used to set the rotation center and rotation matrix to rotate the point cloud data, and the point to be filtered out is rotated above the target surface; wherein if the rotation is selected around the x-axis by θ angle, the x-coordinate remains unchanged, and two-dimensional rotation is performed on the yOz plane composed of y and z, then the formula is:

[0114] x' = x

[0115] y' = y cosθ - z sinθ

[0116] z' = y sinθ + z cosθ

[0117] and further obtain the rotation matrix;

[0118]

[0119] S4: The local region of the target is optimized by comprehensively using the dilation, erosion three-dimensional morphological operator and the opening operation and the closing operation of the dilation-erosion combination, and the local surface feature region of the complete and smooth point cloud data is extracted; the step S4 is specifically: the three-dimensional dilation operator d is used to connect the interruption of the target surface, the three-dimensional erosion operator e is used to separate the adhered objects and shrink the region, the opening operation is used to smooth the boundary of the target region, the closing operation is used to connect the boundary of the object, fill the disconnected connected region, and eliminate the boundary points and noise points;

[0120] wherein the dilation is a morphological operation for finding the local maximum value, and the effect is to reconnect the disconnected surface and expand the target region; the definition of A being dilated by S is: all points of the structure element S are displaced to the input set Z, and then the union is calculated to obtain the dilated image, and the calculation formula is:

[0121]

[0122] wherein A is the set of the subspace Z 2 , S is the structure element of the subspace Z 2 , A being dilated by S is the set of all displacements, and the dilation operation operator is

[0123] Erosion operation is a morphological operation to find local minimum, the main function is to separate the adhesion of objects, the boundary line is inwardly retracted; using S to erode A is to translate the set of points contained in S to A with z, the expression of erosion is:

[0124]

[0125] wherein, is an erosion operator;

[0126] Opening operation is to erode the image first and then dilate, this kind of operation can separate the adhesion of regions, and can also smooth the boundary of the image without changing the original area, its expression is:

[0127]

[0128] wherein, is an opening operator;

[0129] Closing operation is to dilate the image first and then erode, which can connect the boundary between two adjacent objects and fill the disconnected connected regions, its expression is:

[0130]

[0131] wherein, · is a closing operator;

[0132] Specific embodiment: select the three-dimensional point cloud data of teeth, the method for extracting local surface features of point cloud based on morphology, specifically comprising the following steps:

[0133] (1) delete duplicate points; according to the obvious curvature change characteristics of the junction between teeth and gums, the segmentation line of teeth is extracted, as shown in Figure 2 ; three-dimensional point cloud editing and processing software CloudCompare is used to screen out and delete duplicate points according to the minimum distance between points;

[0134] (2) uniformly sample the tooth point cloud data; the sampled points are uniformly distributed, but the position and basic contour of the point cloud are not changed, and have good precision; and by reducing the number of points in the point cloud, the processing speed of subsequent operation is improved;

[0135] (3) remove outliers of tooth point cloud data by statistical filtering; this method is only used to eliminate very obvious outliers, without setting filtering parameters in advance, and can better suppress noise of point cloud; the pretreated tooth point cloud data is as shown in Figure 3 ;

[0136] (4) Calculate the normal value of each point in the tooth point cloud; use NormalEstimation in the PCL library to calculate the normal feature of the surface where the target point is located;

[0137] (5) Estimate the curvature of each point based on the normal value using the PrincipalCurvaturesEstimation method in PCL, and determine the concave-convex condition of the local area according to the sign of the curvature value.

[0138] (6) Combine the curvature range corresponding to the local surface type to screen out the rough tooth segmentation line. The points on the concave surface are mostly negative curvature, and the points on the convex surface are mostly positive curvature, as shown in Figure 4 ; Since this curvature is not the mathematical curvature, but is calculated by least squares and PCA, the average curvature calculated is mapped to the range [0, 1], and finally the points with average curvature in the range [0.05, 0.08] are screened out, obtaining the initial demarcation point of the tooth and the gum, which also includes the edge point of the tooth;

[0139] (7) Gradual morphological filtering is used to filter out the edge points of the tooth. The size of the filter window is increased by linear iteration operation, and the window size increasing speed is adjusted according to the curvature change. The height difference threshold is adaptively adjusted according to the change of the terrain slope and the window size, and the three-dimensional morphological operator is used to filter out the non-ground points and retain the morphological features of the tooth segmentation line, as shown in Figure 5 ;

[0140] (8) Three-dimensional morphological operator is used to optimize the tooth segmentation line, as shown in Figure 6 ; The three-dimensional dilation operator reconnects the disconnected surface and expands the target area. The three-dimensional erosion operator separates the adhered objects, shrinks the boundary line inward, and eliminates the boundary points and noise points, extracting the relatively complete and smooth point cloud local surface feature area, as shown in Figure 7 .

[0141] The above only describes the preferred embodiments of the present application and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for extracting local surface features from point clouds based on morphology, characterized in that, Specifically, the following steps are included: S1: Collect point cloud data and perform data preprocessing on the point cloud data; S2: Search for neighborhood points of each point in the point cloud in 3D space, calculate normal information, and then calculate the curvature of each point based on the normal information to obtain curvature values; then filter out the target surface by referring to the curvature value range corresponding to the local surface features. S3: Progressive morphological filtering is used to continuously adjust the filtering parameters. The size of the filtering window is adjusted according to the curvature change of the selected local area to filter out non-target points in the point cloud data after curvature filtering and extract the points of the target local area. S4: The target local region is optimized by comprehensively using dilation and erosion three-dimensional morphological operators and opening and closing operations of dilation and erosion combination, and extracting local surface feature regions of complete and smooth point cloud data. The specific steps of step S3 are as follows: S31: Progressive morphological filtering is used to distinguish between ground points and non-ground points in point clouds; S32: Use the three-dimensional morphological dilation operator d and the three-dimensional morphological erosion operator e to filter out non-ground points marked on the target local surface by the filtering window, thereby extracting the ground points of the target local surface; The specific steps of step S31 are as follows: S311: Increase the filter window and perform iterative calculations by linearly increasing the window size, and then adjust the size of the filter window based on the curvature change of the selected local area of ​​the filter window. S312: Set the slope value in advance based on the average slope of the point cloud data, and extract the target surface based on the slope change; S313: Adjust the elevation difference threshold using the slope value set in step S312. Simultaneously, adaptively adjust the elevation difference threshold according to changes in terrain slope and window size. Define points with elevation difference values ​​greater than the set elevation difference threshold as non-ground points. Record the slope information of the filtered surface and calculate the elevation difference threshold based on the slope information to maintain dynamic updates of the elevation difference threshold setting.

2. The method for extracting local surface features of point clouds based on morphology according to claim 1, characterized in that, The data preprocessing for point-to-point data in step S1 specifically includes the following steps: S11: Calculate the distance between points in the point cloud data, and filter and delete duplicate points in the point cloud data based on the minimum distance between points; S12: The point cloud space is divided into several spaces using a point cloud downsampling method; then, a point is selected in each space after the division to replace all the points in that space. S13: Remove sparse outliers in the downsampled point cloud data by point cloud filtering, and then obtain the main points by performing neighborhood search on each point in the point cloud.

3. The method for extracting local surface features of point clouds based on morphology according to claim 2, characterized in that, In step S12, the following steps are performed: uniform surface sampling is performed on the point cloud data to divide the point cloud space into several spaces; then, in each space, a point is selected within a sphere with radius R to replace all points within the sphere with that point as its center and radius R.

4. The method for extracting local surface features of point clouds based on morphology according to claim 2, characterized in that, In step S13, after statistical filtering of the point cloud data, spatial topological relationships are established, and the set of neighboring points of the target point is obtained by searching for neighboring points through k nearest neighbors. If the neighboring points of the target point satisfy a Gaussian distribution, the average distance between the target point and k neighboring points is calculated. The standard range of the average distance is determined by the mean and standard deviation of the global distance. Target points whose average distance with neighboring points exceeds the standard range are defined as outliers and removed. If the average distance between the target point and k neighboring points does not exceed the standard range, the target point is retained as the main point of the model.

5. The method for extracting local surface features of point clouds based on morphology according to claim 2, characterized in that, The specific steps of step S2 are as follows: S21: Obtain the normal equation based on the least squares fitting plane, and obtain the normal value of the plane represented by the eigenvector through the normal equation; S22: Based on the acquired point cloud data and the calculated normal values ​​of each point in the point cloud data, the curvature information of the target point is approximated by the surface variation in the k-neighborhood, and the normal curvature of each point in the point cloud data is calculated. S23: Set curvature intervals by combining the curvature value ranges corresponding to different types of local surfaces, and filter out the target local surface features.

6. The method for extracting local surface features of point clouds based on morphology according to claim 5, characterized in that, In step S21, the KdTree nearest neighbor algorithm is used. KdTree is used as a data structure based on binary space partitioning to divide the entire k-dimensional vector space into several regions. By continuously performing branching judgments during the search process, the position of each local space is determined. Then, the normal vector estimation in the point cloud library is used to calculate the normal features of the surface where the target point is located. Specifically, the plane expression is transformed from the general form to the normal form equation. The optimal solution of the equation is obtained by introducing the Lagrange multiplier method. The partial derivatives of the parameters in the expression are calculated and the resulting expression is unified. The covariance matrix in the expression is decomposed into eigenvalues ​​to obtain each eigenvalue. The eigenvector corresponding to the smallest eigenvalue is the normal vector of the fitted plane. The normal vector represents the normal value of the point.

7. The method for extracting local surface features of point clouds based on morphology according to claim 5, characterized in that, In step S22, the point cloud curvature estimation method in the point cloud library is used to search for the neighborhood points of each point in the point cloud based on KdTree, and the principal curvature of each point in the point cloud is calculated. The calculation formula is as follows: Where σ represents the surface variation of the target point in the k-neighborhood; λ0 represents the surface variation along the normal vector; λ1 and λ2 represent the distribution of the target point on the tangent plane, respectively. The surface variation of the target point in the k-neighborhood represents the curvature information, and the concavity and convexity of the surface in the local region are determined according to the sign of the curvature value.

8. The method for extracting local surface features of point clouds based on morphology according to claim 1, characterized in that, The formula for increasing the filter window in step S311 when performing iterative calculations by linearly increasing the window size is: In k =nk+2; Among them, w k is the window size for the k-th filtering iteration; k is the iteration number, k = 0, 1, 2, ..., m; n is the initial window size.

9. The method for extracting local surface features of point clouds based on morphology according to claim 8, characterized in that, In step S313, assuming the slope value s remains constant, the current elevation difference threshold, window size, and slope value s are related by the following formula: Among them, h max(t)，k w represents the maximum elevation difference threshold at time t during the k-th iteration. k w is the window size for the k-th filtering iteration. k-1 h is the window size for the (k-1)th filtering iteration. k This represents the height difference threshold under the current filtering window; h0 represents the initial height difference threshold, and c is the grid cell; h max This represents the maximum height difference threshold.

10. The method for extracting local surface features of point clouds based on morphology according to claim 8, characterized in that, The formula for the three-dimensional morphological dilation operator d in step S32 is as follows: The formula for calculating the corrosion operator e is: Where w is the structuring element window; z(x, y) represents the elevation value of the x-th row and y-th column of the generated image after calculation.

11. The method for extracting local surface features of point clouds based on morphology according to claim 10, characterized in that, In step S32, if the location of a non-target point cannot be filtered out using the filtering window, the `rotate` function from the Open3d library is used to rotate the point cloud data by setting the rotation center and rotation matrix, rotating the points to be filtered out above the target surface. If a rotation of angle θ around the x-axis is chosen, with the x-coordinate remaining unchanged, and a two-dimensional rotation is performed on the yOz plane composed of y and z, then the formula is: x′=x y′=y cosθ-z sinθ z′=y sinθ+z cosθ; Then obtain the rotation matrix; 12. The method for extracting local surface features of point clouds based on morphology according to claim 10, characterized in that, The specific steps of step S4 are as follows: using the three-dimensional dilation operator d to connect the interruption points of the target surface, using the three-dimensional erosion operator e to separate the adhered objects and shrink the region, using the opening operation to smooth the boundary of the target region, using the closing operation to connect the boundary of the object, filling the disconnected connected region, and eliminating boundary points and noise points. Dilation is a morphological operation that finds local maxima, reconnecting broken points on a surface and expanding the target region. The definition of A being dilated by S is: The dilated image is obtained by inputting the set of displacements z of all points in the structuring element S, and then calculating their union. The calculation formula is: Where A is the subspace Z 2 The set of subspaces, S is a subspace of Z. 2 The structuring element, A expanded by S, is the set of all displacements. This is the expansion operator; Erosion is a morphological operation that finds local minima. Its main function is to separate adhered objects and shrink the boundary lines inward. Eroding A with S involves translating the set of points in A containing all structuring elements S by using the z-axis. The erosion expression is: in, This is the erosion operation operator; The opening operation first erodes and then dilates the image. This operation can separate contiguous regions and smooth image boundaries without changing the original area. Its expression is: in, Represents the opening operator; Closing operations first dilate the image and then erode it to connect the boundaries between two adjacent objects and fill broken connected regions. The expression is: Here, · represents the closing operator.

Images