Context geometric global optimization-based point cloud building contour line extraction method

Abstract

The invention provides a point cloud building contour line extraction method based on context geometric global optimization, and relates to the technical field of airborne laser radar data processing, and the method comprises the steps: obtaining airborne laser radar point cloud data; processing the point cloud data by using a progressive triangulation network encryption filtering algorithm, and separating the point cloud data into ground points and non-ground points; obtaining an initial seed point, and extracting a building point; an Alpha Shapes algorithm is adopted to process the building points to obtain building contour points; combining the contour point groups to obtain building contour line segments; constructing an undirected weight map, segmenting the undirected weight map, and distributing a direction angle for each contour line segment; and generating a final building contour line according to the adjusted contour line segment. According to the invention, the collected airborne laser radar point cloud data of the building area is processed in a data driving mode to extract the building contour line, so that subsequent various applications based on the building contour line are facilitated.

Description

Technical Field

[0001] The present invention relates to the technical field of airborne laser radar data processing, and in particular to a point cloud building contour extraction method based on contextual geometry global optimization. Background Art

[0002] Automatically extracting accurate building contours from remote sensing data is widely used in fields such as building 3D reconstruction, urban map updates, and urban disaster monitoring. Airborne LiDAR (Light Detection And Ranging, LiDAR) can actively collect and obtain high-density and high-precision point cloud data. Therefore, LiDAR-based building contour extraction methods are widely used. However, due to the complex shape of buildings and the discrete and irregular characteristics of point clouds, the building contours extracted directly from point clouds have serious jagged problems and cannot accurately depict the true outline of the building. Therefore, how to accurately and automatically extract building contours is a difficult problem. At the same time, the shape of buildings has significant artificial design features, such as potential parallel and perpendicular geometric relationships between contour segments. How to integrate these features into the contour extraction process to ensure that the extracted contours meet human visual cognition is also a difficult problem.

[0003] Currently, many researchers have proposed various building outline extraction methods based on point cloud data. These methods generally include steps such as point cloud data preprocessing, building point extraction, and outline generation. Traditional methods mainly rely on local feature analysis of point cloud data, such as normal vector calculation, curvature analysis, and boundary detection, to identify and extract building outlines. With the development of technology, some methods based on multi-scale analysis and machine learning have also been introduced into building outline extraction to improve the accuracy and automation of extraction. However, these methods often focus on the local geometric features of the point cloud and ignore the global geometric relationships between building outline segments, resulting in extraction results that do not conform to the actual structural characteristics of the building.

[0004] Chinese invention patent 202411148205.4 discloses a building contour extraction and analysis method based on three-dimensional laser point cloud data. The method improves the accuracy of contour extraction through technical means such as smoothing, voxel grid filtering, normal vector calculation, and multi-scale analysis. However, this method has the following technical defects: First, the method mainly relies on local geometric features and boundary detection algorithms to identify building contours, and fails to fully consider the global geometric relationship between building contour segments, such as parallel and perpendicular constraints. Second, although the method introduces a free point repair mechanism, the repair process is only based on the geometric information of neighboring points. It does not transform the contour regularization problem into a global optimization problem and is prone to falling into local optimal solutions. In addition, the method does not have a specific extraction mechanism for the main direction of the building and cannot accurately grasp the overall structural characteristics of the building. As a result, the extracted contour line, although it fits the point cloud data, may not conform to human visual cognition habits, especially when dealing with complex building structures. Summary of the Invention

[0005] In view of this, the present invention proposes a point cloud building contour line extraction method based on contextual geometry global optimization. By considering the potential parallel, perpendicular and other geometric relationships between contour line segments, an undirected weighted graph containing data penalty terms and smoothness penalty terms is constructed and global optimization is performed based on the minimum cut principle. The extracted building contour lines not only meet the accuracy requirements, but also conform to human visual cognition, effectively solving the contour line extraction problem under complex building shapes, and providing more accurate building contour line data for applications such as building three-dimensional reconstruction, urban map updating, and urban disaster monitoring.

[0006] The technical solution of the present invention is achieved as follows:

[0007] The present invention provides a point cloud building contour extraction method based on contextual geometry global optimization, comprising:

[0008] S1. Obtain airborne lidar point cloud data;

[0009] S2. Processing the point cloud data using a progressive triangulation encryption filtering algorithm to separate the point cloud data into ground points and non-ground points;

[0010] S3, using a semi-suppressed fuzzy C-means clustering algorithm to obtain initial building seed points from the non-ground points, and performing region growing based on the initial building seed points to extract building points;

[0011] S4, processing the building points using the Alpha Shapes algorithm to obtain building contour points; performing straight line fitting on the building contour points and performing region growing clustering grouping, and determining the main direction of the building by calculating the angle between adjacent points in the contour point group; merging the contour point groups using a hierarchical clustering strategy to obtain building contour line segments;

[0012] S5. Construct an undirected weighted graph using the building contour line segments and the main directions of the building, and segment the undirected weighted graph based on the minimum cut principle, assigning a direction angle that conforms to the contextual geometric relationship to each contour line segment; perform straight line intersection calculation based on the adjusted contour line segments to obtain the building contour corner points, and connect the building contour corner points to generate the final building contour line.

[0013] Preferably, step S2 includes:

[0014] S21, dividing the entire point cloud data into regular grids, and defining the lowest point in each grid as a ground point;

[0015] S22, constructing a triangulated network using the acquired ground points and using it as a reference ground;

[0016] S23. For each triangle, count the unmarked points that fall within the triangulated area, and calculate the distance dist of each point to the nearest point in the triangulated area, as well as the angle α between the line connecting the point and the vertex of the triangle and the triangle face;

[0017] S24, if the distance dist of the unmarked point is greater than the distance threshold dist th And the angle α is greater than the angle threshold α th , then mark the unmarked point as a non-ground point, otherwise mark it as a ground point;

[0018] S25. Repeat steps S23-24 to update the triangulated network until all points are marked as ground points or non-ground points.

[0019] Preferably, in step S3, the process of the semi-suppressed fuzzy C-means clustering algorithm includes:

[0020] For each non-ground point, five features are calculated based on its cylindrical neighboring points, including the height difference Δh between the highest point and the lowest point, the average height h mean , plane roughness σ, planarity planness and point density density;

[0021] The feature set of all non-ground points is used as the input of the fuzzy C-means classifier for unsupervised clustering. The non-ground point set is divided into two categories according to the following principle:

[0022]

[0023] Where, d ij is the sampling point x j The distance to the i-th cluster center, u ij is the fuzzy membership matrix, n is the number of non-ground points, and c is the number of categories;

[0024] According to the geometric characteristics of building points, the category that meets the building characteristics in the two divided categories is used as the initial seed point of the building.

[0025] Preferably, in step S3, performing region growing based on the initial building seed point to extract building points includes:

[0026] A1. Select the initial seed points of the building as the seed set S and establish an empty result point set R;

[0027] A2. Calculate the normal vector n for each seed point M. M and height h M , and fit the plane through its neighboring points to obtain the plane equation ax+by+cz+d=0;

[0028] A3. Search the set T of k nearest neighbors of point M.

[0029] A4. For each point N in T, calculate its normal vector n N and height h N , and the distance d from the point N to the fitting plane of the seed point M;

[0030] A5. If point N satisfies the following three conditions:

[0031] The distance d from point N to the fitting plane M is less than the preset threshold d th ;

[0032] Normal vector n of point N N Normal vector n to the seed point M M The angle θ between MN Less than the preset threshold θ th ;

[0033] Height difference|h M -h N | is less than the preset threshold h th ;

[0034] Then add point N to the result point set R and to the seed set S;

[0035] A6. Repeat steps A2-A5 until all points in the seed set S are processed and finally R contains all identified building points.

[0036] A7. Cluster the identified building points according to the connectivity between points, remove the point sets with too small an area, and obtain the building points.

[0037] Preferably, in step S4, processing the building points using the Alpha Shapes algorithm includes:

[0038] Project the building points onto the xy plane to form a two-dimensional point set;

[0039] For any building point P(x P ,y P ), use KD-tree to search for the nearest neighbor point NP within 2α, select any point Q(x Q ,y Q )Calculate the center coordinates of the two circles O1(O x1 ,O y1 )、O2(O x2 ,O y2 );

[0040] For the remaining points in NP, if none of them fall within a circle calculated above, then point P is defined as a contour point;

[0041] Extract the contour points from the remaining building points to form a building contour point set.

[0042] Preferably, in step S4, performing straight line fitting on the building contour points and performing region growth clustering grouping includes:

[0043] B1. Use the least squares method to perform a straight line fitting on each building contour point to obtain the straight line fitting error of each contour point;

[0044] B2. Select the contour point with the smallest straight line fitting error as the initial seed point;

[0045] B3, establishing the initial contour point group L where the seed point is located;

[0046] B4. In the set of k nearest neighbor points of the seed point, select the point whose vertical distance to the fitting line of the midpoint of the current group L is less than the threshold D. th The points are added to group L;

[0047] B5. Update the straight line fitting parameters of group L;

[0048] B6. Repeat steps B4-B5 until no new points are added to group L;

[0049] B7. Reselect the point with the smallest fitting error among the ungrouped points as a new seed point;

[0050] B8. Repeat steps B3-B7 until all contour points are assigned to corresponding groups.

[0051] Preferably, in step S4, determining the main direction of the building by calculating the angle between adjacent points in the contour point group includes:

[0052] Group the contour points N = {p1(x1,y1,z1),…,p n (x n ,y n ,z n )}, calculate the angle α between any adjacent contour points i , and calculate the angle variance δ:

[0053] α i =acos(n1,n2)

[0054] n1=(x i-1 -x i ,y i-1 -y i )

[0055] n2=(x i+1 -x i ,y i+1 -y i )

[0056]

[0057] Where acos represents the arc cosine function, which is used to calculate the angle between two vectors; n1 is the vector pointing from the i-1th point to the ith point; n2 is the vector pointing from the i+1th point to the ith point; is the average value of all angles; n is the number of grouped contour points; δ is the angle variance;

[0058] A group of contour points with an angle variance less than a set threshold is selected, and the fitted straight line is defined as the main direction of the building.

[0059] Preferably, in step S4, merging the contour point groups using a hierarchical clustering strategy includes:

[0060] C1. Calculate the straight line fitting error after merging any two adjacent groups of contour points;

[0061] C2. If the minimum merge error is less than the set threshold, the two groups of contour points corresponding to the minimum merge error are merged;

[0062] C3, repeat steps C1 to C2 until the minimum merging error is greater than the set threshold, then stop merging;

[0063] C4. Perform straight line fitting on each final merged contour point group to obtain the building contour segment.

[0064] Preferably, step S5 includes:

[0065] S51, establishing an undirected weighted graph structure based on the geometric relationships of the acquired building contour segments, where each node of the graph represents a contour segment, and the edges between the nodes represent the connection relationship between the segments;

[0066] S52. For each contour segment, calculate its data penalty value under the two constraints of the main direction and the original direction, and construct the node weight of the undirected weight graph;

[0067] S53. For any two adjacent contour line segments, calculate a smoothing penalty value of their direction angle difference, and construct an edge weight of an undirected weighted graph;

[0068] S54. Based on the constructed undirected weighted graph, the minimum cut principle is used to perform graph segmentation optimization and assign the optimal direction angle to each contour segment;

[0069] S55. Refit each contour line segment according to the optimized direction angle to obtain an adjusted contour line segment;

[0070] S56. Calculate the intersection points of adjacent contour line segments as building contour corner points, and connect all corner points to generate a closed building contour line.

[0071] Preferably, the calculation formula of the data penalty value in step S52 is:

[0072]

[0073] Where D p (l p ) is the data penalty value, l p is a straight line fitted by least squares using the main directions and grouped contour points, p i is the i-th contour point belonging to the contour segment p, dist(p i ,l p ) is point p i To the straight line l p The vertical distance, a lp 、b lp 、c lp is a straight line l p The parameter x i ,y i For p i The plane coordinates of , n is the number of contour points contained in the contour segment p;

[0074] The calculation formula of the smoothing penalty value in step S53 is:

[0075] V p,q (l p ,l q)=exp(-|α p -α q |)

[0076] Where, l p and l p For two adjacent contour segments, α p and α q are line segments l p and l q direction angle.

[0077] The present invention has the following beneficial effects compared to the prior art:

[0078] (1) The point cloud building contour line extraction method based on contextual geometry global optimization proposed in the present invention realizes the accurate extraction of building contour lines by transforming contour line regularization into a global optimization problem, constructing an undirected weighted graph containing data penalty terms and smoothness penalty terms, and performing segmentation based on the minimum cut principle. This effectively solves the technical difficulties of poor accuracy and serious jaggedness problems of building contour lines directly extracted from point clouds in traditional methods, so that the extracted contour lines meet the accuracy requirements while conforming to human visual cognition, providing more accurate building contour line data for applications such as building three-dimensional reconstruction, urban map updating, and urban disaster monitoring.

[0079] (2) The present invention proposes a building main direction detection method based on grouped contour point angle analysis. The main direction of the building is determined by calculating the angle between adjacent points in the contour point group and their variance. Compared with the method of determining the main direction based on length constraints, this method can accurately extract the main direction of the building and is not limited by the number of main directions.

[0080] (3) The present invention adopts the semi-suppressed fuzzy C-means clustering algorithm to obtain the initial seed points of buildings from non-ground points, and extracts building points by region growth based on the seed points. The two parameters of normal vector angle and height difference are constrained to effectively improve the accuracy of building point recognition;

[0081] (4) The present invention adopts a hierarchical clustering strategy to merge grouped contour points. By calculating the straight line fitting error after merging any two adjacent groups of contour points and merging them based on the error, the contour line segments that meet the structural characteristics of the building are adaptively obtained, avoiding the jagged problem caused by excessive segmentation of the contour line and improving the smoothness and continuity of the contour line.

[0082] (5) The present invention takes into account the potential parallel, perpendicular and other geometric relationships between contour segments during the contour regularization process. By constructing an undirected weight graph and combining the data penalty term with the smoothness penalty term, the global optimization of the contour segments is achieved. The final extracted building contour line not only maintains consistency with the original point cloud data, but also meets the requirements of the building contour line regularity, thereby improving the geometric accuracy of the contour line. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0084] Figure 1 is a flow chart of the method of the present invention;

[0085] Figure 2 It is a flowchart of the technology implementation of the present invention;

[0086] Figure 3 Selecting and constructing a triangulated network diagram for the initial seed point of the present invention;

[0087] Figure 4 Schematic diagram of angle and distance calculation according to the present invention;

[0088] Figure 5 Schematic diagram of contour point extraction using the Alpha Shapes algorithm of the present invention;

[0089] Figure 6 This is a schematic diagram of contour point grouping according to the present invention;

[0090] Figure 7 This is a schematic diagram of the main direction of the contour line of the present invention;

[0091] Figure 8 This is a schematic diagram of contour segment extraction according to the present invention;

[0092] Figure 9 A schematic diagram of an undirected weighted graph constructed by contour segments and main directions of the present invention;

[0093] Figure 10 This is a schematic diagram of the building contour extraction results of the present invention. DETAILED DESCRIPTION

[0094] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the embodiments described 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 making creative efforts are within the scope of protection of the present invention.

[0095] like Figure 1 As shown, the present invention provides a point cloud building contour line extraction method based on contextual geometry global optimization, comprising:

[0096] S1. Obtain airborne lidar point cloud data;

[0097] S2. Processing the point cloud data using a progressive triangulation encryption filtering algorithm to separate the point cloud data into ground points and non-ground points;

[0098] S3, using a semi-suppressed fuzzy C-means clustering algorithm to obtain initial building seed points from the non-ground points, and performing region growing based on the initial building seed points to extract building points;

[0099] S4, processing the building points using the Alpha Shapes algorithm to obtain building contour points; performing straight line fitting on the building contour points and performing region growing clustering grouping, and determining the main direction of the building by calculating the angle between adjacent points in the contour point group; merging the contour point groups using a hierarchical clustering strategy to obtain building contour line segments;

[0100] S5. Construct an undirected weighted graph using the building contour line segments and the main directions of the building, and segment the undirected weighted graph based on the minimum cut principle, assigning a direction angle that conforms to the contextual geometric relationship to each contour line segment; perform straight line intersection calculation based on the adjusted contour line segments to obtain the building contour corner points, and connect the building contour corner points to generate the final building contour line.

[0101] like Figure 2As shown in the figure, the technical implementation process of the present invention is as follows: first, the airborne lidar point cloud data is introduced as the processing object; second, the progressive triangulation encryption filtering algorithm is used to filter the point cloud data to separate the point cloud data into ground points and non-ground points, providing a basis for subsequent building identification; then, the semi-suppressed fuzzy C-means clustering method is used to obtain the initial building seed points from the non-ground points, and the region growing clustering method is used based on these seed points to identify the complete building points; then, the Alpha The Shapes algorithm processes building points to obtain building contour points, fits these contour points to straight lines and groups them through region growing clustering, determines the main direction of the building by calculating the angle between adjacent points in the contour point group, and merges the contour point groups using a hierarchical clustering strategy to obtain the initial building contour segments; finally, an undirected weighted graph is constructed using the building contour segments and the main directions of the building, and the undirected weighted graph is segmented based on the minimum cut principle, assigning a direction angle that conforms to the contextual geometric relationship to each contour segment, and performing a straight line intersection calculation based on the adjusted contour segments to obtain the building contour corner points, which are then connected to generate the final building contour line, so that the extracted contour line meets the accuracy requirements while conforming to the geometric characteristics of the building and human visual cognition.

[0102] Specifically, in one embodiment of the present invention, step S1 includes:

[0103] By equipping an aerial vehicle with a complete LiDAR system, which includes components such as an inertial measurement unit (IMU), differential GPS (DGPS), a laser scanning ranging system, and an imaging device, the system actively scans building scenes and collects point cloud data. During flight, the airborne LiDAR system emits laser pulses to the ground. By measuring the time difference between the emission of the laser pulse and the reception of the reflected echo, it calculates the distance the laser pulse has traveled. Combined with the carrier position and attitude information provided by the IMU and DGPS, the three-dimensional spatial coordinates of the ground target point are obtained. The airborne LiDAR point cloud data is stored in the form of a set of three-dimensional coordinate points, each point containing XYZ coordinate information as well as attribute information such as intensity and number of echoes. These point cloud data are characterized by high density and high precision, and can accurately reflect the geometric characteristics of ground buildings.

[0104] Specifically, in one embodiment of the present invention, step S2 includes:

[0105] S21, dividing the entire point cloud data into regular grids, and defining the lowest point in each grid as a ground point;

[0106] S22, constructing a triangulated network using the acquired ground points and using it as a reference ground;

[0107] S23. For each triangle, count the unmarked points that fall within the triangulated area, and calculate the distance dist of each point to the nearest point in the triangulated area, as well as the angle α between the line connecting the point and the vertex of the triangle and the triangle face;

[0108] S24, if the distance dist of the unmarked point is greater than the distance threshold dist th And the angle α is greater than the angle threshold α th , then mark the unmarked point as a non-ground point, otherwise mark it as a ground point;

[0109] S25. Repeat steps S23-24 to update the triangulated network until all points are marked as ground points or non-ground points.

[0110] In this embodiment, the progressive triangulation encryption filtering algorithm is used to filter the airborne lidar point cloud data obtained in step S1, and the point cloud data is separated into ground points and non-ground points. First, the point cloud data is preprocessed by rough filtering to remove obvious noise points and outliers. Then, the lowest point is selected from the point cloud data as the seed point, and the initial triangulation is constructed by performing Delaunay triangulation on all seed points to generate the initial ground model. The schematic diagram of initial seed point selection and triangulation construction is shown in FIG. Figure 3 shown.

[0111] The point cloud data is then classified using a progressive triangulation filtering algorithm. The core of this algorithm is to iteratively add unclassified points that meet specific criteria to the ground model, thereby separating ground and non-ground points in the point cloud data. Each iteration consists of the following three main steps:

[0112] (1) Calculate the distance and normal vector angle from the unclassified point to the current triangulation. For each unclassified point, calculate its vertical distance dist and normal vector angle α to the current triangulation. The vertical distance dist represents the distance from the point to the triangle plane where its projection is located, and the normal vector angle α represents the angle between the point and the normal vector of the triangle plane where its projection is located. The schematic diagram is as follows: Figure 4 Specifically, the calculation formulas for the distance dist and the angle α are as follows:

[0113]

[0114] Where, is the vector formed by the unmarked point and the triangulation vertex, is the normal vector corresponding to the triangulated surface, the triangle plane equation is ax+by+cz+d=0, and arccos() is the inverse cosine function.

[0115] (2) According to the set distance threshold dist th and angle threshold αth , determine whether the unclassified point is a ground point. The judgment conditions are:

[0116] If dist <dist th And α<α th , then the point is classified as a ground point;

[0117] If dist≥dist th or α≥α th , then the point is classified as a non-ground point.

[0118] (3) Re-perform Delaunay triangulation on the newly added ground points and the existing ground points to update the triangulated network model.

[0119] The above three steps are repeated until all points are classified as ground points or non-ground points, or the preset maximum number of iterations is reached. In practical applications, the distance threshold dist is set according to different terrain characteristics and building distribution. th and angle threshold α th It can be adjusted appropriately to obtain more accurate classification results.

[0120] The advantage of the progressive triangulation encryption filtering algorithm is that it can adaptively process areas with complex terrain and has a good filtering effect on point cloud data in areas with dense urban buildings.

[0121] Specifically, in one embodiment of the present invention, step S3 is to identify building points from non-ground points, specifically including using a semi-suppressed fuzzy C-means clustering method to obtain initial building seed points from non-ground points, and then using a region growing clustering method based on the obtained initial seed points to identify building points from non-ground points.

[0122] Specifically include:

[0123] S3.1 Semi-suppressed fuzzy C-means clustering algorithm to extract initial building seed points

[0124] First, the semi-suppressed fuzzy C-means clustering algorithm is applied to the non-ground points obtained in step S2 to extract the initial seed points of the building. The process includes:

[0125] For each non-ground point, five features are calculated based on its cylindrical neighboring points, including the height difference Δh between the highest point and the lowest point, the average height h mean , plane roughness σ, planarity planness and point density density; the calculation formula of the five features is as follows:

[0126] Δh=h max -h min

[0127]

[0128] Where h max With h min h is the height of the highest and lowest points among the neighboring points of the cylinder; i is the elevation of a point among the neighboring points, d i is the distance from the point to the surface, n is the number of neighboring points of the cylinder, s is the bottom area of the cylinder, λ1, λ2, and λ3 are the three eigenvalues corresponding to the decomposition of the covariance matrix of the neighboring points, where λ1≤λ2≤λ3.

[0129] The feature set of all non-ground points is used as the input of the fuzzy C-means classifier for unsupervised clustering. The non-ground point set is divided into two categories according to the following principle:

[0130]

[0131] Where, d ij is the sampling point x j The distance to the i-th cluster center, U is the fuzzy membership matrix, and its element u ij is the membership of the jth sample point to the ith cluster, V represents the cluster center set, that is, V=(v1,v2,....,v c ), where v i is the coordinate of the center point of the i-th cluster. In the feature space, each cluster center is composed of five feature dimensions (elevation difference, average elevation, surface roughness, planarity, and point density); n is the number of non-ground points, m is the fuzziness parameter that controls the fuzziness of the cluster, and c is the number of categories. In this embodiment, c = 2. The fuzzy C-means algorithm updates the membership matrix U and cluster center V through an iterative optimization process until the objective function J is obtained. m (U,V) converges to the minimum value.

[0132] Based on the geometric characteristics of building points, the points that meet the building characteristics in the two divided categories are selected as initial building seed points. Building points have the following geometric characteristics: large elevation difference, high average elevation, low surface roughness, good planarity, and high point density. Points that meet these geometric characteristics are selected as initial building seed points.

[0133] S3.2 Extracting complete building points based on the region growing algorithm of the initial seed point

[0134] After obtaining the initial building seed points, the region growing algorithm is used to extract the complete building points. , including:

[0135] A1. Select the initial seed points of the building as the seed set S and establish an empty result point set R;

[0136] A2. Calculate the normal vector n for each seed point M.M and height h M , and fit the plane through its neighboring points to obtain the plane equation ax+by+cz+d=0;

[0137] A3. Search for the set of k nearest neighbors T of point M. Set the value of k between 8 and 15 to ensure that there are enough points for judgment.

[0138] A4. For each point N in T, calculate its normal vector n N and height h N , and the distance d from the point N to the fitting plane of the seed point M;

[0139] A5. If point N satisfies the following three conditions:

[0140] The distance d from point N to the fitting plane M is less than the preset threshold d th ;

[0141] Normal vector n of point N N Normal vector n to the seed point M M The angle θ between MN Less than the preset threshold θ th ;

[0142] Height difference|h M -h N | is less than the preset threshold h th ;

[0143] Then add point N to the result point set R and to the seed set S;

[0144] A6. Repeat steps A2-A5 until all points in the seed set S are processed and finally R contains all identified building points.

[0145] A7. Cluster the identified building points based on inter-point connectivity, remove the point sets that occupy too small an area, obtain building points, and complete building detection.

[0146] In this embodiment, the distance d and the angle θ MN The calculation formula is as follows:

[0147]

[0148] in, is the normal vector of the point to be judged N and the seed point M.

[0149] This region-growing algorithm, based on initial seed points, combines three criteria: distance from the point to the plane, normal vector angle, and height difference. It effectively identifies building points and distinguishes them from other non-ground points (such as vegetation and vehicles). The core idea of the region-growing algorithm is to start with a reliable seed point and gradually incorporate neighboring points with similar geometric characteristics into the building point set, thereby achieving complete building recognition.

[0150] The plane distance constraint ensures that the selected points lie on the same plane or a near-plane, the normal angle constraint ensures the points' consistent orientation, and the height difference constraint prevents misclassification due to sudden changes in height. The combined application of these three constraints enables the algorithm to accurately distinguish building roofs from other non-ground objects.

[0151] The final connectivity clustering and area filtering further improve the accuracy of building recognition, eliminating point sets that may be formed by noise or small non-building targets, ensuring that the final extracted building point cloud has high completeness and accuracy.

[0152] Specifically, in one embodiment of the present invention, step S4 is the process of extracting building contour points and analyzing building contour lines after obtaining the building point cloud. This step uses the Alpha Shapes algorithm to process the building points in step S3 to obtain contour points, and then uses region growing and hierarchical clustering methods to process these contour points to obtain the main direction of the building and the initial building contour line segment. It specifically includes:

[0153] S4.1 uses the Alpha Shapes algorithm to extract building contour points, including:

[0154] Project building points onto the xy plane to form a two-dimensional point set; simplify the complexity of contour extraction.

[0155] For any building point P(x P ,y P ), use KD-tree to search for the nearest neighbor point NP within 2α, where α is the preset radius parameter to control the precision of contour point extraction;

[0156] Select any point Q(x Q ,y Q )Calculate the center coordinates of the two circles O1(O x1 ,O y1 )、O2(O x2 ,O y2 ), both circles have radius α, and both circles pass through points P and Q:

[0157] O x1 =x p +0.5(x Q -xP )-H(y Q -y P )

[0158] O y1 =y p +0.5(y Q -y P )-H(x P -x Q )

[0159] O x2 =x p +0.5(x Q -x P )+H(y Q -y P )

[0160] O y2 =y p +0.5(y Q -y P )+H(x P -x Q )

[0161]

[0162] Where x p ,y P is the coordinate of building point P, x Q ,y Q is the coordinate of any point Q in the point set NP, and S represents the distance between two points;

[0163] Extract the outer boundary edges of the remaining triangles to form a set of building contour points. Specifically, for the remaining points in NP, if none of them fall within a circle calculated above, then define point P as a contour point. The Alpha Shapes algorithm extracts contour points as shown in the following figure. Figure 5 shown.

[0164] S4.2 performs straight line fitting on the contour points to calculate the fitting error and uses the region growing method to group the contour points, including:

[0165] B1. Use the least squares method to perform a straight line fitting on each building contour point to obtain the straight line fitting error of each contour point;

[0166] B2. Select the contour point with the smallest straight line fitting error as the initial seed point;

[0167] B3, establishing the initial contour point group L where the seed point is located;

[0168] B4. In the set of k nearest neighbor points of the seed point, select the point whose vertical distance to the fitting line of the midpoint of the current group L is less than the threshold D.th The points are added to group L;

[0169] B5. Update the straight line fitting parameters of group L;

[0170] B6. Repeat steps B4-B5 until no new points are added to group L;

[0171] B7. Reselect the point with the smallest fitting error among the ungrouped points as a new seed point;

[0172] B8. Repeat steps B3-B7 until all contour points are assigned to corresponding groups.

[0173] In this embodiment, for any contour point and its neighboring points, a linear fit is performed using the least squares method, and the linear fitting error σ is obtained:

[0174]

[0175] Where a and b are the parameters of the fitted line y = ax + b, and n is the number of points of the fitted line. These parameters describe the local linear characteristics of the contour points and provide a basis for subsequent contour point grouping.

[0176] Based on the straight line fitting error, the contour points are grouped using the region growing method: (1) The point with the smallest fitting error is selected as the initial seed point, and the linear characteristics of the region where the point is located are the most obvious. (2) The region growth is realized based on the vertical distance from the point to the straight line, that is, the points whose vertical distance to the straight line where the current seed point is located is less than the threshold are selected to be added to the group. (3) The above process is repeated until all contour points are assigned to different groups, and finally the contour points are grouped. The specific diagram is as follows Figure 6 shown.

[0177] S4.3 Calculate the angle characteristics of the contour point group to determine the main direction of the building, including:

[0178] Group the contour points N = {p1(x1,y1,z1),…,p n (x n ,y n ,z n )}, calculate the angle α between any adjacent contour points i , and calculate the angle variance δ:

[0179] α i =acos(n1,n2)

[0180] n1=(x i-1 -x i ,y i-1 -y i )

[0181] n2=(xi+1 -x i ,y i+1 -y i )

[0182]

[0183] Where acos represents the arc cosine function, which is used to calculate the angle between two vectors; n1 is the vector pointing from the i-1th point to the ith point; n2 is the vector pointing from the i+1th point to the ith point; is the average value of all angles; n is the number of grouped contour points; δ is the angle variance;

[0184] Select the contour point group whose angle variance is less than the set threshold, and define the fitting straight line as the main direction of the building. Figure 7 As shown in Figure 3, this method can accurately extract the main directions of buildings and is not limited by the number of main directions.

[0185] S4.4 merges contour point groups to form contour line segments through hierarchical clustering, including:

[0186] C1. Calculate the straight line fitting error after merging any two adjacent groups of contour points;

[0187] C2. If the minimum merge error is less than the set threshold, the two groups of contour points corresponding to the minimum merge error are merged;

[0188] C3. Repeat steps C1 to C2 until the minimum merging error is greater than the set threshold, then stop merging; this bottom-up hierarchical clustering method can effectively reduce the number of contour segments while maintaining the accuracy of the contour.

[0189] C4. Perform straight line fitting on each final merged contour point group to obtain the building contour line segment. Each group of merged contour points corresponds to a building contour line segment, such as Figure 8 As shown, each line segment is fitted using the least squares method to obtain the parameter representation of the line segment.

[0190] Through the above steps, this method successfully extracted the building's contour points and analyzed the building's main orientation and initial contour segments. These results laid the foundation for subsequent contour optimization, significantly reducing jagged edges and making the extracted contour smoother and more accurate, better reflecting the building's true shape.

[0191] Specifically, in one embodiment of the present invention, step S5 is a process of constructing a target optimization function using the main building directions and the initial building contour segments, implementing contour regularization, and generating building contours. This step constructs an undirected weight graph, combines data penalty terms and smoothing penalty terms, and uses a minimum cut technique to achieve global optimization of the contours. It specifically includes:

[0192] S51, establishing an undirected weighted graph structure based on the geometric relationships of the acquired building contour segments, where each node of the graph represents a contour segment, and the edges between the nodes represent the connection relationship between the segments;

[0193] S52. For each contour segment, calculate its data penalty value under the two constraints of the main direction and the original direction, and construct the node weight of the undirected weight graph;

[0194] S53. For any two adjacent contour line segments, calculate a smoothing penalty value of their direction angle difference, and construct an edge weight of an undirected weighted graph;

[0195] S54. Based on the constructed undirected weighted graph, the minimum cut principle is used to perform graph segmentation optimization and assign the optimal direction angle to each contour segment;

[0196] S55. Refit each contour line segment according to the optimized direction angle to obtain an adjusted contour line segment;

[0197] S56. Calculate the intersection points of adjacent contour line segments as building contour corner points, and connect all corner points to generate a closed building contour line.

[0198] In a specific example, the implementation process of step S5 is as follows:

[0199] S5.1 Calculate the data penalty term between the initial contour line and the main direction based on its direction angle

[0200] First, calculate the data penalty term D between each initial contour segment and the main direction p (l p ), which measures the fitting error after the contour segment is adjusted to the main direction:

[0201]

[0202]

[0203] Where D p (l p ) is the data penalty value, l p is a straight line fitted by least squares using the main directions and grouped contour points, p i is the i-th contour point belonging to the contour segment p, dist(p i,l p ) is point p i To the straight line l p The vertical distance, a lp 、b lp 、c lp is a straight line l p The parameter x i ,y i For p i The plane coordinates of the contour segment p, n is the number of contour points contained in the contour segment p; D p (l p ) is smaller, the greater the fitting error after the contour segment is adjusted to the main direction, which means that the cost of the adjustment is higher.

[0204] S5.2 Calculate the smoothing penalty term between adjacent line segments based on the direction angles of adjacent line segments

[0205] Then, the smoothness penalty term V between adjacent contour segments is calculated p,q (l p ,l q ), which encourages smooth transitions between adjacent segments:

[0206] V p,q (l p ,l q )=exp(-|α p -α q |)

[0207] Where, l p and l q For two adjacent contour segments, α p and α q are line segments l p and l q The smaller the difference in direction angles, the larger the smoothing penalty value, indicating that the directions of adjacent line segments are closer and the smoother the image. This design helps maintain smooth transitions between building outline segments and avoids jagged or unnatural turns in the outline.

[0208] S5.3 Construct an undirected weighted graph and segment it using the minimum cut technique

[0209] like Figure 9 As shown, based on the data penalty term and smoothing penalty term calculated above, an undirected weighted graph G{V,E} is constructed. The node V in the graph represents the building contour segment, and each node can be assigned to any main direction; the edge E in the graph connects adjacent contour segments, and the weight of the edge is determined by the smoothing penalty term V. p,q (l p ,l q ) is determined; each node is connected to a special source point s and sink point t, and the connection weight is determined by the data penalty term Dp (l p )Decide.

[0210] The undirected weighted graph is segmented using the minimum cut technique to determine the optimal orientation angle for each contour segment. The minimum cut algorithm seeks a cut in the graph that minimizes the total weight of the cuts, thereby achieving a globally optimal orientation distribution for building contour segments. This global optimization method based on graph segmentation simultaneously considers data fitting errors and inter-segment smoothness, achieving more regularized building contour extraction.

[0211] S5.4 Adjust the building outline according to the new direction angle

[0212] Each contour segment is fitted using the least squares method based on the new orientation angle obtained after minimum cut segmentation and the corresponding contour points to obtain the adjusted contour segment. This step adjusts the initial contour segment to be consistent with or perpendicular to the main direction of the building, making the extracted building contour more regular and consistent with the building's geometric characteristics.

[0213] S5.5 Calculate building outline corner points and generate final outline

