ENHANCEMENT OF GENERIC BUILDING MODELS BY RECOGNITION AND ENFORCEMENT OF GEOMETRIC CONSTRAINTS

: Many buildings in 3D city models can be represented by generic models, e.g. boundary representations or polyhedrons, without expressing building-speciﬁc knowledge explicitly. Without additional constraints, the bounding faces of these building reconstructions do not feature expected structures such as orthogonality or parallelism. The recognition and enforcement of man-made structures within model instances is one way to enhance 3D city models. Since the reconstructions are derived from uncertain and imprecise data, crisp relations such as orthogonality or parallelism are rarely satisﬁed exactly. Furthermore, the uncertainty of geometric entities is usually not speciﬁed in 3D city models. Therefore, we propose a point sampling which simulates the initial point cloud acquisition by air-borne laser scanning and provides estimates for the uncertainties. We present a complete workﬂow for recognition and enforcement of man-made structures in a given boundary representation. The recognition is performed by hypothesis testing and the enforcement of the detected constraints by a global adjustment of all bounding faces. Since the adjustment changes not only the geometry but also the topology of faces, we obtain improved building models which feature regular structures and a potentially reduced complexity. The feasibility and the usability of the approach are demonstrated with a real data set.


INTRODUCTION 1.1 Motivation
Scenes of arbitrary complexity can be described by generic models, e.g.boundary representations or polyhedrons respectively, in context of modeling man-made objects.Unfortunately, representations obtained by purely data-driven approaches cannot distinguish between buildings and non-building objects since they contain no building-specific knowledge.Actually, data-driven approaches are unable to recognize man-made objects at all but feature the flexibility to capture objects of arbitrary shape.Specific models, e.g.parametric models, can be directly related to buildings.Still, being too specific has the disadvantage of restricting the number of representable buildings.In general we have the classical dilemma of choosing a generic, but too unspecific, or a specific, but too restrictive model (Heuel and Kolbe, 2001).
Several concepts exist to introduce building shape knowledge without being too restrictive.In (Xiong et al., 2015) for instance, the topological graphs of roof areas are analyzed to combine and instantiate predefined building primitives.Another way is to recognize man-made structures by a hypothesis generationand-testing strategy.Geometric and topological relations such as orthogonality or parallelism are found by statistical hypothesis testing and then established by an adjustment (Meidow, 2014).
To perform the adjustment, a set of consistent and non-redundant constraints is required.Here, consistency refers to non-contradicting constraints and redundancy to independent constraints.Thus selection strategies are required to obtain such minimal sets reliably.Conceivable approaches are automatic theorem proving (Loch-Dehbi andPlümer, 2009, Meidow andHammer, 2016) or decisions drawn on numerical criteria within a greedy algorithm (Pohl et al., 2013).* Corresponding author The pyramidal broach roof has been reconstructed as a hipped roof with a very short ridge line.The adjustment process enforces the concurrence constraint and therefore changes the topology of the roof to a pyramidal roof.

Contribution
Given the generic model of a building in a boundary representation, we recognize geometric relations such as orthogonality or parallelism automatically by hypothesis testing and enforce the recognized constraints by an efficient adjustment procedure.As a result, we obtain enhanced building models with regular structures and potentially simplified topology.
An example is shown in Figure 1: A pyramidal broach roof is initially represented by a hipped roof at which the ridge line is very short.After recognizing that four planes share a common point with high probability, an adjustment enforcing this constraint has been carried out.This leads to a pyramid as a more appropriate roof description.
For the execution of the hypothesis testing and the subsequent adjustment we have to specify the uncertainty of the model geometry.We assume that no information about the geometric uncertainty comes along with the boundary representation.Therefore, we perform a point sampling for each bounding face to obtain a For the recognition of constraints only the data-independent significance level, i.e. the probability of correctly not rejecting a null hypothesis, has to be specified.The usability and the feasibility of the approach are demonstrated by the analysis of a generic boundary representation for a building captured by airborne laser scanning.The reconstruction has been carried out according to the approach proposed in (Xiong et al., 2014), cf. Figure 2.

