Robust and fast reconstruction of complex roofs with active sampling from 3D point clouds

This article proposes a novel method for the 3D reconstruction of LoD2 buildings from LiDAR data. We propose an active sampling strategy which applies a cascade of filters focusing on promising samples at an early stage, thus avoiding the pitfalls of RANSAC‐based approaches. Filters are based on prior knowledge represented by (nonparametric) density distributions. In our approach samples are pairs of surflets—3D points together with normal vectors derived from a plane approximation of their neighborhood. Surflet pairs provide parameters for model candidates such as azimuth, inclination and ridge height, as well as parameters estimating internal precision and consistency. This provides a ranking of roof model candidates and leads to a small number of promising hypotheses. Building footprints are derived in a preprocessing step using machine learning methods, in particular support vector machines.


Abstract
This article proposes a novel method for the 3D reconstruction of LoD2 buildings from LiDAR data. We propose an active sampling strategy which applies a cascade of filters focusing on promising samples at an early stage, thus avoiding the pitfalls of RANSAC-based approaches. Filters are based on prior knowledge represented by (nonparametric) density distributions. In our approach samples are pairs of surflets-3D points together with normal vectors derived from a plane approximation of their neighborhood. Surflet pairs provide parameters for model candidates such as azimuth, inclination and ridge height, as well as parameters estimating internal precision and consistency. This provides a ranking of roof model candidates and leads to a small number of promising hypotheses. Building footprints are derived in a preprocessing step using machine learning methods, in particular support vector machines. approach integrating additional bottom-up information into the top-down process is proposed.
In a former article (Henn, Gröger, Stroh, & Plümer, 2013), we applied the random sample consensus (RANSAC) method (Fischler & Bolles, 1981) to fit semantic roof models in 3D point clouds. RANSAC is a robust estimator and produces good estimations as long as the noise is randomly distributed and unstructured, which is not guaranteed in the context of roofs due to highly structured noise such as dormers or chimneys. In fact, RANSAC nowadays represents a large family of algorithms, the most popular being MLESAC (Torr & Zisserman, 2000), but most of the members of this family are hardly able to cope with situations where the observations are explained by different similar, but competing, hypotheses (such as parallel planes with small distances or planes with the same azimuth but slightly different inclinations as in the case of dormers) and the searched model is supported only by a small proportion of the observations (see Section 3 for more details).
In RANSAC-based approaches the size of the support set (or the score as in MLESAC which maximizes the likelihood rather than the number of inliers) has been used for selecting the best model. As will be illustrated in Section 2, this approach fails if the observations can best be explained by structured models aggregated from several primitives, such as adjacent, parallel gabled roofs or roofs with large shed dormers.
The contribution of this article is a new method for deriving LoD2 buildings with larger roof structures from LiDAR data, which overcomes these and similar problems arising in other approaches. We propose an active sampling strategy that excludes "bad;" samples (which lead to models not supported by the observations) from the beginning. Prior knowledge about models and good samples, represented by probability distributions of parameters, is used to focus on promising hypotheses at a very early stage, in contrast to the random sampling of RANSAC. This approach has been inspired and motivated by probabilistic robotics where this kind of importance sampling has successfully been applied (Thrun, Burgard, & Fox, 2005). We use the term "active sampling" in contrast to "random sampling" typical of RANSAC. Active sampling is similar to importance sampling but not the same. It uses prior knowledge in order to assess samples at a very early stage in order to reject weak samples as early as possible and simultaneously avoid early commitment to any hypothesis. In this way, we actively focus on promising samples and corresponding hypotheses. To the best knowledge of the authors, such methods have not been employed in the context of building reconstruction and especially roof reconstruction until now, with the exception of Nguatem, Drauschke, and Mayer (2013) and Schmittwilken and Plümer (2010) who used an informed sampling for façade reconstruction.
With regard to roof models, we distinguish between generic and extension parameters. The extension parameters (bounding geometry) are given by a polygon, with a rectangle as a special case, whereas the generic parameters may be specified by plane equations. For a generic symmetric gabled roof model, four points in general position are sufficient to generate hypotheses (Henn et al., 2013). In this article, we use pairs of surflets (Wahl, Hillenbrand, & Hirzinger, 2003) instead. An observed surflet consists of a 3D point (from the LiDAR point cloud) and its normal representing a plane approximation of its k-neighborhood in the point cloud. Figure 1 summarizes the surflet-based classification of gabled roofs. Figure 1a depicts observed surflets in a LiDAR point cloud that represent a pair. Such pairs are filtered based on roof properties such as planarity, curvature, symmetry, and opposite azimuth, and consequently represent preliminary model hypotheses which can be easily verified or falsified afterwards. Figure 1b shows a top view of the point cloud from Figure 1a which has been segmented according to the azimuth of the underlying surflet points. The surflet pairs restrict the model hypotheses space enabling the reconstruction of the roof models. Figure 1c depicts a top view of the reconstructed 3D roof model. The corresponding roof model surflets are drawn in the same color as the roof surface.
The verification of the prediction for gabled roofs with regard to the observations is reduced to the verification of a ridge prediction. Verification or falsification of hypotheses for gabled roofs can be based on predicted ridges and the 3D points observed in the vicinity of the predicted ridges. Hence, our method uses a hybrid approach integrating additional bottom-up information into the top-down process, as proposed by Huang et al. (2013). Pairs of surflets integrate the bottom-up (i.e., data-driven) approach in a very tight way with the top-down (i.e., model-driven) approach. They are tailored to specific roof models such as gabled roofs and both model parameters (azimuth, declination, ridge height) and measures of consistency (mean square error, curvature, symmetry).
Many approaches for automatically deriving LoD2 building models without larger roof structures; (e.g., Henn et al., 2013;Zhang, Zang, Agam, & Chen, 2014) require building footprints as additional input. Such footprints typically exist in most highly developed countries. However, even if such data are available, there might be frequent inconsistencies since often the footprints and the other input data (LiDAR, aerial images) have been collected at different points in time, resulting in buildings represented in the LiDAR/image data, but missing in the footprint data, or vice versa.
There are also relevant disaster scenarios where current footprints are not available. Consider, for example, the Haiti earthquake in 2010, where the damage assessment could not be based on reliable footprint information.
F I G U R E 1 Efficient classification of gabled roofs by surflets. Pairs of surflets in a point cloud (observed surflets) are filtered by using roof properties (planarity, curvature, symmetry, opposite azimuth): (a) A surflet pair/ model hypothesis (green, blue) derived by filtering, front view of the point cloud; (b) Segmentation of the noisy point cloud according to the azimuth of surflets, top view; and (c) Reconstructed roof model and corresponding roof model surflets Our approach does not require building footprints as input. They are derived in a preprocessing step using machine learning methods, in particular support vector machines (SVMs), in combination with an α-shape based boundary determination. However, if they are available, they can be incorporated and used to improve the result.
The remainder of this article is organized as follows. The next section gives an overview of related approaches to generating buildings from LiDAR data and their shortcomings. Our new method is presented in Section 3. Section 3.1 is devoted to the classification of LiDAR points by methods from machine learning as well as the estimation and regularization of footprints from the points classified as building points. The methods for detection of roofs and larger roof structures using active sampling and surflets are introduced in Section 3.2. Section 3.3 deals with the 3D reconstruction of structured roof models and discusses the reconstruction results. Section 4 provides concluding remarks and an outlook on the next steps in building reconstruction.

