Efﬁcient Algorithm for Constructing Order K Voronoi Diagrams in Road Networks

: The order k Voronoi diagram (OkVD) is an effective geometric construction to partition the geographical space into a set of Voronoi regions such that all locations within a Voronoi region share the same k nearest points of interest (POIs). Despite the broad applications of OkVD in various geographical analysis, few efﬁcient algorithms have been proposed to construct OkVD in real road networks. This study proposes a novel algorithm consisting of two stages. In the ﬁrst stage, a new one-to-all k shortest path ﬁnding procedure is proposed to efﬁciently determine the shortest paths to k nearest POIs for each node. In the second stage, a new recursive procedure is introduced to effectively divide boundary links within different Voronoi regions using the hierarchical tessellation property of the OkVD. To demonstrate the applicability of the proposed OkVD construction algorithm, a case study of place-based accessibility evaluation is carried out. Computational experiments are also conducted on ﬁve real road networks with different sizes, and results show that the proposed OkVD algorithm performed signiﬁcantly better than state-of-the-art algorithms.


Introduction
The Voronoi diagram (VD) is one of the fundamental geographical constructions for partitioning a geographical space through a competitive means [1][2][3][4]. Specifically, the VD partitions the geographic space into a set of Voronoi regions based on a given set of points of interest, also known as generators. Each Voronoi region contains all the locations that share the same nearest generator, allowing for space tessellation and providing an adjacency structure for topology generation [5]. In addition, a Voronoi region can be used as container for geographical entities, making it an effective indexing structure for organizing spatial data. Due to these unique characteristics, the VD has been widely used in various geographical applications, such as mobile phone data analysis [6], facility service area delimitations [7], spatial data mining [8,9], and etc.
The order k Voronoi diagram (OkVD) is a generalization of the VD that has also received significant attention in the literature [1,10,11]. Given a set of generators, the OkVD is meant to partition the geographical space into a set of order k Voronoi regions (OkVR), where all locations within an OkVR share the same k nearest generators. These generators are determined by their order of proximity to the location, i.e., the first nearest generator, the second nearest generator, and so on until the kth nearest generator. This OkVD provides an effective tool to determine k nearest generators for any location, which is a critical step for numerous spatial analysis tasks, such as spatial interpolations [12], accessibility evaluations [13], and k nearest neighbor queries [14][15][16]. Given the broad applications of the OkVD in the geographical analysis, it is therefore necessary to develop efficient algorithms for constructing OkVDs in real-world geographical spaces.
In the literature, efficient algorithms for constructing Voronoi diagrams (VDs) have been the subject of much research. Early studies focused on the construction of VDs in homogeneous and isotropic planar spaces, where the distance between any two locations is measured by Euclidean distance [3,17]. However, recognizing that movements in urban areas is generally constrained by road networks, Erwig [18] first proposed a network Voronoi diagram (NVD) model for partitioning the urban geographical space. The NVD model quantifies distances between any two locations in the road network explicitly by using shortest path algorithms to compute the network distance [19]. Building on this work, Okabe et al. [1] developed an efficient algorithm for constructing NVDs by modifying conventional shortest path algorithms. Ohsawa [20] applied the concept of Order k Voronoi region for solving the problem of k nearest neighbor queries. As an alternative approach, Ai et al. [2] proposed an NVD construction algorithm by discretizing each road link into a set of linear units of equal length and using the stream flowing approach to construct NVDs. Chen et al. [21] proposed several new NVD models for partitioning the urban geographic space constrained by multi-mode public transport networks and developed effective algorithms for constructing NVDs in public transport networks.
To the best of our knowledge, there are few algorithms in the literature have been developed to construct order k NVD (OkNVD) in real road networks, with the notable exception of Okabe et al. [1]. For convenience, this algorithm is hereafter referred to as Okabe's algorithm. This algorithm consists of two stages. In the first stage, the k nearest generators for each node are determined by performing the one-to-all shortest path algorithm (e.g., Dijkstra's algorithm) |F| times to construct completed one-to-all shortest path trees rooted at all generators, where |F| is the number of generators. The k nearest generators for each node can be determined by sorting |F| shortest paths from all generators to the node. In the second stage, the OkNVD is constructed by determining all nodes and links belonging to each order k network Voronoi region (OkNVR). The idea of a lower-bound envelope is introduced to divide boundary links that belongs to at least two OkNVRs. Although Okabe's algorithm is straightforward and easy to implement, the first stage of Okabe's algorithm tends to have computational overheads, because it calculates too many unnecessary shortest paths, given that |F| is generally much larger than the k value. Furthermore, no detailed procedure was provided for implementing the idea of a lower-bound envelope method in the second stage.
Along the line of previous work [1], this study proposes an effective and efficient algorithm for constructing OkNVD in road networks. The proposed algorithm also consists of two stages, which extends the previous algorithm in the following aspects.

1.
A new one-to-all k shortest path procedure is proposed in the first stage to efficiently determine the k nearest generators at each node. As an extension of Okabe's algorithm, the proposed procedure simultaneously constructs |F| shortest path trees in a single path searching process, sharing the same priority queue to coordinate the construction of shortest path trees from different generators. Unlike Dijkstra's algorithm [22], the proposed procedure closes a node after it has been selected from the priority queue k times from different generators, allowing it to exactly determine the shortest paths from the k nearest generators. This results in partial shortest path trees being constructed, improving Okabe's algorithm. Furthermore, the proposed procedure can directly determine the order of k nearest generators according to their order of selection from the priority queue without the need for shortest path sorting at each node.

