ROOF RECONSTRUCTION FROM POINT CLOUDS USING IMPORTANCE SAMPLING

: We propose a novel fully automatic technique for roof ﬁtting in 3D point clouds based on sequential importance sampling (SIS). Our approach makes no assumption of the nature (sparse, dense) or origin (LIDAR, image matching) of the point clouds and further distinguishes, automatically, between different basic roof types based on model selection. The algorithm comprises an inherent data parallelism, the lack of which has been a major drawback of most Monte Carlo schemes. A further speedup is achieved by applying a coarse to ﬁne search within all probable roof conﬁgurations in the sample space of roofs. The robustness and effectiveness of our roof reconstruction algorithm is illustrated for point clouds of varying nature.


INTRODUCTION
The generation of accurate 3D city models from point clouds is a very active research topic in computer graphics, photogrammetry and computer vision, (Wahl et al., 2008;Huang et al., 2013;Lafarge et al., 2013).The main focus of these related contributions has been to achieve (a) a significant amount of data reduction, (b) the automation of object recognition and reconstruction from unstructured 3D point clouds, and (c) the integration of semantics into building and 3D city models.However, most of the efforts in this domain are based on the analysis of terrestrial and airborne LIDAR data sets.The two major reasons for this are (1) the maturity of the LIDAR technology and (2) the resulting high accuracy and very dense 3D data.
However, due to the general availability of digital cameras, there is an immense interest in the generation of 3D point clouds from multiple images (Snavely et al., 2006;Frahm et al., 2010;Bartelsen et al., 2012).Though this data acquisition pipeline is being improved due to the developments of new algorithms for dense matching such as (Strecha et al., 2008;Hirschmueller, 2008;Kuhn et al., 2013), the acquired 3D points are still inferior to data acquired from LIDAR as far as accuracy is concerned.However, the 3D data are good enough to reliably recognize objects and their parts.Furthermore, these algorithms also benefit from the inherent flexibility using a standard consumer camera.We see a trend in an increasing quality of the 3D point clouds from matching and thus a necessity for devising methods for roof reconstruction from point clouds, which are robust with respect to the different data acquisition process.The goal of our work is to find, classify and fit basic roof types above a given quadrilateral building footprint.Thus, we aim at a level-of-detail (LOD) 2 model, consisting of the basic 3D shape such as a gable or a mansard roof.This will act as an intermediate state in a hierarchical modeling of buildings.A wide variety of applications in robotics, automotive navigation, tourism or city planning can benefit from such an automated building roof reconstruction scheme.We illustrate our algorithm with experimental results for data from matching images.We employ terrestrial images combined with images obtained from unmanned aircraft system collected in Southern Germany.Additionally, LI-DAR data sets and synthetic data sets generated from sampling a CAD model with additive noise are used in our experiments.The paper is structured as follows.Section 2 discusses existing approaches for reconstructing roofs in 3D point clouds.Beginning with an overview, we present our algorithm in Section 3.This is followed by the discussion of experimental results in Section 4. In Section 5, we conclude and give directions for future work.

RELATED WORK
There is a significant amount of work on the reconstruction of building roofs from point clouds.Most often, this is done based on domain knowledge, e.g., using shape priors of the objects of interest (Huang et al., 2013).Sampath and Shan (2010) take curvature as a feature within a fuzzy k-Means framework to cluster points defining roof segments in aerial LIDAR data.A similar region growing approach (Rottensteiner, 2010) detects planar surfaces by combining information from multiple images and point clouds.These surfaces can further be classified into various roof types.A Markov Chain Monte Carlo (MCMC) approach to roof reconstruction is presented in (Huang et al., 2013).Here, samples of possible roof constellations are proposed from within a predefined library of planar shapes.Proposals may or may not be accepted within a Reversible Jump MCMC (RJMCMC) setting depending on their computed likelihoods.However, the generative nature of the algorithm makes it very sensitive to noise and the RJMCMC sampler means a bad scalability of the algorithm.In (Taillandier and Deriche, 2004), a bayesian approach is used to fit polyhedral models from aerial images.The resulting planar patches are used to construct graphs from which buildings are reconstructed.Poullis and You (2009) presented a generic parametrization for the reconstruction of basic roof types.The likelihood for a specific roof type is assumed to underlie a Gaussian mixture model (GMM) whose parameters are inferred.However, the input data is assumed to be in perfect raster form.This is achieved by preprocessing including local energy minimization which can be highly error prone and is sensitive to noise.Moreover, the GMM assumption for a roof type implies perfect symmetry which will often not hold in real world.Recently, Henn et al. (2013)   that can cope with very sparse data sets.It assumes an underlying rectangular footprint.M-estimator random consensus (MSAC) proposed by Torr and Zisserman (2000) is used as a robust estimator to fit basic roof models.To differentiate between competing roof types, a support vector machine was used for classification.However, the approach is not suited for point clouds from image matching which are often highly irregular distributed and more complex roofs (mansard, gambrel etc.).Additionally, the assumption of a perfect rectangular footprint may be too strong a constraint.Besides the above mentioned deficits the majority of the mentioned algorithms are tailored for LIDAR data sets.
Our algorithm contains an extension of the parametrization of (Poullis and You, 2009), combined with a likelihood model similar to (Henn et al., 2013) within a Monte Carlo (MC) setting similar to (Huang et al., 2013), to optimally explore the search space for roof types above a detected general quadrilateral footprint.