Related Work
We consider the approach proposed in (Xiong et al., 2014) as one of the most recent methods to derive polyhedral building models from airborne captured point clouds.Formulations for constraints as multivariate polynomials can be found in (Brenner, 2005).Constraints and statistical tests for geometric entities in homogeneous representation are provided by (Heuel, 2004, Förstner et al., 2000).In (Pohl et al., 2013), a greedy algorithm is used to select a set of independent and consistent constraints automatically.An adjustment model to solve constrained problems formulated in homogeneous coordinates efficiently can be found in (Meidow, 2014).In (Loch-Dehbi andPlümer, 2011, Loch-Dehbi andPlümer, 2009) independent constraints are found by automatic theorem proving using Wu's method (Wu, 1986).The feasibility is demonstrated but no real data sets have been analyzed.

THEORETICAL BACKGROUND
After the explanation of the point sampling process and the determination of adjacent point groups (Subsection 2.1) we explain the utilized parametrization of planes and the corresponding parameter estimation (Subsection 2.2) and formulate geometric constraints (Subsection 2.3).The recognition of man-made structures by hypothesis testing (Subsection 2.4) and the enforcement of the recognized constraints by an adjustment (Subsection 2.5) are then the core of the reasoning process.Within this approach, the deduction of a set of consistent and non-redundant constraints is essential (Subsection 2.6).

Point Sampling and Adjacencies
The uncertainties of the geometric entities are required for the complete reasoning and adjustment process.This refers to the statistical testing of geometric relations (Subsection 2.4) as well as to the adjustment of the planes (Subsection 2.5).Thus the uncertainties of bounding planes have to be specified by covariance matrices.Since we assume that the model does not come along with this information, we perform a point sampling for each face of the boundary representation.This sampling simulates the result of the laser scanning initially used for the acquisition of the given model.
Figure 3 shows an exemplary part of the point sampling.For each face, points on a grid with given spacing have been generated.Then noise has been added to the point coordinates.The simulated point cloud of a model face is used to estimate the values of the plane parameters (Subsection 2.2).Furthermore, the point clouds can easily be used to determine the adjacencies of bounding faces which are then used for hypothesis generation.Note that in this way relations can be established which are not contained in the boundary representation at hand.

Representation of Planes and Parameter Estimation
A convenient representation of a plane is its normal vector N with N = 1 and the signed distance D of the plane to the origin of the coordinate system.The corresponding homogeneous 4-vector is with its homogeneous part A h and its Euclidean part A0.Since ( 1) is an over-parametrization, we need to introduce additional constraints for the plane parameters to obtain unique results within the adjustment.We use the spherical normalized vector with corresponding covariance matrix Σ AA = J Σ AA J T according to error propagation, the covariance matrix ΣAA of the unconstrained plane parameters and the Jacobian For the adjustment we need the spherical normalization (2) as a constraint for the parameters which reads with the factor 1/2 for convenience.
The plane parameters can be estimated in a statistically optimal manner for 3D points which are mutually uncorrelated, i.e. each constraint polynomials m Table 1: Potential geometric relations between planes A, B, C and D expressed by multivariate polynomials and the number of independent equations (degrees of freedom m).
point Xi has covariance matrix σ 2 i I 3 with the 3 × 3 identity matrix I3 and standard deviation σ i for the coordinates of the point.The plane passes through the weighted centroid and its normal is the eigenvector belonging to the smallest eigenvalue of the empirical covariance or moment matrix M of the point group.
The covariance matrix of the plane parameters can be derived from the eigenvalues of M , too (Förstner and Wrobel, 2017).The 4×4 covariance matrix is full and singular, nevertheless the information is considered completely by the adjustment model presented in Subsection 2.5.

Formulation of Constraints
With the homogeneous representation of planes, the formulation of constraints is simple.Table 1 summarizes the considered relations between planes, see (Heuel, 2004) for details.The concurrence constraint states that a 3D point is incident to four planes.
Further specific constraints such as identical slopes for roof areas or horizontal roof ridges are conceivable.But often these constraints are implicitly already given, which limits their practical relevance, depending on the specific application.
For the recognition of constraints by hypothesis testing and for the enforcement of constraints by adjustment, we need the Jacobians of the constraints listed in Table 1.These are simple and can be provided analytically.For the concurrence constraint ∂ det(M )/∂vec(M ) = adj(M ) holds, where adj(M ) is the adjugate of the matrix M = [A, B, C, D] and vec

