CN113205529A - Method for segmenting top surface of building based on airborne LiDAR point cloud - Google Patents

Method for segmenting top surface of building based on airborne LiDAR point cloud Download PDF

Info

Publication number
CN113205529A
CN113205529A CN202110418234.8A CN202110418234A CN113205529A CN 113205529 A CN113205529 A CN 113205529A CN 202110418234 A CN202110418234 A CN 202110418234A CN 113205529 A CN113205529 A CN 113205529A
Authority
CN
China
Prior art keywords
plane
points
point
segmentation
building
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.)
Granted
Application number
CN202110418234.8A
Other languages
Chinese (zh)
Other versions
CN113205529B (en
Inventor
季顺平
王瑄
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202110418234.8A priority Critical patent/CN113205529B/en
Publication of CN113205529A publication Critical patent/CN113205529A/en
Application granted granted Critical
Publication of CN113205529B publication Critical patent/CN113205529B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/66—Analysis of geometric attributes of image moments or centre of gravity
    • 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/20172—Image enhancement details
    • G06T2207/20192—Edge enhancement; Edge preservation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Geometry (AREA)
  • Image Analysis (AREA)

Abstract

The invention provides a method for partitioning the top surface of a building based on airborne LiDAR point cloud, which well realizes high-precision extraction of the top surface of the building and has strong robustness for different data types. Firstly, carrying out European clustering on original building point clouds to segment the point clouds of each building; and then performing multi-scale segmentation based on L0 gradient minimization on each building respectively to obtain an initial plane segmentation result. Specific to the extraction work of the top surface of the house, the L0 gradient minimization algorithm is improved in a targeted mode, the application of strategies comprising sequencing and region expansion and additional geometric constraints are included, and the accuracy and the robustness of the L0 gradient minimization algorithm in plane segmentation are improved well. And finally, refining the initial result by adopting a graph cut-based post-processing optimization algorithm, and solving the problems of excessive segmentation and jagged edges possibly existing in the initial result to obtain a final house top surface extraction result.

Description

Method for segmenting top surface of building based on airborne LiDAR point cloud
Technical Field
The invention relates to a method for extracting a house top surface based on airborne LiDAR point cloud by utilizing L0 gradient minimization and graph cut global optimization, which can be used for LiDAR point cloud classification segmentation, three-dimensional building model reconstruction and the like.
Background
Three-dimensional reconstruction of buildings is an important research topic in the fields of computer vision, photogrammetry, remote sensing and the like. With the rapid development of laser Detection and Ranging (Light Detection and Ranging) systems, it becomes possible to rapidly and accurately obtain large-scale urban three-dimensional point clouds, which brings great convenience to the reconstruction work of three-dimensional buildings. Existing building three-dimensional modeling methods based on airborne LiDAR point clouds can be roughly divided into two categories: a model driven method and a data driven method. The model-driven approach is a top-down approach and the building models are considered to be present in a pre-defined model library. The data driving method is to model the building from bottom to top, firstly, the top surface of the building needs to be extracted from the point cloud, and then the topological structure of the top surface of the building is analyzed to finally construct a complete building model. The data driving method starting from data is more effective because the model driving method cannot accurately restore the house model structure which does not exist in the model base. For the data-driven method, how to accurately and efficiently extract the top surface of the building from the point cloud becomes the key of the three-dimensional modeling work. The method aims at the problem of building top surface extraction, related research layers at home and abroad do not exist at present, but due to different data point cloud densities and different house complexity degrees, stable and high-precision house top surface extraction results cannot be obtained by the existing method, and therefore a high-precision and robust house top surface extraction algorithm needs to be designed to further promote building model reconstruction work.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides the house top surface extraction method based on the L0 gradient minimization algorithm and the graph cut optimization algorithm, and the method can obtain the high-precision house top surface extraction result when facing different types of data.
The technical scheme of the invention is a method for partitioning the top surface of a building based on airborne LiDAR point cloud, which comprises the following steps:
step 1, carrying out European clustering on original building point clouds to segment the point clouds of each building;
step 2, respectively carrying out multi-scale segmentation of each building based on an L0 gradient minimization algorithm of region expansion to obtain an initial plane segmentation result;
and 3, refining the initial plane segmentation result by adopting a graph-segmentation-based post-processing optimization algorithm, and extracting the top surface of the high-precision building.
Further, the region expansion-based L0 gradient minimization algorithm in step 2 is an improvement on the L0 gradient minimization algorithm, and the specific principle is as follows;
the energy function of the L0 gradient minimization algorithm based on regional dilation is defined as:
Figure BDA0003026826680000021
wherein M represents the length of the input point cloud, IiAnd SiRespectively represent ith point p in point cloudiInput and output normal vectors; n is a radical ofiRepresents piThe k nearest neighbor neighborhood points or the neighborhood points in a sphere with R as the radius in the three-dimensional space, i and j respectively represent the indexes of the neighborhood points, and the parameter lambda is the weight of the L0 norm;
the equation (2) is solved approximately by first solving for a pair of adjacent points piAnd pjFor example, their contribution F to the overall energy function F is:
Figure BDA0003026826680000022
to minimize f, two different strategies are adopted: non-fusion of piAnd pjAt this time Si≠SjOr fusion of piAnd pjAt this time Si＝SjAnd the L0 norm between two neighboring points is eliminated;
for Si≠SjThe specific method comprises the following steps:
Figure BDA0003026826680000023
for Si＝SjThe specific method comprises the following steps:
Figure BDA0003026826680000024
combining the two conditions to obtain:
Figure BDA0003026826680000025
geometric constraints are added to further improve the stability of the algorithm, and a plane distance concept is introduced firstly, so that a pair of adjacent points piAnd pjFor example, the planar distance Pd between themi,jThe definition is as follows:
Pdi,j＝abs[pij·Ii] (6)
wherein P isijRepresents from piTo pjVector of (a), IiRepresents piThe normal vector of (2), abs]Taking an absolute value; thus Pdi,jRepresents two points in a normal vector IiProjection distance in direction, when the plane distance between two points is greater than a prescribed threshold value TpdNon-zero gradients between them are no longer penalized;
equation 6 is therefore rewritten as:
Figure BDA0003026826680000031
next, the above formula is extended from neighboring points to a set of neighboring points, in GiRepresents a set of points, YiRepresenting its normal vector, the number of points in the set of points is wiSet of points GiAnd GjThe number of the connecting points between is ci,jCombining the energy formula between the adjacent point sets of the geometric constraint condition as follows:
Figure BDA0003026826680000032
the above formula is solved as:
Figure BDA0003026826680000033
equation 10 is the expansion condition of the L0 gradient minimization algorithm based on the region expansion.
Further, the flow of the region expansion based L0 gradient minimization algorithm of step 2 is as follows,
1) calculating local smooth characteristics of each point in the point cloud to be segmented, and sequencing the operation sequence of the points according to the local smooth characteristic values;
2) initializing a point set and setting a preset parameter beta, wherein each point is regarded as a point set with only one point;
3) introducing an auxiliary parameter lambda, and increasing the value of lambda by k times in each iteration; the initialization of the auxiliary parameters is accomplished by the following formula:
Figure BDA0003026826680000034
wherein mean { } represents taking the median;
4) in each iteration process, all point sets are processed respectively in sequence, and for the point set GiDetermining whether to use the point set G according to the expansion conditioniAll neighborhood point sets GjIs added to GiIf G isjAdding GiCorresponding to wiWill be updated, GjWill be added to NiAfter completion of one expansion, GjThe related information is deleted, and the number of the point sets to be processed is reduced by 1;
5) and circularly processing all the point sets until the auxiliary parameter lambda is greater than the preset parameter beta, terminating iteration, enabling the rest point sets to correspond to the segmented planes, and enabling the points in each point set to have consistent plane normal vectors.
Further, the local smooth feature is the maximum value of the angle difference value between the point and the normal vector of other neighborhood points in the neighborhood.
Furthermore, the multi-scale segmentation in step 2 comprises the following steps,
(1) selecting an initial segmentation scale parameter beta, and taking point cloud clustering obtained by Euclidean clustering as input, wherein the segmentation scale parameter beta is the same as a preset parameter;
(2) performing plane verification on each cluster, if the plane index is greater than a threshold value, directly outputting the plane as a legal plane, otherwise, performing plane segmentation under the current segmentation scale through an L0 gradient minimization algorithm based on region expansion to obtain a new cluster;
the plane index is defined as: for a cluster, fitting a plane by PCA (principal component analysis), and calculating the distance d from each point in the cluster to the fitted planeiThe plane index is the average distance from all points to the fitting plane;
(3) when all the clusters are processed in the step (2), if a new cluster is generated, the new generated cluster is used as input to carry out next iteration, and the segmentation scale parameter is reduced; if no new clusters are generated, i.e., all input clusters are legal planes, the algorithm terminates.
Further, in the post-processing flow based on graph cut in step 3, the multi-scale division result is used as an initial plane, and an objective function is optimized through a graph cut algorithm to obtain an optimal corresponding relation from each point to the plane, wherein the optimized objective function is composed of three parts, which are respectively: data items, smoothing items, and tag items, defined as follows:
Figure BDA0003026826680000041
where L is a plane, P is a set of points for all points on the plane L, LpI.e. a plane to which the representative point p belongs, wherein p, q epsilon and N in the formula represent two adjacent points;
wherein the data item is the distance from the measuring point to the plane, and the calculation formula is as follows:
Figure BDA0003026826680000042
wherein
Figure BDA0003026826680000043
Is LpThe plane equation of (a);
the smooth term is used to penalize the number of adjacent points belonging to different planes, δ () is an indicator function, when L isp≠LqThe function value is 1, and ω is a weight;
the last item is a label item, which is used for punishing the number of different planes corresponding to all the points finally, and k is a weight.
Further, the energy function minimization problem in equation 12 is solved by α -expansion.
The invention provides an algorithm for extracting the top surface of a house by utilizing L0 gradient minimization and graph cut global optimization based on airborne LiDAR point cloud, which well realizes high-precision extraction of the top surface of the house and has strong robustness for different data types. Firstly, carrying out European clustering on original building point clouds to segment the point clouds of each building; and then performing multi-scale segmentation based on L0 gradient minimization on each building respectively to obtain an initial plane segmentation result. Specific to the extraction work of the top surface of the house, the L0 gradient minimization algorithm is improved in a targeted mode, the application of strategies comprising sequencing and region expansion and additional geometric constraints are included, and the accuracy and the robustness of the L0 gradient minimization algorithm in plane segmentation are improved well. And finally, refining the initial result by adopting a graph cut-based post-processing optimization algorithm, and solving the problems of excessive segmentation and jagged edges possibly existing in the initial result to obtain a final house top surface extraction result.
Drawings
FIG. 1 is an overall algorithm flow diagram of the present invention.
FIG. 2 is a flow chart of the multi-scale segmentation algorithm of the present invention.
Detailed Description
The technical solution of the present invention is further described below by way of examples with reference to the accompanying drawings.
The technical scheme adopted by the invention is as follows: firstly, carrying out European clustering on original building point clouds to segment the point clouds of each building; then, performing multi-scale segmentation based on L0 gradient minimization on each building to obtain an initial plane segmentation result; and finally, refining the initial result by adopting a graph cut-based post-processing optimization algorithm, and solving the problems of excessive segmentation and jagged edges possibly existing in the initial result to obtain the final building top surface extraction result.
The Euclidean clustering method is a common method in point cloud example segmentation, and the point cloud of each building can be obtained through Euclidean clustering, so that the time efficiency of the algorithm is improved, and the stability of the algorithm is improved.
The above-described L0 gradient minimization algorithm is the fundamental method of our multi-scale segmentation, which was originally applied in signal processing and used to control the number of non-zero gradients between adjacent signals. The invention extends an L0 gradient minimization algorithm to a plane segmentation algorithm according to the aim of the invention, and provides an L0 gradient minimization algorithm based on region expansion.
The energy function of the L0 gradient minimization algorithm is defined as:
Figure BDA0003026826680000051
i denotes the input signal and S denotes the output signal. S represents the gradient of the output signal. I | · | purple wind0Representing the L0 norm. The parameter λ controls the weighting of the L0 norm to determine the degree of smoothing of the final output signal.
After being expanded into the point cloud and discretized, equation 1 can be expressed as:
Figure BDA0003026826680000061
wherein M represents the length of the input point cloud, IiAnd SiRespectively represent ith point p in point cloudiInput and output normal vectors. N is a radical ofiRepresents piThe neighborhood of (a) may be k nearest neighbor neighborhood points or neighborhood points in a sphere with a radius of R in a three-dimensional space, i and j respectively represent indexes of the neighborhood points, and the parameter λ is a weight of L0 norm.
Since directly solving the minimum value of equation 2 is an NP-hard problem, it is necessary to solve this problem approximately.
First with a pair of adjacent points piAnd pjFor example, their contribution F to the overall energy function F is:
Figure BDA0003026826680000062
to minimize f, two different strategies may be adopted: non-fusion of piAnd pjAt this time Si≠SjOr fusion of piAnd pjAt this time Si＝SjAnd the L0 norm between two neighboring points is eliminated.
For Si≠SjThe specific method comprises the following steps:
Figure BDA0003026826680000063
for Si＝SjThe specific method comprises the following steps:
Figure BDA0003026826680000064
combining the two cases we can get:
Figure BDA0003026826680000065
meanwhile, in consideration of the particularity of the three-dimensional plane segmentation problem, geometric constraints are added to further improve the stability of the algorithm. The concept of a plane distance is introduced here, asA pair of adjacent points piAnd pjFor example, the planar distance Pd between themi,jThe following can be defined:
Pdi,j＝abs[pij·Ii] (19)
wherein P isijRepresents from piTo pjVector of (a), IiRepresents piThe normal vector of (2), abs]The absolute value is taken. Thus Pdi,jIt represents two points in the normal vector IiThe projected distance in the direction, when this projected distance is too large, we consider that they do not belong to the same plane, although they may have similar normal vectors. So that when the planar distance between two points is greater than a prescribed threshold value TpdWe no longer penalize the non-zero gradients between them.
Equation 6 can therefore be rewritten as:
Figure BDA0003026826680000071
next, the above formula is extended from neighboring points to a set of neighboring points. With GiRepresents a set of points, YiRepresenting its normal vector, the number of points in the set of points is wiSet of points GiAnd GjThe number of the connecting points between is ci,j(ii) a The energy formula between the adjacent point sets in combination with the above geometric constraint condition is:
Figure BDA0003026826680000072
the plane distance is also calculated for the set of neighboring points by equation 7, where PijRepresentative of a set of points GiAnd GjConnects the formed vectors.
The above formula is solved as:
Figure BDA0003026826680000073
equation 10 is the dilation condition for our region dilation based L0 gradient minimization algorithm. We now describe the overall flow of the L0 gradient minimization algorithm.
1) Calculating local smooth characteristics of each point in the point cloud to be segmented, and sequencing the operation sequence of the points according to local smooth characteristic values, wherein the local smooth characteristics are the maximum values of the angle difference values of the normal vectors of the point and other neighborhood points in the neighborhood of the point;
2) initializing a point set and setting a preset parameter beta, wherein each point is regarded as a point set with only one point;
3) introducing an auxiliary parameter lambda, and increasing the value of lambda by 1.2 times in each iteration; the initialization of the auxiliary parameters is accomplished by the following formula:
Figure BDA0003026826680000081
where mean { } represents taking the median.
4) In each iteration process, all point sets are processed respectively in sequence, and for the point set GiWhether to set the points G is determined according to the expansion condition (equation 10)iAll neighborhood point sets GjIs added to GiIf G isjAdding GiCorresponding to wiWill be updated, GjWill be added to NiAfter completion of one expansion, GjThe related information is deleted, and the number of the point sets to be processed is reduced by 1;
5) and circularly processing all the point sets until the auxiliary parameter lambda is greater than the preset parameter beta, terminating iteration, enabling the rest point sets to correspond to the segmented planes, and enabling the points in each point set to have consistent plane normal vectors.
The algorithm flow is as follows:
Figure BDA0003026826680000082
Figure BDA0003026826680000091
because the segmentation effect of the single L0 gradient minimization algorithm is greatly influenced by the parameter lambda, in order to enable the algorithm to have self-adaptability and higher robustness, the invention provides a multi-scale segmentation process, firstly, a plane index concept is introduced, and a plane index c is used for measuring whether a cluster is a legal plane; for a cluster, fitting a plane by Principal Component Analysis (PCA), and calculating the distance d from each point in the cluster to the fitting planeiThe plane index is the average distance from all points to the fitting plane; only if the plane index of one cluster is greater than a preset threshold value TcWe will consider it to be a legal plane.
1) And performing Euclidean clustering to obtain point cloud clusters of different houses.
2) Selecting a larger segmentation scale parameter beta as an initial beta, and using point cloud clustering obtained by Euclidean clustering as input.
3) And firstly carrying out plane verification on each cluster, if the plane index is larger than a threshold value, directly outputting the plane as a legal plane, and otherwise, carrying out plane segmentation through an L0 gradient minimization algorithm based on region expansion under the current segmentation scale to obtain a new cluster.
4) And (3) when all the clusters are processed in the step (3), if a new cluster is generated, the algorithm takes the new generated cluster as input to carry out the next iteration, the segmentation scale parameter is reduced, and the parameter beta is reduced by 1.5 times. If no new clusters are generated, i.e., all input clusters are legal planes, the algorithm terminates.
The multi-scale algorithm is robust to the starting lambda and can detect as many potential planes as possible. In order to refine the segmentation result obtained by multi-scale segmentation, a post-processing flow based on graph segmentation is adopted to solve the problems of excessive segmentation and jagged edges in the multi-scale segmentation.
The post-processing flow based on graph cut mainly uses a multi-scale segmentation result as an initial plane to optimize an objective function through a graph cut algorithm so as to obtain the optimal corresponding relation from each point to the plane. The optimized objective function consists of three parts, namely: data items, smoothing items, and tag items. The definition is as follows:
Figure BDA0003026826680000101
where L is a plane, P is a set of points for all points on the plane L, LpI.e. a plane to which a representative point p belongs, wherein p, q ∈ N in the formula represents two adjacent points;
the data item mainly measures the distance from a point to a plane, and the calculation formula is as follows:
Figure BDA0003026826680000102
wherein
Figure BDA0003026826680000103
Is LpThe plane equation of (c).
The smooth term is used to penalize the number of adjacent points belonging to different planes, δ () is an indicator function, when L isp≠LqThe function is 1, and ω is the weight of this cost.
The last item is a label item, which mainly penalizes the number of different planes corresponding to all final points, and if the weight k is set to be smaller, the number of the final planes is more. It is noted that the resulting planes here are a subset of the input planes.
The energy function minimization problem in equation 12 can be solved by alpha-expansion.
After the post-processing flow is finished, a refined plane segmentation result is finally obtained.
To verify the role of the present invention in building roof extraction, we evaluated the method of the present invention on two sets of data, one from Hangzhou City, China and one from Vaihingen, Germany. The Vaihingen dataset has three zones, and the following table details the details of the two sets of datasets.
Table 1 details of the two sets of data sets
Figure BDA0003026826680000111
The main evaluation indexes are as follows: integrity (C)m) Accuracy (correction, C)r) Quality of split (Q)l) Over-segmentation index (R)o) An under-segmentation index (R)u) Edge precision (B)p) Edge recall (B)r) Edge F-measure (F-measure)m). These eight indices are evaluation indices commonly used in plane segmentation.
The integrity reflects the percentage of correctly extracted planes in the truth value, the accuracy reflects the percentage of correctly extracted planes in the extracted planes, and the segmentation quality is a comprehensive index and can be regarded as a main index. The over-segmentation index and the under-segmentation index are used to measure the degree of over-segmentation and under-segmentation of the extracted result. The edge F index reflects the edge accuracy of the extraction plane, as well as the composite index of the edge recall.
The results of the experiment are shown in the following table:
TABLE 2 results of the experiment
Figure BDA0003026826680000112
Figure BDA0003026826680000121
Example (b):
referring to fig. 1 and 2, the present invention extracts building roof surfaces from airborne LiDAR point clouds with high efficiency and high accuracy through a multi-scale L0 gradient minimization algorithm and a graph cut optimization algorithm.
Take the measurement area number one of the published data set Vaihingen of the ISPRS as an example.
The method comprises the following steps: european clustering
First we calculate the resolution of the point cloud by counting the average distance of all points from the nearest point. And setting a minimum distance threshold of Euclidean clustering according to the point cloud resolution, and proving that the Euclidean clustering segmentation effect can be well obtained by triple point cloud resolution through experiments. The specific flow of Euclidean clustering is as follows:
establishing a Kd-tree according to the minimum distance threshold, and determining the neighborhood points of each point through the Kd-tree, wherein the distances from the neighborhood points to the central point are not more than the distance threshold; then, starting from a point randomly, forming an initial cluster by the point, adding the neighborhood point of the point into the cluster according to a Kd-tree, then, for each newly added point, further adding the neighborhood points into the cluster through the operation, repeating the operation until no new point is added, obtaining a cluster at the moment, removing the point forming the cluster from the original point cloud, and repeating the process from the random point selection until all the points are added into a cluster.
Step two: multi-scale segmentation
Referring to fig. 2, a clustering result of the euclidean clustering is used as an input, parameter setting is performed according to data characteristics, experiments prove that a multi-scale segmentation process is insensitive to an initial lambda, and the initial lambda is set to be 0.1 in order to reduce unnecessary iteration. Then determining a threshold T for plane verificationcExperiment proves TcA selection of 0.05 gives good results. For the L0 gradient minimization algorithm, a plane distance threshold T is setpdIs 0.1. After all parameter settings are determined, we perform multi-scale segmentation according to the flow described in fig. 2 to obtain an initial plane segmentation result.
Step three: based on graph cut post-processing
Firstly, determining the parameter settings of graph cut optimization, namely the weights of a smooth item and a label item. For the data characteristics, we set the weight of the smoothing term to be 0.05 and the weight of the tag term to be 5. And then obtaining initial labels according to the initial plane segmentation result, wherein each plane model is used as one label, and one plane label comprises the gravity center and a plane normal vector. And then, solving the global optimization problem by utilizing an alpha-expansion algorithm to obtain labels corresponding to all points, determining a point set corresponding to each surface according to the labels, and finally obtaining the extraction result of the top surface of the house.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.