| REL ATED WORK
In the field of automatic building and especially roof reconstruction from LiDAR data, two approaches are mainly used: bottom-up (nonparametric) and top-down (parametric). A comparison of both approaches can be found in Haala and Kada (2010). Model-driven approaches (e.g. Henn et al., 2013;Kada & McKinley, 2009;Poullis & You, 2009) define a catalog of parametric 3D roof models a priori and fit the roof models to the data, followed by a model selection (e.g. minimum description length, Akaike information criterion, SVM). The approach of Kada and McKinley (2009) uses normal vectors for the determination of the models which best fit the point cloud. Huang, Brenner, and Sester (2011) and Lafarge, Descombes, Zerubia, and Pierrot-Deseilligny (2010) use reversible jump Markov chain Monte Carlo to fit the primitives. Gaussian mixture models and the EM-algorithm are used by Poullis and You (2009) for this task. Besides the detection of locally fitted primitives, Zhou and Neumann (2012) detect relations between roof and wall polygons, such as orientation, placement, and parallelism as well as orthogonality and equality of roof planes. Top-down approaches regularly have to deal with the problem of the flexibility. If, however, roof primitives are aggregated in a regular, constrained way as specified by formal grammars, for instance, they are able to model buildings by aggregating primitives and using specific constraints based on frameworks such as attribute grammars and Markov logic networks (Dehbi, Gröger, & Plümer, 2016;Dehbi, Hadiji, Gröger, Kersting, & Plümer, 2017).
In data-driven approaches, geometric primitives such as planes are segmented in a first step and aggregated afterwards. Among others, popular methods are RANSAC (Tarsha-Kurdi et al., 2008), 3D Hough transform (Overby, Bodum, Kjems, & Ilsø, 2004;Vosselman & Dijkman, 2001), graph matching (Oude Elberink, 2009) and region growing (Rottensteiner, 2006;Verma, Kumar, & Hsu, 2006). A further data-driven approach has been presented by Poullis (2013). Laser points are clustered first according to the point normal and height variance (based on an 8-neighborhood). Points are clustered (unsupervised) into patches, and neighboring patches are aggregated based on a likelihood function. Finally, the boundary is derived and refined, by employing an energy minimizing function using graph cuts. Data-driven approaches are usually faced with data perturbation such as occlusions or reflections, leading to disturbances such as incompletely reconstructed roofs. Furthermore, geometric constraints have to be satisfied. Hence, a postprocessing step, as described in Rottensteiner (2006), is needed which may lead to inconsistency with the previous fitting step (Huang et al., 2013).
A hybrid approach for building reconstruction from LiDAR data combining bottom-up and top-down strategies has been proposed by Satari, Samadzadegan, Azizi, and Maas (2012). Main roof planes are detected first in a bottom-up manner, using point-based classification methods, clustering and the α-shape algorithm. Afterwards, the derived roof planes are verified by supervised machine learning methods (SVM with polynomial kernel); features are the angle and gradient differences as well as the Euclidean distance between LiDAR points and the plane. In the second, model-based step, dormers are reconstructed using dormer models and constraints (e.g., identical azimuth and inclination difference of more than 20 • between main roof plane and dormer roof). For the classification of dormer types, again an SVM is employed. Another hybrid approach was proposed by Lafarge and Mallet (2012). In a first step, building points are classified using energy minimization and graph cut based methods, based on local point features. In a model-based way, geometric primitives (planar, spherical, cylindrical, and conoidal shapes) are fitted to the points classified as building points. This segmentation is derived by iterative nonlinear optimization. For points not covered by primitives, meshes are determined in a data-driven manner. In the third step, the primitives and meshes are combined/aggregated using an adjacency graph. Similar approaches target also the reconstruction of façades (Lafarge, Keriven, Bredif, & Vu, 2013). For the derivation of building footprints from LiDAR points, Huang and Sester (2011) apply a hybrid approach. A segmentation (data-driven) of the points into planes is achieved by a 3D Hough transformation. Based on the segmentation, primitives (rectangles/triangles) are fitted to the point cloud, allowing for overlaps along the primitives. The method used is reversible jump Markov chain Monte Carlo.
In order to estimate geometric primitives, bottom-up as well as top-down approaches make use of RANSAC (Beder & Förstner, 2006;Schnabel, Wessel, Wahl, & Klein, 2008;Tarsha-Kurdi et al., 2008). In Henn et al. (2013) RANSAC is used to directly fit roof models. Geometric primitives such as planes together with the underlying constraints are estimated in a model-driven way. In order to derive roof model parameters needed for RANSAC, this method requires an appropriate decomposition of the building footprints into rectangles which is not given in many cases (e.g., for complex buildings). RANSAC is a robust estimator even though the proportion of model outliers (observations which do not support the model) is more than 50% of the observations. Furthermore, RANSAC produces good estimates as long as the noise is randomly distributed and unstructured, which is not guaranteed in the context of roofs due to highly structured noise such as dormers or chimneys. However, RANSAC fails to deliver good results if the searched model is supported only by a small proportion of the observations. Examples are given in Figure 2 where the gabled roofs are not detected since the support of bigger, false models is dominant.
This problem is also present for dominant dormers (Figure 2c).
These problems are a consequence of the sampling and scoring strategy of RANSAC. It searches randomly for samples, generates models, and then scores them using the sum of the residuals. This is done without taking the structure and the distribution of the residuals into account. For instance, four points in general position are needed in order to sample (construct) a gabled roof (Henn et al., 2013). However, this choice does not guarantee a good model prediction. We make use of surflets in order to enable a 3D building reconstruction including larger roof structures. Wahl et al. (2003) construct surflet-pair relation histograms as an object descriptor for classification F I G U R E 2 RANSAC: problems in estimating the roof model for observations (observations supporting the model are depicted in red, those that do not support the model in yellow); (a) Expected models (gabled roofs) are contained in the model space of RANSAC. Since the model is supported only by a part of the observations, it delivers a false estimation; (b) A false model (flat roof covering the whole observations) obtains more support than the smaller gabled roofs; and (c) Large dormers mislead RANSAC into making a false estimation of a gable roof in 3D point clouds. Rusu, Blodow, and Beetz (2009) modify the previous descriptor in order to reduce their computation times, whereby the surflets in histograms are solely considered in their local neighborhood. However, Papazov and Burschka (2011) take a fraction of surflets into account in order to detect their underlying models.
This certainly leads to a loss of information about the object shape.
Strategies for efficient sampling have been introduced in the context of robotics (Thrun et al., 2005). This kind of sampling is known as importance sampling and computes for every point a fitness value, a "probability of being a good point," using nonparametric density estimation techniques. In the field of building reconstruction, Schmittwilken and Plümer (2010) and Schmittwilken (2012) estimate stairs and windows from terrestrial laser scans. Sampling strategies for RANSAC have been adapted in order to search for samples in a terrestrial laser point cloud that most likely represent the object of interest. In particular, stochastic universal sampling (Baker, 1987) has been used.

