ROAD SURFACE DETECTION FROM MOBILE LIDAR DATA

The accurate three-dimensional road surface information is highly useful for health assessment and maintenance of roads. It is basic information for further analysis in several applications including road surface settlement, pavement condition assessment and slope collapse. Mobile LiDAR system (MLS) is frequently used now a days to collect detail road surface and its surrounding information in terms three-dimensional (3D) point cloud. Extraction of road surface from volumetric point cloud data is still in infancy stage because of heavy data processing requirement and the complexity in the road environment. The extraction of roads especially rural road, where road-curb is not present is very tedious job especially in Indian roadway settings. Only a few studies are available, and none for Indian roads, in the literature for rural road detection. The limitations of existing studies are in terms of their lower accuracy, very slow speed of data processing and detection of other objects having similar characteristics as the road surface. A fast and accurate method is proposed for LiDAR data points of road surface detection, keeping in mind the essence of road surface extraction especially for Indian rural roads. The Mobile LiDAR data in XYZI format is used as input in the proposed method. First square gridding is performed and ground points are roughly extracted. Then planar surface detection using mathematical framework of principal component analysis (PCA) is performed and further road surface points are detected using similarity in intensity and height difference of road surface pointe in their neighbourhood. A case study was performed on the MLS data points captured along wide-street (two-lane road without curb) of 156 m length along rural roadway site in the outskirt of Bengaluru city (South-West of India). The proposed algorithm was implemented on the MLS data of test site and its performance was evaluated it terms of recall, precision and overall accuracy that were 95.27%, 98.85% and 94.23%, respectively. The algorithm was found computationally time efficient. A 7.6 million MLS data points of size 27.1MB from test site were processed in 24 minutes using the available computational resources. The proposed method is found to work even for worst case scenarios, i.e., complex road environments and rural roads, where road boundary is not clear and generally merged with road-side features.


