Compound 3D building modeling with structure-aware partition and primitive assembly from airborne laser scanning point clouds

ABSTRACT

Urban digital twins and realistic three-dimensional (3D) scenes have further enhanced the ever-increasing demand for building models. Airborne Laser Scanning (ALS) point clouds are the predominant data source for 3D building reconstruction owing to their highly detailed spatial features. However, information on the reconstruction of compound buildings (i.e. buildings with diverse plane structures or various roof structures) is lacking. Existing methods often produce building models that suffer from topological inconsistencies and/or geometric errors. To address these issues, this study proposes an automatic building reconstruction approach based on ALS point cloud that offers both modeling flexibility and geometric consistency. Specifically, a structure-aware partitioning algorithm that partitions a complex building roof into several regular roof patches is incorporated. Subsequently, a primitive model library comprising representative primitives suitable for building reconstruction is created. Suitable primitive model selection and reconstruction optimization are employed to ensure the accuracy and compactness in the final model. To comprehensively evaluate the performance, two typical ALS point cloud datasets are tested. Experiments demonstrate the generated models can achieve a boundary consistency of 92% while ensuring Root Mean Square Error below 0.09 m. Comparative experiments with state-of-the-art methods further validate the effectiveness of the proposed approach in compound building reconstruction.

KEYWORDS:

• Airborne laser scanning point cloud
• building modeling
• structure-aware partitioning algorithm
• primitive model library
• reconstruction optimization

1. Introduction

Urban digital twins and realistic three-dimensional (3D) scenes are integral components of the digital earth that have further enhanced the ever-increasing demand for 3D building structure reconstruction (Guo et al. 2021). Geographic Information System (GIS)-based 3D building models can represent the spatiotemporal information required for smart city applications (Guo et al. 2024). In particular, structural building models are essential in various multiple geospatial applications such as urban planning (Liu et al. 2021), intelligent transportation (Anand et al. 2019), disaster assessment (Vetrivel et al. 2018), and deformation monitoring (Awrangjeb, Gilani, and Siddiqui 2018). Over the past few decades, Airborne Laser Scanning (ALS) technology has been widely used to efficiently capture accurate and highly detailed spatial data of urban scenes. ALS point cloud-based building structure reconstruction has attracted considerable attention from the photogrammetry, remote sensing, and computer vision communities (Bauchet and Lafarge 2019; Qin, Hu, and Dai 2018). However, owing to the limitations of existing methods, achieving robust and accurate building structure reconstruction remains challenging.

ALS point clouds can describe the details of building roofs well; however, they fail to fully capture the facades of buildings (Huang et al. 2022; Wu et al. 2017). Therefore, previous studies have primarily focused on roof reconstruction. Nevertheless, several significant obstacles can hinder the task: (1) Noise or outliers, missed points caused by occlusions, and varying point densities can greatly impact modeling precision (Chen et al. 2022). (2) Numerous compound buildings exist in urban scenes, presenting challenges during crucial steps (e.g. segmentation, assembly, and optimization) in the modeling process (Li et al. 2022; Weixiao et al. 2023). (3) Urban buildings have various roof structures (e.g. flats, gables, hip roofs), requiring more flexible and generic primitive model libraries to achieve precise modeling (Jain and Susan 2013; Shao et al. 2021).

Extensive studies have attempted to address these challenges and enhance the effectiveness, robustness, and accuracy of ALS data-based roof reconstruction. Existing methods can be primarily categorized into data-, model-, and hybrid-driven. Data-driven methods extract geometric features directly from point clouds and assemble them into complete models by considering the topological relationships between them (Chen et al. 2020; Poullis 2013; Wang et al. 2020). Although these methods preserve the building details well, the results lack geometric consistency and are susceptible to noise and missing points. Model-driven methods create primitive model libraries that cover multiple building roofs. A primitive model can match point clouds by changing its size and/or shape (Dehbi, Koppers, and Plümer 2019; Li, Zhang, and Shan 2020; Zhang, Li, and Shan 2021). Model-driven methods are more robust and geometrically consistent; however, they may miss building details. The key to hybrid-driven methods lies in achieving the optimal reconstruction result by striking a balance between data and primitive models (Malihi et al. 2016; Li, Rottensteiner, and Heipke 2019). Hybrid-driven methods achieve more accurate models than those achieved using model-driven methods; however, for datasets with dozens of buildings, the processing time is longer, and the model volume is significantly larger. Considering the above, in this study, a model-driven method is adopted for the urban scene dataset.

This study focuses on a novel model-driven method for generating frame-structured rooftop models, particularly for urban scenes with a multitude of compound buildings. The main contributions of this study are summarized as follows: (1) A novel structure-aware partitioning method for handling compound buildings with diverse plane structures is introduced; regular partitioning of complex structures is implemented, and good results are achieved. (2) A primitive model library comprising representative primitives in the urban scene is established, balancing the diversity and effectiveness of reconstruction. (3) An energy function for selecting the optimal primitive model for each partitioned regular roof patch is introduced; the function measures the degree of matching between the roof patch and each primitive model transformed by the non-rigid Coherent Point Drift (CPD) algorithm.

The remainder of this paper is organized as follows: Section 2 presents an overview of recent studies on ALS point cloud-based building modeling and summarizes their advantages and disadvantages. Section 3 presents the reconstruction method, including a regular structure-aware partitioning of the building roof, suitable primitive model selection, and reconstruction optimization. Section 4 presents the experimental results and relevant analyses. Finally, Section 5 concludes the study.

2. Related works

ALS point cloud-based building modeling has been a highly researched topic over the past few decades. This section presents a perspective on the core techniques, including primitive model library creation and assembly, and reviews existing building modeling methods.

2.1. Creation of primitive model library

