AUTOMATED ROAD INFORMATION EXTRACTION FROM HIGH RESOLUTION AERIAL LIDAR DATA FOR SMART ROAD APPLICATIONS

Automatic extraction of road features from LiDAR data is a fundamental task for different applications, including asset management. The availability of updated and reliable models is even more important in the context of smart roads. One of the main advantages of LiDAR data compared with other sensing instruments is the possibility to directly get 3D information. However, the task of deriving road networks form LiDAR data acquired with Airborne Laser Scanning (ALS) may be quite complex due to occlusions, low feature separability and shadowing from contextual objects. Indeed, even if roads elements can be identified in the ALS point cloud, the automated identification of the network starting form them can be involved due to large variability in the size of roads, shapes and presence of connected off-road features such as parking lots. This paper presents a workflow aimed at partially solving the automatic creation of a road network from high-resolution ALS data. The presented method consists of three main steps: (i) labelling of road points; (ii) a multi-level voting scheme; and (iii) the regularization of the extracted road segments. The developed method has been tested using the “Vaihingen”, “Toronto” and “Tobermory” data set provided by the ISPRS. * Corresponding author


INTRODUCTION
Nowadays, smart roads are more and more important and relate to a large number of other topics such as digital transformation of infrastructures (Meijer et al, 2018), autonomous driving (Varaiya, 1993), connected vehicles (Lu et al, 2014), etc. In particular, the possibility to combine traditional GIS products with BIM-based outputs and especially the possibility to parameterize some road geometric features can create new workflows for road management, design and for vehicle-toinfrastructure communication (Ndashimye et al, 2017). For those purposes, the availability of reliable and updated road information systems is of paramount importance. LiDAR (Light Detection and Ranging) data can provide reliable, accurate and repeatable information to feed smart road systems. One of the main advantages of LiDAR data compared with other sensing instruments is the possibility to directly get 3D information, to be used to separate buildings and trees from roads. In addition, due to its relatively narrow scanning angle (Rottensteiner and Clode, 2008) Airborne laser scanning (ALS) data are generally free of serious occlusions of the road surface. Finally, the possibility to have full waveform information and study different reflections allows to detect road in rural areas under tree coverage. Although LiDAR data presents a high potential for smart road application, the original point cloud cannot be directly used for further analysis and integration with other data sources. However, the transformation of raw LiDAR data into exploitable products and their integration into BIM/GIS environments is not a trivial task and still requires expert personnel. Clode at al., (2007) presented one of the first approaches for automatic road extraction from LiDAR data. This method identifies road candidate points by filtering from a given distance from a Digital Surface Model (DSM) derived by LiDAR data and integrating laser intensity (see Scaioni et al., 2018) of points to remove bare grounds. Starting from road points, road patches were firstly extracted and later connected into a road network. Template matching was developed to identify the road network in Zhao and You (2012). In this framework a two-stage procedure was developed. In the first stage, LiDAR points are labelled as "ground" and "off-ground" points. In the second step, a template scheme is used to search for roads on ground intensity images. Road widths and orientations are determined by a subsequent voting scheme. The main trend in recent years is associated with the combination of different algorithms. For example, Hu et al, (2014) presented the combination of three algorithms: Mean Shift Clustering (MSC), Tensor Voting, and Hough Transform. First, MSC is used to detect road centre points. In the second step, Tensor Voting is implemented to highlight the main linear features. Finally, Hough Transform allows to detect road centrelines. Similarly, Li et al. (2015) developed a multi-step semiautomated procedure. A manual election of road seed points is carried out in the first step. Then road areas are detected by using a region growing method starting from the manually identified seeds. In the last step, fast parallel thinning is used to extract road centrelines. A proper selection of the seeds can be carried out to remove incorrectly extracted roads. Another multi-step approach was developed by Hui et al. (2016). Firstly, Skewness (see Crosilla et al., 2013) balancing is used to define an optimal intensity threshold for road points. Secondly, narrow streets are removed with Rotating Neighbourhood algorithm. Finally, road centrelines are derived by using a Hierarchical Fusion and optimization technique. Another combination of algorithms for road centreline extraction is presented in Li et al. (2016). Here the three main steps are: (i) the detection of the road centre point by using adaptive MSC, (ii) local principal component analysis for extracting linear distributed points, and (iii) hierarchical grouping for connecting primitives into a complete road network. Tejenaki et al. (2019) detect road centrelines by combining intensity data and normalized DSM. MSC is used to filter intensity data, that are combined with different normalized DSM products to minimize the effects of large parking lots and the like. In the final stage, road centrelines are extracted by using a Voronoi-diagram-based approach and then by removing dangle lines. Wen et al. (2019) present a Deep Learning approach aimed at detecting road markings by using deep learning and, a modified U-net framework starting from intensity images. Similarly, Jung et al. (2019) developed a multi-step procedure for lane marking extraction using a rasterized version of the point cloud. The procedure is quite involved and the involved steps include image segmentation, morphological filtering, and lane association. Some general aspects may be highlighted from this overview of the existing literature. First, although LiDAR intensity is a fundamental aspect for identifying roads, the separation between road points and bare soil points is still problematic. In addition, different roads may present multiple intensity values connected with the multiplicity of road construction materials (e.g., asphalt, concrete, etc.). A second issue for an automated method is related to the influence of attached areas (e.g., parking lots) on road extraction. Irregular point distribution and variations of road pattern/width may also affect the efficiency of road detection. To partially cope with this aspect, in this paper we are presenting an automated procedure for road information extraction, i.e., road centreline and road width, that can be used in combination with a BIM/GIS framework. The paper is organized as follows: Section 2 presents an overview of the proposed method and the possible integration of road centreline extraction in a BIM/GIS framework. Details on the developed road centreline methodology are discussed in Section 3. Section 4 shows some results of an experimental study carried out to validate the developed method. In particular, the "Vaihingen", "Toronto" and "Tobermory" data, provided by the International Society of Photogrammetry and Remote Sensing (ISPRS), were used in the experiments. The last Section 5 draws some conclusions and future work.