INTRODUCTION
Roads are often the single largest publicly owned national asset.Road network reduces the distance between people, markets, services and knowledge, which facilitate transport, trade, social integration and economic development.Road network is growing rapidly in order to accommodate the increasing traffic load.Therefore, effective management and planning of road network is essential.This requires comprehensive and accurate information (Kavzoglu et al., 2009), i.e. geometric and radiometric details of road including conditions of road.Road surface threedimensional (3D) information collection is among the most important steps for pavement condition assessment.Information about the road and its surface is becoming important for the road maintenance (Mc Elhinney et al., 2010;Yang et al., 2013).The assessment of road for maintenance and safety measures is highly essential on periodical basis.The in-situ measurements along with visual examination and interpretation are traditionally utilized by many for road surface evaluation (Schnebele et al., 2015).In developing nations like India the road environment is quite complex and heterogeneous in nature.Therefore, collection and accurate processing of roadway data along tens of thousands kilometre long roadway is a big challenge due to heterogeneous surrounding.Therefore the researchers have been documenting the road information using different types of remote sensing data collected over the corridor since past few decades.Satellite images, aerial photographs and point cloud data are main data sources.Most of the proposed methods (Ferchichi and Wang, 2005;Wan et al., 2007;Mokhtarzade and Valadan Zoej, 2007;Mena and Malpica, 2005;Mohammadzadeh et al., 2006;Wang et al., 2005) using satellite images and aerial photographs only provide road pixels and its twodimensional (2D) location information.Accurate 3D road information is possible using point cloud data captured by airborne (Jiangui and Guang, 2011) and mobile LiDAR system (MLS).High-precision road terrain obtained in terms of point cloud data is the basic data for further analysis for applications such as road surface settlement, pavement, and slope collapse (Liu et al., 2013).MLS is gaining popularity in 3D LiDAR mapping applications along various corridors because of the extreme ease in capturing comprehensive high resolution 3D topographic data at highway speed (Yadav et al., 2014).A MLS observes roadway and nearby objects in the form of their dense coordinates at high speed (Yen et al., 2011) which makes it possible to construct a detailed 3D model of the road environment (Lehtomäki et al., 2010).A MLS achieve the highest data quality and completeness among the traditional roadway data collection methods (Jalayer et al., 2014).MLS data are used for accurate road geometry computation over existing conventional methods due to significant savings in time and human resources while yielding high accuracy.MLS has advantage over traditional roadway data collection methods (Puente et al., 2013) in terms of high speed data capture (time and cost reduction), high density of collected point cloud, comprehensive topographic survey and minimization of the erroneous/ questionable data, remote acquisition and measurement (increases survey efficiency and safety), and deliverable being coloured 3D point cloud that accurately represent roadway objects in the scene.Only a few researchers have tried to extract the road surface using 3D point cloud data captured along the corridor using MLS.Available literature on road surface LiDAR points extraction are categorized into two parts: (I) road curb based road surface detection, (II) road surface detection using global property of roads (e.g., topology and smoothness) and uniform road LiDAR points intensity value.The methods proposed in the first category use curb information, where height jump is noticeable along the scan line.In (Ibrahim and Lichti, 2012) a method to extract the road surface LiDAR points using different density datasets having straight and curved curb is proposed.Point cloud segmentation into ground and non-ground is performed.Refinement of ground segment using morphological neighbourhood property is done followed by edge detection for curb and constructing closed polygon to extract the street floor.Topology and smoothness of road and the local shape features of point clouds are used in (Yang et al., 2013) for detecting and tracking the curb.Point clouds are partitioned into road cross sections using GPS time of points and moving window operator is used for extracting candidate road areas.Curb points are detected by geometrically analysing local points.Angular distance to ground normal map is computed in (Hervieu and Soheilian, 2013) using LiDAR data.It is used for generating feature map and distinguishing 3D curb patterns from other features, such as ground, facades or stairs of higher heights.Road edges are detected by applying prediction/estimation model on feature map.Two stages algorithm for extracting the road edge from LIDAR and navigation data is proposed in (Mc Elhinney et al., 2010).The first stage of this algorithm creates a set of road cross sections.In the second stage these cross sections are processed into 2D lines.These lines are then analysed based on the slope, returned intensity, returned pulse width and proximity to the vehicle to determine the road edges.In (Serna and Marcotegui, 2014) a method for automatic analysis of curb accessibility from MLS data is proposed.First, input point cloud is mapped to range images.Second, the image is interpolated in order to avoid connectivity problems and the quasi-flat zones algorithm is used to segment the ground (road+sidewalk).Then, curb candidates are selected using height and elongation criteria, and close curbs are reconnected using Be źier curves.Three-dimensional virtual grid, multi-scale neighborhood iterative analysis and local slope filtering ideas are applied in (Liu et al., 2013).They have extracted expressway surface points from mobile laser scanning point clouds.In (Zhang, 2010) a method to identify road regions and road-edges using LiDAR range data is discussed.The road segment and road-edge points are first identified using the elevation information extracted from the range data.Input point cloud is first separated into sequences of points that represent scan lines in (Miyazaki et al., 2014).Further a line-based region growing method is applied in order to detect planar structures with precise boundary from point clouds with uneven distribution density of points.In category II, a rough classification of point cloud into ground and non-ground is executed (Lien et al., 2012).Cubic curve least squares fitting using lowest point, second lowest point away from it by road range and intermediate points is performed for selecting the road points.Major limitation of existing methods is to extract rural road surface which are mostly without curb.It is difficult for the method proposed in (Yang et al., 2013) to deal with curbs with boundaries that are characterized as asphalt/soil, asphalt/vegetation, or asphalt/grassy bank.It is difficult to mark exact edge location, if point density is less and sufficient elevation jump is not present at boundary of carriageway (Mc Elhinney et al., 2010).Method proposed in (Wu et al., 2013) is dependent on the road range, which cannot be calculated correctly if height difference between the road surface and the shoulder is not obvious.Computation time efficiency is not discussed in the existing methods for extracting the 3D road surface points.It is very important parameter and need to be addressed, because data sizes of gigabytes/km are generally available from these MLS surveys.The objective of this paper is to develop a computationally time efficient method for extracting road surface without curb.Further objective is to eliminate the limitations of existing methods as described in the previous paragraph.A test site from rural locality is selected for algorithm testing.

