URBAN ROAD DETECTION IN AIRBORNE LASER SCANNING POINT CLOUD USING RANDOM FOREST ALGORITHM

The objective of this research is to detect points that describe a road surface in an unclassified point cloud of the airborne laser scanning (ALS). For this purpose we use the Random Forest learning algorithm. The proposed methodology consists of two stages: preparation of features and supervised point cloud classification. In this approach we consider ALS points, representing only the last echo. For these points RGB, intensity, the normal vectors, their mean values and the standard deviations are provided. Moreover, local and global height variations are taken into account as components of a feature vector. The feature vectors are calculated on a basis of the 3D Delaunay triangulation. The proposed methodology was tested on point clouds with the average point density of 12 pts/m that represent large urban scene. The significance level of 15% was set up for a decision tree of the learning algorithm. As a result of the Random Forest classification we received two subsets of ALS points. One of those groups represents points belonging to the road network. After the classification evaluation we achieved from 90% of the overall classification accuracy. Finally, the ALS points representing roads were merged and simplified into road network polylines using morphological operations. * Corresponding author


INTRODUCTION
The classification of airborne laser scanning point clouds aims typically at identification of ground points and points belonging to various objects above the topographic surface.Over the past years the high resolution ALS data have been also increasingly used for topographic object detection and modelling.
Information about road networks is a significant component of topographic data basis.It should be up to date because of its importance for many applications, e.g.incident and emergency responses, the transportation management, etc. Various data sources can be utilised to obtain the geometrical information about road networks.
Aerial imagery and the satellite imaging are well established data sources for the road network extraction (Shackelford et al., 2003), that is performed using solely two dimensional geometrical and radiometric information provided by the imagery.
The accurate and dense ALS point clouds can be explored in order to extract the road network for topographic data basis as well as for the road surface inspection.Different approaches have been proposed for extracting of a road network from ALS data (e.g.Hu et al., 2004, Clode et al, 2004, Choi et al., 2008, Ferraza et al, 2016).The approaches that rely solely on ALS data base on the creating of a digital surface model (DSM) or a digital terrain model (DTM) and then performing transformations and filtering on the reduced dataset.In order to filter the point cloud, Clode et al. (2004) use only intensity and height differences in DSM.Another approach is proposed by Choi et al. (2008).The authors combine approaches based on clustering using similar height and pulse intensity with additional constraints for height differences.Although the results are good, as authors point, the method is characterized by a high level of the false positive rate, especially on roofs of big, flat, and relatively low buildings (Choi et al., 2008).
In the paper of Zhao et al. (2011), a method is proposed, that bases on transformation of ALS data into binary image using point's intensity and height differences from DTM/DSM.Afterwards, the raster analysis is deployed to extract streets and streets junctions.The disadvantage of the method is the main assumption that the road network takes the form of a regular grid and most roads bend at right angle (Zhao, 2011).
A method that combines two data sources, ALS point clouds and imagery has been also developed (Hu, et al., 2004).ALS point cloud processing and filtering is supported in this method by an analysis of a high resolution aerial imagery.In described work authors use ALS data and imaginary as the separate datasets and analyze them separately in order to extract different features and then to combine the results (Hu et al., 2004).
The topic of the road detection using ALS data is mostly discussed in the context of urban areas, but there is also approach that addresses it to rural or forest areas.According to Ferraza et al. (2016), some additional assumptions have to be set up in this case.The authors of discussed work base their method on the classification of DTM utilising the random forest algorithm executed on the calculated geometry parameters of DTM triangles.Then, after filtering and processing the results, the road graph is generated, that provides very good representation of the road network, even on the areas with high canopy (Ferraza et al., 2016).This research addresses the identification of ALS points that are reflections on the road surface.The objective is to detect road points using ALS features and areal imagery features simultaneously.This novel approach is in the contrary to the work of (Hu et al., 2004), where these data sources were explored separately.We classify points into two subsets, road points and non-road points.No filtering is performed before the classification.In order to achieve the objective we use the random forest algorithm that is a classifier based on the decision trees.The feature vector required by this algorithm consists of components derived from ALS data and R, G, B taken from imagery.We consider ALS points that represent the last echo only.In order to simplify the road axes we use a set of morphological operations.The proposed approach is tested on the high urbanised city scene.

STUDY AREA
The test site with an area of 22 km 2 is located in Warsaw, the capitol of Poland.This area stands out for a large diversity of the road networks.There are single lane roads and multiple lane roads with median, alleys in residential areas, property access roads (driveways), shopping centre parking lots, roundabouts, viaducts.Moreover, parks, low and high buildings including large trading halls with flat roofs exist in the study area.Therefore, the study area is representative for an urban environment of a large city.The location of the study area on the city map is presented in Figure 1.

Figure 1. Location of the study area
The ALS data was acquired using the Riegl LMS-Q680 scanner.The average point density is about 12 points per square meter.Beside the coordinates, the following parameters were provided for each point: intensity, the echo number and the total number of echoes.Additionally, point cloud was coloured utilising orthoimage created from images acquired by the IGI DigiCAM60 camera.

METHODOLOGY
The extraction of the road network axes from an ALS point cloud is a complex issue, which requires several steps of the data processing.The main task is the separation of the point cloud into two sub datasets: road points and non-road points.To achieve this goal, pre-established constraints and filtering of the data transformed to the DSM are mostly used (e.g.Clode et al., 2004;Choi et al., 2008).
In this work the random forest classification is used to filter out points not belonging to roads.A number of point attributes is determined for this purpose.The result of the classification is then transformed to the binary raster, which is than reduced to the representation of the road axes using a set of the morphological algorithms.

Random Forests
Random forest algorithm is the statistical, robust method of homogeneous, large datasets classification.Random forests are the combination of tree predictors, where each tree depends on the values of a random vector sampled independently and with the same probability distribution for all trees in the forest.The generalization error of a forest of tree classifiers depends on the strength of the individual trees in the forest and the correlation between them.Internal estimates monitor: error, strength, and correlation.These parameters are used to show the response to increasing number of features used in the splitting.Moreover, the internal estimates are used to measure the variable importance (Breiman et al., 1984, Breiman, 2001).2013) integrated the random forest classifier into the conditional random field and applied this framework to classify ALS data in dense urban areas.The approach uses an extended set of contextual information for each point in the dataset.Each 3D point is classified in order to enable the detection of small objects as cars in the scene.Nevertheless, the approach shows a weakness in detection of trees.In this study we investigated this issue during road surface detection under tree canopies.

Training Data
The accuracy and correctness of a supervised classification depend strongly on the quality of training data.In this research we use 3 separate training fields, each representing different features that are commonly occurring in urban areas (Figure 2).The field 1 (570 m x 530 m) is an area covering the estate of tall apartment buildings, trees, small, single-lane roads and two-lane road separated by a median.The field 2 (570 m x 530 m) is a typical residential area.Roads occurring in this area are narrow, surrounded by trees and parked cars.In addition, there are gaps in the dataset due to water bodies.The field 3 (490 m x 470 m) is an area covering industrial district, where large and relatively low buildings occur.From the top view, those buildings are confusingly similar to parking lots.This training field should help to eliminate incorrect classification of building roofs as roads.
For the accuracy assessment the field 4 (570 m x 530 m) was selected.Similarly to other three test fields it was classified manually, but was used only for classification evaluation.It was not used as a training dataset.

Feature Calculation
The features selection is an essentials step in the data classification using random forest.Since in our approach 3D points are classified, the primary information available from the point cloud is mainly explored.The point cloud is initially prefiltered and only points that represent last echo are considered.To compose the feature vector, the information available directly from the point cloud data record (point intensity, number of echoes) and RGB are used.Furthermore, additional features are calculated for each point analysing geometrical configuration of surrounding points in the direct neighbourhood.The 3D Delaunay triangulation is built for this purpose and in each point P(x,y,h) the mean normal vector is calculated.
The extended vector of the features is composed of the following components: 1. intensity, 2. number of echoes, 3. Red, 4. Green, 5. Blue, 6. the mean value  of the deviation angles between the normal vectors and the z-axis, 7. the standard deviation of the deviation angles between the normal vector and the z-axis, 8. the mean value of the deviation angles between the normal vectors and the z-axis, calculated based on triangles within the 3 meter radius, 9. the standard deviation of the deviation angles between the normal vector and the z-axis, calculated based on triangles within the 3 meter radius, 10. information whether the point represents local minimum, 11. deviation between the mean height of the whole scene and the height of the considered point, 12. deviation between the mean height of the points within the 3 meter radius and the height of the considered point.
A given point P can be the vertex of several adjacent triangles of the 3D Delaunay triangulation.Each of the triangles has own normal vector n(x,y,h).The angle α between n(x,y,h) and z-axis is calculated for each triangle (Figure 3).Then the mean value of the α angles for adjacent triangles is calculated and associated with the point P.Moreover, the standard deviation for  is also calculated.This parameter can be considered as the terrain roughness in the point P.
In a similar manner the mean value  and its standard deviation are computed for the 3 meter neighbourhood.