OVERVIEW OF THE DEVELOPED METHODOLOGY
The proposed workflow for extraction of road centreline starting from ALS data is presented in Figure 1. The first step of the proposed approach is the classification of the acquired point cloud by using Random Forests (Liaw and Wiener, 2002). This first classification allows for the definition of two main groups of points: "ground points" (including roads, bare soil, grasslands, parking lots, etc.) and "off-ground points" (including trees, buildings, cars, people, etc.). Starting from "ground points" the road extraction workflow is composed of five main steps, as shown in Figure 1. The main assumption supporting the multiscale linearity voting is that points belonging to road elements are aligned along a line. Since there are multiple typologies of rods in a road network system in terms or road width and length, a multiscale approach is used to identify points of a road element. In particular, the probability that a point belongs to a road is evaluated at different scales and a final voting array is determined. In the second step, a first set of raw centrelines are extracted. Indeed, starting from the votes received in the voting phase, LiDAR intensity and point planarity, MSC is used to classify "road points." Then α-shape (Fayed and Mouftah, 2009) is used to build road polygons and subsequently to derive a first set of raw centrelines. Regularization of centrelines is the aim of next processing step. Centreline smoothing is obtained by clustering centrelines according to proximity and centreline orientation. Finally, since the extracted centrelines may not be connected to each other a cell complex is created starting from the regularized centrelines and a connected road network is derived by using a minimization criterion.