Recognition of Man-Made Structures
For each bounding face, a point cloud is sampled simulating the scanning results.Two point clouds A and B are considered to be adjacent if the Euclidean distance between two points-one from cloud A and one from cloud B-is shorter than a threshold T d , say T d = 0.5 m.Choosing a large value for T d increases the number of established neighborhoods.The cliques of the neighborhood graph can be determined by the Bron-Kerbosch algorithm (Bron andKerbosch, 1973, Cazals andKarande, 2008); they constitute the set of potential relations which have to be checked.

Hypothesis Checking
To check an assumed relation, we formulate a classical hypotheses test.We test the null hypothesis H0 that an assumed relation is valid against the alternative hypothesis Ha that the relation does not hold.Formally, we test (5) with the expectation values µ d for the distance measures d as provided by Table 1.
The corresponding test statistic Figure 4 illustrates the distribution of the test statistic ( 6) checking the parallelism of two planes obtained by simulation.We sampled two planes with the rather small number of n1 = 4 and n2 = 6 3D points, estimated the two sets of plane parameters independently, and computed the test statistic d , see Table 1.The simulation has been performed 10,000 times with uniformly distributed random coordinates X and Y in the interval [0, 1] and normal distributed Z-coordinates with σ Z = 0.01 and mean 0.
As one can see from the histogram, the empirical data fits well to the Fisher-distribution with m = 2 and n = 6 + 4 − 2 • 3 = 4 degrees of freedom.For larger values of σZ , the effect of linearization errors become visible.
For the test, we choose the data-independent significance number of α = 0.05, i.e. the probability of erroneously rejecting the null hypothesis in case it actually holds (type I error).Then the critical value F1−α;m,n of the corresponding distribution is computed which limits the regions of acceptance and rejection for the test statistic.The hypothesis will not be rejected if T < F1−α;m,n holds for the value of the test statistic (6).
Note that the power of a test, i.e., the probability of not rejecting a null hypothesis while the alternative hypothesis is true, depends also on the sample size.

Pre-Check
As one can see from equation ( 6), the value of the test statistic becomes small if the variance of the distance measure is large.In this case, the null hypothesis tends to be not rejected although the alternative hypothesis is correct (type II error).Therefore, we perform a pre-examination of the data quality by checking the acceptability of the achieved precision for the plane parameters (McGlone et al., 2004).To do so, we provide a reference covariance or criterion matrix CEE for a plane in canonical form E = [0, 0, 1, 0] T , e.g.
which specifies the required precision for the plane parameters, see Figure 5. Thus, we require a standard deviation of at most 3 • for the two principal normal directions N E = [0, 0, 1] T of a The planes in general orientation are moved to the reference position [0, 0, 1, 0] T by a 3D motion accompanied by error propagation.For this purpose, the rotation and translation can be extracted from the estimated covariance The estimate is acceptable if the largest eigenvalue of cf. (McGlone et al., 2004).

