Building contour extraction and analysis method based on 3D laser point cloud data

Abstract

The present invention discloses a building contour extraction and analysis method based on three-dimensional laser point cloud data, and relates to the technical field of data analysis. The present invention performs smoothing processing on three-dimensional laser point cloud data, removes noise points, and performs sparse processing through a voxel grid filter to retain structural information; applies a ground point classification algorithm to separate ground points from building points in point cloud data; calculates the normal vector of each point, and analyzes local geometric features based on this, and uses a point cloud segmentation algorithm to divide the point cloud into different areas; identifies edge points in the building point concentration through a boundary detection algorithm, and generates preliminary building contours by combining multi-scale analysis and curve fitting; identifies and classifies free points, determines whether they are part of the building contour, and repairs them; smoothes the updated contour, compares it with reference data, verifies accuracy and optimizes it, and obtains the target building contour.

Description

Technical Field

[0001] The present invention relates to the technical field of data analysis, and in particular to a building contour extraction and analysis method based on three-dimensional laser point cloud data. Background Art

[0002] With the rapid development of 3D laser scanning technology, backpack 3D laser scanners have been widely used in the fields of building 3D modeling, urban planning, cultural relics protection, etc. due to their high efficiency, flexibility and high precision. However, after obtaining a large amount of 3D laser point cloud data, how to accurately and quickly extract the contours of buildings from these data has become an attractive problem to be solved. Traditional building contour extraction and analysis methods mostly rely on manual operations, such as extracting the edges of buildings through 2D image processing technology, or outlining using vector data in GIS (Geographic Information System). These methods are not only time-consuming and costly, but also easily affected by human factors, resulting in inaccurate and inconsistent extraction results. In addition, with the increase in the complexity of buildings, such as the presence of multiple curved surfaces, irregular shapes or occlusions, the limitations of traditional methods become more and more obvious.

[0003] In recent years, automated extraction methods based on 3D laser point cloud data have gradually attracted attention. These methods use the 3D spatial characteristics of point cloud data to automatically identify and extract the contours of buildings through algorithms. However, existing automated extraction methods still have some shortcomings. For example, in the process of extracting building contours based on 3D laser point cloud data, free points often appear, and existing technologies may not fully consider the local features around the free points, resulting in some free points with potential information being mistakenly deleted, which may cause the extracted contours to be inaccurate, affecting subsequent analysis and application. Summary of the invention

[0004] The purpose of the present invention is to provide a building contour extraction and analysis method based on three-dimensional laser point cloud data to solve the problems raised in the above background technology.

[0005] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0006] A building contour extraction and analysis method based on three-dimensional laser point cloud data, the method comprising the following steps:

[0007] Step S100. Acquire three-dimensional laser point cloud data, smooth the three-dimensional laser point cloud data, and remove noise points; use a voxel grid filter to sparse the smoothed three-dimensional laser point cloud data, so as to retain structural information; use a ground point classification algorithm to classify ground points and building points in the three-dimensional laser point cloud data into a ground point set and a building point set;

[0008] Step S200. For the building point set, calculate the normal vector of each point in the building point set, and analyze the local geometric features through the normal vector; apply the point cloud segmentation algorithm to segment the building point set into several different regions based on the local geometric features, and each region represents a type of point set with similar features;

[0009] Step S300. Using the normal vector and curvature boundary detection algorithm, identify the initial edge points in the building point set; using the multi-scale analysis method, dynamically adjust the size of the detection window to identify contour segments of different scales; combine the contour segments of different scales and the identified initial edge points to perform curve fitting to generate a preliminary building contour line;

[0010] Step S400. Based on the preliminary building outline, identify free points, and classify the free points using local geometric features and machine learning models to determine whether the free points are part of the building outline; repair the free points that belong to the building outline, and integrate them into the preliminary building outline using the geometric information of neighboring points, thereby obtaining the target building outline.

[0011] Furthermore, step S100 includes:

[0012] S101. Collect three-dimensional laser point cloud data, and store the three-dimensional laser point cloud data in XYZ coordinate format; for each point Pi (xi, yi, zi) in the three-dimensional laser point cloud data, calculate the corresponding neighborhood point set N (Pi), and the neighborhood point set is defined as: N (Pi) = {Pj|||Pi-Pj||≤r}, where xi, yi, zi represent the coordinate values ​​of point Pi on the X, Y and Z axes respectively; i and j both represent the number of the point in the three-dimensional laser point cloud data, and i≠j; ||Pi-Pj|| represents the Euclidean distance between point Pi and point Pj; r represents the radius of the neighborhood;

[0013] Based on the neighborhood point set N(Pi), the smoothed coordinate value is calculated for each point in the 3D laser point cloud data, and the corresponding calculation formula is:

[0014] hi'=[Σ Pj∈N(Pi) wi j×hj] / Σ Pj∈N(Pi) wij;