Primitive model libraries are a fundamental requirement for model-driven methods. The creation of a flexible primitive model library is crucial for achieving effective and accurate modeling. In general, more diverse primitive models enable more accurate results. Therefore, the development of primitive model libraries that accommodate complex building roofs in urban scenes is an important research topic in the domain of model-driven methods.

Initial studies on primitive model libraries focused only on common individual roof structures. Zhang and Fan (2022) integrated LiDAR point clouds with optical imagery to primitive models with clear edges and proposed simple primitive models with cubic or hip roof structures. Subsequently, a parametric approach was employed to standardize the geometric characteristics of the primitive models. Huang, Brenner, and Sester (2013) designed three primitive sets consisting of a limited number of planes, achieving building modeling with only a few parameters. However, this approach is applicable only in simple urban scenarios. Eventually, additional primitive models were added. For instance, Gooding, Crook, and Tomlin (2015) defined a primitive with a ridge length that was less than 35% of the footprint length of the hip primitive, whereas a length that was greater than 68% of the footprint length was considered as the long hip primitive. Zheng, Weng, and Zheng (2017) categorized roof types based on the orientation of the roof surface. A gable roof is defined as a roof with two upward-sloping sides that converge at the ridge. A hip roof is characterized by all sides sloping towards the ridge, albeit usually with gentler slopes compared to those of a gable roof. In recent years, primitive model libraries have been further expanded. The BuildingReconstruction software (Buyukdemircioglu, Kocaman, and Isikdag 2018) includes 32 CityGML-format roof types for automatic reconstruction. In addition, Zhang, Li, and Shan (2021) developed a secondary primitive library for model bonding to address the challenges encountered during roof reconstruction, such as chimneys and attics, thereby facilitating the handling of more complex scenarios.

These methods can perform well in scenes with individual buildings and simple horizontal structures. However, there is a lack of primitive model libraries that consider ubiquitous compound buildings that do not meet the accuracy requirements for reconstruction. Therefore, it is necessary to develop more general and detailed primitives as well as a flexible structure partitioning strategies to expand the applicability of modeling.

2.2. Primitive model assembly

Substantial research focuses also on the primitive model assembly, involving two key techniques: primitive recognition and geometric fitting. Primitive recognition employs traditional or deep learning-based methods to make the primitive model match the point cloud. Geometric fitting ensures that the shape, size, and position of the primitives are consistent with the point cloud.

For primitive recognition, deep learning-based methods have attracted the attention of researchers (Chen et al. 2021; Ponciano et al. 2021; Yang, Hou, and Li 2023; Zhang et al. 2019). Bittner et al. (2019) conducted a study on the classification of roofs, including flat and non-flat roofs, utilizing a 50-cm resolution satellite image within a conflict network. In contrast, Pohle-Fröhlich et al. (2019) evaluated the performances of PointNet, PointNet++, and U-Net in 3D roof segmentation based on slope orientation and found that PointNet++ and U-Net outperform other methods in terms of accuracy. To further improve the automation and efficiency of 3D roof segmentation, Zhang et al. (2022) proposed a two-branch network comprising an instance and a semantics branch. The extracted instance and semantics features are aggregated to obtain a more discriminative feature, thereby improving the accuracy of semantic predictions. Learning-based methods offer the advantage of a high recognition accuracy when trained on sufficient labeled samples. However, the labeling process is labor-intensive and time-consuming.

Considering the difficulties associated with deep learning-based methods, in traditional methods, Vo et al. (2015) applied a coarse-to-fine strategy to detect plane primitives based on the similar geometric properties of adjacent points and improved the efficiency by octree-based voxelization. The topological information between the geometric elements of the roof is also typically utilized to recognize the roof. Wang, Peethambaran, and Chen (2018) applied the Roof Topology Graph (RTG) to describe the topological information of building roofs. Li et al. (2022) further expanded the four types of sloping roofs, namely shed, gable, hip, and pyramid, and recognized them by selecting the primitives matching the most with the point cloud. The primitive model with the highest score was determined as the best match by introducing three score functions and measuring the matching results. Dey et al. (2023) promoted 3D building reconstruction by segmenting different feature points (such as boundaries, folded edges, and roof points) on the building roof and then used the obtained feature points to establish the relationship between feature lines and roof patches.

After primitive recognition, geometric fitting entails further determination of the parameters of each primitive, such as length, width, height, center position, and rotation angle. Zhang, Li, and Shan (2021) simplified the parameters of primitives into two distinct categories, namely shape and position parameters, and used PointNet++ for their determination. This method can also be extended to construct superstructures, such as attics and chimneys. Similarly, to optimize the processing time, non-learning studies focusing on parameter determination were also conducted. Dehbi et al. (2021) proposed an active sampling strategy that pairs data samples to form surfaces, enabling the extraction of essential parameters, including azimuth, inclination, ridge height, internal accuracy, and consistency. Li et al. (2022) employed an improved RANSAC multi-primitive segmentation framework for roof patch segmentation. This method involved an energy function to facilitate parameter determination, effectively addressing the challenges of sparse and partially missing point clouds. Tarsha Kurdi et al. (2023) achieved cylindrical modeling by determining the axis and cross-section. A new method for automatically generating cross-sections was introduced. In addition, potential points related to the building decoration were used to improve the building models. These methods only consider the geometric inconsistencies between primitives, resulting in poor detail accuracy in the final model.

