3d Building Reconstruction from Lidar Point Clouds by Adaptive Dual Contouring

This paper presents a novel workflow for data-driven building reconstruction from Light Detection and Ranging (LiDAR) point clouds. The method comprises building extraction, a detailed roof segmentation using region growing with adaptive thresholds, segment boundary creation, and a structural 3D building reconstruction approach using adaptive 2.5D Dual Contouring. First, a 2D-grid is overlain on the segmented point cloud. Second, in each grid cell 3D vertices of the building model are estimated from the corresponding LiDAR points. Then, the number of 3D vertices is reduced in a quad-tree collapsing procedure, and the remaining vertices are connected according to their adjacency in the grid. Roof segments are represented by a Triangular Irregular Network (TIN) and are connected to each other by common vertices or-at height discrepancies-by vertical walls. Resulting 3D building models show a very high accuracy and level of detail, including roof superstructures such as dormers. The workflow is tested and evaluated for two data sets, using the evaluation method and test data of the " ISPRS Test Project on Urban Classification and 3D Building Reconstruction " (Rottensteiner et al., 2012). Results show that the proposed method is comparable with the state of the art approaches, and outperforms them regarding undersegmentation and completeness of the scene reconstruction.


INTRODUCTION 1.1 Motivation
For more than two decades, 3D building reconstruction has been an active research topic of remote sensing, photogrammetry, and computer vision (Rottensteiner et al., 2012;Wang, 2013;Haala and Kada, 2010;Lafarge and Mallet, 2012).Continuing research is driven by the increasing demand for accurate, automatically produced, and detailed 3D city models (Wang, 2013).City models are used for urban planning (Verma et al., 2006), change detection (Rau and Lin, 2011) and environmental or telecommunication simulations (Geibel and Stilla, 2000;Rau and Lin, 2011).Today's utilization of city models expands to everyday user-driven mobile applications, such as location based services (Wang, 2013;Brenner, 2005), 3D Geographic Information Systems (GIS) for navigation, driver assistance systems, virtual tourism (Zhou and Neumann, 2010), and augmented reality.The effort for keeping 3D city models up-to-date depends on the level of automation in building reconstruction.
LiDAR point clouds are well suited for automatic building reconstruction.In comparison to optical stereo imagery, where stereo matching is needed to obtain 3D geometry, LiDAR data contains directly measured, and thus very accurate 3D information (Meng et al., 2010;Haala and Kada, 2010).With continuously increasing LiDAR sensor capacities and point densities, research on building reconstruction has set a focus on LiDAR point clouds (Geibel and Stilla, 2000;Haala and Kada, 2010).

Related Work
Building reconstruction requires the extraction of individual buildings' points from a LiDAR scene.Once buidlings are extracted, there are two main approaches to reconstruction, i.e. model-and data-driven approaches.
Model-driven approaches select for each building point cloud, or parts of it the best fitting parametric model and its corresponding parameters from a predefined catalogue (Maas and Vosselman, 1999;Vosselman and Dijkman, 2001;Kada and McKinley, 2009;Haala and Kada, 2010;Zhang et al., 2012).Model-driven reconstruction is robust, effective and fast, because regularization constraints, such as parallelity and orthogonality, are already inherent in the parametric models.However model-driven approaches are limited to the beforehand defined model portfolio and are therefore not flexible to model all roof shapes.Data-driven approaches connect individual roof segments, which are constructed according to a preliminary segmentation of the building point cloud.Even though data-driven approaches require a high effort for subsequent regularization, they are widely used (e.g.Rottensteiner et al., 2012;Wang, 2013).The advantages of these approaches are a high fit to the input data and flexibility in modeling complex roof shapes.Roof segmentation can be achieved by surface-fitting techniques such as RANSAC (Sohn et al., 2008;Tarsha-Kurdi et al., 2008;Brenner, 2000) or Hough transform (Vosselman and Dijkman, 2001;Sohn et al., 2012;Vosselman et al., 2004), or using region growing methods (Rottensteiner, 2003;Oude Elberink and Vosselman, 2009;Perera et al., 2012;Verma et al., 2006;Nurunnabi et al., 2012;Dorninger and Pfeifer, 2008;Lafarge and Mallet, 2012).Typically, each segment is delimited by a polygonal segment boundary, which is created by using e.g.Alpha-shapes (Dorninger and Pfeifer, 2008;Kada and Wichmann, 2012;Sampath and Shan, 2007;Wang and Shan, 2009), the Voronoi neighborhood (Maas and Vosselman, 1999;Matei et al., 2008;Rottensteiner, 2003) or using a 2D-grid-cell projection (Sun and Salvaggio, 2013;Zhou and Neumann, 2008).Polyhedral 3D models are commonly constructed on the basis of heuristics for extracting and connecting 3D lines along the segment boundaries (Dorninger and Pfeifer, 2008;Vosselman and Dijkman, 2001;Sohn et al., 2008;Rau and Lin, 2011;Rottensteiner, 2003).Structural modeling procedures estimate the coordinates of the building model's 3D vertices by error propagation techniques (Lafarge and Mallet, 2012) or lo-cal error minimization (Fiocco et al., 2005) are less frequently used.The latter technique is also used by Zhou and Neumann (2010), who apply a 2.5D dual contouring algorithm on a point cloud which is segmented into different height layers.The asset of their method is the outstanding flexibility to model complex roof shapes, including non-planar roof segments.However, the algorithm cannot create step edges between roof segments connecting within one roof height layer, which results in a deficiency for modeling superstructures.