[0015] Where hi' represents the coordinate value of xi or yi or zi after smoothing of point Pi (xi, yi, zi), Pj represents the point in the neighborhood point set N (Pi); hj represents the coordinate value of point Pj (xj, yj or zj); wij represents the Gaussian weight;

[0016] S102. Use a voxel grid filter to thin out the smoothed 3D laser point cloud data. The specific analysis process is: divide the entire point cloud data space into a grid with a voxel size of v. For all points in each voxel, calculate the position coordinates of the centroid G of the points in the voxel. The specific calculation formula is:

[0017] G=((1 / |v|)Σ Pk∈v xk,(1 / |v|)Σ Pk∈v yk,(1 / |v|)Σ Pk∈v zk);

[0018] Wherein, |v| represents the number of points in voxel v; Pk represents the point in voxel v, and k represents the number of the point in the voxel; xk, yk and zk represent the coordinate values ​​of point Pk on the X, Y and Z axes respectively;

[0019] S103. Using a ground point classification algorithm, the ground points and building points in the three-dimensional laser point cloud data are classified into a ground point set and a building point set.

[0020] Furthermore, the specific process of classifying the ground points and building points in the three-dimensional laser point cloud data in S103 is as follows:

[0021] Obtain a smoothed and sparse 3D laser point cloud data set, preliminarily distinguish ground points and building points according to a preset height threshold Tg, compare the Z-axis coordinate value of each point in the smoothed and sparse 3D laser point cloud data set with the height threshold Tg, if the Z-axis coordinate value of the corresponding point is greater than the height threshold Tg, it is preliminarily classified as a building point set, if the Z-axis coordinate value of the corresponding point is less than or equal to the height threshold Tg, it is preliminarily classified as a ground point set;

[0022] Randomly select m points from the initially classified ground point set, fit a plane model based on the selected m points, and the corresponding fitting plane equation is: ax+by+cz+d=0, where a, b, c and d represent plane parameters, and x, y and z represent coordinate variables in three-dimensional space; for each point in the three-dimensional laser point cloud data set, calculate the distance to the fitting plane. If the distance from a point to the fitting plane is less than or equal to the preset threshold L, it is a ground point, otherwise it is a building point;

[0023] The preliminary classification result is compared with the classification result after plane fitting. If the two classifications are consistent, the original classification is maintained; if the two classifications are inconsistent, the following steps are performed: if the preliminary classification is a building point, but after plane fitting, it is found that the distance of the point is less than or equal to the threshold L, it is reclassified as a ground point; if the preliminary classification is a ground point, but after plane fitting, it is found that the distance of the point is greater than the threshold L, it is reclassified as a building point.

[0024] Further, step S200 includes:

[0025] S201. For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), and construct the covariance matrix C according to the neighborhood point set N(Pi), and C=(1 / n)Σ u∈[1,n] (Pu-p)(Pu-p) T , where n represents the number of points in the corresponding neighborhood point set N(Pi); Pu = (xu, yu, zu) represents the u-th point in the neighborhood point set N(Pi), and u is a positive integer from 1 to n; p represents the mean point of the point set N(Pi); according to the covariance matrix C, the eigenvalue and the corresponding eigenvector are obtained, and the eigenvector corresponding to the minimum eigenvalue of the covariance matrix is ​​used as the normal vector n, and the inverse of the minimum eigenvalue of the covariance matrix is ​​used as the curvature ρ;

[0026] S202. Calculate the angle θ between the normal vector n and the global normal vector n1, and the specific calculation formula is:

[0027] cosθ=(n·n1) / (|n|||n1|);

[0028] Among them, the global normal vector n1 is calculated by plane fitting based on the building point set; a point P0 is randomly selected from the building point set as a seed point, and the point P0 is added to a region. Starting from each point in the current region, its neighboring points are checked. If the normal vector and curvature of the neighboring point are similar to those of the current region, it is added to the current region, and the above process is repeated until the points in the region no longer meet the segmentation conditions;

[0029] S203. According to S201 and S202, the building point set is divided into several different regions, and each region represents a type of point set with similar characteristics.

[0030] Furthermore, the specific process of determining whether the normal vector and curvature of the neighboring point are similar to the normal vector and curvature of the current area in step S202 is as follows:

[0031] The measure value for judging whether the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region is defined as similarity S. The calculation formula of similarity S is: S=αS1+(1-α)S2, where α represents the weight coefficient; S1 represents the similarity between the normal vector of point P0 and the normal vector of the current region, and S1=n_P0·n_R, n_P0 represents the normal vector of point P0, and n_R represents the normal vector of the current region; S2 represents the similarity between the curvature of point P0 and the curvature of the current region, and S2=1 / (1+|ρ_P0-ρ_R|), ρ_P0 represents the curvature of point P0, and ρ_R represents the curvature of the current region; when the similarity S is greater than or equal to the similarity threshold S0, the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region.