ROAD CENTRELINE EXTRACTION
Starting from the original ALS point cloud the first step is the classification of the point cloud into three sets of object classes: "ground", "building" and "tree". In particular, the classification is carried out by using Random Forest (RF) classifier. A RF classifier produces multiple decision trees using a randomly selected subset of training samples and variables. In other words, RF randomly and iteratively samples the selected variables to generate a large set (named as "forest") that represents the statistical behaviour of numerous decision trees. To combine the votes over the constructed trees a majority vote is generally used while a bagging strategy is generally used to create a training set from the original data. The bagging randomly selects about two thirds of the samples from a training data to train these trees. Then, the remaining samples, generally called "out-of-bag" (OOB), are used for cross-validation and to estimate the classification error. Features used to perform the classification are subdivided into three main categories: • Height-based; • Eigenvalue-based; and • Local-plane-based.
Height based features considered in this work are: • Δz: defined as the height difference between the point and the lowest point found in a cylinder of a radius of 7 m. This parameter allows to discriminate between "ground" and "off-ground" points; and • σz 2 : the height variance computed for 50-nearest neighbouring points.
Where λ1> λ2> λ3 are the eigenvalue of the variance-covariance matrix computed considering the 50-nearest neighbouring points. Eigenvalue-based features allow to discriminate between man-made objects and trees. Finally, local-plane-based features are: • Nz: the deviation angle of the normal vector of the fitted plane from the vertical direction; • Ri: residuals of the local estimated plane; and • N 2 : the variance of the point normal with respect to the normal vector of the fitted plane.
Local-plane features are computed on the basis of Least-Squares fitting of the 50-nearest neighbouring points. Training sets are manually selected to train the RF according to the previously defined classes: "tree", "ground" and "building". Once the classification is carried out the layer named as "ground" contains roads, parking lots, bare ground, low-land grass, etc. To distinguish among the categories in the previous listed points, the main idea exploited in this paper is that the geometrical distribution of points into a road is different compared with the one we can observe in parking lots and bare ground. In particular along road, points are more densely sampled and present smoother surfaces than in correspondence of bare soil and under trees. In addition, roads present points distributed along a line while parking lots and bare ground show an equally-distributed point density. A specific multiscale linearity voting was designed to exploit this idea (Figure 2). A road network may be characterized by roads with different width and lengths. The idea developed in this paper is that at an appropriate scale road points can be identified as a dense liner segment forming a smooth surface. In this paper a scale is defined as the size of the side of a square centred in a road point candidate. The definition of the scales for the analyse depends on: • road's width: the minimum scale corresponds to the minimum road width; • road's average length: influences the definition of minimum and maximum scales; and • road density: the maximum scale should be defined not to include multiple roads inside a single scale.
Once a specific scale is defined (Figure 2) a square around each "ground" point (named as "seed") is investigated. By using Random Sample Consensus (RANSAC - Fischler and Bolles, 1981), a line is firstly fitted to the local set of points. In RANSAC a tolerance value equal to 1/3 of the side of the clustering window (Figure 3a) is considered for linear fitting. A minimum number of support points is also used for line acceptance, that is equal to 1/3 of the points in the square (Figure 3b). To prevent possible over-segmentation, i.e., to identify those linear features that are not really existing, region growing is then performed for all points in the clustering window using as new seeds the ones detected by RANSAC. The linearity of the grown region is then evaluated by computing the linearity coefficient: where λ1> λ2 are the 2D eigenvalues of the variance-covariance matrix of the points in the search window. If Lλ> Llim the feature is classified as linear and all points belonging to it are assigned with a vote. Llim is a user-defined threshold to discriminate between linear and no-linear elements at the defined scale. If Lλ< Llim no votes are assigned to these points and the original "seed" point is penalized. If the line is accepted, the local surface smoothness (planarity) is computed as: ( 2) where zj is the elevation of a point pj and ̅ z is the average elevation of all points in the line. Similarly, the planarity of any remaining points pi in the search window is computed. Once the linearity voting is carried out, the separation among different point classes can be performed (Figure 4). Indeed, the The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) specific spatial point distribution associated with each object's class directly reflects into a different linearity vote. At this stage, besides road points, the detected "ground" points contain parking lots, bare ground, and low grass. However, compared to parking lots, bare ground, and low-grass road points are assigned with a higher linearity vote. Therefore, the subsequent process is finalized to recognize road points by votes and to extract the primitives of road centrelines. To separate "road" points from "non-road" points, MSC is applied.
The following characteristics are evaluated per each point: the number of votes received at each specific scale, the laser intensity, and the point smoothness. MSC allows for the detection of group of points having similarities in terms of these three parameters: • points characterized by a lower number of votes and lower smoothness at each level can be recognized as either "off-ground points" or "bare ground/low grass;" • points characterized by a lower number of votes and high smoothness can be categorized as "parking lot;" and • points with higher votes and smoothness at least in one level are assigned to the "road" class.
After voting points classified as "roads" can be transformed into a polygon representation by using the α-shapes. An α-shape is defined as a frontier, which is a linear approximation of the original shape. Such α-shapes can be used to perform the boundary reconstruction from an irregular point cloud. The parameter α controls the precision of the boundary. In particular, the value of α represents the radius of a rolling around the point cloud. The rolling track of the circle form the boundary of the point-set, which also allows for the regularization of the shape. In particular, the α parameter is chosen equal to the average point density in the data set. Some spurious points still remain in the "road" data set, mainly in road shoulder, that makes this estimate of the road boundaries quite inaccurate. Starting from the polygon of the road, the centreline can be extracted by considering only the boundary points and by determining the centreline with the methods based on Thiessen polygon (Brassel and Reif, 1979).
Due to the noise in the original polygon the outcomes of the centreline may result quite jagged. In addition, some road intersection areas are missing. Indeed, due to the characteristics of the voting scheme, intersection areas tend to be classified as "parking lot". For this reason, a regularization is taken out to form a complete road network. In particular, to reconstruct the final road network a hierarchical grouping method is adopted. Firstly, adjacent primitive with small collinear thresholds are connected into longer road segments. In a second stage connectivity is established by building a cell complex by connecting and intersecting these smoothed road lines in a way similar to the one presented in Previtali et al., (2018). The main aim of this procedure is to create a connected set of lines by minimizing the length of "guessed" connections needed to establish a closed network among really "detected" smoothed centrelines.
Finally, a removal of short lines which are likely not meaningful roads is carried out.

