Hierarchical Segmentation Method for Generating Road Intersections from Crowdsourced Trajectory Data

: Maintaining the data freshness and completeness of road intersection information is the key task of urban road map production and updating. Compared to professional surveying methods, crowdsourced trajectory data provide a low-cost, wide-coverage and real-time data resource for road map construction. However, there may exist the problems of spatio-temporal heterogeneity and uneven density distribution in crowdsourced trajectory data. Hence, in light of road hierarchies, the paper proposes a hierarchical segmentation method to generate road intersections from crowdsourced trajectories. The proposed method ﬁrstly implements an adaptive density homogenization processing on raw trajectory data in order to decrease the uneven density discrepancy. Then, a hierarchical segmentation strategy is developed to extract multi-level road intersection elements from coarse scale to ﬁne scale. Finally, the structural models of road intersections are delineated by an iterative piecewise ﬁtting method. Experimental results show that the proposed method can accurately and completely extract road intersections of different shapes and scales, with an accuracy of about 87–90%. Particularly, the precision and recall of road intersection detection are obviously increased by about 7% and 20% by adaptive density homogenization, indicating the advantages of dealing with uneven trajectory data.


Introduction
High Definition (HD) routable road map is an indispensable data infrastructure for many Intelligent Transportation System (ITS) applications, such as self-driving, Intelligent Connected Vehicles (ICVs) and Mobility as a Service (MaaS) [1,2].As critical hubs of interconnecting roadways, road intersections play a significant role in route guiding [3], traffic controlling [4] and road safety [5].Generally, road intersections have more complex network structures and more diverse traffic rules than common roadways [6].Nevertheless, the road intersection information in current map products is limitedly represented on both the level of detail and data freshness, leading to hardly meeting the demands of many ITS applications [7].Hence, it becomes one of the crucial topics in the field of road map production and data updating to generate the up-to-date and high-detailed models of road intersections [8].
In the last decades, numerous targeted approaches, such as map vectorization, mobile mapping and remote sensing have been developed for commercial or official road map construction [9][10][11].However, these professional-surveying approaches may be confronted with some vulnerabilities, such as long updating cycle, outdated data source (for map vectorization), expensive data acquisition/processing cost or ticklish object-occlusion problem [9,12,13].Presently, in the new data-collection paradigm of "citizens as sensors" [14], some crowdsourced mapping platforms (e.g., OpenStreetMap) make important data supplements for map construction and updating in spite of submissions of varying quality [15,16].Particularly, the increasingly amount of crowdsourced trajectory data not only prompt the research on spatiotemporal dynamic analysis [17][18][19], but also provide a continuousdynamic, wide-coverage and low-cost way for constructing routable road map [20][21][22].However, some unpredictable bias and observation errors may occur to crowdsourced trajectories.On the one hand, limited by the equipped GNSS (Global Navigation Satellite System) devices, unstable positioning precision and uneven sampling intervals make crowdsourced trajectories possibly delineate inconsistent shapes of true moving path (see Figure 1a).On the other hand, impacted by uneven road traffic flows, the collected trajectories at different grade roads may show largely varied densities.As shown in Figure 1b, the trajectories collected at high-grade roads are generally highly dense, whereas those trajectories collected at low-grade streets are relatively sparse.
ppl.Sci.2022, 12, x FOR PEER REVIEW 2 of 16 map vectorization), expensive data acquisition/processing cost or ticklish object-occlusion problem [9,12,13].Presently, in the new data-collection paradigm of "citizens as sensors" [14], some crowdsourced mapping platforms (e.g., OpenStreetMap) make important data supplements for map construction and updating in spite of submissions of varying quality [15,16].Particularly, the increasingly amount of crowdsourced trajectory data not only prompt the research on spatiotemporal dynamic analysis [17][18][19], but also provide a continuous-dynamic, wide-coverage and low-cost way for constructing routable road map [20][21][22].However, some unpredictable bias and observation errors may occur to crowdsourced trajectories.On the one hand, limited by the equipped GNSS (Global Navigation Satellite System) devices, unstable positioning precision and uneven sampling intervals make crowdsourced trajectories possibly delineate inconsistent shapes of true moving path (see Figure 1a).On the other hand, impacted by uneven road traffic flows, the collected trajectories at different grade roads may show largely varied densities.As shown in Figure 1b, the trajectories collected at high-grade roads are generally highly dense, whereas those trajectories collected at low-grade streets are relatively sparse.To address the above-mentioned challenges, the paper proposed a hierarchical segmentation method for generating road intersections from crowdsourced trajectories.The main contributions of the proposed method can be summarized as the following aspects.
 An adaptive density homogenization operator is performed on raw track data to decrease the uneven density distribution of crowdsourced trajectories. A hierarchical segmentation process is implemented to targetedly extract multilevel road intersection elements implied on crowdsourced trajectory data. An iterative piecewise fitting method is adopted to approximately delineate the geometric structures of road intersections.
The remainder of the paper is organized as follows.The state of art road network extraction is systematically reviewed in Section 2. After that, the road intersection extraction method is detailedly elaborated in Section 3. Experimental results and analysis are then presented in Section 4. Finally, Section 5 concludes the paper and outlooks the future works.

Literature Review
In the early stage, optical remote sensing images provided the primary data sources for fast acquiring road network data.The existing methods can mainly be categorized into four kinds, e.g., road detection methods based on template matching [23,24], road extraction models based on prior knowledge [25], object-oriented road segmentation methods [26] and road semantic labelling methods based on deep learning [13].Generally, the image-based road extraction methods have the advantages of high efficiency and mature To address the above-mentioned challenges, the paper proposed a hierarchical segmentation method for generating road intersections from crowdsourced trajectories.The main contributions of the proposed method can be summarized as the following aspects.

•
An adaptive density homogenization operator is performed on raw track data to decrease the uneven density distribution of crowdsourced trajectories.

•
A hierarchical segmentation process is implemented to targetedly extract multi-level road intersection elements implied on crowdsourced trajectory data.