[0032] Furthermore, step S300 includes:

[0033] S301. Using the boundary detection algorithm based on normal vector and curvature, identify the initial edge points in the building point set; set a multi-scale detection window, and the different detection window scales are W1, W2, ..., WM; for each detection window Wf (f is a positive integer from 1 to M), calculate the average edge strength B_Wf of all points in the window, and the specific calculation formula is:

[0034] B_Wf=(1 / |Wf|)Σ Pi∈Wf Bi;

[0035] Where |Wf| represents the number of points in the detection window Wf, Pi is the point in the detection window Wf, Bi represents the edge strength of the point in the detection window Wf; and each detection window corresponds to an edge strength average threshold Bf;

[0036] S302. Compare the edge intensity average value B_Wf with the corresponding edge intensity average threshold value Bf, mark the points in the window where the edge intensity average value B_Wf is greater than the corresponding edge intensity average threshold value Bf as edge points of the corresponding scale, and extract the corresponding contour line segment Ls according to the edge points of the corresponding scale at different detection window scales;

[0037] The contour segments Ls of all detection windows are obtained, and the contour segments of different scales and the identified initial edge points are combined for curve fitting to generate preliminary building contours.

[0038] Furthermore, the specific process of identifying the initial edge points in the building point set in S301 is as follows:

[0039] For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), normal vector ni and curvature ρi, calculate the edge strength Bi, and the specific calculation formula is:

[0040] Bi=β|Δni|+(1-β)|ρi|;

[0041] Among them, β represents the weight coefficient, |Δni| represents the change of the normal vector of point Pi, and |ρi| represents the absolute value of the curvature of point Pi; the edge intensity Bi is compared with the corresponding edge intensity threshold B0, and the point whose edge intensity Bi is greater than the edge intensity threshold B0 is marked as the initial edge point.

[0042] Furthermore, step S400 includes:

[0043] S401. Obtain the preliminary building outline and the points around the preliminary building outline, calculate the distance D between the points around the preliminary building outline and the preliminary building outline, and if D is greater than the set threshold D0, mark it as a free point, thereby forming a free point set;

[0044] S402. For each point in the free point set, calculate the corresponding normal vector and curvature; use the points on the preliminary building outline as a labeled data set to train the machine learning model, input the normal vector and curvature of each point in the free point set as the input local geometric features into the trained machine learning model, and output the category label of each free point belonging to the preliminary building outline, and the category label is 0 and 1, where 0 indicates that it does not belong to the preliminary building outline, and 1 indicates that it belongs to the preliminary building outline;

[0045] S403. For the free point Py with the category label 1, obtain the nearest contour segment L, calculate the projection point P_1 on the contour segment L, and the specific calculation formula is:

[0046] P_1=P_s+[(Py-P_s)·(P_e-P_s)] / |P_e-P_s| 2 ]·(P_e-P_s);

[0047] Among them, P_s represents the starting point of the contour segment L, P_e represents the end point of the contour segment L, |P_e-P_s| 2 Represents the square of the direction vector of the contour segment L, P_e-P_s represents the direction vector of the contour segment L; determine whether the projection point P_1 is on the contour segment L. If the projection point P_1 is on the contour segment L, the repair point Q is the projection point P_1. If the projection point P_1 is not on the contour segment L, the repair point Q is set to the point between P_s and P_e that is closest to the free point Py; for each free point Py with a category label of 1, calculate the distance from the point in the corresponding neighborhood point set to the repair point Q, and then perform weighted averaging to obtain the final position Q0 of the repair point, add the final position Q0 of the repair point to the preliminary building contour line, and generate the target building contour line.

[0048] Furthermore, the specific process of determining whether the projection point P_1 is on the contour segment L in S403 is as follows:

[0049] Calculate the vector P_11 from the projection point P_1 to the starting point P_s of the contour segment L, and express it as: P_11=P_1-P_s; calculate the vector P_12 from the projection point P_1 to the end point P_e of the contour segment L, and express it as: P_12=P_1-P_e; express the direction vector of the contour segment L as P_L; calculate the proportional coefficient t, and t=(P_11·P_L) / |P_L| 2 ; When the proportional coefficient 0≤t≤1, the projection point P_1 is on the contour segment L, otherwise the projection point P_1 is not on the contour segment L.