2.
A new recursive procedure is introduced in the second stage to effectively divide boundary links into different OkNVRs. This study rigidly proves the hierarchical tessellation property of OkNVD. Specifically, each order j NVR is covered by a set of order j + 1 NVRs, and these order j + 1 NVRs are mutually exclusive except at the boundary points (j = 0, 1, . . . , k − 1). Using this hierarchical tessellation property, the introduced procedure, therefore, can recursively divide every boundary link into the order first, second, . . . , until the kth NVRs.

3.
To demonstrate the applicability of the proposed OkNVD construction algorithm, a comprehensive case study of evaluating place-based accessibility is carried out in Shenzhen, a mega-city in China. Computational experiments are also conducted using five real road networks with various settings of k and |F| values. Results show that the proposed OkVD algorithm performed significantly better than state-of-theart algorithms.
The remainder of this article is organized as follows. Section 2 gives the model formulation of OkNVD and proves its hierarchical tessellation property. Section 3 describes the proposed OkNVD construction algorithm. Section 4 reports the numerical example of OkNVD application in place-based accessibility, and computational experiments using five real road networks. Section 5 presents the conclusions together with recommendations for future research in the area.

Model Formulation of Order k Network Voronoi Diagrams
The list of notations used in this article is summarized in Abbrevations. A road network can be represented by a directed graph, G(N, A), consisting of a set of nodes N, and a set of links A. Each link a ij A has a tail node n i N and a head node n j N, and a positive link cost t ij (length or travel time). Each node n j has a set of successor nodes SUCC n j and a set of predecessor nodes PRED n j . Any location q in the road network can be represented as a ij , θ q using the linear reference technique [23], where θ q ∈ [0, 1] is the relative position on the link a ij . The values of θ q = 0 and θ q = 1 refer to the tail and head nodes of the link respectively. With location q, the link a ij can be divided into two partial links a iq and a qj , and their link costs are assumed to be proportional to relative positions as: Let p oq be the shortest path from origin o to destination q. The path cost, denoted by t oq , can be calculated by the summation of corresponding link costs along the path: where δ oq ij is the binary variable of link-path incidence relationship, δ oq ij = 1 indicates link a ij on the path p oq and otherwise δ oq ij = 0. Let F = f 1 , . . . , f i , . . . , f |F| be a set of generators located in the road network, where |F| is the number of generators in F. Generators are typically a certain type of service facility. Given any location q in the network, its first outward nearest generator, denoted by f q,1 out ∈ F, satisfies: t f q Then, the set of k outward nearest generators can be determined and denoted by F q,K out = f q,1 out , . . . , f q,K out . Therefore, the outward OkNVD (order k network Voronoi diagram), denoted by VD K out , is to partition the road network into a set of sub-networks, called outward OkNVRs (order k network Voronoi regions), denoted by . . . , vr K,i Out , . . . , such that all locations within any vr K,i Out share the same set of k nearest generators, F q,K out . The outward OkNVR, vr K,i Out , can be expressed as where ∧ represents that all conditions are satisfied simultaneously. Because the road network is directed, the inward OkNVD can also be defined by using the shortest path in a reverse direction from location q to the corresponding generator. Let be the set of k inward nearest generators. The jth inward nearest generator satisfies following condition: where t q f j is the cost of the shortest path from location q to the jth inward nearest generator f q,j in . Then, the inward OkNVD, denoted by VD K in , consists of a set of inward OkNVRs, VD K in = . . . , vr K,i in , . . . . Each inward OkNVR, vr K,i in , delimits all locations sharing the same set of k inward nearest generators, F q,K in , and it can be mathematically expressed as The outward and inward OkNVDs are applicable for different types of generators and different applications. The outward OkNVDs are commonly for emergency facilities (e.g., fire stations) and applications whose outward travel cost t f q i is critical; while the inward OkNVDs are usually for service facilities (e.g., supermarkets and shopping malls) and applications whose inward travel cost t q f i is critical. Figures 1 and 2 illustrate the outward and inward OkNVDs in a well-known small road network, namely the Sioux Falls network. Figure 1a,b respectively give the outward and inward OkNVDs when k = 1, in which OkNVDs degrade to the classical outward and inward NVDs. In both two figures, generators are given in star shapes, and their Voronoi regions are shown in the same color. As can be seen, the outward NVDs are not identical to the inward counterparts, since the road network is directed. It can also be seen that both outward and inward NVDs form a tessellation on the road network with two observations: (1) All NVDs cover the whole network (i.e., collectively exhaustive); and (2) Voronoi regions are mutually exclusive except at the boundary points. This is because the Sioux Falls network (as most real road networks) is fully connected, i.e., there exists at least one path between any two locations in the road network.  and . As can be seen, the outward (or inward) order second NVDs not only form a tessellation on the whole network but also on each outward (or inward) order first NVR. Thereby, the order-second NVD together with the order first NVD form a hierarchical tessellation on the road network. This hierarchical tessellation can also be represented by the tree structure shown in Figure 3a,b, in which the Voronoi regions were kept the same color as Figure 2a,b. The original road network is represented as a single root node at level 0 of the tree. The order first NVD is to partition the original network into a set of order first NVRs, which are represented as child nodes of nodes at level 1. Each order first NVR is further divided into a set of order second NVRs representing as its child nodes, i.e., nodes at level 2.   and . As can be seen, the outward (or inward) order second NVDs not only form a tessellation on the whole network but also on each outward (or inward) order first NVR. Thereby, the order-second NVD together with the order first NVD form a hierarchical tessellation on the road network. This hierarchical tessellation can also be represented by the tree structure shown in Figure 3a,b, in which the Voronoi regions were kept the same color as Figure 2a,b. The original road network is represented as a single root node at level 0 of the tree. The order first NVD is to partition the original network into a set of order first NVRs, which are represented as child nodes of nodes at level 1. Each order first NVR is further divided into a set of order second NVRs representing as its child nodes, i.e., nodes at level 2. Figure 2a,b respectively shows the outward and inward OkNVDs when k = 2 in the same Sioux Falls network. In these two OkNVDs, six OkNVRs are given in different colors. Each OkNVR delimits all locations sharing the same k nearest generators. For example, all locations within in vr 2,1 out have the same first and second nearest generators, f 1 and f 2 . As can be seen, the outward (or inward) order second NVDs not only form a tessellation on the whole network but also on each outward (or inward) order first NVR. Thereby, the order-second NVD together with the order first NVD form a hierarchical tessellation on the road network. This hierarchical tessellation can also be represented by the tree structure shown in Figure 3a,b, in which the Voronoi regions were kept the same color as Figure 2a,b. The original road network is represented as a single root node at level 0 of the tree. The order first NVD is to partition the original network into a set of order first NVRs, which are represented as child nodes of nodes at level 1. Each order first NVR is further divided into a set of order second NVRs representing as its child nodes, i.e., nodes at level 2.   Since both outward and inward OkNVDs have an identical tessellation property, we only present the outward case for simplicity. Given a OkNVR, , , it is said to be a child region of order − 1 NVR, , , if , = , ∪ { , } holds. All child NVRs of , form the child region set, denoted by CRS , = … , , , … ⊂ . Conversely, this , is said to be the parent region of any , ∈ CRS , . The tessellation property of OkNVD can be proved rigidly as below.

