AUTOMATIC 3D RECONSTRUCTION OF COMPLEX BUILDINGS FROM INCOMPLETE POINT CLOUDS WITH TOPOLOGICAL-RELATION CONSTRAINTS

: Automatic 3D building reconstruction from laser scanning or photogrammetric point clouds has gained increasing attention in the past two decades. Although many efforts have been made, the complexity of buildings and incompletion of point clouds, i.e., data missing, still make it a challenging task for automatic 3D reconstruction of buildings in large-scale urban scenes with various architectural styles. This paper presents an innovative approach for automatic generation of 3D models of complex buildings from even incomplete point clouds. The approach first decomposes the 3D space into multiple space units, including 3D polyhedral cells, facets and edges, where the facets and edges are also encoded with topological-relation constraints. Then, the units and constraints are used together to approximate the buildings. On one hand, by extracting facets from 3D cells and further extracting edges from facets, this approach simplifies complicated topological computations. On the other hand, because this approach models buildings on the basis of polyhedral cells, it can guarantee that the models are manifold and watertight and avoid correcting topological errors. A challenging dataset containing 105 buildings acquired in Central, Hong Kong, was used to evaluate the performance of the proposed approach. The results were compared with two previous methods and the comparisons suggested that the proposed approach outperforms other methods in terms of robustness, regularity, and accuracy of the models, with an average root-mean-square error of less than 0.9 m. The proposed approach is of significance for automatic 3D modelling of buildings for urban applications.


INTRODUCTION
Three-dimensional (3D) reconstruction from point clouds has been an active topic in photogrammetry and computer vision communities.However, missing data on the vertical surfaces, especially in the aerial laser scanning (ALS) data, could be a serious problem.Therefore, most of the previous methods, including model-driven, data-driven and hybrid-driven methods (Haala and Kada, 2010;Wang et al., 2018), only focused on reconstructing rooftops of buildings and produced 2.5D building models (Zhou and Neumann, 2010).Although dense image matching (DIM) point clouds generated from oblique images through multi-view stereo pipelines (Vu et al., 2012;Wu et al., 2018) provide more information on vertical surfaces than ALS data, the data missing issue still exists where there are occlusions between high-rise buildings that are closely located.
In addition, the architectural styles of buildings can vary greatly with culture, location and time (Qiao et al., 2009), which makes it impossible to fit the point clouds with a predefined model library as in model-driven methods (Huang et al., 2013;Kada and McKinley, 2009;Poullis and You, 2009).The complexity of buildings also increases the difficulties in topological computation in model-driven and hybrid-driven methods (Vosselman and Dijkman, 2001;Sohn et al., 2008;Verma et al., 2006), resulting in crack effects in the final models (Poullis, 2013;Xie et al., 2018) or needing extra work to correct the topological errors in the models manually (Xiong et al., 2014).And even when the topological errors are enforced to be correct, the results might violate some geometric relations with respect to regularity, such as co-planarity, parallelism and orthogonality.Particular methods have been developed to preserve the topological relations while reconstructing building rooftops (Chen et al., 2014;Chen et al., 2017), but they cannot capture the topological relations between two roof components when there is a large height jump.
A recent trend of generating building models is to decompose the 3D space into a set of basic units, e.g., boxes, facets and cells, and then use these basic units to fit the building surfaces (Li et al., 2016;Nan and Wonka, 2017;Verdie et al., 2015).This strategy avoids error-prone topological computations during the reconstruction and is able to produce true 3D building models.Although these methods have made certain achievements, their performances are hindered by their own limitations.For example, the box-based fitting method (Li et al., 2016) is applicable only to scenes that satisfy Manhattan assumption.Polyfit (Nan and Wonka, 2017) is a more general method that generates a set of hypothetical faces by intersecting planes extracted by Random Sample Consensus (RANSAC) (Schnabel et al., 2007) and then selected 3D faces that best approximated the building surface using a binary integer programming.However, PolyFit is only intended for the reconstruction of simple polygonal surfaces, and in practical applications, it might encounter computational bottlenecks with complex objects (Nan and Wonka, 2017).
In this paper, an innovative approach for 3D reconstruction of complex buildings from even incomplete point clouds is developed.This method adopts a space-decomposition-andapproximation strategy; but unlike previous methods, the topological relations between the basic space units are extracted after space decomposition.Then, the topological relations are used as constraints in the approximation step to select the optimal space units to form the building models.The topological-relation constraints guarantee the building models to be watertight, manifold, and, to some extent, regularized.Experiments with point clouds containing 105 buildings in Central, Hong Kong, were carried out to evaluate the performance of the proposed approach.