[0050] Compared with the prior art, the beneficial effects of the present invention are as follows: the present invention combines smoothing processing with voxel grid filters to ensure that the noise of point cloud data is removed while retaining important structural information; this method can effectively reduce the complexity of point cloud data and lay a good foundation for subsequent processing. By calculating normal vectors and applying point cloud segmentation algorithms to process local geometric features, this method can provide more accurate area division and identify the distribution of buildings and other objects. This processing method is very effective for segmenting and identifying different features. Using a boundary detection algorithm based on normal vectors and curvature combined with a multi-scale analysis method, the size of the detection window is dynamically adjusted to adapt to contours of different scales; this method can more comprehensively capture the contour features of the building and improve the accuracy of edge point recognition. Based on the preliminary building contour line, free points are identified and repaired, and classified and repaired in combination with local geometric features and machine learning models; this processing method can effectively retain potential information, thereby ensuring the accuracy of the extracted target building contour line. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0052] Figure 1 It is a flow chart of a building contour extraction and analysis method based on three-dimensional laser point cloud data of the present invention. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0054] See also Figure 1 , the present invention provides a technical solution:

[0055] A building contour extraction and analysis method based on three-dimensional laser point cloud data, the method comprising the following steps:

[0056] Step S100. Acquire three-dimensional laser point cloud data, smooth the three-dimensional laser point cloud data, and remove noise points; use a voxel grid filter to sparse the smoothed three-dimensional laser point cloud data, so as to retain structural information; use a ground point classification algorithm to classify ground points and building points in the three-dimensional laser point cloud data into a ground point set and a building point set;

[0057] Step S200. For the building point set, calculate the normal vector of each point in the building point set, and analyze the local geometric features through the normal vector; apply the point cloud segmentation algorithm to segment the building point set into several different regions based on the local geometric features, and each region represents a type of point set with similar features;

[0058] Step S300. Using the normal vector and curvature boundary detection algorithm, identify the initial edge points in the building point set; using the multi-scale analysis method, dynamically adjust the size of the detection window to identify contour segments of different scales; combine the contour segments of different scales and the identified initial edge points to perform curve fitting to generate a preliminary building contour line;

[0059] Step S400. Based on the preliminary building outline, identify free points, and classify the free points using local geometric features and machine learning models to determine whether the free points are part of the building outline; repair the free points that belong to the building outline, and integrate them into the preliminary building outline using the geometric information of neighboring points, thereby obtaining the target building outline.

[0060] Step S100 includes:

[0061] S101. Collect three-dimensional laser point cloud data, and store the three-dimensional laser point cloud data in XYZ coordinate format; for each point Pi (xi, yi, zi) in the three-dimensional laser point cloud data, calculate the corresponding neighborhood point set N (Pi), and the neighborhood point set is defined as: N (Pi) = {Pj|||Pi-Pj||≤r}, where xi, yi, zi represent the coordinate values ​​of point Pi on the X, Y and Z axes respectively; i and j both represent the number of the point in the three-dimensional laser point cloud data, and i≠j; ||Pi-Pj|| represents the Euclidean distance between point Pi and point Pj; r represents the radius of the neighborhood;

[0062] In this embodiment, for each point Pi (xi, yi, zi) in the three-dimensional laser point cloud data, all other points Pj are traversed, ||Pi-Pj|| is calculated, and if ||Pi-Pj||≤r, the point Pj is added to the neighborhood point set N(Pi).

[0063] Based on the neighborhood point set N(Pi), the smoothed coordinate value is calculated for each point in the 3D laser point cloud data, and the corresponding calculation formula is:

[0064] hi'=[Σ Pj∈N(Pi) wi j×hj] / Σ Pj∈N(Pi) wij;

[0065] Where hi' represents the coordinate value of xi or yi or zi after smoothing of point Pi (xi, yi, zi), Pj represents the point in the neighborhood point set N (Pi); hj represents the coordinate value of point Pj (xj, yj or zj); wij represents the Gaussian weight;

[0066] S102. Use a voxel grid filter to thin out the smoothed 3D laser point cloud data. The specific analysis process is: divide the entire point cloud data space into a grid with a voxel size of v. For all points in each voxel, calculate the position coordinates of the centroid G of the points in the voxel. The specific calculation formula is:

[0067] G=((1 / |v|)Σ Pk∈v xk,(1 / |v|)Σ Pk∈v yk,(1 / |v|)Σ Pk∈v zk);

[0068] Wherein, |v| represents the number of points in voxel v; Pk represents the point in voxel v, and k represents the number of the point in the voxel; xk, yk and zk represent the coordinate values ​​of point Pk on the X, Y and Z axes respectively;

[0069] S103. Using a ground point classification algorithm, the ground points and building points in the three-dimensional laser point cloud data are classified into a ground point set and a building point set.

[0070] The specific process of classifying the ground points and building points in the three-dimensional laser point cloud data in S103 is as follows:

[0071] Obtain a smoothed and sparse 3D laser point cloud data set, preliminarily distinguish ground points and building points according to a preset height threshold Tg, compare the Z-axis coordinate value of each point in the smoothed and sparse 3D laser point cloud data set with the height threshold Tg, if the Z-axis coordinate value of the corresponding point is greater than the height threshold Tg, it is preliminarily classified as a building point set, if the Z-axis coordinate value of the corresponding point is less than or equal to the height threshold Tg, it is preliminarily classified as a ground point set;