In general, there are traditional geometry-based and deep learning-based methods for modeling 3D building roof models from ALS point cloud. The geometry-based methods can be further divided into model-driven (top-down) and data-driven (bottom-up) approaches. Model-driven methods attempt to approximate building structures using a combination of predefined simple primitives. The key challenges lie in constructing the basic roof shapes into appropriate primitives and combining them correctly based on geometric positioning (Huang, Brenner, and Sester 2013; Li and Shan 2022; Li, Zhang, and Shan 2020; Xiong et al. 2015). Data-driven methods, on the other hand, first extract geometric features (points, lines, planes) from the point cloud using segmentation, and then analyze the topologies of these features to reconstruct the roof model (Li et al. 2017; Nan and Wonka 2017; Ochmann, Vock, and Klein 2019; Valero, Bosché, and Forster 2018). For deep learning-based methods, the neural networks are employed to assist various aspects of 3D building reconstruction, such as detecting vertices and edges, predicting roof types, and constructing topological graphs. However, these approaches are often limited by the availability of labeled 3D building model data (Jiang et al. 2023; Li et al. 2022; Yu, Liu, and Zhang 2022; Zhang, Li, and Shan 2021).

3. Method

The framework is shown in Figure 1 and comprises three main components: (1) regular structure-aware partitioning, (2) primitive model assembly, and (3) reconstruction optimization. Specifically, the input data needs to be preprocessed to remove ground points and retain roof. This paper performs clustering on roofs and classifies them as regular or irregular roof patches using the Minimum Bounding Rectangle (MBR) method. For irregular roof patches, a structure-aware partitioning algorithm is employed to further partition them into regular roof patches. Subsequently, a primitive model library that encompasses common building roof types in the urban scene is created. Primitive models are selected and parameterized to best match the regular roof patches. Finally, reconstruction optimization is performed on the primitives to obtain a watertight model by eliminating geometric errors and addressing topological inconsistencies.

Figure 1. Framework of the proposed method.

3.1. Regular structure-aware partitioning of building roofs

A regular structure-aware partitioning algorithm is proposed to adaptively partition a building roof into several regular roof patches. The partitioning algorithm consists of two parts: irregular roof patch recognition (Section 3.1.1) and irregular roof patch partitioning (Section 3.1.2).

3.1.1. Irregular roof patch recognition using MBR

A roof patch is a contiguous region in a point cloud. In the partitioning algorithm, roof patches are utilized as fundamental units. To extract roof patches, the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm (Kumar and Reddy 2016; Zhang et al. 2022) is applied because of its ability to aggregate proximate and densely populated points into groups. This algorithm is based on density clustering and the clustering of areas with sufficient space density into a point set. For compound buildings, in order to ensure the extraction accuracy of the roof patches, we process the rough and fine extraction through multiple clustering.