Overview of the Algorithm
The input to our system is an unstructured 3D point cloud, D, containing a building as viewed from outside and above.The output is a set of connected polygonal surfaces above a general quadrilateral footprint with a defined roof type, e.g., mansard, gable or pinnacle.We suppose the 3D region above the footprint defines the most plausible search space for a roof.Appropriately, we transform D to a coordinate frame where the direction vector (0, 0, 1) points upwards.Starting at the highest point above the footprint and descending in discrete steps, we estimate probability distributions over possible roof configurations at several height levels.Making no assumptions about the form of these probability distributions over roofs at each level, our objective is to find the configuration that "best" explains D.
In Fig. 1, we show a gable roof illustration of this probabilistic roof fitting scheme using four roof configurations (coloured black, green, blue and yellow) and three height levels (l0, l1, l2).First, we build initial roof configurations (level l0) and estimate their likelihoods.Then, iteratively a transition is made to the next level and the roof configurations likelihood is updated.At each level, the configuration with the highest likelihood (marked bold in Fig. 1) explains our data better than all other configurations present.At height level l2, configuration θ1 (black in Fig. 1) provides a better fit than all other configurations.The problem is defined in a Bayesian setting and a sequential Monte Carlo (SMC) approach is used to search for a solution.

Notation
Random points p j are considered in the subspace Φ ∈ R 3 of the 3D Euclidean space.An ordered sequence θ ≡ (p 1 • • • p q ) ∈ Φ q uniquely defines a 3D surface in some standard way, e.g., in our experiments the p's are the vertices of a 3D polygonal chain.Thus, θ represents a single roof configuration.The symbol l represents a discrete height level index.P l defines a plane with normal vector (0, 0, 1), passing through height level l.The symbol ∼ means "distributed according to" as in x ∼ N (m, σ 2 ) defining the Gaussian random variable x with mean m and standard deviation σ and y ∼ U(a, b) refers to a uniform random variable y with support [a, b].Draws are collected, unless otherwise stated, in an independent and identically distributed (i.i.d.) way.The same notation will be used for a random variable and its realization to reduce notational clutter.
Being in an MC setting, we use a set of discrete roof configurations to approximate a distribution over roofs.For example, at level l, {θ l } M i=1 are a finite set of M such 3D polygonal chains used as a discrete representation of p(θ l |D), the distribution over polygonal chains at level l.The likelihood of a configuration, L(θ (i) l ), denotes how the data is explained by that polygonal chain.The normalized likelihood defines the weight of a configuration, w (i) l .This is used as an approximation of the posterior probability, p(θ (i) l |D).At each height level, the configuration with the highest weight is reckon on as the maximum a-posteriori probability (MAP) configuration, denoted θMAP (bold lines in Fig. 1).