[0072] Randomly select m points from the initially classified ground point set, fit a plane model based on the selected m points, and the corresponding fitting plane equation is: ax+by+cz+d=0, where a, b, c and d represent plane parameters, and x, y and z represent coordinate variables in three-dimensional space; for each point in the three-dimensional laser point cloud data set, calculate the distance to the fitting plane. If the distance from a point to the fitting plane is less than or equal to the preset threshold L, it is a ground point, otherwise it is a building point;

[0073] The preliminary classification result is compared with the classification result after plane fitting. If the two classifications are consistent, the original classification is maintained; if the two classifications are inconsistent, the following steps are performed: if the preliminary classification is a building point, but after plane fitting, it is found that the distance of the point is less than or equal to the threshold L, it is reclassified as a ground point; if the preliminary classification is a ground point, but after plane fitting, it is found that the distance of the point is greater than the threshold L, it is reclassified as a building point.

[0074] Step S200 includes:

[0075] S201. For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), and construct the covariance matrix C according to the neighborhood point set N(Pi), and C=(1 / n)Σ u∈[1,n] (Pu-p)(Pu-p) T , where n represents the number of points in the corresponding neighborhood point set N(Pi); Pu = (xu, yu, zu) represents the u-th point in the neighborhood point set N(Pi), and u is a positive integer from 1 to n; p represents the mean point of the point set N(Pi); according to the covariance matrix C, the eigenvalue and the corresponding eigenvector are obtained, and the eigenvector corresponding to the minimum eigenvalue of the covariance matrix is ​​used as the normal vector n, and the inverse of the minimum eigenvalue of the covariance matrix is ​​used as the curvature ρ;

[0076] S202. Calculate the angle θ between the normal vector n and the global normal vector n1, and the specific calculation formula is:

[0077] cosθ=(n·n1) / (|n|||n1|);

[0078] Among them, the global normal vector n1 is calculated by plane fitting based on the building point set; a point P0 is randomly selected from the building point set as a seed point, and the point P0 is added to a region. Starting from each point in the current region, its neighboring points are checked. If the normal vector and curvature of the neighboring point are similar to those of the current region, it is added to the current region, and the above process is repeated until the points in the region no longer meet the segmentation conditions;

[0079] S203. According to S201 and S202, the building point set is divided into several different regions, and each region represents a type of point set with similar characteristics.

[0080] The specific process of determining whether the normal vector and curvature of the neighboring point are similar to the normal vector and curvature of the current area in step S202 is as follows:

[0081] The measure value for judging whether the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region is defined as similarity S. The calculation formula of similarity S is: S=αS1+(1-α)S2, where α represents the weight coefficient; S1 represents the similarity between the normal vector of point P0 and the normal vector of the current region, and S1=n_P0·n_R, n_P0 represents the normal vector of point P0, and n_R represents the normal vector of the current region; S2 represents the similarity between the curvature of point P0 and the curvature of the current region, and S2=1 / (1+|ρ_P0-ρ_R|), ρ_P0 represents the curvature of point P0, and ρ_R represents the curvature of the current region; when the similarity S is greater than or equal to the similarity threshold S0, the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region.

[0082] In this embodiment, it is assumed that the normal vector of the current area is n_R = (0.5, 0.5, 0.7), the curvature is ρ_R = 0.2, the normal vector of point P0 is n_P0 = (0.4, 0.6, 0.7), the curvature is ρ_P0 = 0.25, and α = 0.7; S1 = 0.99, S2 = 0.952; S = 0.979, assuming the similarity threshold S0 = 0.5, since S = 0.979 > S0 = 0.5, the point P0 has a high similarity with the current area and can be added to the current area.

[0083] Step S300 includes:

[0084] S301. Using the boundary detection algorithm based on normal vector and curvature, identify the initial edge points in the building point set; set a multi-scale detection window, and the different detection window scales are W1, W2, ..., WM; for each detection window Wf (f is a positive integer from 1 to M), calculate the average edge strength B_Wf of all points in the window, and the specific calculation formula is:

[0085] B_Wf=(1 / |Wf|)Σ Pi∈Wf Bi;

[0086] Where |Wf| represents the number of points in the detection window Wf, Pi is the point in the detection window Wf, Bi represents the edge strength of the point in the detection window Wf; and each detection window corresponds to an edge strength average threshold Bf;

[0087] S302. Compare the edge intensity average value B_Wf with the corresponding edge intensity average threshold value Bf, mark the points in the window where the edge intensity average value B_Wf is greater than the corresponding edge intensity average threshold value Bf as edge points of the corresponding scale, and extract the corresponding contour line segment Ls according to the edge points of the corresponding scale at different detection window scales;