Both regular and irregular roof patches are observed in compound buildings after clustering. Regular roof patches are defined as roof patches with rectangular horizontal contours and no hollow regions inside, whereas irregular roof patches are not defined as such. To improve the modeling accuracy, an irregular roof patch recognition method was designed based on the MBR (Kwak and Habib 2014). Specifically, the points of each roof patch are divided into equal-sized grids along a rectangle constructed using the minimum and maximum coordinates of the point cloud, as depicted in Figure 2. The grid resolution has been determined through multiple experiments and is set to 1 m in this study. The accumulated area of grids containing falling points is calculated as S_grid, and the rectangular area is S_MBR, as shown in Eq. 1(1) $\begin{array}{rl}& S_{\mathrm{grid}}=\sum _{i=1}^N(\mathrm{gri}{d}_i\cdot {\mathrm{\Delta }}^2)\\ & \mathrm{gri}{d}_i=\{\begin{array}{c}1,\mathrm{if}\mathrm{card}\{p/p\in \mathrm{gri}{d}_i\}\ne 0\\ 0,\mathrm{otherwise}\end{array}\\ & S_{\mathrm{MBR}}=l\cdot w\end{array}$(1) . A grid with no falling points is not counted. (1) $\begin{array}{rl}& S_{\mathrm{grid}}=\sum _{i=1}^N(\mathrm{gri}{d}_i\cdot {\mathrm{\Delta }}^2)\\ & \mathrm{gri}{d}_i=\{\begin{array}{c}1,\mathrm{if}\mathrm{card}\{p/p\in \mathrm{gri}{d}_i\}\ne 0\\ 0,\mathrm{otherwise}\end{array}\\ & S_{\mathrm{MBR}}=l\cdot w\end{array}$(1) where N denotes the number of all grids, Δ denotes the grid resolution, p denotes a point in grid, card denotes the cardinality of grids, and l and w denote the length and width of the MBR of a roof patch, respectively.

Figure 2. Irregular roof patch recognition based on the Minimum Bounding Rectangle (MBR) method (top view).

Considering the geometric characteristics of roof patches, the area of their horizontal contours is always less than or equal to the MBR area. The MBR ensures that all points are included while maintaining the smallest possible area, with the area S_MBR being usually larger than the grid area S_grid. In particular, regular roof patches exhibit small area differences, whereas irregular roof patches often have significant area variations and can thus be readily distinguished. Based on this, this study calculates the discriminant factor IFR indicating whether a roof patch is regular or not, as shown in Eq. 2(2) $\mathrm{IFR}=\{\begin{array}{c}\mathrm{yes},\mathrm{if}|S_{\mathrm{MBR}}-S_{\mathrm{grid}}|\ge S_{\mathrm{MBR}}\cdot \gamma \\ \mathrm{no},\mathrm{otherwise}\end{array}$(2) . (2) $\mathrm{IFR}=\{\begin{array}{c}\mathrm{yes},\mathrm{if}|S_{\mathrm{MBR}}-S_{\mathrm{grid}}|\ge S_{\mathrm{MBR}}\cdot \gamma \\ \mathrm{no},\mathrm{otherwise}\end{array}$(2) where γ denotes the partition threshold (γ = 0.3 is suggested; the discussion on the specific value is presented in Section 4.4). If the area difference between S_MBR and S_grid is larger than S_MBR • γ, the horizontal contour of the roof patch differs significantly from a rectangle, i.e. the roof patch is irregular.

3.1.2. Irregular roof patch partitioning

To facilitate subsequent primitive-based modeling and ensure the modeling accuracy of compound buildings, it is necessary to further partition each irregular roof patch into regular roof patches that are easy to reconstruct. A regular structure-aware partitioning method is proposed that partitions an irregular roof patch into multiple regular roof patches by analyzing the corner points on the roof contour and ridge (only for sloping roofs). The main steps are as follows:

• Step 1: For each irregular roof patch, the alpha-shape method of Li, Rottensteiner, and Heipke (2019) is used to extract the outer contour, the value of alpha is determined as 1 m based on the average point density of the datasets selected in this paper and the region-growing algorithm of Wang et al. (2020) is used to detect the planes.

• Step 2: From the detected planes, the plane closest to the highest point is identified, and the included angle between the normal and vertical vectors is calculated. If the included angle is less than a certain threshold (i.e. 10°), the roof patch is treated as a flat roof; otherwise, it is treated as a sloping roof.

• Step 3: For each flat roof, if the outer contour is a concave polygon, the building typically exhibits either an L-shaped (as shown in Figure 3a) structure roof. The RANSAC algorithm is employed to detect each outer contour line and select one passing through the concave point as the partition line. For an outer contour with multiple concave points like C-shaped of E-shaped, the line passing through two concave points is prioritized as the partition line. If the outer contour is a convex polygon, the buildings are usually ring-shaped (as shown in Figure 3b). The inner contour lines are extracted as partition lines.

• Step 4: For each sloping roof, if the outer contour is a concave polygon (L-shaped, as shown in Figure 3c), the corner points of the roof contour and ridge are extracted by intersecting the planes. The concave point and its nearest contour and ridge corners are selected to form the partition plane. In the case of a convex polygon (as shown in Figure 3d), the corner points are divided into ridge corners and inner and outer contour corners. A ridge corner, its nearest inner contour corner, and the nearest outer contour corner are selected to form a partition plane.

• Step 5: The irregular roof patches are partitioned into several regular roof patches based on the formed partition planes. It is possible to partition multiple times until all results meet the criteria for determining regular roof patches.

Figure 3. L – and ring-shaped structure roofs (a. L-shaped flat roof; b. ring-shaped flat roof; c. L-shaped sloping roof; d. ring-shaped sloping roof).

The partitioning procedure is shown in Figure 4, taking a C-shaped sloping roof as an example.

Figure 4. Partitioning procedure.

3.2. Primitive model assembly for reconstruction

Considering the common characteristics of urban building structures, a novel primitive model library encompassing nine types of primitive models is designed. This study attempts to use as few entries as possible to build them, thereby improving the adaptability and flexibility of the library. Primitives have rectangular footprints that align with the common architectural design of most buildings.

3.2.1. Primitive model library creation

Owing to the complexity and diversity of roof types, it is essential to create a primitive model library. Complex primitives can enhance the modeling accuracy; however, they also introduce additional parameters, which may reduce the generality. Therefore, an ideal primitive model library should strike a balance between diversity and the use of as few parameters as possible.

To achieve this, nine primitive models are designed in the library. These primitives cover the common roof forms of buildings in the urban scene. The nine primitive models and their control parameters are shown in Figure 5. For each primitive model, the control parameters define the spatial attributes. These parameters include shape and position parameters. The former describe the shape and size of the primitive model, including the wall length l_w, wall width w_w, wall height h_w, roof height h_r, roof length l_r and roof width w_r; the latter describe the precise orientation and position of the primitive model, which are defined by the center point position o and the rotation angle k around the Z-axis.

Figure 5. Primitive model library and control parameters.

3.2.2. Suitable primitive model selection

Multiple regular roof patches are obtained after partitioning. A suitable primitive model from the primitive model library is selected to best match each regular roof patch. However, the transformation between the roof patch and the primitive model is non-rigid. In contrast, the non-rigid CPD algorithm (Myronenko and Song 2010) determines a transformation between the reference point cloud and the possible deformed point cloud, guaranteeing optimal registration between them. The algorithm exhibits robust and accurate performance in processing point clouds with noise, point density variation, and conclusions. Compared with traditional registration methods, it shows robust and accurate performance in processing point clouds with noise, point density variation and conclusions. Therefore, this study applies the non-rigid CPD algorithm to achieve the best match between any primitive model and the roof patch.

For each primitive model, the sampled point cloud is marked as T, and the point cloud of the roof patch is marked as Q. The transformation is shown in Eq.3: (3) $P(x_p,y_p,z_p)=T(x_t,y_t,z_t)+\mathrm{GW}$(3) where P is the transformed point cloud of T, G is a Gaussian distribution matrix of T, W is a transformation coefficient matrix, and GW is the displacement matrix. The Expectation Maximization (EM) optimization is used to calculate the GW. Specifically, the E-step is used to construct the probability distribution of p_mn. The M-step updates the displacement, as shown in Eq. 4(4) $\begin{array}{rl}& p_{\mathrm{mn}}={\displaystyle \frac{p(m)}{{\sum }_{k=1}^{M}p(k)+{\displaystyle \frac{w}{1-w}{\displaystyle \frac{{(2\pi {\sigma }^2)}^{3/2}M}{N}}}}}\\ & p(i)=\exp (-{\displaystyle \frac{1}{2{\sigma }^2}\Vert Q_j-(T_i+G_{i}W){\Vert }^2})\end{array}$(4) : (4) $\begin{array}{rl}& p_{\mathrm{mn}}={\displaystyle \frac{p(m)}{{\sum }_{k=1}^{M}p(k)+{\displaystyle \frac{w}{1-w}{\displaystyle \frac{{(2\pi {\sigma }^2)}^{3/2}M}{N}}}}}\\ & p(i)=\exp (-{\displaystyle \frac{1}{2{\sigma }^2}\Vert Q_j-(T_i+G_{i}W){\Vert }^2})\end{array}$(4) where M and N represent the total numbers of points in T and Q, respectively, T_i and Q_j represent points in T and Q, respectively, σ² is the covariance, w is the amount of noise (0 < w < 1), and G_i is an exclusive matrix for T_i. After transforming T, we obtain P, which aligns closely with Q. It is noteworthy that P is not identical to Q because of the differences in the number of points and positions of adjacent points that do not completely overlap.

An energy function is then devised to compute the degree of matching between P and Q. Each primitive model in the primitive model library is evaluated, and the one with the highest score is designated as the optimal primitive model, as shown in Eq. 5(5) $E=\mathrm{EOverlapping}+\mathrm{\lambda EVarying}$(5) : (5) $E=\mathrm{EOverlapping}+\mathrm{\lambda EVarying}$(5) where E_Overlapping is the overlapping term, E_Varying is the varying term, and λ is the weight coefficient for balancing them. The computational result of the energy function can be regarded as a measure of the degree of matching. A higher value of E indicates a better match between the primitive model and roof patch. The overlapping term is used to quantify the degree of overlap between P and Q and is calculated using Eq. 6(6) $\mathrm{EOverlapping}={\displaystyle \frac{1}{N}\sum _{i=1}^N\exp (-||\mathrm{ai}-\mathrm{bj}|{|}^2)}$(6) : (6) $\mathrm{EOverlapping}={\displaystyle \frac{1}{N}\sum _{i=1}^N\exp (-||\mathrm{ai}-\mathrm{bj}|{|}^2)}$(6) where a_i is a point in Q, and b_j is its nearest point in P. E_Overlapping amplifies the difference between a_i and b_j through the distance: When the two points are closer, point a_i contributes a higher score to E_Overlapping. This process is repeated for each primitive model in the library, resulting in an E_Overlapping value for that specific primitive model with respect to the original point cloud. The varying term primarily focuses on the roof parameters that influence the roof structure within the control parameters. These parameters include roof height h_r, roof length l_r, and roof width w_r. The variations in these three parameters directly affect the roof type, leading to changes in the primitive model type. E_Varying is calculated using Eq. 7(7) $\begin{array}{c}\mathrm{EVarying}=\mathrm{\alpha D}(h_r)+\mathrm{\beta D}(l_r)+D(w_r)\\ D(h_r)=\exp (-|h_r^Q-h_r^P|)\\ D(l_r)=\exp (-|l_r^Q-l_r^P|)\\ D(w_r)=\exp (-|w_r^Q-w_r^P|)\end{array}$(7) : (7) $\begin{array}{c}\mathrm{EVarying}=\mathrm{\alpha D}(h_r)+\mathrm{\beta D}(l_r)+D(w_r)\\ D(h_r)=\exp (-|h_r^Q-h_r^P|)\\ D(l_r)=\exp (-|l_r^Q-l_r^P|)\\ D(w_r)=\exp (-|w_r^Q-w_r^P|)\end{array}$(7) where α and β are weight coefficients. D(•) denotes the parameter difference calculated using P and Q. A smaller difference between the roof parameters results in a greater contribution to E_Varying.

In this algorithm, each primitive model in the library is iterated, and its degree of matching with the original point cloud is computed using the energy function. Finally, the primitive model with the highest matching degree is selected as the suitable primitive model. This procedure is illustrated in Figure 6.

Figure 6. Procedure of the suitable primitive model selection algorithm.

3.3. Reconstruction optimization

Once multiple primitive models are obtained, it is necessary to eliminate geometric errors and address topological inconsistencies. To quantify the geometric discrepancies between adjacent primitives, an improved Hausdorff distance (Taha and Hanbury 2015) is employed according to Eq. 8(8) $\begin{array}{c}\mathrm{HD}(A,B)=min\{\mathrm{hd}(A,B),\mathrm{hd}(B,A)\}\\ \mathrm{hd}(A,B)=\mathrm{avg}\{d(p_i,p_b)|p_i\in A,p_b\in B\},p_b\mathrm{is}\mathrm{the}\mathrm{nearest}\mathrm{point}\mathrm{of}{p}_i\\ \mathrm{hd}(B,A)=\mathrm{avg}\{d(p_j,p_a)|p_j\in B,p_a\in A\},p_a\mathrm{is}\mathrm{the}\mathrm{nearest}\mathrm{point}\mathrm{of}{p}_j\end{array}$(8) : (8) $\begin{array}{c}\mathrm{HD}(A,B)=min\{\mathrm{hd}(A,B),\mathrm{hd}(B,A)\}\\ \mathrm{hd}(A,B)=\mathrm{avg}\{d(p_i,p_b)|p_i\in A,p_b\in B\},p_b\mathrm{is}\mathrm{the}\mathrm{nearest}\mathrm{point}\mathrm{of}{p}_i\\ \mathrm{hd}(B,A)=\mathrm{avg}\{d(p_j,p_a)|p_j\in B,p_a\in A\},p_a\mathrm{is}\mathrm{the}\mathrm{nearest}\mathrm{point}\mathrm{of}{p}_j\end{array}$(8) where A and B represent the point clouds sampled from the surfaces adjacent to the two primitives, hd(A, B) and hd(B, A) denote the one-way Hausdorff distances from A to B and B to A, respectively, and d(,) denotes the Euclidean distance. The distance between the two point clouds is calculated by averaging the point-to-point distances. A smaller value of HD(,) indicates smaller geometric discrepancies between the adjacent primitives.

In the case of compound buildings, the potential geometric errors between primitive models can be categorized into four types: ① rotation angle error of adjacent primitives (Figure 7a); ② crack error (Figure 7b); ③ intersection error (Figure 7c); ④ alignment error of adjacent primitives (Figure 7d).

Figure 7. All possible geometric errors (a. rotation angle error; b. crack error; c. intersection error; d. alignment error).

Processing is not only applicable to individual errors but also to complex errors that involve multiple errors within adjacent primitives. The specific steps are as follows:

• Step 1: Rotation angle error adjustment. Since adjacent primitives have been approximately assembled, a slight angle adjustment is performed within a specified range (e.g. −5° to +5°) to achieve rotational consistency within the model. The rotation center is set as the center of the source primitive, and the target primitive is incrementally rotated by a small angle (e.g. 0.1°). The position yielding the smallest discrepancies is determined as the final position of the target primitive.

• Step 2: Crack or intersection error elimination. The average distance between the adjacent surfaces of the target and source primitives is calculated as the correction width. Using the direction and correction width, a displacement vector is determined and applied to the target primitive to resolve cracks or intersection errors.

• Step 3: Alignment error correction. The two nearest corner points from the target and source primitives, namely P_target and P_source, respectively, are selected. A horizontal displacement vector ${\mathrm{V}}_{\mathrm{horizontal}}$ is calculated from P_source to P_target, as shown in Eq.9. The target primitive is then adjusted by applying the horizontal displacement vector to align the adjacent primitives. (9) $\begin{array}{rl}& {\mathrm{V}}_{\mathrm{horizontal}}={\mathrm{V}}_{\mathrm{ST}}-[{\mathrm{V}}_{\mathrm{ST}}\cdot (0,0,1)]\cdot (0,0,1)\\ & {\mathrm{V}}_{\mathrm{ST}}=P_{\mathrm{target}}-P_{\mathrm{source}}\end{array}$(9)

Following the pairwise optimization, the next two primitives also undergo pairwise optimization, and this process continues until all primitives have been considered. This results in a compact model in the horizontal direction. Finally, the wall and roof heights of the primitive models partitioned from the same roof patch are unified, eliminating geometric errors in the vertical direction of the model.

4. Experiments

4.1. Dataset description

To test the proposed method, two ALS point cloud datasets are selected as experimental datasets. Dataset Ⅰ comprises over 100 individual buildings within the campus of Nanjing University of Information Science & Technology (NUIST). Dataset Ⅱ is AHN3 – an openly available ALS point cloud dataset that covers the entire Netherlands. AHN3 includes numerous buildings with diverse plane structures and is therefore suitable for verifying the performance of modeled buildings with complex horizontal layouts. Figure 8 provides an overview of the raw point clouds in the two datasets and Table 1 lists a detailed description of the dataset parameters.

Figure 8. Overview of raw point clouds.

Table 1. Detailed descriptions of dataset parameters.

Parameters	Dataset Ⅰ	Dataset Ⅱ
Name	NUIST campus	AHN3
Location	Nanjing, China	Delft, Netherland
Type	ALS point cloud	ALS point cloud
Number of Points	24.3 M	20 M
Average Point Span/m	0.15	0.20
Collection Time	2022/10	2021/8
Area Size/m^2	1,904 × 1,428	1,600 × 1,300
Characteristics	Buildings with various roof structures; dataset exhibits obvious terrain undulation.	Numerous buildings with diverse plane structures.

4.2. Reconstruction result evaluation

4.2.1. Visual results

To visually present the global distribution of the building models and showcase the reconstruction details, the reconstruction results (depicted as orange models) are displayed alongside the background images in Google Earth in Figure 9. Additionally, three magnified views are included below each main scene to highlight specific details.

Figure 9. Three-dimensional (3D) views of the reconstruction models and their corresponding background images.

Figure 9 shows the reconstruction results and distribution of all buildings in the two datasets. In the magnified views of Dataset Ⅰ, various roof types, including cube, hip, and mansard, have been successfully reconstructed. Additionally, the method demonstrates its ability to produce favorable results for buildings with complex horizontal layouts, such as L-shaped, H-shaped, and even more intricate roof structures, as shown in the magnified views of Dataset Ⅱ.

4.2.2. Procedure results from modeling typical buildings

A selection of typical buildings from the two datasets is presented in Figure 10 to demonstrate the procedure results of the proposed method. These buildings represent various compound building types and exhibit distinct characteristics, rendering them suitable for evaluating the performance of the proposed method. Specifically, the buildings in the second, fourth, and fifth rows demonstrate complex horizontal layouts, including L-shaped, C-shaped, ring-shaped, and other configurations. The buildings in the first and third rows show intricate roof structures, such as stepped roofs and non-standard shapes. Furthermore, the structures in the last two rows are intricate compound buildings characterized by both complex horizontal layouts and roof structures.

Figure 10. Procedure results (a. google earth images; b. building point clouds; c. extraction results of roof patches; d. partitioning results in red box; e. primitive model assembly in green box; f. final models).

Figure 10(c) illustrates the successful extraction of multiple roof patches with varying shapes or heights, including large continuous surfaces such as ceilings and intricate structures such as attics (e.g. first row in Figure 10). This indicates that the method efficiently preserves not only the primary building structure but also intricate details. Additionally, the method partitions irregular roof patches by considering the roof structures, as shown in the red box in Figure 10(d). The robustness of the primitive selection algorithm is demonstrated even for roof patches with small areas.

The optimized models exhibit closely joined primitives, and there are no noticeable geometric errors in the final models. This is particularly evident in cases where multiple primitives belong to a single irregular roof patch, as depicted in the green box in Figure 10(d) and (e). The absence of jagged incisions resulting from geometric errors indicates that the method can withstand the impact of intermediate operations and generate watertight models.

4.2.3. Reconstruction accuracy assessment

A set of representative buildings was selected from the two datasets to highlight the details, as shown in Figure 11. The first row shows the reconstruction models, the second row shows the corresponding aerial images, and the third row shows the models overlaid with the building point cloud. The accuracy of the reconstruction is demonstrated by combining the point clouds with the reconstruction models. The colors in the visualization are assigned based on the distances between the points and their nearest surfaces in the model. Blue indicates a small modeling error, whereas red indicates a large modeling error.

Figure 11. Set of models for accuracy assessment (a-c. Dataset I; d-f. Dataset II).

The proposed method successfully reconstructs roof structures with different patterns and eliminates geometric errors between primitives. The final models presented in Figure 11 are generally compact. Importantly, all primitives exhibit large blue areas, indicating a high reconstruction accuracy for multiple primitives. Notably, in the third row, the area occupied by blue is significantly larger than those occupied by the other colors. Besides efficiently reconstructing complex structures and various roof types, the method also captures and displays additional superstructures, which are shown in blue.

The quantitative results of the buildings shown in Figure 11 are provided in Table 2. Regression analysis is conducted between the Root Mean Square Error (RMSE) (Huang et al. 2022), Maximum Error (ME), and Intersection over Union (IoU) (Zhang et al. 2022) and the processing time. The diagrams are shown in Figure 12.

Figure 12. Distribution of the RMSE, ME and IoU (a. RMSE/ME of buildings in two datasets, b. IoU of buildings in two datasets).

Table 2. Accuracy evaluation of reconstruction models. (Building NO. corresponds to Figure 11).

Dataset	Building NO.	Points	RMSE/m	ME/m	IoU	Time/s
Dataset Ⅰ	1	28,134	0.0772	0.1049	94%	20
	2	36,742	0.0709	0.1158	96%	18
	3	74,405	0.0811	0.1181	94%	65
	4	44,092	0.0634	0.0912	93%	42
	5	96,562	0.0765	0.1457	94%	12
	6	70,280	0.0806	0.1461	92%	23
	7	82,305	0.0782	0.1453	94%	50
	8	78,095	0.0583	0.1441	93%	53
	9	98,544	0.0651	0.1080	93%	74
Dataset Ⅱ	10	94923	0.0539	0.1053	96%	56
	11	60024	0.0491	0.0905	93%	60
	12	68044	0.0513	0.0825	95%	43
	13	62976	0.0457	0.0972	94%	49

Table 2 and Figure 12 show that the errors of buildings with hip structures (Buildings 5-6) and mansard roofs (Buildings 7-8) are slightly higher than those of buildings comprising only cube primitives (Buildings 1-2). This is primarily owing to the increased computational complexity involved in the reconstruction of such complex roof structures. Nevertheless, even for buildings with multiple different primitives, the RMSE remains below 0.09 m, and IoU exceeds 92%, as indicated in Table 2. This demonstrates that the proposed method maintains a high level of reconstruction accuracy even in the presence of diverse primitive compositions.

Furthermore, for buildings with diverse plane structures (Buildings 10-13), Table 2 shows that the RMSE and ME values of the reconstruction models are relatively low. This indicates that the method is effective in handling buildings with diverse plane structures while still achieving accurate reconstructions. Nevertheless, buildings with diverse plane structures may result in a more complex footprint, leading to an IoU decrease.

4.3. Performance assessment and hyper-parameter discussion

A compound building point cloud is sampled to generate a set of point clouds with varying densities for performance evaluation. A parameter, namely Subdivision Level (SL), is introduced to represent the degree of subdivision. A higher SL corresponds to reduced point density. The modeling results for different point densities are presented in Figure 13a, and a summary of the results is provided in Table 3.

Figure 13. Performance assessment (a. varying point density data reconstruction; b. Gaussian noise data reconstruction).

Table 3. Reconstruction performance of data with varying point density.

Property	Original data	SL = 15	SL = 20	SL = 25	SL = 30
Points	74,405	8,275	2,446	1,324	431
Point span/m	0.1276	0.5371	0.7278	1.0412	1.6492
RMSE/m	0.0811	0.1061	0.1125	0.2053	0.4731

Figure 13a demonstrates that when the point span is less than 1 m, even after significant down-sampling, the reconstruction models can retain all essential features of the sparse data samples, such as stepped roofs and attics. However, when the point span exceeds 1 m, a portion of the points exceeds the minimum clustering density set by the DBSCAN algorithm. Consequently, the clustering algorithm may fail to extract the correct roof patches, leading to a loss of certain features in the reconstruction results.

Gaussian noise is introduced into the original point cloud to assess the robustness of the proposed method in relation to noise. In Figure 13b, the noise level σ represents the point density, with a higher noise level corresponding to sparser points. When σ is below 1 m, the reconstruction models can preserve all the important features, including stepped roofs and superstructures. However, when the noise level σ exceeds 1 m, an excessive number of noise points fill the gaps in the stepped roofs, causing the clustering algorithm to treat them as a single cluster. This can result in loss of detail and reduced accuracy in the reconstruction.

Table 4 provides an overview of the performance of the proposed method for different noise levels in the input data. Notably, increasing the noise level does not alter the total number of point clouds, but only affects the point span. When the noise level exceeds the minimum distance between two actual surfaces in a building, detecting the appropriate roof planes becomes challenging, leading to incorrect model generation. However, when the noise level σ is below 1 m, the model’s RMSE remains below 0.12 m. This demonstrates the robustness of the reconstruction process regarding handling different noise levels in the point clouds.

Table 4. Reconstruction performance of data with different noise levels.

Property	Original data	σ = 0.5m	σ = 0.7m	σ = 1.0m	σ = 1.5m
Points	74405	74405	74405	74405	74405
Point span /m	0.1276	0.2022	0.2429	0.3137	0.4437
RMSE/m	0.0811	0.1052	0.1176	0.2132	0.3518

To evaluate the impact of the partition threshold on the reconstruction, a C-shaped building is selected, and the partition threshold is varied multiple times while keeping the other steps of the reconstruction process unchanged. This allows the generation of a set of reconstruction models for comparison.

This study includes a discussion regarding the selection of the partition threshold in Eq. 2(2) $\mathrm{IFR}=\{\begin{array}{c}\mathrm{yes},\mathrm{if}|S_{\mathrm{MBR}}-S_{\mathrm{grid}}|\ge S_{\mathrm{MBR}}\cdot \gamma \\ \mathrm{no},\mathrm{otherwise}\end{array}$(2) , and a value of γ = 0.3 is chosen. This value is determined through multiple tests on finding the most suitable threshold. As shown in Figure 14, when the partition threshold is set below 0.3, redundant partitions are generated. As highlighted by the red box, the mansards are not fully integrated within the main C-shaped building. This can be attributed to redundant partitions that cause inaccuracies during determining the dimensions of the main structure. These errors affect the overall reconstruction of the model structure. Conversely, if the partition threshold is set too high, incomplete partitioning occurs. The Table 5 shows the reconstruction accuracy under different threshold conditions. Therefore, γ = 0.3 is a satisfactory partition threshold that strikes a balance between redundancy and completeness in the partitioning process.

Figure 14. Partition threshold comparison.

Table 5. Reconstruction performance of data with different partition threshold.

Property	γ = 0.2	γ = 0.3	γ = 0.4	γ = 0.5
RMSE/m	0.1107	0.0782	0.1964	0.1997
ME/m	0.1771	0.1453	0.2279	0.2306
IoU	92%	94%	78%	78%

4.4. Comparison experiment

To provide a comprehensive evaluation, the proposed method is compared with several state-of-the-art building reconstruction methods, namely the Polyfit (Nan and Wonka 2017), Optimal Model Fitting (OMF) (Zhang, Li, and Shan 2021), and RANSAC-MPR (Li et al. 2022) methods. This experiment assesses the strengths and weaknesses of the proposed method in relation to those of these existing methods. The reconstruction models generated by each method are visually represented by combining colored points with the corresponding reconstruction models, as depicted in Figure 15. This visual representation allows a qualitative assessment of the reconstruction accuracy achieved using each method (Table 6).

Figure 15. Comparative evaluation of compound building reconstruction results (a. building with stepped roof, b. C-shaped building, c. ring-shaped building).

Table 6. Performances of different methods.

NO.	Property	Polyfit	OMF	RAN-MPR	Proposed
Building (a)	Points	74405	74405	74405	74405
	RMSE/m	0.3443	0.2852	0.2063	0.0811
	ME/m	0.3203	0.2397	0.1846	0.1181
	IoU	87%	83%	91%	94%

The second row in Figure 15 shows that the Polyfit method can generate consistent surfaces on a large scale. However, it does not strictly adhere to geometric relationships, as evidenced by instances in which the ceiling is not parallel to the horizontal plane. Table 6 further supports this observation, as the Polyfit method exhibits the highest RMSE, indicating relatively large errors in the reconstruction. The OMF method employs deep learning for roof-type recognition; however, it is applicable only to individual buildings. As depicted in the third row of Figure 15, OMF underperforms in case of compound buildings, such as C-shaped or ring-shaped structures, mistakenly filling empty regions and leading to significant errors in the reconstruction process. The RANSAC-MPR method can handle complex structures, and Table 6 shows that it achieves the lowest RMSE among the three methods; however, it lacks a reconstruction optimization step. As depicted in the fourth row of Figure 15, there are noticeable misalignments and cracks in the reconstruction models. Notably, the RANSAC-MPR method employs a primitive model library with only four primitives, resulting in certain deficiencies in the reconstruction accuracy for some buildings. For instance, in Figure 15b, the mansard roof is incorrectly reconstructed as a cube, and the corresponding point cloud is highlighted in red, indicating a modeling error.

5. Conclusion

This study presents a fully automatic approach for modeling compound 3D buildings using ALS point clouds. The proposed method reconstructs compound building point clouds into multiple primitives through regular structure-aware partitioning and suitable primitive model selection. The final model is obtained through reconstruction optimization. Extensive experiments conducted on various datasets demonstrate that the reconstructed model achieves RMSE below 0.09 m, ME below 0.15 m, and IoU exceeding 92%. This method exhibits strong robustness in relation to both sparse point clouds and noise, ensuring the accurate appearance and reconstruction of both individual and compound buildings. Overall, this method achieves a fast and effective compound building reconstruction using ALS point clouds.

Importantly, the current primitive model library includes only roof primitives, which are sufficient for reconstructing the majority of buildings in urban areas. However, the method cannot effectively handle curved surfaces. In the future, expanding the primitive model library to include other primitives, such as spheres and cylinders, could further enhance its versatility and capabilities.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Additional information

Funding

This work was jointly supported by the National Natural Science Foundation of China [grant Nos. 42171433, 42101449, 42071440], the Postgraduate Research & Practice Innovation Program of Jiangsu Province [grant number KYCX22_1213], the Natural Science Foundation of Hubei Province, China [grant number 2022CFB773], the National Key Research and Development Program of China [grant number 2023YFB3906102], the Science and Technology Program of Southwest China Research Institute of Electronic Equipment [grant number JS20200500114] and the Chutian Scholar Program of Hubei Province.