Automated 3D Scene Reconstruction from Open Geospatial Data Sources: Airborne Laser Scanning and a 2D Topographic Database

Abstract

Open geospatial data sources provide opportunities for low cost 3D scene reconstruction. In this study, based on a sparse airborne laser scanning (ALS) point cloud (0.8 points/m²) obtained from open source databases, a building reconstruction pipeline for CAD building models was developed. The pipeline includes voxel-based roof patch segmentation, extraction of the key-points representing the roof patch outline, step edge identification and adjustment, and CAD building model generation. The advantages of our method lie in generating CAD building models without the step of enforcing the edges to be parallel or building regularization. Furthermore, although it has been challenging to use sparse datasets for 3D building reconstruction, our result demonstrates the great potential in such applications. In this paper, we also investigated the applicability of open geospatial datasets for 3D road detection and reconstruction. Road central lines were acquired from an open source 2D topographic database. ALS data were utilized to obtain the height and width of the road. A constrained search method (CSM) was developed for road width detection. The CSM method was conducted by splitting a given road into patches according to height and direction criteria. The road edges were detected patch by patch. The road width was determined by the average distance from the edge points to the central line. As a result, 3D roads were reconstructed from ALS and a topographic database.

1. Introduction

Three-dimensional (3D) technologies are gaining popularity in industrial, consumer, entertainment, healthcare, education, and governmental applications. According to a new market research report [1], 3D modeling accounted for the highest market share in 2013. In the next five years, 3D modeling and mapping market is expected to grow from $1.1 billion in 2013 to $7.7 billion by 2018. 3D modeling solutions enable users to rapidly construct 3D maps of surrounding areas that are suitable for professional visualization systems. Smartphone companies, such as Google, Microsoft, Apple, and Samsung, have indicated substantial interest in 3D map applications. 3D-enabling devices, such as cameras, scanners, GPS components, and other acquisition devices, have created vast potential for the construction of 3D environments [1]. 3D maps offer some advantages over 2D maps, including providing superior navigation, decision making and information visualization in urban planning and being highly applicable to many smart city concepts and location-based services.

3D reconstruction is the process of determining the shape and appearance of objects using this information to produce virtual replicas of, for example, natural environments, old towns, and archaeological elements. Virtual reality refers to computer-simulated environments that are replicas of the real world or imaginary worlds, but these models include rendered, textured planes and surfaces. Examples of 3D modeling pipelines can be found in [2], and an early example of virtual reality can be found in [3].

The reconstruction of 3D environments requires not only advanced technology but also available data sources. Major data sources used in reconstruction include photogrammetric images, laser scanning and existing map data. Photogrammetry is the technology of deriving 3D data from 2D images by mono-plotting (single-ray back projection), stereo-imagery interpretation or multi-imagery block adjustment. Laser scanning is based on laser (lidar) range measurements from a carrying platform and the precise positioning and orientation of the platform. The costs of data acquisition are usually fairly high and include the costs of equipment, planning, design and measurement. Open data sources, which provide numerous opportunities for exploring the potential of data, are becoming increasingly available. Data are described as “open” if they are freely accessible by anyone to use and republish as they wish, without restrictions from copyright, patents or other mechanisms of control [4]. The advantages of using open data sources are apparent: (i) they provide an economic cost savings; (ii) they maximize the use of data for various purposes; (iii) they stimulate research projects; (iv) they foster the connection between developers; and (v) they maximize achievements.

Open geospatial data have gained substantial popularity. In the past few years, some national governments have provided geospatial data as open data sources on websites to share and explore the potential of data through the development of applications to address public and private demands. In spring 2010, the UK government allowed a significant number of datasets to be freely accessible to the general public via a program named ShareGeo Open [5]. These datasets included many core datasets held by the Office of National Statistics, the Central Government and the Ordnance Survey. Since 1 May 2012, the National Land Survey of Finland (NLS) has allowed its topographic datasets to be freely available to the public and to companies at no cost. According to the NLS’s agreement, the open data product can be used without compensation and with extensive and permanent rights of use [6]. In March 2014, the commercial software company ESRI launched the ArcGIS Open Data site, which enables organizations to create custom open data websites.

As more open geospatial databases become available, the trend to update the geospatial databases from 2D to 3D became evident due to the increasing need for 3D applications, such as flood risk modeling, flight path planning, and environment and coastal protection. In 2012, the Netherlands established the 3D national standard for large-scale topography [7]. In December 2013, Singapore launched a plan for the development and maintenance of a 3D topographic database. In Finland, the need to update a 2D topographic database to a 3D topographic database is urgent because airborne laser scanning (ALS) data continuously proliferate throughout Finland. Therefore, one of our motivations was to investigate the NLS open datasets and to advance the technology for updating 2D topographic databases to 3D.

Building reconstruction from ALS has been continuously developed over the decades since the commercialization of ALS systems. The proposed methods are generally categorized into two groups according to the data source: methods utilizing multiple data sources for building reconstruction, which typically use both ALS data and ground plans, and methods utilizing ALS data alone. The former is less affected by the density of the ALS point cloud due to the use of the supplementary data. For instance, the ground plan may define the outer edges of the buildings, and the inner roof edges may be obtained from the intersections of the neighboring roof patches. However, when using ALS data alone, the density of the points plays an important role during the development of reconstruction methods. Figure 1 depicts building roof point clouds of different densities: 0.8 points/m² (left) and 8 points/m² (right). It can be seen that the sparse data do not necessarily capture the building corners. Connecting neighboring edge points from sparse data will result in jagging line pieces. Meanwhile, the dense data clearly define the edges of the roof. Methods developed for a dense dataset may not be suitable for a sparse dataset. The current trend shows that open geospatial data sources will be widely used in the coming years, especially when considering its free cost and convenience. Therefore, it is urgent in developing methods for 3D scene reconstruction from such datasets.

Figure 1. Comparison of the densities of the ALS point clouds. (a) 0.8 points/m²; (b) 8 points/m².

Figure 1. Comparison of the densities of the ALS point clouds. (a) 0.8 points/m²; (b) 8 points/m².

The main goals of this paper are as follows: (1) to explore the use of open geospatial datasets for 3D scene reconstruction; (2) to develop a pipeline of 3D building reconstruction, including voxel-based roof patch segmentation, extraction of the key-points representing the roof patch outline, step edge identification and adjustment, and CAD building model generation; (3) to conduct 3D road detection and reconstruction from ALS and topographic databases. A constrained search method (CSM) is proposed. The open datasets used in this paper include the ALS and 1:10,000 topographic databases. In the following section, we briefly review the methods proposed in previous studies.

2. Related Work

The previous approaches related to building reconstruction can be grouped into three categories: data-driven, model-driven and a combination of data-driven and model-driven. According to Haala and Kada [8], the building reconstruction methods were classified into three types: (i) reconstruction with parametric shapes; (ii) reconstruction with segmentation; and (iii) reconstruction with digital surface model (DSM) simplification; methods in the first category are model-driven, whereas methods in the latter two categories are data driven. The strengths and weaknesses of the data- and model-driven methods have been discussed in previous studies [9]. For example, data-driven methods are more flexible and do not require prior knowledge; however, the density of the data has a significant effect on the resulting models. Model-driven approaches predefine parametric shapes or primitives, such as simple roof prototypes (e.g., gable, hip, gambrel, mansard, shed and dormer). Complex models are reconstructed by the combination of different primitives. One of the advantages of a model-driven approach is that a complete building roof model can be constructed according to predefined shapes when some building roof data are missing (e.g., due to reflection or an obstacle). However, failure is possible when reconstructing complex buildings and building models that are excluded in predefined shapes [8].