| DE TEC TI ON AND RECON S TRUC TI ON OF COMPLE X 3D ROOFS
This section gives an overview on the surflet based active sampling of roof hypotheses as the first step for reconstructing complex 3D roofs. Section 3.1 deals with the classification of LiDAR 3D point clouds, and the estimation and regularization of footprints from the classified points. Section 3.2. describes the approach for detecting roofs and roofs structures which will be reconstructed afterwards. A summary of the whole method is given in Algorithm 1. This is a synthesis of the methods which will be demonstrated in the following sections.
In order to overcome the deficits related to random sampling, we propose a novel sampling strategy which applies a cascade of filters and rejects bad samples as early as possible. This approach covers different roof types and will be demonstrated for gabled roofs. Based on the point normal of each surflet s i , points of flat roofs have been determined via a threshold of 10 • for the angle of the normal to the z-axis.
The first property of roofs which is used for the selection of samples is planarity. Thus a point is a good candidate if its neighborhood with regard to the observations is approximately planar. One measure of planarity is the mean square error (MSE) of its k nearest neighbors with regard to their approximating plane. A second, slightly different measure is the curvature κ (Pauly, Gross, & Kobbelt, 2002), which is calculated from the eigenvalues i of the principal component analysis of that neighborhood via: The second roof property is that ridges are parallel to the ground (x-y plane). Ridges are guaranteed to be horizontal if the azimuths of the opposite roof planes differ by exactly π. We consider points together with their normals, as mentioned in the introduction. This pair (point and normal vector) is well known in computer graphics as a surflet (Wahl et al., 2003). We assume that the normal vectors are normalized (unit length) and oriented upward.
Upward orientation has the immediate consequence that roof pitches are positive and roof planes are oriented.
Our approach starts with a given LiDAR 3D point cloud AllPts and delivers a reconstructed boundary representation Brep of the reconstructed roof models. During the first preprocessing step (see Section 3.1) each point is classified as "building", "vegetation", or "ground". Subsequently building areas and their (preliminary) Outlines are detected and regularized. Also, the classified points Pts ⊂ ℝ 3 and the probability  that a given point is a part of a building are provided. In the next step, surflet-based active sampling, demonstrated in Algorithm 3, generates preliminary roof hypotheses  together with best-ranked surflets  b ⊂ ℝ 3 × ℝ 3 as well as best-ranked surflet pairs . Subsequently, a model selection of 3D roofs is performed, leading to an initial set of roof models  0 (see Algorithm 4). Afterwards the previous predictions of roofs and roof parts are refined in order to discriminate between and combine roof types and roof parts, resulting in refined roof and dormer models  and , respectively.
Finally, the roofs detected from the previous step as well as their parts such as dormers are reconstructed in a topologically correct manner.