METHOD
The proposed workflow (Fig. 1) adapts the method of Zhou and Neumann (2010) for modeling superstructures.The algorithm is modified for a situation-adaptive estimation of the 3D building model's vertices from a detailed roof segmentation.
Input to the procedure are LiDAR data, clustered into sets of Li-DAR points representing different buildings, hereafter referred to as building point clouds.First, the roof points are segmented on the basis of Triangulated Irregular Network (TIN) of the data.Second, a boundary polygon is created for each segmented cluster.Third, vertices of the 3D building model are estimated and connected using an adaptive 2.5D Dual Contouring procedure.Regularization for enhancing model simplicity is not included in this is work.4)).

Roof Segmentation
The goal of roof segmentation is to cluster the points according to the roof segments they belong to.While most region growing (RG) segmentation techniques assume roof segments to be planar, the proposed segmentation is designed to allow any continuous shape.For this purpose, a robust TIN-based RG technique is proposed.Moreover, in contrast to most RG techniques (e.g.Dorninger and Pfeifer, 2008;Oude Elberink and Vosselman, 2009;Sampath and Shan, 2010;Sun and Salvaggio, 2013), the proposed method does not assign one segment label to each LiDAR point, but to each triangle of the TIN.Thereby, LiDAR points, which are part of differently labelled triangles, have more than one label.The labelling of triangles minimizes gaps between adjacent intersecting roof segments and allows accurate boundary determination.
The RG procedure iteratively starts at seed triangles defined by a minimum Local Unevenness Factor (LU F ), defined as where A k is equal to the area and n k = [n x,k , n y,k , n z,k ] T is the normal of the k-th of Kt triangles which are in the neighborhood of triangle t; nx/y/z are the means of all n x/y/z,k .For RG, each triangle is tested for a fixed threshold on the local angular deviation and for two adaptive thresholds.These two robust adaptive thresholds are built according to Nurunnabi et al. (2012) and are using the LU F and the LiDAR points' distances to the current best fitting segment plane.

Segment Boundaries
For each point cloud representing a roof segment, a polygonal boundary is created in an iterative convex-hull collapsing procedure.Iteratively, each line segments of the convex hull is refined by the LiDAR point with a minimum distance measure.The refinement stops when the line segment is shorter than a directionally dependent threshold, which is created by considering the LiDAR point spacings in across-track and along-track sampling directions.

Building Modeling
The proposed modeling algorithm estimates and connects the 3D vertices of the building model using an adaptive 2.5D Dual Contouring procedure.In section 2.3.1, the 2.5D dual contouring principle is introduced.A 2D grid is overlain on the data, and grid data is computed (section 2.3.2).In an iterative quadtree collapsing procedure, a Quadratic Error Function (QEF) is constructed from the grid data of each four adjacent grid cells.The minimization of the situation-adaptive QEF results in the coordinates of one or more 3D vertices of the building model, depending on whether the vertices represent a step-or an intersection edge (section 2.3.3).