Various researchers have made significant contributions to building reconstruction research, e.g., Vosselman et al. [10,11,12,13], Brenner et al. [14,15,16,17], Elberink [18,19,20,21,22], Haala et al. [8,23,24,25,26], and Rottensteiner et al. [27,28,29,30,31,32]. In addition, more recent studies on building reconstruction have been presented by Bulatov et al. [33], Chen et al. [34], Dorninger and Pfeifer [35], Hron and Halounová [36], Huang et al. [37], Jochem et al. [38], Kada et al. [39,40], Kim and Shan [41], Lafarge et al. [42,43], Perera et al. [44], Rau and Lin [45], Sampath and Shan [46], Seo et al. [47], Xiong et al. [48,49], Yan et al. [50], Yang et al. [51], and Zhou and Neumann [52,53]. The papers selected for this brief review are introduced below.

Vosselman [10] proposed a building reconstruction method using ALS data, in which a 3D Hough transformation was utilized to segment roof patches. In this method, edges are identified using the intersection of faces and an analysis of height discontinuities. Geometric constraints were used to enforce building regularities. In later iterations of this approach [11,12], both ALS and ground plan were used, and a subsequent study [13] applied multiple data sources: ALS, maps and aerial photographs.

Brenner and Haala [14] used DSMs and 2D ground plans as data sources in an automatic and/or semiautomatic reconstruction process. Ground plans are divided into rectangular primitives using a heuristic algorithm. For each 2D primitive, a number of different 3D parametric primitives from a fixed set of standard types are instantiated, and their optimal parameters are estimated. The final object representation is obtained by merging all 3D primitives. In subsequent approaches [15], regularized DSMs and ground plans were utilized. A random sampling consensus (RANSAC) was applied for roof plane segmentation. In a subsequent study [16], the author suggested the use of weak constructive solid geometry (CSG) primitives to model 3D objects. In a later work, Brenner reviewed the methods of building reconstruction from images and laser scanning [17].

Elberink and Vosselman [20] noted problems with using dense ALS data for automated building reconstruction and discussed problems with model-driven methods and the combination of data- and model-driven approaches. Specifically, the following problems were noted: uneven distributions of laser points, inaccurate determinations of the parameters for roof plane segmentation, inconsistencies between point clouds and ground plans when using ground plans, misclassification or incompleteness of laser point data when ground plans were not used, errors in building outline detection due to the missing laser points, the need to rely on assumptions about 3D building shapes, challenges in reconstructing complex building shapes and certain small details when using model-driven methods, and conflicts when applying the thresholds for the final shape of the model when using a combination of data- and model-driven approaches. Elberink and Vosselman [21] employed building vectors and an ALS point cloud with a density of 20 points/m² for building reconstruction. A target-based graph matching approach was developed. A data-driven method was used in tandem with a model-driven method in the process. The test data corresponded to residential areas of Dutch cities with villa and apartment house architectural styles. Complete shape matching was observed for 72% of the 728 buildings. The author noted that reconstruction was not feasible when the buildings contained complex height jumps and flat roofs because the proposed algorithm could not reliably locate all edges of flat roof segments, and the locations of corner points inside the polygon were not detected.

Haala and Brenner [23,24] used laser scanning data and ground plans for building reconstruction. The first step was to acquire a DSM from laser point cloud data. Next, the DSM was simplified to reduce the number of presented points. The ground plan was then decomposed according to the DSM normal. An interactive editing tool was developed to refine the initial reconstruction. Finally, 3D CAD models were reconstructed. In a subsequent approach by Haala et al. [26], the authors presented a cell decomposition method for both roof and facade reconstruction from input data from ground plans, ALS and TLS.

Rottensteiner and Briese [27,28] used laser data (for regularized DSM) and aerial images (to perform segmentation of aerial image grey levels and expansion by region-growing algorithms). In this method, planar roof segments were detected using the DSM normal vectors integrated with the segments from the aerial images. Plane intersections and step edges were detected, and a polyhedral model was derived. Rottensteiner et al. [29] only used laser data. Roof planes were detected using the surface normal vectors. The plane intersects and step edges were then detected. Finally, all step edges and intersection lines were combined to form the polyhedral models. The author’s recent contributions to the ISPRS benchmark for urban object detection and 3D building reconstruction can be found in Rottensteiner et al. [27,31,32]. More information about the ISPRS benchmark will be introduced later.

In a recent study, Dorninger and Pfeifer [35] presented a mean shift segmentation method for roof segmentation and 2D α-shape outline extraction using an angular criterion for final polygon generalization. Sampath and Shan [46] utilized the eigenvalues of the covariance matrix to separate the lines and planes. The planar components were grouped into a small patch. The normal vector of the small patch was determined. The same normal vectors were clustered (using k-means clustering) to form roof patches. Buildings were reconstructed by estimating the intersection lines using a plane adjacency matrix. Lafarge and Mallet [43] proposed a hybrid representation for building reconstruction using geometric 3D-primitives (e.g., planes, cylinders, spheres or cones) to represent standard roof sections and mesh-patches to describe more irregular roof components. Xiong et al. [48] presented a building reconstruction method from a set of noisy point clouds by utilizing roof layer division and global optimized model geometry. In this method, the roof points are decomposed into layers. Each layer is reconstructed independently. The nodes of the roof patches are computed based on the roof planes and their intersections. Zhou and Neumann [53] utilized the inter-element similarities and relationships of the building structures (e.g., direction and placement similarities between planar elements) and detected their locally fitted plane primitives and global regularities to reconstruct the building models. Yang et al. [51] presented a method for building outline extraction. The resulting building outlines consisted of a set of edge points, not key-point-represented outlines. As is evident, many researchers have proposed the use of shape decomposition or the detection of plane intersection lines to segment different roof patches instead of patch outline extraction.

The International Society for Photogrammetry and Remote Sensing (ISPRS) benchmark in urban object detection and 3D building reconstruction was launched in 2012. The results were published in July 2014 [27]. In this benchmark, two test areas with five datasets were provided: the first test contained georeferenced 8-cm ground sample distance (GSD) aerial color infrared images and an ALS with a density of 4–7 points/m² from an area in Vaihingen, Germany; the second test contains a set of georeferenced 15-cm GSD RGB color images and ALS data with a density of approximately 6 points/m² from an area in Vaihingen, Germany. According to the results, 14 different building reconstruction methods were submitted. Ten methods were only based on ALS points, two methods only employed images, one method was based on a raster DSM from ALS and one method used both images and ALS data. The results indicated that ALS data were preferred to images during the building reconstruction. Eight methods were based on generic building models, five methods employed adaptive predefined models, and one method was based on primitives. Data-driven methods are more prevalent in this benchmark. One primitive-based reconstruction method was suitable for large roofs but produced many errors in small superstructures on roofs. During the reconstruction process, under-segmentation was the dominant error type for areas with small buildings, whereas over-segmentation errors were common for areas with large roofs. Thus, Rottensteiner et al. [27] concluded that “for future research, the focus can be on smaller roof structures and a better treatment of step edges in complex flat roof buildings”.