Roof Fitting using a Sequential Monte Carlo Method
Having the complete set of data, D at our disposal, we define the sequence as height levels above the foot print.Starting at the highest point in D above the footprint and descending in discrete ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W3, 2013 CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey levels along the vertical direction, we employ sequential importance sampling (SIS) to consistently keep track of roof configurations at these height levels.This iteratively approximates the distribution over roof configurations.The goal is to move the starting configurations based on a forward transition distribution, p(θ l |θ l−1 ), which retains expected properties of the 3D polygonal chains to be fitted and on an observation model, p(D|θ l ), that provides evidence about whether a measured 3D point in D is in the close vicinity of the "true" 3D polygonal chain.The complete algorithm summarized in the flow chart given in Fig. 2 can be described as follows: 1.In initialization, we generate M configurations for every considered roof type (mansard, gable, flat, etc.) 2. In forward prediction, these M configurations per roof type are propagated using the transition probability distribution, p(θ l |θ l−1 ).This results in a change of configurations from one level to another 3. Within an update stage, likelihoods of the configurations, L θ , are estimated and normalized, thus giving a discrete approximation of the posterior probability, p(θ l |D) 4. Steps 2 and 3 are iterated until a minimum level well above the ground plane (denoted by lmax index) is reached 5.At the minimum level, we are confronted with three problems: (1) Which of the considered roof types explains our data better than the others?(2) Amongst the levels, which posterior distribution, p(θ l |D), should be chosen?(3) Considering all the M configurations accounting for the chosen posterior distribution, which configuration explains our data "best"?
We divide the above algorithm into four major parts: Initialization, forward prediction, weight estimation (roof configuration likelihood estimation) and model selection.Details of these major parts of our algorithm are described in the following.

Initialization
The generation of initial configurations means to choose an appropriate parametrization for the various basic roof types and corresponding prior distributions for them.We could extend the generic mansard roof parametrization of (Poullis and You, 2009) to include all symmetric and asymmetric variants of these basic roof types.However, this would be too computational intensive for a scene with few mansard roofed buildings.Therefore, we have devised more distinctive parametrization for different types of roofs.We use non-informative priors for the parameters defining the building heights, p 0 , p 1 , p 2 , p 3 in Fig. 3 and Fig. 4.This prior compensates for uncertainty in the estimates of heights of the buildings walls.Therefore, for the points p 0 , p 1 , p 2 , p 3 , we make four i.i.d.draws given by where k ∈ {0, 1, 2, 3} stands for the respective points p k , h d is the height tolerance as depicted in Figs. 3 and 4 and i is the index of the i-th of M configurations.Since the direction vector (0, 0, 1) has been defined as the vertical direction of the scene, equation ( 1) means a 1-dimensional search along the zaxis.These four draws representing a prior over the height of the walls are fundamental for all roof types.In addition to equation (1), we use domain knowledge and derive informative priors over the ridges, peaks, hips etc. to give a complete prior for a certain roof type.Gable Roof Prior Fig. 3 illustrates the parametrization of the gable roof.A prior over the ridge combined with equation (1) defines our informative prior model for the gable roof type.A realization of the ridge is defined by the line segment between points p ra and p rb projected on plane P l , where p (i) ra ∼N p ma, ht and p with p ma = 0.5 • (p 0 + p 1 ) and p mb = 0.5 • (p 2 + p 3 ), ht defines the tolerance in the ridge and i is the index of the i-th configuration.These are 1-dimensional draws along the direction of the walls.Since Equations ( 1) and (2) are i.i.d.samples, this prior model for the gable roof captures the complete family of gable roofs (symmetric and asymmetric).It is constructed twice such that in each case, the two sloping roof planes are oriented according to wall pairs.This ensures the right orientation of the roof planes to the pair of opposite walls.
Pinnacle Roof Prior Similarly, for the pinnacle roof shown in This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.75 tion of the peak is defined by with p mid = 0.25 • (p 0 + p 1 + p 2 + p 3 ) the centroid of the four points defining the footprint, hc is the tolerance of the intersection point of all four hips and i is the configuration index.p peak is projected on P l to get the true peak.
Mansard Roof Similarly, Fig. 4 shows the parametrization of the mansard roof.
where hn and hm are tolerance values defining the lengths of the ridges and again, i is the index representing the i-th configuration.
The true ridges are a projection of the line segments h mwidth and h (i) mwidth on P l .If the outcome of Equation ( 4) lies within a circle with diameter hc (cf.Fig. 4d),i.e., a highly probable pinnacle roof candidate, it is rejected.The final valid search space for mansard roof is marked gray in Fig. 4d.
Other Roof Types The most basic roof types are the flat and shed roofs.Their parametrization is obtained by varying the four points defining the building height, i.e., equation (1).Figs.4e  and 4f shows our template for the general family of hipped roofs, symmetric, asymmetric and (half)-hipped.Obviously, from this template the parametrization is obtained by a combination of gable and mansard roofs.To increase the effective number of samples and avoid "inverse" and thus invalid roof types, a constraint is added to all parameterizations that guarantees the ridge being higher than the building walls.