Proposition 1. If the road network is fully connected, the outward order first NVD forms a tessellation on the road network.
Proof. See Proposition A1 in the Appendix B. □

Proposition 2.
If the road network is fully connected, , forms a tessellation on , .
Proof. See Proposition A2 in the Appendix B. □ According to Propositions 1 and 2, OkNVRs are stacked under their parent NVRs (i.e., order − 1 NVRs, …, until the order first NVRs), forming a level hierarchical tessellation on the road network. This hierarchical tessellation property of OkNVDs can be useful for many applications to organize spatial data. Such a property is fully utilized in the next section to construct OkNVDs.

Proposed Algorithm for Constructing OkNVDs
This section presents the proposed algorithms for constructing outward and inward OkNVDs in the road network. We first describe the algorithm for constructing outward OkNVD, called ConstructOutwardOkNVD. As shown in Algorithm 1, this ConstructOut-wardOkNVD algorithm consists of two stages. The first stage is to efficiently determine the shortest paths to , at each node . To this end, a new one-to-all k shortest path finding procedure, called FindKNearestPOIs, is proposed. The second stage is to construct the outward OkNVD. Each node only belongs to a OkNVR and can be easily determined using the calculated , . However, a link may belong to multiple OkNVRs. In the second stage, a recursive procedure, called AddLinkToOkNVRs, is introduced to effectively divide the link into corresponding OkNVRs using the above hierarchical tessellation property. The detailed steps of FindKNearestPOIs and AddLinkToOkNVRs procedures are described in the following sections. Since both outward and inward OkNVDs have an identical tessellation property, we only present the outward case for simplicity. Given a OkNVR, vr K,i out , it is said to be a child region of order holds. All child NVRs of vr K−1,i out form the child region set, denoted by CRS K−1,i out = . . . , vr K,i out , . . . ⊂ VR K out . Conversely, this vr K−1,i out is said to be the parent region of any vr K,i out ∈ CRS K−1,i out . The tessellation property of OkNVD can be proved rigidly as below.

Proposition 1.
If the road network is fully connected, the outward order first NVD forms a tessellation on the road network.
Proof. See Proposition A1 in the Appendix A.

Proposition 2.
If the road network is fully connected, CRS K−1,i out forms a tessellation on vr K−1,i out .
Proof. See Proposition A2 in the Appendix A.
According to Propositions 1 and 2, OkNVRs are stacked under their parent NVRs (i.e., order k − 1 NVRs, . . . , until the order first NVRs), forming a k level hierarchical tessellation on the road network. This hierarchical tessellation property of OkNVDs can be useful for many applications to organize spatial data. Such a property is fully utilized in the next section to construct OkNVDs.

Proposed Algorithm for Constructing OkNVDs
This section presents the proposed algorithms for constructing outward and inward OkNVDs in the road network. We first describe the algorithm for constructing outward OkNVD, called ConstructOutwardOkNVD. As shown in Algorithm 1, this ConstructOutwar-dOkNVD algorithm consists of two stages. The first stage is to efficiently determine the shortest paths to F j,K out at each node n j . To this end, a new one-to-all k shortest path finding procedure, called FindKNearestPOIs, is proposed. The second stage is to construct the outward OkNVD. Each node n j only belongs to a OkNVR and can be easily determined using the calculated F j,K out . However, a link a ij may belong to multiple OkNVRs. In the second stage, a recursive procedure, called AddLinkToOkNVRs, is introduced to effectively divide the link into corresponding OkNVRs using the above hierarchical tessellation property. The detailed steps of FindKNearestPOIs and AddLinkToOkNVRs procedures are described in the following sections.