With regard to road reconstruction, many studies have provided evidence that an ALS point cloud is suitable for road detection and reconstruction. Clode et al. [54] proposed a method for classifying roads using both the intensity and range of LIDAR data. Jan et al. [55] introduced a temporal filtering method based on conditional random fields for the detection and reconstruction of the curbs and road surfaces. Zhang [56] utilized the height information of the point cloud and pattern recognition techniques for road candidate selection. The line representation of the candidates was then determined. Finally, road surfaces and road edges were extracted.

Recently, mobile laser scanning (MLS) has been widely for fast and accurate road modeling. Recent studies based on MLS point clouds for road modeling can be found in Boyko and Funkhouser [57], Yang et al. [58], and Sukhan and Park [59]. Boyko and Funkhouser [57] utilized a road map and a large-scale unstructured 3D point cloud for 3D road reconstruction. Yang et al. [58] extracted road markings from MLS data using georeferenced point cloud images and the strength of the reflectance. Sukhan and Park [59] presented a drivable road model using four-layered laser scanners and one camera, and the positions of the vehicles and the obstacles were identified.

Although some studies of 3D road modeling have been implemented, there is still a room for further development. In this paper, we developed approaches tailored for upgrading 2D roads from the NLS topographic database (presented as the central lines of the road carriageways) to 3D road models (featuring the road edges and heights) and investigated the availability of ALS point clouds of different densities.

3. Data Sources

The open geospatial data file download service has been provided by the NLS since May 2012. Customers are able to select their desired product, format and coordinate system as well as their areas of interest within Finland. The data sources that we used include two sets of ALS point clouds with different densities, a 2D topographic database and orthophoto.

3.1. Airborne Laser Scanning

ALS data provide 3D coordinates of point cloud and laser echo information. The coordinate system is ETRS-TM35FIN. ALS is currently available in some parts of Finland. Data acquisition for the entire country is planned for 2019 [6]. The program of utilizing ALS is based on the technology transfer from the Finnish Geodetic Institute to the NLS. From the NLS, the concept was also modified to Swedish National Laser Scanning.

Open ALS data from NLS have a density of 0.8 points/m². Its vertical accuracy is typically 15 cm, whereas its planimetric accuracy is approximately 35 cm. Quality control for the Finnish national ALS data has been addressed by Ahokas et al. [60].

The second set of ALS data, which was provided by the local municipality, has a density of 8 points/m², a vertical accuracy of 10 cm, and a planimetric accuracy of 25 cm. The objective of using two sets of ALS with different densities is to investigate the feasibility of road edge extraction from each set and determine their applications.

3.2. Topographic Database

The Finnish topographic database contains vectors of buildings, transportation networks, administrative borders, geographic names, land use, and waterways. The positional accuracy of the topographic database corresponds to the scale range of 1:5000–1:10,000, which indicates 1 m for the position accuracy and 0.25 m for the vertical accuracy. The traffic road network and geographic names are continuously updated, the buildings, constructions and administrative borders are annually updated and the remaining elements are updated approximately every five to 10 years. The topographic database is available as a version that encompasses the entire country [6]. The formats for the topographic database include Esri shape, MapInfo MIF, and also MAAGIS/XL and GML.

3.3. Orthophoto

NLS orthophotos comprise a color aerial photo data set with a terrain resolution of 0.5 m, which encompasses the entire country of Finland. NLS orthophotos are updated every three to 10 years. The product is a part of the NLS raster material. Color orthophotos can be used as textures for terrain presentations. The format is the jpeg2000 format (. jp2).

4. Developed Methods for 3D Building Reconstruction

This section presents the developed methods for 3D CAD building model reconstruction from a sparse ALS point cloud (0.8 points/m²). We assume that the ALS point cloud has been classified into buildings and ground. The procedure for building reconstruction is illustrated in Figure 2. “CAD” is the abbreviation of “computer-aided design”. In this paper, a CAD model mainly refers to an optimized model. It is represented by a 3D wireframe model. In this section, we also used the terms such as “key point” and “roof patch”. A key point of a building model is the point in which direction or height changes. A key point heavily affects the description of the shape of the building. “Roof patch” is a plane which contains all coplanar points. In a simple case, one roof has one patch. In most of cases, a roof contains several roof patches, which are possibly intersected (e.g., gable and hip roofs) or parallel with different height patches.

During the CAD building model reconstruction, the core issue is the extraction of the key points representing the roof patch outlines. The following steps are suggested for obtaining the final building models: (i) voxel-based building roof patch segmentation and the extraction of the key-points representing the outline of each patch; (ii) establishment of the data structure for the roof patches (the building direction, plane normal, patch neighboring relations and inclusive relationships); (iii) the acquisition and adjustment of step edges; and (iv) building base height acquisition and meshing.

Figure 2. Procedure for 3D CAD building model reconstruction.

Figure 2. Procedure for 3D CAD building model reconstruction.

4.1. Roof Patch Segmentation

The input for the roof patch segmentation algorithm was implemented based on the classified individual buildings (Figure 3a). In the planar patch segmentation, we utilized voxel-based processing initially developed for the efficient processing of mobile laser scanning [61] point clouds, which has also been found to be feasible for the processing of ALS point clouds. First, the point cloud was divided into a regular voxel structure, and only occupied voxels (henceforth “voxels” for brevity) were retained and saved into an array. Each element of the array contained the $x$-, $y$-, and $z$-coordinates of the points of one voxel. For each voxel, the indices of the voxels in the 124-neighborhood were retrieved and saved in a separate table. Second, locally planar voxels were classified via principal component analysis (PCA) to obtain the variances ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3$ in the three principal directions [62]. A voxel was classified as planar if the points in its 124-neighborhood were located approximately within a plane. We used a flatness measure $F=1-{\lambda }_3/\left({\lambda }_1+{\lambda }_2+{\lambda }_3\right)$ [63] in the classification: if $F$ was greater than 0.95, the voxel was labelled as planar. Third, planar seed surfaces were extracted from the planar points (Figure 3b) using a modified version of the USF range image segmentation algorithm [64]. The modifications were as follows: (i) instead of pixels, we used voxels; (ii) a voxel’s normal was calculated via PCA using all points in the voxel’s 124-neighborhood; (iii) we chose the seed voxel randomly and calculated the initial plane parameters from the seed voxel’s 124-neighborhood; (iv) we used 124-connectivity; (v) we did not compare angles between the normal in the growing phase; (vi) the perpendicular distance between a new candidate voxel and the plane of the surface grown so far was calculated between the candidate voxel’s center of mass and the plane (we used a threshold of 0.2 m); (vii) we did not calculate the point-to-point distance but accepted all 124-neighbors that satisfied perpendicular distance condition (vi); (viii) we re-estimated the plane parameters of the surface after five new voxels had been added; and (ix) patches smaller than 2 voxels were discarded. Fourth, the seed patches were used to obtain connected planar surface patches from the original point cloud (Figure 3c). The fourth phase was necessary because, for example, all points on the planes are not classified as a planar patch due to noise and roof crests. The fourth phase began by selecting one seed patch, fitting a plane to it and selecting all points of the building point cloud that resided on the seed plane. Then, selected points were merged with the seed points, duplicates were removed, and connected components were retrieved from the points. If the largest connected component that overlapped the seed patch was sufficiently large, it was classified as a planar patch. The points on the planar patch were removed from the building point cloud and the remaining seed patches. Next, a new seed patch was selected, and the same procedure was repeated for the remaining points. This process continued until all seed patches were treated. An application of the algorithm to a more complex roof structure is shown in Figure 4. The height jumps of different patches are displayed on the right side of Figure 4. Different colors represent different patches.