| Preprocessing and detection of buildings in 3D point clouds
A coarse knowledge of the class of each point ("building", "vegetation", and "ground") improves our sampling method, as will be described in Section 3.2. The class information (category and probability) of the points labeled "building" (henceforth called "building candidates") reduces the search space in sampling. This allows surflets to be solely sampled from roof points. Furthermore, the building outlines are used to bound building model hypotheses to a local area. A refined bounding geometry will be acquired based on the inliers of the selected model and used for the reconstruction of the entire building geometry. For that purpose, we applied a point based classification to the 3D point cloud with a subsequent determination of boundaries for the building candidates. These can be the boundaries of detached buildings or of connected building complexes. The algorithm for our preprocessing is given in Algorithm 2.

| Classification of 3D point clouds
An SVM was used for the classification of points. We applied features from Chehata, Guo, andMallet (2009), Niemeyer, Rottensteiner, andSörgel (2012), and Zhou and Neumann (2008)  single point pt i by its position, its normal and the distribution of positions and normals of its k-neighborhood. Since the information on the ground surface (e.g., a digital terrain model) is not available in many cases, the height of the ground has to be determined from the point cloud. Therefore, a feature for the approximation of height difference to the ground surface, as proposed by Chehata et al. (2009) has been calculated. They assume that the point with the lowest z-value in a vertical cylinder centered at the point pt i with radius r is the point on the ground surface.
In contrast, we sort the points by their z-component and select the subset L of the lowest 5% of the points. Our approximated height h gr of the ground surface is the median height h (owing to its robustness against outliers) of the points in this subset: Furthermore, the fast point feature histogram (FPFH), proposed by Rusu et al. (2009), has been included in the feature set. It has been proven to be discriminatory for object detection in point clouds. FPFH is a 3D feature descriptor for point classification based on the surflets in the local neighborhood . Each combination of two surflets (surflet pair) in  gives a set of three angular relations (e.g., the azimuthal angle between the two surflets). The angular variations are discretized and graphed as a histogram for each relation. An overview over the features used is given in Table 1.
In order to reflect the impact of the features used in the classification task, an excerpt from the point cloud is visualized and segmented according to single chosen features in Figure 3. The relevance and contribution of the individual features with regard to the identification of roofs and their outlines can be shown.
In order to train a classification model, a set of labeled instances is needed. Therefore a set of approximately 80,000 points (including ≈ 21,000 building points) has been labeled manually as "building", "vegetation", or "ground". Given the labels and corresponding features (see Table 1) for the training data set, a (soft margin) SVM (Schölkopf & Smola, 2002), which was also used in Zhou and Neumann (2008) and Henn et al. (2013), has been learned. We used both the radial basis function (RBF) kernel and the linear kernel as kernel functions. After optimization of the kernel parameters, the accuracy of the model was derived using 10-fold cross-validation. For the linear kernel an accuracy of 82.83% and for the RBF kernel of 87.56% was achieved.
Despite its lesser accuracy the linear SVM model was chosen for the classification of the point cloud: first, the classification process is much faster using a linear kernel function; and second, for the generation of the gabled roof model hypotheses only a rough estimate of the class is necessary. building candidates nevertheless include single scattered points, which are regarded as "outliers". In order to eliminate the latter and to smoothen the labeling, we applied local outlier probabilities (LoOP) outlier detection F I G U R E 3 Excerpt from the point cloud segmented according to a subset of the features from Table 1. The relevance and contribution of the individual features to the identification of roofs and their outlines can be shown. Each figure separately illustrates one feature (Kriegel, Kröger, Schubert, & Zimek, 2009) followed by DBSCAN clustering (Ester, Kriegel, Sander, & Xu, 1996).
The advantage of LoOP is the combination of local, density-based outlier scoring with a probabilistic, statistically oriented approach. Therefore an upper limit for the outlier probability (e.g., 95%) can be set as a standardized measure. The points marked as "outliers" were excluded from the DBSCAN clustering. DBSCAN has two significant advantages over standard clustering algorithms, such as k-means: first, the number of clusters (here, the number of single detached buildings and connected building complexes) need not be known a priori; and second, there is no assumption on the shape of the clusters required. For this reason, DBSCAN can identify L-or U-shaped clusters which are beyond the scope of k-means. Its parameters eps and minPts define the radius of a region search for a query point and the minimum number of points which must belong to a single cluster.