EXPERIMENTS AND EVALUATION
To assess the efficiency of the proposed approach, an experimental study was performed taking into consideration different detests. In particular, three different data sets were evaluated: 1. "Vaihingen" Data set is provided by the "ISPRS Test Project on Urban Classification and 3-D Building Reconstruction." This data set represents a typical situation of a European small town with low rise buildings and irregular road networks in the city centre. LiDAR data were acquired with a Leica ALS50 system at point density of approx. 4 points/m 2 ; 2. "Toronto" Data set is provided by the same ISPRS Test Project od Data set 1. In this case, a typical North American city with high-rise buildings and regular road network with main roads aligned along two orthogonal directions is concerned. Data set 2 was provided by Optech Inc. at point density of approx. 6 points/m 2 ; and 3. "Tobermory" Data set is delivered by the ISPRS TC III/WG V and represents a typical small town in North America with a quite small downtown and a significant extension of secondary roads (dirty roads). Data set 2 was captured with an Optech Titan LiDAR system at point density of approx. 10 points/m 2 .
As it can be observed these data sets are covering quite different typical scenarios. In addition, they present a large number of elements that may prevent an efficient extraction of roads such as: parking lots, trees covering a significant portion of the road, grasslands, squares, roads with different widths and directions. The list of parameters used for the different data sets are presented in Table 1.
To evaluate the efficiency of the extracted roads ( Figure 5) results have been compared to the ones achieved on the basis of OpenStreetMap (OSM) data, which have been used as reference. The evaluation of these outcomes has been carried out by comparing the distances between the "extracted" and the "refences" lines. A preliminary global alignment of the two data sets ("extracted" and "reference" lines) has been carried out to take into consideration a possible co-registration bias. After the alignment, the distances between corresponding features in both data sets have been computed. A statistical analysis (Scaioni, 2010) of the results is presented in Table 2 while a visual representation of the discrepancies is shown in Figure 6. As it can be observed, the results for the "Vaihingen" and "Toronto" Data sets are similar and outperforms the results obtained from the "Tobermory" Data set. This difference may be due to the fact that while the "Vaihingen" and "Toronto" Data sets refer to urban environments, the "Tobermory" Data set is mainly focusing on a rural area. In this case, the definition of the road extremes is more complicated for the presented algorithm.