Dual Contouring principle
For building reconstruction, a 2D grid is overlain to the segmented LiDAR points in the x-yplane.The Dual Contouring principle can be illustrated by the example of estimating the vertices of boundary polygons, which separate the segments in the 2D-plane.In each grid cell, a polygon vertex is estimated by minimizing its distances from local boundary lines, which separate the LiDAR points of different segments.Then, a polygon is created by connecting the vertices of adjacent cells (Fig. 2).
The purpose of 2.5D Dual Contouring is to estimate vertices of the 3D building polygon in each grid cell, which are described by so-called hyperpoints.Depending on whether a hyperpoint X = x y zv zv+1 ... zV T describes a step edge or an intersection edge, it contains two or more 3D vertices with the same x-y-coordinates, but with different z-coordinates.Each zcoordinate defines a 3D vertex X3D,v = [x, y, zv], in which all segments S k within a local height layer Hv intersect (Fig. 3).
A hyperpoint's optimal x and y-coordinates minimize the 2Ddistances E2D(X) to local boundary lines (LBL) between the respective segments in the 2D plane.E2D(X) is computed as (2) where ni (k,l) is the normal of the i-th local boundary line LBL k,l between the segments S k and S l , and p i (k,l) is point on this line.
Additionally, an optimal 3D vertex X3D,v corresponding to Hv, v = (1, ..., V ) minimizes the 3D-distances E3D(X) to the local surface planes (LSP), which are fitted to the LiDAR points belonging to each segment S k , k = (1, ..., M ).E3D(X) is computed as where mj k is the j-th normal of the local surface plane LSP k on segment S k , and q j k is a point on this plane.

Local grid data
For each vertex of the 2D grid, a local surface plane LSP = [m, q] is determined, where m = [mx, my, mz] T is the plane's normal, and q = [qx, qy, qz] T is a point on the plane.q xy = [qx, qy] T is equal to the grid vertex.
Each LSP is associated with a segment label lLSP , according to the segment, which is the closest to q xy (Fig. 4 b).Vector m is determined by averaging the normals of the K nearest TINtriangles belonging to S k .
For each grid cell, a local boundary line LBL k,l = [n, p] is estimated for each pair of segments S k and S l using a Least Squares approach (Fig. 4 a).LBL k,l is estimated from all LiDAR points belonging to S k and S l , which are within a buffer zone around the grid cell.LBL which have no intersection point p with the grid cell's border are discarded.