| Approximation and regularization of building model outlines
In the next step, the polygonal outline of each point cluster was determined. The same method was applied for the calculation of the boundary of non-rectangular roof models (e.g., two intersecting gabled roofs) and dormers determined by We used a method based on a combination of the methods of Sampath and Shan (2007) and Lee, Han, Byun, and Kim (2011).
In a first step, the rough building or model outline was estimated. Following Edelsbrunner, Kirkpatrick, and Seidel (1983), α-shapes were used to determine an approximation of the boundary. Irregular and rough outlines, however, are quite untypical of building models. In most cases the outlines are linear and the intersection angles of two adjacent bounding segments are about 90 • . Consequently, the outlines had to be regularized, taking parallelisms and orthogonalities into account.
Analogous to Sampath and Shan (2007), segments of the α-shape were grouped based on the angle difference of subsequent segments. For each group, initial values were computed which were needed in a subsequent estimation of final boundaries. Therefore the dominant direction of the segments was needed which was derived by the method of Lee et al. (2011). Given the dominant directions, the initial values, and the points of the grouped segments, a Gauss-Helmert model estimation with angle constraints was performed which provided the estimated parameters of the bounding lines. The intersection of these lines led to the regularized building footprint.

| Detection of roofs and roof structures
This section describes active sampling for the detection of roofs and their parts from a LiDAR point cloud. This approach uses regularized approximations of predicted building outlines obtained in a previous step as described in Section 3.1. For the sake of clarity and conciseness, our sampling based detection method will be demonstrated for gabled roofs.