Enforcement of Constraints
For completeness, we provide a brief summary of the utilized adjustment model.For details please refer to (Meidow, 2014).The solution of the proposed adjustment is equivalent to the adjustment with constraints for observations only (Koch, 1999, Mc-Glone et al., 2004).In these models, the existence of a regular covariance matrix for the observations is assumed as its inverse serves as a weight matrix within the adjustment.For the treatment of singular covariance matrices, we specialize the general adjustment model introduced in (Meidow et al., 2009).
We compile the unconstrained, spherically normalized parameters of K given planes in a parameter vector x = [A T 1 , . . ., A T K ] T whose singular covariance matrix Σxx features block-diagonal shape.The adjusted or constrained parameters x are related to the unconstrained parameters x by additive unknown corrections which have to be estimated, thus x = x + holds.
The functional model subsumes the constraints to be enforced by the adjustment.The intrinsic constraints are the spherical normalizations (2) of the homogeneous vectors needed to obtain unique results.They are formulated as K restrictions k( x) = 0 on the parameters.The geometric constraints are the H restrictions h( x) = 0 on the plane parameters, e.g.orthogonality or parallelism, cf.Table 1.
The stochastic model reflects the uncertainty of the unconstrained parameters, described by an initial block-diagonal covariance matrix Σ (0) xx which results from the parameter estimation (Subsection 2.2).The covariance matrix Σxx is singular due to overparametrization.The explicit constraints for the observed planes, i.e., the normalization of the homogeneous entities, and the corresponding implicit constraints comprised in the covariance matrix Σxx must be consistent.
For the derivation of the normal equation system we linearize the constraints.The sought constrained parameters x are given by the unconstrained parameters adjusted by estimated corrections , or by approximate values x0 and estimated updates ∆x.Thus x = x + = x0 + ∆x holds and the linearization of the constraints reads with the Jacobians K and H and the auxiliary variables For the intrinsic constraints K = Diag(A1, A2, . . ., AK ) and KK T = I holds due to the constraints (A T k A k − 1)/2 = 0 for each entity indexed with k.The Jacobian H can be extracted from the equations listed in Table 1.
Minimizing the sum of weighted squared residuals yields the Lagrangian with the Lagrangian vectors λ and µ.Setting the partial derivatives of ( 13) to zero yields the necessary conditions for a minimum and the normal equation system reads Solving ( 15) constitutes a solution to the problem at hand but this implies the computation of the pseudo inverse for Σxx, a matrix whose size depends on the number of parameters, i.e., the number of planes.Rewriting (15) leads to the reduced system With ΣxxK T = O and KK T = I the relation holds due to the definition of pseudo inverses and the normal equation matrix can be inverted explicitly: Thus the substitution of the estimate = −ΣxxH T µ − K T k0 in (10) yields the Lagrangian vector µ = Σ −1 hh h0 − HK T k0 with the regular covariance matrix of the small deviations (12) due to variance propagation.Finally we obtain the estimate for the corrections.Note that the size of the matrix to be inverted depends only on the number of constraints and that the pseudo inverse of Σxx has not to be computed explicitly.Furthermore, no special treatment of entities not being affected by any extrinsic constraint is necessary.In this case the corresponding estimates for the corrections (20) will simply be zeros.
The estimates are x = x + and one has to iterate until convergence.Note that the covariance Σxx must be adjusted during the iterations since the linearization points for K change and ΣxxK T = O must be fulfilled.This can be achieved by error propagation for the spherical normalization (2).

Deduction of Required Constraints
For the adjustment presented in the previous subsection, a set of consistent, i.e. non-contradicting, and non-redundant constraints is mandatory.Redundant constraints will lead to singular covariance matrices ( 19) for the small deviations.Since we are dealing with imprecise and noisy observations, we have to face the possibility of non-rejected hypotheses which are contradictory, too.
Greedy algorithms are common approaches to select such sets of constraints-either by algebraic methods, e.g.(Loch-Dehbi andPlümer, 2011, Meidow andHammer, 2016), or by numerical methods, e.g.(Pohl et al., 2013).We adopt the approach proposed in (Pohl et al., 2013) to select sets of constraints automatically: After adding a constraint, we consider the estimated rank and the estimated condition number of the resulting covariance matrix (19).In case of a rank deficiency and/or a small condition number, we neglect the additional constraint.

EXPERIMENTS
After introducing the input data, i.e. the utilized point cloud and the derived generic building representation, we illustrate the point sampling, summarize the deduced constraints, and show the results of the adjustment process which enforces the man-made structures.

Input Data and Point Sampling
We start with the derivation of a generic boundary representation for buildings.Figure 6 shows a detail of an airborne laser scan which represents a rural scene with roofs of farmhouses as dominant man-made objects.The RANSAC-based shape detection (Schnabel et al., 2007) as provided by the 3D point cloud and mesh processing software CLOUDCOMPARE has been used to extract and group points which are likely to lie on planar facets of the buildings.
Next we derived a generic polyhedral building representation by the procedure proposed in (Xiong et al., 2014).Figure 2 shows the resulting boundary representation which is the input for the proposed reasoning process.Due to the specific reconstruction process, the walls of the buildings are vertical and the bottom area of the building is horizontal.Apart from that, the boundary representation of the roof is unconstrained.
For each facet of the boundary representation we carried out a point sampling by equally spacing points with added noise to the point coordinates (Figure 3).The point spacing of D = 0.10 m and the noise of σ = 0.03 m for point-plane distances simulate the assumed scanning process.Thresholds influence the degree of generalization for the reconstruction result.However, small facets cannot be avoided in general and a few rather uncertain facets can be created for the boundary representation.The consideration of the uncertainties of the planes with the pre-check proposed in Subsection 2.4.2 allows for this.