Claims (7)

1. A method for segmenting a building roof based on airborne LiDAR point clouds, comprising the steps of:
step 1, carrying out European clustering on original building point clouds to segment the point clouds of each building;
step 2, respectively carrying out multi-scale segmentation of each building based on an L0 gradient minimization algorithm of region expansion to obtain an initial plane segmentation result;
and 3, refining the initial plane segmentation result by adopting a graph-segmentation-based post-processing optimization algorithm, and extracting the top surface of the high-precision building.
2. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 1, wherein: the region expansion-based L0 gradient minimization algorithm in the step 2 is an improvement on the L0 gradient minimization algorithm, and the specific principle is as follows;
the energy function of the L0 gradient minimization algorithm based on regional dilation is defined as:
Figure FDA0003026826670000011
wherein M represents the length of the input point cloud, IiAnd SiRespectively represent ith point p in point cloudiInput and output normal vectors; n is a radical ofiRepresents piThe k nearest neighbor neighborhood points or the neighborhood points in a sphere with R as the radius in the three-dimensional space, i and j respectively represent the indexes of the neighborhood points, and the parameter lambda is the weight of the L0 norm;
the equation (2) is solved approximately by first solving for a pair of adjacent points piAnd pjFor example, their contribution F to the overall energy function F is:
Figure FDA0003026826670000012
to minimize f, two different strategies are adopted: non-fusion of piAnd pjAt this time Si≠SjOr fusion of piAnd pjAt this time Si＝SjAnd the L0 norm between two neighboring points is eliminated;
for Si≠SjThe specific method comprises the following steps:
Figure FDA0003026826670000013
for Si＝SjThe specific method comprises the following steps:
Figure FDA0003026826670000021
combining the two conditions to obtain:
Figure FDA0003026826670000022
geometric constraints are added to further improve the stability of the algorithm, and a plane distance concept is introduced firstly, so that a pair of adjacent points piAnd pjFor example, the planar distance Pd between themi,jThe definition is as follows:
Pdi,j＝abs[pij·Ii] (6)
wherein P isijRepresents from piTo pjVector of (a), IiRepresents piThe normal vector of (2), abs]Taking an absolute value; thus Pdi,jRepresents two points in a normal vector IiProjection distance in direction, when the plane distance between two points is greater than a prescribed threshold value TpdNon-zero gradients between them are no longer penalized;
equation 6 is therefore rewritten as:
Figure FDA0003026826670000023
next, the above formula is extended from neighboring points to a set of neighboring points, in GiRepresents a set of points, YiRepresenting its normal vector, the number of points in the set of points is wiSet of points GiAnd GjThe number of the connecting points between is ci,jCombining the energy formula between the adjacent point sets of the geometric constraint condition as follows:
Figure FDA0003026826670000024
the above formula is solved as:
Figure FDA0003026826670000025
equation 10 is the expansion condition of the L0 gradient minimization algorithm based on the region expansion.
3. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 2, wherein: the flow of the region expansion based L0 gradient minimization algorithm of step 2 is as follows,
1) calculating local smooth characteristics of each point in the point cloud to be segmented, and sequencing the operation sequence of the points according to the local smooth characteristic values;
2) initializing a point set and setting a preset parameter beta, wherein each point is regarded as a point set with only one point;
3) introducing an auxiliary parameter lambda, and increasing the value of lambda by k times in each iteration; the initialization of the auxiliary parameters is accomplished by the following formula:
Figure FDA0003026826670000031
wherein mean { } represents taking the median;
4) in each iteration process, all point sets are processed respectively in sequence, and for the point set GiDetermining whether to use the point set G according to the expansion conditioniAll neighborhood point sets GjIs added to GiIf G isjAdding GiCorresponding to wiWill be updated, GjWill be added to NiAfter completion of one expansion, GjThe related information is deleted, and the number of the point sets to be processed is reduced by 1;
5) and circularly processing all the point sets until the auxiliary parameter lambda is greater than the preset parameter beta, terminating iteration, enabling the rest point sets to correspond to the segmented planes, and enabling the points in each point set to have consistent plane normal vectors.
4. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 3, wherein: the local smooth feature is the maximum value of the angle difference value between the point and the normal vector of other neighborhood points in the neighborhood.
5. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 3, wherein: the flow of the multi-scale segmentation in the step 2 is that,
(1) selecting an initial segmentation scale parameter beta, and taking point cloud clustering obtained by Euclidean clustering as input, wherein the segmentation scale parameter beta is the same as a preset parameter;
(2) performing plane verification on each cluster, if the plane index is greater than a threshold value, directly outputting the plane as a legal plane, otherwise, performing plane segmentation under the current segmentation scale through an L0 gradient minimization algorithm based on region expansion to obtain a new cluster;
the plane index is defined as: for a cluster, fitting a plane by PCA (principal component analysis), and calculating the distance d from each point in the cluster to the fitted planeiThe plane index is the average distance from all points to the fitting plane;
(3) when all the clusters are processed in the step (2), if a new cluster is generated, the new generated cluster is used as input to carry out next iteration, and the segmentation scale parameter is reduced; if no new clusters are generated, i.e., all input clusters are legal planes, the algorithm terminates.
6. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 1, wherein: in the post-processing flow based on graph cut in the step 3, the multi-scale segmentation result is used as an initial plane, and an objective function is optimized through a graph cut algorithm to obtain the optimal corresponding relation from each point to the plane, wherein the optimized objective function consists of three parts, namely: data items, smoothing items, and tag items, defined as follows:
Figure FDA0003026826670000041
where L is a plane, P is a set of points for all points on the plane L, LpI.e. a plane to which the representative point p belongs, wherein p, q epsilon and N in the formula represent two adjacent points;
wherein the data item is the distance from the measuring point to the plane, and the calculation formula is as follows:
Figure FDA0003026826670000042
wherein
Figure FDA0003026826670000043
Is LpThe plane equation of (a);
the smooth term is used to penalize the number of adjacent points belonging to different planes, δ () is an indicator function, when L isp≠LqThe function value is 1, and ω is a weight;
the last item is a label item, which is used for punishing the number of different planes corresponding to all the points finally, and k is a weight.
7. The method of segmenting a building roof based on airborne LiDAR point cloud of claim 6, wherein: the energy function minimization problem in equation 12 is solved by alpha-expansion.
CN202110418234.8A 2021-04-19 2021-04-19 Method for segmenting top surface of building based on airborne LiDAR point cloud Active CN113205529B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110418234.8A CN113205529B (en) 2021-04-19 2021-04-19 Method for segmenting top surface of building based on airborne LiDAR point cloud

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110418234.8A CN113205529B (en) 2021-04-19 2021-04-19 Method for segmenting top surface of building based on airborne LiDAR point cloud

