RECONSTRUCTION OF 3 D MODELS FROM POINT CLOUDS WITH HYBRID REPRESENTATION

The three-dimensional (3D) reconstruction of urban buildings from point clouds has long been an active topic in applications related to human activities. However, due to the structures significantly differ in terms of complexity, the task of 3D reconstruction remains a challenging issue especially for the freeform surfaces. In this paper, we present a new reconstruction algorithm which allows the 3D-models of building as a combination of regular structures and irregular surfaces, where the regular structures are parameterized plane primitives and the irregular surfaces are expressed as meshes. The extraction of irregular surfaces starts with an oversegmented method for the unstructured point data, a region growing approach based the adjacent graph of super-voxels is then applied to collapse these super-voxels, and the freeform surfaces can be clustered from the voxels filtered by a thickness threshold. To achieve these regular planar primitives, the remaining voxels with a larger flatness will be further divided into multiscale supervoxels as basic units, and the final segmented planes are enriched and refined in a mutually reinforcing manner under the framework of a global energy optimization. We have implemented the proposed algorithms and mainly tested on two point clouds that differ in point density and urban characteristic, and experimental results on complex building structures illustrated the efficacy of the proposed framework.


INTRODUCTION
Buildings are the dominant type of man-made objects and these 3D building models are necessary to visualize and analyse the urban environment (Brook et al., 2013).Its reconstruction and modelling, due to its vast applications in many areas, have long been an active topic in many research communities, from city planning, computer vision, geomatics, photogrammetry, to remote sensing.During past decades, many algorithms and systems have been proposed to reconstruct building models, based on the type of input data: images (Bulatov et al., 2014;Toschi et al., 2017) vs LiDAR (Rychard and Borkowski, 2016;Xu et al., 2017), or their combinations (Aringer and Roschlaub, 2014;Sun and Salvaggio, 2013), the degree of automation: interactive (Nan et al., 2010;Xiong et al., 2014) automatic (Lafarge and Mallet, 2011) the reconstruction scale: single building (Oesau et al., 2016) vs block/city (Duan and Lafarge, 2016;Verdie et al., 2015), and the output models: 2.5D (Chauve et al., 2010;Rouhani et al., 2017), or full 3D (Rychard and Borkowski, 2016), which make the production of 3D models better and faster.Even though much progress has been made, the reliable generation of detailed building models remains to be a challenging issue (Musialski et al., 2013), and existing algorithms meeting present requirements are still in the phase of development (Rottensteiner et al., 2014), especially for the buildings with complex structures (e.g.irregular surfaces).

RELATED WORK
The problem of building modelling has received considerable attention in recent years.Given the large body of work in this broad topic, areas of 3D reconstruction and primitives extraction related to our approach are reviewed.