Results
Figure 7 illustrates all relations recognized by the hypothesisgeneration-and-testing strategy.The utilized constraints for the planes are identity, verticality, orthogonality, parallelism, and concurrence.The latter considers four planes meeting in a common 3D point-a situation often encountered in building polyhedrons.Relations with the bottom area of the building are omitted in Figure 7 for the sake of clarity.Obviously, many relations are redundant, however, it is not clear which can be omitted without loss of information.
Table 2 summarizes the number of recognized and deduced constraints itemized according to their type.The set of deduced constraints is minimal in the sense of consistency and redundancy but the actual number of constraints depends on the order of constraints as processed by the greedy algorithm.In general, the number of actually required constraints is considerably smaller than the number of recognized constraints.
Figure 8 shows the result of the adjustment process with the strictly enforced constraints listed in Table 2. Eye-catching is the change of topology for the building with the pyramidal broach roof: In the given boundary representation, the roof is modeled as hipped roof with a very short ridge line (Figure 2).In the reasoning process, the four planes concurring in a single point have been recognized and the adjustment changed the roof representation to a pyramid.2. The reasoning process led to a more adequate representation of the pyramidal broach roof.

CONCLUSIONS AND OUTLOOK
The feasibility and usability of the proposed approach has been demonstrated with a first real data set.The only input is a generic boundary representation of a building.Data-dependent thresholds are avoided by the exploitation of statistical methods which take the uncertainty of the bounding planes into account.
For the future, several improvements and investigations should be carried out: • For the simulation of point observations as a result of laser scanning, more sophisticated methods can be envisioned, which take sensor-specific characteristics into account and consider sensor trajectories and scene geometries.
• The formulation of hypotheses for geometric relations should be expanded to express cartographic rules for the generalization of boundary representations as a further application.
• The estimation of matrix ranks and condition numbers is inexact and becomes unrealizable for large-scale problems.Therefore, algebraic methods should be considered for the deduction of consistent and non-redundant constraints as they provide exact results.Some progress in this regard has been achieved, see (Meidow and Hammer, 2016).
• The set of deduced constraints depends on the order of the constraints within the greedy algorithm.Thus, general rules for the ordering of constraints are advantageous.
• The results depend on the applied point sampling and the adding of noise.This dependency should further be investigated.
Figure1: The pyramidal broach roof has been reconstructed as a hipped roof with a very short ridge line.The adjustment process enforces the concurrence constraint and therefore changes the topology of the roof to a pyramidal roof.

Figure 2 :
Figure 2: Detail of a point cloud obtained by airborne laser scanning and the derived generic polyhedral building reconstruction.point cloud which simulates the result of an airborne laser scanning.Estimates for the plane parameters are obtained from these point clouds which are then the basis for the subsequent reasoning process.

Figure 3 :
Figure 3: Result of the point sampling for each face of the boundary representation (bird's eye view) with sampling distance 20 cm and σ = 3 cm for each coordinate.The points of the bottom area have been omitted for clarity.

Figure 4 :
Figure 4: Distribution of the test statistic (6) for the parallelism of two planes.Two point sets with n1 = 4 and n2 = 6 points have been sampled 10,000 times.is Fisher-distributed with m and n degrees of freedom.Here, m denotes the degrees of freedom for the specific relation, e.g.m = 2 for the test of parallelism, and n is the redundancy of the parameter estimation, i.e. the number of equations used to estimate the plane parameters involved in the relation minus the number of unknown plane parameters.The covariance matrix Σ dd of the contradictions d is derived by error propagation.

Figure 5 :
Figure 5: Plane in canonical form E = [0, 0, 1, 0] T with normal vector N = [0, 0, 1] T .The confidence region is a hyperboloid of two sheets, specified by the covariance matrix (7).plane and a maximum standard deviation of 0.05 m for distance DE = 0 of the plane to the origin.

Figure 8 :
Figure 8: Constrained boundary representation as the result of the reasoning and adjustment process.38 functionally independent constraints have been recognized and enforced, cf.Table 2. The reasoning process led to a more adequate representation of the pyramidal broach roof.

Table 2 :
). Constraints with the bottom area of the building are omitted for clarity.Recognized and actually required constraints for the adjustment.