Overview of the Approach
As shown in Figure 1, the proposed approach includes three main steps.In the first step, planar primitives are extracted from the point clouds using RANSAC and refined based on geometric relation constraints.The refined planar primitives are then used to partition the 3D space into a set of polygonal cells with a designed partitioning order through a half binary space partition (BSP) tree (Sohn et al., 2008).In the second step, facets are first extracted based on the 2D topological relations between the cells and then edges are extracted based on the 1D topological relations between facets.The last step selects a set of space units, i.e., cells, facets and edges, that best approximate the building, and this binary labelling problem is solved by integer linear programming (ILP) (Johnson et al., 2000).The fidelity measurements derived from cells and facets and regularity measurements derived from edges formulate the objective function of the ILP problem, and the topological relations between these space units formulate the constraints of the ILP problem.
Figure 1.Overview of the proposed method.

3D Arrangement based on Planar Primitives
The space decomposition generates a set of 3D cells, called cell complex, to present the 3D space that is occupied by the bounding box of a building.These cells have simple and convex geometries and are compactly connected with each other.A subset of the cell complex will later facilitate the watertight and manifold building models without any intersecting computation to determine the edges or corner points.

Extraction and Refinement of Planar Primitives:
A planar primitive P is defined as a bordered plane with a normal vector n and a distance coefficient d.In this paper, a set of initial planar primitives P = {P1, P2, …, Pn} are extracted using RANSAC (Schnabel et al., 2007) and their boundaries are extracted using alpha-shape (Liang et al., 1998).To enhance the regularities between the planar primitives, the initial primitives P are refined based on a set of rules for geometric relationships as shown in Table 1, where nz denotes the unit vector in z-axis.
n xy denotes the projection of normal n on the xy-plane,  and d' are angle and distance thresholds.

Co-planarity
For two parallel plane primitives Pi and Pj, Pi and Pj are co-planar if |di -dj| < d'.
Table 1.Geometric relations used to refine planar primitives extracted from buildings.
To handle the possible regularity conflicts between the planar primitives during the refinement, a geometric relation priority rank is defined as horizontality = verticality > parallelism > orthogonal > z-symmetry = xy-parallelism > co-planarity.
Based on this priority rank, the initial planar primitives are refined on the basis of three plane clusters.They are horizontal planes, vertical planes and oblique planes.

3D Arrangement with a Half-Cut BSP Tree:
After the extraction and refinement of the planar primitives, the bounding box of the building is then partitioned into a set of convex cells with 3D Boolean operations.By each cut, the existing cells, which are parent cells, are split into zero or two child cells, which can be recorded by a BSP tree (Sohn et al., 2008).Instead of partitioning the entire 3D space with each cut as in the studies of Verdie et al. (2015) and Nan and Wonka (2017), which result in redundant cell complex and lead to large computational cost, this method uses a half-cut BSP tree that only partitions the parent cells occupying the planar primitives during each cut, as shown in Figure 2. Because the half-cut BSP tree only partitions parent cells containing the planar primitive by each cut, it is noted that the final cell complex is related to the partitioning order.Thus, a partitioning order is defined as follows.First, vertical planar primitives are designed to have higher priority than horizontal or oblique ones to avoid incomplete partitioning caused by missing data on vertices surfaces, which is actually a common problem in point clouds.Second, in the same priority class, planar primitives with larger areas are considered to have higher priority than smaller ones, aiming to minimize the size of the final cell complex.Note that, planar primitives being too small are considered to be error-prone and are removed during the 3D arrangement.

Extraction of 3D Facets and Edges based on Topological Relations
After 3D arrangement, a set of simple and convex cells can be obtained and is annotated as C = {C1, C2, … Cm}.In this stage, the 3D space is further decomposed into 3D facets and edges based on the topological relations between the cell units.

Extraction of 3D Facets:
The 3D facets refer to the interfaces between the 3D cells.To extract the 3D facets, the bounding faces of each cell are first obtained from their geometries and regarded as initial 3D facets F 1 , which denotes 3D facets that are connected to single cells.For each pair of facets Fi, Fj  F 1 that are co-planar, a new facet Fij  F 2 (F 2 denotes 3D facets that are connected to two cells) is decided whether or not to be created based on the 2D topological relation between Fi and Fj, as shown in Figure 3.If the topological relation is not separated nor connecting, a new facet Fij is created and Fi and Fj are updated.This process is repeated until there are no facets in F 1 producing new Fij  F 2 .

Extraction of 3D Edges:
The extraction of 3D edges uses a similar strategy as the extraction of 3D facets.First, the bounding edges of each 3D facet are obtained and annotated as E 1 , which denotes edges connected to single facets.Then 3D edges E 2 (E 2 denotes edges connected to two facets) are determined based on the 1D topological relations between colinear edges Ei and Ej  E 1 , as shown in Figure 4, through an iterative update process.

