CN110246092B - Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope - Google Patents
Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope Download PDFInfo
- Publication number
- CN110246092B CN110246092B CN201910365632.0A CN201910365632A CN110246092B CN 110246092 B CN110246092 B CN 110246092B CN 201910365632 A CN201910365632 A CN 201910365632A CN 110246092 B CN110246092 B CN 110246092B
- Authority
- CN
- China
- Prior art keywords
- point
- neighborhood
- distance
- points
- distribution
- 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 description 26
- 230000001186 cumulative effect Effects 0.000 claims description 14
- 238000005315 distribution function Methods 0.000 claims description 12
- 238000004364 calculation method Methods 0.000 claims description 7
- AYFVYJQAPQTCCC-GBXIJSLDSA-N L-threonine Chemical compound C[C@@H](O)[C@H](N)C(O)=O AYFVYJQAPQTCCC-GBXIJSLDSA-N 0.000 claims description 3
- 230000000717 retained effect Effects 0.000 claims description 3
- 238000010845 search algorithm Methods 0.000 claims description 2
- 238000000926 separation method Methods 0.000 claims 1
- 238000001914 filtration Methods 0.000 abstract description 33
- 230000000694 effects Effects 0.000 abstract description 21
- 230000002159 abnormal effect Effects 0.000 abstract description 6
- 238000004422 calculation algorithm Methods 0.000 abstract description 6
- 238000010276 construction Methods 0.000 abstract description 6
- 230000006870 function Effects 0.000 description 12
- 238000010586 diagram Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 238000010146 3D printing Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- FFBHFFJDDLITSX-UHFFFAOYSA-N benzyl N-[2-hydroxy-4-(3-oxomorpholin-4-yl)phenyl]carbamate Chemical compound OC1=C(NC(=O)OCC2=CC=CC=C2)C=CC(=C1)N1CCOCC1=O FFBHFFJDDLITSX-UHFFFAOYSA-N 0.000 description 1
- 230000002146 bilateral effect Effects 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000000428 dust Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- 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
-
- 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/20—Special algorithmic details
- G06T2207/20024—Filtering details
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Complex Calculations (AREA)
Abstract
Aiming at the series problems that the threshold value of the number of the current adjacent points cannot be ensured by experience setting, the fitting capability of a common filtering algorithm is limited, certain limitation exists and the like, the invention provides a denoising strategy taking the distance and the slope of the adjacent points into consideration for the ground three-dimensional Laser point cloud (TLS), and realizes the construction of a one-dimensional Gamma denoising function taking the distance and the slope of the adjacent points into consideration independently and the construction of a two-dimensional Gamma denoising function taking the distance and the slope of the adjacent points into consideration under the support of the construction of a k-d tree index. The two constraint conditions of the filter provided by the invention consider the conditions of the neighborhood quantity around the query point, the distance mean value between the neighborhood point and the central point, the distribution density and the like from different directions, and the filter is effective for TLS data abnormal values, no matter the TLS data abnormal values are sparse abnormal values or noise types such as isolated noise points and the like, and can achieve the technical effect of automatically filtering the noise points.
Description
Technical Field
The invention belongs to the field of surveying and mapping and geographic information science.
Background
In recent years, computer technology and three-dimensional digital information launching stations are rapidly advanced, three-dimensional point model data of a target object are obtained in a three-dimensional silk wonderful mode and are subjected to reverse modeling, and the research hotspot is at present. Three-dimensional point modeling is particularly used in the fields of 3D printing, virtual reality, computer modeling, and the like. However, in the process of obtaining point data through the three-dimensional scanner, random errors always occur in the measured three-dimensional data, and certain noise information exists in the measured data information due to the irregularity of the surface of the measured object, so that in the modeling process, the accuracy detection and denoising work on the obtained data plays an important role in the accuracy of the modeling result. In the prior art, a plurality of denoising methods for point cloud models exist, for example, 1) a moving surface fitting point cloud denoising method fully utilizes the continuity of ground objects or terrain surfaces in a micro range, namely, the method has a second-order conductibility and can be approximated by a quadric surface. 2) Point cloud denoising methods based on mathematical morphology theory extract image features with the size of features in a "reduced" (eroded) or "increased" (dilated) image. 3) The bilateral filtering point cloud denoising method can effectively reduce image noise and maintain the algorithm of edge detail information to the maximum extent. 4) An isolated noise point eliminating method based on statistical analysis is characterized in that point cloud data of ground buildings are organized by using a k-d tree firstly, and then a threshold value is automatically set by using a statistical analysis method based on a normal distribution principle so as to filter isolated noise points.
Aiming at the series problems that the threshold value of the number of the current adjacent points cannot be ensured by experience setting, the fitting capability of a common filtering algorithm is limited, certain limitation exists and the like, the invention provides a denoising strategy taking the distance and the slope of the adjacent points into consideration for the ground three-dimensional Laser point cloud (TLS), and realizes the construction of a one-dimensional Gamma denoising function taking the distance and the slope of the adjacent points into consideration independently and the construction of a two-dimensional Gamma denoising function taking the distance and the slope of the adjacent points into consideration under the support of the construction of a k-d tree index. The two constraint conditions of the filter consider the conditions of the number of the neighborhood around the query point, the distance mean value between the neighborhood point and the central point, the distribution density and the like from different directions, the filter is effective for TLS data abnormal values, whether the TLS data abnormal values are sparse abnormal values or noise types such as isolated noise points, and the technical effect of automatically filtering the noise points can be achieved.
Disclosure of Invention
The invention relates to a three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope, which comprises the following steps:
(1) reading original point cloud data set P ═ P 1 ,P 2 ,…,P n Storing the obtained data into a three-dimensional array, performing variance calculation on each dimension, selecting a dimension t (t is 1,2 and 3) with the largest variance, selecting a median mu of all values in the dimension t as a first comparison object, namely a dividing axis, obtaining two subsets, and storing the subsets into a tree node N (t, mu);
(2) repeating the process (1) in step 1 for the two subsets obtained until the original data set P ═ P 1 ,P 2 ,…,P n It cannot be subdivided; if a subset cannot be subdivided in the process, the data in the subset is saved to the leaf node.
(1) data point P to be inquired i (x, y, z) comparing with the mu value in each tree node N (t, mu) from the root node, and entering a right subtree query if the value in the dimension t is greater than mu; otherwise, entering left sub-tree query, and calculating P after reaching leaf node i (x, y, z) the distance d between the data points stored on this leaf node and the corresponding target point coordinates P are recorded o (x,y,z)；
(2) Backtracking, i.e. determining whether there is a distance P in an unaccessed data point i (x, y, z) is closer, and P is judged i Whether the distance between (x, y, z) and the unaccessed data point under the parent node to which the (x, y, z) belongs is larger than d, if so, the distance between the (x, y, z) and P does not exist in the branch i (x, y, z) closer points; otherwise, there is a departure P i (x, y, z) closer point, then visit that point, repeat (1) in step 2, replace point P retained by (1) in step 2 o (x, y, z) until backtracking to the absence of a AND P in the root node i (x, y, z) closer points;
(1) for original point cloud data set P ═ { P ═ P 1 ,P 2 ,…,P n Each data point P in i (x, y, z), performing m-neighborhood search using the k-d tree indexAsking, obtaining the nearest neighbor set Q ═ { P ═ P 1 ,P 2 ,…,P j …,P m }；
(2) Calculating P i Euclidean distance d from point (x, y, z) to each point in Q ij Then, as shown in formula I, the distance mean value to all neighborhood points is calculated
As shown in equation two:
(3) calculating point set P ═ P according to formula I and formula II 1 ,P 2 ,…,P n Mean of m neighborhoods of each point in the structure
(1) constructing neighborhood distance means of all points
The Gamma distribution probability density function of (1);
set of points P ═ P 1 ,P 2 ,…,P n Midpoint P i Performing gamma distribution fitting on the m neighborhood mean values of (x, y, z) to obtain gamma distribution parameter alpha 1 And beta 1 The probability distribution density function of the random variables of the Gamma distribution is shown as the formula (III):
in the formula of alpha 1 ,β 1 Delta is the shape parameter, proportion parameter and position parameter of gamma distribution, where delta is the lower limit of model distribution, the distance mean is greater than zero, and random variable is
Actually distributed on the positive half axis of the x axis, the real equivalent equation of the probability density function is given as follows:
(2) establishing all-point neighborhood mean distance
Cumulative probability of (F) 1 The functional expression of (a);
all point neighborhood mean distance
Cumulative probability of (F) 1 The function expression of (c) is represented by formula (v):
(3) the cumulative probability of selection reaches a certain threshold value
As a condition for preliminarily determining the noise point, the cumulative probability F is exceeded in the noise point determination process 1 Represents a query point P i The mean of m neighborhood points around (x, y, z) is large, i.e. larger than
Such a central query point may be initially identified as a noise point;
step 5, solving each query point P i M around (x, y, z)Distance change slope k of neighborhood points i Intercept b i And all point slopes k i ′：
(1) Euclidean distance d of m nearest neighbor points around query point ij Linear fitting is carried out, such as formula (I), each query point P is solved i (x, y, z) distance change slope k of m neighborhood points around i ：
y i ＝k i d ij +b i ⑥
(2) Calculating the point set P ═ P according to the formula 1 ,P 2 ,…,P n The slopes k of all points in } i ′；
Step 6, all the point slopes k i ' fitting gamma curve distribution:
(1) set of points P ═ P 1 ,P 2 ,…,P n Point P in i Slope k of (x, y, z) i ' fitting gamma curve distribution, as shown in formula (c), solving alpha of gamma distribution 2 And beta 2 ：
In a similar manner, wherein α 2 ,β 2 Respectively are a shape parameter and a proportion parameter of gamma distribution;
(2) establishing all-point neighborhood mean slope k i Cumulative probability of F 2 The function expression of (2), cumulative probability F 2 Is shown as (b):
(3) computing a query point P i Euclidean distance change slope k of neighborhood points around (x, y, z) i If the cumulative probability of (2) exceeds the cumulative probability F 2 Represents a query point P i The distance of m neighborhood points around (x, y, z) varies greatly, i.e. the slope is greater than a given threshold k' thre Indicating the neighborhood point distribution around the query pointVery uneven, such query points are considered noise points;
(1) because the two constraint conditions are independent of each other, the two constraint conditions are combined to construct a probability density distribution function, so that noise points can be better removed, and the method is based on
And k i The established two-dimensional Gamma joint probability density function is shown as the formula ninthly:
wherein alpha is 1 ,β 1 ,α 2 ,β 2 Are all greater than zero, are called
Obedience parameter is alpha 1 ,β 1 ,α 2 ,β 2 And (3) two-dimensional Gamma joint distribution of the distribution, the two-dimensional Gamma joint distribution function satisfies the formula (R):
in addition, the binary function in the distribution function satisfies the formula
(2) According to the property of two-dimensional Gamma joint distribution, setting an inner point which is a space region between a certain plane G and a curved surface f (x, y) as an inner point, and otherwise, setting an outer point which is a noise point;
and 8, removing all the noise points which do not meet the step (7), and obtaining a filtered point set.
Drawings
FIG. 1 raw point cloud data;
FIG. 2 shows a histogram of distance means and slope k, a fitting curve, and a joint probability density of the two, where the number of query point neighborhoods is 10;
FIG. 3 shows a histogram of distance means and slope k, a fitting curve, and a joint probability density of the two, where the number of query point neighborhoods is 30;
FIG. 4 shows a histogram of distance means and slope k, a fitted curve, and a joint probability density of the two, for a neighborhood number of query points of 40;
FIG. 5 is a histogram of distance means and slope k, a fitted curve, and a joint probability density of the two, for a neighborhood number of query points of 50;
FIG. 6 is a histogram of distance means and slope k, a fitted curve, and a joint probability density of the two, for a neighborhood number of query points of 60;
FIG. 7 illustrates the variation of each parameter of the mean value of the neighborhood distance for different numbers of neighborhoods;
FIG. 8 shows the variation of each parameter of the slope for different numbers of neighborhoods;
FIG. 9 shows the number of points and variance after filtering by the Gamma filter for different neighborhood numbers;
FIG. 10 is a graph of the results after filtering when m is 20;
FIG. 11 is a graph of the results after filtering when m is 30;
FIG. 12 is a graph of the results after filtering when m is 40;
FIG. 13 is a flow chart of a method.
Wherein, in fig. 1, diagrams (a) and (b) are respectively:
(a) an original point cloud data front view; (b) and (5) a side view of the original point cloud data.
In fig. 10, graphs (a) and (b) are respectively:
(a) when m is 20, the front view effect graph is obtained after filtering;
(b) when m is 20, the effect graph is viewed from the side after filtering;
in fig. 11, the graphs (a) and (b) are:
(a) when m is 30, the front view effect graph is obtained after filtering;
(b) when m is 30, filtering and then side-looking at the effect graph;
in fig. 12, the diagrams (a) and (b) are respectively:
(a) when m is 40, the front view effect graph is obtained after filtering;
(b) when m is 40, the effect graph is viewed from the side after filtering;
Detailed Description
In order to verify the effectiveness of the method, point cloud data on the front face of the administrative building of the university of science and engineering in Jiangxi are collected, as shown in fig. 1(a) and 1(b), 820025 data points are contained, the original point cloud contains elements such as buildings, cars, surrounding low trees and the like, the scene is complex, the types of ground objects are more, and the complex performance of the research scene effectively verifies the robustness of the filter constructed by the method. The original point cloud data is analyzed, and a large amount of non-ground and non-ground object point clouds are generated due to factors such as air dust and the like, and the noise points bring great interference to three-dimensional modeling of a geographic scene. The experimental software and hardware environment of the invention is as follows: intel (R) Xeon (R) CPU E3-1231v3 @3.40GHz and 8GB of memory; the Gamma filter experimental algorithm is realized by combining with PCL method library programming under MATLAB R2017 and Microsoft visual Studio 2017 compiling environment.
According to the method, the number of ideal neighborhood points (m) is determined according to multiple experimental results, when m is 10, 20, 30, 40, 50 and 60, experimental analysis is carried out, and meanwhile, the effect of the method is compared with the effect of a normal distribution experiment in the prior art IV. The distance mean value and neighborhood slope histogram, the Gamma function fitting curve, the normal distribution function fitting curve, and the Gamma combined probability density and the normal distribution combined probability density experimental effects of the distance mean value and the neighborhood slope histogram, the Gamma function fitting curve and the normal distribution function fitting curve under different neighborhood values of the algorithm and the normal distribution algorithm of the prior art are respectively shown in fig. 2 to fig. 6.
In order to better embody the denoising effect of the Gamma filter provided by the invention, the fitting curves of the distance mean value and the histogram of the slope k and the joint probability density of the distance mean value and the slope k are compared with the normal distribution filtering effect of the prior art IV, and when the cumulative probability of both the Gamma distribution curve and the normal distribution curve reaches 95%, the corresponding neighborhood distance mean value and slope value are shown in fig. 2 to fig. 6.
From fig. 2 to 6, when different numbers of neighborhoods are taken, threshold values of the neighborhood mean value and the slope of the query point change, when m takes 10, 30, 40, 50, and 60, values of the neighborhood mean value are 0.13573, 0.1979, 0.2962, 0.32793, and 0.36866, respectively, and values of the slope are 0.0038785, 0.0031125, 0.0024233, 0.0022197, and 0.0020282, respectively. When the neighborhood mean value and the slope histogram are respectively fitted, the change rule of the Gamma curve is consistent with that of the histogram, when the cumulative probability reaches 95%, the prior normal distribution is converged, the values of the horizontal axis of the Gamma filter are all larger than zero, and the fitting effect is more in line with the actual situation that the neighborhood mean value and the slope are larger than zero. As can be seen from fig. 7, when the same number of neighborhoods is selected, for the single constraint condition of the neighborhood mean, the number of noise filtering points when filtering is performed by using the Gamma curve is greater than the number of noise filtering points when filtering is performed by using the normal distribution curve, which indicates that the Gamma filter can filter more noise points under the same condition. As can be seen from fig. 8, when the same number of neighborhoods is selected, for the slope single constraint condition, the number of noise filtering points by filtering using the Gamma curve is greater than the number of noise filtering points by filtering using the normal distribution curve. Therefore, regardless of whether the single condition is the neighborhood mean or the slope around the query point, the effect of filtering by using the Gamma curve is better than that of filtering by using the normal distribution curve.
As can be seen from the fitting curves of the joint distribution in fig. 2 to fig. 6, when different numbers of neighborhoods are adopted, the Gamma joint probability density maps are all converged, the high information point regions are all in the central region, the distribution condition is more favorable for selecting a threshold value to perform filtering processing, the high information region of the two-dimensional joint normal distribution function is not converged in a certain central region, and when a threshold value of the neighborhood mean value or the slope is set, most information content cannot be contained. Meanwhile, as can be seen from fig. 2 to 6, when different numbers of neighborhoods are selected, the neighborhood mean value is calculatedWhen the sum slope is fitted by adopting two-dimensional combined Gamma distribution, the center of the curve is a region with more query points belonging to the mean value of a certain neighborhood and more query points belonging to a certain slope, and points outside the threshold range (namely the region in a red frame) are noise points. The neighborhood mean of the query point is too large, which means that the query point is far away from other points, and can be considered as a noise point, as shown in fig. 2, when m is 10, the query point is observed along the horizontal axis of the Gamma joint distribution, and when m is 10
When the query point is a noise point, the query point is considered as a noise point; the larger the slope of the query point is, the more uneven the distribution of neighborhood points around the query point is illustrated, as in fig. 2, when m is 10, the slope k 'of the query point is observed along the vertical axis of the Gamma joint distribution'>k′ thre When 0.0020282, the query point is considered to be a noise point.
When a Gamma filter is selected for noise point removal, different neighborhood point numbers are selected, the filtering effects are different, and in order to achieve a better filtering effect, the filtering effect with a better effect is selected from 6 different neighborhood point numbers. As can be seen from the relationship between the number of points after Gamma filtering and the number of neighborhood points in fig. 7, when the neighborhood mean is selected as the constraint condition, the number of noise points to be filtered is large when the number of neighborhood points is 30. As can be seen from the relationship between the number of points after Gamma filtering and the number of neighboring points in fig. 8, when the slope is selected as the constraint condition, the number of noise points to be filtered is large when m is 20, but the number of noise points to be filtered is not large when m is 30, and β is a constraint condition 2 In the distribution of (1), after m is 30, beta 2 The change is relatively gradual. As can be seen from fig. 9, when a two-dimensional joint Gamma probability density distribution is selected, that is, when two constraints work together, the number of noise points to be filtered is the largest when m is 30.
Therefore, in summary, when m is 30, the combined Gamma filter can effectively filter noise points under two constraints of neighborhood distance mean and slope. The experimental effect is shown in fig. 10, fig. 11, and fig. 12, and it can be seen from the diagrams that when m is 20, some noise points are not filtered (red frame), when m is 40, some details on the wall surface are filtered (blue frame), and when m is 30, the noise points are completely filtered, and the wall surface detail information can be better retained, and the filtering effect is the best.
The invention can be realized by any object-oriented programming language, and the realization scheme taking C + + language as an example is as follows:
(1) calling an open source Point Cloud Library (PCL) to create a Point Cloud data read-write class: PCReadClass and PCWriteClass;
(2) creating a point cloud k-d tree index creation class and a neighborhood point search class: kdTreeClass and NPFindClass;
(3) creating a query point and neighborhood point mean distance calculation class: avgditcalclass;
(4) creating a calculation class of the slope average value of the query point and the neighborhood point: avgckcalcclass;
(5) creating a calculation class of distance probability density Gamma distribution function of all point neighborhood mean values: avgdistrgammacls;
(6) creating a calculation class of the slope probability density Gamma distribution function of the neighborhood mean value of all the points: avgkrammass class;
(7) creating a two-dimensional Gamma distribution function calculation class of all point neighborhood mean distance and neighborhood mean slope: DistAndKGamaclass;
(8) calling AvgDistGamma class to judge whether the mean distance of the current query point is greater than a threshold value, and if so, storing the mean distance into a preliminary candidate point set Q; if the distance of the mean value of the current query point is smaller than the threshold value, calling AvgKGama class to judge whether the mean value of the neighborhood slope of the current query point is larger than a given threshold value, and if the mean value of the neighborhood slope of the current query point is larger than the threshold value, storing the mean value into a preliminary candidate point set Q;
(9) and calling DistAndKGamaclass to further delete the noise points of the candidate point set Q.
The above description is only for the preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (1)
1. A three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope is characterized by comprising the following steps:
step 1, constructing an index of an original point cloud k-d tree:
(1) reading original point cloud data set P ═ P 1 ,P 2 ,…,P n Storing the obtained data into a three-dimensional array, performing variance calculation on each dimension, selecting a dimension t (t is 1,2 and 3) with the largest variance, selecting a median mu of all values in the dimension t as a first comparison object, namely a dividing axis, obtaining two subsets, and storing the subsets into a tree node N (t, mu);
(2) repeating the process (1) in step 1 for the two subsets obtained until the original data set P ═ P 1 ,P 2 ,…,P n It cannot be subdivided; if a certain subset can not be divided in the process, storing the data in the subset to leaf nodes;
step 2, designing a nearest neighbor search algorithm of the k-d tree index:
(1) data point P to be inquired i (x, y, z) comparing with the mu value in each tree node N (t, mu) from the root node, and entering a right subtree query if the value in the dimension t is greater than mu; otherwise, entering left sub-tree query, and calculating P after reaching leaf node i (x, y, z) the distance d between the data points stored on this leaf node and the corresponding target point coordinates P are recorded o (x,y,z)；
(2) Backtracking, i.e. determining whether there is a distance P in an unaccessed data point i (x, y, z) is closer, and P is judged i Whether the distance between (x, y, z) and the unaccessed data point under the parent node to which the (x, y, z) belongs is larger than d, if so, the distance between the (x, y, z) and P does not exist in the branch i (x, y, z) closer points; otherwise, there is a separation of P i (x, y, z) closer point, then visit that point, repeat (1) in step 2, replace point P retained by (1) in step 2 o (x, y, z) until backtracking to the absence of a AND P in the root node i (x, y, z) closer points;
step 3, calculating the distance mean value between the query point and m neighborhood points around the query point and the global point neighborhood distance mean value:
(1) for original point cloud data set P ═ { P ═ P 1 ,P 2 ,…,P n Each data point P in i (x, y, z), performing m-neighborhood search by using the k-d tree index to obtain a nearest neighbor set Q ═ P 1 ,P 2 ,…,P j …,P m }；
(2) Calculating P i Euclidean distance d from point (x, y, z) to each point in Q ij Then, as shown in formula I, the distance average to all neighborhood points is calculated
As shown in equation two:
(3) calculating point set P ═ P according to formula I and formula II 1 ,P 2 ,…,P n Mean of m neighborhoods of each point in the structure
Step 4, neighborhood distance mean value of all points
Performing Gamma distribution fitting:
(1) constructing neighborhood distance means of all points
The Gamma distribution probability density function of (1);
set of points P ═ P 1 ,P 2 ,…,P n Midpoint P i Performing Gamma distribution fitting on the m neighborhood mean values of (x, y, z) to obtain a Gamma distribution parameter alpha 1 And beta 1 The probability distribution density function of the random variables of the Gamma distribution is shown as the formula (III):
in the formula of alpha 1 ,β 1 Delta is the shape parameter, proportion parameter and position parameter of Gamma distribution, where delta is the lower limit of model distribution, the average value is greater than zero, and random variable
Actually distributed on the positive half axis of the x axis, the real equivalent equation of the probability density function is given as follows:
(2) establishing all-point neighborhood mean distance
Cumulative probability of (F) 1 The functional expression of (a);
all point neighborhood mean distance
Cumulative probability of (F) 1 The function expression of (c) is represented by formula (v):
(3) the cumulative probability of selection reaches a certain threshold value
As a condition for preliminarily determining the noise point, inIn the noise point determination process, the cumulative probability F is exceeded 1 Represents a query point P i The mean of m neighborhood points around (x, y, z) is large, i.e. larger than
Such a central query point may be initially identified as a noise point;
step 5, solving each query point P i (x, y, z) distance change slope k of m neighborhood points around i Intercept b i And all point slopes k i '：
(1) Euclidean distance d of m nearest neighbor points around query point ij Linear fitting is carried out, such as formula (I), each query point P is solved i (x, y, z) distance change slope k of m neighborhood points around i ：
y i ＝k i d ij +b i ⑥
(2) Calculating the point set P ═ P according to the formula 1 ,P 2 ,…,P n The slopes k of all points in } i '
Step 6, all the point slopes k i ' Gamma curve distribution fitting:
(1) set of points P ═ P 1 ,P 2 ,…,P n Point P in i Slope k of (x, y, z) i Fitting Gamma curve distribution, as shown in formula (c), solving the alpha of Gamma distribution 2 And beta 2 ：
In a similar manner, wherein α 2 ,β 2 Respectively are a shape parameter and a proportion parameter of Gamma distribution;
(2) establishing all-point neighborhood mean slope k i ' cumulative probability F 2 The function expression of (2), cumulative probability F 2 Is shown as (b):
(3) computing a query point P i Euclidean distance change slope k of neighborhood points around (x, y, z) i If the cumulative probability of (D) exceeds the cumulative probability F 2 Represents a query point P i The distance of m neighborhood points around (x, y, z) varies greatly, i.e. the slope is greater than a given threshold k' thre Indicating that the neighborhood points around the query point are extremely unevenly distributed, and identifying the query point as a noise point;
step 7, establishing neighborhood mean distance of all points
And neighborhood mean distance slope k i ' two-dimensional Gamma joint probability density function:
(1) because the two constraint conditions are independent of each other, the two constraint conditions are combined to construct a probability density distribution function, so that noise points can be better removed, and the method is based on
And k i The established two-dimensional Gamma joint probability density function is shown as the formula ninthly:
wherein alpha is 1 ,β 1 ,α 2 ,β 2 Are all greater than zero, are called
Obedience parameter is alpha 1 ,β 1 ,α 2 ,β 2 And (3) two-dimensional Gamma joint distribution of the distribution, the two-dimensional Gamma joint distribution function satisfies the formula (R):
in addition theretoSatisfaction formula of binary function in distribution function
F(X,Y)＝P{X≤x,Y≤y} 
(2) Setting an inner point which is a spatial region between a certain plane G and a curved surface f (x, y) as an inner point according to the property of two-dimensional Gamma joint distribution, and otherwise, setting an outer point which is a noise point;
and 8, removing all the noise points which do not meet the step (7), and obtaining a filtered point set.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910365632.0A CN110246092B (en) | 2019-05-02 | 2019-05-02 | Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910365632.0A CN110246092B (en) | 2019-05-02 | 2019-05-02 | Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110246092A CN110246092A (en) | 2019-09-17 |
| CN110246092B true CN110246092B (en) | 2022-09-02 |
Family
ID=67883609
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910365632.0A Active CN110246092B (en) | 2019-05-02 | 2019-05-02 | Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110246092B (en) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111612806B (en) * | 2020-01-10 | 2023-07-28 | 江西理工大学 | Building facade window extraction method and device |
| CN111369610B (en) * | 2020-03-05 | 2022-09-06 | 山东交通学院 | Point cloud data gross error positioning and eliminating method based on credibility information |
| CN111524084B (en) * | 2020-05-13 | 2023-08-22 | 曹彬才 | Photon counting laser radar point cloud denoising method based on multi-peak Gaussian fitting |
| CN113031524B (en) * | 2021-02-07 | 2022-06-17 | 南京航空航天大学 | Cubic spline-based press fitting force envelope curve generation method |
| CN112986964B (en) * | 2021-02-26 | 2023-03-31 | 北京空间机电研究所 | Photon counting laser point cloud self-adaptive denoising method based on noise neighborhood density |
| CN113466826A (en) * | 2021-05-12 | 2021-10-01 | 武汉中仪物联技术股份有限公司 | Data denoising method, device, equipment and medium for range radar |
| CN114932562B (en) * | 2021-08-10 | 2024-04-19 | 南京航空航天大学 | Underground cable tunnel inspection robot based on laser radar and implementation method |
| CN116628429B (en) * | 2023-07-26 | 2023-10-10 | 青岛远度智能科技有限公司 | Intelligent control method for stable lifting of unmanned aerial vehicle |
| CN117828737B (en) * | 2024-01-04 | 2024-07-12 | 济南市园林绿化工程质量与安全中心 | Digital twin landscape design method |
| CN117974933B (en) * | 2024-03-28 | 2024-06-11 | 岐山县华强工贸有限责任公司 | 3D printing mold rapid scanning method for disc brake calipers |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102629367A (en) * | 2012-01-17 | 2012-08-08 | 安徽建筑工业学院 | Bilateral filtering de-noising method of point cloud data based on KDTree |
| CN104573705A (en) * | 2014-10-13 | 2015-04-29 | 北京建筑大学 | Clustering method for building laser scan point cloud data |
| CN106918813A (en) * | 2017-03-08 | 2017-07-04 | 浙江大学 | A kind of three-dimensional sonar point cloud chart image intensifying method based on distance statistics |
| CN107169464A (en) * | 2017-05-25 | 2017-09-15 | 中国农业科学院农业资源与农业区划研究所 | A kind of Method for Road Boundary Detection based on laser point cloud |
| CN107392875A (en) * | 2017-08-01 | 2017-11-24 | 长安大学 | A kind of cloud data denoising method based on the division of k neighbours domain |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107223268B (en) * | 2015-12-30 | 2020-08-07 | 中国科学院深圳先进技术研究院 | Three-dimensional point cloud model reconstruction method and device |
-
2019
- 2019-05-02 CN CN201910365632.0A patent/CN110246092B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102629367A (en) * | 2012-01-17 | 2012-08-08 | 安徽建筑工业学院 | Bilateral filtering de-noising method of point cloud data based on KDTree |
| CN104573705A (en) * | 2014-10-13 | 2015-04-29 | 北京建筑大学 | Clustering method for building laser scan point cloud data |
| CN106918813A (en) * | 2017-03-08 | 2017-07-04 | 浙江大学 | A kind of three-dimensional sonar point cloud chart image intensifying method based on distance statistics |
| CN107169464A (en) * | 2017-05-25 | 2017-09-15 | 中国农业科学院农业资源与农业区划研究所 | A kind of Method for Road Boundary Detection based on laser point cloud |
| CN107392875A (en) * | 2017-08-01 | 2017-11-24 | 长安大学 | A kind of cloud data denoising method based on the division of k neighbours domain |
Non-Patent Citations (3)
| Title |
|---|
| local 3D point cloud denoising via bipartite graph approximation & totalvariation;Chinthaka Dinesh 等;《2018 IEEE 20th international workshop on multimedia signal processing》;20181129;全文 * |
| 基于频率直方图的K邻域稀疏离群点移除算法;郭子选等;《计算机应用与软件》;20161215(第12期);147-153页 * |
| 采用密度k-means和改进双边滤波的点云自适应去噪算法;郭进等;《传感器与微系统》;20160720(第07期);169-172页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110246092A (en) | 2019-09-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110246092B (en) | Three-dimensional laser point cloud denoising method considering neighborhood point mean distance and slope | |
| CN105844602B (en) | A kind of airborne LIDAR point cloud three-dimensional filtering method based on volume elements | |
| CN111598780B (en) | Terrain adaptive interpolation filtering method suitable for airborne LiDAR point cloud | |
| CN108961294B (en) | Three-dimensional point cloud segmentation method and device | |
| CN107657659A (en) | The Manhattan construction method for automatic modeling of scanning three-dimensional point cloud is fitted based on cuboid | |
| Pound et al. | A patch-based approach to 3D plant shoot phenotyping | |
| CN115222625A (en) | Laser radar point cloud denoising method based on multi-scale noise | |
| US20150178594A1 (en) | Point cloud simplification | |
| Shen et al. | A Kd-tree-based outlier detection method for airborne LiDAR point clouds | |
| CN114332291B (en) | Method for extracting outline rule of oblique photography model building | |
| CN114266987A (en) | Intelligent identification method for high slope dangerous rock mass of unmanned aerial vehicle | |
| CN110910435B (en) | Building point cloud extraction method and device, computer equipment and readable storage medium | |
| CN111275616B (en) | Low-altitude aerial image splicing method and device | |
| Zhu et al. | 3D reconstruction of plant leaves for high-throughput phenotyping | |
| CN111462017A (en) | Denoising method for tunnel laser point cloud data | |
| CN111524127A (en) | Urban road surface extraction method for low-altitude airborne laser radar data | |
| CN114332134A (en) | Building facade extraction method and device based on dense point cloud | |
| Meyer et al. | CherryPicker: Semantic skeletonization and topological reconstruction of cherry trees | |
| CN118052938A (en) | Building multi-detail level model reconstruction method based on multi-source data fusion | |
| CN111986223B (en) | Method for extracting trees in outdoor point cloud scene based on energy function | |
| CN116186864B (en) | Deep foundation pit model rapid modeling method and system based on BIM technology | |
| CN112884884A (en) | Candidate region generation method and system | |
| CN110222742B (en) | Point cloud segmentation method, device, storage medium and equipment based on layered multi-echo | |
| CN109118565B (en) | Electric power corridor three-dimensional model texture mapping method considering shielding of pole tower power line | |
| CN113655498B (en) | Method and system for extracting cone barrel information in racetrack based on laser radar |
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 |