The FindKNearestPOIs Procedure
To calculate shortest paths from F j,K out at each node n j , a straightforward approach [1] consists of two steps. The first step is to perform the classical one-to-all Dijkstra's algorithm by |F| times to construct the completed one-to-all shortest path trees rooted at all generators. Accordingly, |F| shortest paths from all generators at each node n j are calculated. The second step is to sort the calculated |F| shortest paths at each node to determine F j,K out . Using the F-Heap data structure [24], Dijkstra's algorithm runs in the worst-case complexity of O(|N|Log(|N|) + |A|), where |N| is the number of nodes and |A| is the number of links. Therefore, the step requires O(|F||N|Log(|N|) + |F||A|). Since the second step requires O(|N||F| 2 ), this approach runs in the worst-case complexity of O(|F||N|Log(|N|) + |F||A| + |N||F| 2 ). Such an approach is easy for implementation. However, its computational performance significantly degrades with |F|, i.e., the number of generators. It tends to have computational overheads by calculating too many unnecessary shortest paths, given that |F| is generally much larger than the k value in practice.
In this study, a new one-to-all k shortest path procedure, i.e., FindKNearestPOIs, is proposed to efficiently determine only k shortest paths from F j,K out rather than |F| shortest paths from all generators. The proposed procedure modifies the Dijkstra's algorithm in three aspects. Firstly, the proposed procedure simultaneously constructs |F| shortest path trees in a single path searching process. The Dijkstra's algorithm maintains only one label at each node. The proposed procedure, however, maintains a label list with |F| size at each node to record the shortest paths from different generators. Secondly, the proposed procedure makes use of only one priority queue for constructing all |F| shortest path trees. By sharing a single priority queue, shortest path trees from all generators can be coordinated during the path-searching process. Thirdly, the proposed procedure closes a node when it was selected from the priority queue by k times instead of a single selection used in Dijkstra's algorithm. After the node was closed, k shortest paths from F j,K out were exactly determined at the node according to the monotonic increasing property of the priority queue. In addition, the ascending order of shortest paths from different generators can be directly obtained according to their selection order from the priority queue. Therefore, the proposed procedure can efficiently determine k shortest paths from F j,K out without the requirement of constructing |F| completed shortest path trees from all generators and sorting the shortest paths at each node. Algorithm 2 gives the detailed steps of the proposed FindKNearestPOIs procedure. Unlike Dijkstra's algorithm, each node n j in the proposed procedure maintains three objects, including a label set LS j , a selected generator set SG j and an open status os j . The label set LS j employs an array list structure with the size of |F|. The ith location of the label set, denoted by LS j [i], records the shortest path from the ith generator f i . This way efficiently retrieves the shortest path from a specific generator. Note that it is not necessary to generate |F| labels at each node; and a null object is used if the shortest path from a certain generator has not been generated. The selected generator set, SG j , employs a linked list structure for recording the order of generators whose shortest paths to node n j were selected from the priority queue. Finally, the open status os j is a binary indicator, and os j = true indicates that n j is open.
In the initialization step, LS j and SG j of all nodes are set to be empty, and os j of all nodes are set to be true. A path (or called label) from each generator f i to itself, denoted by p f f , is created and set its path cost t f f as zero. An additional property, denoted by g f f , is introduced at the label to store the generator order i. This created path is added into the priority queue, denoted by SE.
At each iteration, a path p f j at the top of the priority queue SE is selected. Its generator order, g f j , is then added into the selected generator set SG j . Once the number of generators in SG j reached the k value, the node n j is closed by setting the os j attribute to be false. In this case, k shortest paths from F j,K out have been determined at node n j , and thus all paths of LS j still in the priority queue are eliminated without considerations. Subsequently, the selected path p f j is extended to its successor nodes, which are not closed. Through an open successor node n w , a temporal path p f w is created by p f j ⊕ a jw , where ⊕ is a path concatenation operator. The newly created path p f w compares to the existing path p f w as LS w g f w from the same generator. If the existing path p f w = ∅ or t f w < t f w holds, the newly created path is kept and inserted into the priority queue; and it is discarded otherwise. The procedure continues the path selection and extension process until SE is empty. When the procedure terminates, it can determine the k nearest generators of each node n j (i.e., F j,K out ) at the selected generator set SG j and k shortest paths from F j,K out of each node n j at the label set LS j .
We can prove the optimality of the proposed FindKNearestPOIs procedure as below.
Proposition 3. When the proposed FindKNearestPOIs procedure terminates, it can determine shortest paths from K nearest generators F j,K out for all nodes.
Proof. See Proposition A3 in the Appendix A.
The worst-case complexity of the proposed FindKNearestPOIs procedure is analyzed. Using the F-Heap data structure [24], the selection of the minimum cost path from the priority queue requires O(Log(|F||N|)), and both remove and add path into the priority queue requires O(1), where |F||N| is the maximum paths in the priority queue. Since k|N| is the maximum number of path selection from the priority queue, the proposed procedure thus runs in the worst-case complexity of O(k|N|Log(|F||N|) + k|A|). This complexity is significantly better than the previous Okabe's algorithm, i.e., O(|F||N|Log(|N|) + |F||A| + |N||F| 2 ).

The AddLinkToOkNVRs Procedure
In Okabe's algorithm [1], a complicated idea of a lower-bound envelope was introduced to divide boundary links which belong to at least two distinctive OkNVRs. However, no detailed procedure was given to implement such an idea.
In this study, a new recursive procedure, AddLinkToOkNVRs, is introduced to effectively divide boundary links into different OkNVRs using the hierarchical tessellation property, i.e., Propositions 1 and 2. The detailed steps of the introduced procedure is given in Algorithm 3. Figure 4 illustrates this procedure using a simple example of two directed links a jw and a wj . Two end nodes, n j and n w , are given in the figure with the selected generator sets, SG j = {1, 2, 3} and SG w = {4, 5, 1}. Because SG j = SG w holds, these two links are not in the same OkNVR, so called boundary links (their end nodes are shown in different colors). This AddLinkToOkNVRs procedure divides boundary links into a set of partial links through the hierarchical tessellation from the first level until the kth level. When dividing links at the first level, we focus on only SG j [1] (i.e., Generator 1) and SG w [1] (i.e., Generator 4) and their shortest paths, p f j 19: For each successor node n w ∈ SUCC n j 20: If os w = true Then 21: Create path p f w := p f j ⊕ a jw , and Set t f w := t f j + t jw and g f w := g f j .