Forward Prediction
In Fig. 2, we describe the movement of roof configurations from one height level to another as sampling from a forward transition distribution, p(θ l |θ l−1 ) given by the following simple first-order model Having present the initialization and the forward transition probability, we now proceed to present how the configuration likelihoods of the 3D polygonal chains are computed.

Roof Configuration Likelihood Computation
We use MSAC (Torr and Zisserman, 2000) to determine if a point in D is within the close vicinity of a 3D polygonal chain defined by a configuration, θ.MSAC underlies a likelihood function proportional to a truncated Gaussian bounded by the RANSAC threshold, T 2 .Henn et al. (2013) used a similar likelihood model on aerial LIDAR data.Yet we optimally explore the configuration space in a Bayesian setting.Thus at level l, for a certain configuration, θ (i) l , assuming measurements being conditionally independent, the likelihood function is given by with and e 2 j the shortest Euclidean distance from point p j to the surface of the 3D polygonal chain defined by the configuration θ (i) l .The summation index j runs over all the points in D.

Model Selection
Equation ( 6) defines the likelihood for a single configuration of a specific roof type.We define the configuration weight w, as the normalized likelihood, given by Therefore, the posterior probability at level l for a roof type is approximated by where δ() defines the Dirac function.So far, three major questions are unanswered: (1) the roof type to choose (2) the height level to select and (3) within the chosen height level and roof type, the configuration that fits D "best".For the first two questions, we compare weights of the MAP configurations at each level for all the competing roof types and select the roof type whose θMAP has the highest weight.The height level within which the selected, θMAP resides defines the selected height level and is used to answer the third question.For this, we adopt the M-opened modelling perspective (Smith and Bernardo, 2000) i.e., we believe the "best" configuration may not be within the M configurations of the selected height level for the chosen roof type, however, through Bayesian model averaging, we compute the minimum mean square error estimator (MMSE) configuration as the inferred "best" configuration.For the chosen roof type within the selected height level, this is given by The MMSE estimator in this setting is analogous to taking a weighted least square fit with the weights proportional to the estimated likelihoods as defined in section 3.3.3.

Implementation Details
The major computational bottleneck is Equation ( 6), since the likelihood has to be estimated M times per roof type for every level.To speed up this compuation, we use a coarse to fine search.First D is downsampled using a voxel grid with leaf size 0.4m × 0.4m × 0.4m defining the input point density.A search is conducted with a large enough inter-height level value to capture the height level where the roof most probably lies.Around the selected level, we apply a finer inter-level value for the final search.We also use a fast version of the polygon hit test from (Schneider and Eberly, 2002) for the various polygon segments of a 3D polygonal chain defining a configuration.Furthermore, since the heights of walls of a building are correlated, a pair-wise sampling of the building height priors was prefered to four i.i.d.draws of Equation (1).In our formulation of SIS, a strong data parallelism is implied due to the 3D space partitioning in discrete height levels.Using alternative sampling schemes such as Gibbs, Metropolis-Hastings or sequential importance resampling samplers, this proper decoupling of samples from one level to another is impossible.

Parameters and Data Sets
For all performed experiments, we assumed the input data is on a metric scale and the following parameters were used: ISPRS Annals of the Photogrammetry, Remote Sensing and Information Sciences, Volume II-3/W3, 2013 CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey   For the mansard roof type, hn and hm were both chosen to be 0.5 meters less than length and width of the footprint.Furthermore, we discretized the search space above the footprint in 10 and 5 equal steps for the coarse and the fine search respectively.
Example outputs from data set obtained from image matching is shown in Figs. 5 and 6.For Fig. 5, a combination of images taken from the ground and a UAV were used for matching and the mansard roof type best explains the data set.Meanwhile Fig. 6, shows a highly unequal point distribution in the data set with very few points on the roof as a result of the matching of images taken from the ground alone.We were still able to get the pinnacle roof type as the best fitted model.Since synthetic data sets are more or less perfect, for all experiments, all points were perturbed using 11) where uj ∼ N (0, 1m) for evaluation purpose.