3D reconstruction
Among this huge variety of 3D building reconstruction methods proposed in the literatures, we can distinguish three basic types, namely, data-driven methods, model-driven methods, and hybrid-driven methods.
Data-driven methods: a common assumption for these approaches is that a building is a polyhedral model, and each individual element is reconstructed directly based on the geometric information (Sampath and Shan, 2010;Zhou, 2012).It usually starts with the extraction of planar patches from LiDAR data using segmentation algorithms, followed then a polyhedral building model is generated from the planar patches by intersection.In additional, some regularization rules are usually applied to enhance these models.The main advantage of the data-driven approach is that it can reconstruct polyhedral buildings with complex shapes, while the drawback is their sensitivity to the incompleteness of data arising from occlusion, shadow and missing information, and so on.
Model-driven approaches: it searches the most appropriate models that are fitted by the corresponding given point clouds among primitive building models contained in a predefined model library (Huang, 2013;Karantzalos and Paragios, 2010), but only relatively simple blocks were allowed to reconstruct in the early stages.The model-driven approaches are robust to reconstruct buildings because some of the previous knowledge like parallel and symmetry can be easily combined, and can create a watertight geometric model with semantic labelling.The main weakness is that they cannot generate a satisfactory result of various building styles except for certain types of models in the predefined primitive library which usually have relatively simple roof styles.Due to the growing needs for the accurate 3D building models and these shortcomings in both of above approaches, an alternative hybrid reconstruction method (Elberink and Vosselman, 2009;Rychard and Borkowski, 2016) is gradually being concerned by the public, which combined the conventional model-driven methods with data-driven methods.
Hybrid reconstruction: the main objective of the hybrid reconstruction, based on the widely used roof topology graphs (RTG), is to recognize building structure and search the best fitting primitives from the predefined library.Such an idea is first seen in (Verma et al., 2006) and the normal vector is added as an attribute mark to the RTG, while the primitive types are limited to some simple ones and significantly reduces its application.(Oude Elberink and Vosselman, 2011) expand the library of parametric primitives and introduce more attributes like convexity or concavity to reconstruct objects by the subgraph matching, and it can be easily prone to errors when mismatches of sub-graphs occur.Similar works were explored to track the improvements in reliability and availability of RTG, and were further developed to distinguish the roof features and to interpret building structures (Perera and Maas, 2014).Recently, (Rychard and Borkowski, 2016) proposed an automated method to avoid multiple matching of the same roof elements and to interpret the semantic knowledge of roof elements based on RTG database.However, an inherent problem to the hybrid reconstruction based on the RTG can be easily suffered error-prone in itself (e.g.roof plane was lost because noise, data missing and more in local plane extraction), and the completeness stored in RTG library are rather limited leading to the fact that all building types in the real world cannot be reconstructed.What's more, such complex 3D models are expressed as some meaningful parametric primitives and difficult to interpret the irregular surfaces.

Primitives extraction
Extensive studies have been done to improve the robustness and efficiency of plane primitives extraction, which can be roughly categorized into four categories: methods based model fitting (Chen et al., 2014;Schnabel et al., 2007), methods based region growing (Deschaud et al., 2010;Vo et al., 2015;Yang and Dong, 2013), methods based feature clustering (Kim et al., 2016;Zhou et al., 2016), and methods based global energy optimization (Pham et al., 2016;Yan et al., 2014).Even though these proposed methods can generally provide satisfactory extraction results, there still exists limitations to extract primitives from point clouds.The results are non-deterministic and sensitive to the choice of parameters for the existing greedy methods as region growing, model fitting, and feature clustering, while approach by global energy optimization is hard to optimise.Approaches for irregular surface extraction is usually performed by Hough Transform (Duda and Hart, 1972) and Random Sample Consensus (Fischler and Bolles, 1981), which heavily required unambiguous formula of objects.Inspired by (Dong et al., 2018), our proposed approach detects the irregular surface firstly, and formulates the energy for the multi-scaled plane voxels not the points, which is much simpler and admits efficient approximate optimisation algorithms.

METHODOLOGY
The proposed method, as illustrated in Figure 1, encompasses three key components that extract potential free-form surfaces, segment the planar primitives, and reconstruct the hybrid models.

Generate super-voxels
Collapse super-voxels based region growing

Global energy optimization
Extract freeform surfaces

Construct hybrid models
Hybrid models combined with mesh and planes Taking the individual point cloud as input, an oversegmentation approach is introduced to divide the unorganized point clouds into super-voxels (Section 3.1.1),and a region growing approach based the curvature and normal difference of the adjacent super-voxels is used to collapse similar supervoxels (Section 3.1.2).The nonplanar points belong to the irregular surfaces are extracted by a thickness threshold (Section 3.1.3).Plane primitives are segmented by two steps: generating candidate planes by multi-scale segmentation (Section 3.2.1),and optimizing these candidates with energy minimization based super-voxels (Section 3.2.2).A hybrid models combined irregular meshes and regular planes are constructed in the Section 3.3.