Publications (2)

Publication Number Publication Date
CN113205529A true CN113205529A (en) 2021-08-03
CN113205529B CN113205529B (en) 2022-04-29

Family

ID=77027371

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110418234.8A Active CN113205529B (en) 2021-04-19 2021-04-19 Method for segmenting top surface of building based on airborne LiDAR point cloud

Country Status (1)

Country Link
CN (1) CN113205529B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113822999A (en) * 2021-09-22 2021-12-21 杭州趣村游文旅集团有限公司 Building segmentation method of digital country three-dimensional model
CN116071530A (en) * 2023-03-07 2023-05-05 青岛市勘察测绘研究院 Building roof voxelized segmentation method based on airborne laser point cloud
CN117576144A (en) * 2024-01-15 2024-02-20 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090310867A1 (en) * 2008-06-12 2009-12-17 Bogdan Calin Mihai Matei Building segmentation for densely built urban regions using aerial lidar data
CN105139379A (en) * 2015-07-30 2015-12-09 滁州学院 Airborne Lidar point cloud building top surface gradual extraction method based on classifying and laying
CN106815841A (en) * 2017-01-13 2017-06-09 合肥师范学院 Image object partitioning algorithm based on T junction point clue
CN107944383A (en) * 2017-11-21 2018-04-20 航天科工智慧产业发展有限公司 Building roof patch division method based on three-dimensional Voronoi diagram
WO2020072865A1 (en) * 2018-10-05 2020-04-09 Interdigital Vc Holdings, Inc. A method and device for encoding and reconstructing missing points of a point cloud
CN111815776A (en) * 2020-02-04 2020-10-23 山东水利技师学院 Three-dimensional building fine geometric reconstruction method integrating airborne and vehicle-mounted three-dimensional laser point clouds and streetscape images

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090310867A1 (en) * 2008-06-12 2009-12-17 Bogdan Calin Mihai Matei Building segmentation for densely built urban regions using aerial lidar data
CN105139379A (en) * 2015-07-30 2015-12-09 滁州学院 Airborne Lidar point cloud building top surface gradual extraction method based on classifying and laying
CN106815841A (en) * 2017-01-13 2017-06-09 合肥师范学院 Image object partitioning algorithm based on T junction point clue
CN107944383A (en) * 2017-11-21 2018-04-20 航天科工智慧产业发展有限公司 Building roof patch division method based on three-dimensional Voronoi diagram
WO2020072865A1 (en) * 2018-10-05 2020-04-09 Interdigital Vc Holdings, Inc. A method and device for encoding and reconstructing missing points of a point cloud
CN111815776A (en) * 2020-02-04 2020-10-23 山东水利技师学院 Three-dimensional building fine geometric reconstruction method integrating airborne and vehicle-mounted three-dimensional laser point clouds and streetscape images

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
赵瑞斌等: "机载LiDAR点云中精细建筑物顶面的渐进提取", 《计算机辅助设计与图形学学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113822999A (en) * 2021-09-22 2021-12-21 杭州趣村游文旅集团有限公司 Building segmentation method of digital country three-dimensional model
CN113822999B (en) * 2021-09-22 2024-02-23 杭州趣村游文旅集团有限公司 Building segmentation method of digital rural three-dimensional model
CN116071530A (en) * 2023-03-07 2023-05-05 青岛市勘察测绘研究院 Building roof voxelized segmentation method based on airborne laser point cloud
CN116071530B (en) * 2023-03-07 2023-07-28 青岛市勘察测绘研究院 Building roof voxelized segmentation method based on airborne laser point cloud
CN117576144A (en) * 2024-01-15 2024-02-20 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment
CN117576144B (en) * 2024-01-15 2024-03-29 湖北工业大学 Laser point cloud power line extraction method and device and electronic equipment