•
An iterative piecewise fitting method is adopted to approximately delineate the geometric structures of road intersections.
The remainder of the paper is organized as follows.The state of art road network extraction is systematically reviewed in Section 2. After that, the road intersection extraction method is detailedly elaborated in Section 3. Experimental results and analysis are then presented in Section 4. Finally, Section 5 concludes the paper and outlooks the future works.

Literature Review
In the early stage, optical remote sensing images provided the primary data sources for fast acquiring road network data.The existing methods can mainly be categorized into four kinds, e.g., road detection methods based on template matching [23,24], road extraction models based on prior knowledge [25], object-oriented road segmentation methods [26] and road semantic labelling methods based on deep learning [13].Generally, the image-based road extraction methods have the advantages of high efficiency and mature technologies, which are suitable for large-scale scenarios.However, the image-based methods may suffer from some difficulties in approaching poor illumination, severe weather condition and ground object occlusion, as well as obtaining dynamic traffic information.
After that, Light Detection And Ranging (LiDAR) technologies produced a new data source for extracting high detailed road features, such as 3D road boundary, road markings and affiliated road facilities (e.g., street lamp, traffic signs and street trees).The road boundary extraction methods were mostly based on the assumption that the road surfaces are approximately flat and sharp elevation fluctuations probably occur at road curbs [11,27].The road marking extraction methods generally utilized the intensity difference of point clouds and specific shapes of road markings to obtain typical road marking information [28,29].Additionally, numerous studies attempted to apply template matching, super voxel and deep learning in extracting pole-like road facilities [30,31].In brief, the LiDARbased methods can extract highly precise and very detailed road features under different environmental conditions, but meanwhile, high data cost together with the vulnerable occlusion problem may limit their application in large-scale areas.
Presently, the emerging of crowdsourced mapping (i.e., Volunteer Geographic Information, VGI) prompts a paradigm shift from professional surveying to ubiquitous mapping [32].As a low-cost, wide-coverage and real-time data resource, crowdsourced trajectories collected by vehicles or pedestrians can be processed to extract road networks [20,[33][34][35][36].The trajectory-based road extraction approaches can mainly be summarized as rasterization methods and vectorization methods.The rasterization methods first converted vector trajectories to raster images and then made use of various image processing algorithms, such as morphological operators and skeletonization algorithms, to generate vector road networks [37,38].The rasterization methods are proved efficient and suitable to handle wide-range construction of road maps, but it may be challenging to obtain road semantic information and topological relations for their neglecting the trajectory movement characteristics [39].The vectorization methods mainly consist of two kinds of strategies.One was devoted to automated road map construction by exploring the clustering patterns of moving trajectories [40], and the other one adopted a trajectory-aligning strategy to incrementally construct road networks [41,42].Generally, the vectorization methods can generate detailed road structures and infer some road semantic information (e.g., road width and turning rules); however, these methods are still not robust enough to deal with large trajectory density disparities.Particularly, Zhang et al. (2020) proposed a hybrid road extraction method by considering the road changing patterns implied in multi-temporal trajectory data [12].Furthermore, some studies attempted to construct lane-level road maps from high frequency/precision trajectory data.Most of them first computed the lateral distances of trajectory points to road central lines and inferred the lane number or location from the distance distributions based on the Gaussian Mixture Model (GMM), Bayesian Classifier and Kernel Density Estimation (KDE) [1,[43][44][45].
Particularly, complete road intersection information is very important for route planning and traffic management.Hence, more and more attention has been paid to extracting road intersections from crowdsourced trajectory data in recent years [4,46].The related methods often supposed road intersections as the convergent places of two or more roadways, where vehicle trajectories probably manifest some specific driving behaviors, such as changing lane, changing speed and steering off [7].Based on such driving characteristics, the general workflow of road intersection extraction consisted of three main steps, i.e., road intersections detection, turning modes identification and road geometries derivation [47,48].Fathi and Krumm (2010) developed a trace classifier trained by localized shape descriptors to distinguish road intersections from non-intersections [49].Wang et al. (2015) introduced G-statistics hotspot analysis to mine the agglomeration patterns of large turning behaviors at road intersections [3].Xie et al. (2017) computed the longest common subsequence between two GNSS traces to detect the connecting points where the candidates of road intersections [50].Recently, Yang et al. (2021) applied a mask region convolutional neural network for automatically detecting road intersections [51].The turning modes within road intersections can be further identified based on trajectory clustering [48] or entrances/exits matching [46,52].Based on turning modes identification, two widely-used strategies of curve fitting [47] and piecewise approximation [1] are adopted to derive the road central lines of different turning modes.Although previous methods [44][45][46][47][48][49][50][51] made much beneficial work to obtain road intersection information, there are still some improvement spaces in the robustness to approach heterogeneous trajectory data and the reliability to extract road intersections with different shapes and structures.The paper thus implements a density homogenization process for the purposes of alleviating raw trajectories' density unevenness, and meanwhile, decreases huge computational cost of processing all raw data.Particularly, in light of road hierarchies, a hierarchical segmentation strategy is proposed to extract multi-level road intersection elements from the heterogeneous trajectory data.

Materials and Methods
As illustrated in Figure 2, the proposed method mainly consists of three key steps: adaptive density homogenization, hierarchical trajectory segmentation and road structure delineation.In detail, a data preprocessing is first performed on raw trajectory data to filter out the abnormal points and traces.Then, an adaptive density homogenization is utilized to reduce the impacts of various traffic volume on trajectory sampling density.After that, a hierarchical trajectory segmentation strategy is designed to gradually extract multi-level road intersection elements according to their direction and position characteristics.Finally, the geometric structures of road intersections are delineated by an iterative piecewise fitting algorithm.
turning modes within road intersections can be further identified based on trajectory clustering [48] or entrances/exits matching [46,52].Based on turning modes identification, two widely-used strategies of curve fitting [47] and piecewise approximation [1] are adopted to derive the road central lines of different turning modes.Although previous methods [44][45][46][47][48][49][50][51] made much beneficial work to obtain road intersection information, there are still some improvement spaces in the robustness to approach heterogeneous trajectory data and the reliability to extract road intersections with different shapes and structures.The paper thus implements a density homogenization process for the purposes of alleviating raw trajectories' density unevenness, and meanwhile, decreases huge computational cost of processing all raw data.Particularly, in light of road hierarchies, a hierarchical segmentation strategy is proposed to extract multi-level road intersection elements from the heterogeneous trajectory data.

Materials and Methods
As illustrated in Figure 2, the proposed method mainly consists of three key steps: adaptive density homogenization, hierarchical trajectory segmentation and road structure delineation.In detail, a data preprocessing is first performed on raw trajectory data to filter out the abnormal points and traces.Then, an adaptive density homogenization is utilized to reduce the impacts of various traffic volume on trajectory sampling density.After that, a hierarchical trajectory segmentation strategy is designed to gradually extract multi-level road intersection elements according to their direction and position characteristics.Finally, the geometric structures of road intersections are delineated by an iterative piecewise fitting algorithm.

Iterative piecewise fitting
Figure 2. The workflow of the proposed method.

Adaptive Density Homogenization
As listed in Table 1, raw trajectory data are in .txtdata format, which records the car ID, time, longitude, latitude and speed for each trajectory point.The trajectory of one vehicle is reorganized as a sequential trajectory point set grouped by car ID and sorted by time.Limited by the equipped GNSS devices, unstable positional drift and varying sampling intervals often occurred to raw trajectory points.Hence, on the one hand, the trajectory points, which have no displacements with the former and next points, are first removed for avoiding unnecessary calculations.On the other hand, the vehicle trajectories, with some sub-segments (i.e., the distance between two successive points) longer than the predefined threshold, will be split at those sub-segments.In this paper, we refer to Deng et al. [44] and set the threshold as 80km/h*30s by considering the most-frequent sampling

Adaptive Density Homogenization
As listed in Table 1, raw trajectory data are in .txtdata format, which records the car ID, time, longitude, latitude and speed for each trajectory point.The trajectory of one vehicle is reorganized as a sequential trajectory point set grouped by car ID and sorted by time.Limited by the equipped GNSS devices, unstable positional drift and varying sampling intervals often occurred to raw trajectory points.Hence, on the one hand, the trajectory points, which have no displacements with the former and next points, are first removed for avoiding unnecessary calculations.On the other hand, the vehicle trajectories, with some sub-segments (i.e., the distance between two successive points) longer than the predefined threshold, will be split at those sub-segments.In this paper, we refer to Deng et al. [44] and set the threshold as 80 km/h × 30 s by considering the most-frequent sampling interval of track points and the highest speed limit in urban areas.After that, we select those trajectories which have more than 5 track points to extract road intersections in view of data reliability.
Moreover, as mentioned before, vehicle trajectory data are, actually, a sparse sensing way of road information due to uneven traffic flows.Generally, dense trajectories are obtained on high-grade roads and large road intersections while sparse trajectories are collected at secondary roads and small street blocks.Hence, we propose an adaptive density homogenization method to decrease the trajectory density discrepancies at different grade roads.The adaptive density homogenization process encompasses three key steps, described as follows.As mentioned before, vehicle trajectories probably manifest some specific behavior characteristics, such as changing lane, changing speed and steering off.However, due to uneven sampling density of raw trajectory data, the number of trajectories with such specific moving characteristics may greatly vary at different regions, which could make it difficult to make a uniform rule or threshold for detecting road intersection points from raw trajectory data.To acquire a density-even trajectory data, we should first recognize high-density trajectory cells and do a screening processing to ensure the uniform sampling in the whole study area.So, we partition the covering area of raw trajectory data into regular cells with the fixed size of w * h.The cell width w and cell height h are calculated by the following formula.
where dx i and dy i denote the absolute values of X and Y coordinate differences between two successive track points.k is a density parameter, uniformly set as 1000 here.
Then, the proportion of track points in each cell G j is computed as follows.
where Num(j) is the number of trajectory points in the cell G j .The cell G j will be identified as a high-density trajectory cell (renumbered as HG j ) when Prop j is larger than a predefined threshold N grid .

Calculating Neighborhood Circles for High-Density Trajectory Cell
The identified high-density cells indicate that there are a large number of trajectories collected at these cells.However, the covering range of high-density trajectories at each high-density cell is uncertain in spite of having a fixed cell size.To determine the covering range of high-density trajectories, we compute a neighborhood circle for each high-density cell.As shown in Figure 3, for each HG j , we suppose its center point is c j .Then, the neighborhood circle of HG j is calculated as follows.
Step1: Build the buffer region for center c j by an increasing radius ε i = ε 0 + i * s, where the initial radius ε 0 and the step length s are both set as 10 m and i denotes the iterative number.
Step2: Count the ratio of trajectory points in two adjacent buffer regions as Formula (3), where Num(ε i ) is the number of trajectory points falling in the ε i -buffer region.
Step3: Construct the neighborhood circle for HG j using c j and the circle radius r j defined as Formula (4), where N circle is a given threshold and the maximum buffer radius is set as 200 m because of the maximum road intersection size.
where First{ε i . . .} means the first element in the set {ε i . . .}. Step1: Build the buffer region for center cj by an increasing radius εi=ε0 +i*s, where the initial radius ε0 and the step length s are both set as 10m and i denotes the iterative number.
Step2: Count the ratio of trajectory points in two adjacent buffer regions as formula (3), where Num(εi) is the number of trajectory points falling in the ε  -buffer region.
Step3: Construct the neighborhood circle for HGj using cj and the circle radius rj defined as formula (4), where Ncircle is a given threshold and the maximum buffer radius is set as 200m because of the maximum road intersection size.
where First{εi…} means the first element in the set {εi…}.

Trajectory Density Homogenization of High-Density Trajectory Cell
To retain as much road changing information implied in raw trajectory data as possible, we select ten-days trajectory data (e.g., 1 st , 4 th , …, 28 th ) at equal intervals from the raw one-month trajectory data and perform a density homogenization processing on the remained trajectory data.In detail, the density homogenization processing is presented as follows.
Step1: The remaining trajectories are traversed one by one to judge whether it falls in the neighborhood circle of a high-density cell.
Step2: If the trajectory does not completely fall in the neighbor circle of any highdensity cell, it will be reserved and return Step1 to traverse the next trajectory; otherwise, the original trajectory will be split into several trajectory segments by the neighborhood circles overlaying with it.
Step3: Among the trajectory segments split by Step2, those trajectory segments falling in neighborhood circles will be removed and other ones will be reserved, return Step1 to traverse the next trajectory.
Finally, the initially-selected ten-days trajectory data together with those reserved trajectory data after density homogenization are used in the following step of categorylevel segmentation.

Hierarchical Segmentation Model for Crowdsourced Trajectories
The trajectory data in road intersection areas may show specific movement patterns such as changing lane, changing speed and steering off.Based on this, many studies identified the trajectory segments of lane changing or the trajectory points with large turns by hotspot analysis [48], decision tree learning [53] or dynamic programming [50].These

Trajectory Density Homogenization of High-Density Trajectory Cell
To retain as much road changing information implied in raw trajectory data as possible, we select ten-days trajectory data (e.g., 1st, 4th, . . ., 28th) at equal intervals from the raw onemonth trajectory data and perform a density homogenization processing on the remained trajectory data.In detail, the density homogenization processing is presented as follows.
Step1: The remaining trajectories are traversed one by one to judge whether it falls in the neighborhood circle of a high-density cell.
Step2: If the trajectory does not completely fall in the neighbor circle of any highdensity cell, it will be reserved and return Step1 to traverse the next trajectory; otherwise, the original trajectory will be split into several trajectory segments by the neighborhood circles overlaying with it.
Step3: Among the trajectory segments split by Step2, those trajectory segments falling in neighborhood circles will be removed and other ones will be reserved, return Step1 to traverse the next trajectory.
Finally, the initially-selected ten-days trajectory data together with those reserved trajectory data after density homogenization are used in the following step of categorylevel segmentation.

Hierarchical Segmentation Model for Crowdsourced Trajectories
The trajectory data in road intersection areas may show specific movement patterns such as changing lane, changing speed and steering off.Based on this, many studies identified the trajectory segments of lane changing or the trajectory points with large turns by hotspot analysis [48], decision tree learning [53] or dynamic programming [50].These methods are restrictedly applicable for high-quality trajectory data, and furthermore, may result in the trajectories at gas stations, parking lots and residential areas being misidentified as the representative trajectories in road intersection areas because the trajectories at these regions show similar turning behaviors with those at road intersections.As previous studies examined [54], urban road networks are characterized by hierarchical organization, which provides an important inspiration for extracting road intersections from trajectory data.The key idea of generating road intersections from trajectory data is to hierarchically segment massive heterogeneous trajectory points or lines into different clusters, and then geometrically delineate the trajectory points or lines of each cluster so as to form graphical structures of road intersections.Therefore, on the basis of trajectory density homogenization, this paper proposes a hierarchical segmentation method to identify multi-level road intersection elements from coarse scale to fine scale.In detail, the first category-level seg-mentation is to distinguish the trajectory points of road-intersection and non-intersection according to the indicators of speed change, direction change and turning distance ratio.
The second object-level segmentation is to use Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm to further partition the trajectory points, previously, categorized as road intersection into different clusters that correspond to various road intersection objects.Finally, the third turning-modes segmentation is to extract the internal traffic rules by analyzing the starting and ending direction of each trajectory segments within road intersection.

Category-Level Segmentation
Previous studies proposed to use the velocity change, direction change and straight curve ratio of trajectory points to recognize the trajectory points at road intersections.However, vehicles that pass through residential areas, gas stations and parking lots may have the similar moving characteristics with those passing through road intersections.Hence, based on the trajectory data after adaptive density homogenization, we try to use three indicators, namely speed change ∆v, direction change ∆θ and turning distance ratio ∆R to distinguish the road-intersection trajectory points from non-intersection trajectory points.For one trajectory consisting of m points, i.e., T = {p 1 . . .p m }, the speed change ∆v i and direction change ∆θ i of one point p i are calculated according to Formula (5), respectively.
Particularly, the turning distance ratio ∆R is an indicator to measure the turning degree of one trajectory point, and its calculation process is shown in Figure 4. Firstly, the trajectory T = {p 1 . . .p m } was segmented as a bend hierarchical tree according to the method of [55], and the turning distance h i of each trajectory point p i and the baseline length bend i of the bend segment were recorded, respectively.As shown in Figure 4, for trajectory point p 2, the baseline length of its belonging bend segment was the Euclidean distance between p 1 and p 10 , i.e., bend 2 = |p 1 p 10 |, and the turning distance h 2 of trajectory point p 2 is equal to the perpendicular distance from p 2 to the baseline p 1 p 10 .Then, the turning distance ratio of trajectory point p 2 was calculated as ∆R 2 = h 2 /bend 2 , so did other trajectory points.When the three indicators of speed change ∆v, direction change ∆θ and turning distance ratio ∆R are all greater than the given thresholds T v , T α and T dr , the track point will be identified as a candidate of road-intersection trajectory point.from trajectory data.The key idea of generating road intersections from trajectory data is to hierarchically segment massive heterogeneous trajectory points or lines into different clusters, and then geometrically delineate the trajectory points or lines of each cluster so as to form graphical structures of road intersections.Therefore, on the basis of trajectory density homogenization, this paper proposes a hierarchical segmentation method to identify multi-level road intersection elements from coarse scale to fine scale.In detail, the first category-level segmentation is to distinguish the trajectory points of road-intersection and non-intersection according to the indicators of speed change, direction change and turning distance ratio.The second object-level segmentation is to use Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm to further partition the trajectory points, previously, categorized as road intersection into different clusters that correspond to various road intersection objects.Finally, the third turning-modes segmentation is to extract the internal traffic rules by analyzing the starting and ending direction of each trajectory segments within road intersection.

Category-Level Segmentation
Previous studies proposed to use the velocity change, direction change and straight curve ratio of trajectory points to recognize the trajectory points at road intersections.However, vehicles that pass through residential areas, gas stations and parking lots may have the similar moving characteristics with those passing through road intersections.Hence, based on the trajectory data after adaptive density homogenization, we try to use three indicators, namely speed change ∆v, direction change ∆θ and turning distance ratio ∆R to distinguish the road-intersection trajectory points from non-intersection trajectory points.For one trajectory consisting of m points, i.e., T = {p1…pm}, the speed change ∆vi and direction change ∆θi of one point pi are calculated according to formula (5), respectively.
Particularly, the turning distance ratio ∆R is an indicator to measure the turning degree of one trajectory point, and its calculation process is shown in Figure 4. Firstly, the trajectory T= {p1…pm} was segmented as a bend hierarchical tree according to the method of [55], and the turning distance hi of each trajectory point pi and the baseline length bendi of the bend segment were recorded, respectively.As shown in Figure 4, for trajectory point p2, the baseline length of its belonging bend segment was the Euclidean distance between p1 and p10, i.e., bend2 = |p1p10|, and the turning distance h2 of trajectory point p2 is equal to the perpendicular distance from p2 to the baseline p1p10.Then, the turning distance ratio of trajectory point p2 was calculated as ∆R2= h2/bend2, so did other trajectory points.When the three indicators of speed change ∆v, direction change ∆θ and turning distance ratio ∆R are all greater than the given thresholds Tv, Tα and Tdr, the track point will be identified as a candidate of road-intersection trajectory point.As illustrated in Figure 5, take the blue trace for example.The identified roadintersection trajectory points based on speed change, direction change and turning distance ratio are marked as red dots in Figures 5a, 5b and 5c, respectively.It can be seen that only using the indicator of speed change or direction change may lead to misidentifying some non-intersection trajectory points that suddenly change the speed or direction (see the cyan ellipse in Figure 5) and some local-turning trajectory points that pass through residential areas, gas stations and parking lots (see the cyan rectangle in Figure 5) as the candidates of road-intersection trajectory points.Comparatively, the turning distance ratio can efficiently distinguish irregular-turning trajectory points at specific places (e.g., gas stations) and regular-turning trajectory points at actual intersections.
As illustrated in Figure 5, take the blue trace for example.The identified road-inter-section trajectory points based on speed change, direction change and turning distance ratio are marked as red dots in Figure 5a, Figure 5b and Figure 5c, respectively.It can be seen that only using the indicator of speed change or direction change may lead to misidentifying some non-intersection trajectory points that suddenly change the speed or direction (see the cyan ellipse in Figure 5) and some local-turning trajectory points that pass through residential areas, gas stations and parking lots (see the cyan rectangle in Figure 5) as the candidates of road-intersection trajectory points.Comparatively, the turning distance ratio can efficiently distinguish irregular-turning trajectory points at specific places (e.g., gas stations) and regular-turning trajectory points at actual intersections.

Object-Level Segmentation
The category-level segmentation of trajectory data is to extract the coarse information about whether the trajectory points are road intersection or non-intersection.The next step is to perform object-level segmentation on the candidate trajectory points of road-intersection, so as to determine the location and spatial coverage of each road intersection.As the DBSCAN algorithm does not need to set the number of clusters in advance, and has certain advantages of finding different shaped clusters and handling noise points, this paper uses the DBSCAN algorithm [56] to further segment road-intersection trajectory points.Suppose the road-intersection trajectory points are P = {p1…pn}, the object-level segmentation process based on two parameters of Eps and MinPts is carried out as follows.
Step1: Calculate the neighborhood point set of each point pi in P as is in formula (6), where Dis(pi, pj) is the Euclidean distance between pi and pj, and then choose those roadintersection trajectory points with |NEps(pi)| larger than MinPts as the core seed point.
Step2: Randomly select one unprocessed core seed point and establish a new cluster for it, and meanwhile, add the neighborhood points of the current core seed point to the newly-created cluster.
Step3: Sequentially traverse each unprocessed core seed point in the current cluster, and add the road-intersection trajectory points density-connected with the core seed point to the current cluster until all core seed points in the current cluster are processed.
Step4: Return to Step2 to traverse the remaining unprocessed core seed points until all core seed points have been processed.
Based on DBSCAN, each obtained cluster C of road-intersection trajectory points actually corresponds to a road intersection object.After that, the location of the corresponding road intersection is defined as the cluster center c and the spatial coverage of the corresponding road intersection is calculated as the maximum distance from the cluster points to the cluster center according to formula (7).Finally, we select the trajectory Legend Road-intersection   Example trace Non-intersection

Object-Level Segmentation
The category-level segmentation of trajectory data is to extract the coarse information about whether the trajectory points are road intersection or non-intersection.The next step is to perform object-level segmentation on the candidate trajectory points of roadintersection, so as to determine the location and spatial coverage of each road intersection.As the DBSCAN algorithm does not need to set the number of clusters in advance, and has certain advantages of finding different shaped clusters and handling noise points, this paper uses the DBSCAN algorithm [56] to further segment road-intersection trajectory points.Suppose the road-intersection trajectory points are P = {p 1 . . .p n }, the object-level segmentation process based on two parameters of Eps and MinPts is carried out as follows.
Step1: Calculate the neighborhood point set of each point p i in P as is in Formula ( 6), where Dis(p i , p j ) is the Euclidean distance between p i and p j , and then choose those roadintersection trajectory points with |N Eps (p i )| larger than MinPts as the core seed point.
Step2: Randomly select one unprocessed core seed point and establish a new cluster for it, and meanwhile, add the neighborhood points of the current core seed point to the newly-created cluster.
Step3: Sequentially traverse each unprocessed core seed point in the current cluster, and add the road-intersection trajectory points density-connected with the core seed point to the current cluster until all core seed points in the current cluster are processed.
Step4: Return to Step2 to traverse the remaining unprocessed core seed points until all core seed points have been processed.
Based on DBSCAN, each obtained cluster C of road-intersection trajectory points actually corresponds to a road intersection object.After that, the location of the corresponding road intersection is defined as the cluster center c and the spatial coverage of the corresponding road intersection is calculated as the maximum distance from the cluster points to the cluster center according to Formula (7).Finally, we select the trajectory segments in the spatial coverage of each road intersection as the input of the following step of turning-modes segmentation.d = max

Turning-Modes Segmentation
According to the importance assessment of trajectory features conducted by the literature [57], where in the headings angle of starting and ending points are better than the heading difference of starting points, the heading difference of ending points and the average heading angle of intermediate points to express the road-intersection trajectory characteristics.Therefore, this paper uses the headings angle of starting and ending points to extract the internal traffic mode within road intersection.The specific process of turning-modes segmentation is presented as follows: Step1: For each trajectory segment of one road intersection, we first partition the heading angle of its starting or ending point into an eight directions spider code [58].As illustrated in Figure 6, when the heading angle of one starting or ending point belongs to [337.5 • ,360 • ] ∪ [0 • ,22.5 • ), the heading angle is coded as H, or when the heading angle is between [22.5 • ,67.5 • ), it is coded as A, and so on.

Turning-Modes Segmentation
According to the importance assessment of trajectory features conducted by the literature [57], where in the headings angle of starting and ending points are better than the heading difference of starting points, the heading difference of ending points and the average heading angle of intermediate points to express the road-intersection trajectory characteristics.Therefore, this paper uses the headings angle of starting and ending points to extract the internal traffic mode within road intersection.The specific process of turning-modes segmentation is presented as follows: Step1: For each trajectory segment of one road intersection, we first partition the heading angle of its starting or ending point into an eight directions spider code [58].As illustrated in Figure 6, when the heading angle of one starting or ending point belongs to [337.5°,360°]∪[0°,22.5°), the heading angle is coded as H, or when the heading angle is between [22.5°,67.5°), it is coded as A, and so on.Step2: Based on the above eight directions spider code, each trajectory segment likely has 8×8=64 heading code combinations of the starting and ending points, which can be further judged to belong to an element in the heading code matrix shown in Figure 7a.
Step3: According to the heading code matrix, we calculate the proportion of trajectory segments belonging to each matrix element and select the trajectory segments in those elements with the proportion larger than a threshold Ratio as one typical turning mode.
As shown in Figure 7a, the circle size denotes the number of trajectory segments belonging to each matrix element, and the black rectangles highlight the matrix elements with the trajectory proportion larger than the given threshold.Particularly, Figure 7b separately elaborates the trajectory segments of four typical turning modes in the rectangles (i)-(iv) of Figure 7a, which correctly correspond to the actual traffic modes at the road intersection.Step2: Based on the above eight directions spider code, each trajectory segment likely has 8 × 8 = 64 heading code combinations of the starting and ending points, which can be further judged to belong to an element in the heading code matrix shown in Figure 7a.
Step3: According to the heading code matrix, we calculate the proportion of trajectory segments belonging to each matrix element and select the trajectory segments in those elements with the proportion larger than a threshold Ratio as one typical turning mode.
As shown in Figure 7a, the circle size denotes the number of trajectory segments belonging to each matrix element, and the black rectangles highlight the matrix elements with the trajectory proportion larger than the given threshold.Particularly, Figure 7b separately elaborates the trajectory segments of four typical turning modes in the rectangles (i)-(iv) of Figure 7a, which correctly correspond to the actual traffic modes at the road intersection.

Road Structure Delineation
As shown in Figure 8, for the trajectory segment set of one turning mode TS = {Ti|Ti = {p1 (i) …pti (i) }, i=1…n}, the road central line of the turning mode is delineated by the following iterative piecewise fitting method: Step1: Select the longest trajectory segment in TS as the initial central line TC = {C1…Cm}.

Road Structure Delineation
As shown in Figure 8, for the trajectory segment set of one turning mode TS = {T i |T i = {p 1 (i) . . .p ti (i) }, i = 1. . .n}, the road central line of the turning mode is delineated by the following iterative piecewise fitting method: Step1: Select the longest trajectory segment in TS as the initial central line TC = {C 1 . . .C m }.
Step2: Build the buffer regions G Step3: Search the trajectory points falling in each buffer region, and calculate their average geometric centers as the new positions of C 1 . . .C m to update the current central line.
Step4: Repeat Step2 and Step3 until the Hausdorff distance of the fitted central lines between two iterations is less than a given threshold ε.

Road Structure Delineation
As shown in Figure 8, for the trajectory segment set of one turning mode TS = {Ti|T = {p1 (i) …pti (i) }, i=1…n}, the road central line of the turning mode is delineated by the follow ing iterative piecewise fitting method: Step1: Select the longest trajectory segment in TS as the initial central line TC = {C1…Cm}.
Step3: Search the trajectory points falling in each buffer region, and calculate thei average geometric centers as the new positions of C1...Cm to update the current central line Step4: Repeat Step2 and Step3 until the Hausdorff distance of the fitted central line between two iterations is less than a given threshold ε.

Experimental Data and Parameter Setting
In order to verify the effectiveness of the proposed method, we used trajectory data from Changsha City, Hunan Province, collected in October 2018 for comprehensive anal ysis.The used trajectory data covers an area of about 7.6×7.0km, which is characterized by heavy traffic flow, rich commercial facilities and a variety of road intersections with different shapes and complex structures.The experimental dataset contains a total o

Experimental Data and Parameter Setting
In order to verify the effectiveness of the proposed method, we used trajectory data from Changsha City, Hunan Province, collected in October 2018 for comprehensive analysis.The used trajectory data covers an area of about 7.6 × 7.0 km, which is characterized by heavy traffic flow, rich commercial facilities and a variety of road intersections with different shapes and complex structures.The experimental dataset contains a total of 14,839,457 trajectory points collected by 5439 vehicles.The positional precision is about 5-10 m and the most-frequent sampling interval is about 30 s.The experimental parameters for density homogenization, hierarchical trajectory segmentation and road structure delineation are listed in Table 2. N grid is the proportion threshold for identifying high-density trajectory cells (Section 3.1.1).N circle is the proportion threshold for calculating the neighborhood circles of high-density cells (Section 3.1.2).T v , T α and T dr are the thresholds of speed change, direction change and turning distance ratio for identifying road-intersection trajectory points (Section 3.2.1).Eps and Minpts are the clustering parameters of DBSCAN for detecting object-level road intersections (Section 3.2.2).Ratio is the proportion threshold for determining the turning modes within road intersection (Section 3.2.3).ε is used to terminate the iterative piecewise fitting process for delineating road intersection structures (Section 3.3).These parameters are mostly determined through multiple experiments, of which some parameters are empirically predefined according to the actual road width (10-70 m) and trajectory data quality (positional precision: 5-10 m, time sampling: 30 s).These used parameters could be resized for different data sources or study areas.

Result Analysis of Road Intersection Detection
Based on the density homogenization and coarse-to-fine trajectory segmentation strategy, we hierarchically extracted the trajectory points of road-intersection, the spatial coverages of different road intersections and the internal turning modes within road intersections.Figure 9a,b illustrates the results of road intersections detection with and without adaptive density homogenization.To evaluate the accuracy of road intersection detection, we referred to the national platform for common geospatial information services (https://map.tianditu.gov.cn/) and checked the identified road intersections one by one.By manual inspection, the road intersections correctly identified by the proposed method are classified as true positive (red triangle), the road intersections wrongly identified by the proposed method are classified as false positive (green circle) and, meanwhile, the road intersection not identified by the proposed method are marked as false negative (blue rectangle).
tion threshold for determining the turning modes within road intersection (Section 3.2.3).ε is used to terminate the iterative piecewise fitting process for delineating road intersection structures (Section 3.3).These parameters are mostly determined through multiple experiments, of which some parameters are empirically predefined according to the actual road width (10m-70m) and trajectory data quality (positional precision: 5m-10m, time sampling: 30s).These used parameters could be resized for different data sources or study areas.Based on the density homogenization and coarse-to-fine trajectory segmentation strategy, we hierarchically extracted the trajectory points of road-intersection, the spatial coverages of different road intersections and the internal turning modes within road intersections.Figure 9a and Figure 9b illustrates the results of road intersections detection with and without adaptive density homogenization.To evaluate the accuracy of road intersection detection, we referred to the national platform for common geospatial information services (https://map.tianditu.gov.cn/) and checked the identified road intersections one by one.By manual inspection, the road intersections correctly identified by the proposed method are classified as true positive (red triangle), the road intersections wrongly identified by the proposed method are classified as false positive (green circle) and, meanwhile, the road intersection not identified by the proposed method are marked as false negative (blue rectangle).It can be seen from Figure 9 that more correctly-detected and fewer wrongly-identified road intersections are found in the result of road intersection detection with adaptive density homogenization.In particular, as shown in the two cyan ellipses B of Figure 9, through trajectory density homogenization, the proposed method identifies some road intersections with sparse trajectory density, indicating our method can efficiently overcome large, varied trajectory densities caused by uneven traffic flow.In addition, as depicted in the right of cyan rectangle A of Figure 9a, local movement indicators, such as speed change and direction change, may misidentify some trajectory points passing through parking lots and residential areas as the trajectory points at road intersection.We introduce the turning distance ratio to distinguish the trajectory points of road intersection and nonintersection, appropriately avoiding the road intersection misidentification based on local movement indicators.
After manual inspection, three accuracy indicators of Precision, Recall and F-Value are calculated to quantitatively evaluate the performance of road intersection detection, which is stated in Table 3.It was found that 132 road intersections were identified when directly using raw trajectory data without density homogenization, of which 109 were correctly identified and 23 were wrongly identified, as well as 50 road intersections which were not detected.Overall, Precision was about 82.6%, Recall was about 68.5%, and F-Value was about 74.9%.Through trajectory density homogenization, 166 road intersections were identified, of which 149 were correctly identified, 17 were incorrectly identified as well as only 15 missed road intersections.It can be seen that the accuracy indicators of road intersection detection by trajectory density homogenization were significantly improved, of which Precision increased to 89.8%, Recall increased to 90.8%, and F-Value reached 90.3%.

Result Analysis of Road Intersection Extraction
Figure 10 shows the results of internal turning mode recognition and road intersection structure extraction.The red lines in Figure 10a are the structural models of road intersection extracted by our method.Although the study area has various kinds of road intersections with different shapes and structures, the proposed method can correctly extract complex and diverse road intersection structures based on density homogenization and hierarchical segmentation.In particular, Figure 10b, Figure 10d and Figure 10c, Figure 10e illustrate the enlarged views of the generated geometric models and extracted turning modes of the two road intersections in the green and blue circles in Figure 10a, respectively.It can be seen that the generated geometric models and extracted turning modes have high consistency with the actual remote sensing images, demonstrating the high efficiency and accuracy of the hierarchical segmentation method on extracting road intersection information for routable road map construction.
Furthermore, turning mode segmentation is the key step of road intersection extraction.To quantitatively evaluate the effectiveness of turning mode segmentation, we randomly selected 16 road intersections (including five cloverleaf intersection, five cross intersection, three T-shaped junctions and three Y-shaped junctions) and manually counted the internal turning modes within them for the ground truth.By comparing with the ground truth of manual recognition, the accuracy statistics results based on formulas (8) are also listed in the last row of Table 3, showing that the Precision indicator of turning modes recognition is about 85.7%, the recall indicator is about 88.4%, and F value is about 87.0%.That further indicates that the proposed method by combining trajectory density homogenization and hierarchical trajectory segmentation can efficiently deal with spatio-temporal heterogeneity and uneven density distribution of crowdsourced trajectory data and achieve good performance on generating structural road intersection map.
(8) are also listed in the last row of Table 3, showing that the Precision indicator of turning modes recognition is about 85.7%, the recall indicator is about 88.4%, and F value is about 87.0%.That further indicates that the proposed method by combining trajectory density homogenization and hierarchical trajectory segmentation can efficiently deal with spatiotemporal heterogeneity and uneven density distribution of crowdsourced trajectory data and achieve good performance on generating structural road intersection map.

Discussion and Conclusions
Road intersections are critical parts of traffic road networks, which are characterized by complex structures and diverse traffic rules.Due to the limited road intersection information in current road map products, numerous studies attempt to extract road intersections from crowdsource trajectory data.However, crowdsourced trajectory data may be faced with spatio-temporal heterogeneity and uneven density distribution, leading to larger improvement spaces in the robustness to extract road intersections with different shapes and structures.This paper proposes a road intersection extracting method based on adaptive density homogenization and hierarchical trajectory segmentation.Firstly, raw trajectory data are processed to identify high-density trajectory cells and a trajectory density homogenization is performed on high-density trajectory cells in order to effectively reduces the time consumption of direct processing all raw trajectory data, and simultaneously averting the misidentification of non-intersection trajectory points and the omission of true road-intersection trajectory points caused by uneven density distribution.Moreover, in light of road hierarchies, a hierarchical trajectory segmentation strategy of category-level, object-level and turning-modes is then adopted to extract multi-level road intersection elements from coarse scale to fine scale.Finally, the iterative piecewise fitting method is utilized to generate the structural models of road intersections.A practical trajectory dataset is used for comprehensive analysis and the results show that the proposed method can efficiently deal with crowdsourced trajectory data with low sampling frequency and uneven density distribution, achieving high accuracy and recall rate of road intersection extraction.Especially, the proposed indicator of turning distance ratio can effectively reduce the misidentification of irregular turning tracks at parking lots, gas stations and residential areas.Nevertheless, there are also some aspects to be studied in the

Discussion and Conclusions
Road intersections are critical parts of traffic road networks, which are characterized by complex structures and diverse traffic rules.Due to the limited road intersection information in current road map products, numerous studies attempt to extract road intersections from crowdsource trajectory data.However, crowdsourced trajectory data may be faced with spatio-temporal heterogeneity and uneven density distribution, leading to larger improvement spaces in the robustness to extract road intersections with different shapes and structures.This paper proposes a road intersection extracting method based on adaptive density homogenization and hierarchical trajectory segmentation.Firstly, raw trajectory data are processed to identify high-density trajectory cells and a trajectory density homogenization is performed on high-density trajectory cells in order to effectively reduces the time consumption of direct processing all raw trajectory data, and simultaneously averting the misidentification of non-intersection trajectory points and the omission of true road-intersection trajectory points caused by uneven density distribution.Moreover, in light of road hierarchies, a hierarchical trajectory segmentation strategy of category-level, object-level and turning-modes is then adopted to extract multi-level road intersection elements from coarse scale to fine scale.Finally, the iterative piecewise fitting method is utilized to generate the structural models of road intersections.A practical trajectory dataset is used for comprehensive analysis and the results show that the proposed method can efficiently deal with crowdsourced trajectory data with low sampling frequency and uneven density distribution, achieving high accuracy and recall rate of road intersection extraction.Especially, the proposed indicator of turning distance ratio can effectively reduce the misidentification of irregular turning tracks at parking lots, gas stations and residential areas.Nevertheless, there are also some aspects to be studied in the future research.On the one hand, the calculation of neighborhood circle and the process of trajectory density homogenization should be improved to ensure enough trajectories of different traffic modes being reserved.On the other hand, the turning-modes segmentation method should introduce more trajectory features because eight directions spider code will be limited to extract large road intersections with more than 64 turning modes.Finally, tedious parameter settings should also be relieved, especially for decreasing the parameter number and the parameter sensitivity.

Figure 1 .
Figure 1.Spatio-temporal heterogeneities and uneven density distribution of crowdsourced trajectory data: (a) inconsistent trajectory shape and (b) varied trajectory densities.

Figure 1 .
Figure 1.Spatio-temporal heterogeneities and uneven density distribution of crowdsourced trajectory data: (a) inconsistent trajectory shape and (b) varied trajectory densities.

Figure 2 .
Figure 2. The workflow of the proposed method.

Figure 3 .
Figure 3. Calculating the neighborhood circles for high-density cell

Figure 3 .
Figure 3. Calculating the neighborhood circles for high-density cell.

Figure 4 .Figure 4 .
Figure 4. Calculating the turning distance ratio based on bend hierarchical tree Figure 4. Calculating the turning distance ratio based on bend hierarchical tree.

Figure 5 .
Figure 5. Three indicators of speed change, direction change and turning distance ratio for categorylevel segmentation: (a) speed charge, (b) direction change and (c) turning distance ratio.

Figure 5 .
Figure 5. Three indicators of speed change, direction change and turning distance ratio for categorylevel segmentation: (a) speed charge, (b) direction change and (c) turning distance ratio.

Figure 6 .
Figure 6.Eight directions spider code of heading directions of trajectory points

Figure 6 .
Figure 6.Eight directions spider code of heading directions of trajectory points.

Figure 7 .
Figure 7. Turning-modes segmentation based on heading code matrix of starting and ending points: (a) the heading code matrix of starting and ending points and (b) the trajectory segments of typical turning modes.

Figure 7 .
Figure 7. Turning-modes segmentation based on heading code matrix of starting and ending points: (a) the heading code matrix of starting and ending points and (b) the trajectory segments of typical turning modes.

Figure 7 .
Figure 7. Turning-modes segmentation based on heading code matrix of starting and ending points (a) the heading code matrix of starting and ending points and (b) the trajectory segments of typica turning modes.

Figure 8 .
Figure 8. Delineating road intersection structures by iterative piecewise fitting algorithm

Figure 8 .
Figure 8. Delineating road intersection structures by iterative piecewise fitting algorithm.

Figure 9 .
Figure 9.Comparison of road intersections detection without/with density homogenization.Figure 9. Comparison of road intersections detection without/with density homogenization.(a) without density homogenization, (b) with density homogenization.

Figure 9 .
Figure 9.Comparison of road intersections detection without/with density homogenization.Figure 9. Comparison of road intersections detection without/with density homogenization.(a) without density homogenization, (b) with density homogenization.

Figure 10 .
Figure 10.Results of identifying turning modes and extracting structural models of road intersection.

Figure 10 .
Figure 10.Results of identifying turning modes and extracting structural models of road intersection.(a) the structural models of road intersection.(b-e) the enlarged views of the generated geo-metric models and extracted turning modes of the two road intersections in the green and blue circles in (a).

Table 1 .
The main information in raw trajectory data.

Table 2 .
Experimental parameters setting

Table 3 .
Accuracy statistics of road intersections and turning modes identification.