Adaptive
QEF In contrast to Zhou and Neumann (2010), whose building point cloud is segmented into different height layers, the presented method works on a detailed segmentation of the roof into diffent segments.While Zhou and Neumann (2010) estimate one z-coordinate for each global roof height layer, (eq.4), the proposed method requires that the LiDAR points within one cell are grouped into local height layers.The advantage of local height layers is that step edges can be created between segments from one global roof height layer.This allows to model complex roof structures such as dormers and shed roof segments.Grouping the segments into local height layers Hv is achieved by estimating a step edge probability SEP for each local boundary line LBL (Eq.5).Within one cell, all segments with an SEP < 0.5 are grouped into one local height layer Hv.Assuming the z-coordinates of Li-DAR points to be normally distributed around their true value, the equation for computing SEP is designed to use the minimum step edge height Tstep as standard deviation: where Tstep is a fixed step edge threshold and dz a measure expressing the local height difference of the two segments.
In Zhou and Neumann (2010), all local boundary lines LBL represent step edges, as their input point cloud is only segmented into global roof height layers.When constructing the QEF, it has to be considered that in the proposed method, LBL can also represent intersection edges.In case of a step edge, ideally only the local boundary lines LBL are considered for estimating a hyperpoint's horizontal position [x, y].The distances to local surface planes LSP should not be considered.As LSP cannot be omitted, a balancing weight w k,l is computed for each group of ni k,l , i = 1, ..., I k,l using the corresponding SEP (k,l) .Each w k,l ranges from [0, ..., 1, ..., wmax], corresponding to a SEP (k,l) of [0, ..., 0.5, ..., 1], where wmax is the maximum weight.
If the number of LSP for each roof segment represented in a QEF is not equally distributed, QEF minimization will result in a distortion.For instance, minimizing the E3D(X) in a QEF containing three LSP k and one LSP k+1 will not minimize the point's distances to both roof segments S k and S k+1 .In order to weight each roof segment equally, each normal mj k is scaled by the number N k of corresponding LSPk.The weighted and scaled QEF is defined as 2.3.4QEF solution and quadtree collapsing Each QEF is solved by least squares adjustment after composing a matrix equation Xw = argmin where A is the model matix, b is the vector of observations and W is the vector containing the weightings of each QEF line, w (k,l) , k, l = 1, ..., M , k = l and 1 N k , k = 1, ..., M .Additional solution constraints ensure that the QEF solution lies inside the quadtree cell.
The number of hyperpoints for a building model should me minimized.Therefore, the grid cells are treated as leaf cells of a quadtree, which is iteratively collapsed.For this purpose, the grid is designed to have 2 n cells.By this, an iterative collapsing of groups of 4 adjacent cells into one larger quadtree cell is possible.For deciding whether to collaps a group of four quadtree cells, a combined QEF is constructed from the LSP and LPL of these cells.If the non weighted residual error RQEF = AX − b is larger than a threshold Rmax, the four quadtree cells are collapsed to a larger one.The procedure iterates until there is no group of four quadtree cells which can further be collapsed.

Building polygon creation
Each hyperpoint vertex carries labels for the segments of the local height layer it is estimated from.Therefore, within each pair of adjacent hyperpoints, there is a pair of hyperpoint vertices sharing at least one segment label.Those hyperpoint vertices are connected to form a 3D edge.After creating a 3D edge between all adjacent hyperpoints, the building roof is represented by to two types of connections, i.e. 3D triangles and 3D quads of 3D edges (Fig. 5).In order to represent the roof by a triangulation, each quad is separated into two triangles.If two possibilities for separation are possible (Fig. 5 a), the separation resulting in the best fit to the input point cloud is chosen.Where a 3D edge is only part of one 3D triangle (single edge), a vertical wall has to be created between this 3D edge and another single edge from the same pair of hyperpoints.If no such other single edge is found, a vertical wall is created to ground.The result is a watertight 3D building model, represented by a large number of roof triangles and vertical wall elements.Subsequent regularization procedures are recommended for increasing model simplicity, but are not in the scope of this paper.

TESTS AND EVALUATION
In this section, tests of the proposed method, data and results are described.The results are evaluated using the method applied in the "ISPRS Test Project on Urban Classification and 3D Building Reconstruction" (ISPRS benchmark project, Rottensteiner et al., 2012) and compared to the results of other state of the art methods.

Test data and parameters
Two different datasets are used for testing, a smaller scene from Munich, Germany, and a larger scene from Vaihingen, Germany.The Vaihingen scene corresponds to the test scene "Vaihingen Area 1" of the ISPRS benchmark project.

Test results
All the buildings in both datasets were reconstructed.Figure 7 a presents the reconsted city scene of Vaihingen, consisting of 21 buildings.Figure 7b presents the reconstructed Munich scene with 8 buildings.

Evaluation
The results of the reconstructed buildings are evaluated according to the evaluation of the the ISPRS benchmark project (Rotten-    steiner et al., 2012).The ground truth used by the ISPRS benchmark project was not made available.Therefore, the results are evaluated against manually extracted 2D reference segment polygons.For the Vaihingen data set, ground truth was extracted from the ortho image delivered with the test data for the ISPRS benchmark project.For the Munich data set, ground truth was extracted from high resolution ortho image.The following eight evaluation parameters are calculated: ) or less than (F N ) 50 % with estimated segment polygons.Cm,10 is analogously for segment areas ≥ 10m 2 .
• Correctness Cr = T Pe T Pe + F P , where T P d are true positive and F P are false positive estimated segments, whose area is ≥ 2.5 m 2 and which are overlapping by at least (T P d ) or less than (F P ) 50 % with reference polygons.Cr,10 is computed analogously for segment areas ≥ 10m 2 .
• RMSE computed from the 2D distances dxy of the estimated segment outline vertices to their reference segment, while dxy > 3m are neglected.
• NO: Number of oversegmented reference segments, i.e. those corresponding to more than one T P d .
• NU: Number of undersegmenting estimated segsments, i.e. those corresponding to more than one T Pr.
• NC: Number of references which are both under-and oversegmented.

