CN111340723B - Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method - Google Patents
Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method Download PDFInfo
- Publication number
- CN111340723B CN111340723B CN202010109945.2A CN202010109945A CN111340723B CN 111340723 B CN111340723 B CN 111340723B CN 202010109945 A CN202010109945 A CN 202010109945A CN 111340723 B CN111340723 B CN 111340723B
- Authority
- CN
- China
- Prior art keywords
- points
- point
- interpolation
- elevation
- point cloud
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 83
- 238000001914 filtration Methods 0.000 title claims abstract description 34
- 230000003044 adaptive effect Effects 0.000 claims abstract description 11
- 230000014509 gene expression Effects 0.000 claims description 24
- 238000004364 calculation method Methods 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 15
- 230000011218 segmentation Effects 0.000 claims description 13
- 230000000877 morphologic effect Effects 0.000 claims description 12
- 230000009466 transformation Effects 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 11
- 230000007797 corrosion Effects 0.000 claims description 5
- 238000005260 corrosion Methods 0.000 claims description 5
- 238000012216 screening Methods 0.000 claims description 5
- 230000003628 erosive effect Effects 0.000 claims description 4
- 230000006740 morphological transformation Effects 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 241000760358 Enodes Species 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 238000013459 approach Methods 0.000 claims description 3
- 238000005452 bending Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 238000003708 edge detection Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 238000000638 solvent extraction Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 6
- 230000008569 process Effects 0.000 description 5
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 230000000750 progressive effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000002360 preparation method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 230000001502 supplementing effect Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- G06T5/70—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration by the use of local operators
- G06T5/30—Erosion or dilatation, e.g. thinning
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
Abstract
The invention discloses a terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method, which comprises the following steps: step 1, inputting original airborne laser radar three-dimensional point cloud data; step 2, denoising the three-dimensional point cloud data by adopting a method combining a kd tree and an elevation histogram; step 3, acquiring an interpolation reference point for the denoised three-dimensional point cloud data by using a two-dimensional region growing method; step 4, optimizing the interpolation reference points based on DSM line characteristics; and 5, acquiring the elevation fluctuation degree of foot points in a local area based on the relative variation coefficient for the optimized interpolation reference points, performing terrain-adaptive thin plate spline interpolation, acquiring an approximation model of a local terrain surface as a filtering classification surface, segmenting ground points and non-ground points, and outputting filtered point clouds. The filtering method is less influenced by noise, better maintains the detail characteristics of the terrain, and can maintain more stable filtering effect under various terrain data.
Description
Technical Field
The invention relates to the field of surveying and mapping science and technology, in particular to a terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method.
Background
As a novel earth surface information acquisition means, the airborne laser radar has the advantages of being all-weather, strong in penetrability, capable of acquiring earth surface three-dimensional geometric information quickly and the like compared with a traditional photogrammetry method. However, the obtained point cloud data is only the discretization representation of the ground surface space information, lacks semantic information (the type information of ground objects corresponding to the laser foot points), and increases the difficulty for subsequent application. Therefore, the filtering of the point cloud data is an important step in the application of point cloud comprehensive processing.
The existing filtering methods for airborne laser radar point clouds are divided into the following categories: mathematical morphology-based algorithms, gradient-based algorithms, triangulation-based algorithms, interpolation or fitting-based algorithms, cluster segmentation-based algorithms, and the like. The morphological method is difficult to give consideration to small and complex terrains in the area with large topographic relief, and even if a progressive window is adopted, the problem is still difficult to solve; furthermore, how to set the threshold value in the progressive window, especially the threshold value adapted to the local terrain gradient, remains a difficult point to study. The difficulty of the gradient-based algorithm is how to judge the gradient threshold of the point cloud and how to obtain the gradient threshold in a self-adaptive manner in an area with large terrain gradient change; while the reference to the slope of the terrain fracture is almost worthless. The method based on the triangulation network has a good effect under most terrain conditions, but also faces reasonable seed point selection; correctness at a steep slope, ridge or cliff; the efficiency of the construction, the search and the encryption of the triangulation network and the like. The method based on the curved surface can better simulate the terrain surface in a local area, but has the problem of large calculation amount, especially when global optimization and progressive iteration are needed; and simultaneously, the effect on a part of terrain fracture area is to be improved. Based on the segmentation method, the segmentation process is sensitive to parameters, and the segmentation effect under different parameters is greatly different; the segmentation calculation amount of the three-dimensional point cloud is large; meanwhile, the accuracy of the segmentation is particularly important, and the segmentation error is difficult to eliminate in the subsequent processing process.
The existing method mainly has three problems: the method has the advantages that firstly, the terrain adaptability problem is solved, most of traditional algorithms can achieve satisfactory effects aiming at specific terrains, but the algorithms are not good in performance in complex and mixed terrains; few algorithms can maintain the stability of the filtering effect under various complex terrain and ground feature conditions. Secondly, parameter self-adaptation problem, the algorithm with stronger terrain adaptability needs to modify parameters manually by a user according to conditions and filtering results in different terrains, the parameter modification difficulty is higher, and the accuracy of the filtering results depends on complicated parameter setting. And thirdly, the calculation efficiency problem is solved, a plurality of classical and improved algorithms are added in an iterative encryption process, the algorithm complexity is further increased, the hardware requirement on processing of massive point clouds is high, the time consumption is long, and the method is difficult to be used in practical production application.
Disclosure of Invention
The technical problem to be solved by the invention is to provide a terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method aiming at the defects in the prior art, and solve the problems of terrain adaptability and parameter self-adaptation in airborne laser radar point cloud filtering.
The invention uses the regularized thin plate spline function interpolation as a basic filtering classification surface acquisition method, adopts a two-dimensional region growing method with small calculation amount to acquire enough interpolation reference points, adopts the relative variation coefficient to acquire the undulation degree of local terrain, and controls the strictness degree of the thin plate spline interpolation according to the value of the relative variation coefficient, thereby realizing self-adaptation to the terrain with different undulation degrees. Compared with the prior art, the thin plate spline interpolation method adopted by the invention has the advantages that the retaining of the terrain details is better, the filtering error rate is lower, and the better stability can be kept under various terrain conditions; meanwhile, the redundant iterative encryption process is avoided, unnecessary interpolation calculation is reduced, and the calculation efficiency is higher. Therefore, the method has important practical value and wide application prospect. .
The technical scheme adopted by the invention for solving the technical problems is as follows:
the invention provides a terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method, which comprises the following steps:
step 1, inputting original airborne laser radar three-dimensional point cloud data, including three-dimensional coordinate information of laser foot points;
step 2, denoising the three-dimensional point cloud data by adopting a method combining a kd tree and an elevation histogram, wherein the denoising comprises removing outlier noise points and cluster noise points;
step 3, acquiring an interpolation reference point for the denoised three-dimensional point cloud data by using a two-dimensional region growing method;
step 4, optimizing the interpolation reference points based on DSM line characteristics;
and 5, acquiring the elevation fluctuation degree of foot points in a local area based on the relative variation coefficient for the optimized interpolation reference points, performing terrain-adaptive thin plate spline interpolation, acquiring an approximation model of a local terrain surface as a filtering classification surface, segmenting ground points and non-ground points, and outputting filtered point clouds.
Further, the noise point removing method combining the kd-tree and the mathematical morphology transformation described in step 2 of the present invention specifically includes:
step 2.1, constructing a three-dimensional kd tree for the point cloud, partitioning the point cloud if the number of points of the point cloud is greater than a certain threshold value, determining the size of a block according to the point density and the overall distribution of the point cloud, and then constructing the kd tree for each block respectively to accelerate the point query speed;
2.2, traversing all points of the kd tree based on the three-dimensional point cloud organized by the kd tree, and removing the current point as an isolated noise point if the number of points in the sphere range is less than a set threshold of the number of adjacent points within the sphere range for each point, wherein the radius of each point is greater than or equal to a certain threshold;
step 2.3, performing virtual gridding on the points obtained in the step 2.2, and taking the lowest elevation as the elevation of each grid unit;
and 2.4, performing morphological transformation on the gridded points, and removing cluster noise points based on the result of bottom-cap transformation.
Further, the morphological transformation of the lattice-formed points described in step 2.4 of the present invention, and the specific method of removing the cluster-like noise points based on the result of the bottom-hat transformation is as follows:
step 2.4.1, performing morphological opening operation on the gridded points, namely performing operation of corrosion first and then expansion; the concrete formula is as follows:
Gopen＝dilate(erode(Gorigin))
wherein dilate and enode refer to the expansion and erosion operations, respectively, of the grid, GoriginRepresenting the original grid, GopenRepresenting the grid after the open operation;
step 2.4.2, according to the result of the morphological open operation of the step 2.4.1, performing the morphological close operation, namely the operation of firstly expanding and then corroding; the expansion and corrosion operations are the same as in step 2.4.1; closing the result of the open operationCalculating to obtain the opening and closing grid GcloseThe concrete formula of (1) is as follows:
Gclose＝dilate(erode(Gopen))
step 2.4.3, the difference is calculated between the result of the opening and closing operation and the original grid to obtain a result G of bottom-cap transformationbottomHat(ii) a The concrete formula is as follows:
Gbottom Hat＝Gclose-Gorigin
step 2.4.4, traversing all points of the three-dimensional point cloud, finding out the grid unit of each point, and if the result grid G of the bottom cap transformation is obtainedbottom HatThe corresponding grid unit is larger than a certain threshold value, and the elevation of the current point is lower than the grid G after the opening and closing operationcloseAnd if the corresponding grid unit elevation has a certain threshold value, judging the point as a noise point and removing the noise point so as to remove the cluster-shaped low noise point.
Further, the method for obtaining the interpolation reference point by using the two-dimensional region growing method in step 3 of the present invention specifically comprises:
3.1, obtaining an elevation grid according to the denoised three-dimensional point cloud data, and acquiring grid points with the lowest elevation in a set range by using a local minimum method to serve as initial growth seed points;
step 3.2, according to the growth seed points obtained in the step 3.1, taking the elevation difference value or the gradient difference value of adjacent grids as a judgment standard, carrying out two-dimensional region growth, and bringing the growth seed points and the region growth points meeting the conditions into an interpolation reference point set for subsequent interpolation judgment; the concrete formula is as follows:
elevation growth conditions: z-zadj＜zt
The gradient increase condition is as follows:
wherein, (x, y, z) is the three-dimensional coordinate of the current point, (x)adj，yadj，zadj) Is a three-dimensional coordinate of 8 adjacent points in the neighborhood, ZtAnd StRespectively, a height difference threshold and a gradient threshold;
and 3.3, removing the appeared elevation abnormal points and isolated points from the interpolation point set by using the outlier noise removing method again according to the interpolation reference points obtained in the step 3.2 to obtain alternative interpolation reference points.
Further, the interpolation reference point optimization based on the DSM line feature described in step 4 of the present invention specifically includes:
step 4.1, according to the alternative interpolation reference point obtained in the step 3.3, DSM is constructed, and for the cavity area, a Kriging interpolation method is used for filling;
step 4.2, carrying out normalization processing on the DSM elevation, quantizing the DSM elevation into a gray map, extracting elevation fracture lines by using Canny and Scharr edge detection operators, and marking the DSM grids in which elevation fracture areas are located;
4.3, optimizing the detected fracture line, removing the line with the length less than a certain threshold value, connecting and filling two or more disconnected lines, and marking and updating the fracture area mark;
step 4.4, scanning the DSM from multiple directions, and screening out grid units which are positioned between the two broken lines and have the elevations higher than the elevations on the two sides of the broken lines; the alternate points located in these areas are removed from the interpolated reference points, and the point with the lowest elevation is taken as the set of interpolated reference points in the remaining grid of reference points.
Further, the specific method of step 4.3 of the present invention is:
step 4.3.1, traversing the end points of the lines of the extracted fracture line graph, taking the current end point as the center, indexing other end points of the non-identical lines in a window with a certain size, and processing the following conditions:
(1) if the distance between the two end points is smaller than the threshold value and the elevation difference is smaller than the threshold value, connecting the two lines;
(2) if the distance between the two end points is smaller than the threshold value and the direction difference between the two fracture lines is smaller than the threshold value, connecting the two lines;
(3) if the distance between the two end points is smaller than the threshold value and one end point is on the extension line of the other line, the two lines are connected;
and 4.3.2, after the elevation fracture lines are connected, updating the grid marks of the fracture line graphs, calculating the lengths of the fracture lines on the basis of the central points of the grids, removing the isolated lines with the lengths smaller than the threshold value, and updating the grid marks again after the processing is finished.
Further, the method for obtaining the classification surface in step 5 of the present invention is to obtain a filtering classification surface by using a thin plate spline interpolation method of a thin plate adaptive terrain, and the specific method is as follows:
step 5.1, searching the lowest point of the global elevation according to the alternative reference point set obtained in the step 4.3, and subtracting the lowest point elevation from the elevation of all alternative reference points to prevent numerical overflow;
step 5.2, calculating a global variation coefficient, namely obtaining a global mean value and a global standard deviation, and dividing the mean value by the standard deviation to obtain the variation coefficient for measuring the fluctuation condition of data; the calculation formula is as follows:
wherein, the calculation of the variation coefficient uses the relative height value z of the point to be interpolated as zi-zminCalculating relative elevation mean
And the standard deviation σ;
and 5.3, traversing all the points, acquiring a classification surface by using a thin plate spline interpolation method, and filtering the three-dimensional points.
Further, in step 5.3 of the present invention, the method of using thin-plate spline interpolation is used to obtain the classification surface and filter the three-dimensional points, and the specific method is as follows:
step 5.3.1, obtaining points in a certain size window of the neighborhood of the current point as candidate points;
step 5.3.2, sorting the candidate points according to the elevation from small to large, selecting the first N points with the minimum elevation, wherein N is more than 9, and taking an integer between 9 and 18; if the number of points in the window is less than N due to the point cloud holes or the point cloud holes in the boundary area, expanding the search window and repeating the previous step;
step 5.3.3, calculating the coefficient of variation of the N points, and then calculating the relative coefficient of variation lambda, wherein the specific calculation method is to use the coefficient of variation eta calculated herelocalDividing by the global coefficient of variation eta obtained in step 5.2globalThe concrete formula is as follows:
step 5.3.4, according to the relative variation coefficient calculated in the step 5.3.3, if the value of the relative variation coefficient approaches to 0, the elevation fluctuation of the adjacent area is very weak, which indicates that the area is very flat and close to the horizontal state, and at this time, in order to improve the calculation efficiency, the elevation average value of the N points is directly used as the horizontal classification surface of the area where the current point is located; if the value of the relative variation coefficient is not 0, entering the next step;
step 5.3.5, according to step 5.3.4, if the relative coefficient of variation is not 0, interpolating by using a regularized thin plate spline interpolation function; taking the relative variation coefficient as a regularization term coefficient of a thin plate spline function, constructing an interpolation equation and solving; the function of the ordinary thin-plate spline interpolation is as follows:
wherein, P is a trend function and expresses an approximate first-order linear trend surface, and the specific expression is as follows:
p(x，y)＝a0+a1x+a2y
wherein x and y are the plane coordinates of the current point to be interpolated, aiCalculating a constant coefficient;
summation of the second half
A basic function representing the interpolation expression is used for controlling the bending degree of the interpolation surface and ensuring the minimum curvature of the thin plate spline surface obtained by interpolation; wherein, ω isiThe coefficient is to be calculated; di＝(x-xi)2+(y-yi)2Expressing the plane distance square terms of the current point to be interpolated and each interpolation reference point; n is the number of reference points used for interpolation of the current point;
the coefficients to be solved in the interpolation function comprise n +3 items in total of a and omega, so that at least n +3 equations are needed to be solved; on this basis, 3 constraints are added:
converting the original equation into a matrix form, wherein the expression is as follows:
wherein for K in the coefficient matrixn×nThe elements of (A) are:
di,j＝(xj-xi)2+(yj-yi)2
Ki，jthe horizontal distance between each point to be interpolated; the expressions for the remaining terms are:
the regularization thin plate spline function interpolation method introduces a regularization item in a common thin plate spline interpolation equation; the general expression of the regularization term for thin-plate splines is as follows:
Kn＝λ·γ2·I
wherein λ is a regularization term coefficient, I is an identity matrix, γ is a scalar, and the expression is as follows:
computing a regularization term matrix KnThen adding the K of the original thin plate spline interpolation matrix formula to obtain a regularized thin plate spline interpolation formula, solving the coefficient obtained by the formula to obtain an exact expression of the interpolation formula, substituting the exact expression into the x and y coordinates of the point to be interpolated to obtain the interpolation Z coordinate of the point to be interpolated;
5.3.6, calculating the interpolation height of the current point position according to the solved thin plate spline function model, taking the interpolation height as a segmentation elevation, determining points higher than the segmentation elevation by a certain height as non-ground points, and classifying the rest points as ground points; and traversing all the points according to the mode to finally obtain the filtered point cloud.
The invention has the following beneficial effects:
1) compared with other methods for obtaining the terrain surface by fitting and interpolating, the method has the advantages of simple calculation and better terrain approximation degree by using a flexible thin plate spline function with strong adaptability as a fitting method of the terrain surface;
2) the method has the advantages that the alternative interpolation reference points are rapidly acquired by using a mode based on region growing, the constraint conditions of region growing and the detection of region growing results are increased, and the uniformity and the integrity of the interpolation points are ensured;
3) further screening the alternative points based on DSM line characteristic constraint conditions aiming at the growth result of the two-dimensional area to ensure the accuracy of the final interpolation reference point;
4) the relative variation coefficient is used as the regularization coefficient of the interpolation function, so that the terrains with different fluctuation degrees are ensured to use thin plate spline interpolation means with different strictness degrees, the universality of the interpolation on various terrains is enhanced, and meanwhile, the method has better noise resistance and noise resistance.
Drawings
The invention will be further described with reference to the accompanying drawings and examples, in which:
FIG. 1 is a flow chart of an embodiment of the present invention;
FIG. 2 is a schematic diagram of an embodiment of the present invention for obtaining an interpolated reference using region growing method engagement;
FIG. 3 is a schematic diagram of approximating a local topographic surface using thin-plate splines, in accordance with an embodiment of the present invention
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
As shown in FIG. 1, the terrain adaptive onboard LiDAR point cloud regularization thin-plate spline interpolation filtering method in the embodiment of the invention comprises the following steps:
step 1, data preparation. The data required to be input by the method comprises the following steps: the original airborne laser radar three-dimensional point cloud comprises information such as three-dimensional coordinates of laser foot points; the filtering parameters and the numerical values of the parameters input by a user are flexibly set according to the actual condition of the point cloud, the requirements of point cloud filtering precision, hardware environment and the like.
Step 2, removing outlier noise points and cluster noise points by combining the kd tree and the elevation histogram;
step 3, obtaining an interpolation reference point by using two-dimensional region growth;
step 4, optimizing an interpolation reference point based on DSM line characteristic constraint;
and 5, acquiring the elevation fluctuation degree of foot points in a local area based on the relative variation coefficient, performing thin plate spline interpolation of a self-adaptive terrain, acquiring an approximation model of a local terrain surface as a filtering classification surface, segmenting ground points and non-ground points, and outputting a filtered point cloud.
The details of each step are as follows:
the step 1 process belongs to the preparation step and is not described in detail here.
And step 2 is a noise point removing method combining a kd tree and mathematical morphology transformation.
Step 2.1, constructing a three-dimensional kd tree for the point cloud, if the number of points of the point cloud is huge, partitioning the point cloud (the size of a block can be determined according to the point density and the overall distribution of the point cloud), and then constructing the kd tree for each block respectively to accelerate the point query speed;
2.2, traversing all points of the tree based on the three-dimensional point cloud organized by the kd tree, and removing the current point as an isolated noise point if the number of points in the sphere range is less than 3-10 (the threshold value of the number of adjacent points can be self-determined) in the sphere range with the radius R of more than or equal to 5m (the radius can be automatically increased and decreased according to the actual situation of the point cloud) for each point;
step 2.3, performing virtual gridding on the points obtained in the step 2.2, and taking the lowest elevation as the elevation of each grid unit;
step 2.4, morphologically transforming the grid points, and removing cluster noise points based on the result of bottom-hat transformation, specifically comprising the following steps:
and 2.4.1, performing morphological opening operation on the gridded points, namely performing operation of corrosion first and expansion later.
The specific formula is as follows:
Gopen＝dilate(erode(Gorigin))
where the dilate and the enode refer to the expansion and erosion operations, respectively, on the grid.
And 2.4.2, performing morphological closing operation according to the morphological opening operation result of the step 2.4.1, namely performing the operation of firstly expanding and then corroding. The expansion and erosion operations are the same as the previous step. The formula for obtaining the open-close grid by performing the close operation on the open operation result is as follows:
Gclose＝dilate(erode(Gopen))
and 2.4.3, performing difference calculation on the result of the opening and closing operation (firstly performing morphological opening operation and then performing morphological closing operation) and the original grid to obtain a result of bottom-cap transformation. The specific formula is as follows:
Gbottom Hat＝Gclose-Gorigin
step 2.4.4, traversing all points of the three-dimensional point cloud, finding out the grid unit of each point, and if the result grid G of the bottom cap transformation is obtainedbottom HatThe corresponding grid unit is larger than a certain threshold value, and the elevation of the current point is lower than the grid G after the opening and closing operationcloseAnd if the corresponding grid unit elevation is a certain threshold value, judging the point as a noise point and removing the noise point. In this way, clustered low noise spots can be removed.
And 2.4.5, updating the elevation grid.
And 3, obtaining an interpolation reference point and performing post-processing screening by using two-dimensional region growth.
3.1, on the basis of obtaining the elevation grid, obtaining grid points with the lowest elevation in a larger range by using a local minimum method, and taking the grid points as initial growth seed points;
and 3.2, according to the growing seed points obtained in the step 3.1, taking the elevation difference value or the gradient difference value of the adjacent grids as a judgment standard, carrying out two-dimensional region growing, and bringing the seed points and the region growing points meeting the conditions into an interpolation reference point set for subsequent interpolation judgment.
Elevation growth conditions: z-zadj＜zt
The gradient increase condition is as follows:
wherein, (x, y, z) is the three-dimensional coordinate of the current point, (x)adj，yadj，zadj) Three-dimensional coordinates of adjacent points in the 8 neighborhoods. ZtAnd StRespectively, a head and a grade threshold.
And 3.3, removing possible elevation abnormal points and isolated points from the interpolation point set by using an outlier noise removing method again according to the interpolation reference points obtained in the step 3.2 to obtain alternative interpolation reference points.
And 4, optimizing the acquired reference points by using DSM line characteristics to delete some ground object points close to the elevation of the ground point in a specific direction, such as ground entrances of the viaduct.
Step 4.1, according to the alternative interpolation reference point obtained in the step 3.3, DSM is constructed, and for the cavity area, a Kriging interpolation method is used for filling;
step 4.2, carrying out elevation normalization on the DSM, quantizing the elevation normalization into a gray scale map, extracting elevation fracture lines by using edge detection operators such as Canny and Scharr, and marking the DSM grids where elevation fracture areas are located;
4.3, optimizing the detected fracture line, removing the shorter line, connecting and supplementing two or more disconnected lines, and marking and updating the fracture area mark; the method comprises the following specific steps:
step 4.3.1, traversing the end points of the extracted fracture line graph, taking the current end point as the center, indexing other end points of the non-identical line in a window (5 x 5) with a certain size, and processing the following conditions:
(5) if the distance between the two end points is close and the elevation difference is small, the two lines are connected;
(6) the two end points are close to each other, and the two fracture lines are close to each other in direction, so that the two lines are connected;
the two end points are close to each other, and one end point is arranged on the extension line of the other line, so that the two lines are connected;
step 4.3.2, after connecting elevation broken lines, updating the grid marks of the broken line graphs, calculating the lengths of the broken lines on the basis of the central points of the grids, removing isolated lines (a large-scale buffer area is constructed on the basis of the whole line, and if no end points of other lines exist in the whole buffer area, the isolated lines are judged to be the isolated lines), enabling the lines to be short in length (the size of a grid unit needs to be considered in the setting of a length threshold, and the calculated value of the common length is 1-2 meters so as to ensure that small-sized lines in residences and the like cannot be deleted), and updating the grid marks again after the processing is finished;
and 4.4, scanning the DSM in multiple directions, and screening out grid units which are positioned between the two broken lines and have the elevations which are obviously higher than the elevations on the two sides of the broken lines. The alternative points in these areas are removed from the interpolation reference points, and in the remaining reference point grids, the point with the lowest elevation is included in the interpolation reference point set (if there is a higher requirement for the number of interpolation reference points, the number of interpolation reference point sets included in each qualified reference point grid can be increased as appropriate, provided that the elevation range of all the reference points in a single grid is within a certain range, and there is no large difference).
And 5, obtaining a filtering classification surface by using a thin plate spline interpolation method of the thin plate adaptive terrain.
Step 5.1, searching the lowest point of the global elevation according to the alternative reference point set obtained in the step 4.3, and subtracting the elevation of the lowest point from the elevation of all alternative reference points (preventing numerical overflow);
and 5.2, calculating the global variation coefficient, namely obtaining a global mean value and a global standard deviation, and dividing the mean value by the global standard deviation to obtain the variation coefficient which can be used for measuring the fluctuation condition of the data. The calculation formula is as follows:
the coefficient of variation is calculated using the relative height z ═ z of the point to be interpolatedi-zminCalculating relative elevation mean
And the standard deviation σ.
And 5.3, traversing all the points, acquiring a classification surface by using a thin plate spline interpolation method, and filtering the three-dimensional points. The method comprises the following specific steps:
step 5.3.1, obtaining points (<25) in a window of 5 x 5 of a current point neighborhood as candidate points;
step 5.3.2, sorting the candidate points according to the elevation from small to large, and selecting the first N points with the minimum elevation (N is more than 9, and can be an integer between 9 and 18 generally); if the point cloud holes or the point cloud holes in the boundary area cause that the number of points in the 5 x 5 window is less than N, expanding the search window and repeating the previous step;
step 5.3.3, calculating the coefficient of variation of the N points, wherein the method is the same as the step 5.1 and the step 5.2; and calculating the relative variation coefficient, wherein the specific calculation method is to divide the calculated variation coefficient by the global variation coefficient obtained in the step 5.2, and the specific formula is as follows:
step 5.3.4, according to the relative variation coefficient calculated in the step 5.3.3, if the value of the relative variation coefficient approaches to 0, the elevation fluctuation of the adjacent area is very weak (the area is very flat and is close to the horizontal state), and at this time, in order to improve the calculation efficiency, the elevation average value of the N points is directly used as the horizontal classification surface of the area where the current point is located; if the value of the relative variation coefficient is not 0, entering the next step;
step 5.3.5, according to step 5.3.4, if the relative coefficient of variation is not 0, interpolating by using a regularized thin plate spline interpolation function; and (5) taking the relative variation coefficient as a regularization term coefficient of a thin plate spline function, constructing an interpolation equation and solving. The function of the ordinary thin-plate spline interpolation is as follows:
wherein, P is a trend function and expresses an approximate first-order linear trend surface, and the specific expression is as follows:
p(x，y)＝a0+a1x+a2y
x and y are the plane coordinates of the current point to be interpolated, aiThe constant coefficient is to be calculated.
Summation of the second half
A basic function representing the interpolation expression is used for controlling the bending degree of the interpolation surface and ensuring the minimum curvature of the thin plate spline surface obtained by interpolation; wherein, wiThe coefficient is to be calculated; di＝(x-xi)2+(y-yi)2Expressing the plane distance square terms of the current point to be interpolated and each interpolation reference point; n is the number of reference points used for interpolation of the current point;
the coefficients to be solved in the interpolation function comprise n +3 items in total of a and omega, so that at least n +3 equations are needed to be solved; on this basis, 3 constraints are added:
converting the original equation into a matrix form, wherein the expression is as follows:
wherein for K in the coefficient matrixn×nThe elements of (A) are:
di,j＝(xj-xi)2+(yj-yi)2
Ki，jthe horizontal distance between each point to be interpolated; the expressions for the remaining terms are:
the regularization thin plate spline function interpolation method introduces a regularization item in a common thin plate spline interpolation equation; the general expression of the regularization term for thin-plate splines is as follows:
Kn＝λ·γ2·I
wherein λ is a regularization term coefficient, I is an identity matrix, γ is a scalar, and the expression is as follows:
computing a regularization term matrix KnThen adding the K of the original thin plate spline interpolation matrix formula to obtain a regularized thin plate spline interpolation formula, solving the coefficient obtained by the formula to obtain an exact expression of the interpolation formula, substituting the exact expression into the x and y coordinates of the point to be interpolated to obtain the interpolation z coordinate of the point to be interpolated;
5.3.6, calculating the interpolation height of the current point position according to the solved thin plate spline function model, taking the interpolation height as a segmentation elevation, determining points higher than the segmentation elevation by a certain height as non-ground points, and classifying the rest points as ground points; and traversing all the points according to the mode to finally obtain the filtered point cloud.
It will be understood that modifications and variations can be made by persons skilled in the art in light of the above teachings and all such modifications and variations are intended to be included within the scope of the invention as defined in the appended claims.
Claims (4)
1. A terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method is characterized by comprising the following steps:
step 1, inputting original airborne laser radar three-dimensional point cloud data, including three-dimensional coordinate information of laser foot points;
step 2, denoising the three-dimensional point cloud data by adopting a method combining a kd tree and an elevation histogram, wherein the denoising comprises removing outlier noise points and cluster noise points;
step 3, acquiring an interpolation reference point for the denoised three-dimensional point cloud data by using a two-dimensional region growing method;
the method for obtaining the interpolation reference point by using the two-dimensional region growing method in the step 3 comprises the following specific steps:
3.1, obtaining an elevation grid according to the denoised three-dimensional point cloud data, and acquiring grid points with the lowest elevation in a set range by using a local minimum method to serve as initial growth seed points;
step 3.2, according to the growth seed points obtained in the step 3.1, taking the elevation difference value or the gradient difference value of adjacent grids as a judgment standard, carrying out two-dimensional region growth, and bringing the growth seed points and the region growth points meeting the conditions into an interpolation reference point set for subsequent interpolation judgment; the concrete formula is as follows:
elevation growth conditions: z-zadj＜zt
The gradient increase condition is as follows:
wherein, (x, y, z) is the three-dimensional coordinate of the current point, (x)adj，yadj，zadj) Is a three-dimensional coordinate of 8 adjacent points in the neighborhood, ZtAnd SlRespectively, a height difference threshold and a gradient threshold;
3.3, removing the appeared elevation abnormal points and isolated points from the interpolation point set by using an outlier noise removing method again according to the interpolation reference points obtained in the step 3.2 to obtain alternative interpolation reference points;
step 4, optimizing the interpolation reference points based on DSM line characteristics;
the interpolation reference point optimization based on the DSM line characteristics in the step 4 comprises the following specific steps:
step 4.1, according to the alternative interpolation reference point obtained in the step 3.3, DSM is constructed, and for the cavity area, a Kriging interpolation method is used for filling;
step 4.2, carrying out normalization processing on the DSM elevation, quantizing the DSM elevation into a gray map, extracting elevation fracture lines by using Canny and Scharr edge detection operators, and marking the DSM grids in which elevation fracture areas are located;
4.3, optimizing the detected fracture line, removing the line with the length less than a certain threshold value, connecting and filling two or more disconnected lines, and marking and updating the fracture area mark;
step 4.4, scanning the DSM from multiple directions, and screening out grid units which are positioned between the two broken lines and have the elevations higher than the elevations on the two sides of the broken lines; removing the alternative points in the areas from the interpolation reference points, and taking the point with the lowest elevation in the residual reference point grids into the interpolation reference point set;
step 5, acquiring the elevation fluctuation degree of foot points in a local area based on the relative variation coefficient for the optimized interpolation reference point, performing terrain-adaptive thin plate spline interpolation, acquiring an approximation model of a local terrain surface as a filtering classification surface, segmenting ground points and non-ground points, and outputting filtered point clouds;
the method for obtaining the classification surface in the step 5 is to obtain a filtering classification surface by using a thin plate spline interpolation method of a thin plate self-adaptive terrain, and the specific method is as follows:
step 5.1, searching the lowest point of the global elevation according to the alternative reference point set obtained in the step 4.3, and subtracting the lowest point elevation from the elevation of all alternative reference points to prevent numerical overflow;
step 5.2, calculating a global variation coefficient, namely obtaining a global mean value and a global standard deviation, and dividing the mean value by the standard deviation to obtain the variation coefficient for measuring the fluctuation condition of data; the calculation formula is as follows:
wherein, the calculation of the variation coefficient uses the relative height value z of the point to be interpolated as zi-zminCalculating relative elevation mean
And the standard deviation σ;
step 5.3, traversing all the points, obtaining a classification surface by using a thin plate spline interpolation method, and filtering the three-dimensional points;
the method for obtaining the classification surface by using the thin plate spline interpolation in the step 5.3 is used for filtering the three-dimensional points, and the specific method is as follows:
step 5.3.1, obtaining points in a certain size window of the neighborhood of the current point as candidate points;
step 5.3.2, sorting the candidate points according to the elevation from small to large, selecting the first N points with the minimum elevation, wherein N is more than 9, and taking an integer between 9 and 18; if the number of points in the window is less than N due to the point cloud holes or the point cloud holes in the boundary area, expanding the search window and repeating the previous step;
step 5.3.3, calculating the coefficient of variation of the N points, and then calculating the relative coefficient of variation lambda, wherein the specific calculation method is to use the coefficient of variation eta calculated herelocalDividing by the global coefficient of variation eta obtained in step 5.2globalThe concrete formula is as follows:
step 5.3.4, according to the relative variation coefficient calculated in the step 5.3.3, if the value of the relative variation coefficient approaches to 0, the elevation fluctuation of the adjacent area is very weak, which indicates that the area is very flat and close to the horizontal state, and at this time, in order to improve the calculation efficiency, the elevation average value of the N points is directly used as the horizontal classification surface of the area where the current point is located; if the value of the relative variation coefficient is not 0, entering the next step;
step 5.3.5, according to step 5.3.4, if the relative coefficient of variation is not 0, interpolating by using a regularized thin plate spline interpolation function; taking the relative variation coefficient as a regularization term coefficient of a thin plate spline function, constructing an interpolation equation and solving; the function of the ordinary thin-plate spline interpolation is as follows:
wherein, P is a trend function and expresses an approximate first-order linear trend surface, and the specific expression is as follows:
p(x，y)＝a0|a1x|a2y
wherein x and y are the plane coordinates of the current point to be interpolated, aiCalculating a constant coefficient;
summation of the second half
A basic function representing the interpolation expression is used for controlling the bending degree of the interpolation surface and ensuring the minimum curvature of the thin plate spline surface obtained by interpolation; wherein, wiThe coefficient is to be calculated; di＝(x-xi)2+(y-yi)2Expressing the plane distance square terms of the current point to be interpolated and each interpolation reference point; n is the number of reference points used for interpolation of the current point;
the coefficients to be solved in the interpolation function comprise n +3 items of a and w, so that at least n +3 equations are needed to be solved; on this basis, 3 constraints are added:
converting the original equation into a matrix form, wherein the expression is as follows:
wherein for K in the coefficient matrixn×nThe elements of (A) are:
di,j＝(xj-xi)2+(yj-yi)2
Ki，jthe horizontal distance between each point to be interpolated; the expressions for the remaining terms are:
the regularization thin plate spline function interpolation method introduces a regularization item in a common thin plate spline interpolation equation; the general expression of the regularization term for thin-plate splines is as follows:
Kn＝λ·γ2·I
wherein λ is a regularization term coefficient, I is an identity matrix, γ is a scalar, and the expression is as follows:
computing a regularization term matrix KnThen adding the K of the original thin plate spline interpolation matrix formula to obtain a regularized thin plate spline interpolation formula, solving the coefficient obtained by the formula to obtain an exact expression of the interpolation formula, substituting the exact expression into the x and y coordinates of the point to be interpolated to obtain the interpolation z coordinate of the point to be interpolated;
5.3.6, calculating the interpolation height of the current point position according to the solved thin plate spline function model, taking the interpolation height as a segmentation elevation, determining points higher than the segmentation elevation by a certain height as non-ground points, and classifying the rest points as ground points; and traversing all the points according to the mode to finally obtain the filtered point cloud.
2. The terrain adaptive onboard LiDAR point cloud regularization thin-plate spline interpolation filtering method according to claim 1, wherein the noise point removing method combining kd-trees and mathematical morphology transformation in step 2 is specifically:
step 2.1, constructing a three-dimensional kd tree for the point cloud, partitioning the point cloud if the number of points of the point cloud is greater than a certain threshold value, determining the size of a block according to the point density and the overall distribution of the point cloud, and then constructing the kd tree for each block respectively to accelerate the point query speed;
2.2, traversing all points of the kd tree based on the three-dimensional point cloud organized by the kd tree, and removing the current point as an isolated noise point if the number of points in the sphere range is less than a set threshold of the number of adjacent points within the sphere range for each point, wherein the radius of each point is greater than or equal to a certain threshold;
step 2.3, performing virtual gridding on the points obtained in the step 2.2, and taking the lowest elevation as the elevation of each grid unit;
and 2.4, performing morphological transformation on the gridded points, and removing cluster noise points based on the result of bottom-cap transformation.
3. The terrain adaptive onboard LiDAR point cloud regularization thin-plate spline interpolation filtering method according to claim 2, wherein the morphological transformation of the gridded points in step 2.4 is performed by a specific method for removing cluster-like noise points based on the result of the bottom-hat transformation:
step 2.4.1, performing morphological opening operation on the gridded points, namely performing operation of corrosion first and then expansion; the concrete formula is as follows:
Gopen＝dilate(erode(Gorigin))
wherein dilate and enode refer to the expansion and erosion operations, respectively, of the grid, GoriginRepresenting the original grid, GopenRepresenting the grid after the open operation;
step 2.4.2, according to the result of the morphological open operation of the step 2.4.1, performing the morphological close operation, namely the operation of firstly expanding and then corroding; the expansion and corrosion operations are the same as in step 2.4.1; performing a closing operation on the opening operation result to obtain an opening and closing grid GcloseThe concrete formula of (1) is as follows:
Gclose＝dilate(erode(Gopen))
step 2.4.3, the difference is calculated between the result of the opening and closing operation and the original grid to obtain a result G of bottom-cap transformationbottomHat(ii) a The concrete formula is as follows:
GbottomHat＝Gclose-Gorigin
step 2.4.4, traversing all points of the three-dimensional point cloud, finding out the grid unit of each point, and if the result grid G of the bottom cap transformation is obtainedbottomHatThe corresponding grid unit is larger than a certain threshold value, and the elevation of the current point is lower than the grid G after the opening and closing operationcloseAnd if the corresponding grid unit elevation has a certain threshold value, judging the point as a noise point and removing the noise point so as to remove the cluster-shaped low noise point.
4. The terrain adaptive onboard LiDAR point cloud regularization thin plate spline interpolation filtering method according to claim 1, characterized in that the specific method of step 4.3 is:
step 4.3.1, traversing the end points of the lines of the extracted fracture line graph, taking the current end point as the center, indexing other end points of the non-identical lines in a window with a certain size, and processing the following conditions:
(1) if the distance between the two end points is smaller than the threshold value and the elevation difference is smaller than the threshold value, connecting the two lines;
(2) if the distance between the two end points is smaller than the threshold value and the direction difference between the two fracture lines is smaller than the threshold value, connecting the two lines;
(3) if the distance between the two end points is smaller than the threshold value and one end point is on the extension line of the other line, the two lines are connected;
and 4.3.2, after the elevation fracture lines are connected, updating the grid marks of the fracture line graphs, calculating the lengths of the fracture lines on the basis of the central points of the grids, removing the isolated lines with the lengths smaller than the threshold value, and updating the grid marks again after the processing is finished.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010109945.2A CN111340723B (en) | 2020-02-23 | 2020-02-23 | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010109945.2A CN111340723B (en) | 2020-02-23 | 2020-02-23 | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111340723A CN111340723A (en) | 2020-06-26 |
| CN111340723B true CN111340723B (en) | 2022-04-15 |
Family
ID=71186968
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010109945.2A Active CN111340723B (en) | 2020-02-23 | 2020-02-23 | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111340723B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111783721B (en) * | 2020-07-13 | 2021-07-20 | 湖北亿咖通科技有限公司 | Lane line extraction method of laser point cloud and electronic equipment |
| CN112231848B (en) * | 2020-11-09 | 2023-04-07 | 北京理工大学 | Method and system for constructing vehicle spraying model |
| CN113239136B (en) * | 2021-05-17 | 2023-12-19 | 北京车和家信息技术有限公司 | Data processing method, device, equipment and medium |
| CN115100363B (en) * | 2022-08-24 | 2022-11-25 | 中国科学院地理科学与资源研究所 | Underground abnormal body three-dimensional modeling method and device based on ground penetrating radar |
| CN115861571B (en) * | 2023-01-18 | 2023-04-28 | 武汉大学 | Semantic perception triangle network model building entity reconstruction method |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103093466A (en) * | 2013-01-21 | 2013-05-08 | 武汉大学 | Building three-dimensional change detection method based on LiDAR point cloud and image |
| CN103745436A (en) * | 2013-12-23 | 2014-04-23 | 西安电子科技大学 | LiDar point cloud data morphological filtering method based on area prediction |
| CN105026955A (en) * | 2012-12-28 | 2015-11-04 | 诺基亚技术有限公司 | A method and apparatus for de-noising data from a distance sensing camera |
| CN106154247A (en) * | 2016-06-24 | 2016-11-23 | 南京林业大学 | A kind of multiple dimensioned Full wave shape laser radar data optimizes decomposition method |
| CN110119438A (en) * | 2019-04-23 | 2019-08-13 | 东华理工大学 | Airborne LiDAR point cloud filtering method based on Active Learning |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10302768B2 (en) * | 2016-05-09 | 2019-05-28 | Microsoft Technology Licensing, Llc | Multipath signal removal in time-of-flight camera apparatus |
-
2020
- 2020-02-23 CN CN202010109945.2A patent/CN111340723B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105026955A (en) * | 2012-12-28 | 2015-11-04 | 诺基亚技术有限公司 | A method and apparatus for de-noising data from a distance sensing camera |
| CN103093466A (en) * | 2013-01-21 | 2013-05-08 | 武汉大学 | Building three-dimensional change detection method based on LiDAR point cloud and image |
| CN103745436A (en) * | 2013-12-23 | 2014-04-23 | 西安电子科技大学 | LiDar point cloud data morphological filtering method based on area prediction |
| CN106154247A (en) * | 2016-06-24 | 2016-11-23 | 南京林业大学 | A kind of multiple dimensioned Full wave shape laser radar data optimizes decomposition method |
| CN110119438A (en) * | 2019-04-23 | 2019-08-13 | 东华理工大学 | Airborne LiDAR point cloud filtering method based on Active Learning |
Non-Patent Citations (2)
| Title |
|---|
| Automatic and Unsupervised Water Body Extraction Based on Spectral-Spatial Features Using GF-1 Satellite Imagery;Yongjun Zhang 等;《IEEE GEOSCIENCE AND REMOTE SENSING LETTERS》;20190630;第927-931页 * |
| 自适应移动盒子的机载LiDAR点云去噪算法;李炼 等;《测绘科学》;20160430;第144-147页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111340723A (en) | 2020-06-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111340723B (en) | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method | |
| CN106529469B (en) | Unmanned aerial vehicle-mounted LiDAR point cloud filtering method based on self-adaptive gradient | |
| CN106199557B (en) | A kind of airborne laser radar data vegetation extracting method | |
| CN110363861B (en) | Laser radar point cloud-based field crop three-dimensional reconstruction method | |
| Zhang et al. | Filtering airborne LiDAR data by embedding smoothness-constrained segmentation in progressive TIN densification | |
| CN112784403B (en) | Numerical simulation method for establishing discrete element model of jointed rock mass based on point cloud data | |
| CN110443810B (en) | Point cloud plane segmentation method based on quick adjacent voxel query | |
| CN114332366A (en) | Digital city single house point cloud facade 3D feature extraction method | |
| CN110335352B (en) | Double-element multi-resolution hierarchical filtering method for airborne laser radar point cloud | |
| CN111598780B (en) | Terrain adaptive interpolation filtering method suitable for airborne LiDAR point cloud | |
| Hui et al. | A mean shift segmentation morphological filter for airborne LiDAR DTM extraction under forest canopy | |
| CN116310849B (en) | Tree point cloud monomerization extraction method based on three-dimensional morphological characteristics | |
| CN111950589B (en) | Point cloud region growing optimization segmentation method combined with K-means clustering | |
| CN112669333A (en) | Single tree information extraction method | |
| CN113379919A (en) | Vegetation canopy height rapid extraction method based on unmanned aerial vehicle RGB camera | |
| CN115049925A (en) | Method for extracting field ridge, electronic device and storage medium | |
| CN114862715A (en) | TIN (triangulated irregular network) progressive encryption denoising method fusing terrain feature semantic information | |
| CN115063555A (en) | Method for extracting vehicle-mounted LiDAR point cloud street tree growing in Gaussian distribution area | |
| CN113409332B (en) | Building plane segmentation method based on three-dimensional point cloud | |
| CN111861946B (en) | Adaptive multi-scale vehicle-mounted laser radar dense point cloud data filtering method | |
| CN114463338A (en) | Automatic building laser foot point extraction method based on graph cutting and post-processing | |
| CN113177897A (en) | Rapid lossless filtering method for disordered 3D point cloud | |
| CN112945196A (en) | Strip mine step line extraction and slope monitoring method based on point cloud data | |
| Omidalizarandi et al. | Segmentation and classification of point clouds from dense aerial image matching | |
| CN116597199A (en) | Point cloud tree classification method and system based on airborne LiDAR |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |