Roof Reconstruction of Aerial Point Cloud Based on BPPM Plane Segmentation and Energy Optimization

In this article, a novel reconstruction method for aerial LiDAR point cloud building models is proposed to obtain valid roof building models. There are two problems in the reconstruction. First, in the process of primitive segmentation, due to the uneven density of the point cloud, there is the problem of over- or undersegmentation of the plane, resulting in the inability to extract concise and suitable building planes, which, in turn, affects the reconstruction topology. Second, in the model construction process, due to the variety and complexity of building structures, obtaining regular, compact, and topologically correct surface models from sparse and noisy point clouds is still a challenge. To address the first problem, first, the initial primitives are obtained using an improved multiresolution supervoxel-based region growing segmentation algorithm. Then, a new progressive primitive fusion algorithm Belief Propagation Primitive Merge is proposed to optimize the fragmented primitives. For the second problem, first, the Corner KLine regularization algorithm is proposed to obtain the building footprints. Then, the height map is constructed from the point cloud to extract the polyline of the building boundaries and deduce the vertical planes. Finally, a new energy function is proposed to encourage the selection of the recommended combination of model planes to obtain a compact and valid roof reconstruction model. Experiments are performed on different roof point clouds to quantitatively evaluate the proposed method, and qualitative experiments are conducted with comparative experiments to confirm the effectiveness of the algorithm.

The point clouds obtained by different LiDAR devices are different in processing and application scenarios, for example, terrestrial laser scanning and mobile laser scanning only obtain point clouds on the outside of street buildings without being able to obtain the overall structure of the building. An aerial laser scanner (ALS) is capable of achieving the roof structure and covers large areas; thus, the ALS point cloud data are selected for modeling in this article. However, ALS data have sparse, noisy, unevenly dense, and lack of topology, resulting in the loss of planar primitives and edge information during reconstruction, which makes it more challenging to reconstruct a reliable building model. Therefore, this article proposes a new roof reconstruction method, which is based on the data-driven approach. The goal of the algorithm is to construct an automated reconstruction method to generate a valid 3-D model that is watertight and regular, and the main contributions are as follows.
1) In the plane segmentation process, a novel rough-to-fine ALS point cloud clustering algorithm based on supervoxel is proposed. This includes the progressive primitive fusion model Belief Propagation Primitive Merge (BPPM), which transforms the primitive merging problem into a probabilistic inference problem to eliminate the primitive fragmentation produced in rough clustering. 2) In the model boundary extraction process, the Corner KLine (CKL) regularization algorithm is proposed to obtain the results of the building boundary segmentation, and then, the building boundary is regularized using the structural laws of the building (parallelism, collinearity, and orthogonality) to adjust the orientation as a constraint for generating the model sidewalls. 3) In the model reconstruction process, a new energy term is introduced to enhance the model topology and details based on a hypothesis and selection-based framework for building reconstruction [5]. Then, the roof building model of LOD2 [6] is obtained by selecting hypothetical surfaces. The rest of this article is organized as follows. Section II reviews the three popular algorithms for ALS roof reconstruction and analyzes the advantages and drawbacks of each algorithm. Section III describes the algorithm of this article in detail from the two parts of primitive segmentation and model reconstruction. The experimental results, the qualitative and quantitative evaluation, and discussion are described in Section IV. Finally, Section V concludes this article.

II. RELATED WORK
Currently, ALS building reconstruction consists of datadriven, model-driven, and hybrid-driven approaches. In this section, the current state of the research is presented in terms of the three approaches mentioned above.

A. Building Reconstruction Based on the Data-Driven Approach
The data-driven approach is a bottom-up reconstruction method, which first uses planar segmentation to obtain the topology of the primitive. Geometric elements, such as key points and internal and external lines of the primitive, are then extracted and finally triangulated to form the building model. Zhou and Neumann [7] proposed an algorithm to generate a 2.5-D building model, which extends the dual contour to the 2.5-D method to realize the optimization of surface and boundary. Although this method can generate building models with arbitrary shapes, the model quality needs to be optimized. Sampath and Shan [8] introduced a fuzzy k-means clustering algorithm to segment the roof planes and then used the distance threshold to obtain the adjacent matrix of the roof topology. Finally, a regularization method is proposed in the modeling process to generate topologically building models. In order to optimize the 2-D topology of the roof, Yan et al. [9] adopted a novel snake algorithm, which imposes multigeometric constraints in the energy function to optimize the topology, enforce parallelism on the 2-D topology, and adjust the edges of the twisted building model. The final model showed that the method is effective. Kim and Shan [10] proposed a roof reconstruction method that includes planar segmentation and roof modeling, which used a level set function to segment the roof primitives and estimate the topological connections, and then obtained roof models by intersecting lines of adjacent roofs. Poullis and You [11] introduced a large-scale point cloud reconstruction method, which includes plane segmentation, fitting, boundary regularization, and 3-D modeling. In addition, a P2C planar primitive clustering was proposed in [12], and an energy optimization algorithm for boundary constrained building modeling was proposed to obtain a nonoverlapping and watertight 3-D model. Then, in 2019, Poullis [13] proposed a novel large-scale reconstruction framework, which assumes the similarity distance as an adaptive computational parameter of the Weibull distribution and combines attribute information into tensor clustering, demonstrating the validity of the model results. Zhang et al. [14] proposed a framework combining classification and reconstruction, which first classifies point clouds using a training-free neural network (ReLU-NN). The point sets are then integrated by edge-aware point resampling. Finally, the reconstruction model is obtained using the 2.5-D dual-contour method. Cao et al. [15] introduced a parallel roof reconstruction framework, which uses a density clustering algorithm to separate individual buildings, then segment and regularize the building roof surface, and finally reconstruct the roof using topological and geometric message. Chen et al. [16] proposed a multi-LOD reconstruction framework, which first uses Voronoi subgraphs to obtain the primitive boundaries and then optimizes the roof boundaries. Finally, the building is modeled using a topologically aware method to get the model. In order to obtain a complete model of the building, Awrangjeb et al. [17] proposed a new roof reconstruction method, which uses adjacency matrix intersection line analysis to identify and fill in missing planes. Li et al. [18] proposed a TIN-based region growing method for clustering planes in order to solve the problem of uneven point cloud density affecting model quality, which projected to the 2-D grid map to construct the labeled graph optimized and regularized boundaries. Finally, compact building models were obtained. Hu et al. [19] proposed a method for the 3-D reconstruction of buildings combining semantic and topological information, using energy function to extract planar primitives, then constructing primitive topology based on neighbor distance and similarity. Finally, a progressive decomposition method is introduced to reconstruct the building to obtain the effective building models with semantic structure. Huang et al. [20] proposed a novel method for reconstructing watertight 3-D building models, which is based on extended assumptions and a selected polygonal surface reconstruction framework to obtain the final model. Xie et al. [21] proposed a rule-and hypothesis-based reconstruction method that introduces pairwise constraints, triplet constraints, and proximity constraints to remove the fuzzy region topology. Liu et al. [22] proposed a GUIDED-based planar segmentation, topology optimization, and corner point optimization method for building reconstruction, which improves the geometric accuracy of the model. Li and Wu [23] proposed a CSG-BREP tree-based approach to generate geometric building models. And a new optimization strategy for 3-D spatial cell selection was developed to extract 3-D surfaces from topological relations.
Although the data-driven approach is commonly used in modeling, it has drawbacks: the model result is sensitive to the density and noise of the point cloud and relies on the results of primitive segmentation; the computational process consumes relatively high time costs.

B. Building Reconstruction Based on the Model-Driven Approach
The model-driven approach [24] is a top-down reconstruction method, which first builds a library of model primitives (flat, shed, gabled, hipped, and pyramidal) and then processes the building point cloud to find the model and parameters that match the library. The final model is approximated by the point cloud.
Poullis and You [25] proposed a novel interactive building reconstruction method, which reconstructs linear and nonlinear surfaces using the symmetry constraints of parameterized geometric primitive. Huang et al. [26] introduced a generation statistical-based model library method for roof reconstruction to address the problem of model irregularity and incompleteness, which utilized the reversible jump Markov chain Monte Carlo to select primitives. Henn et al. [24] proposed a model-driven automatic reconstruction method, which combines RANSAC and support vector machine algorithms, and the extracted features improve the accuracy of model selection. Zheng and Weng [27] proposed a model-driven reconstruction method for irregular planes, which uses a decision tree based on physical morphological features to classify point clouds into seven classes and uses static moment equations to calculate the principal directions to reconstruct the roof. Zhang et al. [28] proposed an optimal model fitting approach to address the problem of overreliance on the roof plane segmentation accuracy in the modeling process, which transforms the reconstruction problem into an optimization problem that minimizes the distance between the point cloud and the TIN model to generate semantic building models. Li and Shan [29] proposed a multiprimitive reconstruction method to get regular building models; the method utilizes RANSAC to perform primitive segmentation and then combines the holistic primitive fitting and 3-D Boolean operations to generate the building model.
Although the model-driven approach produces a more regular roof model than the data-driven approach, it is limited by the completeness of the model categories in the model library and has difficulty matching building roofs when the number of models in the model library is insufficient.

C. Building Reconstruction Based on the Hybrid-Driven Approach
The hybrid-driven approach is a combination reconstruction method of data-driven and model-driven approaches, which divides complex roofs into subsets, and each subset gets a matching parametric model in the model library. Another approach is to construct a roof topology graph (RTG), which transforms the reconstruction problem into a matching problem between the topology graph and the model library objects. This method not only combines the advantages of the two reconstruction methods mentioned above, but also reduces the shortcomings of each.
Xu et al. [30] proposed a hierarchical roof topology to obtain the relationship between primitives, which was first employed to classify planes. Then, parallel perpendicular colinear relationships within and outside the model are detected and weak connections are distinguished by a hypothesis-verification approach, resulting in a 3-D model. Lafarge et al. [31] proposed an object-based method to reconstruct the model, which first obtains the footprint of the building. Then, the adjacent matrix and height are utilized to regularize the building, and the building models are obtained. Xiong et al. [32] introduced a graph edit dictionary method to automatically find and correct graphs with wrong roof topologies, which puts the manual-edited topology graph into the dictionary, and then obtained the best model with a heuristic method. Then, in 2015, Xiong et al. [33] decomposed the semantic structure building according to the basic elements and physical meaning of the RTG and obtained the building model by combining the buildings according to CSG-style pairs of subgraphs. Lin et al. [34] proposed a reconstruction method for low-rise buildings, which first partitions the point cloud semantics into different classes and then transforms the reconstruction problem into a combination problem of multiple levels of small structural blocks. Hu et al. [35] introduced a roof attribute graph method to represent roof topological relationships and proposed a decomposition and refinement process that follows gestalt principles to reconstruct roofs generating a complete and visually distinctive roof model. However, the computational complexity and the automatization of the reconstruction algorithm are still the difficulties and challenges of the hybrid-driven reconstruction method.

A. General Overview of Method
The workflow design of roof reconstruction proposed in this article is generally divided into two stages: roof primitive plane extraction and roof model reconstruction. And the overall process results are shown in Fig. 1.
1) Primitive plane extraction: First, multiresolution supervoxels are constructed based on aerial LiDAR building roof point cloud, and rough segmentation is performed using an improved region growing algorithm. Second, a progressive roof primitive merging optimization algorithm is proposed. Finally, the primitive boundaries are optimized using the energy function. See Section III-B for details. 2) Roof model generation: First, the footprint of each building is extracted and regularized using CKL. Second, the vertical side walls of the buildings are inferred. Finally, an energy optimization function is used to generate the building models. The details are provided in Section III-C.

B. Extraction of Roof Plane Primitives
In general, in order to obtain accurate roof planar primitives from airborne LiDAR building roof point clouds for subsequent reconstruction, a novel roof planar primitive segmentation optimization method is proposed, which consists of the following steps. First, a boundary constrained multiresolution supervoxel segmentation algorithm is applied to get the initial supervoxels (shown in our previous work [36]), Second, an improved supervoxel region growing algorithm is proposed to obtain initial planes. Finally, a novel BPPM model is proposed to progressively merge primitives belonging to the same plane. The pseudocode for roof primitive extraction is shown in Algorithm 1 , and the details are described as follows.
1) Multiresolution Supervoxel Generation: Due to the disorder of point clouds, it is necessary to construct spatial topology. In our previous work [36], multiresolution supervoxels have the advantage of being capable of preserving boundaries and providing a highly approximate representation with less data. Therefore, supervoxels are used as the basic unit of the planar primitive extraction.
The initial supervoxel is obtained by BPSS and energy function optimization, where the energy function in (1) includes distance constraint and representative point constraint. Then, the energy optimal solution (NP-hard problem) is solved by the fusion and exchange algorithm, and the detailed description of algorithm can be found in [37] min where E(Z) is the energy function, z ij means point p j is represented by p i , which p i is a representative point, and p j is the neighboring points of the representative point. D(p i , q j ) represents the distance between the representative point and neighboring points. And λ is a parameter to trade off distance and the number of supervoxels. C(Z) is the number of representative points, which penalize the number of supervoxels. However, when this method is applied to ALS point cloud, this method cannot detect low-density planar area, resulting in inaccurate clustering results. Therefore, the local features and residuals of each supervoxel are calculated to obtain the lowdensity supervoxels, as shown in (2). Finally, the low-density supervoxels are resegmented with low resolution, and the results are used as the basic unit for planar clustering where P λ , L λ , and S λ represent local geometry features.

2) Supervoxel-Based Region Growing Rough Segmentation:
In this study, the improved region growing algorithm is applied to cluster supervoxels. First, the supervoxel with the smallest residual is selected as the initial seed. Then, if the adjacent supervoxels satisfy the combined feature conditions, merge them. Finally, when there are no more supervoxels that satisfy the conditions, merging is terminated. Since the roof point cloud has parallel planes, only relying on the angle threshold to segment the point cloud will cause undersegmentation in parallel planes with inconsistent heights. In order to segment two planes that are parallel yet not connected, as shown in Fig. 2, the proximity between the supervoxels is included in where S dist represents the similarity distance between adjacent supervoxels. D ij and A ij are the proximity and angle difference of supervoxel i and j, respectively. λ 1 and λ 2 are weights to Algorithm 1: Roof Plane Segmentation Algorithm.
measure proximity and smoothness where S i and S j are the angles formed by the supervoxel normal vector and the centroid connection vector, respectively. θ is determined by the sigmoid function of the difference between S i and S j [38].
3) BPPM Algorithm: After rough segmentation, the initial plane primitive patches are obtained. Nevertheless, many primitives belonging to the same plane remain unmerged during the region growing of supervoxels due to uneven point cloud density and difficulties in setting the supervoxel resolution, resulting in incomplete planar primitives for reconstruction. To fuse primitives in the same plane and improve the segmentation accuracy, a new progressive fusion algorithm is proposed. Inspired by [39], the fusion problem of primitive fragmentation is transformed into a probabilistic inference problem in this article. The probability distribution obtained by each primitive node in the model is iteratively fused by passing messages to neighboring primitive nodes, making the probability distribution of each node eventually maintain a stable state. The detailed description of the BPPM algorithm is presented from the following three points: graph model building, similarity message propagation, and model calculation principle. 1) Graph model building: In fact, belief propagation is a special Markov model. The primitives obtained by coarse segmentation correspond to each node of the graph model, which consists of an explicit node and an implicit node. Concretely, explicit nodes correspond to internal features in the primitive plane, and implicit nodes correspond to similarity features between primitives. In essence, the beliefs are the combined internal and external features evaluated between the current primitive and the neighboring primitives. Meanwhile, according to the primitive size (the number of points contained in the primitive), this article classifies the primitives into two categories, complete (large size) and fragment (small size), corresponding to "c-node" and "f-node," respectively. The mapping relationship between the model and the segmentation result to be optimized is shown in Fig. 3. 2) Similarity message propagation: Substantially, the belief is to quantify the similarity between primitives in the same plane, aiming to realize the fusion of the same planes and retain the differences of different planes. In each fusion, the algorithm transforms the planar information of the primitive and local feature information into a message passed to the neighboring primitives, affecting the belief of the other primitives. The propagation of node messages is shown in Fig. 4, which is mainly divided into propagation between c-node and c-node and between c-node and f-node. Note that the propagation terminal of the message ends at the terminal primitive of the single building.

3) Model calculation principle:
In the proposed algorithm, the belief of current primitive is first calculated (product of all messages). Then, the neighboring primitive j is prefused [see (9)], and the mean squared error (MSE) of the fused primitive is calculated. Eventually, if the new primitive satisfies the condition, the prefusion result is accepted and updates the graph model to iterate the fusion again. After each primitive belonging to the same plane is merged, the iteration is terminated. Essentially, the maximum a posteriori estimate is obtained in each iteration of the BPPM model to merge the adjacent primitive with the strongest feature similarity. The belief obtained by the primitive is the product of messages, which represents the similarity of local neighboring primitives and is updated in each iteration. Assuming i is the current primitive, the belief is calculated as where φ i (x i ) is likelihood between internal planar features and external similarity, N is the number of adjacent primitives of the ith node, m ij (x i ) is the message passed into primitive i, and k represents normalization constant where p represents the points contained in the ith primitive, P g is the fitting plane of the points in the ith primitive. D(p m , P g ) and σ i represent the distance and variance from point m to the corresponding primitive plane, respectively. Then, the formula of the incoming message from the node is as follows: where m ij (x j ) represents the propagation message from primitive j to i during fusion and ϕ ij is feature similarity between primitives where G i j , D ij , and A ij represent feature similarity, proximity, and center normal vector angle between primitives i and j, respectively arg min In this way, the oversegmentation primitives are merged by the BPPM algorithm, and the innovations of BPPM model are as follows. 1) BPPM divides primitives by size and progressively merges primitives belonging to the same plane using internal and external features of the primitives. 2) A premerge step is determined using MSE to determine whether the merged primitives conform to the same plane.
3) The iterative termination condition for the progressive optimization is that all primitives belonging to the same plane are merged.

4) Optimization of Planar Primitives Based on Energy Minimization:
Although the BPPM algorithm merges fragmented primitives and reduces the oversegmentation rate, the primitive planes still have irregular boundaries. Hence, the roof planes are optimized by the energy minimization algorithm mentioned in [40], and the energy function is shown in (10), including distance cost and boundary cost where D(A) represents the distance term and G(A) is the boundary term where d(x i ) represents the distance from point x i to plane P k and g(x i ) represents the smoothness term of the point x i in the region of the neighbor N i . The label of point x i in term g(x i ) is assumed to be k. If the label of A xj is equal to k, then η = 1; otherwise, η = 0. And the graph-cut-based algorithm is used for energy minimization. The results of the building segmentation process are shown in Fig. 5.

C. Building Model Reconstruction
After obtaining building primitives, a novel roof reconstruction method is proposed to obtain the polygon model of roof buildings that conform to human intuitive vision. The overall flowchart for the reconstruction is shown in Fig. 6, which consists of three main steps: building footprint extraction, building vertical wall extraction, and building model generation.

1) Extraction of Building Corners:
In order to obtain the regular building model, the alphashape algorithm is utilized to obtain the initial building boundaries. However, due to the uneven density of the point cloud and the alphashape radius setting, the corner point extraction process results in underdetection [41]. In addition, using curvature alone can only locate ambiguous corners, but cannot obtain the precise number of corners, because the angle detection result will change with the curvature value. To obtain the simplified boundary corners, a combination of curvature points and line intersection angle is first utilized. The boundary points between each of the two corner points are then determined sequentially using the line residuals to determine whether they contain missing corner points. Finally, an improved missing corner detection algorithm is introduced to find the missing corners. a) Initial corner detection: The initial corner consists of the curvature points and the line intersection angle points. Fig. 7(b) displays the building roof boundary points obtained from alphashape, and Fig. 7(e) shows the extracted initial As the boundary points obtained are sequential, each boundary point is first traversed, and its first three and last three neighbor points are selected. Afterward, the first three and last three points are fitted with PCA to find the intersection angle, which is obtained by fitting the line from the first three points to intersect the line from the last three points. Note that the intersection angle is not the local internal angle of the boundary, as shown in Fig. 7(d). The green point in Fig. 7(c) is the result of the corner obtained from the intersection angle. Besides, curvature is an intuitive local feature in corner point detection. Since most of the roof point clouds have high curvature at the corner points of the boundaries, the corner points are chosen to be obtained by calculating the local curvature of the boundary points and by threshold selection, as shown in Fig. 7(c). And the curvature calculation is shown as where C r represents the curvature at a point on the boundary and λ min is the minimum eigenvalue obtained from PCA. b) Missing corner detection: Due to the boundary point extraction threshold setting and the unevenness of the point cloud density [see Fig. 7(f)], there are still many undetected corners in the initial extracted results. Therefore, in order to find the missing corners, it is necessary to first detect the line segment containing the undetected corners. The detailed procedures are as follows.
1) Calculate the residuals and the number of points of the boundary points between adjacent corners; if the number of points is greater than 3 and the residual is greater than the threshold, then proceed to steps 2 and 3. 2) Line segments that are detected as potentially containing missing corner points are divided into two categories. The first category is (x length > y length ), where the length of x is greater than the length of y between corners containing boundary points; then, x min and y min are subtracted and projected onto the coordinate axis to obtain the new scatter point E 1 . The second type is (x length < y length ), where x min and y min are subtracted in each point, and the coordinates of each point are dropped to obtain the scatter E 2 . The final fitted curve is obtained by fitting the above two types of scatter points separately using the smooth spline (smoothed cubic spline) method. As shown in Fig. 8.
3) The spacing between adjacent extremes obtained from the fitted curve should be greater than 0.8 in the horizontal direction and greater than d max in the vertical direction, where d max is determined by the density of the point cloud. For the detailed analysis of the threshold statistics, see [41]. 2) Building Footprint Regularization: Extracting the roof corners, which are directly used as external boundary points of the building, usually suffers from jagged boundaries leading to model artifacts. Therefore, in order to obtain a building model with regular boundaries, the boundaries need to be regularized. Therefore, this article proposes a CKL-based boundary point regularization method, which includes two steps of clustering and adjustment.
More specifically, for the boundary B p of a building, CKL first obtains the number of building boundary clusters based on the number of corner points. Then, the coordinates of the midpoint of each two corner point clusters are calculated clockwise as the seed points of the line segmentation, where the line segment is indicated by L = {L 1 , L 2 , . . ., L k }, and each set of lines contains points p = {p n 1 , p n 2 , . . ., p n m } belonging to the line. The line segmentation process assigns the points according to distance and angle; if the distance and angle to the fitted line are less than the threshold, then they are accepted. If there are no unassigned points in the initial set of boundary points, then the line set L is accepted; otherwise, the line set L is cleared and the seed points are reselected for line set assignment. CKL is an iterative algorithm of line segmentation where the condition for ending the loop is whether the initial set of boundary points is empty. The pseudocode of CKL is shown in Algorithm 2.
Note that regularization and orientation adjustments are made after obtaining multiple clusters of line segments. For each line segment cluster after clustering, its average direction is calculated, and then, the line segment directions of the clusters are adjusted so that they have the same average direction. The method is to construct the minimum boundary rectangle by the MBR algorithm [42]. Then, traverse each line segment cluster after CKL segmentation, and if the line segment cluster is parallel (i.e., the difference in angle is less than 10 • ) or orthogonal to the two directions of the smallest rectangular box, the clusters of line segments that are colinear and the clusters of line segments that are approximately perpendicular are approximated and adjusted. The final regularized building footprint effect obtained after the adjustment is shown in Fig. 9.
3) Building Vertical Wall Extraction: The airborne LiDAR point cloud data are obscured during the pushing and scanning process, making it impossible to obtain information on areas such as vertical walls. However, the vertical walls can visually reflect the connection relationship with the roof, where the vertical walls include not only the vertical walls generated by the building footprint but also the internal vertical walls of the building. Therefore, in this article, the algorithm proposed by Huang et al. [20] is used to extract the building polyline by constructing the height map through the external footprint constraint of the building and obtain the internal and external walls of the building to obtain the topologically correct building model [43]. The detailed description is as follows.
First, the point cloud is constrained to the extent of the building's outer boundary footprint, and the 3-D roof point cloud data are projected into 2-D. The 2-D Delaunay triangle is used to construct a triangular irregular grid to generate the 2-D raster height map. After that, the raster height map is morphologically operated to fill the holes. And find contour points where there is variation in height values, which means that a set of discrete contour pixel points is extracted from the height map using the Canny detector as an initial estimate of the vertical plane. Finally, the regular vertical wall contour lines are obtained using the regularization method for the obtained contour points. And extrude each contour line to get the inferred vertical walls to constrain the plane selection. The building vertical plane generation process is illustrated in Fig. 10(b)-(e). It can be seen from the figure that the deduced vertical walls are able to solve the model reconstruction problem caused by incomplete building planes and finally obtain valid 3-D building models.

4) Building Model Generation:
The building model generation method proposed in this article is based on the hypothetical selection framework [5]; however, the model generation method in this article is based on the selection of planar elements generated by intersecting the vertical and segmentation planes with each other. Unlike [5] and [20], this article uses BPPM to generate the roof planes to calculate the intersection, while no manual footprint is required as input. This not only ensures the realistic connection between the roof model sidewalls and the roof model planes, but also will make the model more regular.  In addition, to encourage the combination of realistic building models that match the physical world, new energy terms are introduced. The overall energy function of the model face is shown as where E d is the fit term for the roof model surface and E c is the model complexity term for the planar structure. The details are described in [5]. The new resource term E n proposed in this article is described as follows. For roof point clouds, roof planes with complex structures usually have large differences. The reconstruction results show multiple roofs; therefore, a new elevation-preferred resource formulation E n is introduced to physically weigh the higher roofs and the other roof planes. In this article, we use Gurobi to minimize the energy formulation to obtain the final plan combination of the building, and solving the hard constrained plan combination is a binary linear program problem, and the results of the comparison of the energy optimization are shown in Fig. 11.
where |F | denotes the number of hypothetical planes; z i represents the z-axis coordinates of the hypothetical plane center; z mean , z max , and z min represent the average height, highest height, and lowest height of the roof hypothetical plane, respectively; and α, β, and γ are the weighting factors.
In summary, the reconstruction algorithm in this article differs from the methods of [5] and [20] as follows.
1) In the boundary regularization procedure, a new CKLbased algorithm is proposed in this article to obtain the regular building footprint. 2) To guarantee the model topology correctness, a new energy term is proposed to constrain the building model plane selection. 3) Without inputting footprint as the model prior, the model is directly reconstructed.

A. Experimental Setup
Our proposed reconstruction framework is effective for rooftop buildings. Experiments were conducted on two public datasets, and the experimental results were evaluated quantitatively. The results of the algorithm are also compared qualitatively with the comparison algorithms in the planar primitive extraction process and the model reconstruction process, respectively. This section focuses on the experimental datasets, the evaluation indicators, and the parameter settings in each step. The experiment is carried out on a computer with Intel(R) Core(TM) i7-4510U CPU @ 2.00 GHz and 16-GB RAM hardware. The algorithm is implemented in C++ through the CGAL and PCL libraries.

B. Description of Datasets
In this article, three datasets are used. The specific parameters of the datasets are described in Table I, which includes specific information about the point clouds, including dataset name, collection sites, density, equipment, and the acquisition website. Figs. 12-14 show the point clouds for the three datasets. The datasets are described as follows.
1) The first dataset is the AHN3 dataset, which was collected in the Netherlands with an average point set density of 8. This dataset contains semantic labels for buildings, trees,   etc., which can be directly obtained by obtaining individual building point clouds. Both the raster data and the raw point cloud data can be downloaded through PDOK and NationalGeoregister.
2) The second dataset is the Dayton Annotated LiDAR Earth Scan (DALES) dataset [44], which was obtained by Riegl Q1560. The location of the captured site is the city of Surrey, Canada, which includes categories such as ground, vegetation, buildings, and power lines.
3) The third dataset was captured in Vaihingen, Germany, courtesy of the German Society for Photogrammetry, Remote Sensing and Geoinformation, and was acquired by the Leica ALS50 unit flying at an altitude of 500 m. The average density of the point cloud is 4 points/m 2 . Three areas are selected for reconstruction experiments, which are all in the ISPRS benchmark for "Test Project on Urban Classification and 3D Building Reconstruction" [45], [46].

C. Evaluation Metrics
To quantitatively evaluate the method proposed in this article, two procedures are evaluated from primitive segmentation and model reconstruction. According to [47], planar segmentation uses the following evaluation indicators: completeness (C m ), correctness (C r ), and quality (Q l ), i.e., The RMSE of the boundaries is calculated in the reconstruction process. The reconstructed model topology is evaluated on the ISPRS data using N 1:M (oversegmentation), N N :1 (undersegmentation), and N N :M (the plane both over-and undersegmented). In addition, the models are evaluated using completeness (Comp), correctness (Corr), and quality (Q).

D. Parameter Settings
The parameters of the algorithm are described in detail in terms of both segmentation and reconstruction processes, and we illustrate how their values were chosen. The values of the algorithm parameters are shown in Table II. There are six significant parameters in the planar segmentation procedure, i.e., R s , S dist , λ 1 /λ 2 , T f , T M , and Δd. R s is the supervoxel resolution for the supervoxel segmentation process; S dist is the discrimination threshold for the clustering of similarity features between supervoxels in the rough segmentation procedure. λ 1 and λ 2 are the weight values for proximity and smoothness, respectively. T f is the number of points' threshold for the BPPM model to classify primitive types. T M is the MSE that limits the fusion of the same plane during planar segmentation. Δd is the distance threshold from the noise point to the plane for the planar energy optimization. Subsequently, the reconstruction process has four important parameters: the residual threshold from point to the fitted line in the missing corner point detection process and the CKL line segmentation process (i.e., ε), the horizontal threshold for missing corner filtering process (i.e., d max ), and the distance threshold from point to line in the CKL process (i.e., T D ). And in the modeling process, λ n is the weighting factor for the new energy term.
Notice that the parameters in the table with empty values require different settings due to the different point cloud densities. During the segmentation process, the supervoxel resolution R S of the AHN3 and DALES datasets is set to 2, and to 3 in the Vaihingen dataset. In the reconstruction process, the parameters d max and T D are set to 0.3 for both the AHN3 and DALES datasets and to 0.6 for the Vaihingen dataset. The other remaining parameters are all set equal in three datasets and are shown in Table II.

1) Experimental Results of Planar Segmentation Comparison:
In order to compare the performance of the planar segmentation algorithm proposed in this article with other state-ofthe-art methods, the following algorithms are used. QTPS [48] is a novel ALS point cloud planar segmentation method based on supervoxel and quasi-a-contrario theory. HCBR [40] is a hierarchical clustering algorithm with boundary relabeling. SVGS [49] is a multifeature graph segmentation algorithm based on supervoxels. RG is a classical algorithm for planar segmentation. Thus, these comparative algorithms are chosen for evaluation.
The results of the segmentation algorithm in this article and the comparative segmentation algorithm are shown in Fig. 15, where the four rows represent four different buildings, and the different columns represent the results of the different algorithms for building segmentation. Specifically, the first column shows the experimental results of this article, the second column shows the results of the SVGS segmentation algorithm, the third column shows the results of the RG segmentation algorithm, the fourth column shows the results of the QTPS algorithm, and the last column shows the results of the HCBR algorithm segmentation. Hereby, the advantage of this article's algorithm is that it balances completeness and segmentation accuracy, and is able to effectively handle oversegmented regions, and guarantees the accuracy of the boundaries. In the column of Fig. 15(b), it can be observed that the SVGS algorithm has undersegmentation in the complex region of the planar structure and cannot obtain well results in the region of planar primitive boundaries. In Fig. 15(c), it can be observed that the RG algorithm suffers from some point cloud data loss in the planar intersection region, which affects the completeness of the point cloud data. In the black box in Fig. 15(d), it can be observed that the QTPS algorithm also has some under-and oversegmented areas. Moreover, the results obtained by the HCBR algorithm in Fig. 15(e) have weak segmentation completeness and holes in the building point clouds within the black boxes, which, in turn, influence the topological accuracy of the subsequent reconstruction steps. Thus, the above results show that the segmentation algorithm in this article guarantees accuracy and completeness.  Compared to traditional algorithms, the proposed method performs better in areas of buildings with uneven density. As shown in Fig. 16, the algorithm in this article does not suffer from serious oversegmentation in the buildings with uneven density in the red box, which emphasizes the advantage of the BPPM algorithm.
In order to emphasize the impact of the planar segmentation algorithm proposed in this article on the reconstruction process and to highlight the advantages, the ablation experiments are designed. The segmentation results are first obtained using different algorithms in the planar segmentation process, and then, the obtained planes are put into the reconstruction methods of this article separately to obtain model results for comparison, as shown in Fig. 17.
Typically, Fig. 17(a) shows the planar segmentation results and model reconstruction results of the algorithm in this article. It can be seen that the BPPM algorithm in this article is able to eliminate oversegmented planes and has the advantage of providing well-defined topological planes for the reconstruction process, because, on the one hand, the BPPM algorithm is a probability reasoning optimization model that guarantees the merging of planes with similar features; on the other hand, the MSE of the planes is used as the bound for iterative fusion. Fig. 17(b) shows results of the region growing algorithm, which has undersegmentation in some narrower planes [in the red box of the segmentation results in Fig. 17(b)], thereby causing the reconstruction process to loss topological planes. Fig. 17(c) shows the results of SVGS, a supervoxel-based planar segmentation algorithm. It can be seen that the existing supervoxel-based segmentation algorithm, when applied to the ALS roof point cloud, fails to achieve satisfactory segmentation results, with serious over-and undersegmentation, resulting in the incorrect topology of the roof plane during the reconstruction process, leading to invalid building models. Fig. 17(d) shows the result of planar segmentation based on the P-linkage [50]. From the segmentation result, it can be seen that the planar oversegmentation is serious, and part of the point cloud is missing, which makes the topology of the model reconstruction result also invalid, and only a few incorrect roof planes can be obtained. Fig. 17(e) shows the result of the planar segmentation and model of the HCBR algorithm. From the result, the completeness of the point cloud after segmentation is poor, and the segmentation makes the point cloud missing and the primitive missing, which makes the hypothetical planes of the roof provided before the plane selection process is incorrect and, thus, makes the model obtained in Fig. 17(e) invalid. Frankly speaking, the reconstruction step proposed in this article places further demands on the planar segmentation algorithm, requiring both completeness and accuracy of the segmented point cloud data. However, if the segmented data are not complete, the reconstruction process will affect the regularity of the external boundaries of the model and may also cause the reconstructed model to lose topology at the roof. The model topology error in Fig. 17(c) occurs if the plane is undersegmented, which means that the model is partially missing in the red box. Oversegmentation of the plane also leads to model topology errors.
In summary, the segmentation algorithm in this article balances completeness and accuracy and obtains superior model results. It is demonstrated that compared with other existing algorithms on ALS point cloud data, the planar segmentation algorithm in this article has the advantage of being able to provide reliable topological planes for the reconstruction process.
2) Experimental Results of Model Reconstruction: In order to analyze the model reconstruction algorithms, model reconstruction results on three ALS point cloud datasets are shown and analyzed in this section, including model reconstruction results for large-scale point clouds and individual building point clouds, as well as model results for comparative algorithms.
The large-scale roof reconstruction results for the DALES dataset are shown in Fig. 18, where Fig. 18(a) shows the raw LiDAR roof point cloud, Fig. 18(b) shows the result of the reconstruction of the roof model, and Fig. 18(c) shows the result of superimposing the point cloud with the line model in DALES. The blue area in Fig. 18(c) is the roof point cloud, and the orange line is the line of the building model. Fig. 18(d) shows the side view of the results of Fig. 18(c). Consequently, the results of the large-scale DALES model reconstruction show that the reconstruction algorithm in this article is able to satisfy the model detail requirements and obtain valid results.
To clearly observe the details of the building model, Fig. 19 magnifies the reconstruction results of the algorithm in this article for individual buildings, which includes the fused representation of the building model line generated by the algorithm in this article and the point cloud. Each number corresponds to three graphs; from top to bottom, the first row shows the original building point cloud, the second row shows merging the point cloud with the line model, and the third row is the reconstruction model of the proposed algorithm, which combines building lines, models, and point clouds to display. In Fig. 19, (1)- (10) are the building point cloud reconstruction results for the AHN3 dataset and (11)- (18) are the building point cloud reconstruction results for the DALES dataset. Although the regularization process direction adjustment makes the point cloud in the area inside the black box fit differently from the building model, in general, the buildings obtained by the algorithm in this article are compact and can be well fused into the point cloud. It is because the energy function of the proposed reconstruction algorithm is able to constrain the model selection process so that the details of the model fit the point cloud.
Since the original ISPRS data contain the corresponding optical images of the buildings, the optical images are added to Fig. 20 as comparison to observe and compare the correctness of the model details. The results of four typical buildings in the Vaihingen data are selected for magnification, including the original point cloud, the model, and the corresponding optical image. Although there are minor missing building details, overall, the algorithm in this article obtains compact and guaranteed topological accuracy and boundary adaptation of the building models on the low-density ISPRS dataset. Because the CKL regularization process incorporates the buildings structural laws, enabling to obtain accurate building footprints.
Moreover, the proposed reconstruction algorithm is compared with the current state-of-the-art algorithm, that is, the results are compared with the algorithm of this article on AHN3, DALES, and Vaihingen datasets using three comparison methods, City3D in 2022 [20], TINS in 2019 [18], and 2.5-D dual contour in 2010 [7], respectively. Fig. 21 shows the performance of the proposed method versus other methods, where the first two rows from top to bottom in the first column (a) are AHN3 roof data, the middle two rows are DALES roof data, and the last row is Vaihingen data. The second column (b) shows the reconstruction results of the algorithm in this article, the third column (c) shows the reconstruction results of City3D, the fourth column (d) shows the reconstruction results of 2019TINS, and the fifth column (e) shows the reconstruction results of 2.5-D dual contour. The last column (f) is the manual model result. From the comparison results, we can observe that the reconstruction algorithm proposed in this article is the strongest in terms of boundary regularity. Meanwhile, it can get the correct model topology while ensuring less number of planes, and the model is compact and effective; as we can see from Table VI, our reconstructed model has the  highest boundary accuracy. 3) Model Evaluation: In order to obtain the pixel-level model precision in the Vaihingen dataset, the DSM results are obtained in Fig 22, where Fig. 22(a)-(c) shows the reference DSM for the Vaihingen dataset and Fig. 22(d)-(f) shows the DSM results obtained by the algorithm in this article. As can be seen from the results, the building models obtained in the Vaihingen data are compact and have well-defined boundaries and better accuracy. However, some inevitable errors (false positives and false negatives) are present at the boundaries of individual buildings. Because in the regularization process, the outer contours of the building are shifted to satisfy the orientation constraints to make it more regular. Overall, the algorithm obtained superior results for the reconstructed model. In addition, the successful implementation of the 3-D model by the algorithm in this article during the reconstruction process can be observed from the pixel-level model accuracy of the Vaihingen dataset in Fig. 23, which shows the 3-D model converted to the label graph. As the   by the reference DSM and superior precision is achieved in the pixel-level model (see Table V).
Due to the lack of ground truth in AHN3 and DALES datasets, the common method for quality assessment of model results is to calculate the RMS of the surface fit error. To evaluate the accuracy of the reconstructed model, the point cloud with error-related colors is superimposed on the building models in Fig. 24, where the color bars indicate the error magnitude. Three buildings are included in Fig. 24, corresponding to models (8), (13), and (15) in Table IV. The results show that the flatness of the algorithm in this article fitted well.