| Gabled roofs
Sampling is based on surflets  ⊂ ℝ 3 × ℝ 3 . A surflet s = (pt, n) ∈  is a combination of a 3D point pt and its normal vector n derived from the plane approximation of its k-neighborhood.
From our experience, k = 5 is a good choice for specifying this neighborhood. In some cases surflets derived from broader homogeneous regions will be preferred. Normal vectors, which in general are unique up to sign, are The best supporting local plane and its MSE are derived as usual. The curvature κ is defined as the ratio between the smallest eigenvalue and the sum of the three eigenvalues (see Equation 1).
Based on this information, a surflet s i for pt i with s i = (pt i , n i ) is generated. In the next step the curvature and MSE are used to define a filter characterizing the fitness of single surflets with respect to roof samples. Based on the smoothness criterion, the surflets are ranked in order to get a pre-filtered surflet set .
Note that small values of curvature and MSE suggest high fitness of single surflets with respect to roof sampling. Ranking is based on the derivation of a cumulative density function (cdf). Since neither curvature nor MSE can be assumed to be normally distributed, the calculation of the cdf is based on a kernel density estimation or (for the sake of efficiency) histograms. The weighted sum of the ranks corresponding to both features provides a first estimate of the fitness of a single surflet. The next step is to calculate the fitness of pairs  ′ of surflets, which are given as the cartesian product of the surflets which passed through the first filter. Ranking of pairs of surflets is done by the determination of the declinations Δϕ of the normal vectors of s k and s l for each pair sp = (s k , s l ) ∈  � (symmetry). Furthermore, the sum of the azimuths of normal vectors of s k and s l is computed (horizontality of ridges). The cdf represents a ranking criterion for the surflet pairs  ′ , resulting in a set of best-ranked surflet pairs  which will be used to verify the ridge prediction. We found that neither the size of the support set nor the score as applied in most RANSAC-based approaches (Henn et al., 2013;Tarsha-Kurdi et al., 2008) is a good measure for the appropriateness of a hypothesis. Instead, in the spirit of Karl Popper's "falsificationism" (Popper, 1959) and adapting an idea of Schmittwilken and Plümer (2010), we try to falsify inappropriate hypotheses based on their predictions. Particularly well suited for falsification are ridges (i.e., the points which should be observed in their neighborhood). For this purpose, the ridge line is computed for each pair sp i = (s k , s l ) ∈  by intersecting the corresponding planes. The occurrence of 3D points is predicted in the neighborhood of the ridge calculated, and then it is verified or falsified whether these points are inliers with regard to the ridge following a scoring strategy (Torr & Zisserman, 2000). If the ridge line can be verified with respect to the observations Pts, the corresponding roof hypotheses are added to a set of hypotheses . The whole algorithm for the surflet based active sampling for gabled roof hypothesis generation is summarized in pseudo-code in Algorithm 3.

F I G U R E 4
Roof model surflets of gabled roof. Roof properties are reflected by using surflet pairs (orange): (a) The horizontality of ridges is represented by the fact that the azimuths 1 and 2 of the normals ⃗ n 1 and ⃗ n 2 of the surflets differ by exactly π; and (b) Symmetric gabled roof: The declination φ must be identical Based on the set of hypotheses  acquired from the surflet based active sampling step, a model selection is performed in order to obtain a preliminary set of roof models. Algorithm 4 summarizes the procedure of model selection of 3D roofs using pseudo-code. For this purpose, the competitive hypotheses within  which support the same roof are grouped in a group g i ∈ . The model selection derives the best unbounded roof model m i for each group g i ∈  using model ranking based on prediction of ridges (number of observed points in the immediate neighborhood of the predicted ridge), general support (size of support set and score as in Torr and Zisserman (2000)) and average inlier density. These unbounded roof models are then collected in a set ' in order to be refined in a next step. Preliminary boundaries for each m i ' ∈ ' are determined in order to build a set of initial roof models.

