Automatically Tracking Road Centerlines from Low-Frequency GPS Trajectory Data

High-quality digital road maps are essential prerequisites of location-based services and smart city applications. The massive and accessible GPS trajectory data generated by mobile GPS devices provide a new means through which to generate maps. However, due to the low sampling rate and multi-level disparity problems, automatically generating road maps is challenging and the generated maps cannot yet meet commercial requirements. In this paper, we present a GPS trajectory data-based road tracking algorithm, including an active contour-based road centerline refinement algorithm as the necessary post-processing. First, the low-frequency trajectory data were transferred into a density estimation map representing the roads through a kernel density estimator, for a seeding algorithm to automatically generate the initial points of the road-tracking algorithm. Then, we present a template-matching-based road-direction extraction algorithm for the road trackers to conduct simple correction, based on local density information. Last, we present an active contour-based road centerline refinement algorithm, considering both the geometric information of roads and density information. The generated road map was quantitatively evaluated using maps offered by the OpenStreetMap. Compared to other methods, our approach could produce a higher quality map with fewer zig-zag roads, and therefore more accurately represents reality.


Introduction
With the development of location-based services, digital road maps are increasingly needed, making road network extraction (RNE) an active research topic. Especially in dense metropolitan cities, people frequently access navigation, route planning, and real-time traffic services, increasing the need for accurate and up-to-date digital road maps. Compared to the time-consuming and labor-intensive conventional map survey, various RNE methods were developed based on remote sensing (RS) images [1,2], GPS trajectory [3,4], and light detection and ranging (LiDAR) point-cloud data [5], since the 1970s. These semi-automatic and automatic algorithms not only lower the cost but also shorten the time of producing and updating road maps. However, due to the inherent drawbacks of the input data, automatically generating high-quality road maps remains challenging [1,3].
Among the existing methods, extracting road networks from RS images is rather mature, and various methods were developed, such as classification-based [6], mathematic morphology [7], road tracking [8][9][10], and active contour [11]. However, these methods are susceptible to the occlusion and shadow in RS images, therefore discontinuous road segments are easy to appear. At present, semi-automatic methods, which require human input, generally perform better than automatic methods. Road extraction from GPS trajectory data started in 1990; since then, it underwent two stages-map augmentation and map generation. Map augmentation requires an existing map as prior information; conversely, map generation constructs a road map with a topological structure from an empty map. Based on the widely deployed GPS sensors, algorithms can use the collective mobility pattern behind trajectory data to extract semantic information, such as turning restrictions and obtaining GPS trajectory data [12], which is much cheaper compared to RS images or LiDAR data. Given the volume of GPS trajectory data, most of these methods are automatic.
To overcome these issues, some methods adopt an intersection-first approach in constructing road maps [14,18]. As intersections are defined as the convergence of three or more road branches, they possess many features that can be detected with a designed algorithm. Detected intersections can serve as valuable information for identifying the road segments connecting them. There has been a recent surge in road intersection detection methods [18,21], which provide a better foundation for this approach. However, few GPS trajectory data-based methods can fully use known intersections, whereas RS-images-based road tracking is mostly driven by points like intersection-centers falling within the region of the road surface.
Therefore, we propose a GPS trajectory data-based road-tracking algorithm, including an active contour-based road centerline refinement algorithm, to extract smoothed and interconnected road centerlines from low-frequency trajectories. By adopting an automatic seeding technique, this map generation method is automatic and thus capable of handling massive trajectory data. The major contributions of our method include.
(1) We enrich the intersection-first approach in extracting road network from GPS trajectory data, making road tracking an available tool to process low-frequency trajectories, which can make better use of known intersections.
(2) We propose a template-matching-based seeding algorithm to automatically provide starting points for the road-tracking algorithm, using the density information of GPS points. To overcome the disparity problem between the heat zone of intersection region and road branches, a circular sampling approach is introduced.
(3) We propose a template-matching-based road direction extraction algorithm to perform correction when tracking road centerlines, which adapts the road-tracking algorithm to GPS trajectory data. To avoid the high GPS noise problem, the real-time correction is simple but robust.
(4) We propose an active contour-based road centerline refinement algorithm to conduct further correction of extracted road centerlines. We include road orientation restrictions in the active contour model to make better use of intersections. The extended post-processing stage makes use of geometric properties of roads, providing a general solution to zig-zag roads.
Our paper is organized as follows. Section 2 reviews the related work on extracting road networks from GPS trajectory data and RS images. Section 3 describes the proposed method for tracking road centerlines. Section 4 presents a set of experimental results and analyses. Finally, Section 5 discusses the conclusions.