Building Approximation with Topological-relation Constraints
Taking each cell, facet and edge as a node, the approximation of the buildings can be considered as a labelling problem that finds the best nodes that form an optimal model.In solving this problem, fidelity energies are encoded in the cell nodes and facet nodes, and regularity energies are encoded in the edge nodes.The global energy function is given as: (1) where l is the binary labelling configuration that assigns each cell, facet and edge a label lC, lF and lE, respectively ( lC, lF , lE  {0, 1}; 1 denotes that the cell/facet/edge is selected, and 0 denotes that it is not).

Cell Fidelity Energies:
Cells inside the building are defined as occupied and those outside the building are defined as empty.As all of the cells are convex and supposed not to bestride the building surfaces, the determination of occupied cells can be simplified as to determine whether their centroids are inside the building.If a ray starting from the centroid of a cell has an odd number of intersection points with the building surfaces, the cell is occupied (Figure 5   For each of cell C  C, a number of rays starting from its centroid are drawn and used to determine the possibility of the cell being occupied based on the number of rays having an odd number of intersection points with the building planar primitives.The cell energies are formulated as: (2) where N(C) is the number of cells in C, lC  {0, 1} is the binary label assigned to C, lC = 1 denotes C is occupied and 0 denotes empty, and pC is the ratio of the number of rays having odd numbers of intersection points with building surfaces to the total number of rays.

Facet Fidelity Energies:
For a 3D facet F  F, if it is supported by points in the building point cloud, it is defined as occupied, otherwise, it is empty.The facet fidelity energies are given as follows: (3) where N(F) is the number of facets in F, lF  {0, 1} is the binary label assigned to F, lF = 1 denotes F is occupied and 0 denotes empty, and pF is the coverage ratio of the area of supporting points to the area of facet F, as computed in Nan and Wonka (2017).

Edge Regularity Energies:
The edge energies are used to introduce regularity constraints about the buildings into the optimal approximation.As man-made objects are generally believed to conform to planarity and orthogonality, the edge energies favour flat and right angles between corresponding facets by setting it a regularity value of 0, and penalize other angles by setting it a regularity value of 1.The edge regularity energies have the following formulation: (4) where N(E) denotes the number of the edges in E, and A(E) is the regularity value set to an edge E  E.

Energy Optimization with Topological-Relation
Constraints: With the cell, facet and edge energies described above, the global energy in Equation ( 1) is presented as a sectional-continuous function.To simplify the ILP problem, the probabilities of the cells and facets being occupied are binarilized into {0, 1} based on two user defined thresholds thresholds C and F (e.g., C = 0.5 and F = 0.3), as shown in Equation ( 5) and ( 6). ( The global energy in Equation ( 1) is therefore turned into the quadratic format as follows: (7) Simultaneously, due to the topological relationships between the cells, facets and edges, the binary variables lC, lF and lE (C  C, F  F and E  E) also meet the following constraints: (8) The first constraint means that for a facet F  F, it is occupied (lF = 1) when only one of its connected cells is occupied (lC = 1), otherwise it is empty (lF = 0).The second constraint means an edge E  E is occupied (lE = 1) when both of its two connected facets are occupied (lF = 1), otherwise it is empty (lE = 0).
With the objective function in Equation ( 7) and the constraints in Equation ( 8), the ILP problem can be solved by the Gurobi solver (Gurobi, 2015).Finally, the building model is formed by the occupied cells, facets and edges labelled as 1

Dataset Description
The dataset used for experimental evaluation is the DIM point clouds acquired in Central, Hong Kong, with a point density of 8 points/m 2 , containing 105 buildings with various architecture styles.Some of the buildings have serious problems of missing data, as shown in Figure 6.

Qualitative Evaluation of the Reconstruction Results
Figure 7 gives an overview of the reconstruction results of the 105 buildings in the test dataset.It can be seen that although the architectural styles of the buildings are very complicated and vary significantly, the reconstruction results are generally consistent with the appearances of the buildings in the point clouds, indicating a promising performance of the proposed method.Figure 8 shows the detailed reconstruction results of four challenging types of buildings, i.e., buildings with complex structures, with missing data, with curved surfaces, and with true 3D structures.Taking Figure 7 and Figure 8 together, it can be noted that the proposed approach performs promisingly in modelling buildings with various architectural styles.Even with the complexity of buildings and poor-quality input data, the proposed approach produces satisfactory building models.

Quantitative Evaluation and Comparisons with Previous Methods
For quantitative evaluation of the geometric accuracy of the models generated from point clouds, the root-mean-square error (RSME) formulated as in Equation ( 9) is used.Figure 9 and 10 show the quantitative evaluations of the reconstruction results of several complex buildings.The results from our approach are compared with the results by three previous methods: Poisson, 2.5D D-C (Zhou and Neumann, 2010) and PolyFit (Nan and Wonka, 2017) and.Poisson and 2.5D D-C methods produced 3D models presented by triangular mesh, while PolyFit and our approach produced polygonal 3D models, which have regularized boundaries, require smaller storage space and can be easily converted into CityGML standards.In general, according to Table 2, the Poisson results have the highest accuracy, but they totally ignored the regularity of the buildings.The enlarged views in Figure 9 shows the details of the reconstructed building rooftops.Results from our approach have the smallest RSME values of 0.67 m and 0.58 m (as shown in Table 2) for the two buildings, indicating a higher modelling accuracy than PolyFit and 2.5D D-C.The main defects of the models generated by our approach occurred at linear and small structures on the building roofs, which can be eliminated during the extraction of the planar primitives of the buildings.Points carved into the building, which might be generated by mismatching during the MVS pipeline, could also cause inaccuracies on the building façades (such as Building-1 shown in Figure 9).However, this problem is consistent with the results of PolyFit and 2.5D D-C.
Although PolyFit shows competitive performance in reconstructing simple buildings with high regularity, it tends to generate inaccurate models for buildings with complex structures (such as Building-3 shown in Figure 10) or just failed to output the building models (such as Building-4 and 5 shown in Figure 10).In contrast, both our approach and the 2.5D D-C show high robustness in the reconstruction of buildings with various complex architectural styles.However, again, 2.5D D-C only produces triangular mesh models with highly irregular boundaries, mostly with relatively lower geometric accuracy compared with our appraoch (as indicated by the RSME measurements of Building-1~4 shown in Table 2).With respect to the modelling of small structures on the building rooftops, 2.5D D-C shows better performance than our approach (as illustrated by Building-3 in Figure 10), because of the failure to present such small structures with planar primitives.However, 2.5D D-C only focuses on the reconstruction of non-vertical structures (e.g., building rooftops), whereas our approach considers the structures of the entire building and has better performance in modelling building façades (as illustrated by Building-4 in Figure 10).
In general, our approach, Poisson and 2.5D D-C have higher robustness than PolyFit, meanwhile, our approach is also able to generate polygonal models with high regularity and true 3D structures, which is beyond the capability of Poisson and 2.5D D-C.Although small structures are likely to be left out in the final models, which can be a common issue for reconstruction methods based on planar segments, our approach, in fact, exhibits outstanding performances considering the limited qualities of the input data.

CONCLUSIONS AND DISCUSSION
In this paper, we propose an innovative and robust approach for automatic generation of 3D building models, regardless the complexity of buildings and the incompletion of point clouds.Qualitative and quantitative experimental evaluations using a challenging dataset in Central, Hong Kong, show promising performances of the proposed approach.The proposed approach provides a practical solution for 3D reconstruction of buildings in large-scale urban scenes, which constitute the spatial data infrastructure for smart cities. Future efforts will be made to investigate the automatic generation of hybrid models, where non-planar structures will be presented as regularized triangular meshes rather than approximated by planes.In addition, converting the polygonal building models into standard format (such as CityGML models) with enriched topological and semantical information will also be investigated, so that the models can be used to support various 3D GIS applications.

Figure 2 .
Figure 2. Space partition with planar primitive based on half-cut BSP tree.

Figure 3 .
Figure 3. Extraction of facets based on 2D topological relations.

Figure 4 .
Figure 4. Extraction of edges based on 1D topological relations.

Figure 5 .
Figure 5. Occupied and empty cells illustrated in 2D.

Figure 6 .
Figure 6.Building point clouds with various complex structures and missing data.

Figure 7 .
Figure 7. Overview of the textured point clouds (first row) and the reconstruction results (second row) of the 105 buildings in Central, Hong Kong.
p -B || denotes the distance between a point p and the geometric model B of a building B, and N (B) is the number of points in the building point cloud .

Figure 8 .
Figure 8. Reconstruction results of four challenging types of buildings.

Figure 9 .
Figure 9. Evaluations of reconstruction results of two simple buildings and comparison with the results of Poisson, 2.5D D-C and PolyFit.Details of the reconstructed building rooftops are shown in the black rectangles, where the input point clouds are overlapped on the building models.

Figure 10 .
Figure 10.Evaluations of reconstruction results of three complex buildings and comparison with the results of Poisson, 2.5D D-C and PolyFit.Details of the reconstructed building rooftops are shown in the black rectangles, where the input point clouds are overlapped on the building models.

Table 2 .
RMSE measurements (m) of building reconstruction results of different methods.