[0088] The contour segments Ls of all detection windows are obtained, and the contour segments of different scales and the identified initial edge points are combined for curve fitting to generate preliminary building contours.

[0089] The specific process of identifying the initial edge points in the building point set in S301 is as follows:

[0090] For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), normal vector ni and curvature ρi, calculate the edge strength Bi, and the specific calculation formula is:

[0091] Bi=β|Δni|+(1-β)|ρi|;

[0092] Among them, β represents the weight coefficient, |Δni| represents the change of the normal vector of point Pi, and |ρi| represents the absolute value of the curvature of point Pi; the edge intensity Bi is compared with the corresponding edge intensity threshold B0, and the point whose edge intensity Bi is greater than the edge intensity threshold B0 is marked as the initial edge point.

[0093] Step S400 includes:

[0094] S401. Obtain the preliminary building outline and the points around the preliminary building outline, calculate the distance D between the points around the preliminary building outline and the preliminary building outline, and if D is greater than the set threshold D0, mark it as a free point, thereby forming a free point set;

[0095] S402. For each point in the free point set, calculate the corresponding normal vector and curvature; use the points on the preliminary building outline as a labeled data set to train the machine learning model, input the normal vector and curvature of each point in the free point set as the input local geometric features into the trained machine learning model, and output the category label of each free point belonging to the preliminary building outline, and the category label is 0 and 1, where 0 indicates that it does not belong to the preliminary building outline, and 1 indicates that it belongs to the preliminary building outline;

[0096] S403. For the free point Py with the category label 1, obtain the nearest contour segment L, calculate the projection point P_1 on the contour segment L, and the specific calculation formula is:

[0097] P_1=P_s+[(Py-P_s)·(P_e-P_s)] / |P_e-P_s| 2 ]·(P_e-P_s);

[0098] Among them, P_s represents the starting point of the contour segment L, P_e represents the end point of the contour segment L, |P_e-P_s| 2 Represents the square of the direction vector of the contour segment L, P_e-P_s represents the direction vector of the contour segment L; determine whether the projection point P_1 is on the contour segment L. If the projection point P_1 is on the contour segment L, the repair point Q is the projection point P_1. If the projection point P_1 is not on the contour segment L, the repair point Q is set to the point between P_s and P_e that is closest to the free point Py; for each free point Py with a category label of 1, calculate the distance from the point in the corresponding neighborhood point set to the repair point Q, and then perform weighted averaging to obtain the final position Q0 of the repair point, add the final position Q0 of the repair point to the preliminary building contour line, and generate the target building contour line.

[0099] The specific process of determining whether the projection point P_1 is on the contour segment L in S403 is as follows:

[0100] Calculate the vector P_11 from the projection point P_1 to the starting point P_s of the contour segment L, and express it as: P_11=P_1-P_s; calculate the vector P_12 from the projection point P_1 to the end point P_e of the contour segment L, and express it as: P_12=P_1-P_e; express the direction vector of the contour segment L as P_L; calculate the proportional coefficient t, and t=(P_11·P_L) / |P_L| 2 ; When the proportional coefficient 0≤t≤1, the projection point P_1 is on the contour segment L, otherwise the projection point P_1 is not on the contour segment L.

[0101] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device.

[0102] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A building contour extraction and analysis method based on three-dimensional laser point cloud data, characterized by: The method comprises the following steps: Step S100. Acquire three-dimensional laser point cloud data, smooth the three-dimensional laser point cloud data, and remove noise points; use a voxel grid filter to sparse the smoothed three-dimensional laser point cloud data, so as to retain structural information; use a ground point classification algorithm to classify ground points and building points in the three-dimensional laser point cloud data into a ground point set and a building point set; Step S200. For the building point set, calculate the normal vector of each point in the building point set, and analyze the local geometric features through the normal vector; apply the point cloud segmentation algorithm to segment the building point set into several different regions based on the local geometric features, and each region represents a type of point set with similar features; Step S300. Using the normal vector and curvature boundary detection algorithm, identify the initial edge points in the building point set; using the multi-scale analysis method, dynamically adjust the size of the detection window to identify contour segments of different scales; combine the contour segments of different scales and the identified initial edge points to perform curve fitting to generate a preliminary building contour line; Step S400. Based on the preliminary building outline, identify free points, and classify the free points using local geometric features and machine learning models to determine whether the free points are part of the building outline; repair the free points that belong to the building outline, and integrate them into the preliminary building outline using the geometric information of neighboring points, thereby obtaining the target building outline.

2. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 1 is characterized in that: The step S100 includes: S101. Collect three-dimensional laser point cloud data, and store the three-dimensional laser point cloud data in XYZ coordinate format; for each point Pi (xi, yi, zi) in the three-dimensional laser point cloud data, calculate the corresponding neighborhood point set N (Pi), and the neighborhood point set is defined as: N (Pi) = {Pj|||Pi-Pj||≤r}, where xi, yi, zi represent the coordinate values ​​of point Pi on the X, Y and Z axes respectively; i and j both represent the number of the point in the three-dimensional laser point cloud data, and i≠j; ||Pi-Pj|| represents the Euclidean distance between point Pi and point Pj; r represents the radius of the neighborhood; Based on the neighborhood point set N(Pi), the smoothed coordinate value is calculated for each point in the 3D laser point cloud data, and the corresponding calculation formula is: hi'=[Σ Pj∈N(Pi) wij×hj] / Σ Pj∈N(Pi) Wij; Where hi' represents the coordinate value of xi or yi or zi after smoothing of point Pi (xi, yi, zi), Pj represents the point in the neighborhood point set N (Pi); hj represents the coordinate value of point Pj (xj, yj or zj); wij represents the Gaussian weight; S102. Use a voxel grid filter to thin out the smoothed 3D laser point cloud data. The specific analysis process is: divide the entire point cloud data space into a grid with a voxel size of v. For all points in each voxel, calculate the position coordinates of the centroid G of the points in the voxel. The specific calculation formula is: G=((1 / |v|)Σ Pk∈v xk,(1 / |v|)Σ Pk∈v yk,(1 / |v|)Σ Pk∈v zk); Wherein, |v| represents the number of points in voxel v; Pk represents the point in voxel v, and k represents the number of the point in the voxel; xk, yk and zk represent the coordinate values ​​of point Pk on the X, Y and Z axes respectively; S103. Using a ground point classification algorithm, the ground points and building points in the three-dimensional laser point cloud data are classified into a ground point set and a building point set.

3. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 2 is characterized in that: The specific process of classifying the ground points and building points in the three-dimensional laser point cloud data in S103 is as follows: Obtain a smoothed and sparse 3D laser point cloud data set, preliminarily distinguish ground points and building points according to a preset height threshold Tg, compare the Z-axis coordinate value of each point in the smoothed and sparse 3D laser point cloud data set with the height threshold Tg, if the Z-axis coordinate value of the corresponding point is greater than the height threshold Tg, it is preliminarily classified as a building point set, if the Z-axis coordinate value of the corresponding point is less than or equal to the height threshold Tg, it is preliminarily classified as a ground point set; Randomly select m points from the initially classified ground point set, fit a plane model based on the selected m points, and the corresponding fitting plane equation is: ax+by+cz+d=0, where a, b, c and d represent plane parameters, and x, y and z represent coordinate variables in three-dimensional space; for each point in the three-dimensional laser point cloud data set, calculate the distance to the fitting plane. If the distance from a point to the fitting plane is less than or equal to the preset threshold L, it is a ground point, otherwise it is a building point; The preliminary classification result is compared with the classification result after plane fitting. If the two classifications are consistent, the original classification is maintained; if the two classifications are inconsistent, the following steps are performed: if the preliminary classification is a building point, but after plane fitting, it is found that the distance of the point is less than or equal to the threshold L, it is reclassified as a ground point; if the preliminary classification is a ground point, but after plane fitting, it is found that the distance of the point is greater than the threshold L, it is reclassified as a building point.

4. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 2 is characterized in that: The step S200 includes: S201. For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), and construct the covariance matrix C according to the neighborhood point set N(Pi), and C=(1 / n)Σ u∈[1,n] (Pu-p)(Pu-p) T , where n represents the number of points in the corresponding neighborhood point set N(Pi); Pu = (xu, yu, zu) represents the u-th point in the neighborhood point set N(Pi), and u is a positive integer from 1 to n; p represents the mean point of the point set N(Pi); according to the covariance matrix C, the eigenvalue and the corresponding eigenvector are obtained, and the eigenvector corresponding to the minimum eigenvalue of the covariance matrix is ​​used as the normal vector n, and the inverse of the minimum eigenvalue of the covariance matrix is ​​used as the curvature ρ; S202. Calculate the angle θ between the normal vector n and the global normal vector n1, and the specific calculation formula is: cosθ=(n·n1) / (|n|||n1|); Among them, the global normal vector n1 is calculated by plane fitting based on the building point set; a point P0 is randomly selected from the building point set as a seed point, and the point P0 is added to a region. Starting from each point in the current region, its neighboring points are checked. If the normal vector and curvature of the neighboring point are similar to those of the current region, it is added to the current region, and the above process is repeated until the points in the region no longer meet the segmentation conditions; S203. According to S201 and S202, the building point set is divided into several different regions, and each region represents a type of point set with similar characteristics.

5. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 4 is characterized in that: The specific process of determining in step S202 that the normal vector and curvature of the neighboring point are similar to the normal vector and curvature of the current area is as follows: The measure value for judging whether the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region is defined as similarity S. The calculation formula of similarity S is: S=αS1+(1-α)S2, where α represents the weight coefficient; S1 represents the similarity between the normal vector of point P0 and the normal vector of the current region, and S1=n_P0·n_R, n_P0 represents the normal vector of point P0, and n_R represents the normal vector of the current region; S2 represents the similarity between the curvature of point P0 and the curvature of the current region, and S2=1 / (1+|ρ_P0-ρ_R|), ρ_P0 represents the curvature of point P0, and ρ_R represents the curvature of the current region; when the similarity S is greater than or equal to the similarity threshold S0, the normal vector and curvature of the neighborhood point are similar to the normal vector and curvature of the current region.

6. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 4 is characterized in that: The step S300 includes: S301. Using the boundary detection algorithm based on normal vector and curvature, identify the initial edge points in the building point set; set a multi-scale detection window, and the different detection window scales are W1, W2, ..., WM; for each detection window Wf (f is a positive integer from 1 to M), calculate the average edge strength B_Wf of all points in the window, and the specific calculation formula is: B_Wf=(1 / |Wf|)Σ Pi∈Wf Yes; Where |Wf| represents the number of points in the detection window Wf, Pi is the point in the detection window Wf, Bi represents the edge strength of the point in the detection window Wf; and each detection window corresponds to an edge strength average threshold Bf; S302. Compare the edge intensity average value B_Wf with the corresponding edge intensity average threshold value Bf, mark the points in the window where the edge intensity average value B_Wf is greater than the corresponding edge intensity average threshold value Bf as edge points of the corresponding scale, and extract the corresponding contour line segment Ls according to the edge points of the corresponding scale at different detection window scales; The contour segments Ls of all detection windows are obtained, and the contour segments of different scales and the identified initial edge points are combined for curve fitting to generate preliminary building contours.

7. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 6 is characterized in that: The specific process of identifying the initial edge points in the building point set in S301 is as follows: For each point Pi(xi,yi,zi) in the building point set, obtain the corresponding neighborhood point set N(Pi), normal vector ni and curvature ρi, calculate the edge strength Bi, and the specific calculation formula is: Bi=β|Δni|+(1-β)|ρi|; Among them, β represents the weight coefficient, |Δni| represents the change of the normal vector of point Pi, and |ρi| represents the absolute value of the curvature of point Pi; the edge intensity Bi is compared with the corresponding edge intensity threshold B0, and the point whose edge intensity Bi is greater than the edge intensity threshold B0 is marked as the initial edge point.

8. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 6 is characterized in that: The step S400 includes: S401. Obtain the preliminary building outline and the points around the preliminary building outline, calculate the distance D between the points around the preliminary building outline and the preliminary building outline, and if D is greater than the set threshold D0, mark it as a free point, thereby forming a free point set; S402. For each point in the free point set, calculate the corresponding normal vector and curvature; use the points on the preliminary building outline as a labeled data set to train the machine learning model, input the normal vector and curvature of each point in the free point set as the input local geometric features into the trained machine learning model, and output the category label of each free point belonging to the preliminary building outline, and the category label is 0 and 1, where 0 indicates that it does not belong to the preliminary building outline, and 1 indicates that it belongs to the preliminary building outline; S403. For the free point Py with the category label 1, obtain the nearest contour segment L, calculate the projection point P_1 on the contour segment L, and the specific calculation formula is: P_1=P_s+[(Py-P_s)·(P_e-P_s)] / |P_e-P_s| 2 ]·(P_e-P_s); Among them, P_s represents the starting point of the contour segment L, P_e represents the end point of the contour segment L, |P_e-P_s| 2 Represents the square of the direction vector of the contour segment L, P_e-P_s represents the direction vector of the contour segment L; determine whether the projection point P_1 is on the contour segment L. If the projection point P_1 is on the contour segment L, the repair point Q is the projection point P_1. If the projection point P_1 is not on the contour segment L, the repair point Q is set to the point between P_s and P_e that is closest to the free point Py; for each free point Py with a category label of 1, calculate the distance from the point in the corresponding neighborhood point set to the repair point Q, and then perform weighted averaging to obtain the final position Q0 of the repair point, add the final position Q0 of the repair point to the preliminary building contour line, and generate the target building contour line.

9. The building contour extraction and analysis method based on three-dimensional laser point cloud data according to claim 8, characterized in that: The specific process of determining whether the projection point P_1 is on the contour segment L in S403 is as follows: Calculate the vector P_11 from the projection point P_1 to the starting point P_s of the contour segment L, and express it as: P_11=P_1-P_s; calculate the vector P_12 from the projection point P_1 to the end point P_e of the contour segment L, and express it as: P_12=P_1-P_e; express the direction vector of the contour segment L as P_L; calculate the proportional coefficient t, and t=(P_11·P_L) / |P_L| 2 ; When the proportional coefficient 0≤t≤1, the projection point P_1 is on the contour segment L, otherwise the projection point P_1 is not on the contour segment L.

Images