Related Work
Since the input data are heterogenous, road network extraction methods based on GPS trajectory data are very different from methods based on RS images. Typically, the output road map is represented by a simple graph G(V, E), where V is a set of nodes representing the terminal points of roads, and E is a set of edges denoting roads.
First, we introduce the GPS trajectory data-based method. As its input, GPS trajectory data consist of positions of a group of mobile targets across a span of time, sometimes including speed and heading direction. Existing GPS trajectory-based methods can be roughly divided into three categories-(1) clustering-based, (2) KDE, and (3) incremental track insertion.
Clustering-based methods first extract the basic node or edge elements of a road map using a clustering algorithm such as K-means. Then, the road map is constructed by connecting these extracted elements. Edelkamp and Schrödl [13] first extracted a set of seed points spread along the road centerlines, by employing the K-means algorithm to cluster the input GPS points. Then, the adjacent extracted points are connected, based on the connectivity of the trajectories. Karagiorgou [14] first extracted road intersections by employing an agglomerate clustering algorithm to the turn-points, which were defined as where the heading direction changed dramatically. Chen [15] first extracted road segments by employing a graph-based clustering algorithm to a sampling set of input points; then, the intersections were detected by training a support vector machine (SVM) classifier based on their proposed traj-SIFT feature. In conclusion, clustering-based methods generally use the trajectory motion information and therefore suffer more loss from low frequencies.
KDE methods begin with transforming GPS points into a density map through a kernel density estimator, with Gaussian kernel. In most cases, this process occurs in conjunction with discretizing the plane space into the grid of pixels; therefore, imageprocessing algorithms can be applied to that density map regarding it as a gray-scale image. Biagioni and Eriksson [16] applied an image thinning algorithm on a binarized density map to produce a set of road skeletons using various binary thresholds. These skeletons are merged into a gray-scale skeleton, where an edge's confidence of being an actual road is represented by its gray scale. Wang [17] proposed a Morse-theory-based road centerlines extraction method, which can capture both heavily-and sparsely-sampled roads with respect to local density. Both methods can produce road maps with different levels of coverage by tuning the parameters and adopting a post-processing approach to better remove the spurious roads. In conclusion, KDE methods are robust to low frequency, but tend to generate zig-zag roads and spurious roads, which requires an effective map refinement algorithm as post-processing.
Incremental track insertion methods construct road maps by incrementally inserting trajectories into an existing map with the help of a map-matching algorithm. Cao and Krumm [18] first grouped similar input trajectories by simulating the physical attraction between them, then incrementally inserted every trajectory by checking whether its node should merge with an existing graph node. Ahmed and Wenk [19] employed Fréchet distance, allowing a trajectory to be partially matched to the map. Different from the first two categories, incremental track insertion methods support online updating, but require high-quality input trajectory data.
Two of the main challenges in extracting road networks from GPS trajectory data are the low frequency and the multi-level disparity of trajectory data; both are ubiquitous among the feasible trajectory data. First, low sampling rate (average sampling interval over 20 s) greatly increases the uncertainty of trajectory data, posing a serious limitation to clustering-based and incremental track insertion methods. Second, the disparity impacts different levels of the extracted road network. At the road level, the staying points of trajectories lead to abnormal local heat zone and distorted road segments. At the network level, the sheer traffic volume difference between different roads makes algorithm tend to mistake roads being travelled by very few trajectories as spurious roads. Finding a global threshold to separate false clusters or roads from correct ones is infeasible.
To address these challenges, researchers are increasingly focusing on the intersectionfirst approach. Following this approach, road network extraction is divided into two stepsroad intersection detection and intersection-based map construction. The difficulty in each step is reduced by adopting a divide-and-conquer strategy. Karagiorgou [14] first proposed a turn-points clustering-based road intersection detection method; Deng [21] introduced a local G * statistic in clustering-based detection; Zhang [18] proposed a KDE-based detection method, considering the intersection's heat zone feature. The goal of intersection detection is to achieve both good coverage and precision, under the low-frequency input data.
However, in map construction, most of the existing GPS trajectory-data-based methods only use the intersection's position; information such as orientations of its connected roads are neglected. Therefore, we combined suitable RS images-based road network extraction methods with this intersection-first approach to more effectively use the known intersections. Here, we focused on briefly introducing the road tracking and active contour methods; both accept a set of seed points as the input, and present roads with a similar parameter model.
Road tracking methods trace the road centerlines by repeating two key steps-prediction and correction. Until a stopping criterion is satisfied, the road tracker traces a single road starting from an initial seed point, which is provided manually or by an automatic seeding algorithm. At the prediction step, an algorithm predicts the next point most likely on the road centerline given the current state of the road tracker. At the correction step, the prediction is corrected by considering its nearby image information, then the state of the road tracker is also updated. Vosselman [8] applied Kalman filter to model the tracking process as a linear dynamic system, and a least squares profile matching was chosen to provide measurement input to the Kalman filter at the correction step. Movaghati [10] combined the extended Kalman filter (EKF) with a particle filter (PF) to find potential new road branches after the road tracker reached its stopping criterion. Road tracking methods are generally sensitive to initial seed points, and face challenge in overcoming intersections, abrupt changes in road direction, and obstructions during the tracing process.
Last, we introduce the active contour method-it defines an energy function on a deformable line, which recursively fits a deformable line to a contour or centerline of a road in the direction of minimal energy, starting from an initialization position. This model was first proposed by Kass [22]; also known as snake, its energy function consists of both internal and image forces. Internal force (elasticity and stiffness) models the geometric features of a road, whereas image force (intensity gradient) represents certain features of interest in the image. Fua and Leclerc [23] proposed the ribbon snakes model to include a width parameter in the image force. By tuning the ratio between internal force and image force, active contour methods can avoid the overfitting problem when extracting roads from an image, which often yields zig-zag roads. As a result, it can also be used as a geometric refinement algorithm for road networks, given background image information.
In conclusion, Figure 1 outlines the rough categories of GPS trajectory-data-based and RS-images-based road network extraction methods. Road tracking methods can also benefit from adopting the intersection-first approach-if intersections and their connected road branches can be detected with high precision and high coverage, the difficulty of road tracking would be significantly reduced. However, due to GPS noise and multi-level disparity problems, the RS images-based methods cannot be directly applied to the density estimation map of GPS points.
Our method falls into the intersection-first approach category. Based on previous road intersection detection works, we present a combination of KDE and road tracking methods to overcome low-frequency and multi-level disparity of GPS trajectory data in road network extraction. First, we propose a density-based seeding algorithm and a road direction extraction algorithm, both based on template matching, to adapt road tracking to the density estimation map. We adopt a simple real-time correction in this tracking process to avoid tackling drastic GPS noise with limited local information. Then, we apply an active contour-based refinement algorithm to the extracted road network to recursively redo the correction, considering both the geometric information of an entire road and the density information. We emphasize that the impact of disparity is many-fold; by adopting a divide-and-conquer strategy, this problem can be easier solved at the local level. Our method falls into the intersection-first approach category. Based on previous road intersection detection works, we present a combination of KDE and road tracking methods to overcome low-frequency and multi-level disparity of GPS trajectory data in road network extraction. First, we propose a density-based seeding algorithm and a road direction extraction algorithm, both based on template matching, to adapt road tracking to the density estimation map. We adopt a simple real-time correction in this tracking process to avoid tackling drastic GPS noise with limited local information. Then, we apply an active contour-based refinement algorithm to the extracted road network to recursively redo the correction, considering both the geometric information of an entire road and the density information. We emphasize that the impact of disparity is many-fold; by adopting a divide-and-conquer strategy, this problem can be easier solved at the local level.

Tracking Road Centerlines from Low-Frequency Trajectory Data
Our road network extraction method starts with a set of known intersection positions, which can be automatically provided by previous road intersection detection works [18,21] or can be input manually. The proposed method consists of three key steps for extracting road centerlines from low-frequency GPS trajectories-(1) generating seed points that exclusively define one end of a road; (2) tracking road centerlines with simple correction; and (3) refining road centerlines based on the active contour model. The workflow of our road centerline tracking method is shown in Figure 2.

Tracking Road Centerlines from Low-Frequency Trajectory Data
Our road network extraction method starts with a set of known intersection positions, which can be automatically provided by previous road intersection detection works [18,21] or can be input manually. The proposed method consists of three key steps for extracting road centerlines from low-frequency GPS trajectories-(1) generating seed points that exclusively define one end of a road; (2) tracking road centerlines with simple correction; and (3) refining road centerlines based on the active contour model. The workflow of our road centerline tracking method is shown in Figure 2.