Also Published As

Publication number Publication date
CN113205529B (en) 2022-04-29

Similar Documents

Publication Publication Date Title
CN113205529B (en) Method for segmenting top surface of building based on airborne LiDAR point cloud
Biosca et al. Unsupervised robust planar segmentation of terrestrial laser scanner point clouds based on fuzzy clustering methods
CN101877128B (en) Method for segmenting different objects in three-dimensional scene
CN110443810B (en) Point cloud plane segmentation method based on quick adjacent voxel query
CN110222642B (en) Plane building component point cloud contour extraction method based on global graph clustering
Zhao et al. Indoor point cloud segmentation using iterative gaussian mapping and improved model fitting
CN103106632B (en) A kind of fusion method of the different accuracy three dimensional point cloud based on average drifting
CN105139379A (en) Airborne Lidar point cloud building top surface gradual extraction method based on classifying and laying
Zhao et al. Automated recognition and measurement based on three-dimensional point clouds to connect precast concrete components
Xue et al. A derivative-free optimization-based approach for detecting architectural symmetries from 3D point clouds
Jung et al. A line-based progressive refinement of 3D rooftop models using airborne LiDAR data with single view imagery
CN112241676A (en) Method for automatically identifying terrain sundries
Chen et al. A local tangent plane distance-based approach to 3D point cloud segmentation via clustering
Hamid-Lakzaeian Structural-based point cloud segmentation of highly ornate building façades for computational modelling
Xu et al. A voxel-and graph-based strategy for segmenting man-made infrastructures using perceptual grouping laws: Comparison and evaluation
Wang et al. Roof plane segmentation from lidar point cloud data using region expansion based l 0 gradient minimization and graph cut
Xin et al. Rapid registration method by using partial 3D point clouds
CN110310322A (en) Method for detecting assembly surface of 10-micron-level high-precision device
Lu et al. A lightweight real-time 3D LiDAR SLAM for autonomous vehicles in large-scale urban environment
Demir Automated detection of 3D roof planes from Lidar data
Xin et al. Accurate and complete line segment extraction for large-scale point clouds
CN116452826A (en) Coal gangue contour estimation method based on machine vision under shielding condition
CN115147471A (en) Laser point cloud automatic registration method based on curvature density characteristics
CN114092775A (en) Method, system and equipment for extracting window of supervision type building
Huang et al. A multi-scale point clouds segmentation method for urban scene classification using region growing based on multi-resolution supervoxels with robust neighborhood

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