22:
Retrieve existing path p f w := LS w g f w . 23 :  The procedure consists of three steps.  The procedure consists of three steps. The first step (i.e., Step 2.1) is to divide links a jw and a wj into four partial links. Using the equal travel times between left-bound route  On the same break point, we create two new nodes, namely a left-bound node n u and right-bound node n v . Then, we divide two links a jw and a wj into four newly partial links, including two left-bound partial links, a ju and a uj , and two right-bound partial links, a vw and a wv . For example, partial link a ju = a jw , 0, θ u represents part of link a jw from 0 to θ u using the linear reference.
The second step (i.e., Step 2.2) is to determine the selected generator sets (SG u and SG v ) and label sets (LS u and LS v ) for newly created nodes n u and n v . It is easy to determine SG j [1] (e.g., Generator 1) and SG w [1] (e.g., Generator 4) as the first and second generators of the left-bound node n u , e.g., SG u = {1, 4}. A reverse order is used for the right-bound node n v , i.e., SG v = {4, 1}. Since LS u and LS v have identical path sets, we use LS to represent them, i.e., LS [1]  maintained at n j and n w respectively. All constructed candidate paths are inserted into the priority queue, denoted by CSE. Then, the top k − 2 paths from CSE are identified and stored at LS, and their generators are inserted into SG u and SG v , e.g., SG u = {1, 4, 2} and SG v = {4, 1, 2}. After this step was performed, we divided two input links into four partial links within two NVRs, e.g., vr 1,1 out and vr 1,2 out , at the first level.
The last step (i.e., Step 2.3) is to recursively call the AddLinkToOkNVRs procedure for the link division at the next level by using the newly created left-bound links (a ju and a uj ) and right-bound links (a vw and a wv ) as two input links respectively.
For the link division at the ith level, this procedure focuses on only Scenario 2 of SG j [i] = SG w [i]. This is because we have Scenario 1 of SG j [l] = SG w [l] for ∀l = 1, . . . , i, in which SG j [l] and LS j SG j [l] are kept for the setting (SG u and SG v ) and label sets (LS u and LS u ) at the ith level. The link division of the ith level has the similar three steps as that of the first level. This procedure recursively performed until Scenario 3 of SG u [l] = SG v [l] for ∀l = 1, . . . , k holds. After the procedure terminated, input links were divided into a set of partial links belonging to different OkNVRs, e.g., vr 3,1 out , . . . , vr 3,6 out in the illustrative example. i + t jw θ u and g fu i := g fj i . 09: Create right-bound path p fv i := p fw i ⊕ a wv and Set t fv i := t fw i + t wv θ v and g fv i := g fw i . 10: Create left-bound node n u and partial links a ju := a jw , 0, θ u and a uj := a wj , θ v , 1 . 11: Create right-bound node n v and partial links a vw := a wj , θ u , 1 and a wv := a jw , 0, θ v .
Step 2.2. Set the selected generator set and the label set for two newly created nodes: 12 Create candidate path p fu l := p fj l ⊕ a ju and Set t fu l := t fj l + t jw θ u and g fu l := g fj l . 20: Create candidate path p fv l := p fw l ⊕ a wv and Set t fv l := t fw l + t wv θ v and g fv l := g fw l . We can prove the optimality of the introduced AddLinkToOkNVRs procedure below.

Proposition 4.
When the introduced AddLinkToOkNVRs procedure terminates, it can correctly divide two input links into a set of OkNVRs.

Proof. See Proposition A4 in Appendix A.
The worst-case complexity of the introduced AddLinkToOkNVRs procedure is also analyzed. In the procedure, Step 2.1 runs in O(1); and Step 2.2 runs in O(kLog(k)) using the F-heap data structure. Because the procedure recursively divides input links from the first level until the k level, it runs in the complexity of O k 2 Log(k) .

Complexity Analysis of the Proposed OkNVD Construction Algorithms
With Propositions 3 and 4, we can easily prove the optimality of the proposed Con-structOutwardOkNVD algorithm as follows.

Proposition 5.
When the ConstructOutwardOkNVD algorithm terminates, it can correctly construct outward OkNVDs in the road network.

Proof. It can be easily followed by Propositions 3 and 4.
Based on the worst-case complexity of AddLinkToOkNVRs procedure (see Section 3.2), we can determine that the second stage of the ConstructOutwardOkNVD algorithm runs in O |A|k 2 Log(k) . Since its first stage runs in O(k|N|Log(|F||N|) + k|A|) (see Section 3.1), the ConstructOutwardOkNVD algorithm therefore runs in O k|N|Log(|F||N|) + |A|k 2 Log(k) .
The proposed ConstructOutwardOkNVD algorithm can be modified easily to construct inward OkNVDs, namely the ConstructInwardOkNVD algorithm. Both FindKNearestPOIs and AddLinkToOkNVRs procedures should be modified by using the backward search (i.e., path extension to predecessor nodes) instead of the forward search (i.e., path extension to successor nodes). For example, the modified FindKNearestPOIs procedure should create a new path p w f := a jw ⊕ p f j from each predecessor node n w ∈ PRED n j (see Lines 19 and 22). The ConstructInwardOkNVD algorithm has the same worst-case complexity as the ConstructOutwardOkNVD algorithm.