Automatic Seeding Based on Kernel Density Estimation Map
Seed points are the prerequisite of road-tracking algorithms. In our method, a seed point is defined as ( , , ), where ( , ) ∈ is the start point of road and ∈ [0,2 ] is the orientation of road at that point (anti-clockwise, equals 0 when heading toward Earth's true east). Every seed point is derived from an existing intersection, and exclusively defines one end of a road.
Given a known intersection's center position ( , ), the possibility of having a road

Automatic Seeding Based on Kernel Density Estimation Map
Seed points are the prerequisite of road-tracking algorithms. In our method, a seed point is defined as (x, y, θ), where (x, y) ∈ R 2 is the start point of road and θ ∈ [0, 2π] is the orientation of road at that point (anti-clockwise, equals 0 when heading toward Earth's true east). Every seed point is derived from an existing intersection, and exclusively defines one end of a road.
Given a known intersection's center positio I i (x i , y i ), the possibility of having a road branch with orientation θ is quantitively evaluated by matching the corresponding slice image with a set of ideal road template images. Least squares correlation matching is applied, along with a circular sampling approach, tackling the disparity problem of GPS trajectories. Based on the minimum-matching-distance criterion, multiple seed points x i , y i , θ i,j can be drawn, and then be examined by the seeding algorithm. Here we use a point-based intersection model, which can cover a majority of intersections with cross-shape, T-shape, and Y-shape.

Density Estimation
Based on the well-known KDE algorithm [24], the sparsely sampled GPS points can be transformed into a consecutive density function ρ(x, y) on a rectangular domain R ⊂ R 2 , which is the area of interest for road network extraction.
Since using a consecutive density function is calculation-intensive, a common practice is discretizing R into a grid G R with a cell-length c, and then only computing ρ(x, y) at the grid center points. Therefore, the density estimation result can be viewed as a gray-scale image D(i, j), where D(i, j) stands for the density at grid i, j.
Let P = → p 1 , . . . , → p n be the collection of all GPS points. First, a 2D histogram H(i, j) is produced by counting the number of GPS points falling within each grid. Then, the density image D(i, j) is computed as: where K σ (u, v) is the popular Gaussian kernel and has a bandwidth parameter σ. Following a previous work [16], we used σ = 8.5 m.
It is faster to compute Equation (1) in two rounds of 1D convolution. Last, we discard the marginal information of K σ (u, v), keeping ∑ |u|,|v|≤M K σ (u, v) ≥ 0.995 as the central information. These improvements support using a smaller grid cell-length c and conducting faster template-matching-based algorithms.

Template Matching
The previous section demonstrated a method to generate global density image from input GPS points. However, due to the GPS noise and the disparity problem [16], automatic seeding requires focusing on the local density image (slice image) to apply pre-processing, such as median filters and image binarization.
Given a position (x, y), Figure 3 shows how a slice image with orientation θ is generated. Let S(x, y, θ) denote the slice image after rotation, for the convenience of matching with a set of fixed road template images. Its pixel value s(i, j) stands for the density at the recalculated grid center → x ij , which is computed as: Then, density value s(i, j) is re-estimated using Equation (1) and the nearby GPS points: where K σ (d) is the Gaussian kernel function. N is the image size satisfying N·c = 80 m, which means the slice image covers an area of 80 × 80 m.
Then, density value ( , ) is re-estimated using Equation (1) and the nearby GPS points: where ( ) is the Gaussian kernel function. is the image size satisfying • = 80 m, which means the slice image covers an area of 80 × 80 m.
(a) (b) (c) Figure 3. Generated slice image with arbitrary orientation-(a) the square local area; (b) re-estimating density value on every pixel using the nearby points; and (c) rotating the skewed slice image.
Next, we prepare a set of template images based on an ideal road model. A previous work [25] showed that a Gaussian mixture model (GMM) can model the spread of GPS points within a road, reflecting the number of lanes and their corresponding traffic volume in reality. Here, we use single Gaussian model as a simplification, as shown in Figure  4. Therefore, in a road template image, density of one pixel is solely determined by the distance between the pixel center and road centerline.

Let
( , ) denote a road template image with size × pixels, its pixel value ( , ) can be computed as: Next, we prepare a set of template images based on an ideal road model. A previous work [25] showed that a Gaussian mixture model (GMM) can model the spread of GPS points within a road, reflecting the number of lanes and their corresponding traffic volume in reality. Here, we use single Gaussian model as a simplification, as shown in Figure 4. Therefore, in a road template image, density of one pixel is solely determined by the distance between the pixel center and road centerline.
Then, density value ( , ) is re-estimated using Equation (1) and the nearby GPS points: where ( ) is the Gaussian kernel function. is the image size satisfying • = 80 m, which means the slice image covers an area of 80 × 80 m.
(a) (b) (c) Figure 3. Generated slice image with arbitrary orientation-(a) the square local area; (b) re-estimating density value on every pixel using the nearby points; and (c) rotating the skewed slice image.
Next, we prepare a set of template images based on an ideal road model. A previous work [25] showed that a Gaussian mixture model (GMM) can model the spread of GPS points within a road, reflecting the number of lanes and their corresponding traffic volume in reality. Here, we use single Gaussian model as a simplification, as shown in Figure  4. Therefore, in a road template image, density of one pixel is solely determined by the distance between the pixel center and road centerline.

Let
( , ) denote a road template image with size × pixels, its pixel value ( , ) can be computed as: Let Temp(α, B) denote a road template image with size N × N pixels, its pixel value t(i, j) can be computed as: where α is the orientation of road, and B is the bias parameter. The bandwidth τ was determined by 3τ = N/2 to make the template image cover the entire density descending course. By shifting Temp(α, 0) B pixels' width perpendicular to its road orientation, we obtain Temp(α, B). The sizes of the slice image and template image were set to the same value throughout our method. If the slice image was close enough to a template image, one road was detected in that slice image. When the slice image deviated from the road centerline, the road could still be detected from a broader template set {Temp(α, B)|B = 0, ±1, . . . , ±K}. Since a traffic volume difference might exist between different roads, we could not directly compare pixel value s(i, j) with t(i, j). Therefore, we applied median filter and binarization to both images as pre-processing, before the matching. Figure 5 shows that the median filter could effectively remove the GPS noise, then the traffic volume was normalized by binarization.
If the slice image was close enough to a template image, one road was detected in that slice image. When the slice image deviated from the road centerline, the road could still be detected from a broader template set { ( , )| = 0, ±1, … , ± }. Since a traffic volume difference might exist between different roads, we could not directly compare pixel value ( , ) with ( , ). Therefore, we applied median filter and binarization to both images as pre-processing, before the matching. Figure 5 shows that the median filter could effectively remove the GPS noise, then the traffic volume was normalized by binarization. Last, the matching should be based on a distance measure between two images. We chose the squared distance for calculation convenience. Given an image , let denote the image after the above pre-processing. The matching distance between two images was defined as: The smaller the matching distance, the more similar were the two images.

Seeding Based on Template Matching and Circular Sampling
After these preparations, we can now quantitively evaluate the possibility of a road connected to an intersection from arbitrary directions. The search was performed by circular sampling to avoid different road branches, overlapping within the intersection region.
As shown in Figure 6, we categorized the slice image into three classes-(1) road segment, (2) near intersection, and (3) highly noisy area. Experiments showed that the template-matching algorithm worked for road segments and highly noisy areas, but suffered a significant drop in accuracy for the areas near intersections. Detailed experimental results are provided in Section 4. Therefore, we proposed a circular sampling-based seeding algorithm, as shown in Figure 7, to avoid the drawback of the near-intersection areas. Given an intersection's center position ( , ), we used the slice image ( + • , + • , ) to detect whether Last, the matching should be based on a distance measure between two images. We chose the squared distance for calculation convenience. Given an image S, let S denote the image after the above pre-processing. The matching distance between two images was defined as: The smaller the matching distance, the more similar were the two images.

Seeding Based on Template Matching and Circular Sampling
After these preparations, we can now quantitively evaluate the possibility of a road connected to an intersection from arbitrary directions. The search was performed by circular sampling to avoid different road branches, overlapping within the intersection region.
As shown in Figure 6, we categorized the slice image into three classes-(1) road segment, (2) near intersection, and (3) highly noisy area. Experiments showed that the template-matching algorithm worked for road segments and highly noisy areas, but suffered a significant drop in accuracy for the areas near intersections. Detailed experimental results are provided in Section 4. that slice image. When the slice image deviated from the road centerline, the road could still be detected from a broader template set { ( , )| = 0, ±1, … , ± }. Since a traffic volume difference might exist between different roads, we could not directly compare pixel value ( , ) with ( , ). Therefore, we applied median filter and binarization to both images as pre-processing, before the matching. Figure 5 shows that the median filter could effectively remove the GPS noise, then the traffic volume was normalized by binarization. Last, the matching should be based on a distance measure between two images. We chose the squared distance for calculation convenience. Given an image , let denote the image after the above pre-processing. The matching distance between two images was defined as: The smaller the matching distance, the more similar were the two images.

Seeding Based on Template Matching and Circular Sampling
After these preparations, we can now quantitively evaluate the possibility of a road connected to an intersection from arbitrary directions. The search was performed by circular sampling to avoid different road branches, overlapping within the intersection region.
As shown in Figure 6, we categorized the slice image into three classes-(1) road segment, (2) near intersection, and (3) highly noisy area. Experiments showed that the template-matching algorithm worked for road segments and highly noisy areas, but suffered a significant drop in accuracy for the areas near intersections. Detailed experimental results are provided in Section 4. Therefore, we proposed a circular sampling-based seeding algorithm, as shown in Figure 7, to avoid the drawback of the near-intersection areas. Given an intersection's center position ( , ), we used the slice image ( + • , + • , ) to detect whether Therefore, we proposed a circular sampling-based seeding algorithm, as shown in Figure 7, to avoid the drawback of the near-intersection areas. Given an intersection's center position (x, y), we used the slice image S(x + r·cosθ, y + r·sinθ, θ) to detect whether there exists a road branch with orientation θ. By matching S with a set of template images, a matching distance function f x,y (θ) with the orientation variable was defined as: where µ represents how far a deviation we accept from the road centerline; here, we set µ = 0.125. Radius r was set to 60 m through a series of experiments. Intuitively, smaller f x,y (θ) represents a higher chance of having a road branch with orientation θ. Figure 8 shows that each road branch corresponds to a local minimum of f x,y (θ). As a result, by finding all local minimum points {θ k } of f x,y (θ), the corresponding seed points {(x, y, θ k )} were generated. Since the number of intersections was limited, we adopted a simple enumerating solution in finding these local minimum points. We only computed f x,y (θ) on a discretized set {θ i = 2iπ/n| i = 0, . . . , n}, an angle θ i that satisfied f x,y (θ i ) = min |θ j −θ i |≤π/12 f x,y θ j would be selected by our seeding algorithm. a matching distance function , ( ) with the orientation variable was defined as: where represents how far a deviation we accept from the road centerline; here, we set = 0.125 . Radius was set to 60 m through a series of experiments. Intuitively, smaller , ( ) represents a higher chance of having a road branch with orientation . Figure 7. Detecting road branches based on template matching and circular sampling. Figure 8 shows that each road branch corresponds to a local minimum of , ( ). As a result, by finding all local minimum points { } of , ( ) , the corresponding seed points {( , , )} were generated. Since the number of intersections was limited, we adopted a simple enumerating solution in finding these local minimum points. We only computed , ( ) on a discretized set { = 2 ⁄ | = 0, … , }, an angle that satisfied would be selected by our seeding algorithm. However, as shown in Figure 9, the fluctuation of , ( ) caused by GPS noise led to many false seed points. By applying Gaussian smoothing on , ( ), this problem could only be mitigated; therefore, we proposed density-based connectivity to verify the generated seed points. Given two planar points ⃗ , ⃗ , let ⃗ , ⃗ denote the connectivity between these points.   Figure 8 shows that each road branch corresponds to a local minimum of , ( ). As a result, by finding all local minimum points { } of , ( ) , the corresponding seed points {( , , )} were generated. Since the number of intersections was limited, we adopted a simple enumerating solution in finding these local minimum points. We only computed , ( ) on a discretized set { = 2 ⁄ | = 0, … , }, an angle that satisfied would be selected by our seeding algorithm. However, as shown in Figure 9, the fluctuation of , ( ) caused by GPS noise led to many false seed points. By applying Gaussian smoothing on , ( ), this problem could only be mitigated; therefore, we proposed density-based connectivity to verify the generated seed points. Given two planar points ⃗ , ⃗ , let ⃗ , ⃗ denote the connectivity between these points.
, ( , ) = ⃗ + ⃗ − ⃗ . However, as shown in Figure 9, the fluctuation of f x,y (θ) caused by GPS noise led to many false seed points. By applying Gaussian smoothing on f x,y (θ), this problem could only be mitigated; therefore, we proposed density-based connectivity to verify the generated seed points. Given two planar points ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 10 of 25 The proposed connectivity ranged from 0 to 1. Assuming ⃗ is the intersection's center position, when ⃗ refers to a false road branch, the density value drops quickly along ⃗ , which results in low connectivity. Figure 10 shows that true road branches tend to have larger integral area under curve ( , ), therefore, a connectivity threshold can effectively separate the false road branches from the true ones. Figure 10 also shows that, in most cases, the start point ⃗ has a much higher density value for traffic from different road branches overlapping at the intersection, which results in a lower ⃗ , ⃗ in general. Therefore, when verifying a seed point ( , , ) by computing density-based connectivity, the segment falling within the intersection region should be truncated, as follows: The proposed connectivity ranged from 0 to 1. Assuming → P is the intersection's center position, when → PQ refers to a false road branch, the density value drops quickly along → PQ, which results in low connectivity. Figure 10 shows that true road branches tend to have larger integral area under curve ρ(x t , y t ), therefore, a connectivity threshold can effectively separate the false road branches from the true ones.
where ⃗ = ( , ) + • ( , ) and is a predefined connectivity threshold. We set = 0.1 through a series of experiments. The other two parameters = 30 m and = 60 m were set through visualization, as shown in Figure 11. Seed points that satisfied Equation (8) were considered true, and were fed to the road-tracking algorithm in the next step.

Tracking Road Centerlines Based on Kernel Density Estimation Map
After the seeding algorithm generated a set of seed points from the density map, the tracking algorithm set up multiple road trackers to incrementally generate estimation curves of road centerlines. Each tracker corresponded to a single seed point ( , , ). The estimation curves were further refined in the next post-processing course.
Here, we used a simple graph ( , ) to model the output road network, with known intersections as its initial vertex set . For every road tracker, the prediction and correction two-step process was repeated-at the prediction step, the algorithm calculated the position of a new vertex based on the current state of the road tracker (using a predefined road  Figure 10 also shows that, in most cases, the start point → P has a much higher density value for traffic from different road branches overlapping at the intersection, which results in a lower c → P, → Q in general. Therefore, when verifying a seed point (x, y, θ) by computing density-based connectivity, the segment falling within the intersection region should be truncated, as follows: where → P d = (x, y) + d·(cosθ, sinθ) and ϕ is a predefined connectivity threshold. We set ϕ = 0.1 through a series of experiments. The other two parameters l = 30 m and L = 60 m were set through visualization, as shown in Figure 11. Seed points that satisfied Equation (8) were considered true, and were fed to the road-tracking algorithm in the next step. The proposed connectivity ranged from 0 to 1. Assuming ⃗ is the intersection's center position, when ⃗ refers to a false road branch, the density value drops quickly along ⃗ , which results in low connectivity. Figure 10 shows that true road branches tend to have larger integral area under curve ( , ), therefore, a connectivity threshold can effectively separate the false road branches from the true ones. Figure 10 also shows that, in most cases, the start point ⃗ has a much higher density value for traffic from different road branches overlapping at the intersection, which results in a lower ⃗ , ⃗ in general. Therefore, when verifying a seed point ( , , ) by computing density-based connectivity, the segment falling within the intersection region should be truncated, as follows: where ⃗ = ( , ) + • ( , ) and is a predefined connectivity threshold. We set = 0.1 through a series of experiments. The other two parameters = 30 m and = 60 m were set through visualization, as shown in Figure 11. Seed points that satisfied Equation (8) were considered true, and were fed to the road-tracking algorithm in the next step.

Tracking Road Centerlines Based on Kernel Density Estimation Map
After the seeding algorithm generated a set of seed points from the density map, the tracking algorithm set up multiple road trackers to incrementally generate estimation curves of road centerlines. Each tracker corresponded to a single seed point ( , , ). The estimation curves were further refined in the next post-processing course.
Here, we used a simple graph ( , ) to model the output road network, with known intersections as its initial vertex set . For every road tracker, the prediction and correction two-step process was repeated-at the prediction step, the algorithm calculated the position of a new vertex based on the current state of the road tracker (using a predefined road

Tracking Road Centerlines Based on Kernel Density Estimation Map
After the seeding algorithm generated a set of seed points from the density map, the tracking algorithm set up multiple road trackers to incrementally generate estimation curves of road centerlines. Each tracker corresponded to a single seed point (x, y, θ). The estimation curves were further refined in the next post-processing course.
Here, we used a simple graph G(V, E) to model the output road network, with known intersections as its initial vertex set V. For every road tracker, the prediction and correction two-step process was repeated-at the prediction step, the algorithm calculated the position of a new vertex based on the current state of the road tracker (using a predefined road model), and the successive positions formed the edges of the graph G. At the correction step, the algorithm measured the local density map around the newest position to refine the prediction and update the current state of the road tracker. As such, new edges were incrementally added to the shared graph G, until a stopping criterion was satisfied. Figure 12 shows the flow chart of our road-tracking algorithm.
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 11 of 25 model), and the successive positions formed the edges of the graph . At the correction step, the algorithm measured the local density map around the newest position to refine the prediction and update the current state of the road tracker. As such, new edges were incrementally added to the shared graph , until a stopping criterion was satisfied. Figure  12 shows the flow chart of our road-tracking algorithm.

Road Model
We first introduced a road model that was used by previous methods [10], which models one road as a series of states formed a polyline representing the road centerline. Figure 13 shows that this definition of road state contained the necessary information to predict the most likely next point on the road centerline, along with other crucial road attributes. Previous methods also modeled the tracking process as a stochastic process [8,10]. As a result, the prediction step could be written in an equation form:

Road Model
We first introduced a road model that was used by previous methods [10], which models one road as a series of states {X 0 , X 1 , . . . , Here, attribute θ k represented the orientation of the road at position (x k , y k ), . θ k represented the curvature of road at (x k , y k ), and {(x k , y k )} 0≤k≤n formed a polyline representing the road centerline. Figure 13 shows that this definition of road state contained the necessary information to predict the most likely next point on the road centerline, along with other crucial road attributes. model), and the successive positions formed the edges of the graph . At the correction step, the algorithm measured the local density map around the newest position to refine the prediction and update the current state of the road tracker. As such, new edges were incrementally added to the shared graph , until a stopping criterion was satisfied. Figure  12 shows the flow chart of our road-tracking algorithm.

Road Model
We first introduced a road model that was used by previous methods [10], which models one road as a series of states formed a polyline representing the road centerline. Figure 13 shows that this definition of road state contained the necessary information to predict the most likely next point on the road centerline, along with other crucial road attributes. Previous methods also modeled the tracking process as a stochastic process [8,10]. As a result, the prediction step could be written in an equation form: Previous methods also modeled the tracking process as a stochastic process [8,10]. As a result, the prediction step could be written in an equation form: where s is the tracking step length parameter and W k represents the random variation in the road state between successive positions. The next step was measuring the local density map to support the correction step, which is referred to as another equation called the measurement equation: where V k represents the measurement noise. Assuming that both Equations (9) and (10) were linear and both W k and V k were Gaussian noise, the classical correction based on the Kalman filter could be derived [8]. However, since precisely modeling W k and V k was infeasible for GPS trajectory data under complex urban environments, we adopted a simple correction in the tracking process by replacing W k and V k with a zero matrix. The underlying assumption was that roads have a fixed curvature except for a few turn-points.

Road Orientation Extraction Based on Template Matching
To support the correction step, we measured two attributes (θ k and B k ) from the local density map of the predicted state X k|k−1 = f (X k−1 ), based on the template matching previously introduced in Section 3.1.2. Attributeθ k represented the road orientation at the predicted point x k|k−1 , y k|k−1 , and B k represented the grid deviation from the road centerline at a predicted point.
Since the tracking process required excessive measurements, we introduced another matching mode requiring far less computation. Recall that in Section 3.1.3, generating a slice image S(x + r·cosθ, y + r·sinθ, θ) by redoing the density estimation placed a heavy computation burden on the seeding algorithm. Therefore, in this section, we only computed a slice image with fixed θ = 0, and let T(∆θ, µ) denote the enlarged template image set, based on discretized orientations: Using the same minimum-matching-distance criterion, the measurement result was obtained as:θ where Temp(θ k , B k ) belonged to the pre-generated set T(∆θ, µ). Here, we used ∆θ = 0.25 • and µ = 0.125, which could satisfy the simple correction needed. Figure 14 demonstrates the difference between the two matching modes. By allowing the slice image to have an arbitrary orientation, the algorithm could use more local density information, but the computation cost was not equivalent to the gain. where is the tracking step length parameter and represents the random variation in the road state between successive positions. The next step was measuring the local density map to support the correction step, which is referred to as another equation called the measurement equation: where represents the measurement noise. Assuming that both Equations (9) and (10) were linear and both and were Gaussian noise, the classical correction based on the Kalman filter could be derived [8].
However, since precisely modeling and was infeasible for GPS trajectory data under complex urban environments, we adopted a simple correction in the tracking process by replacing and with a zero matrix. The underlying assumption was that roads have a fixed curvature except for a few turn-points.

Road Orientation Extraction Based on Template Matching
To support the correction step, we measured two attributes ( and ) from the local density map of the predicted state | = ( ), based on the template matching previously introduced in Section 3.1.2. Attribute represented the road orientation at the predicted point ( | , | ), and represented the grid deviation from the road centerline at a predicted point.
Since the tracking process required excessive measurements, we introduced another matching mode requiring far less computation. Recall that in Section 3.1.3, generating a slice image ( + • , + • , ) by redoing the density estimation placed a heavy computation burden on the seeding algorithm. Therefore, in this section, we only computed a slice image with fixed = 0, and let (∆ , ) denote the enlarged template image set, based on discretized orientations: Using the same minimum-matching-distance criterion, the measurement result was obtained as: where ( , ) belonged to the pre-generated set (∆ , ). Here, we used ∆ = 0.25° and = 0.125, which could satisfy the simple correction needed. Figure 14 demonstrates the difference between the two matching modes. By allowing the slice image to have an arbitrary orientation, the algorithm could use more local density information, but the computation cost was not equivalent to the gain.
(a) (b) Figure 14. Difference between two template-matching modes-(a) slice images with variable orientation allow for using more image information; and (b) fixed slice image is less computationally intensive.

Road Tracking Algorithm
In this section, we combined the road model and the measuring algorithm to present our road-tracking algorithm. Starting with the density map ( , ) and a set of seed points, this algorithm incrementally built a road network ( , ) with a topological structure. The tracking algorithm consisted of the following five steps: Figure 14. Difference between two template-matching modes-(a) slice images with variable orientation allow for using more image information; and (b) fixed slice image is less computationally intensive.

Road Tracking Algorithm
In this section, we combined the road model and the measuring algorithm to present our road-tracking algorithm. Starting with the density map D(i, j) and a set of seed points, this algorithm incrementally built a road network G(V, E) with a topological structure. The tracking algorithm consisted of the following five steps: (1) Initialization: For every seed point (x, y, θ), create a road tracker with an initial state X 0 := [x, y, θ, 0] T , and let distinguished seed positions {(x, y)} be the initial vertex set of road network G.
(2) Prediction: For every road tracker, based on the current state X k−1 and Equation (9), compute the predicted state X k|k−1 = f (X k−1 ).
(3) Measurement: Based on the predicted state X k|k−1 and Equation (12), obtain measurement resultθ k and B k from density map D.
(4) Correction: Based on the measurement result and the following equation, compute the corrected state X k|k : where c is the grid cell length and s is the tracking step length.
If Equation (14) is satisfied, then → v k|k is absorbed by the closer vertex → P or → Q of the road network, and a new edge from → v k−1 to the absorb-vertex is added to the road network. Then, the corresponding road tracker stops. Otherwise, examine the density-connectivity criterion using Equation (7).
where ϕ is the same connectivity threshold in Equation (8). If Equation (15) is satisfied, it is likely that the newly extracted road segment → v k|k crossed a road edge or ran into a dead end (zero density area). The corresponding road tracker stopped and the new edge was nullified in this situation. Otherwise, a new vertex and edge were added to the road network, and the road tracker continued tracking with the newest state X k = X k|k , by repeating steps (2) to (5).
The tracking step length played a key role in our road-tracking algorithm. Intuitively, when the step length was too large, the road tracker could not handle the dramatic change in road orientation well. Conversely, when the step length was too small, the road tracker would be trapped or distorted by the local heat zone of the density maps, such as the intersection regions. Here, we chose s = 50 m through a series of experiments.
As shown in Figure 15, after the input of ground-truth intersections, our road-tracking algorithm could generate continuous road centerlines interconnected to form a network, and achieved good coverage of existing roads (travelled by at least three trajectories). This provided a solid foundation for post-processing, to further enhance the quality of the road network.
In conclusion, the proposed template-matching-based road-orientation-extraction algorithm adapted the road-tracking algorithm for the GPS trajectory data. Although the correction step was simple, it considered 80 × 80 m area of local density information and was robust enough. A road-intersection-detection algorithm with high recall could help the road-tracking algorithm to avoid tackling the challenging intersection regions in advance. In conclusion, the proposed template-matching-based road-orientation-extraction a gorithm adapted the road-tracking algorithm for the GPS trajectory data. Although th correction step was simple, it considered 80 × 80 m area of local density information an was robust enough. A road-intersection-detection algorithm with high recall could hel the road-tracking algorithm to avoid tackling the challenging intersection regions in ad vance.

Road Centerline Refinement Based on Active Contour
Due to the high GPS noise, most existing KDE-based road-network-extraction meth ods produced a number of zig-zag roads in the extraction result, as shown in Figure 1 The proposed road tracking algorithm faced the same problem-due to the discretize output of the measurement algorithm, the extracted road centerlines often have fluctua ing orientations on consecutive points, resulting in zig-zag roads, especially when th tracking step length is small.
Developing a density-aware road-centerline-refinement algorithm provided a gen eral solution to this problem. We chose an active contour model for it defined an energ function on a deformable line, including a pair of structural internal force and externa image force; therefore, it could balance between generating straight roads and over-fittin zig-zag roads. The refinement algorithm took a road centerline from a road-tracking algo rithm as its input; based on a predefined energy function, the geometric refinement of th road centerline was done by minimizing its energy, yielding a smoothed road.

Active Contour Model
In this section, we introduce the definition of the original active contour model b Kass [22], along with the minimization procedure, through variational calculu Let ⃗( ) = ( ), ( ) denote a deformable line in plane space, its energy function is de fined as: represents the internal energy, which controls the shape of curve ⃗( ) with i first-order term and second-order term; represents the image energy when ⃗( ) placed on a background image. Parameter ( ) is also known as the elasticity paramete

Road Centerline Refinement Based on Active Contour
Due to the high GPS noise, most existing KDE-based road-network-extraction methods produced a number of zig-zag roads in the extraction result, as shown in Figure 16. The proposed road tracking algorithm faced the same problem-due to the discretized output of the measurement algorithm, the extracted road centerlines often have fluctuating orientations on consecutive points, resulting in zig-zag roads, especially when the tracking step length is small. To determine , since orientation restrictions were independent from image information (we could assume that the image force equaled to zero), the following equation held true for the optimization result: Focusing on one side of the road, one could make a reasonable assumption that [ ][ ⃗ ⃗ ⃗ ⃗ ⃗ ] = 0 also held true, where ⃗ = ⃗ − ℎ ⃗ 0 was defined by prolonging the road in the opposite orientation. Based on these assumptions, ⃗ and ⃗ could now be solved using the first two rows of coefficient matrix : where ℎ ⃗ = ( • , • ) was computed by using the road orientation and tracking step length. By similarly computing ⃗ and ⃗ , yielded the complete restriction ma- Previous works demonstrated that road centerlines could be viewed as mountain ridges of the density field [17]; therefore, we define the image force as a negative density Developing a density-aware road-centerline-refinement algorithm provided a general solution to this problem. We chose an active contour model for it defined an energy function on a deformable line, including a pair of structural internal force and external image force; therefore, it could balance between generating straight roads and over-fitting zig-zag roads. The refinement algorithm took a road centerline from a road-tracking algorithm as its input; based on a predefined energy function, the geometric refinement of the road centerline was done by minimizing its energy, yielding a smoothed road.

Active Contour Model
In this section, we introduce the definition of the original active contour model by Kass [22], along with the minimization procedure, through variational calculus. Let → v (s) = (x(s), y(s)) denote a deformable line in plane space, its energy function is defined as: where E int represents the internal energy, which controls the shape of curve → v (s) with its first-order term and second-order term; E img represents the image energy when → v (s) is placed on a background image. Parameter α(s) is also known as the elasticity parameter, where setting a larger α(s) results in an optimal curve closer to the shortest path; and parameter β(s) is also known as the stiffness parameter, where setting a smaller β(s) allows the optimal curve to have a more drastic change in direction.
Based on Equation (16), the optimal curve was defined and obtained by minimizing the total energy. Kass proposed a numerical method using the Euler-Lagrange differential → v ss ds, then the following equation was satisfied: Turning Equation (16) into discrete form yields an easier solution method, especially when α(s) and β(s) were not constant. Following Kass' approach, two independent Euler-Lagrange equations could be drawn from Equation (17), and could be written in matrix form as: Solving (X, Y) from Equation (18) yields the optimal curve, which could be achieved by recursively computing the following equations, starting from the initial curve (X 0 , Y 0 ): where γ is a step parameter and I is an n × n identity matrix. The inverse matrix of (A + γI) could be computed by LU decomposition. When the total energy change was less than a threshold, the curve (X t , Y t ) converged to a local minima.

Geometric Refinement of the Extracted Road Centerlines
Based on the active contour model, the geometric refinement of a road centerline could be defined as minimizing its energy, yielding a smoothed road. The refinement algorithm took its input → v n from the road-tracking algorithm. To solve the optimal road centerlines, first, we computed the coefficient matrix A. In our case, curve → v (s) was unclosed and had two fixed ends, resulting in the boundary con- → v (n) were fixed, we could only optimize the remaining points; therefore, A was an (n − 1) × (n − 1) pentadiagonal banded matrix: where a = 2α + 6β, b = −(α + 4β), c = β by Kass' approach [22]. Following the common assumption made by road-tracking methods that road centerlines generally have low curvatures, we could assume that α(s) = α and β(s) = β were constant, except for a few turn-points. If → v i is a turn-point, then β i is set to 0, and the ith row of A equals [· · · 0 −α 2α −α 0 · · ·]. The adaptive selection of turn-points is discussed later.
Next, we introduced orientation restrictions was satisfied, as we want orientation restrictions to only directly affect the first points of a road from two ends.
To determine B, since orientation restrictions were independent from image information (we could assume that the image force equaled to zero), the following equation held true for the optimization result: Focusing on one side of the road, one could make a reasonable assumption that prolonging the road in the opposite orientation. Based on these assumptions, → b 1 and → b 2 could now be solved using the first two rows of coefficient matrix A: where → h 0 = (s·cosθ 0 , s·sinθ 0 ) was computed by using the road orientation and tracking step length. By similarly computing Previous works demonstrated that road centerlines could be viewed as mountain ridges of the density field [17]; therefore, we define the image force as a negative density value P → v (s) = −ρ → v (s) . However, due to the disparity problem, the image force varied largely from curve to curve, making it infeasible to apply a set of global parameters. Therefore, the image force must be normalized at the local level. Combining the orientation restrictions, Equation (18) was rewritten as: where λ is the normalization factor, which was computed using the density value along The new equation could still be solved by Kass' approach (see Equation (19)). Last, we identified four tuning strategies through a series of experiments: (1) Set a rather small α; therefore, the image force plays a priority role in refining road centerlines; (2) Set a rather large β/α; here, we use β/α = 10; (3) Set a rather large step parameter γ/α; here, we use γ/α = 5; (4) Allow the curve to have multiple adaptively selected turn-points, where β = 0 is satisfied.
We briefly introduced these strategies. Since the input curve from the road-tracking algorithm was close enough to the real road centerline, a large step parameter could produce a robust iteration course. The impact of disparity was multi-level; local heat zones existed in the road due to the stay-points in the trajectories. Figure 17 shows that these heat zones might distort an active contour, as the image force dragged its points into a density peak. Setting a larger β could mitigate this problem; however, it produced a new problem-if α(s) and β(s) were constant, the active contour could not accurately fit a road with drastic changes in orientation (Figure 18a). A previous work [11] showed that introducing turn-points could solve this problem (Figure 18b). At a turn-point, was set to 0, allowing the optimal curve to drastically change in orientation. Turn-points could be adaptively selected by checking whether the orientation change | − | in the initial curve was larger than a given threshold (e.g., 6 ⁄ ). Experiments showed that the adaptive turn-points could considerably improve the geometric refinement of road centerlines, which significantly reduced the zigzag roads in the extracted road network.

Datasets and Evaluation of the Extracted Road Network
Real trajectory data of taxis in the Shenzhen city [12] were used to verify the effectiveness of our road network extraction method using low-frequency GPS trajectory data. The Shenzhen dataset contained over 75 million GPS points from 14,716 taxis, with an average sample rate of 26.1 s [26], collected in one day, which represented feasible but A previous work [11] showed that introducing turn-points could solve this problem (Figure 18b). At a turn-point, was set to 0, allowing the optimal curve to drastically change in orientation. Turn-points could be adaptively selected by checking whether the orientation change | − | in the initial curve was larger than a given threshold (e.g., 6 ⁄ ). Experiments showed that the adaptive turn-points could considerably im prove the geometric refinement of road centerlines, which significantly reduced the zig zag roads in the extracted road network.

Datasets and Evaluation of the Extracted Road Network
Real trajectory data of taxis in the Shenzhen city [12] were used to verify the effec tiveness of our road network extraction method using low-frequency GPS trajectory data The Shenzhen dataset contained over 75 million GPS points from 14,716 taxis, with an A previous work [11] showed that introducing turn-points could solve this problem (Figure 18b). At a turn-point, β i was set to 0, allowing the optimal curve to drastically change in orientation. Turn-points could be adaptively selected by checking whether the orientation change |θ i − θ i−1 | in the initial curve was larger than a given threshold (e.g., π/6). Experiments showed that the adaptive turn-points could considerably improve the geometric refinement of road centerlines, which significantly reduced the zig-zag roads in the extracted road network.

Datasets and Evaluation of the Extracted Road Network
Real trajectory data of taxis in the Shenzhen city [12] were used to verify the effectiveness of our road network extraction method using low-frequency GPS trajectory data. The Shenzhen dataset contained over 75 million GPS points from 14,716 taxis, with an average sample rate of 26.1 s [26], collected in one day, which represented feasible but low-quality GPS trajectory data. In our experiments, only a subset covering an area of 2.0× 2.0 km was used for quantitative evaluation (Figure 19). The ground-truth map was provided by OpenStreetMap (OSM) [27], and we only kept roads that were travelled by at least three trajectories. The partial Shenzhen taxi trajectory data and manually processed OpenStreetMap can be found in supplementary material. To compare the extracted road network with ground-truth data, we adopted a graphsampling-based distance measure proposed by Biagioni and Erikson [3]. This measure considered both the geometric and topological similarity of maps. The sample points on the extracted road network were called marbles, and the sample points on the groundtruth road network were called holes. If a marble was close enough to a hole (judging from a predefined distance threshold), a matched marble and a matched hole were yielded. Based on this intuition, two metrics could be computed by counting the matched marbles and holes: (1) . Combing these two metrics, the F-score could be computed as: where = 1 − and = 1 − . The F-score represents the proximity of the extracted road network to the groundtruth map, under a specific distance threshold. Changing the distance threshold could reflect different standards for evaluating a generated road network, from satisfying commercial requirements to visualization needs. In our experiments, a set of distance thresholds {5 m, 10 m, 20 m} was used.

Parameter Tuning
Our method consisted of three steps-automatic seeding, road tracking, and road centerline refinement. The tracking algorithm was sensitive to the generated seed points, therefore automatic seeding is a key step, which involves two critical parameters-the circular sampling radius and the density-connectivity threshold .
The parameter tuning was based on the ground-truth intersections offered by Open-StreetMap. The generated seed points were evaluated by computing the error of road ori- To compare the extracted road network with ground-truth data, we adopted a graphsampling-based distance measure proposed by Biagioni and Erikson [3]. This measure considered both the geometric and topological similarity of maps. The sample points on the extracted road network were called marbles, and the sample points on the ground-truth road network were called holes. If a marble was close enough to a hole (judging from a predefined distance threshold), a matched marble and a matched hole were yielded. Based on this intuition, two metrics could be computed by counting the matched marbles and holes: (1) spurious = spurious_marbles/(spurious_marbles+ matched_marbles) and (2) missing = empty_holes/(empty_holes + matched_holes). Combing these two metrics, the F-score could be computed as: where precision = 1 − spurious and recall = 1 − missing.
The F-score represents the proximity of the extracted road network to the ground-truth map, under a specific distance threshold. Changing the distance threshold could reflect different standards for evaluating a generated road network, from satisfying commercial requirements to visualization needs. In our experiments, a set of distance thresholds {5 m, 10 m, 20 m} was used.

Parameter Tuning
Our method consisted of three steps-automatic seeding, road tracking, and road centerline refinement. The tracking algorithm was sensitive to the generated seed points, therefore automatic seeding is a key step, which involves two critical parameters-the circular sampling radius r and the density-connectivity threshold ϕ.
The parameter tuning was based on the ground-truth intersections offered by Open-StreetMap. The generated seed points were evaluated by computing the error of road orientation using known road branches of the ground-truth intersections. Using a fixed degree threshold 7.5 • , Table 1 shows the precision and recall of the generated seed points (before density connectivity check), using different circular sampling radius r. It could be seen that as sampling radius increased, the number of generated seed points increased, however, the recall first increased and then decreased. Both precision and recall reached its maximum value when r was 60 m. Based on this parameter setting and a loosened degree threshold of 15 • , Table 2 shows the precision and recall after the density-connectivity check, using different threshold ϕ. The density-connectivity check tried to remove the false seed points without misclassifying the true ones. About 5.2% of the generated seed points could not match with any ground-truth seed points, even when the degree threshold was greatly loosened. In this case, a connectivity threshold of 0.1 was able to remove half of the false seed points, while preventing the recall from decreasing too much.
The proposed road-tracking algorithm involved one key parameter, the tracking step length s. Since the road centerlines would be further refined in post-processing, here, we focus on the recall and complexity of the road network generated by road tracking. We use ground-truth seed points as input to exclude road tracking's sensitivity to seed points. The distance threshold in evaluation of the road network was set to 20 m, causing every extracted road be counted in the recall.
As shown in Table 3, given the correct seed points, the proposed road-tracking algorithm showed the robustness for covering the existing roads. Either too small or too large s led to redundant vertexes and edges, resulting in more spurious roads. Lastly, Table 4 lists the parameter settings for full experiments.

Evaluation of the Extracted Road Network
Since previous intersection detection methods [18,21] lacked a positioning error evaluation of the detection result, the first part of our experiments was based on the ground-truth intersections. The simulation result using intersections with different coverages and noise levels is discussed in the next section. Table 5 shows that our method achieved both relatively high precision and high recall when extracting road networks, if the input intersection was high quality. In the recent road intersection detection method proposed by Deng [21], a high precision (94.52%) was achieved using low-frequency trajectory data in Wuhan, but the recall (55.65%) was rather low. Intuitively, a recall-focus intersection detection method could better meet our method's needs to cover more roads in the extraction result.  Table 5 also shows that the quality of automatically generated road network was still inadequate for commercial applications such as navigation and automated driving, a more realistic usage could be identifying changes and missing roads in the existing map, at a daily level (noticed that the Shenzhen dataset we used was collected on one day). Figure 20a shows that our method could generate smooth road centerlines from the noisy and unevenly distributed density map. Figure 20b demonstrates the performance of geometric refinement, which shows that the refined roads achieved a balance between straight road and zig-zag roads, connecting the local peaks of the density field. However, the refined road was still susceptible to the heat zone of intersection region near its two ends.
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 20 of 25 achieved using low-frequency trajectory data in Wuhan, but the recall (55.65%) was rather low. Intuitively, a recall-focus intersection detection method could better meet our method's needs to cover more roads in the extraction result.  Table 5 also shows that the quality of automatically generated road network was still inadequate for commercial applications such as navigation and automated driving, a more realistic usage could be identifying changes and missing roads in the existing map, at a daily level (noticed that the Shenzhen dataset we used was collected on one day). Figure 20a shows that our method could generate smooth road centerlines from the noisy and unevenly distributed density map. Figure 20b demonstrates the performance of geometric refinement, which shows that the refined roads achieved a balance between straight road and zig-zag roads, connecting the local peaks of the density field. However, the refined road was still susceptible to the heat zone of intersection region near its two ends. In conclusion, quantitative evaluation and visual inspection demonstrated our method's effectiveness as a semi-automatic method, the automatic one would continuously benefit from the development of the road-intersection-detection method. However, the position of intersections was not refined by our method. There were also a small proportion of spurious roads in the extracted road network, derived from the false seed points, which required a more in-depth research on an active contour-based road network refinement.

Analysis of Template Matching
To verify the template-matching algorithm, we randomly sampled 176 slice images In conclusion, quantitative evaluation and visual inspection demonstrated our method's effectiveness as a semi-automatic method, the automatic one would continuously benefit from the development of the road-intersection-detection method. However, the position of intersections was not refined by our method. There were also a small proportion of spurious roads in the extracted road network, derived from the false seed points, which required a more in-depth research on an active contour-based road network refinement.

Analysis of Template Matching
To verify the template-matching algorithm, we randomly sampled 176 slice images of a road with ground-truth orientation offered by OpenStreetMap, and then manually labeled them as road segments (total of 111 slices), near intersection (total of 40 slices), or highly noisy area (total of 25 slices). Figure 21 shows that based on the minimummatching-distance criterion, our orientation extraction algorithm worked for the dominant road segment area; therefore, it could satisfy the need for simple correction of the roadtracking algorithm. Figure 21 also shows that adopting a circular sampling approach in the automatic seeding algorithm was crucial to preventing the areas near intersections from undermining the matching distance function. algorithm. Figure 21 also shows that adopting a circular sampling approach in the automatic seeding algorithm was crucial to preventing the areas near intersections from undermining the matching distance function.

Analysis of Automatic Seeding
First, we verified the effectiveness of density connectivity check using ground-truth intersections as input. Figure 22a shows the statistical result of matching the generated seed points with ground-truth road branches, from which we observed a proportion of obvious false seed points having an unacceptable orientation error (>60°). Figure 22b shows that after the density connectivity check, these obvious false seed points were removed, increasing the precision of the generated seed points at a small cost of recall.  Table 6 lists the detail precision, recall, and the F-score of the generated seed points (after density connectivity check), using different degree thresholds. As shown in Table 6, our seeding algorithm could balance precision and recall.

Analysis of Automatic Seeding
First, we verified the effectiveness of density connectivity check using ground-truth intersections as input. Figure 22a shows the statistical result of matching the generated seed points with ground-truth road branches, from which we observed a proportion of obvious false seed points having an unacceptable orientation error (>60 • ). Figure 22b shows that after the density connectivity check, these obvious false seed points were removed, increasing the precision of the generated seed points at a small cost of recall. algorithm. Figure 21 also shows that adopting a circular sampling approach in the automatic seeding algorithm was crucial to preventing the areas near intersections from undermining the matching distance function.

Analysis of Automatic Seeding
First, we verified the effectiveness of density connectivity check using ground-truth intersections as input. Figure 22a shows the statistical result of matching the generated seed points with ground-truth road branches, from which we observed a proportion of obvious false seed points having an unacceptable orientation error (>60°). Figure 22b shows that after the density connectivity check, these obvious false seed points were removed, increasing the precision of the generated seed points at a small cost of recall.  Table 6 lists the detail precision, recall, and the F-score of the generated seed points (after density connectivity check), using different degree thresholds. As shown in Table 6, our seeding algorithm could balance precision and recall.   Table 6 lists the detail precision, recall, and the F-score of the generated seed points (after density connectivity check), using different degree thresholds. As shown in Table 6, our seeding algorithm could balance precision and recall. Next, we analyzed the sensitivity to input intersection positions under the fixed degree threshold of 7.5 • . Since the previous road-intersection-detection methods lacked the fundamental positioning error model, we modeled the positioning error as 2D Gaussian distribution, and introduced circular error probable (CEP) to represent the noise level. When the CEP was L meters, 50% of the positioning error were within L meters. Figure 23 shows that the highly noisy intersection input eventually lowered the performance of the seeding algorithm, posing a moderate requirement on the precision of the road intersection detection method. distribution, and introduced circular error probable (CEP) to represent the noise level. When the CEP was meters, 50% of the positioning error were within meters. Figure 23 shows that the highly noisy intersection input eventually lowered the performance of the seeding algorithm, posing a moderate requirement on the precision of the road intersection detection method.

Input Sensitivity Analysis
Last, we tested the stability of our road network extraction method using simulation intersections with different coverages and noise levels (CEP). We modeled three types of input intersection-(1) ground-truth, (2) precision-focused, and (3) recall-focused. Simulation intersections were generated by randomly picking from ground-truth intersections and then adding 2D Gaussian noise to the center position; the detail CEP and recalls of different inputs are listed in Table 7. Ideal prerequisite The ground-truth input represents the manually selected input, provided by human operators or OpenStreetMap. The precision-focused input represents previous precisionfocused road intersection detection methods; here, the recall was determined from the recent work by Deng [21]. The recall-focused input represents ideal methods that could better meet our method's need.
As shown in Figure 24, using the precision-focused input, our method could more accurately generate parts of the road network. However, it would eventually miss about 20% of the existing roads. Using recall-focused input, this automatic approach could cover over 90% of the existing roads in the extraction result, regardless of the moderate intersection position error. Its ability to cover existing roads was very close to the semi-automatic approach. Finally, if the vertex of extracted road network could also be refined, which was what our method lacked, the map quality could be further improved.

Input Sensitivity Analysis
Last, we tested the stability of our road network extraction method using simulation intersections with different coverages and noise levels (CEP). We modeled three types of input intersection-(1) ground-truth, (2) precision-focused, and (3) recall-focused. Simulation intersections were generated by randomly picking from ground-truth intersections and then adding 2D Gaussian noise to the center position; the detail CEP and recalls of different inputs are listed in Table 7. The ground-truth input represents the manually selected input, provided by human operators or OpenStreetMap. The precision-focused input represents previous precisionfocused road intersection detection methods; here, the recall was determined from the recent work by Deng [21]. The recall-focused input represents ideal methods that could better meet our method's need.
As shown in Figure 24, using the precision-focused input, our method could more accurately generate parts of the road network. However, it would eventually miss about 20% of the existing roads. Using recall-focused input, this automatic approach could cover over 90% of the existing roads in the extraction result, regardless of the moderate intersection position error. Its ability to cover existing roads was very close to the semi-automatic approach. Finally, if the vertex of extracted road network could also be refined, which was what our method lacked, the map quality could be further improved.

Conclusions
In this paper, a low-frequency GPS trajectory data-based road network extraction method was proposed. The feasibility of GPS trajectory data provided a new means to automatically generating road maps. However, the prevalence of low sampling rate and multi-level disparity complicated the achievement of both high precision and high coverage in the generated map. Among various road-network-extraction methods, the intersection-first approach used the valuable information from road intersections in a divide-andconquer manner, and it could be further benefited from the development of intersectiondetection methods. Therefore, we followed this approach to propose a combined method that applied a road-tracking algorithm on the density map estimated by the KDE algorithm, where information of road intersections was represented as automatically generated seed points.
We found that the divide-and-conquer strategy was the key to tackling the multilevel disparity problem. By ruling out the near intersection region in advance, the template matching-based road-orientation-extraction algorithm performed better in general. The road-tracking algorithm could trace a road connecting two known intersections with more robustness. We pointed out that a recall-focus road-intersection-detection method could better suffice the need of road network extraction in this approach. By adopting an active contour-based road-centerline-refinement algorithm, our method yielded smoothed road centerlines than the traditional KDE-based methods. Compared to the previous GPS-trajectory-data-based road-network extraction methods, our method achieved a good balance between precision and recall using low-frequency trajectories. Moreover, the proposed road centerline refinement algorithm could serve as a general solution to zig-zag roads, considering both road structure and image information.
However, our method lacks the ability to refine the position of road intersections and to handle complicated intersection structures (e.g., highways and roundabouts). There also exists a small proportion of spurious roads due to false seed points, which requires an in-depth road-network-refinement method. Therefore, our future work will focus on applying the active contour for optimizing the position of a known intersection, and developing its ability to extract the precise boundary of city blocks from a density map, which is helpful in identifying the dead-end roads.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, the partial Shenzhen taxi trajectory data and OpenStreetMap we used in experiment is added in the supplementary file.

Conclusions
In this paper, a low-frequency GPS trajectory data-based road network extraction method was proposed. The feasibility of GPS trajectory data provided a new means to automatically generating road maps. However, the prevalence of low sampling rate and multi-level disparity complicated the achievement of both high precision and high coverage in the generated map. Among various road-network-extraction methods, the intersection-first approach used the valuable information from road intersections in a divide-and-conquer manner, and it could be further benefited from the development of intersection-detection methods. Therefore, we followed this approach to propose a combined method that applied a road-tracking algorithm on the density map estimated by the KDE algorithm, where information of road intersections was represented as automatically generated seed points.
We found that the divide-and-conquer strategy was the key to tackling the multi-level disparity problem. By ruling out the near intersection region in advance, the template matching-based road-orientation-extraction algorithm performed better in general. The road-tracking algorithm could trace a road connecting two known intersections with more robustness. We pointed out that a recall-focus road-intersection-detection method could better suffice the need of road network extraction in this approach. By adopting an active contour-based road-centerline-refinement algorithm, our method yielded smoothed road centerlines than the traditional KDE-based methods. Compared to the previous GPStrajectory-data-based road-network extraction methods, our method achieved a good balance between precision and recall using low-frequency trajectories. Moreover, the proposed road centerline refinement algorithm could serve as a general solution to zig-zag roads, considering both road structure and image information.
However, our method lacks the ability to refine the position of road intersections and to handle complicated intersection structures (e.g., highways and roundabouts). There also exists a small proportion of spurious roads due to false seed points, which requires an indepth road-network-refinement method. Therefore, our future work will focus on applying the active contour for optimizing the position of a known intersection, and developing its ability to extract the precise boundary of city blocks from a density map, which is helpful in identifying the dead-end roads.
Supplementary Materials: The following are available online at https://www.mdpi.com/2220-996 4/10/3/122/s1, the partial Shenzhen taxi trajectory data and OpenStreetMap we used in experiment is added in the supplementary file.