4) Time Consumption:
The time consumption of the algorithm proposed in this article and the comparison algorithm in the primitive segmentation process is as follows. On average, region growth consumes the shortest time and QTPS the longest. The algorithm in this article consumes less time than QTPS and greater than HCBR because of the increased time consumption in the multiresolution supervoxel segmentation process. The time complexity of the clustering process is O (N ).
The time complexity of the proposed algorithm is O(N 2 ); and as seen from Table IV, this article consumes longer time than 2.5-D DC and 2019TINS, and comparable with 2022City3D time, because in the process of generating the building model, the method in this article infers that the process of combining the intersection of assumed planes is a combination optimization problem. However, since the algorithm proposed in this article obtains a model that balances compactness and regularity, it is effective despite the large time spent.

5) Discussion:
Herein, a roof point cloud reconstruction framework is proposed, which aims to obtain an effective roof reconstruction model. From the above segmentation and reconstruction experimental results, there are three other points relevant to the discussion of the proposed method: selection of parameters, comparison analysis, and limitations. a) Selection of parameters for the planar segmentation and reconstruction: During the supervoxel construction process, the lower initial supervoxel resolution leads to oversegmentation as the features of the clustering process deviate from the actual feature values, while the larger supervoxel resolution leads to undersegmentation in the clustering procedure, such as dormers in smaller buildings. Moreover, the BPPM planar segmentation procedure proposed in this article requires the parameter T M to control the degree of merging of planar primitives. If T M is set too large, the planes are merged excessively, as shown in Fig. 25(d), resulting in nonidentical planar primitives being merged, causing the results deviate from the actual segmentation results, thus affecting the results of the reconstructing process. If T M is set too small, as shown in Fig. 25(b), it will also result in incomplete fusion of primitives, which should be in the same plane. To trade off the over-and undersegmentation of the roof plane and to obtain accurate planar segmentation results, the supervoxel resolution and T M still need to be obtained through iterative experiments.
In the process of the reconstruction, the setting of the threshold value affects the corner extraction results and, thus, the reconstruction quality. The higher threshold value will result in excessive corner point extraction and redundancy. The lower threshold value will produce the situation of missing corner points. Therefore, the choice of threshold value will directly affect the subsequent linear segmentation results. Thus, in this article, the residual threshold from point to the fitting line is set to 0.1 by trial and error, and the intersection angle range is 45-135 • .
b) Comparison analysis of the segmentation and reconstruction: The planar segmentation algorithm proposed in this article uses multiresolution supervoxels as the basic unit instead of using a combination of scatter and supervoxel clustering, because such combination affects the efficiency of planar clustering, for example, [51], which directly affects the overall reconstruction efficiency, although it segments the point cloud planes well. Compared with the point-based segmentation algorithm, the supervoxel-based method has lower undersegmentation rate and higher completeness in the region with uneven density. Since the BPPM algorithm is used in this article to merge the original segmentation, the common oversegmentation problem in supervoxel-based planar segmentation algorithms is solved. Specifically, the progressive segmentation algorithm proposed in this article is able to handle fragmentation to optimize the reconstruction results due to the inclusion of supervoxel neighboring features in the coarse segmentation process and the termination condition of BPPM model fusion that all the primitives are smaller than the planar MSE threshold. The segmentation of small roofs like dormers can obtain better experimental results on AHN3 and DALES data, and it is necessary to resize the supervoxel resolution on the sparse density in Vaihingen data  and rely on the cost function optimization to get the building details, but the segmentation algorithm in this article is superior to the current segmentation algorithms.
The building roof reconstruction model obtained in this article has fewer planes than the 2.5-DDC method and the City3D  method, and the model is lighter, as shown in Table VI. Although quantitatively more than the 2019TINS method, the experimental results show that the algorithm in this article outperforms the 2019TINS method in terms of topological accuracy and regularity of the model, because this article considers the regularity of the boundary during the reconstruction process. Also, the reconstruction process in this article requires no input of the building footprint. c) Limitations: The algorithm proposed in this article is able to reconstruct both accurate and watertight building models on ALS point cloud datasets. However, there are still limitations of the algorithm that need to be addressed. Experimental results show that the algorithm proposed in this article is not automatically parameterized and requires the adjustment of parameters to obtain accurate building models on different density datasets. Furthermore, the sparse or missing roof point cloud data lead to incorrect model reconstruction. In addition, the point cloud reconstruction of buildings with large curvature is not considered in this article, and the above issues will be considered in future research to reconstruct more compact and regular building models.

V. CONCLUSION
In this article, a reconstruction algorithm was proposed and evaluated on three public ALS datasets. Three main contributions were included: BPPM-based planar primitive segmentation, CKL-based building footprint extraction, and energyoptimized roof reconstruction method. The proposed method was compared with the state-of-the-art reconstruction algorithms in the literature. Qualitative results showed that the model obtained in this article is the most regular and compact and effective. In general, this article achieved the following objectives.
1) The BPPM algorithm obtained accurate and complete planes.
2) The regularity of the building boundary was improved.
3) The model polygon corner points were simplified, and the topology and accuracy of the model were guaranteed.
However, there are still areas for improvement in this article. Robustness for low-density point cloud data needs to be enhanced, and surface roof reconstruction will be considered for research. In future research, more refined 3-D models will be reconstructed in combination with deep learning theory.