Figure 3. Procedure of plane segmentation.

Figure 3. Procedure of plane segmentation.

Figure 4. Example of roof patch segmentation with height jumps. (a) ALS building points; (b) detected roof patches shown in different colors; (c) roof patches at different heights.

Figure 4. Example of roof patch segmentation with height jumps. (a) ALS building points; (b) detected roof patches shown in different colors; (c) roof patches at different heights.

4.2. Building Roof Data Structure Establishment

After roof patch segmentation, the data structure for the building roofs should be established. In this section, an example is given to show how a building roof data structure is treated. Figure 5 displays the numbering of the roof patches.

Figure 6 shows an example of an inclusive relationship. First, patch 1 and patch 2 are neighboring relations. If the centroid of patch 1 (red point) is located on patch 2 and the height of patch 1 is greater than the height of patch 2, $patch 1\in patch2$. These patches have an inclusive relationship if the following conditions are satisfied: (i) two planar patches are neighboring relations; (ii) the centroid of one patch is located inside another neighboring plane (in 2D space); and (iii) one plane is higher than another plane.

Figure 5. Segmented building roof patches with marked patch numbers.

Figure 5. Segmented building roof patches with marked patch numbers.

Figure 6. The determination of an inclusive relationship: $Patch 1\in Patch 2$.

Figure 6. The determination of an inclusive relationship: $Patch 1\in Patch 2$.

Table 1 shows the data structure for the roof patches. The data structure in Table 1 contains multiple feature descriptions: the consistency in the patch direction and building direction, intersecting angle of the patch normal and vertical direction, inclusive relations, neighboring relations and height of the patch centroid. The patch and building directions can be obtained by PCA of the points. The normal of a patch is defined here as the mean value of the normal of all points of the patch. The normal angle is the intersection angle of the normal and vertical directions. The relationship between two neighboring planar patches can be determined by calculating the 2D shortest distances (−xy plane) between the boundaries of two patches. When the shortest distance is less than a threshold, e.g., 1.5 m, we refer to them as neighboring relations.

Table 1 reveals that all patches have consistent directions in an entire building. The patch normal angle is defined as a zero horizontal plane and a 90° vertical plane. As shown in Table 1, six smaller patches contain angles greater than 10°. Two patches (P12 and P15) contain angles greater than 20°, which indicates that the large patches are more stable than the smaller planes in normal angle estimation. Table 1 also shows that three pairs of patches satisfy the criteria of an inclusive relationship: P7 and P4, P15 and P2, P16 and P3.

Table 1. Data structure of building roofs.

Plane (Number of Points >10)	Number of Points	Plane and Building Direction (Consistent = 1)	Plane Normal Angle (degree)	Inclusive Relationship $\left(Pi\in Pj\right)$	Neighboring Relations	The Height of the Centroid of the Plane (m)
P1	1652	1	1.0896	0	P2, P3, P5, P8, P10, P11	14.87
P2	1550	1	−5.5742	0	P1, P3, P6, P12, P15	22.2
P3	976	1	3.3197	0	P2, P1, P12, P16	26.09
P4	648	1	3.8889	0	P6, P7, P13	22.23
P5	386	1	3.7306	0	P1, P9, P10, P11, P14	22.91
P6	243	1	−8.8889	0	P2, P4, P5	11.35
P7	237	1	7.5949	$P7\in P4$	P4	25.65
P8	74	1	14.5946	0	P1, P10	11.66
P9	56	1	−12.8571	0	P5	19.69
P10	48	1	−7.5	0	P1, P5, P8	19.60
P11	53	1	6.7925	0	P1, P5	18.36
P12	66	1	21.8182	0	P2, P3	15.49
P13	51	1	−14.1176	0	P4	7.1
P14	19	1	−18.9474	0	P5	18.2
P15	14	1	−25.7143	$P15\in P2$	P2, P6	26.81
P16	30	1	12	$P16\in P3$	P3	27.24

Table 1. Data structure of building roofs.

Plane (Number of Points >10)	Number of Points	Plane and Building Direction (Consistent = 1)	Plane Normal Angle (degree)	Inclusive Relationship $\left(Pi\in Pj\right)$	Neighboring Relations	The Height of the Centroid of the Plane (m)
P1	1652	1	1.0896	0	P2, P3, P5, P8, P10, P11	14.87
P2	1550	1	−5.5742	0	P1, P3, P6, P12, P15	22.2
P3	976	1	3.3197	0	P2, P1, P12, P16	26.09
P4	648	1	3.8889	0	P6, P7, P13	22.23
P5	386	1	3.7306	0	P1, P9, P10, P11, P14	22.91
P6	243	1	−8.8889	0	P2, P4, P5	11.35
P7	237	1	7.5949	$P7\in P4$	P4	25.65
P8	74	1	14.5946	0	P1, P10	11.66
P9	56	1	−12.8571	0	P5	19.69
P10	48	1	−7.5	0	P1, P5, P8	19.60
P11	53	1	6.7925	0	P1, P5	18.36
P12	66	1	21.8182	0	P2, P3	15.49
P13	51	1	−14.1176	0	P4	7.1
P14	19	1	−18.9474	0	P5	18.2
P15	14	1	−25.7143	$P15\in P2$	P2, P6	26.81
P16	30	1	12	$P16\in P3$	P3	27.24

4.3. Extraction of Roof Patch Outlines

The roof patch outlines are extracted by the following procedure:

(a) Rotate the planar patch to be horizontal (Figure 7a,b), retrieve the minimum bounding box (MBB) and rotate the patch again by setting the sides of the MBB as the new coordinate axes;

The planar patch was projected into the plane in which it lies, and a MBB was retrieved for the projection in the following way. First, a MBB with an arbitrary rotation was retrieved, and its area was calculated. Next, the rotation was varied in small steps, and the MBB and its area were retrieved in each step. The rotation that provided the MBB with the smallest area was selected. Finally, the projections were rotated such that the sides of the MBB were parallel to coordinate axes.

(b) Transform the horizontal patch to a binary image and obtain the MBB of the image (Figure 7c);