[0214] Finally, using the adjusted contour segments, the building contour corners are calculated by straight line intersection, and these corners are connected to generate the final building contour line (such as Figure 10 The formula for calculating the coordinates of the intersection of adjacent straight lines is as follows:

[0215]

[0216] Where a1, b1 and a2, b2 are the line parameters of adjacent lines. In this way, the corner points of the building's outline can be accurately located, and the generated outline more accurately expresses the building's geometric shape.

[0217] The proposed method for extracting building outlines from point clouds, based on global optimization of contextual geometry, considers potential parallel and perpendicular relationships between outline segments during the contour regularization process. This ensures that the extracted outlines meet extraction accuracy while also being consistent with human visual perception. This global optimization strategy effectively addresses the aliasing problem that often occurs when extracting building outlines directly from point clouds, improving both the accuracy and visual quality of building outline extraction.

[0218] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A point cloud building contour extraction method based on contextual geometry global optimization, characterized in that: include: S1. Obtain airborne lidar point cloud data; S2. Processing the point cloud data using a progressive triangulation encryption filtering algorithm to separate the point cloud data into ground points and non-ground points; S3, using a semi-suppressed fuzzy C-means clustering algorithm to obtain initial building seed points from the non-ground points, and performing region growing based on the initial building seed points to extract building points; S4, using the Alpha Shapes algorithm to process the building points to obtain building outline points; Performing straight line fitting on the building contour points and performing region growth clustering grouping, and determining the main direction of the building by calculating the angle between adjacent points in the contour point group; A hierarchical clustering strategy is used to merge the contour point groups to obtain building contour segments; S5. Constructing an undirected weighted graph using the building contour segments and the main directions of the building, and segmenting the undirected weighted graph based on a minimum cut principle, assigning a direction angle that conforms to a contextual geometric relationship to each contour segment; A straight line intersection calculation is performed based on the adjusted contour line segments to obtain building contour corner points, and the building contour corner points are connected to generate a final building contour line.

2. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 1, characterized in that: Step S2 includes: S21, dividing the entire point cloud data into regular grids, and defining the lowest point in each grid as a ground point; S22, constructing a triangulated network using the acquired ground points and using it as a reference ground; S23. For each triangle, count the unmarked points that fall within the triangulated area, and calculate the distance dist of each point to the nearest triangle, as well as the angle α between the line connecting the point and the vertex of the triangle and the triangle face; S24, if the distance dist of the unmarked point is greater than the distance threshold dist th And the angle α is greater than the angle threshold α th , then mark the unmarked point as a non-ground point, otherwise mark it as a ground point; S25. Repeat steps S23-24 to update the triangulated network until all points are marked as ground points or non-ground points.

3. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 1, characterized in that: In step S3, the process of the semi-suppressed fuzzy C-means clustering algorithm includes: For each non-ground point, five features are calculated based on its cylindrical neighboring points, including the height difference Δh between the highest point and the lowest point, the average height h mean , plane roughness σ, planarity planness and point density density; The feature set of all non-ground points is used as the input of the fuzzy C-means classifier for unsupervised clustering. The non-ground point set is divided into two categories according to the following principle: Where, d ij is the sampling point x j The distance to the i-th cluster center, u ij is the fuzzy membership matrix, n is the number of non-ground points, and c is the number of categories; According to the geometric characteristics of building points, the category that meets the building characteristics in the two divided categories is used as the initial seed point of the building.

4. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 3, characterized in that: In step S3, region growing is performed based on the initial building seed points to extract building points, including: A1. Select the initial seed points of the building as the seed set S and establish an empty result point set R; A2. Calculate the normal vector n for each seed point M. M and height h M , and fit the plane through its neighboring points to obtain the plane equation ax+by+cz+d=0; A3. Search the set T of k nearest neighbors of point M. A4. For each point N in T, calculate its normal vector n N and height h N , and the distance d from the point N to the fitting plane of the seed point M; A5. If point N satisfies the following three conditions: The distance d from point N to the fitting plane M is less than the preset threshold d th ; Normal vector n of point N N Normal vector n to the seed point M M The angle θ between MN Less than the preset threshold θ th ; Height difference|h M -h N | is less than the preset threshold h th ; Then add point N to the result point set R and to the seed set S; A6. Repeat steps A2-A5 until all points in the seed set S are processed and finally R contains all identified building points. A7. Cluster the identified building points according to the connectivity between points, remove the point sets with too small an area, and obtain the building points.

5. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 1, characterized in that: In step S4, processing the building points using the Alpha Shapes algorithm includes: Project the building points onto the xy plane to form a two-dimensional point set; For any building point P(x P ,y P ), use KD-tree to search for the nearest neighbor point NP within 2α, select any point Q(x Q ,y Q )Calculate the center coordinates of the two circles O1(O x1 ,O y1 )、O2(O x2 ,O y2 ); For the remaining points in NP, if none of them fall within a circle calculated above, then point P is defined as a contour point; Extract the contour points from the remaining building points to form a building contour point set.

6. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 5, characterized in that: In step S4, performing straight line fitting on the building outline points and performing region growth clustering grouping includes: B1. Use the least squares method to perform a straight line fitting on each building contour point to obtain the straight line fitting error of each contour point; B2. Select the contour point with the smallest straight line fitting error as the initial seed point; B3, establishing the initial contour point group L where the seed point is located; B4. In the set of k nearest neighbor points of the seed point, select the point whose vertical distance to the fitting line of the midpoint of the current group L is less than the threshold D. th The points are added to group L; B5. Update the straight line fitting parameters of group L; B6. Repeat steps B4-B5 until no new points are added to group L; B7. Reselect the point with the smallest fitting error among the ungrouped points as a new seed point; B8. Repeat steps B3-B7 until all contour points are assigned to corresponding groups.

7. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 6, characterized in that: In step S4, determining the main direction of the building by calculating the angle between adjacent points in the contour point group includes: Group the contour points N = {p1(x1,y1,z1),…,p n (x n ,y n ,z n )}, calculate the angle α between any adjacent contour points i , and calculate the angle variance δ: a i =acos(n1,n2) n1=(x i-1 -x i ,y i-1 -y i ) n2=(x i+1 -x i ,y i+1 -y i ) Where acos represents the arc cosine function, which is used to calculate the angle between two vectors; n1 is the vector pointing from the i-1th point to the ith point; n2 is the vector pointing from the i+1th point to the ith point; is the average value of all angles; n is the number of grouped contour points; δ is the angle variance; A group of contour points with an angle variance less than a set threshold is selected, and the fitted straight line is defined as the main direction of the building.

8. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 6, characterized in that: In step S4, merging the contour point groups using a hierarchical clustering strategy includes: C1. Calculate the straight line fitting error after merging any two adjacent groups of contour points; C2. If the minimum merge error is less than the set threshold, the two groups of contour points corresponding to the minimum merge error are merged; C3, repeat steps C1 to C2 until the minimum merging error is greater than the set threshold, then stop merging; C4. Perform straight line fitting on each final merged contour point group to obtain the building contour segment.

9. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 1, characterized in that: Step S5 includes: S51, establishing an undirected weighted graph structure based on the geometric relationships of the acquired building contour segments, where each node of the graph represents a contour segment, and the edges between the nodes represent the connection relationship between the segments; S52. For each contour segment, calculate its data penalty value under the two constraints of the main direction and the original direction, and construct the node weight of the undirected weight graph; S53. For any two adjacent contour line segments, calculate a smoothing penalty value of their direction angle difference, and construct an edge weight of an undirected weighted graph; S54. Based on the constructed undirected weighted graph, the minimum cut principle is used to perform graph segmentation optimization and assign the optimal direction angle to each contour segment; S55. Refit each contour line segment according to the optimized direction angle to obtain an adjusted contour line segment; S56. Calculate the intersection points of adjacent contour line segments as building contour corner points, and connect all corner points to generate a closed building contour line.

10. The method for extracting building contours from point clouds based on contextual geometry global optimization according to claim 9, characterized in that: The calculation formula of the data penalty value in step S52 is: Where D p (l p ) is the data penalty value, l p is a straight line fitted by least squares using the main directions and grouped contour points, p i is the i-th contour point belonging to the contour segment p, dist(p i ,l p ) is point p i To the straight line l p The vertical distance, a lp 、b lp 、c lp is a straight line l p The parameter x i ,y i For p i The plane coordinates of , n is the number of contour points contained in the contour segment p; The calculation formula of the smoothing penalty value in step S53 is: V p,q (l p ,l q )=exp(-|α p -a q |) Where, l p and l q For two adjacent contour segments, α p and α q are line segments l p and l q direction angle.

Images