| Hip roofs and dormers
The preliminary set of initial gabled roof models  0 acquired in the previous model selection step have to be refined in order to associate them to a most specific roof model type (in particular, hip or pavilion roof) if necessary. These steps are explained in detail in Algorithm 5. As stated in Henn et al. (2013), hip or pavilion roofs can be regarded as a special case of gabled roofs. Again, surflets are used. With regard to a given gable roof model, outliers which lie under the predicted roof model (negative residuals) may be explained by a hip roof/pavilion roof hypothesis, whereas outliers above the roof (positive residuals) may be explained by dormer hypotheses. In the case of a hip roof/pavilion roof, the azimuth of both surflets must be orthogonal to the azimuths of the roof.
Further, in order to classify dormers, a set of surflet pairs  ′ is constructed such that for each sp = (s k , s l ) ∈  ′ with s k ∈  i , s l ∈  i , the residuals of both surflets in sp have a positive sign. Three types of dormers can be discriminated. Surflet pairs from  ′ belonging to a gable dormer are characterized by an azimuth which is orthogonal to the azimuth of the corresponding roof model m i . Furthermore, a surflet from S i belongs to a flat roof dormer if the corresponding azimuth is undefined, whereas a surflet belongs to a shed dormer if the corresponding azimuth is the same as the azimuth of the roof model m i modulo π. Analogously to the roof model hypotheses, the competitive dormer hypotheses which support the same model are collected together in a group. Then a model selection is performed within each group in order to obtain the best unbounded model which will subsequently be refined by the determination of its boundary. This leads to a resulting set of dormer models . The latter are, in addition to the previous detected roof models, the basis for the reconstruction task which will be addressed in the next section. | 127 DEHBI Et al.

| Reconstruction of roofs and roof structures and experimental results
This section deals with the reconstruction of 3D roofs based on the roof models M with dormers  acquired from the refinement algorithm (Algorithm 5). Furthermore, the reconstruction results are presented and discussed.
The regularized approximations Outlines are determined following the method described in Section 3.1. Then the reconstruction method outputs a boundary representation Brep of the reconstructed 3D roof models. Figure 5 shows reconstructed 3D roof models based on our active sampling. Sophisticated roofs which ordinarily cause detection difficulties are reconstructed and highlighted on the margin of the figure. In particular, the otherwise challenging dormer identification and reconstruction turns out to be successful. Figure 6 illustrates a specific building where the roof is characterized by dormers of different types (shed and gabled). Although the single dormers represent a small fraction of the roof surface and the dormer clusters are small and close to each other, a separation using the active sampling is possible and leads to a good reconstruction.
The reconstruction is based on a test data set from Vaihingen in Germany provided by the DGPF which contains airborne laserscanner (ALS) data consisting of 10 ALS strips (Rottensteiner et al., 2014). The median point density is 6.7 points per square meter in overlapping strip areas. In regions covered by only one strip the mean point density is 4 points/m 2 . The number of buildings in the training phase amounted to about 140.
The problems of RANSAC-based sampling approaches described in the previous sections (cf. Figure 2) have been overcome due to both active sampling and the falsification of inappropriate hypotheses via ridge predictions.
For comparison, an outline of the resulting models derived with our approach is shown in Figure 7. The pseudo-code in Algorithm 6 summarizes the reconstruction method for the 3D roofs.
Roof reconstruction utilizes the boundaries of the roof model which are determined by clustering the inliers with the DBSCAN algorithm (Ester et al., 1996) and calculating its outline with alpha-shapes (Edelsbrunner et al., 1983). In the case of adjacent roofs, plane intersections are calculated in order to derive border lines. Furthermore, the topological connection between the neighboring roof models as well as partonomies such as the relation between dormers and their corresponding roofs are identified. In the next step, areas in the previously regularized approximation of predicted building outlines, OutLines, which are not covered by the roof models from  are detected and neighboring roof models must be extended, if necessary. Furthermore, buildings or roof structures which are too small (false positives) are eliminated. Finally, a regularization of the boundaries takes place in addition to an adjustment of neighboring ridges in order to produce the resulting reconstructed roof models Brep. F I G U R E 5 3D roof models as the result of an active sampling based reconstruction. Selected buildings with sophisticated roof models are highlighted. Data stem from the ISPRS benchmark test on urban object detection and reconstruction (Rottensteiner, Sohn, Gerke, & Wegner, 2014) (test area 3 from Vaihingen, Germany) F I G U R E 6 Reconstructed building based on active sampling with different dormer types (right). The roof and dormers have been reconstructed despite small dormer clusters which are close to each other (left) In order to evaluate our reconstruction results we submitted our Breps to the ISPRS test project on urban classification and 3D building reconstruction. Our results were compared to a reference data set for area 3 in Vaihingen (Germany). The completeness and correctness of our reconstructed models were evaluated on a perarea (Cm ar ∕Cr ar ) and on a per-object basis (Cm ob ∕Cr ob ). Furthermore, the completeness and correctness of objects larger than 55 m 2 (Cm 55 /Cr 55 ) were calculated. The thematic accuracy was determined as described in Rutzinger, Rottensteiner, and Pfeifer (2009). The geometrical accuracy was evaluated using the root mean square (RMS) error of distances from our building outlines to the reference building outlines. While Cr ar amounted to 93.9%, Cr ob was found to be 97.9% (cf. Table 2). The per-area completeness Cm ar was 67.5% and amounted to 85.5% for objects larger than 55 m 2 . This is attributed to the fact that we gave attention to structured roofs and did not perform a reconstruction of flat roofs. Table 3 gives detailed evaluation on a per-roof plane basis. The RMS error came to 1.1 m.
The computation time is dominated by the classification part using an SVM. In general, the learning phase of the classifier, comprising the optimization of the SVM parameters (c and γ) can take several minutes for several building roof models. However, this is only needed once. Application of the classifier is linear, and can be done in F I G U R E 7 Roof models estimated with our approach based on an active sampling. The problems occurred by RANSAC (cf. Figure 2) have been overcome TA B L E 2 Overview of the evaluation of our method using the ISPRS benchmark (Rottensteiner et al., 2014)  a fraction of a second for several roofs. For the implementation of DBSCAN, the dominating factor is the ɛ-neighborhood graph which is computed using a sweep-line algorithm in O(n log n+kn) without an index structure, where n is the number of the input points and k is the average number of points contained in a square-shaped neighborhood of size 2ɛ × 2ɛ centered around the points. Using a data structure (e.g., kd-tree) that supports queries for all points in the ɛ-neighborhood, the running time amounts to O(T setup (n) + nT query (n)). Likewise, the setup of the data structure is performed only once. All in all, in practice, the application of the learned model on a specific building is performed in a fraction of second. The experiments were performed on an Intel(R) Core(TM) i7-3770K processor.
The machine is clocked at 3.4 GHz and has 16 GB RAM. With regard to dormers which occupy a large portion of the roof, our method scores well with other approaches.