METHODOLOGY
In the proposed method first square gridding and rough ground classification are performed and then road surface points are detected.The first stage divides the input MLS data points into square grid of predefined size.Then planar ground points are extracted from each grid.These planar ground points also include road points.Finally road surface points are detected from the planar ground points using thresholding criteria on parameters computed by principal component analysis (PCA), 2D point density, range of intensity values and intensity standard deviation.

2.1
Square gridding  (Yadav et al., 2015).Square grid size (m×m) is defined by the user on the basis of the slope of the ground within the grid.Square gridding is performed on the projected MLS data points on the XY plane.X-axis and Y-axis serve as grid lines.
Figure 1.Rough ground classification by clustering.Points lie within range of minimum Z to min Z h .

Rough ground classification
Rough ground classification is performed on the basis of grouping of data points lying in the specific range of Z values starting from minimum Z ( ).Initially MLS data points of each grid are rearranged in increasing order of their Z values.Points lying within and + ℎ are selected.These points are from the first slice of height h (Figure 1).The first slice includes mostly the LiDAR points that belong to the ground surface and bottom of the other non-ground objects (Wu et al., 2013).Therefore, MLS data points within the first slice contain road surface points.The effect of outliers is avoided by considering the average of 100 data points from top of rearranged data file as minimum Z.It is verified through several preliminary experiments