Numerical Experiments
This section reports the numerical experiments using real road networks. The proposed algorithm was implemented using the C# programming language. Section 4.1 reports their applications in place-based accessibility to supermarkets in Shenzhen, China; and Section 4.2 investigates their computational performance in several real road networks with different sizes.

Place-Based Accessibility to Supermarkets
In this study, Shenzhen, one of the most developed mega-cities in China, was selected as the study area. As shown in Figure 5a, Shenzhen is located in southern China, bordering Hong Kong to the south. It is the first Chinese special economic zone and experienced rapid socioeconomic development during the past decades. By the end of 2013, Shenzhen covered an area of approximately 1996 km 2 and had a total population of approximately 10.5 million, of which more than 70% were immigrants. Shenzhen has ten administrative districts with diverse land-use patterns. Among these districts, Nanshan, Futian and Luohu are core urban regions with dense residential and commercial areas; Bao'an, Longhua, Longgang, and Yantian are suburban regions with several new towns and industrial areas; and other three districts (Guangming, Pingshan, and Dapeng) are rural regions with large-scale agriculture and hilly lands as well as industrial areas. Three datasets were collected for the case study. The first dataset was 39 larg markets in Shenzhen shown by star symbols in Figure 5a. The second dataset was network of Shenzhen city. As shown in Figure 5b, it consisted of 32,065 nodes an directed links. The third dataset was trajectories of 17,406 taxis collected on a typic day, i.e., 23 September 2014. These trajectories were matched onto the road netwo a map matching algorithm [17], and then used to estimate and link travel time info during the evening peak hour (18:00-19:00) by using the method of Shi et al. [25] Using supermarkets as generators, Figure 6a shows the constructed inward when = 1. The generators are shown in star shape. The proposed algorithm con 39 inward OkNVRs, i.e., , , for all generators. Each OkNVR, shown by a di color, delimits a spatial region in which all locations share the same nearest gener can be seen, the sizes of OkNVRs are spatially heterogeneous in Shenzhen city w ferent supermarket densities, i.e., smaller at core urban and suburban regions (e 1-5, 9-17 and 18-22) while larger at rural regions (e.g., Nos. 38 and 39). It can also that all OkNVRs formed a tessellation on the Shenzhen network. They covered th network and they were mutually exclusive except at the boundary points. Figure 6b shows OkNVD for large supermarkets in Shenzhen when = 2. A seen, the proposed algorithm constructed 149 inward OkNVRs, i.e., , , in th area. In this case, each constructed OkNVR delimits a spatial region in which all l share the same set of the first and second nearest generators, , . Compared to of = 1, the sizes of , are smaller, especially for those areas with dense su Three datasets were collected for the case study. The first dataset was 39 large supermarkets in Shenzhen shown by star symbols in Figure 5a. The second dataset was the road network of Shenzhen city. As shown in Figure 5b, it consisted of 32,065 nodes and 81,618 directed links. The third dataset was trajectories of 17,406 taxis collected on a typical workday, i.e., 23 September 2014. These trajectories were matched onto the road network using a map matching algorithm [17], and then used to estimate and link travel time information during the evening peak hour (18:00-19:00) by using the method of Shi et al. [25].
Using supermarkets as generators, Figure 6a shows the constructed inward OkNVD when k = 1. The generators are shown in star shape. The proposed algorithm constructed 39 inward OkNVRs, i.e., vr 1,i in , for all generators. Each OkNVR, shown by a distinctive color, delimits a spatial region in which all locations share the same nearest generator. As can be seen, the sizes of OkNVRs are spatially heterogeneous in Shenzhen city with different supermarket densities, i.e., smaller at core urban and suburban regions (e.g., Nos. 1-5, 9-17 and 18-22) while larger at rural regions (e.g., Nos. 38 and 39). It can also be seen that all OkNVRs formed a tessellation on the Shenzhen network. They covered the whole network and they were mutually exclusive except at the boundary points.    Figure 6b shows OkNVD for large supermarkets in Shenzhen when k = 2. As can be seen, the proposed algorithm constructed 149 inward OkNVRs, i.e., vr 2,i in , in the study area. In this case, each constructed OkNVR delimits a spatial region in which all locations share the same set of the first and second nearest generators, F q,2 in . Compared to the case of k = 1, the sizes of vr 2,i in are smaller, especially for those areas with dense supermarkets. This result is obvious because all inward OkNVRs in the case of k = 2 form a tessellation in the case of k = 1. Each vr 1,i in is covered by one or more vr 2,i in , which are mutually exclusive except at the boundary points. Figure 6c shows the corresponding OkNVD when k = 3. In this case, the whole study area was partitioned into 320 inward OkNVRs, i.e., vr 3,i in . The constructed OkNVRs when k = 3 further form a tessellation when k = 2. Table 1 summarizes the characteristics of OkNVRs under different k values. As can be seen, the average size of OkNVRs (in terms of cumulative link lengths) was reduced with the increase of k values. The total size of all OkNVRs under different k values keep constant, i.e., 5530 km. This is because all OkNVRs at three levels (k = 1 to 3) are stacked, forming a three-level hierarchical tessellation on the whole road network. , for any location q in the road network were calculated. These calculated travel times were employed to quantify the place-based accessibility to large supermarkets in Shenzhen and shown in Figure 7.  Figure 7a shows the accessibility pattern using the least travel time to the nearest supermarket, i.e., calculated in the order first NVD. This is a classical indicator for measuring place-based accessibility in the literature [26]. As can be observed, the level of accessibility to large supermarkets is high for areas surrounded by at least one supermarket but low for those areas far from supermarkets. This indicator assumes that residents always go to the nearest supermarket for food purchases. However, several empirical  Figure 7a shows the accessibility pattern using the least travel time to the nearest supermarket, i.e., t q f 1 calculated in the order first NVD. This is a classical indicator for measuring place-based accessibility in the literature [26]. As can be observed, the level of accessibility to large supermarkets is high for areas surrounded by at least one supermarket but low for those areas far from supermarkets. This indicator assumes that residents always go to the nearest supermarket for food purchases. However, several empirical studies have found that the most frequently used supermarkets by most residents are not the nearest supermarket to home [27]. As an extension to this classical indicator, Chen et al. [28] recently suggested using the average travel time (denoted by t q f k ) of the k nearest supermarkets to better quantify accessibility to food retailers. Specifically, the accessibility level at location q is calculated as follows: A smaller t q f k value implies the easer of a location for accessing supermarkets. This indicator can be directly obtained by using the results of OkNVD. Figure 7b shows the accessibility result using the indicator of t q f k when k = 2. Compared to the classical indicator of t q f 1 , this indicator can well reveal the rural-urban disparity of accessibility to supermarkets caused by heterogeneous supermarket densities. The accessibility level was relatively high in the core urban regions (Nanshan, Futian, and Luohu) and relatively low in rural regions (Guangming, Pingshan, and Dapeng). For example, such indicators can well identify the low accessibility level for areas nearby Supermarket No. 39 in Dapeng district, since only one large supermarket is available in those areas. Similarly, a more distinct rural-urban disparity of accessibility can be observed when k = 3 is used, see Figure 7c. As the value of k increases, the accessibility level of a location with low supermarket density would decrease rapidly, as they lack alternative supermarket supplies. However, a large k value could cover infeasible supermarkets with too large travel times. The selection of a suitable k parameter is critical for this accessibility measure. As suggested by Chen et al. [28], k = 3 is appropriate for supermarket accessibility studies.