| CON CLUS I ON S AND OUTLOOK
The contribution of this article is a new method for the automatic derivation of LoD2 building models with larger roof structures from LiDAR data. Knowledge about roofs and roof structures such as dormers is used to select promising samples ("active sampling") for roofs and roof structures from the LiDAR point cloud and to reject 'bad' samples from this set as early as possible. The concepts of surflets and surflet pairs has been proven to be appropriate to perform this task. Surflets are ranked according to their smoothness and planarity. Pairs derived from the best-ranked surflets which do not satisfy roof properties (such as symmetry, horizontal ridges, or ridges of dormers being orthogonal to the ridge of the roof) are rejected as early as possible and very efficiently. From the remaining pairs, roof hypotheses are generated and selected by applying state-of-the-art model selection methods. Hence, our method implements a hybrid approach combining data-and model-driven approaches, reaping the benefits of both approaches and avoiding the disadvantages. In particular, we have shown that problems which occur when applying the RANSAC algorithm are solved. Currently our approach covers (symmetrical) gabled roofs, hip and pavilion roofs, as well as flat roofs, gabled and shed dormers.
In contrast to many other methods, our approach does not depend on building footprints as input. In a preprocessing step, points of the LiDAR point cloud are classified into building and other points using machine learning methods. The building points are clustered and a set of footprints which fulfills orthogonality and parallelism constraints is derived. This method can also be used for change detection tasks, in particular to detect buildings which are not present in cadastral data sets. The incorporation of cadastral footprints or even from volunteered geographic information data, however, will be of great benefit to improve reconstruction accuracy. A further integration of land use data and elevation models as background knowledge could also enhance the scalability of the proposed method.
The method is presented in this article in detail by giving the pseudo-code algorithms for all steps (Algorithms 1-5). The output of the algorithms is a 3D city model (detail level LoD2) in CityGML, a standard for 3D city models which is used internationally for data modeling and data exchange.
The next steps will be the extension of the catalog of roof types and roof structure types (e.g., chimneys or other dormer types) and the extension of the method to cope with 3D points derived from photogrammetry, in particular to points derived from semi-global matching methods. The geometrical characteristics of these points are different than those from LiDAR points. Furthermore, additional color information must be integrated. Hence, the methods must be adapted accordingly.
Further, intensity and RGB values can be investigated in order speed up the classification and hence the reconstruction process. The investigation of transferability of the approach can be studied for photogrammetric point clouds.