Results
We tested our algorithm on a data sets containing 152 quadrilateral footprints of 3D points of varying density and origin.The results is summarized in the Table 1.All wrongly fitted roofs resulted in a flat roof type.This was due to the very low inclination angle of the roof planes.However, by tuning the parameters one can certainly get rid of these inaccuracies in the fits.

Clutter
The data sets on Figs. 5 and 6 contains dormers and a chimney respectively.However, these are fairly small roof superstructures and thus their effects on a roof fitting algorithm may be considered negligible compared to the highly unequal point distribution in the data set.Meanwhile, we demonstrate the robustness of our algorithm on a data set with substantial clutter through roof superstructures of reasonable sizes.In Fig. 7, results from a LIDAR data set containing a building with T-shaped footprint.
We decompose the footprint into two separate non-overlapping quadrilaterals and apply our roof fitting algorithm 2 on each decomposed footprint.In both cases, gable roof type provided the best fit.The robustness against substantial clutter can be substantiated in twofold.First, our distinctive designed priors accounts for good discission making amongst competing models.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W3, 2013 CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey Secondly, the choice of the likelihood and the proper exploration of the search space within a Bayesian setting provides the optimal fit in MMSE sense.

CONCLUSION
We presented a sampling based method for automatically fitting and classifying building roof models from 3D point clouds.Assuming a general quadrilateral footprint, the method starts with the construction of parametric models of different roof types.These are represented by an ordered sequence of 3D polygonal surfaces.A likelihood function is chosen that considers the closeness of the 3D Data to the 3D polygonal chains.These likelihoods are used within a Bayesian setting to build probability distributions for various roof types in several possible discrete height levels above the footprint.Using model averaging, a model selection scheme is chosen to get the "best" model in MMSE sense.We achieve a computational speedup by applying a coarse to fine search.The method yields geometrically accurate parametric building roof models and is robust concerning the source and nature of the input data.For the generic case of modeling 3D buildings with a quadrilateral footprint, our approach robustly produces reliable results from which further refinements can be made, progressively increasing the level of detail towards detailed model (LOD3).A parametric representation is chosen which directly encodes to state-of-the art visualization and representation environments such as CityGML.Practical problems arise particularly due to extreme occlusions and wrong priors.To overcome this problem, it would be desirable to extend the likelihood function.

Figure 1 :
Figure1: At each height level, the configuration with the highest likelihood is marked in bold.The black arrows describe the transition of the configurations from one level to the next and the gray dots are data points.Different colors correspond to different configurations.Please note that the space between height levels is widely exaggerated for a better comprehension.

Figure 3 :
Figure 3: Top row: gable roof model including variance in building heights (h d ) and ridge position (ht).Bottom row: Pinnacle roof model showing hip intersection point tolerance circle with diameter hc.

Figure 4 :
Figure 4: Top row (a), (b) and (c): Illustration of mansard roof parametrization.Bottom row: In (d) the gray region illustrates the actual search space for the mansard roof excluding pinnacle roof candidates.(e) and (f) illustrate the symmetric and asymmetric (half)-hipped roof parametrization respectively.
Fig. 3, a combination of Equation (1) and a prior over the peak location, p peak , defines the informative prior model.A realiza-ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W3, 2013 CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey

Figure 5 :
Figure 5: MMSE fit for flat, pinnacle, gable and mansard roof types on a mansard roofed building data set obtained from image matching.The yellow points are the inliers to the "best" 3D polygonal chain.As expected, the mansard roof model explains the data best by providing the best fit.

Figure 6 :
Figure 6: Top row: Image and MMSE fit using a gable roof type.The model is viewed from inside the building.Bottom row: MMSE fit using a pinnacle roof model.The model is viewed from outside (left) and inside (right) the building.

Figure 7 :
Figure 7: Data from LIDAR with a T-shaped building footprint decomposed into two separate non-overlapping quadrilaterals.The yellow points are the inliers of the MMSE fitted roof type over each footprint.Best fit in spite of substantial clutter

Table 1 :
Results of our experiments conducted on different data sets.The number of correct versus wrong fits given the footprints.