Evaluation results
Tables 2 and 3 show the evaluation parameters for segmentation.2D segment outlines are evaluated against the reference segments.Results for the Vaihingen scene are compared to those of Awrangjeb and Fraser (2014), who apply the same evaluation method (section 3.3) to their segmentation results.Awrangjeb and Fraser (2014).

Method
Tables 4 and 5 show the evaluation results for reconstruction.The outer boundaries of each group of adjacent roof triangles with equal segment label are projected to the x-y-plane and evaluated against the reference segments.Results for the Vaihingen scene are compared to those of the ISPRS benchmark project.
It has to be considered that the manually extracted ground truth might differ from the ground truth used by the ISPRS benchmark project.

DISCUSSION
The proposed method creates building models of high detail by triangulation.Rottensteiner et al., 2012).
Depending on the chosen level of detail (by setting the parameters C and R), points for building representation are reduced to a minimum at continuities.Second, step edges are always represented as vertical walls, as vertices of a hyperpoint are always vertically arranged.
The novelty of the proposed reconstruction method is using a detailed segmentation as input to a 2.5D Dual Contouring approach.
In contrast to Zhou and Neumann (2010), the proposed situationadaptive QEF construction allows to model accurate step edges between segments from one height layer (Fig. 8).For each pair of local height layers within each grid cell, a decision is made whether to connect them either by an intersection edge or by a step edge.The step edge threshold Tstep, the residual for quadtree collapsing R, and the grid cell size C are of equal importance for deciding whether a step edge or an intersection edge is locally created.If C is chosen too high, roof segments which require both intersection edges and step edges to be modelled accurately (e.g.dormers) will only fall into one height layer and are "smoothed out".The same smoothing effect happens if R or Tstep are chosen too high, or if roof segments within one roof heigh layer are not segmented accurately, i.e. undersegmentation occurs (Fig. 9).(c) Reconstruction, where no step edge can be modelled, because no different segments are identified.
Figure 9. Smoothed roof effect due to undersegmentation.If the roof height layer is undersegmented no step edges can be created, because only one local height layers is identified.However, the roof is approximated to fit the input point cloud in detail.more roof triangles.However, if a regularized building model is required the large number of roof triangles increases the effort for postprocessing.
The building model's 3D edges are constructed from connections of adjacent hyperpoints in the quadtree.In certain cases, this can lead to erroneous edges (Fig. 10b, right arrow).If hyperpoints from non-adjacent quadtree cells are connected to an edge, deformations or "stand-alone walls" in the model can occur (Fig. 10b, left arrow).

CONCLUSION
A novel method for creating detailed building models with complex roof shapes from LiDAR point clouds is proposed in this paper.The 2.5D Dual Contouring method of Zhou and Neumann (2010) is used and adapted in a way that step edges and intersection edges can be created between roof segments.A main contribution of this work is the modification and weighting of the Quadratic Error Function (QEF) for modeling step edges and intersection edges.The modeling depends on the step edge probabilities of local height layers.A prerequisite for adaptive 2.5D Dual Contouring is a roof segmentation technique which stops at smooth edges.The applied robust TIN-based region growing reliably stops at smooth edges.Consequently, undersegmentation is significantly reduced.The resulting building models show a very high fit to the input LiDAR points.Each roof segment is represented by a triangulation, thus also non-planar roof shapes can be modelled.Subsequent model regularization is recommended, because buildings are represented by a large number of vertices.Errors in reconstruction result mostly from wrong or missing connections of the vertices.Thus, the way the connections of the vertices to the building model should be more robust.Wrong connections could be avoided by checking for the consistency of the model with the building footprint.Under assumption that building edges are mostly orthogonal or parallel to the main build-   ing direction, the missing connections could be avoided by aligning the building point cloud to the main building direction before the modeling procedure.The proposed workflow was tested and evaluated using two data sets with different characteristics and varying building complexity.Evaluation of the results has shown that both segmentation of LiDAR point clouds and reconstruction of buildings are comparable to the state-of-the-art methods.
(a) 2D-grid overlain to segmented point cloud.Each colour represents one segment.(b) Blue lines represent local boundary lines for each pair of segments in a cell.(c)QEF minimization (blue cell) using all local boundary lines from adjacent cells (pink lines).The yellow star is the QEF solution, i.e. one optimal vertex of the polygonal boundary line.(d)QEF solutions (yellow stars) are computed for all cells and connected to the polygonal boundary line (black line) according to their adjacency in the grid.