Computational Experiments
This section reports the computational performance of the proposed algorithm. Previous Okabe's algorithm [1] was also implemented for comparisons. Since Okabe's algorithm did not provide the detailed procedure for the second stage, we employed that of our proposed algorithm. Both algorithms were coded in the same C# programming language and used the same F-heap data structure [24] as the priority queue. Since constructing inward and outward OkNVDs requires the same computational performance, we only report that of inward OkNVD constructions for simplicity. All experiments were conducted on a desktop with a four-core Intel 3.3 GHz CPU (only a single processor was used) and 8GB RAM running on Windows 10 operation system. Table 2 reports the computational performance of both OkNVD construction algorithms in Shenzhen network under a different number of generators (i.e., |F| values). Network nodes were randomly selected as new generators. The k value was set as 3. As can be seen, the computational performance of Stage 1 dominates that of both algorithms. Therefore, we only examined the computational performance of Stage 1 hereafter. In addition, computational cost was required to sort |F| shortest paths at each node.
As shown, the proposed algorithm has a much lower rate of degradation with the increase of |F|. When |F| increases from 100 to 1000, its computational time raised by only 91.7% from 2.40 s to 4.60 s. This is because the proposed algorithm calculates only k shortest paths at each node by constructing |F| partial shortest path trees. Therefore, the proposed algorithm performed significantly faster than Okabe's algorithm by calculating much less shortest paths, especially for the scenarios of large |F|. As shown, the proposed algorithm ran 10.45 (i.e., 27.49/2.40 − 1) times faster than Okabe's algorithm when |F| = 100. This computational advantage further increased to 65.34 times when |F| = 1000. Such computational advantage is important for real applications commonly with large |F|, for example |F| = 6899 in Chen et al. [6]. Figure 8 shows the computational performance of both OkNVD construction algorithms in Shenzhen network under different k values. The number of generators was set as 500. As can be seen, the computational performance of Okabe's algorithm is stable under different k values. It is obvious, since Okabe's algorithm calculates |F| shortest paths at each node regardless of different k settings. It can be seen that the computational performance of the proposed algorithm degraded linearly with the increase of k value. This is because the proposed algorithm calculates k shortest paths at each node. The larger the k value was, the more shortest paths were calculated. In practice, the k value is usually low, e.g., k ≤ 5 [28]. It should be noted that the proposed algorithm ran still significantly better than Okabe's algorithm various k values under various k values. Table 3 reports the computational performance of both testing algorithms in five real road networks with different sizes. The number of generators was selected as 500 and the k value was set as 3. Generators were randomly selected from the network nodes. As shown, the computational times of both testing algorithms increased with the network size. This result is expected, since the computational performance of shortest path algorithms themselves degraded with the network size [19]. It can be seen that the proposed algorithm performed consistently better than Okabe's algorithm in all road networks with different sizes. For example, in a small network of Xiamen, the proposed algorithm consumed 0.19 s to construct inward OkNVD, which was 36.73 times (i.e., 7.17/0.19 − 1) faster than Okabe's algorithm. This computational improvement remained relatively stable in all networks ranging from 31 to 37 times. formance of the proposed algorithm degraded linearly with the increase of value is because the proposed algorithm calculates shortest paths at each node. The larg k value was, the more shortest paths were calculated. In practice, the value is u low, e.g., ≤ 5 [28]. It should be noted that the proposed algorithm ran still signifi better than Okabe's algorithm various values under various values.  Table 3 reports the computational performance of both testing algorithms in fiv road networks with different sizes. The number of generators was selected as 500 a k value was set as 3. Generators were randomly selected from the network nod shown, the computational times of both testing algorithms increased with the ne size. This result is expected, since the computational performance of shortest path rithms themselves degraded with the network size [19]. It can be seen that the pro algorithm performed consistently better than Okabe's algorithm in all road network