Extract Freeform Surfaces
The freeform surfaces are always founded in the urban environments, increasing the complexity of 3D modelling, and its models are usually reconstructed by the split small plane primitives.In this paper, these irregular surfaces are directly extracted and clustered as an individual object not a series of small plane patches.Three steps are included: generating supervoxels, collapsing the similar voxels based region growing, and grouping the super-voxels with a thickness threshold.

Generation super-voxels
As the reconstructed point cloud contains several millions of points, we represent the point cloud by a collection of small super-voxels to reduce the computational cost and cluster the points with the similar attributes.We used the supervoxel segmentation approach proposed by (Papon et al., 2013), which over-segments the input point cloud as a set of small patches . Each contains a normal vector i n , a centroid i c , and the curvature i f .In additional, the method can also define an adjacent graph between super-voxels, and E is a collection of edges between adjacent surface patches.We use the resulting graph G for the further region growing base voxels.In this work, we are only interested in geometric features, so we use the spatial distance, and normal vector deviation for super-voxel generation.

Collapse super-voxels
To group the over-segmented super-voxels with the similar features, a region growing approach is applied to the edges of G .It merges two adjacent super-voxels (shared a valid edge E ) that within a fixed threshold of curvature and normal difference.The optimal features i sv of the merged are calculated by the Least Square Estimation and Principal Component Analysis.
After merging two voxels, the adjacency graph G as well as the features i sv are updated.It should be noted that this update is local as only edges with the super-voxels adjacent to the two merged are operated.