Figure 2 .
Figure 2. 2D Illustration of the Dual Contouring principle.The input is segmented point cloud (a) and the QEF solutions are the polygonal boundary lines (d).

Figure 3 .
Figure 3. Hyperpoints at intersection edges (one z-coordinate) and step edges (more than one z-coordinate).
Local boundary lines plotted over the segmented LiDAR points.Each color represents one segment.The lines' colors indicate diffferent step edge probabilities, ranging from zero (blue) to one (red).
Points q and normal vectors m of the local surface planes, plotted over the segment boundaries.Colours correspond to the segment labels l SP .A side view of combined grid data, i.e. local boundary lines LBL and local surface planes LSP.

Figure 4 .
Figure 4. Local boundary lines and local surface planes in a grid.
Figure 5. Quads (yellow areas) and triangles (red areas) resulting from connecting adjacent hyperpoint vertices (dark blue points) to 3D edges (blue lines, top view) point cloud of the Munich test scene.

Figure 6 .
Figure 6.LiDAR point clouds of the test scenes.
Side view of the reconstructed Vaihingen scene, wall segments (grey areas) are connecting the outer roof boundaries to ground or to other roof segments.Side view of the reconstructed Munich scene, wall segments (grey areas) are connecting the outer roof boundaries to ground or to other roof segments.

Figure 8 .
Figure 8. Detailed building model in a local coordinate system.Tstep = 0.2 m, grid cell size C = 2.5 LiDAR points per cell, residual threshold R = 0.8 m.Segments from one height layer are separated by both, step edges and intersection edges.
where T Pr are true positive and F N are false negative reference polygons, whose area is ≥ 2.5 m 2 and which are overlapping by at least (T Pr ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W4, 2015 PIA15+HRIGI15 -Joint ISPRS conference 2015, 25-27 March 2015, Munich, Germany This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.doi:10.5194/isprsannals-II-3-W4-157-2015 One hyperpoint (dots) computed for each cell of the collapsed quadtree.Erroneously connected hyperpoints of neighboring quadtree cells (right arrow), and of non-neighboring hyperpoints (left arrow).
3D model resulting from erroneous hyperpoint connection.An additional (wrong) rooftop triangle due to wrong hyperpoint connection (right part), and a "standalone wall" due to missing hyperpoint connection (left).

Figure 10 .
Figure 10.Errors due to connectivity of hyperpoints from (non)neighboring quadtree cells.

Table 1 .
Table 1 shows the data characteristics of the two test scenes and the applied reconstruction parameters.Characteristics of the two test scenes and applied parameters for testing the proposed method.

Table 2 .
Awrangjeb and Fraser (2014)s of the segmentation.Segmentation results of the Vaihingen Scene are compared to the results ofAwrangjeb and Fraser (2014).

Table 3 .
Under-and oversegmentation and accuracy of the segmentation.Segmentation results of the Vaihingen Scene are compared to the results of

Table 4 .
Rottensteiner et al., 2012)g the input LiDAR point cloud, adaptive 2.5D Dual Contouring has two main advantages: Completeness and correctness of the reconstruction.Results of the Vaihingen Scene are compared to the results of other methods (* seeRottensteiner et al., 2012).

Table 5 .
Under-and over-segmentation and accuracy of the reconstruction.Results of the Vaihingen Scene are compared to the results of other methods (* see