Conclusions
The order k network Voronoi diagram (OkNVD) is an effective geometric construction to partition the network space into a set of order k network Voronoi regions (OkNVRs) such that all locations within an OkNVR share the same k nearest generators. Despite the broad applications of OkNVD, few effective and efficient algorithms had been developed to construct OkNVD in real road networks. This study proposed a novel OkNVD construction algorithm consisting of two stages. In the first stage, a new one-to-all k shortest path procedure was proposed to efficiently determine the k nearest generators for each network node in a single shortest path search process. This proposed procedure improves the worst-case complexity from previous O(|F||N|Log(|N|) + |F||A| + |N||F| 2 ) [1] to O(k|N|Log(|F||N|) + k|A|). In the second stage, a new recursive procedure was introduced to effectively divide boundary links into different OkNVRs using the hierarchical tessellation property. The introduced procedure can effectively solve the link division issue that had not been addressed by the previous study [1].
To demonstrate the applicability of the proposed OkNVD construction algorithm, a comprehensive case study in Shenzhen city was carried out. Results of the case study demonstrated that the proposed OkNVD construction algorithm can well quantify placebased accessibility in a fine-gained spatial resolution. Results of computational experiments showed that the proposed algorithm was significantly faster than previous Okabe's al-gorithm for all five testing networks under different k and |F| settings. For example, it performed about 65 times faster than Okabe's algorithm in the Shenzhen network when k = 3 and |F| = 1000. This computational advantage would be more obvious in real large-scale applications, in which the number of generators |F| are large and the k value is small.
Several directions for future research are worth pursuing. Firstly, this study assumed that link costs are static and deterministic. Traffic conditions in real road networks are dynamic and stochastic [29]. How to incorporate such dynamic and stochastic travel times into the proposed algorithm needs further investigation. Secondly, the proposed algorithm considers only the road transport mode. The extension of this proposed algorithm into the multi-mode transport networks, such as metro and bus, is warranted for further study. Last but not least, finding k nearest generators is a critical step for many network analyses, such as spatial interpolation, spatial clustering, outlier detection, and k nearest neighbor query, etc. [12,14,28]. Subsequent studies should investigate how to incorporate the proposed algorithm in these network analysis methods.

Data Availability Statement:
The data that support the findings of this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Notations used in this article.
A a set of links in the road network |A| number of links in the road network a ij a network link from tail node n i to head node n j a iq a partial link from tail node n i to location q a qj a partial link from location q to head node n j CRS K−1,i out the child region set of vr K−1,i out F the set of generators k key parameter, such as k nearest generators LS j label set of node n j LS j [i] a label for recording the shortest path from the ith generator f i N a set of nodes in the road network |N| number of nodes in the road network n i a network node o a network location os j open status p oq the shortest path from a location o from another location q p f w a temporal path from generator f i to node n w PRED n j a set of predecessor nodes of node n j q a network location SE priority queue for sorting shortest paths SG j the set of determined generators of node n j SG j [i] the ith determined generator of node n j SUCC n j a set of successor nodes of node n j t ij cost of link a ij t oq cost of path p oq t Proof. We first prove that Voronoi regions of the order first outward NVD are collectively exhaustive as follows. Because the road network is fully connected, we can calculate shortest paths from all generators to a location q. Then, we can determine the location's first outward nearest generator, f q,1 out ∈ F, satisfying t xq out . According to Equation (6), the location q belongs to a Voronoi region vr 1,i Out (F q,1 out ). Therefore, any location q in the road network is covered by VR 1 out . Thus, the collectively exhaustive property of VR 1 out is satisfied. We then prove that Voronoi regions of the order first outward NVD are mutually exclusive except at the boundary points as below. As proved, we can determine the first outward nearest generator of location q, i.e., f q,1 out ∈ F satisfying t xq 1 ≤ t xq i , ∀ f i ∈ F − f q,1 out , which can further divide into two scenarios. The first scenario is t out ) for ∀ f j ∈ . . . , f j,1 out , . . . . Therefore, Voronoi regions of the order first outward NVD are mutually exclusive except at the boundary points. Therefore, the order first outward VND forms a tessellation on the road network. Similarly, we also can prove that the order first inward VND forms a tessellation on the road network by using the shortest paths from a location q to all generators. Proposition A2. If the road network is fully connected, CRS K−1,i out forms a tessellation on vr K−1,i out .
Proof. We first prove that Voronoi regions of CRS , satisfying t f q . According to Equation (6), the location q belongs to a Voronoi region vr K,i out ⊂ CRS K−1,i out . Therefore, any location q in vr K−1,i out is covered by CRS K−1,i out . Thus, the collectively exhaustive property of CRS K−1,i out is satisfied. We then prove that Voronoi regions of the CRS K−1,i out are mutually exclusive except at the boundary points as below. As proved, we can determine the kth outward nearest generator of location q, i.e., f q,K out ∈ F − F w,K−1 out satisfying t f q out , which can further divide into two scenarios. The first scenario is t f q According to Equation (6), the location q belongs to only one Voronoi region, i.e., vr K,i out . Without the loss of generality, the other scenario is t f q