(c) Calculate the solidity of the image. The solidity of an image refers to the proportion of pixels in the convex hull that are also in the region, that is, image actual area divided by its convex area. If the solidity value is less than the threshold, e.g., 0.94, the image contains irregular regions. Thus, its complementary image and the image’s MBB are obtained (Figure 7d).

(d) Repeat (c) for all images inside the MBBs, that is, divide the image into sub-images until the solidities of all sub-images are greater than a threshold (e.g., 0.94) (Figure 7e). When the area of the sub-image is less than a threshold, it is not considered for the MBB extraction.

(e) Combine all MBBs and sort the order of the vertices of the polygon;

Figure 7. Outline extraction from binary images: (a) roof patch point cloud; (b) roof patch rotation to horizontal or vertical direction; (c) extract the MBB of the patch; (d) obtain the MBB of the complementary image of (c); (e) acquire the MBB of the complementary image of (d); (f) combine the boundaries of (c), (d) and (e) to extract the outline of the roof patch; (g) achieve the final outline by rotating the outline of (f) to the original direction.

Figure 7. Outline extraction from binary images: (a) roof patch point cloud; (b) roof patch rotation to horizontal or vertical direction; (c) extract the MBB of the patch; (d) obtain the MBB of the complementary image of (c); (e) acquire the MBB of the complementary image of (d); (f) combine the boundaries of (c), (d) and (e) to extract the outline of the roof patch; (g) achieve the final outline by rotating the outline of (f) to the original direction.

(f) Transform the result from (e) to 3D points (the reverse processing of (b)) (Figure 7f);

(g) Rotate the result from (f) to the original position (the reverse rotation processing of (a)) (Figure 7g);

(h) Overlap the key points from (g) with all points of the roof patch and remove the key points that have no neighboring points within a certain distance (e.g., 2 m).

(I) Extract the final outline by sorting the vertices as a closed polygon.

The building in the example contains 16 planar patches. Figure 8 illustrates all outlines of the extracted roof patches. However, the neighboring patches have inconsistent edges. Additional processing is needed to achieve consistent edges.

Figure 8. Superimpose the outlines of all planes on the building point cloud.

Figure 8. Superimpose the outlines of all planes on the building point cloud.

4.4. Acquisition and Adjustment of Step Edges

We obtained individual roof patch outlines from Section 3.3. However, the neighboring patches have inconsistent edges. The step edge is an abruptly height change roof edge connecting its neighboring patches (see the left image of Figure 9). From the patch neighboring relationship list in Table 1, each patch must address the edges of its neighboring patches. The neighboring relationships between two edges of neighboring patches must be identified. The shorter edge is projected onto the longer neighboring edge to ensure consistency in the edges. In this stage, the edges are still considered in 2D. The right picture of Figure 9 shows the results after edge adjustment.

Figure 9. (Left) Step edge represented by red lines; (Right) After edge adjustment.

Figure 9. (Left) Step edge represented by red lines; (Right) After edge adjustment.

4.5. Building Base Height Acquisition and Model Generation

Previous sections have presented the procedure for extracting the roof patch outlines and achieving consistent edges between two neighboring patches. We now have the outlines of the roof patches. The base height of a building is its ground-level height. ALS building points only include roof parts, and ALS data do not include building facade data. The building base height can be obtained from ALS ground points by projecting the building key points onto the ground. However, not all planar patch key points are projected to the ground because some planes are located on top of other planes (the inclusive relationship). Upper patch key points are vertically projected onto its lower planar patch. Then, the ALS ground points are triangulated to form a triangulated irregular network (TIN). Building base heights are calculated by vertically projecting the key points of the buildings onto the TIN model. The ground heights of the building key points are obtained via bilinear interpolation. An average ground height is applied to achieve the leveling building base lines. Building models are constructed by meshing the key points in each vertical plane. A vertical plane contains both roof points and building base points. Figure 10 illustrates the wireframe models and solid models of the buildings.

Figure 10. Building model generation: wireframe model and solid model. Roofs are represented by the red areas.

Figure 10. Building model generation: wireframe model and solid model. Roofs are represented by the red areas.

5. Developed Methods for 3D Road Detection and Reconstruction

Road networks provide abundant information for a 3D scene. A 3D road network can be reconstructed by employing 2D carriageway central lines from a topographic database and ALS data. A topographic database provides both 2D carriageway central lines and the road attributes, e.g., the “class”. The road “class” was defined by the surface materials: motorway, main suburban road, railway, and pedestrian way. Each “class” described the range of the road or carriageway widths e.g., class “12111”: motorway, each carriageway with a width of 10–11 m (refer to Table 2). We utilize the width information to obtain the location of the road edges from the ALS point cloud. For example, for a motorway with a carriageway width range of 10–11 m, the edge search range (R) can be extended to a width of e.g., 5 m in each side. If one side shows slope change and another is flat, then the neighboring carriageway of the flat side is searched, until slope changes are detected. If there is no height change within the search range, it is accepted as a central carriageway. Thus, neighboring carriageways in the both sides need to be retrieved to find the slope change edges.

Table 2. Transportation data (the width of a carriageway) from the NLS topographic database.

Class	Legend	Width (m)
12111	Motorway	10–11
12112	Highway	8–10
12121	High-capacity urban roads	6.5–8
12122	Main suburban roads	5–6.5
12131	Suburban roads and entry	4–5
12132	Small suburban roads with gravel surface	3–4
12141	Roadway	≤3
12313	Pedestrian and bicycle way with gravel surface	≤2
12314	Pedestrian and bicycle way with asphalt surface	≤2
12316	Footpath	1
14111	Electrified railway	1.52 m

Table 2. Transportation data (the width of a carriageway) from the NLS topographic database.

Class	Legend	Width (m)
12111	Motorway	10–11
12112	Highway	8–10
12121	High-capacity urban roads	6.5–8
12122	Main suburban roads	5–6.5
12131	Suburban roads and entry	4–5
12132	Small suburban roads with gravel surface	3–4
12141	Roadway	≤3
12313	Pedestrian and bicycle way with gravel surface	≤2
12314	Pedestrian and bicycle way with asphalt surface	≤2
12316	Footpath	1
14111	Electrified railway	1.52 m

Figure 11. Procedure of 3D road network reconstruction.

Figure 11. Procedure of 3D road network reconstruction.

Figure 11 illustrates the 3D road network reconstruction procedure. The method was developed according to the following steps: (i) obtain the elevations of the carriageway or road central lines from ALS ground points; (ii) determine the ALS search area for the road edge detection; (iii) road patch separation; (iv) use the discrete Laplacian method for road surface detection; (v) calculate the width of the road and create parallel road lines on both sides of the central lines according to the detected road width; and (v) mesh the roads to form a 3D road network.

5.1. Obtaining the Elevations of 2D Central Lines of the Roads or Carriageways from ALS Ground Points

First, the ALS ground points (also called terrain points) are triangulated to form a TIN. The 2D carriageway nodes from topographic databases are projected onto the TIN model, and each node is located in a certain triangle (see Figure 12). The bilinear interpolation method is used to calculate the node height.

Figure 12. Road node projected onto a ground TIN model.

Figure 12. Road node projected onto a ground TIN model.

5.2. Determination of the ALS Search Area for Detecting the Road Edges