Figure 3. Visualisation of the deviation angles between the normal vectors and the z-axis
The component "local minimum" of the feature vector is calculated utilizing points that are vertices of the adjacent triangles.This value is "true" for the point P if ) min( i p h p h  , i=1,2,..,n.
To calculate the 11 th component of the feature vector the mean value of h for the whole scene is calculated.In order to remove ghost points, only points within the interval <h min , h max > are taken into account: (2) where σ is the standard deviation of the heights on the whole scene.Alternatively to this approach, the M-estimation can be applied to determine optimal mean plane for the scene (Chechata et al., 2009).Finally, the height differences between h and heights of the ALS points are calculated.In a similar manner the local values, within the 3 meter neighbourhood, are determined.

Feature Evaluation and Accuracy Assessment
In order to build classifier, the data representing all training fields were used simultaneously.From each training data set 100 000 points representing both road and no-road classes were taken randomly.In the first pass 100 trees were grown utilising all components of the feature vector.This number of trees is sufficient to stabilise the mean squared error of generalisation (Figure 4).As a result of the first pass calculation, the variable importance of the particular features (Figure 5) and the generalisation error of 7% were estimated.The plot (Figure 5) shows the increase in mean squared error averaged over all trees and divided by their standard deviations, calculated for each variable separately.
Figure 5. Feature importance Based on the results of the classifier training, the vector of the features was reduced.Due to low value of the importance of the 8 th variable, the mean value of the deviation angles within 3 meter radius was removed from the feature vector.For the reduced vector of the features we received the generalization error of 7%.Therefore, the reduced set of features is as well effective as the original vector of features.The classifier built in the second pass was used in the further investigations.
All data representing the test fields were classified using constructed classifier.The level of certainty was set up on 85 %.It means that a point was classified as the road point if 85% of the trees voted for that.The results of classification were compared with the results of the manual classification.This comparison is given as the confusion matrix for the particular training field (Table 1, 2, 3) and for the evaluation field (

Road axes detection
As a result of the classification the point cloud is divided into two subsets: road points and non-road points.Since the level of certainty was set up at 85 % all points below this value have been classified as non-road points.The subset "road points" have been then converted to raster with 0.5 m pixel size.In order to achieve road axes, a set of binary morphological operations, implemented in Matlab, is utilised in the following order: majority, clean, thin, clean, diag.

RESULTS
After the classification of the test area, the filtered point cloud representing roads was obtained.The resulting dataset contains only those points from source dataset that were assigned to class roads with 85% level of certainty.The classified points with the assigned level of certainty are shown in Figure 6 for a representative part of the study area.After the generation of the binary image and applying the morphological algorithms, the road axes were produced.The morphological operations were able to reduce small gaps (Figure 7).Nevertheless, the big gaps caused by overpasses still remain (ellipses in the Figure 7).Besides this issue, the overall accuracy and quality of the road detection is satisfactory.
Figure 7. Produced road axes (a part of the study area)

CONCLUSION
In this study, the random forests were successfully applied to detect the urban road in the ALS point cloud.ALS and imagery based parameters were selected as the components of the feature vector.After the feature evaluation by random forests, the local variation of the mean value of the normal vector (feature No 8) has been removed due to low importance.The most relevant features are intensity, number of echoes and height differences.
The RGB features improve the classification in the presence of vegetation.We received the overall classification accuracy of 93%, 96 and 99% for three training datasets, respectively.For the evaluation dataset we received the overall accuracy of 90%.
Random forests can be utilised for the ALS point cloud classification, based on various point attributes, such as heightbased, echo-based, plane-based or full waveform-based.According to Chehata et al. (2009), random forest algorithm can classify ALS data with the accuracy exceeding 90 %.Niemeyer et al. (

Figure 2 .
Figure 2. Training fields: top left: field 1, top right: field 2, bottom left: field 3. Evaluation filed, bottom right -field 4. Blue colour represents road points of ALS point clouds.

Figure 4 .
Figure 4. Mean squared error as a function of the grown tree number

Figure 6 .
Figure 6.Classification results (a part of the study area) Some gaps are present in the data set representing road.These gaps are produced by cars and road overpasses.Unfortunately, overpasses were not present in the training datasets.The detailed analysis of the results shown that involving the RGB in the vector of features improves the classification in the case of the presence of vegetation in the direct neighbourhood of the road.

Table 4 .
Confusion matrix for the field 4The values on the diagonal of the matrix represent values of the correct classification.The off-diagonal values are the numbers of points that were classified incorrectly (e.g.Table1, values 245349 and 84668).Moreover, for the accuracy assessment, the following parameters were estimated: I CP -Overall classification accuracy Kappa -kappa coefficient, overall classification agreement E c1 -Error of commission for class non-roads E c2 -Error of commission for class roads E o1 -Error of omission for class non-roads E o2 -Error of omission for class roads A p1 -producer's accuracy for class non-roads A p2 -producer's accuracy for class roads A u1 -user's accuracy for class non-roads A u2 -user's accuracy for class roads These values for the particular test fields are given in Table6.Since the datasets of the test fields 1, 2 and 3 were used simultaneously for the classifier construction, the values given in the columns 2, 3 and 4 of Table6reflect the internal quality of the classifier.The values collected in the column 5 (Field 4) show the classification accuracy on the evaluation dataset, which can be interpreted as the absolute accuracy.

Table 6 .
Accuracy ratings for training and evaluation fields