Grouping irregular voxels
In order to find an irregular surface from the merged supervoxels in G , a thickness feature from (Viejo and Cazorla, 2014)   is introduced to extract the irregular voxels, which can be measured the dispersion angle for the points in the planar patch neighbourhood, calculate by: ) arctan( Where value of i λ is the singular values in descending order.
Thickness angle gives us an idea about how super-voxel in a neighbourhood fit to a planar surface.We use it to discard patches with a thickness angle greater than 0.04 degree.The obtained irregular voxels will be grouped by K-means clustering to generate the individual freeform objects.

Segment plane primitives
Taking the left planar super-voxels as input, an efficient global energy optimization method is used to precisely segment the plane patches.There are three key steps as multi-scale supervoxels segmentation to obtain the candidate planes, and a global energy optimization to extract the detailed planes.

Generation multi-scale super-voxels
Fine detailed planes may be contained in the collapsing operation or the over-segmentation.Thus, to generate multiscale super voxels, a flatness check for the left planar supervoxels is used to find the potential small planes.Voxels will not be further divided only if the descending singular values are satisfied with the following conditions： Where α s and β s are the scaling factors.Multi-scale super- voxels are generated in an iteration manner by: Where scale R is the initial radius used in the above section of freeform surface extraction.The value k is the number of iteration and is always set to 5.These over-segmented voxels from multi-scaled segmentation can express the regular scene as much as possible.In additional, the super-voxel that no longer to be split can be expressed as a plane by:

Global energy optimization for planes
The candidate planar patches detected from multi-scaled segmented will be optimized by a global energy minimization solution as Eq. ( 5).This problem can be resolved via a popular graph cuts like the extended α -expansion algorithm (Isack and Boykov, 2012), which achieves the balance of geometric errors (data cost), spatial coherence (smooth cost), and the number of planes (label cost).
In the formula (5), the data cost term ∑ − p p L p measures the sum of geometric errors and is calculated as the quadratic perpendicular distance between centroid of super-voxel p and Label p L as Eq. ( 6), especially the use of quadratic distance is equivalent to assuming Gaussian distribution for errors.The second term in Eq. ( 5) is a smoothness prior.It assumes some specific neighbourhood system edge for the adjacent super-voxels and is constructed by Triangulated Irregular Networks.The indicator function (.) δ is selected as Potts Model.The value of (.) δ is 1 only if the specified pair of centroid points p and q to belong to the same plane; while 0 otherwise.The weight pq ω is set inversely proportional to the distance between centroid points of voxel p and q, as Eq. ( 7).
Where closer supervoxels are a priori more likely to fit the same plane.
The label cost term in Eq. ( 5) penalizes the number of planes (Labels).Fewer planes are encouraged to be used to represent the point cloud compactly.We used numbers of inliers point ( i L ) from each super-voxel in the label cost written by: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2, 2018 ISPRS TC II Mid-term Symposium "Towards Photogrammetry 2020", 4-7 June 2018, Riva del Garda, Italy The plane optimization based on global energy minimization is an iterative process and terminates when the global energy is no more decreased, and we refit the plane parameters for the derived labels of the super-voxels.The filtered points from the global energy optimization are attached to the nearest irregular surface.

Reconstruct hybrid models
The final model is a hybrid representation where the extracted irregular surfaces are represented by the meshes (Triangulated Irregular Networks), while the segmented planes are organized by geometric primitives.To simply the problem and get a compact hybrid models, we firstly search the closest point between vertexes of meshes and planes, and make an interpolation operation in the common transition regions to get the new linked lines.

RESULTS AND DISCUSSION
We have implemented our proposed algorithms and mainly tested on two datasets that differ in point density and feature characteristic, as illustrated in Figure 1.The first one (S1) is a stand-alone ALS building point cloud with 17,510 points, 23 points/m2; while the other (S2) is a merged building data from the Multi-View Stereo (MVS) and TLS with 594,871 points, 45 points/m2.In addition, the merged data is in the uneven distribution as the different acquired platforms.To extract the freeform surfaces, a multi-scaled segmentation is used to generate the super-voxels with an iteration manner.In our experiments, an initial scale R is set to ten times of the average resolution; the curvature and normal difference for collapsing super-voxels are set to 0.05 and 10 degree for both of the data.The thickness feature used to extract the irregular super voxels is equal to 0.05 degree for the two data.For both of these data, the numbers of iteration is set to 5, while the values of α   And due to the lack of independent reference data, we have also conducted an evaluation with visual inspection for the reconstructed hybrid models in Figure 3.The basic planar primitives are well reconstructed including the narrow planes in the common transition regions.In additional, the points belong to the planar planes are used to evaluate internal quality measures.The average distance between a point to the reconstructed plane is 3.8 cm (S1) and 6.2 cm (S2), respectively.The larger average bias for S2 is mainly caused by the imprecise point clouds generated by Multi-View Stereo matching.

CONCLUSION AND FUTURE WORK
In this paper we present a new scheme for building modelling from 3D point cloud.By first extracting the points of individual free-form objects, we develop and apply a global energy optimization approach to obtain reliable parametric planes, and the output of the reconstruction is a visually pleasing hybrid models combining meshes and plane primitives.The key ideal is to extract the potential free-form building points and precisely abstract the potential planar planes, which can be useful to obtain a more completed model and avoid the failures of plane segmentation caused by the outliers, noise, occlusion, and clutter in the point clouds.In the future, it would be interesting to improve the robustness in the scheme of freeform points extraction with adaptive parameters, and to smooth the linking between the meshes and planes.Other challenge is to adapt our approach to point clouds generated from internet photo collections which contain more outliers and have spatial distributions highly heterogeneous.

Figure 1 .
Figure 1.The illustration of the proposed method
s and β s are 20 and 5, respectively.These candidate planes are optimized by global energy minimization, where the parameters λ and β are fixed to 1.5 and 20.These parameters produce good results and have been used for all the examples of this paper.As shown in Figure2, points belong to the irregular surfaces are grouped in blue, while the segmented planes are marked respectively in different colors.

Figure. 2 .
Figure.2.Results of extracted freeform surface and planes.Points belong to freeform surfaces are in blue colour, while optimized planes are marked with individual colour.S1 data is shown in the top and S2 in the bottom.

Figure. 3 .
Figure.3.Final output hybrid models.Hybrid models for S1 and S2 data are illustrated in top and bottom, respectively.