According to the road “class” definition, the width range for each carriageway is known. Figure 13 shows one example of a motorway.

For example, consider “12111” motorway with a carriageway width range of 10–11 m. The actual edges could be positioned between the green lines. The extended area for searching the edges is defined as the minimum width plus 10 m; each side is about 5 m. An example of an extended area can be found in Figure 14, in which the red lines represent the extended border for ground point selection.

Figure 13. Example of the width range for a carriageway defined by the road “class” in topographic databases.

Figure 13. Example of the width range for a carriageway defined by the road “class” in topographic databases.

Figure 14. Road patch segmentation. (Left) lines of segmentation (blue line) vertical to the road surfaces through the key points; (Right) ALS ground points segmented into multiple patches.

Figure 14. Road patch segmentation. (Left) lines of segmentation (blue line) vertical to the road surfaces through the key points; (Right) ALS ground points segmented into multiple patches.

5.3. Road Patch Separation

One road typically contains many nodes, and the road may be located on an undulating terrain. When we detect the road edges, ground points are separated into patches according to the detected key points. The algorithm for the key point detection will be described in the text below. The blue lines in the left image of Figure 14 represent the separation lines perpendicular to the road surfaces through the key points. These lines intersect with the extended border lines (red lines) to form the intersection points. The connections between the intersection points create the patch areas. The right image in Figure 14 illustrates ALS ground points located in different patch areas using different colors.

The method for the key point extraction was developed according to two criteria. When one of these criteria is satisfied, we accept the point as a key point.

(i) Direction change in the -xy plane

We consider Figure 15 as an example. The angle between the neighboring lines, e.g., P3P4 and P4P5, is calculated. If the intersection angle is greater than 15°, the intersection point is accepted as a key point. In Figure 15, the angle between P3P4 and P4P5 is 27°; P4 becomes a key point.

(ii) Height change

From the direction change detection, P3 and P4 were selected as key points. Then, we must calculate the height changes between the neighboring points. For example, the road height ranges from P1 to P5: [3.48, 3.34, 3.24, 2.95, 3.42].

For P2, h2−h1 = −0.14, h3−h2 = −0.1; For P4, h4−h3 = −0.29, h5−h4 = 0.47.

From the estimation of the height changes, two criteria are applied for the key point selection:

Identical or different signs in the height differences of the neighboring points:

Different signs (e.g., −0.29, 0.47) indicate that the terrain contains ups and downs (e.g., P4). Consider P4 as an example (although P4 has been selected as a key point during the direction change detection); P4 is also selected as a key point according to this criterion.

The height differences between neighboring points:

If the height difference between neighboring points exceeds the threshold (e.g., 0.3 m), we retain the point as a key point. In the example, h5−h4 = 0.47. Thus, we select h4 as a key point.

From these results, P3 and P4 are selected as key points. The ortholines through P3 and P4 perpendicular to the road central line are created. The ortholines separate the road surfaces into multiple patches.

Figure 15. Example of a road with direction and height changes.

Figure 15. Example of a road with direction and height changes.

5.4. Application of the Discrete Laplacian Method for Road Detection

The road surface was first extracted roughly using Laplace’s differential operator, which can be described as

$$L=\frac{∆V}{4}=\frac{1}{4}\left(\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\partial }^{2}V}{\partial y^2}\right)$$

(1)

where $V$(x,y) is a function of two variables and the interior points of $L$ are obtained by taking the difference between a point in $V$ and the average of its four neighbors. The discrete Laplacian calculates the values on the edges of $L$ by linearly extrapolating the second differences from the interior. We employed the discrete Laplacian method for the height changes with respect to the neighbors. The results shown in Figure 16 are obtained by thresholding the discrete Laplacian. The red points represent the approximate result for the road surface and the blue points represent the ALS points.

Figure 16. Application of the discrete Laplacian method for road detection. (Red) Road surface points detected by the discrete Laplacian operator; (Blue) ALS ground points.

Figure 16. Application of the discrete Laplacian method for road detection. (Red) Road surface points detected by the discrete Laplacian operator; (Blue) ALS ground points.

After the rough extraction of the road surface, histogram analysis was applied to the previous results to remove the majority of the scattered noise. A bivariate histogram (-xz plane) is established to analyze the number of points located in the bins. When the number of points in a bin is less than the threshold, these points are removed.

Figure 17 shows the histogram analysis and results. The road trace is distinct on the left side of the figure. However, the results on the right side of the figure show some scattered points around the road surfaces. Therefore, we employ the range search to estimate the number of nearest neighbor points to remove the scattered points. Figure 18 shows the results after scatter point removal.

Figure 17. Road detection using a histogram analysis. (Left) Histogram image; (Right) the results of the road detection.

Figure 17. Road detection using a histogram analysis. (Left) Histogram image; (Right) the results of the road detection.

Figure 18. Results of the road surface detection.

Figure 18. Results of the road surface detection.

In the case of multi-carriageways motorways or highways, since each carriageway has one central line, a road consists of multiple central lines. The slope change happens on both sides of the road. For a central carriageway, there is no height change on its sides. Therefore, when we use the CSM, if no slope change was detected, the neighboring carriageway was retrieved until the edge with slope change was found.

5.5. Calculation of the Road Width and Estimation of the Road Edges

After extracting the road surface from the patches, edges are generated according to the rule of parallel to the road center line. The road width can be calculated after merging all patches. The road width was determined by the average distance from the edge points to the central line. The mean width can be applied (see Figure 19). For multi-carriageways roads such as motorway or highway, the road width is the widths of multi-carriageways combination.

Figure 19. The road edges are generated by creating parallel lines from the central lines of the road according to the calculated road width.

Figure 19. The road edges are generated by creating parallel lines from the central lines of the road according to the calculated road width.

5.6. Meshing

When the geometries of the roads have been defined by the sets of nodes and their topological relationships, the shapes and surfaces of the roads can be generated by triangulating the nodes, which is also referred to as meshing. The meshed polygons can be visualized in 3D computer graphics and solid models.

6. Experimental Results and Analysis

We utilized the datasets described in Section 2 for 3D scene reconstruction. The test areas were located in Tapiola-Otaniemi-Keilaniemi area, in southeast Espoo, Finland. A set of ALS building points was utilized for building model reconstruction. Two sets of ALS point clouds (densities of 0.8 points/m² and 8 points/m²) and one set of topographic data were applied to test the road detection and reconstruction, and their applicability was investigated. The resulting 3D scene includes 3D terrain model, 3D buildings and 3D roads.

6.1. 3D Terrain Model

The 3D terrain model was generated from ALS ground points. Due to a large number of points from ALS, we have developed a quad-tree algorithm for ground point reduction [65]. After performing this algorithm, the terrain model with representation of reasonable points can be achieved. The orthophoto is applied for the ground texture (Figure 20).

Figure 20. Terrain model for an area of 6 km × 6 km. (Left) terrain model without texture; (Right) terrain model with an orthophoto texture.

Figure 20. Terrain model for an area of 6 km × 6 km. (Left) terrain model without texture; (Right) terrain model with an orthophoto texture.

6.2. 3D Building Models

