In this paper, we proposed a road extraction algorithm that utilises voxel-based skewness and kurtosis balancing in order to update existing road networks accurately and reliably from airborne lidar data. The proposed algorithm consists of initial road extraction followed by a series of refinements including width constraints, spatial interpolations and curve fitting which are successfully implemented In: B. Veenendaal and A. Kealy (Eds.): Research@Locate'15, Brisbane, Australia, 10-12 March 2015, published at http://ceur-ws.org
reduce
join
up
misclosures
and
determine centrelines of the extracted roads. A
numerous sample lidar datasets are tested in
order to assess the proposed algorithm. The
quantified completeness and correctness of the
classification results are 98.20% and 98.54%,
respectively; hence it is concluded that the
proposed
algorithm
works
effectively
in
acquisition of road features from airborne lidar
data.
Keeping road network databases up-to-date is important for many applications e.g. emergency handling,
car navigation, tourism, traffic management and monitoring, intelligent transportation systems, web
map services, location-based services, because of rapid urbanisation in the developing and developed
countries
Height differences, normal vectors, and contextual information are the most widely used
geo-features to extract roads from lidar data. In general, typical feature extraction strategies consist of
segmentation and clustering
Akel et al. (2005) applied a region growing segmentation method based on elevation and normal
vectors to detect road areas.
Because of the shadow effect in high-density vegetation areas, a good level of completeness of extracted roads is difficult to achieve. For the purpose of higher accuracy and completeness, some researchers have utilised aerial images in addition to lidar data. Hu et al. (2004) utilised aerial imagery and lidar data together in combination to extract road points. Their results show that the accuracy and completeness are improved. Zhu et al. (2004) introduced the associated road images extracted from lidar points in order to enhance the results from real road images from aerial mapping. However, road extraction from the fusion of optical images and lidar points is still challenging since the co-registration of the images and the lidar points may become a noticeable error source.
Existing geospatial data can provide a good complement to lidar data as it can be served as a priori
knowledge about roads
This study aims to utilise airborne lidar data and the national infrastructure database in order to update the existing road networks. We employ a voxel-based skewness and kurtosis balancing algorithm which is based on the assumption that ground points exhibit a normal distribution in order to filter out object points from the lidar data. Since multi-classes and various road materials can distort the normal distribution, the existing road database is utilised to segment the lidar point clouds into buffer-like tiles so that the assumption can hold. One of the advantages in doing so is that only a few parameters i.e. voxel size and buffer size are necessary. Also, this method can deal with different types of roads regardless of the road materials (e.g. asphalt or concrete) if the candidate roads in the same tile are made of a single material. It should be noted that not only shadows of objects e.g. trees and buildings but also open areas e.g. parking lots may affect the road extraction results, hence further refinements are essential to improve the accuracy and completeness, once initial extraction is performed. A series of refinements including width constraints, spatial interpolations and curve fitting are implemented in the proposed method in order to reduce false-positives, join up misclosures and determine centrelines of the extracted roads.
Skewness and kurtosis balancing
One of the methods to assess the asymmetry of data distribution is skewness, denoted by sk
sk
1 N
N 3 i1 zi 3 ku 1 N
zi 4 N 4 i1 (1) (2)
Compared with other algorithms such as Expectation-Maximization (EM)
Intensity and elevation are common attributes used in skewness and kurtosis balancing. To test the effectiveness of the two attributes, we apply elevation, intensity, and a combination of the two attributes to a sample data, respectively. The results are shown in Figures 1-3. As seen in Figure 1, the elevation values of some road points in the blue rectangle exceed those of vegetation in the black rectangle because of the slope of the road. As a result, road points in the blue rectangle can be eliminated while the vegetation points in the black rectangle can remain. Hence the elevation-based skewness and kurtosis balancing algorithm (SKBA) may not perform well in mountainous areas. It is observed in Figure 2 that the intensity-based SKBA can remove most of object points, however, the intensity values of some trees may also satisfy the normal distribution. While SKBA using the combination of elevation and intensity shows a reasonable accuracy, a large amount of highly elevated road points can be omitted. Taken all the factors into account, the intensity-based SKBA is applied in this paper. 3.
A framework of the proposed road network update scheme is illustrated in Figure 4. The first step is to segment lidar data according to the buffers generated by the road network. The second step is to partition the buffered points according to geographic coordinates and reorder the points along the main direction of the existing roads. Voxels are created on the basis of geographic coordinates in order to allocate the reordered points to different voxels and remove object points. Once the object points are removed, the intensity-based SKBA is applied to the remaining points. Some of non-road features such as parking spaces can be possibly classified as roads, and the misclosures among road points may occur because of the blockage by tall trees and buildings. Hence refinements including width constraints, spatial interpolations and fitting curves are applied to the initial extraction results.
Skewness and kurtosis balancing does not consider the material characteristics of the road, therefore
different road materials within a road network can disrupt the normal distribution of the intensity
values of the whole dataset. It seems that this is the main reason why the researchers have not used
skewness and kurtosis balancing for road extraction purposes. Since skewness and kurtosis balancing is
largely reliant on the assumption of a normal distribution, preprocessing is required to make the
assumption hold.
Airborne lidar data typically consists of a huge amount of points. Therefore processing the entire points without proper data management is time-consuming. Also, iterations in skewness and kurtosis balancing will take time if the whole dataset is processed at once. To reduce the computation time effectively, our approach is to partition lidar points along the main direction of road segments according to the geographic coordinates. The core steps of partitioning of lidar points are explained as follows: a)
Calculate the coordinate difference between the starting node and the ending node of the line segment in x-direction and y-direction, respectively. b) If the absolute value of the coordinate difference between the starting node and the ending node in x-direction is larger than that in y-direction, set x-direction as the main partitioning direction; otherwise set y-direction as the main partitioning direction. c) d)
Search the point with the minimum coordinate value in the main partitioning direction. For each point, calculate the index number according to the Equation (3):
I INT[(C Cmin ) S] (3) where I is the index number of each point, C is the coordinate value of a point in the main partitioning direction, Cmin is the minimum coordinate value in the main partitioning direction, and S is the step value which is equal to the absolute value of the slope. INT is a mathematical operation that rounds the result to the nearest integer. e)
Cluster the points with the same index number. The partitioned points are indexed for the purpose of efficient data management.
After the voxel-based filtering, the intensity-based SKBA is applied to the remaining points to filter out object points. In fact, points of the maximum intensity value are removed before applying the intensity-based SKBA since the intensity of road points is overall the lowest in the sample datasets. The process can be described as follows: (i) Calculate the skewness and kurtosis values of the dataset. (ii) For a normal distribution, the kurtosis value is 3 and the skewness value is 0. If the absolute skewness is larger than the threshold and the absolute difference between kurtosis and 3 is larger than the threshold, go to step (iii); otherwise, go to step (vi). (iii) The point with the maximum intensity value is regarded as the object point and is removed. (iv) Recalculate the skewness and kurtosis values of the dataset. (v) Repeat steps (ii)-(iv) until all the points are processed. (vi) If there are some remaining points, label them as candidates for road points.
Once the initial road extraction is performed, some open areas can be shown as roads since they present similar characteristics as roads. Also, gaps between road segments may occur because of the blockage by trees, buildings and/or vehicles. Therefore, further refinements are needed to improve the results. The refinement process consists of width constraints, interpolations and fitting a curve to the road points. At first, road sections i.e. cross-sectional points perpendicular to the main direction of roads are generated. Then the width of the road is calculated based on the mean values since the widths of roads within a road network may vary. Then the road section whose width is larger than that of the road segment i.e. a part of the road between two nodes is processed. The points far away from the centre point e.g. points of more than a half of the road width away, are removed. After the false positives are removed, a local inverse distance weighting (IDW) interpolation is applied to connect the misclosures caused by the blockage of trees, buildings, or vehicles along the main direction of the road. The local linear interpolation is used to minimise the bias caused by the curves along the road. Connectivity is enforced to join up the misclosures between two different road segments. In case of the jagged edges, symmetric points to the road network are calculated. To obtain a smooth road surface, a curve fitting is adopted, and the centrelines of roads are calculated. These processes are shown in Figure 5. Because of the buffer distance along and perpendicular to the main direction of the road, the intersections between adjacent roads will appear, especially in conjunction areas. To refine the topology, we force the connecting roads end in the intersection nodes if two roads intersect. In terms of multiple intersection nodes, the mean intersection node is calculated according to the coordinates and we force the connecting roads end in the mean intersect nodes. (a) Road Candidates
(b) Deletion of False Positives (c) Interpolation (d) Curve Fitting
To validate the proposed method, tests with various scenes and lidar data samples are conducted. The first sample is a part of Anzac Parade which is an arterial road and surrounding neighbours near the University of New South Wales (UNSW), consisting of 128,847 points with a point density of 0.94 points/m2. The surroundings of the dataset are complicated as they contain both high-rise buildings and small residential houses. In addition, different types of roads are present in the sample, including an arterial road, a local lane, an avenue and a local street. Also, Anzac Parade consists of two parallel roads with a road separator in the middle. As the study area is near a university campus, many roads are highly occluded by cars. The second sample is a part of Botany Street which is another arterial road and its surrounding residential areas near UNSW, including 69,064 points. Apart from roads, this sample dataset mainly consists of small houses and dense vegetation. The main challenge for this 2 dataset is its sparse point density of 0.74 pt/m . The third sample is a residential area near UNSW, including 144,255 points. In this dataset, the high slopes along the streets cause huge elevation differences up to 50.73 m, which is a great challenge to distinguish ground points from building points. Another problem which may affect the completeness of the result is that roads are heavily blocked by trees in some areas. The fourth sample is a part of Barker Street which is a main road and its surrounding areas near UNSW, consisting of 124,690 points. Residential houses, trees and roads are main features in the dataset and the challenge of this sample is a high slope along the street resulting in the fact that the elevation values of some houses are even lower than those of some roads.
Based on the proposed algorithm, the voxel sizes in x, y, z directions are set to 10 m, 10 m, and 1.5 m, respectively. If the voxel size in x- or y-direction is too small, some building rooftops may split into different voxels, which shortens the height difference between adjacent neighborhoods, therefore misclassifies building rooftops as part of ground points. If the voxel size is too large, some road points in mountainous areas will be removed as well. Figure 7 shows the overlay of the centrelines of the four sample datasets and the aerial image. Figure 8 shows the road extraction results and corresponding centrelines of the four sample datasets, which indicate the feasibility of the proposed algorithm.
Pr
1i01 Di2 10 Ps rN1 Lr Pr2 rN1 Lr (4) (5) where Pr is the positional accuracy of a road segment r, Di is the distance between the ith node of the reference and the extracted road, Ps is the positional accuracy of a sample dataset, N is the total number of road segments in one sample dataset, r=1,2,3, ….N, and Lr is the length of the rth road segment.
Sample 1 Sample 2 Sample 3 Sample 4
Our quantitative analysis results are summarized in Table 2. Although the majority of the extracted roads match with the actual roads, some are slightly deviated from the reference centrelines, some lie on the edge of the roads. The main reason for the deviation is that the elevation change along the road is large, which can lead to the removal of ground points in highly elevated areas. Since there are no points in these areas to indicate which direction to interpolate, an error may occur and the situation becomes worse when the width of the road is large. Although symmetric points corresponding to the road network are calculated to avoid such a situation, existing positional errors in the road networks may result in the positional errors of the symmetric points. Another vulnerable area is the traffic islands. As shown in the extraction results for Samples 1, 3, 4, the approach fails to obtain the shapes of traffic islands. While Sample 1 has the highest correctness, completeness and quality values, the positional accuracy of the test dataset is the lowest among all samples. The unweighted positional accuracy of all the samples are calculated which results in 6.60 m, 4.62 m, 3.78 m and 6.37 m, respectively. This indicates that the road segments in Sample 1 are deviated from the reference centrelines within the tolerance difference. (a) Voxel size of 2 m (b) Voxel size of 8 m (c) Voxel size of 10 m (d) Voxel size of 20 m The proposed method may not perform well on the arcs of traffic islands. The method also yields poor results if road networks are not used because of multi-classes in lidar data. If roads in one buffer consist of various materials, omission error of road extraction may increase. The voxel size in x- or y-direction is crucial in filtering object points in mountainous areas.
Updating road networks is critical in many applications e.g. traffic management, traffic monitoring etc. In this paper, we presented a novel method to update existing road databases with airborne lidar data. The proposed method takes advantage of a voxel-based skewness and kurtosis balancing to distinguish road points from the lidar point clouds. To make the intensity of roads suitable to the proposed SKBA, road networks are utilised to divide the lidar points into different tiles since different types of roads within a road network will distort the normal distribution of intensity values of road points. The lidar points are reordered along the main direction of roads and partitioned into sections according to the coordinates. Voxel-based filtering is applied to filter out remaining object points. Skewness and kurtosis balancing is used to classify the candidate road points. In order to refine the results, width-constraints are applied to eliminate false positives, and misclosures within a road segment are filled by a local interpolation method. Calculating road width automatically based on mean values in the refinement step avoids the problem of different road widths within a road network. The proposed method extracts roads from airborne lidar data with an acceptable completeness and correctness. The quantitative results show that it is a feasible method for updating the road databases regardless of the road types. There are some biases in the extracted centrelines that can be larger than the threshold. The main reason for this is that the removal of ground points in highly elevated areas fails to select the interpolation direction correctly. While the algorithm can fill large misclosures caused by trees and building rooftops, it is still difficult to reconstruct road segments if it is totally blocked by trees or building rooftops.
Akel N.A., Kremeike K., Filin S., Sester M., and Doytsher Y., (2005). Dense DTM generalization aided by roads extracted from LIDAR data. In: IAPRS 36 (Part 3/W19), 54–59.
Baltsavias, E., Zhang C.S., (2005). Automated updating of road databases from aerial images. International Journal of Applied Earth Observation and Geoinformation. Vol 6, pp 199–213. Bao Y.F., Li G.P., Cao C.X., Li X.W., Zhang H., He Q.S., Bai L.Y. and Chang C.Y., (2008). Classification of Lidar Point Cloud and Generation of DTM From Lidar Height and Intensity Data In Forested Area. In: IAPRSIS. Vol. XXXVII. Part B3b. Beijing 2008.
Bartels, M., and Hong W., (2006). Segmentation of Lidar Data Using Measures of Distribution. Bartels, M., and Hong W., (2010). Threshold-free object and ground point separation in LIDAR data.
Pattern Recognition Letters 31(10): 1089-1099. Bartels, M., Hong W, and Mason D.C., (2006). DTM Generation from LIDAR Data using Skewness Balancing. In: ICPR. pp. 566 – 569. Hong Kong.
Zhao, J.P., and You S.Y., (2012). Road Network Extraction from Airborne Lidar Data Using Scene Context. In: IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Zhou, L. and Vosselman G., (2012). Mapping curbstones in airborne and mobile laser scanning data. International Journal of Applied Earth Observation and Geoinformation 18(2012), 293-304. Zhu, P., Lu Z., Chen X., Honda K., and Eiumnoh A., (2004). Extraction of city roads through shadow path reconstruction using laser scanning. PE&RS 70(12), 1433–1440.
Zhu. Q.H. and Mordohai P., (2009). A Minimum Cover Approach for Extracting the Road Network from Airborne Lidar Data. In: IEEE 12th Conference on Computer Vision Workshops.