Road surface points detection
Road surface is piecewise planar structure generally parallel to XY plane (see Figure 2).It is detected in this section.The proposed method uses ground points extracted in the previous section.Therefore, orientation of planar surface is also checked.The angles between and Z axis, and between and Z axis are calculated.The planar surfaces nearly parallel to XY plane are detected by applying thresholds conditions on the previously computed parameters using data set . The same criteria is generalized and used for each disc data set , where represent arbitrary disc (see Equation 1).Total numbers of discs are represented by .Where and are threshold values.represents planar ground surfaces nearly parallel to XY plane and includes non-planar ground surfaces and planar ground surfaces not nearly parallel to XY plane.contains the road surface data points, which are extracted in the following section.Road surface extraction method proposed in this section is based on the geometrical property and radiometric response of road surface.Roads are manmade structures having ribbon like nearly flat pattern.MLS data points from road surface have quite less height variation in their neighborhood, that is, less height standard deviation (see figure 4(a)&(c)).Also, the road MLS data points have nearly similar intensities (I) within narrow range depending on the type of road materials and frequency of pulse energy used in MTLS.Intensity values of MTLS data points have good separability among various objects classes, if the wavelength of the pulse laser is suitable for ground materials.The relative separations between ground features, i.e. asphalt road, grass, building and tree have been compared using intensity data (Hu et al., 2004).It is found that the separability is very high for road vs. grass and road vs. tree (Song et al., 2002).The planar ground points extracted in the previous section contain road points from actual road surface and other non-road planar ground points.The 2D point density of planar ground points decreases away from MLS ground projected trajectory (see Figure 5).The road points have uniformity in intensity and height values in their surrounding (see Figure 4  Road points set is generated using Equation 2, where , and are thresholds for 2D point density and standard deviation of intensity, respectively.and are the lower and upper thresholds of intensity values, which belongs to the road.Above process is initially applied on the seed disc to check its candidature as part of road surface.Further, new seed point is generated by shifting 2D coordinate of the previous seed point along the X-direction by the same user defined radius R. Again the above whole process is repeated.MLS data points are covered by horizontal circle growing along the X and Y directions and road surface points are detected.

Test site and mapping data
The test site was located along rural roadway site in the outskirt of Bengaluru city (South-West of India).It was wide-street (two-lane road without curb) of 156 m length (see Figure 3(a)).Test site was surveyed using StreetMapper 360 (StreetMapper 360, 2014) MLS and 3D point cloud data of site was acquired.StreetMapper 360 is designed for the rapid mapping of highways, infrastructure and buildings (Yadav et al., 2017).It delivers proven accuracies in the most challenging environments.The typical positional accuracy is better than 2 cm (1×σ) and the point-to-point accuracy within the data is 1 cm (1×σ), where σ represents standard deviation.A total of 7,600,988 (of size 27.1MB) data points were acquired from the test site at 2D point density (Pts/m ) of 708.Test site included road-side homogenous trees, bushes and shrubs, utility poles and single storey buildings (see Figure 6(a)).This site was specifically chosen to test the proposed method effectively because the road boundary

Reference data
Reference data for road surface points was collected by visual inspection of colored point cloud data (X-Y-Z-R-G-B format) of test site.Collection was based on the visual identification of road boundary in the colored point cloud data.Once road boundary was identified, all the points within it were segmented out from the input data.Segmented out road surface had total number of 1,866,066 points, which was used as reference data points for validating the result.Average road width was calculated by dividing the road into 5 equal segments.Road width at each 5 segmented location was measured manually and averaged out.It came out to be 7 m.It was used as reference average road width of the test site.

Result
The proposed method was tested on the MLS data points of test site using the threshold values of parameters (see Table 1) the road surface points (see Figure 7) were extracted.First, square gridding and rough ground classification were performed using user defined thresholds listed in Table 1.
Their values were chosen based on the terrain undulation.Road surface behaves like piece-wise planar connected surfaces.The 2D point density of ground points is highest near ground projected MLS trajectory and it decreases away from the trajectory.Large regions formed by considering similarity in intensity and elevation (height) are mainly lies on the road surface and have large number of data points (see Figure 5).In intensity based road surface filtering a range of intensity values from road surface was chosen by analyzing intensity values of different road patches.These road patches were selected from different locations of the road, i.e. along and across the road.Road surface is like ribbon pattern and made up with specific material, so the intensity differences of neighboring points are quite less.These above properties were used to decide threshold values of road surface point through conducting several preliminary experiments.Total extracted road surface points for the chosen test site were 1,798,531.Road boundary of the test site was not defined clearly.It was separated by asphalt/soil, asphalt/vegetation, or asphalt/grassy bank (Yadav and Singh, 2018).Thus, extracted road surface boundary was not smooth (Figure 7).

Discussion
Proposed method uses only XYZ coordinate and intensity I of MLS data points, thus it does not depend on the scanning geometry and neighborhood structure in the data file.

Quantitative evaluation of 3D road extraction
Quantitative evaluation of the extracted road surface by proposed method is shown in

Run time performance
The execution time of the proposed method mainly depends on the parameters: (1) square grid size (m×m), (2) height of first slice (h), (3) radius of circle (R) and point density of dataset ( ). Results were generated on the tuned values of these parameters for getting highest precision, recall, and overall accuracy.Proposed method was implemented on the test site, 7.6 million MLS data points of size 27.1MB, using Matlab2017a installed in standard PC (CPU: Intel Core i5@3.6GHz 8 th Gen., RAM: 16GB).The total execution time of the proposed method was 24 minutes.Discussion on execution time for 3D road surface point's extraction from MLS data points is not clearly emphasised in the existing literature.So time complexity comparative analysis of proposed method with existing methods cannot perform here.

Conclusion and Future Recommendation
A computationally time efficient method is proposed in this paper for extracting road surface from MLS data points.3D road surface point extraction is an important step in the pavement condition assessment.First rough ground classification is performed then planar surfaces having road surface are extracted using mathematical framework of PCA.Further road surface points are extracted based on the similarity in intensity and height difference of road surface pointe in their neighbourhood.Point density achieved in the roughly classified ground points is also highest from road surface thus used as additional criterion to extract road surface points.MLS data points from complex road environment test site were used for testing the proposed method.Precision of 98.85% and recall of 95.27% were achieved.Overall accuracy is 94.23%.Average road width was measured at accuracy of 1.43%.Missing and falsely classified road points were from the road boundary, because the boundary of the road was not clearly geometrically defined.Using the results obtained in the test site and thorough analysis of the existing methods, it can be concluded that the proposed method is more general and removes the limitations of the existing methods by: (i) detecting road surface points without curb (ii) being computationally efficient (iii) being independent of the scanning geometry information and requiring only MLS data points in X-Y-Z-I format.Future work would focus on smoothing of road boundary.
First the 2D version (only X and Y values) of set is organized into partially overlapping circles having fixed radius R. 2D ground points within these partially overlapping circles take the shape of partially overlapping discs of different heights once their elevation values are considered (see Figure3).Mathematical framework of PCA is applied in each disc to extract planar surfaces parallel to XY plane.Initially, 3D coordinates of MLS data points from each disc are selected and the variance-covariance matrix C is computed, that is further used to compute three Eigen values, that is, , and , where ≥ ≥(Yadav et al., 2016).The normalized Eigen values , and are also computed.The Eigen vectors and along maximum variance ( ) and second maximum variance ( ) are computed.The and tell about the extent and direction of maximum spread of MLS data points.Similarly and , and indicate second maximum and minimum, respectively extent of points spread and their direction.In the normalized Eigen values, two of them are nearly equal (0.4 ≤ ≤ 0.6; 0.4 ≤ ≤ 0.6 ) and third one is very small (0 ≤ ≤ 0.1) for planar surfaces (see Figure3).The road surface is associated with the planar surfaces, which are nearly parallel to the XY plane.

Figure 2 .
Figure 2. Point cloud spread and shapes of road surface and other dominant objects present in the MLS data.

Figure 3 .
Figure 3. Eigen values and Eigen vectors, which represent road surface MLS data geometrically.

Figure 4 .
Figure 4. (a) Perspective view of MTLS data cross section across the MTLS path (b) intensity variation of neighbors of road surface, (c) height variation of neighbors of road surface.

Figure 5 .
Figure 5. Showing (a) width of road surface; (b) ground point density decreases across the MLS path; (c) heights and intensity roughness parameters are uniform, respectively on the road surface.
(b) &(c)).These information are used here in designing of a method to extract road points from .Again planar ground points are divided into partially overlapping circles of radius R which take shape of discs once elevation/height values are considered.The 2D point density ( ), average and standard deviation ( ) of intensity values are computed in each disc.
was not structured in the form of raised curbs (see Figure6(b)&(c)) and it was only delineated by non-uniformly distributed shrubs and asphalt/soil bank unlike datasets used in the literature.Even at many locations of road boundary the road and road-side surfaces were not geometrically differentiated (see Figure6(c)) and acted like levelled surfaces.

Figure 6 .
Figure 6.Perspective view of MLS data points: (a) along the roadway of 156 m of length (b) closure view showing the road coverage (c) photograph showing the characteristics of road surface at its boundary, that is, levelled surface.

Figure 7 .
Figure 7. Perspective view of MLS data points of extracted road surface format are used as input, where i and r represent an arbitrary data point and the total number of data points, respectively.XYZ is the coordinate of MLS data point and I represents the intensity value.The input data points are projected onto XY plane.Rectangular minimum bounding box is generated for closing the 2D projected data points.Bounding box is defined by its four vertices.Further, the data points within bounding box are divided into 2D regular square grids

Table 1 .
Parameters values used for the test site.
th I 4

Table 2
TP and TN are the true positive and true negative classified road points; FP and FN are the false positive and false negative classified road points.