The original building point cloud has a density of 0.8 points/m². The building types in the test area have both flat and oblique, regular- and irregular-shaped roofs. The resulting models are shown in Figure 21. The evaluation is conducted based on three aspects: (i) roof patch segmentation; (ii) height difference between models and building point clouds; and (iii) assessment of the distance between the model points and their nearest points in the laser data.

Figure 21. 3D building models of the Tapiola area.

Figure 21. 3D building models of the Tapiola area.

Figure 22. Extracted building roof patches. The different colors represent different patches.

Figure 22. Extracted building roof patches. The different colors represent different patches.

6.2.1. Assessment of roof patch segmentation

Plane extraction is an important step during reconstruction, as it directly affects the number of building reconstructions. Figure 22 illustrates the extracted roof patches. The patches are presented in different colors. The evaluated point cloud contained 236 planar building patches, and the algorithm detected 271 patches. A total of 176 patches, that is, 74.6% of the building patches, were correctly detected. Forty-one extracted segments consisted of noise segments that contained non-planar objects, such as vegetation or parts of one or more roof patches. Fifteen instances of undersegmentation and three instances of oversegmentation were detected; all oversegmented patches were divided into two patches. The extracted patches also included 12 parapets on the edges of the roofs, 16 ambiguous segments, which could be classified as the ground or buildings, and five planar patches that were not part of buildings. The errors from roof patch segmentation were caused by various factors, such as the density of ALS data; the results from building classification and the complex of the building structures and so on.

6.2.2. Assessment of the height difference between the models and original building points

Figure 23 shows the resulting models for the height evaluation. The left image illustrates the reconstructed building models, in which roofs are shown in red and the building facades are shown in gray. In the right image, the building models (in gray) overlap the building point cloud (in green). We selected 15 test locations distributed evenly throughout the test area, where we measured the height of both 3D building models and the original building point cloud. We compared the heights of the models and original points in the locations and estimated the average accuracy of the building models with respect to the original building points. Table 3 reveals the height differences between the models and original ALS roof points, which vary from 0 to 0.28 m. We employed the root mean square error (RMSE) to measure the deviation:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{t=1}^n{\left(Z_t-Z\right)}^2}{n}}$$

where $Z_t$ is the predicted value, $Z$ is the observed value, and $n$ is the number of predictions.

Figure 23. Building model height evaluation. (up) Building models, with roofs in red and facades in gray. (down) Original ALS building roof points in green and building models in gray.

Figure 23. Building model height evaluation. (up) Building models, with roofs in red and facades in gray. (down) Original ALS building roof points in green and building models in gray.

In this context, predicted value is the value obtained from a reconstructed model, and the observed value is the value obtained from the original ALS point cloud. The RMSE was estimated from the difference showed in Table 3 as 0.18 m. The average height difference between the models and the original points was 15 cm. The heights of the model key points were determined by their neighboring points. The referenced height of a roof patch was calculated as the average height of the patch points. Therefore, the difference is reasonable if the deviation between the model and the roof patch points is less than the tolerance of the planar patch detection.

Table 3. Height differences between the ALS building roof points and 3D roof models in 15 test locations.

Test Location Index	Height of ALS Building Points (m)	Height of Building Models (m)	Height Difference (m)
1	31.61	31.59	0.02
2	26.84	26.74	0.10
3	20.04	19.91	0.13
4	20.76	20.52	0.24
5	25.04	24.93	0.11
6	23.81	23.76	0.05
7	22.72	22.44	0.28
8	26.86	26.61	0.25
9	24.91	24.86	0.05
10	25.36	25.10	0.26
11	23.21	23.02	0.19
12	21.81	21.57	0.24
13	25.89	25.81	0.08
14	20.51	20.33	0.18
15	19.72	19.58	0.14
Average			0.15
RMSE			0.18

Table 3. Height differences between the ALS building roof points and 3D roof models in 15 test locations.

Test Location Index	Height of ALS Building Points (m)	Height of Building Models (m)	Height Difference (m)
1	31.61	31.59	0.02
2	26.84	26.74	0.10
3	20.04	19.91	0.13
4	20.76	20.52	0.24
5	25.04	24.93	0.11
6	23.81	23.76	0.05
7	22.72	22.44	0.28
8	26.86	26.61	0.25
9	24.91	24.86	0.05
10	25.36	25.10	0.26
11	23.21	23.02	0.19
12	21.81	21.57	0.24
13	25.89	25.81	0.08
14	20.51	20.33	0.18
15	19.72	19.58	0.14
Average			0.15
RMSE			0.18

6.2.3. Assessment of the distance between the model points and their nearest points in laser data

In our test data, 176 roof patches were reconstructed, with 1142 model key points extracted. These model key points were compared to their nearest laser points, and the distances were calculated. The results showed that 6% of the distances were larger than one pixel (e.g., 1.2 m in the test data), while 2.5% were larger than two pixels. The pixel size used in model reconstruction was determined by the space between the input points. Therefore, the use of denser laser points decreases the deviation between the model point and its nearest laser point.

6.2.4. Discussion concerning quality assessment approaches for reconstructed 3D building models

Elberink and Vosselman [22] presented a quality assessment approach for reconstructed 3D building models. The authors proposed using residuals between model faces and laser data and the nearest distance between model points and laser data to measure the quality of the reconstructed building models. Rottensteiner et al. [27] employed several quality metrics for evaluating the building reconstruction results from ISPRS benchmarks: (a) completeness, correctness and quality on a per-roof-plane basis; (b) completeness, correctness and quality on a per-roof-plane basis for planes larger than 10 m²; (c) over-segmentation, under-segmentation, and both cases; (d) the RMSE of the planimetric distances of the reference roof plane boundary points to their nearest neighbors on the corresponding extracted roof plane boundaries; and (e) the RMSE of the height difference between reference planes and all corresponding extracted planes. The reference planes were manually acquired from aerial stereo-images. The quality indices (a), (b) and (c) describe the quality of the roof plane segmentation, whereas (d) and (e) represent the geometrical accuracy of the roof polygons. Our assessment method is similar to that of Elberink and Vosselman [22] with a notable difference being that they evaluated the deviation of the model from each roof patch point.

6.3. 3D Road Networks

3D road networks were reconstructed based on the ALS point cloud and 2D topographic database. Two sets of ALS point clouds with different point densities (0.8 points/m² and 8 points/m²) were tested. The heights of the carriageway central lines can be acquired from either sparse or dense ALS data. However, the road widths are crucial for the 3D road reconstruction. The solutions for road width acquisition are different based on two sets of different density point clouds. Figure 24 illustrates the difference between the point clouds of different densities. The green lines represent the range of the road widths, which are defined in the topographic database. For the example of a road with a width of 10–11 m, the inner green lines represent the width of 10 m, whereas the outer green lines represent a width of 11 m. As shown in the left image in Figure 24, few points are located between the widths of 10 m and 11 m. It is difficult to accurately detect the edges from the sparse data. A solution for using a sparse dataset for road width acquisition is to utilize the width information from a topographic database. An average width value, e.g., 10.5 m in this case, can be applied for the road reconstruction. The dense point cloud is applied for the automatic edge detection.