Data sets Statistics on discrepancies
Vaihingen  Table 3. Statistics on results obtained for the Vaihingen", "Toronto" and "Tobermory" Data sets.

Orthoimages Segmentation results Extracted road centrelines
" Toronto" Data set " Tobermory" Data set Figure 5. Results of the presented methodology applied to the "Vaihingen" (on the upper row), the "Toronto" (on the middle row), and the "Tobermory" (on the lower row). Data sets: orthoimages of interested areas are shown on the leftmost column; results of the segmentation process are shown in the central column , where ground segments are depicted in green and the remaining classes in blue; and the obtained road centrelines are displayed in the rightmost column.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 6. Comparison between the "extracted" and "reference" centrelines obtained from OpenStreetMap (OSM) for the "Vaihingen" (on the left), "Toronto" (on the centre) and "Tobermory" (on the right) Data sets, respectively. Roads are colorized according to discrepancy with respect to the "reference" data. Figure 7. Cumulated percentage of discrepancy for the "Vaihingen" (on the top), "Toronto" (in the middle) and "Tobermory" (at the bottom) Data sets, respectively. Discrepancies are computed with respect to "reference" centrelines from OSM.
Indeed, some points wrongly classified as "road", belonging in realty to the "bare soil" segment, may influence in a negative way the results of the road extraction. However, in the case of the "Tobermory" Data set it is worth to observe that the developed road extraction methodology has identified a significant number of secondary dirty roads that are not present in the OSM "reference" data. Finally, we have to observe that, even if the maximum discrepancy in the "Tobermory" data set is equal to 53 m (corresponding to an existing road not detected by the developed method) the 80% of the data set presents a discrepancy lower than 3.0 m (Figure 7). In the Vaihingen" and "Toronto" Data sets the maximum discrepancies are 5.5 m and 4.5 m, respectively. In addition, these discrepancies are localized in correspondence of road intersections. Due to the characteristics of the developed methodology road intersections are generally not directly detected (they are recognized as parking lots) and reconstructed only in the final step as result of the reconstruction of the cell complex. This may determine a localization of the largest discrepancies in the correspondence of road intersections. In the case of the "Vaihingen" and "Toronto" Data sets the 80% of the extracted road network presents a discrepancy lower than 0.52 m and 0.92 m respectively (Figure 7).

DISCUSSION AND CONCLUSIONS
Accurate information about road network is fundamental for many practical applications in the field of smart roads. Continuously updated and reliable road network automatically extracted from Airborne Laser Scanning is an important aspect that may contribute to increase applications such as autonomous driving and infrastructure-to-vehicle communication.
The presented methodology has been tested in three different data sets provided by the ISPRS. The "Vaihingen" and "Toronto" Data sets concern two different involved urban environments. The former represents a typical European town and the latter a typical North American downtown with highrise buildings. The third data set, namely "Tobermory," represents a rural area with a large grassland and wooden areas.
The results for this three data sets were compared with OpenStreetMap (OSM) road centrelines. We can observe that road centrelines derived for urban environments ("Vaihingen" and "Toronto") are extremely similar to the OSM ones proving the reliability of the proposed method. In the case of the "Tobermory" data sets discrepancies are larger. This is probably caused by the inaccurate identification of road because points belonging to the road shoulder are erroneously classified as "road" points. This may affect the correct identification of the road centreline. In addition, significant fractions of roads are partially covered by trees, resulting in occlusions that contribute to worsen the final outputs. Further analysis on this issue are needed. On the other hand, it should be noted that a significant number of secondary dirty roads has been identified in the "Tobermory" Data set but they are not reported in OSM. Future work on road centreline extraction is aimed at combining ALS data with drone-based point clouds either derived from photogrammetry or LiDAR systems. Indeed, in those cases the possibility to combine a further information associated to point colour may provide a further parameter to enhance point classification and separation among different classes.