Figure 24. ALS datasets and road with a width range of 10–11 m. (Left) ALS with a density of 0.8 points/m²; (Right) ALS with a density of 8 points/m².

Figure 24. ALS datasets and road with a width range of 10–11 m. (Left) ALS with a density of 0.8 points/m²; (Right) ALS with a density of 8 points/m².

The approach proposed in Section 4 is conducted for the dense dataset. The objective is to robustly achieve road edges. Here, the road edges refer to the height-changing edges. A road typically contains the road surfaces and road shoulders. In this study, the road includes the road surfaces and road shoulders. In order to exhibit the datasets and the detected roads, Figure 25 gives an illustration. Red lines in the left figure represent the 2D central lines of the carriageways, which were from topographic database. The others are ALS point cloud. In the right image of Figure 25, it shows the resulting road surface overlaying on the meshed terrain model. It can be used as visual inspection and accuracy evaluation. The final 3D road networks are shown in Figure 26. The reconstructed road types include motorways, highways, urban roads, suburban roads, pedestrian and bicycle roads and tracks. The road widths vary from 1 m to approximately 40 m. In the right image of Figure 26, the 3D road models are overlapped with the orthophoto.

Figure 25. Road detection and reconstruction from a 2D topographic database and ALS. (Left) red lines are the central lines of the carriageways from topographic database. Others are ALS point cloud; (Right) detected road points overlaid with meshed terrain (dark blue).

Figure 25. Road detection and reconstruction from a 2D topographic database and ALS. (Left) red lines are the central lines of the carriageways from topographic database. Others are ALS point cloud; (Right) detected road points overlaid with meshed terrain (dark blue).

Figure 26. The result of 3D road network reconstruction. (Left) 3D road networks; (Right) 3D roads overlaid on the NLS orthophoto (roads in red).

Figure 26. The result of 3D road network reconstruction. (Left) 3D road networks; (Right) 3D roads overlaid on the NLS orthophoto (roads in red).

The evaluation of the road widths is challenging because the practical road widths are not defined by the height jumping changes. To prevent an inconsistent road edge definition, we evaluate our method by comparing the detected road edges with the interactive measurement from the dense ALS points. The height differences between the ALS points and 3D road models are also included in the assessment process. Eight roads are selected for the assessments. Table 4 presents the accuracy assessment of the widths and heights of the roads. The average width deviation is 22 cm. Its RMSE was estimated as 0.23 m. The average height error is 14 cm, with an RMSE of 0.16 m. The width errors are caused by a) noise points around the roads, such as low vegetation close to the road edges, which produces larger detected road widths compared with the practical widths; and b) road parallel assumption. We use the average width and create parallel lines from the central lines to obtain the road edges. In practice, the edges including the road shoulders exhibit zigzag edges, which produce an inconsistency between the actual parallel edges and ideal parallel edges. The height errors may be attributed to the variations in the road surface. The road heights are extracted from the ALS data. The nearby points of each road node that is located in the ALS ground TIN model are included in the height estimation. An average height value of the nearby points is assigned to the height of the road node. The same height value is applied for the corresponding edge nodes. In most cases, road surfaces are not flat, which creates the potential for errors.

Table 4. Road width and height accuracy evaluation.

Road No.	Road Data from Reference Data (m)		Detected Road Width (m)		Difference (Absolute Value) (m)
	Width	Height	Width	Height	Width	Height
Road 1	39.71	6.26	39.90	6.11	0.19	0.15
Road 2	8.26	8.54	8.14	8.26	0.12	0.28
Road 3	26.45	5.92	26.23	5.82	0.22	0.10
Road 4	32.18	5.14	32.43	5.35	0.25	0.21
Road 5	19.32	8.42	19.50	8.28	0.18	0.14
Road 6	14.89	8.06	14.55	7.98	0.34	0.08
Road 7	9.38	10.17	9.08	10.06	0. 30	0.11
Road 8	6.71	12.08	6.56	12.11	0.15	0.03
Average					0.22	0.14
RMSE					0.23	0.16

Table 4. Road width and height accuracy evaluation.

Road No.	Road Data from Reference Data (m)		Detected Road Width (m)		Difference (Absolute Value) (m)
	Width	Height	Width	Height	Width	Height
Road 1	39.71	6.26	39.90	6.11	0.19	0.15
Road 2	8.26	8.54	8.14	8.26	0.12	0.28
Road 3	26.45	5.92	26.23	5.82	0.22	0.10
Road 4	32.18	5.14	32.43	5.35	0.25	0.21
Road 5	19.32	8.42	19.50	8.28	0.18	0.14
Road 6	14.89	8.06	14.55	7.98	0.34	0.08
Road 7	9.38	10.17	9.08	10.06	0. 30	0.11
Road 8	6.71	12.08	6.56	12.11	0.15	0.03
Average					0.22	0.14
RMSE					0.23	0.16

According to this analysis, we construct the following recommendations:

(i) When only sparse ALS data are available, the average values of the road width ranges from the topographic database can be applied for 3D road reconstruction;

(ii) When dense ALS data are available, the proposed method can be applied for road edge extraction.

7. Conclusions

In this study, we investigated the availability of open geospatial datasets for updating a 2D topographic database to a 3D database and reconstructed a 3D scene from open data sources. We utilized sparse ALS data and developed a building reconstruction pipeline. The pipeline includes voxel-based roof patch segmentation, extraction of the key-points representing the roof patch outline, step edge identification and adjustment, and CAD building model generation. In addition, we developed methods for road detection and reconstruction using a 2D topographic database and ALS point clouds. Our investigation shows that when the density of the point cloud was low (e.g., 0.8 points/m²), it was difficult to acquire the road edges. When a denser point cloud (8 points/m²) was applied, the road surfaces could be successfully detected using a discrete Laplacian operator. The results were validated for both 3D building models and 3D road models, and 74.6% of the roof patches were correctly identified for the data with a low point density of 0.8 points/m². The average height deviation between the models and the laser points was 16 cm. The deviation of the distance between the model points and their nearest laser points was estimated. 6% of these distances were greater than one pixel, while 2.5% were greater than two pixels. The pixel size used in model reconstruction was determined by the space between the input points. The results showed that as the point density increased, the models approached the original laser points more closely. Furthermore, the 3D road reconstruction provided a height deviation of 14 cm and width accuracy of 22 cm. These results provide evidence of the great potential of open geospatial data sources for the low cost 3D scene reconstruction. Compared to previous studies, our contributions lie in the following:

(i)

The availability of open datasets for 3D scene reconstruction has been demonstrated. The use of open datasets reduces the cost of 3D scene reconstruction.

(ii)

The proposed method provided a means of reconstructing 3D buildings from sparse datasets. Our method works for both sparse and dense datasets. A dense point cloud provides higher model accuracy. In contrast, many previous methods developed based on dense datasets may not be applicable for sparse datasets.

(iii)

Our method can produce CAD building models for different roof types, e.g., flat and oblique, regularly and irregularly shaped. Furthermore, no extra step is needed to enforce edge to be parallel or building regularization.

(iv)

Our study has shown the potential of upgrading the NLS 2D topographic database (e.g., buildings and roads) to a 3D database.