A New Algorithm for Euclidean Shortest Paths in the Plane
AA New Algorithm for Euclidean Shortest Paths in the Plane ∗ Haitao WangDepartment of Computer ScienceUtah State University, Logan, UT 84322, USA [email protected]
Abstract
Given a set of pairwise disjoint polygonal obstacles in the plane, finding an obstacle-avoidingEuclidean shortest path between two points is a classical problem in computational geometry andhas been studied extensively. Previously, Hershberger and Suri [SIAM J. Comput. 1999] gave analgorithm of O ( n log n ) time and O ( n log n ) space, where n is the total number of vertices of allobstacles. Recently, by modifying Hershberger and Suri’s algorithm, Wang [SODA 2021] reducedthe space to O ( n ) while the runtime of the algorithm is still O ( n log n ). In this paper, we presenta new algorithm of O ( n + h log h ) time and O ( n ) space, provided that a triangulation of the freespace is given, where h is the number of obstacles. The algorithm, which improves the previouswork when h = o ( n ), is optimal in both time and space as Ω( n + h log h ) is a lower bound on theruntime. Our algorithm builds a shortest path map for a source point s , so that given any querypoint t , the shortest path length from s to t can be computed in O (log n ) time and a shortest s - t path can be produced in additional time linear in the number of edges of the path. Let P be a set of h pairwise disjoint polygonal obstacles with a total of n vertices in the plane. Let F denote the free space , i.e., the plane minus the interior of the obstacles. Given two points s and t in F , we consider the problem of finding a Euclidean shortest path from s to t in F . This is a classicalproblem in computational geometry and has been studied extensively, e.g., [4, 19, 21–23, 25, 26, 33, 35,37, 41–43, 46].To solve the problem, two methods are often used in the literature: the visibility graph and thecontinuous Dijkstra. The visibility graph method is to first construct the visibility graph of the verticesof P along with s and t , and then run Dijkstra’s shortest path algorithm on the graph to find a shortest s - t path. The best algorithms for constructing the visibility graph run in O ( n log n + K ) time [19] or in O ( n + h log (cid:15) h + K ) time [8] for any constant (cid:15) >
0, where K is the number of edges of the visibilitygraph. Because K = Ω( n ) in the worst case, the visibility graph method inherently takes quadratictime. To deal with the case where h is relatively small comparing to n , a variation of the visibilitygraph method was proposed that is to first construct a so-called tangent graph and then find a shortest s - t path in the graph. Using this method, a shortest s - t path can be found in O ( n + h log h + K (cid:48) )time [7] after the free space F is triangulated, where K (cid:48) may be considered as the number of tangentsamong obstacles of P and K (cid:48) = O ( h ). Note that triangulating F can be done in O ( n log n ) time orin O ( n + h log (cid:15) h ) time for any small (cid:15) > O ( n / (cid:15) ) time for any constant (cid:15) >
0. Also using the continuous Dijkstra ∗ A preliminary version will appear in
Proceedings of the 53rd Annual ACM Symposium on Theory of Computing(STOC 2021) . This research was supported in part by NSF under Grant CCF-2005323. a r X i v : . [ c s . C G ] F e b pproach plus a novel conforming subdivision of the free space, Hershberger and Suri [26] presentedan algorithm of O ( n log n ) time and O ( n log n ) space; the running time is optimal when h = Θ( n ) asΩ( n + h log h ) is a lower bound in the algebraic computation tree model (which can be obtained bya reduction from sorting; e.g., see Theorem 3 [15] for a similar reduction). Recently, by modifyingHershberger and Suri’s algorithm, Wang [46] reduced the space to O ( n ) while the running time is still O ( n log n ). All three continuous Dijkstra algorithms [26, 37, 46] construct the shortest path map , denoted by
SPM ( s ), for a source point s . SPM ( s ) is of O ( n ) size and can be used to answer shortest path queries.By building a point location data structure on SPM ( s ) in additional O ( n ) time [17, 31], given a querypoint t , the shortest path length from s to t can be computed in O (log n ) time and a shortest s - t pathcan be output in time linear in the number of edges of the path.The problem setting for P is usually referred to as polygonal domains or polygons with holes inthe literature. The problem in simple polygons is relatively easier [21–23, 25, 33]. Guibas et al. [22]presented an algorithm that can construct a shortest path map in linear time. For two-point shortestpath query problem where both s and t are query points, Guibas and Hershberger [21, 23] built adata structure in linear time such that each query can be answered in O (log n ) time. In contrast, thetwo-point query problem in polygonal domains is much more challenging: to achieve O (log n ) timequeries, the current best result uses O ( n ) space [11]; alternatively Chiang and Mitchell [11] gave adata structure of O ( n + h ) space with O ( h log n ) query time. Refer to [11] for other data structureswith trade-off between space and query time.The L counterpart of the problem where the path length is measured in the L metric alsoattracted much attention, e.g., [2, 10, 12, 13, 34, 36]. For polygons with holes, Mitchell [34, 36] gave analgorithm that can build a shortest path map for a source point in O ( n log n ) time and O ( n ) space;for small h , Chen and Wang [10] proposed an algorithm of O ( n + h log h ) time and O ( n ) space, afterthe free space is triangulated. For simple polygons, Bae and Wang [2] built a data structure in lineartime that can answer each two-point L shortest path query in O (log n ) time. The two-point queryproblem in polygons with holes has also been studied [5, 6, 45]. To achieve O (log n ) time queries, thecurrent best result uses O ( n + h log h/ log log h ) space [45]. In this paper, we show that the problem of finding an Euclidean shortest path among obstacles in P is solvable in O ( n + h log h ) time and O ( n ) space, after a triangulation of the free space F is given. Ifthe time for triangulating F is included and the triangulation algorithm in [3] is used, then the totaltime of the algorithm is O ( n + h log (cid:15) h ), for any constant (cid:15) > With the assumption that thetriangulation could be done in O ( n + h log h ) time, which has been an open problem and is beyond thescope of this paper, our result settles Problem 21 in The Open Problem Project [1]. Our algorithmactually constructs the shortest path map SPM ( s ) for the source point s in O ( n + h log h ) time and O ( n ) space. We give an overview of our approach below.The high-level scheme of our algorithm is similar to that for the L case [45] in the sense that wefirst solve the convex case where all obstacles of P are convex and then extend the algorithm to thegeneral case with the help of the extended corridor structure of P [5, 8–10, 30, 38]. The convex case.
We first discuss the convex case. Let V denote the set of topmost, bottommost,leftmost, and rightmost vertices of all obstacles. Hence, |V| ≤ h . Using the algorithm of Hershberger An unrefereed report [29] announced an algorithm based on the continuous Dijkstra approach with O ( n + h log h log n )time and O ( n ) space. If randomization is allowed, the algorithm of Clarkson, Cole, and Tarjan [14] can compute a triangulation in O ( n log ∗ n + h log h ) expected time. S on the points of V , without considering the obstacleedges. Since |V| = O ( h ), the size of S is O ( h ). Then, we insert the obstacle edges into S to builda conforming subdivision S (cid:48) of the free space. The subdivision S (cid:48) has O ( h ) cells (in contrast, theconforming subdivision of the free space in [26] has O ( n ) cells). Unlike the subdivision in [26] whereeach cell is of constant size, here the size of each cell of S (cid:48) may not be constant but its boundaryconsists of O (1) transparent edges and O (1) convex chains (each of which belongs to the boundary ofan obstacle of P ). Like the subdivision in [26], each transparent edge e of S (cid:48) has a well-covering region U ( e ). In particular, for each transparent edge f on the boundary of U ( e ), the shortest path distancebetween e and f is at least 2 · max {| e | , | f |} . Using S (cid:48) as a guidance, we run the continuous Dijkstraalgorithm as in [26] to expand the wavefront, starting from the source point s . A main challengeour algorithm needs to overcome (which is also a main difference between our algorithm and thatin [26]) is that each cell in our subdivision S (cid:48) may not be of constant size. One critical property ouralgorithm relies on is that the boundary of each cell of S (cid:48) has O (1) convex chains. Our strategy is tosomehow treat each such convex chain as a whole. We also borrow some idea from the algorithm ofHershberger, Suri, and Yıldız [27] for computing shortest paths among curved obstacles. To guaranteethe O ( n + h log h ) time, some global charging analysis is used. In addition, the tentative prune-and-search technique of Kirkpatrick and Snoeyink [32] is applied to perform certain operations related tobisectors, in logarithmic time each. Finally, the techniques of Wang [46] are utilized to reduce thespace to O ( n ). All these efforts lead to an O ( n + h log h ) time and O ( n ) space algorithm to constructthe shortest path map SPM ( s ) for the convex case. The general case.
We extend the convex case algorithm to the general case where obstacles maynot be convex. To this end, we resort to the extended corridor structure of P , which was used beforefor reducing the time complexities from n to h , e.g., [5, 8–10, 30, 38]. The structure partitions the freespace F into an ocean M , O ( n ) bays, and O ( h ) canals. Each bay is a simple polygon that shares anedge with M . Each canal is a simple polygon that shares two edges with M . But two bays or twocanals, or a bay and a canal do not share any edge. A common edge of a bay (or canal) with M iscalled a gate. Thus each bay has one gate and each canal has two gates. Further, M is bounded by O ( h ) convex chains (each of which is on the boundary of an obstacle). An important property relatedto shortest paths is that if both s and t are in M , then any shortest s - t path must be in the union of M and all corridor paths, each of which is contained in a canal. As the boundary of M consists of O ( h )convex chains, by incorporating all corridor paths, we can easily extend our convex case algorithm tocomputing SPM ( M ), the shortest path map SPM ( s ) restricted to M , i.e., SPM ( M ) = SPM ( s ) ∩ M .To compute the entire map SPM ( s ), we expand SPM ( M ) to all bays and canals through their gates.For this, we process each bay/canal individually. For each bay/canal C , expanding the map into C is actually a special case of the (additively)-weighted geodesic Voronoi diagram problem on a simplepolygon where all sites are outside C and can influence C only through its gates. In summary, aftera triangulation of F is given, building SPM ( M ) takes O ( n + h log h ) time, and expanding SPM ( M )to all bays and canals takes additional O ( n + h log h ) time. The space of the algorithm is bounded by O ( n ). Outline.
The rest of the paper is organized as follows. Section 2 defines notation and introducessome concepts. Section 3 presents the algorithm for the convex case. The general case is discussed inSection 4.
For any two points a and b in the plane, denote by ab the line segment with a and b as endpoints;denote by | ab | the length of the segment. 3 Figure 1:
Illustrating the shortest path map
SPM ( s ). The solid (red) curves are walls and the (blue) dotted segmentsare windows. The anchor of each cell is also shown with a black point. For any two points s and t in the free space F , we use π ( s, t ) to denote a shortest path from s to t in F . In the case where shortest paths are not unique, π ( s, t ) may refer to an arbitrary one. Denote by d ( s, t ) the length of π ( s, t ); we call d ( s, t ) the geodesic distance between s and t . For two line segments e and f in F , their geodesic distance is defined to be the minimum geodesic distance between anypoint on e and any point on f , i.e., min s ∈ e,t ∈ f d ( s, t ); by slightly abusing the notation, we use d ( e, f )to denote their geodesic distance. For any path π in the plane, we use | π | to denote its length.For any compact region A in the plane, let ∂A denote its boundary. We use ∂ P to denote theunion of the boundaries of all obstacles of P .Throughout the paper, we use s to refer to the source point. For convenience, we consider s as adegenerate obstacle in P . We often refer to the vertices of P as obstacle vertices and refer to the edgesof P as obstacle edges . For any point t ∈ F , we call the adjacent vertex of t in π ( s, t ) the anchor of t in π ( s, t ) .The shortest path map SPM ( s ) of s is a decomposition of the free space F into maximal regionssuch that all points in each region R have the same anchor [26, 35] in their shortest paths from s ; e.g.,see Fig. 1. Each edge of SPM ( s ) is either an obstacle edge fragment or a bisecting-curve , which is thelocus of points p with d ( s, u ) + | pu | = d ( s, v ) + | pv | for two obstacle vertices u and v . Each bisecting-curve is in general a hyperbola; a special case happens if one of u and v is the anchor of the other,in which case their bisecting-curve is a straight line. Following the notation in [18], we differentiatebetween two types of bisecting curves: walls and windows . A bisecting curve of SPM ( s ) is a wall ifthere exist two topologically different shortest paths from s to each point of the edge; otherwise (i.e.,the above special case) it is a window (e.g., see Fig. 1).We make a general position assumption that for each obstacle vertex v , there is a unique shortestpath from s to v , and for any point p in the plane, there are at most three different shortest pathsfrom s to p . The assumption assures that each vertex of SPM ( s ) has degree at most three, and thereare at most three bisectors of SPM ( s ) intersecting at a common point, which is sometimes called a triple point in the literature [18].A curve in the plane is x -monotone if its intersection with any vertical line is connected; the y -monotone is defined similarly. A curve is xy -monotone if it is both x - and y -monotone.The following observation will be used throughout the paper without explicitly mentioning it again. Usually “predecessor” is used in the literature instead of “anchor”, but here we reserve “predecessor” for otherpurpose. This is usually called bisector in the literature. Here we reserve the term “bisector” to be used later. u Figure 2:
Illustrating an elementary chain (the thicksegments), which contains the vertex v but not u . ppred ( p ) s ppred ( p ) s (a) (b) Figure 3:
Illustrating the predecessor.
Observation 1 n + h log n = O ( n + h log h ) . Proof:
Indeed, if h < n/ log n , then n + h log n = Θ( n ), which is O ( n + h log h ); otherwise, log n = O (log h ) and n + h log n = O ( n + h log h ). (cid:3) In this section, we present our algorithm for the convex case where all obstacles of P are convex. Thealgorithm will be extended to the general case in Section 4.For each obstacle P ∈ P , the topmost, bottommost, leftmost, and rightmost vertices of P are called rectilinear extreme vertices . The four rectilinear extreme vertices partition ∂P into four portions andeach portion is called an elementary chain , which is convex and xy -monotone. For technical reasonthat will be clear later, we assume that each rectilinear extreme vertex v belongs to the elementarychain counterclockwise of v with respect to the obstacle (i.e., v is the clockwise endpoint of the chain;e.g., see Fig. 2). We use elementary chain fragment to refer to a portion of an elementary chain.We introduce some notation below that is similar in spirit to those from [27] for shortest pathsamong curved obstacles.Consider a shortest path π ( s, p ) from s to a point p in the free space F . It is not difficult tosee that π ( s, p ) is a sequence of elementary chain fragments and common tangents between obstaclesof P ∪ { p } . We define the predecessor of p , denoted by pred ( p ), to be the initial vertex of the lastelementary chain fragment in π ( s, p ) (e.g., see Fig. 3 (a)). Note that since each rectilinear extremevertex belongs to a single elementary chain, pred ( p ) in π ( s, p ) is unique. A special case happens if p is a rectilinear extreme vertex and π ( s, p ) contains a portion of an elementary chain A clockwiseof p . In this case, we let pred ( p ) be endpoint of the fragment of A in π ( s, p ) other than p (e.g., seeFig. 3 (b)); in this way, pred ( p ) is unique in π ( s, p ). Note that p may still have multiple predecessorsif there are multiple shortest paths from s to p . Intuitively, the reason we define predecessors as aboveis to treat each elementary chain somehow as a whole, which is essential for reducing the runtime ofthe algorithm from n to h .The rest of this section is organized as follows. In Section 3.1, we compute a conforming subdivision S (cid:48) of the free space F . Section 3.2 introduces some basic concepts and notation for our algorithm. Thewavefront expansion algorithm is presented in Section 3.3, with two key subroutines of the algorithmdescribed in Section 3.4 and Section 3.5, respectively. Section 3.6 analyzes the time complexity ofthe algorithm, where a technical lemma is proved separately in Section 3.7. Using the informationcomputed by the wavefront expansion algorithm, Section 3.8 constructs the shortest path map SPM ( s ).The overall algorithm runs in O ( n + h log h ) time and O ( n + h log h ) space. Section 3.9 reduces thespace to O ( n ) while keeping the same runtime, by using the techniques from Wang [46].5 .1 Computing a conforming subdivision of the free space Let V denote the set of the rectilinear extreme vertices of all obstacles of P . Hence, |V| = O ( h ).Using the algorithm of algorithm of Hershberger and Suri [26] (called the HS algorithm), we build aconforming subdivision S with respect to the vertices of V , without considering the obstacle edges.The subdivision S , which is of size O ( h ), is a quad-tree-style subdivision of the plane into O ( h )cells. Each cell of S is a square or a square annulus (i.e., an outer square with an inner square hole).Each vertex of V is contained in the interior of a square cell and each square cell contains at most onevertex of V . Each edge e of S is axis-parallel and well-covered , i.e., there exists a set C ( e ) of O (1) cellsof S such that their union U ( e ) contains e with the following properties: (1) the total complexity ofall cells of C ( e ) is O (1) and thus the size of U ( e ) is O (1); (2) for any edge f of S that is on ∂ U ( e ) oroutside U ( e ), the Euclidean distance between e and f (i.e., the minimum | pq | among all points p ∈ e and q ∈ f ) is at least 2 · max {| e | , | f |} ; (3) U ( e ), which is called the well-covering region of e , contains atmost one vertex of V . In addition, each cell c of S has O (1) edges on its boundary with the following uniform edge property : the lengths of the edges on the boundary of c differ by at most a factor of 4,regardless of whether c is a square or square annulus.The subdivision S can be computed in O ( h log h ) time and O ( h ) space [26].Next we insert the obstacle edges into S to produce a conforming subdivision S (cid:48) of the free space F . In S (cid:48) , there are two types of edges: those introduced by the subdivision construction (which arein the interior of F except possibly their endpoints) and the obstacle edges; we call the former the transparent edges (which are axis-parallel) and the latter the opaque edges . The definition of S (cid:48) issimilar to the conforming subdivision of the free space used in the HS algorithm. A main difference isthat here endpoints of each obstacle edge may not be in V , a consequence of which is that each cell of S (cid:48) may not be of constant size (while each cell in the subdivision of the HS algorithm is of constantsize). However, each cell c of S (cid:48) has the following property that is critical to our algorithm: Theboundary ∂c consists of O (1) transparent edges and O (1) convex chains (each of which is a portion ofan elementary chain).More specifically, S (cid:48) is a subdivision of F into O ( h ) cells. Each cell of S (cid:48) is one of the connectedcomponents formed by intersecting F with an axis-parallel rectangle (which is the union of a set ofadjacent cells of S ) or a square annulus of S . Each cell of S (cid:48) contains at most one vertex of V . Eachvertex of V is incident to a transparent edge. Each transparent edge e of S (cid:48) is well-covered , i.e.,there exists a set C ( e ) of O (1) cells whose union U ( e ) contains e with the following property: for eachtransparent edge f on ∂ U ( e ), the geodesic distance d ( e, f ) between e and f is at least 2 · max {| e | , | f |} .The region U ( e ) is called the well-covering region of e and contains at most one vertex of V . Note that S (cid:48) has O ( h ) transparent edges.Below we show how S (cid:48) is produced from S . The procedure is similar to that in the HS algorithm.We overlay the obstacle edges on top of S to obtain a subdivision S overlay . Because each edge of S is axis-parallel and all obstacle edges constitute a total of O ( h ) elementary chains, each of which is xy -monotone, S overlay has O ( h ) faces. We say that a face of S overlay is interesting if its boundarycontains a vertex of V or a vertex of S . We keep intact the interesting faces of S overlay while deletingevery edge fragment of S not on the boundary of any interesting cell. Further, for each cell c containinga vertex v ∈ V , we partition c by extending vertical edges from v until the boundary of c . This divides c into at most three subcells. Finally we divide each of the two added edges incident to v into segmentsof length at most δ , where δ is the length of the shortest edge on the boundary of c . By the uniformedge property of S , ∂c has O (1) edges, whose lengths differ by at most a factor of 4; hence dividingthe edges incident to v as above produces only O (1) vertical edges. The resulting subdivision is S (cid:48) .As mentioned above, the essential difference between our subdivision S (cid:48) and the one in the HSalgorithm is that the role of an obstacle edge in the HS algorithm is replaced by an elementary chain.Therefore, each opaque edge in the subdivision of the HS algorithm becomes an elementary chainfragment in our case. Hence, by the same analysis as in the HS algorithm (see Lemma 2.2 [26]), S (cid:48) U ( e ) of each transparent edge e of S (cid:48) is defined in the same way as in the HS algorithm.The following lemma computes S (cid:48) . It should be noted that although S (cid:48) is defined with the help of S overlay , S (cid:48) is constructed directly without computing S overlay first. Lemma 1
The conforming subdivision S (cid:48) can be constructed in O ( n + h log h ) time and O ( n ) space. Proof:
We first construct S in O ( h log h ) time and O ( h ) space [26]. In the following we construct S (cid:48) by inserting the obstacle edges into S . The algorithm is similar to that in the HS algorithm (seeLemma 2.3 [26]). The difference is that we need to handle each elementary chain as a whole.We first build a data structure so that for any query horizontal ray with origin in F , the firstobstacle edge of P hit by it can be computed in O (log n ) time. This can be done by building ahorizontal decomposition of F , i.e., extend a horizontal segment from each vertex until it hits ∂ P . Asall obstacles of P are convex, the horizontal decomposition can be computed in O ( n + h log h ) time and O ( n ) space [28]. By building a point location data structure [17, 31] on the horizontal decompositionin additional O ( n ) time, each horizontal ray shooting can be answered in O (log n ) time. Similarly,we can construct the vertical decomposition of F in O ( n + h log h ) time and O ( n ) space so that eachvertical ray shooting can be answered in O (log n ) time.The edges of S (cid:48) are obstacle edges, transparent edges incident to the vertices of S , and transparentedges subdivided on the vertical segments incident to the vertices of V . To identify the second type ofedges, we trace the boundary of each interesting cell separately. Starting from a vertex v of S , we tracealong each edge incident to v . Using the above ray-shooting data structure, we determine whetherthe next cell vertex is a vertex of S or the first point on ∂ P hit by the ray. As S has O ( h ) edges andvertices, this tracing takes O ( h log n ) time in total. Tracing along obstacle edges is done by startingfrom each vertex of V and following each of its incident elementary chains. For each elementary chain,the next vertex is either the next obstacle vertex on e , where e is the current tracing edge of theelementary chain, or the intersection of e with a transparent edge of the current cell. Hence, tracingall elementary chains takes O ( n ) time in total. The third type of edges can be computed in lineartime by local operations on each cell containing a vertex of V .Finally, we assemble all these edges together to obtain an adjacency structure for S (cid:48) . For this, wecould use a plane sweep algorithm as in the HS algorithm. However, that would take O ( n log n ) timeas there are O ( n ) edges in S (cid:48) . To obtain an O ( n + h log h ) time algorithm, we propose the followingapproach. During the above tracing of elementary chains, we record the fragment of the chain that liesin a single cell. Since each such portion may not be of constant size, we represent it by a line segmentconnecting its two endpoints; note that this segment does not intersect any other cell edges becausethe elementary chain fragment is on the boundary of an obstacle (and thus the segment is inside theobstacle). Then, we apply the plane sweep algorithm to stitch all edges with each elementary chainfragment in a cell replaced by a segment as above. The algorithm takes O ( h log h ) time and O ( h )space. Finally, for each segment representing an elementary chain portion, we locally replace it bythe chain fragment in linear time. Hence, this step takes O ( n ) time altogether for all such segments.As such, the total time for computing the adjacency information for S (cid:48) is O ( n + h log n ), which is O ( n + h log h ). Clearly, the space complexity of the algorithm is O ( n ). (cid:3) Our shortest path algorithm uses the continuous Dijkstra method. The algorithm initially generatesa wavefront from s , which is a circle centered at s . During the algorithm, the wavefront consists of allpoints of F with the same geodesic distance from s (e.g., see Fig. 4). We expand the wavefront until allpoints of the free space are covered. The conforming subdivision S (cid:48) is utilized to guide the wavefront7igure 4: Illustrating the wavefront. The black region are obstacles. The green point is s . The red curves are bisectingcurves of SPM ( s ). The gray region is the free space that has been covered by the wavefront. The boundary between thewhite region and the grey region is the wavefront. The figure is generated using the applet at [24]. expansion. Our wavefront expansion algorithm follows the high-level scheme as the HS algorithm.The difference is that our algorithm somehow considers each elementary chain as a whole, which is insome sense similar to the algorithm of Hershberger, Suri, and Yıldız [27] (called the HSY algorithm).Indeed, the HSY algorithm considers each xy -monotone convex arc as a whole, but the difference isthat each arc in the HSY algorithm is of constant size while in our case each elementary chain maynot be of constant size. As such, techniques from both the HS algorithm and the HSY algorithm areborrowed.We use τ to denote the geodesic distance from s to all points in the wavefront. One may alsothink of τ as a parameter representing time. The algorithm simulates the expansion of the wavefrontas time increases from 0 to ∞ . The wavefront comprises a sequence of wavelets , each emanating froma generator . In the HS algorithm, a generator is simply an obstacle vertex. Here, since we want totreat an elementary chain as a whole, similar to the HSY algorithm, we define a generator as a couple α = ( A, a ), where A is an elementary chain and a is an obstacle vertex on A , and further a clockwiseor counterclockwise direction of A is designated for α ; a has a weight w ( a ) (one may consider w ( a ) asthe geodesic distance between s and a ). We call a the initial vertex of the generator α .We say a point q is reachable by a generator α = ( A, a ) if one can draw a path in F from a to q by following A in the designated direction to a vertex v on A such that vq is tangent to A and thenfollowing the segment vq (e.g., see Fig. 5). The (weighted) distance between the generator α and q is the length of this path plus w ( a ); by slightly abusing the notation, we use d ( α, q ) to denote thedistance. From the definition of reachable points, the vertex a partitions A into two portions and onlythe portion following the designated direction is relevant (e.g., in Fig. 5, only the portion containingthe vertex v is relevant). Henceforth, unless otherwise stated, we use A to refer to its relevant portiononly and we call A the underlying chain of α . In this way the initial vertex a becomes an endpointof A . For convenience, sometimes we do not differentiate α and A . For example, when we say “thetangent from q to α ”, we mean “the tangent from q to A ”; also, when we say “a vertex of α ”, we mean“a vertex of A ”.The wavefront can thus be represented by the sequence of generators of its wavelets. A waveletgenerated by a generator α = ( A, a ) at time τ is a contiguous set of reachable points q such that d ( α, q ) = τ and d ( α (cid:48) , q ) ≥ τ for all other generators α (cid:48) in the wavefront; we also say that q is claimedby α . Note that as A may not be of constant size, a wavelet may not be of constant size either; itactually consists of a contiguous sequence of circular arcs centered at the obstacle vertices A (e.g., seeFig. 6). If a point q is claimed by α , then d ( s, q ) = d ( α, q ) = τ and the predecessor pred ( q ) of q is a ; sometimes for convenience we also say that the generator α is the predecessor of q . If q is on anelementary chain A (cid:48) and the tangent from q to a is also tangent to A (cid:48) , then a new generator ( A (cid:48) , q )8 a A v Figure 5:
Illustrating a generator α = ( A, a ). A consistsof the thick segments, and is designated clockwise directionaround the obstacle. q is a reachable point through vertex v . aA Figure 6:
Illustrating the wavelet generated by α =( A, q ), designated counterclockwise direction. Thewavelet has three pieces each of which is a circular arc. is added to the wavefront (e.g., see Fig. 7(a)). A special case happens when q is the counterclockwiseendpoint of A (and thus q does not belong to A ); in this case, a new generator α (cid:48) = ( A (cid:48) , q ) is alsoadded, where A (cid:48) is the elementary chain that contains q (e.g., see Fig. 7(b)). qa A A (cid:48) AA (cid:48) a q ( a ) ( b ) v Figure 7:
Illustrating a new generator α (cid:48) = ( A (cid:48) , q ). (a) A general case where both A and A (cid:48) are marked with thicksegments. (b) A special case where A is marked with solid (red) segments and A (cid:48) is marked with dashed (blue) segments. As τ increases, the points bounding the adjacent wavelets trace out the bisectors that form theedges of the shortest path map SPM ( s ). The bisector between the wavelets of two generators α and α (cid:48) ,denoted by B ( α, α (cid:48) ), consists of points q with d ( α, q ) = d ( α (cid:48) , q ). Note that since α and α (cid:48) may not beof constant size, B ( α, α (cid:48) ) may not be of constant size either. More specifically, B ( α, α (cid:48) ) has multiplepieces each of which is a hyperbola defined by two obstacle vertices v ∈ α and v (cid:48) ∈ α (cid:48) such that thehyperbola consists of all points that have two shortest paths from s with v and v (cid:48) as the anchors inthe two paths, respectively (e.g., see Fig. 8). A special case happens if α (cid:48) is a generator created bythe wavelet of α , such as that illustrated in Fig. 7(a), then B ( α, α (cid:48) ) is the half-line extended from q along the direction from v to q (the dashed segment in Fig. 7(a)); we call such a bisector an extensionbisector . Note that in the case illustrated in Fig. 7(b), B ( α, α (cid:48) ), which is also an extension bisector, isthe half-line extended from q along the direction from v to q (the dashed segment in Fig. 7(b)), where v is the obstacle vertex adjacent to q in A .A wavelet gets eliminated from the wavefront if the two bisectors bounding it intersect, which iscalled a bisector event . Specifically, if α , α , and α are three consecutive generators of the wavefront,the wavelet generated by α will be eliminated when B ( α , α ) intersects B ( α, α ); e.g., see Fig. 9.Wavelets are also eliminated by collisions with obstacles and other wavelets in front of it. If a bisector B ( α, α (cid:48) ) intersects an obstacle, their intersection is also called a bisector event .Let SPM (cid:48) ( s ) be the subdivision of F by the bisectors of all generators (e.g. see Fig. 10). Theintersections of bisectors and intersections between bisectors and obstacle edges are vertices of SPM (cid:48) ( s ).Each bisector connecting two vertices is an edge of SPM (cid:48) ( s ), called a bisector edge . As discussed before,9 A B ( α, α (cid:48) ) A (cid:48) a (cid:48) v v (cid:48) qq (cid:48) Figure 8:
Illustrating the bisector B ( γ, γ (cid:48) ) defined by two generators γ = ( A, q ) and γ = ( A (cid:48) , q (cid:48) ). The portion between q and q (cid:48) is a hyperbola defined by v and v (cid:48) . αα α B ( α , α ) B ( α, α ) Figure 9:
Illustrating the intersection of two bisectors. a bisector, which consists of multiple hyperbola pieces, may not be of constant size. Hence, a bisectoredge e of SPM (cid:48) ( s ) may not be of constant size. For differentiation, we call each hyperbola piece of e , a hyperbolic-arc . In addition, if the boundary of an obstacle P contains more than one vertex of SPM (cid:48) ( s ), then the chain of edges of ∂P connecting two adjacent vertices of SPM (cid:48) ( s ) also forms anedge of SPM (cid:48) ( s ), called a convex-chain edge .The following lemma shows that SPM (cid:48) ( s ) is very similar to SPM ( s ). Refer to Fig. 11 for SPM ( s ). Lemma 2
Each extension bisector edge of
SPM (cid:48) ( s ) is a window of SPM ( s ) . The union of all non-extension bisector edges of SPM (cid:48) ( s ) is exactly the union of all walls of SPM ( s ) . SPM (cid:48) ( s ) can beobtained from SPM ( s ) by removing all windows except those that are extension bisectors of SPM (cid:48) ( s ) . Proof:
We first prove that each extension bisector edge e of SPM (cid:48) ( s ) is a window of SPM ( s ). Bydefinition, one endpoint of e , denoted by u , is an obstacle vertex, such that e is a half-line extendedfrom u in the direction from v to u for another obstacle vertex v . Hence, for each point p ∈ e , there isa shortest s - p path π ( s, p ) such that vp is a segment of π ( s, p ). As u ∈ vp , p must be on a window w of SPM ( s ). This proves that e ⊆ w . On the other hand, for any point p ∈ w , there is a shortest s - p path π ( s, p ) that contains vp . Since u ∈ vp , p must be on the extension bisector edge e of SPM (cid:48) ( s ),and thus w ⊆ e . Therefore, e = w . This proves that each extension bisector edge of SPM (cid:48) ( s ) is awindow of SPM ( s ).We next prove that the union of all non-extension bisector edges of SPM (cid:48) ( s ) is exactly the unionof all walls of SPM ( s ). • We first show that the union of all non-extension bisector edges of
SPM (cid:48) ( s ) is a subset of theunion of all walls of SPM ( s ). 10 qa Figure 10:
Illustrating the map
SPM (cid:48) ( s ). The red curves are non-extension bisectors. The green and blue segmentsare extension bisectors, where the green and blue ones are Type (a) and (b) as illustrated in Fig. 7. The predecessor ineach cell is also shown with a black point. For example, all points in the cell containing q have a as their predecessor. s Figure 11:
Illustrating the map
SPM ( s ). The dashed segments are windows that are not in SPM (cid:48) ( s ); removing thembecomes SPM (cid:48) ( s ). The anchor of each cell is also shown. Let q be a point on a non-extension bisector B ( α, α (cid:48) ) of two generators α = ( A, a ) and α =( A (cid:48) , a (cid:48) ). Let v and v (cid:48) be the tangent points on A and A (cid:48) from q , respectively. Hence, q hastwo shortest paths from s , one containing qv and the other containing qv (cid:48) . Since B ( α, α (cid:48) ) is notan extension bisector, v (cid:48) (cid:54)∈ qv and v (cid:54)∈ qv (cid:48) . This means that the two paths are combinatoriallydifferent and thus q is on a wall of SPM ( s ). This proves that the union of all non-extensionbisector edges of SPM (cid:48) ( s ) is a subset of the union of all walls of SPM ( s ). • We then show that the union of all walls of
SPM ( s ) is a subset of the union of all non-extensionbisector edges of SPM (cid:48) ( s ).Let q be a point on a wall of SPM ( s ). By definition, there are two obstacle vertices v and v (cid:48) suchthat there are two combinatorially different shortest paths from s to q whose anchors are v and v (cid:48) , respectively. Further, since the two paths are combinatorially different and also due to thegeneral position assumption, v (cid:48) (cid:54)∈ qv and v (cid:54)∈ qv (cid:48) . Therefore, q must be on the bisector edge of thetwo generators α and α (cid:48) in SPM (cid:48) ( s ) whose underlying chains containing v and v (cid:48) , respectively.Further, since v (cid:54)∈ qv (cid:48) and v (cid:48) (cid:54)∈ qv , the bisector of α and α (cid:48) is not an extension bisector. This11roves that the union of all walls of SPM ( s ) is a subset of the union of all non-extension bisectoredges of SPM (cid:48) ( s ).The above proves that the first and second statements of the lemma. The third statement imme-diately follows the first two. (cid:3) Corollary 1
The combinatorial size of
SPM (cid:48) ( s ) is O ( n ) . Proof:
By Lemma 2,
SPM (cid:48) ( s ) is a subset of SPM ( s ). As the combinatorial size of SPM ( s ) is O ( n ) [26, 35], the combinatorial size of SPM (cid:48) ( s ) is also O ( n ). (cid:3) Lemma 3
The subdivision
SPM (cid:48) ( s ) has O ( h ) faces, vertices, and edges. In particular, SPM (cid:48) ( s ) has O ( h ) bisector intersections and O ( h ) intersections between bisectors and obstacle edges. Proof:
By Lemma 2, for any vertex of
SPM (cid:48) ( s ) that is the intersection of two non-exenstion bisectors,it is also the intersection of two walls of SPM ( s ), which is a triple point . Let SPM (cid:48)(cid:48) ( s ) refer to SPM (cid:48) ( s )excluding all extension bisectors. We define a planar graph G corresponding to SPM (cid:48)(cid:48) ( s ) as follows.Each obstacle defines a node of G and each triple point also defines a node of G . For any two nodes of G , if they are connected by a chain of bisector edges of SPM (cid:48) ( s ) such that the chain does not containany other node of G , then G has an edge connecting the two nodes. It is proved in [44] that G has O ( h ) vertices, faces, and edges.It is not difficult to see that a face of SPM (cid:48)(cid:48) ( s ) corresponds to a face of G and thus the numberof faces of SPM (cid:48)(cid:48) ( s ) is O ( h ). For each bisector intersection of SPM (cid:48)(cid:48) ( s ), it must be a triple point andthus it is also a node of G . As G has O ( h ) nodes, the number of bisector intersections of SPM (cid:48)(cid:48) ( s )is O ( h ). For each intersection v between a bisector and an obstacle edge in SPM (cid:48)(cid:48) ( s ), G must has anedge e corresponding a chain of bisector edges and v is the endpoint of the chain; we charge v to e .As e as two incident nodes of G , e can be charged at most twice. Since G has O ( h ) edges, the numberof intersections between bisectors and obstacle edges in SPM (cid:48)(cid:48) ( s ) is O ( h ).We next prove that the total number of extension bisector edges of SPM (cid:48) ( s ) is O ( h ). For eachextension bisection edge e , it belongs to one of the two cases ( a ) and ( b ) illustrated in Fig. 7. ForCase (b), one endpoint of e is a rectilinear extreme vertex. Since each rectilinear extreme vertex candefine at most two extension bisector edges and there are O ( h ) rectilinear extreme vertices, the totalnumber of Case (b) extension bisectors is O ( h ). For Case (a), e is an extension of a common tangentof two obstacles; we say that e is defined by the common tangent. Note that all such common tangentsare interior disjoint as they belong to shortest paths from s encoded by SPM (cid:48) ( s ). We now define aplanar graph G (cid:48) as follows. The node set of G (cid:48) consists of all obstacles of P . Two obstacles have anedge in G (cid:48) if they have at least one common tangent that defines an extension bisector. Since all suchcommon tangents are interior disjoint, G (cid:48) is a planar graph. Apparently, no two nodes of G (cid:48) sharetwo edges and no node of G (cid:48) has a self-loop. Therefore, G (cid:48) is a simple planar graph. Since G (cid:48) has h vertices, the number of edges of G (cid:48) is O ( h ). By the definition of G (cid:48) , each pair of obstacles whosecommon tangents define extension bisectors have an edge in G (cid:48) . Because each pair of obstacles candefine at most four common tangents and thus at most four extension bisectors and also because G (cid:48) has O ( h ) edges, the total number of Case (a) extension bisectors is O ( h ).The bisector intersections of SPM (cid:48) ( s ) that are not vertices of SPM (cid:48)(cid:48) ( s ) involve extension bisectors of SPM (cid:48) ( s ). The intersections between bisectors and obstacle edges in SPM (cid:48) ( s ) that are not in SPM (cid:48)(cid:48) ( s )also involve extension bisectors. Since all extension bisectors and all edges of SPM (cid:48)(cid:48) ( s ) are interiordisjoint, each extension bisector can involve in at most two bisector intersections and at most two12ntersections between bisectors and obstacle edges. Because the total number of extension bisectoredges of SPM (cid:48) ( s ) is O ( h ), the number of bisector intersections involving extension bisectors is O ( h ) andthe number of intersections between extension bisectors and obstacle edges is also O ( h ). Therefore,comparing to SPM (cid:48)(cid:48) ( s ), SPM (cid:48) ( s ) has O ( h ) additional vertices and O ( h ) additional edges. Note thatthe number of convex-chain edges of SPM (cid:48) ( s ) is O ( h ), for SPM (cid:48) ( s ) has O ( h ) vertices. The lemma thusfollows. (cid:3) Corollary 2
There are O ( h ) bisector events and O ( h ) generators in SPM (cid:48) ( s ) . Proof:
Each bisector event is either a bisector intersection or an intersection between a bisector andan obstacle edge. By Lemma 3, there are O ( h ) bisector intersections and O ( h ) intersections betweenbisectors and obstacle edges. As such, there are O ( h ) bisector events in SPM (cid:48) ( s ).By definition, each face of SPM (cid:48) ( s ) has a unique generator. Since SPM (cid:48) ( s ) has O ( h ) faces byLemma 3, the total number of generators is O ( h ). (cid:3) By the definition of
SPM (cid:48) ( s ), each cell C of SPM (cid:48) ( s ) has a unique generator α = ( A, a ), all pointsof the cell are reachable from the generator, and a is the predecessor of all points of C (e.g., see Fig. 10;all points in the cell containing q have a as their predecessor). Hence, for any point q ∈ C , we cancompute d ( s, q ) = d ( α, q ) by computing the tangent from q to A . Thus, SPM (cid:48) ( s ) can also be used toanswer shortest path queries. In face, given SPM (cid:48) ( s ), we can construct SPM ( s ) in additional O ( n )time by inserting the windows of SPM ( s ) to SPM (cid:48) ( s ), as shown in the lemma below. Lemma 4
Given
SPM (cid:48) ( s ) , SPM ( s ) can be constructed in O ( n ) time (e.g., see Fig. 10 and Fig. 11). Proof:
It suffices to insert all windows of
SPM ( s ) (except those that are also extension bisectors of SPM (cid:48) ( s )) to SPM (cid:48) ( s ). We consider each cell C of SPM (cid:48) ( s ) separately. Let α = ( A, a ) be the generatorof C . Our goal is to extend each obstacle edge of A along the designated direction of A until theboundary of C . To this end, since all points of C are reachable from α , the cell C is “star-shaped”with respect to the tangents of A (along the designated direction) and the extension of each obstacleedge of A intersects the boundary of C at most once. Hence, the order of the endpoints of theseextensions is consistent with the order of the corresponding edges of A . Therefore, these extensionendpoints can be found in order by traversing the boundary of C , which takes linear time in the sizeof C . Since the total size of all cells of SPM (cid:48) ( s ) is O ( n ) by Corollary 1, the total time of the algorithmis O ( n ). (cid:3) In light of Lemma 4, we will focus on computing the map
SPM (cid:48) ( s ). To simulate the wavefront expansion, we compute the wavefront passing through each transparentedge of the conforming subdivision S (cid:48) . As in the HS algorithm, since it is difficult to compute a truewavefront for each transparent edge e of S (cid:48) , a key idea is to compute two one-sided wavefronts (called approximate wavefronts ) for e , each representing the wavefront coming from one side of e . Intuitively,an approximate wavefront from one side of e is what the true wavefront would be if the wavefrontwere blocked off at e by considering e as an opaque segment (with open endpoints).In the following, unless otherwise stated, a wavefront at a transparent edge e of S (cid:48) refers to anapproximate wavefront. Let W ( e ) denote a wavefront at e . To make the description concise, asthere are two wavefronts at e , depending on the context, W ( e ) may refer to both wavefronts, i.e., thediscussion on W ( e ) applies to both wavefronts. For example, “compute the wavefronts W ( e )” means“compute both wavefronts at e ”. 13or each transparent edge e of S (cid:48) , define input ( e ) as the set of transparent edges on the boundaryof the well-covering region U ( e ), and define output ( e ) = { g | e ∈ input ( g ) } ∪ input ( e ) . Because ∂ U ( e (cid:48) )for each transparent edge e (cid:48) of S (cid:48) has O (1) transparent edges, both | input ( e ) | and | output ( e ) | are O (1). The wavefront propagation and merging procedures.
The wavefronts W ( e ) at e are computedfrom the wavefronts at edges of input ( e ); this guarantees the correctness because e is in U ( e ) (and thusany shortest path π ( s, p ) must cross some edge f ∈ input ( e ) for any point p ∈ e ). After the wavefronts W ( e ) at e are computed, they will pass to the edges of output ( e ). Also, the geodesic distances from s to both endpoints of e will be computed. Recall that V is the set of rectilinear extreme vertices ofall obstacles and each vertex of V is incident to a transparent edge of S (cid:48) . As such, after the algorithmis finished, geodesic distances from s to all vertices of V will be available. The process of passing thewavefronts W ( e ) at e to all edges g ∈ output ( e ) is called the wavefront propagation procedure , whichwill compute the wavefront W ( e, g ), where W ( e, g ) is the portion of W ( e ) that passes to g through thewell-covering region U ( g ) of g if e ∈ input ( g ) and through U ( e ) otherwise (in this case g ∈ input ( e ));whenever the procedure is invoked on e , we say that e is processed . The wavefronts W ( e ) at e areconstructed by merging the wavefronts W ( f, e ) for edges f ∈ input ( e ); this procedure is called the wavefront merging procedure . The main algorithm.
The transparent edges of S (cid:48) are processed in a rough time order. Thewavefronts W ( e ) of each transparent edge e are constructed at the time ˜ d ( s, e ) + | e | , where ˜ d ( s, e ) isthe minimum geodesic distance from s to the two endpoints of e . Define covertime ( e ) = ˜ d ( s, e ) + | e | .The value covertime ( e ) will be computed during the algorithm. Initially, for each edge e whose well-covering region U ( e ) contains s , W ( e ) and covertime ( e ) are computed directly (and set covertime ( e ) = ∞ for all other edges); we refer to this step as the initialization step , which will be elaborated below.The algorithm maintains a timer τ and processes the transparent edges e of S (cid:48) following the order of covertime ( e ).The main loop of the algorithm works as follows. As long as S (cid:48) has an unprocessed transparentedge, we do the following. First, among all unprocessed transparent edges, choose the one e withminimum covertime ( e ) and set τ = covertime ( e ). Second, call the wavefront merging procedure toconstruct the wavefronts W ( e ) from W ( f, e ) for all edges f ∈ input ( e ) satisfying covertime ( f ) We provide the details on the initialization step. Consider a transparentedge e whose well-covering region U ( e ) contains s . To compute W ( e ), we only consider straight-linepaths from s to e inside U ( e ). If U ( e ) does not have any island inside, then the points of e that can bereached from s by straight-line paths in U ( e ) form an interval of e and the interval can be computedby considering the tangents from s to the elementary chains on the boundary of U ( e ). Since U ( e ) has O (1) cells of S (cid:48) , ∂ U ( e ) has O (1) elementary chain fragments and thus computing the interval on e takes O (log n ) time. If U ( e ) has at least one island inside, then e may have multiple such intervals.As U ( e ) is the union of O (1) cells of S (cid:48) , the number of islands inside U ( e ) is O (1). Hence, e has O (1) such intervals, which can be computed in O (log n ) time altogether. These intervals form thewavefront W ( e ), i.e., each interval corresponds to a wavelet with generator s . From W ( e ), the value We include input ( e ) in output ( e ) mainly for the proof of Lemma 20, which relies on output ( e ) having a cycle enclosing e . overtime ( e ) can be immediately determined. More specifically, for each endpoint v of e , if one of thewavelets of W ( e ) covers v , then the segment sv is in U ( e ) and thus d ( s, v ) = | sv | ; otherwise, we set d ( s, v ) = ∞ . In this way, we find an upper bound for ˜ d ( s, e ) and set covertime ( e ) to this upper boundplus | e | . The algorithm correctness. At the time ˜ d ( s, e ) + | e | , all edges f ∈ input ( e ) whose wavefronts W ( f ) contribute a wavelet to W ( e ) must have already been processed. This is due to the property ofthe well-covering regions of S (cid:48) that d ( e, f ) ≥ · max {| e | , | f |} since f is on ∂ U ( e ). The proof is the sameas that of Lemma 4.2 [26], so we omit it. Note that this also implies that the geodesic distance d ( s, v )is correctly computed for each endpoint v of e . Therefore, after the algorithm is finished, geodesicdistances from s to endpoints of all transparent edges of S (cid:48) are correctly computed. Artificial wavelets. As in the HS algorithm, to limit the interaction between wavefronts fromdifferent sides of each transparent edge e , when a wavefront propagates across e , i.e., when W ( e ) iscomputed in the wavefront merging procedure, an artificial wavelet is generated at each endpoint v of e , with weight d ( s, v ). This is to eliminate a wavefront from one side of e if it arrives at e later thanthe wavefront from the other side of e . Topologically different paths. In the wavefront propagation procedure to pass W ( e ) to all edges g ∈ output ( e ), W ( e ) travels through the cells of the well-covering region U of either g or e . Since U may not be simply connected (e.g., the square-annulus), there may be multiple topologically differentshortest paths between e and g inside U ; the number of such paths is O (1) as U is the union of O (1) cellsof S (cid:48) . We propagate W ( e ) in multiple components of U , each corresponding to a topologically differentshortest path and defined by the particular sequence of transparent edges it crosses. These wavefrontsare later combined in the wavefront merging step at g . This also happens in the initialization step, asdiscussed above. Claiming a point. During the wavefront merging procedure at e , we have a set of wavefronts W ( f, e )that reach e from the same side for edges f ∈ input ( e ). We say that a wavefront W ( f, e ) claims a point p ∈ e if W ( f, e ) reaches p before any other wavefront. Further, for each wavefront W ( f, e ), the set ofpoints on e claimed by it forms an interval (the proof is similar to that of Lemma 4.4 [26]). Similarly,a wavelet of a wavefront claims a point p of e if the wavelet reaches p before any other wavelet of thewavefront, and the set of points on e claimed by any wavelet of the wavefront also forms an interval.For convenience, if a wavelet claims a point, we also say that the generator of the wavelet claims thepoint.Before presenting the details of the wavefront merging procedure and the wavefront propagationprocedure in the next two subsections, we first discuss the data structure (for representing elementarychains, generators, and wavefronts) and a monotonicity property of bisectors, which will be used laterin our algorithm. We use an array to represent each elementary chain. Then, a generator α = ( A, a ) can be representedby recording the indices of the two end vertices of its underlying chain A . In this way, a generatortakes O (1) additional space to record and binary search on A can be supported in O (log n ) time (e.g.,given a point p , find the tangent from p to A ). In addition, we maintain the lengths of the edges ineach elementary chain so that given any two vertices of the chain, the length of the sub-chain between15he two vertices can be obtained in O (1) time. All these preprocessing can be done in O ( n ) time andspace for all elementary chains.For a wavefront W ( e ) of one side of e , it is a list of generators α , α , . . . ordered by the intervalsof e claimed by these generators. Note that these generators are in the same side of e . Formally, wesay that a generator α is in one side of e if the initial vertex of α lies in that side of the supporting lineof e . We maintain these generators by a balanced binary search tree so that the following operationscan be supported in logarithmic time each. Let W be a wavefront with generators α , α , . . . , α k . Insert Insert a generator α to W . In our algorithm, α is inserted either in the front of α or after α k . Delete Delete a generator α i from W , for any 1 ≤ i ≤ k . Split Split W into two sublists at some generator α i so that the first i generators form a wavefrontand the rest form another wavefront. Concatenate Concatenate W with another list W (cid:48) of generators so that all generators of W are inthe front of those of W (cid:48) in the new list.We will show later that each wavefront involved in our algorithm has O ( h ) generators. Therefore,each of the above operation can be performed in O (log h ) time. We make the tree fully persistent bypath-copying [16]. In this way, each operation on the tree will cost O (log h ) additional space but theolder version of the tree will be kept intact (and operations on the old tree can still be supported). Suppose α = ( A , a ) and α = ( A , a ) are two generators on the same side of an axis-parallel line (cid:96) (i.e., a and a are on the same side of (cid:96) ; it is possible that (cid:96) intersects A and A ).Then, the bisector B ( α , α ) intersects (cid:96) at no more than one point. Proof: Consider a generator α = ( A, a ). Recall that a is an endpoint of A . Let a (cid:48) be the otherendpoint of A . We define R ( α ) as the set of points in the plane that are reachable from α (withoutconsidering any other obstacles), and we call it the reachable region of α . The reachable region R ( α )is bounded by A , a ray ρ ( a ) with origin a , and another ray ρ ( a (cid:48) ) with origin a (cid:48) (e.g., see Fig. 12).Specifically, ρ ( a ) is along the direction from z to a , where z is the anchor of a in the shortest pathfrom s to a ; ρ ( a (cid:48) ) is the ray from the adjacent obstacle vertex of a (cid:48) on A (cid:48) to a (cid:48) .We claim that the intersection (cid:96) ∩ R ( α ) must be an interval of (cid:96) . Indeed, since A is xy -monotone,without loss of generality, we assume that a (cid:48) is to the northwest of a . Let ρ (cid:48) ( a ) be the verticallydownward ray from a (e.g., see Fig. 13). Observe that the concatenation of ρ ( a (cid:48) ), A , and ρ (cid:48) ( a ) is xy -monotone; let R (cid:48) be the region of the plane to the right of their concatenation. Since the boundaryof R (cid:48) is xy -monotone and (cid:96) is axis-parallel, R (cid:48) ∩ (cid:96) must be an interval of (cid:96) . Notice that ρ ( a ) must bein R (cid:48) and thus it partitions R (cid:48) into two subregions, one of which is R ( α ). Since (cid:96) intersects ρ ( a ) atmost once and R (cid:48) ∩ (cid:96) is an interval, R ( α ) ∩ (cid:96) must also be an interval.According to the above claim, both (cid:96) ∩ R ( α ) and (cid:96) ∩ R ( α ) are intervals of (cid:96) . Let I denote thecommon intersection of the two intervals. Since B ( α , α ) ⊆ R ( α ) and B ( α , α ) ⊆ R ( α ), we obtainthat B ( α , α ) ∩ (cid:96) ⊆ I .Assume to the contrary that B ( α , α ) ∩ (cid:96) contains two points, say, q and q . Then, both pointsare in I . Let v and u be the tangents points of q on A and A , respectively. Let v and u bethe tangents points of q on A and A , respectively (e.g., see Fig. 14). As I = (cid:96) ∩ R ( α ) ∩ R ( α )and I contains both q and q , if we move a point q from q to q on (cid:96) , the tangent from q to A willcontinuously change from q v to q v and the tangent from q to A will continuously change from q u to q u . Since A and A are in the same side of (cid:96) , either q u intersects q v in their interiors16 aa (cid:48) ρ ( a ) ρ ( a (cid:48) ) zs R ( α ) Figure 12: Illustrating the reachable region R ( α ) of agenerator α = ( A, a ). Aaa (cid:48) ρ (cid:48) ( a ) ρ ( a (cid:48) ) R (cid:48) Figure 13: Illustrating the ray ρ (cid:48) ( a ) and the region R (cid:48) . A A Iq q a = v v u u p a Figure 14: Illustrating the proof of Lemma 5. or q v intersects q u in their interiors; without loss of generality, we assume that it is the formercase. Let p be the intersection of q u and q v (e.g., see Fig. 14). Since q ∈ B ( α , α ), points of q u other than q have only one predecessor, which is a . As p ∈ q u and p (cid:54) = q , p has only onepredecessor a . Similarly, since q ∈ B ( α , α ) and p ∈ q v , a is also p ’s predecessor. We thus obtaina contradiction. (cid:3) Corollary 3 Suppose α = ( A , a ) and α = ( A , a ) are two generators both below a horizontal line (cid:96) . Then, the portion of the bisector B ( α , α ) above (cid:96) is y -monotone. Proof: For any horizontal line (cid:96) (cid:48) above (cid:96) , since both generators are below (cid:96) , they are also below (cid:96) (cid:48) .By Lemma 5, B ( α , α ) ∩ (cid:96) (cid:48) is either empty or a single point. The corollary thus follows. (cid:3) In this section, we present the details of the wavefront merging procedure. Given all contributingwavefronts W ( f, e ) of f ∈ input ( e ) for W ( e ), the goal of the procedure is to compute W ( e ). The algo-rithm follows the high-level scheme of the HS algorithm (i.e., Lemma 4.6 [26]) but the implementationdetails are quite different.We only consider the wavefronts W ( f, e ) and W ( e ) for one side of e since the algorithm for theother side is analogous. Without loss of generality, we assume that e is horizontal and all wavefronts W ( f, e ) are coming from below e . We describe the algorithm for computing the interval of e claimedby W ( f, e ) if only one other wavefront W ( f (cid:48) , e ) is present. The common intersection of these intervalsof all such f (cid:48) is the interval of e claimed by W ( f, e ). Since | input ( e ) | = O (1), the number of such f (cid:48) is O (1).We first determine whether the claim of W ( f, e ) is to the left or right of that of W ( f (cid:48) , e ). To thisend, depending on whether both W ( f, e ) and W ( f (cid:48) , e ) reach the left endpoint v of e , there are twocases. Note that the intervals of e claimed by W ( f, e ) and W ( f (cid:48) , e ) are available from W ( f, e ) and W ( f (cid:48) , e ); let I f and I f (cid:48) denote these two intervals, respectively.17 If both I f and I f (cid:48) contain v , i.e., both W ( f, e ) and W ( f (cid:48) , e ) reach v , then we compute the(weighted) distances from v to the two wavefronts. This can be done as follows. Since v ∈ I f , v must be reached by the leftmost generator α of W ( f, e ). We compute the distance d ( α, v ) bycomputing the tangent from v to α in O (log n ) time. Similarly, we compute d ( α (cid:48) , v ), where α (cid:48) is leftmost generator of W ( f (cid:48) , e ). If d ( α, v ) ≤ d ( α (cid:48) , v ), then the claim of W ( f, e ) is to the left ofthat of W ( f (cid:48) , e ); otherwise, the claim of W ( f, e ) is to the right of that of W ( f (cid:48) , e ). • If not both I f and I f (cid:48) contain v , then the order of the left endpoints of I f and I f (cid:48) will give theanswer.Without loss of generality, we assume that the claim of W ( f, e ) is to the left of that of W ( f (cid:48) , e ).We next compute the interval I of e claimed by W ( f, e ) with respect to W ( f (cid:48) , e ). Note that the leftendpoint of I is the left endpoint of I f . Hence, it remains to find the right endpoint of I , as follows.Let (cid:96) e be the supporting line of e . Let α be the rightmost generator of W ( f, e ) and let α (cid:48) be theleftmost generator of W ( f (cid:48) , e ). Let q be the left endpoint of the interval on e claimed by α in W ( f, e ),i.e., q is the intersection of the bisector B ( α , α ) and e , where α is the left neighboring generator of α in W ( f, e ). Similarly, q be the right endpoint of the interval on e claimed by α (cid:48) in W ( f (cid:48) , e ), i.e., q is the intersection of e and the bisector B ( α (cid:48) , α (cid:48) ), where α (cid:48) is the right neighboring generator of α (cid:48) in W ( f (cid:48) , e ). Let q be the intersection of the bisector B ( α, α (cid:48) ) and e . We assume that the three points q i , i = 0 , , O (log n )time by a bisector-line intersection operation given in Lemma 7. If q is between q and q , then q isthe right endpoint of I and we can stop the algorithm. If q is to the left of q , then we delete α from W ( f, e ). If q is to the right of q , then we delete α (cid:48) from W ( f (cid:48) , e ). In either case, we continue thesame algorithm by redefining α or α (cid:48) (and recomputing q i for i = 0 , , O ((1 + k ) log n ) time, where k is the number of generators thatare deleted. We apply the algorithm on f and other f (cid:48) in input ( e ) to compute the correspondingintervals for f . The common intersection of all these intervals is the interval of e claimed by W ( f, e ).We do so for each f ∈ input ( e ), after which W ( e ) is obtained. Since the size of input ( e ) is O (1), weobtain the following lemma. Lemma 6 Given all contributing wavefronts W ( f, e ) of edges f ∈ input ( e ) for W ( e ) , we can comptethe interval of e claimed by each W ( f, e ) and thus construct W ( e ) in O ((1 + k ) log n ) time, where k isthe total number of generators in all wavefronts W ( f, e ) that are absent from W ( e ) . Lemma 7 (Bisector-line intersection operation) Each bisector-line intersection operation canbe performed in O (log n ) time. Proof: Given two generators α = ( A , a ) and α = ( A , a ) below a horizontal line (cid:96) , the goalof the operation is to compute the intersection between the bisector B ( α , α ) and (cid:96) . By Lemma 5, B ( α , α ) ∩ (cid:96) is either empty or a single point.We apply a prune-and-search technique from Kirkpatrick and Snoeyink [32]. To avoid the lengthybackground explanation, we follow the notation in [32] without definition. We will rely on Theorem 3.6in [32]. To do so, we need to define a decreasing function f and an increasing function g .We first compute the intersection of (cid:96) and the reachable region R ( α ) of α , as defined in the proofof Lemma 5. This can be easily done by computing the intersections between (cid:96) and the boundary of R ( α ) in O (log n ) time. Let I be their intersection, which is an interval of (cid:96) as proved in Lemma 5.Similarly, we compute the intersection I of (cid:96) and the reachable region R ( α ) of α . Let I = I ∩ I . Foreach endpoint of I , we compute its tangent point on A (along its designated direction) to determinethe portion of A whose tangent rays intersect I and let A denote the portion. Similarly, we determinethe portion B of A whose tangent rays intersect I . All above takes O (log n ) time.18 B(cid:96) a f ( a ) B ( α , α ) A B(cid:96) B ( α , α ) bg ( b ) Figure 15: Illustrating the definitions of f ( a ) and g ( b ). We parameterize over [0 , 1] each of the two convex chains A and B in clockwise order, i.e., eachvalue of [0 , 1] corresponds to a slope of a tangent at a point on A (resp., B ). For each point a of A ,we define f ( a ) to be the parameter of the point b ∈ B such that the tangent ray of A at a along thedesignated direction and the tangent ray of B at b along the designated direction intersect at a pointon the bisector B ( α , α ) (e.g., see Fig. 15 left); if the tangent ray at a does not intersect B ( α , α ),then define f ( a ) = 1. For each point b of B , we define g ( b ) to be the parameter of the point a ∈ A suchthat the tangent ray of A at a and the tangent ray of B at b intersect at a point on the line (cid:96) (e.g.,see Fig. 15 right). One can verify that f is a continuous decreasing function while g is a continuousincreasing function (the tangent at an obstacle vertex of A and B is not unique but the issue canbe handled [32]). The fixed-point of the composition of the two functions g · f corresponds to theintersection of (cid:96) and B ( α , α ), which can be computed by applying the prune-and-search algorithmof Theorem 3.6 [32].As both chains A and B are represented by arrays (of size O ( n )), to show that the algorithm canbe implemented in O (log n ) time, it suffices to show that given any a ∈ A and any b ∈ B , we candetermine whether f ( a ) ≥ b in O (1) time and determine whether g ( b ) ≥ a in O (1) time.To determine whether f ( a ) ≥ b , we do the following. We first find the intersection q of the tangentray of A at a and the tangent ray of B at b . Then, f ( a ) ≥ b if and only if d ( α , q ) ≤ d ( α , q ). Noticethat d ( α , q ) = w ( a ) + | (cid:100) a a | + | aq | , where | (cid:100) a a | is the length of the portion of A between a and a .Recall that we have a data structure on each elementary chain C such that given any two vertices on C , the length of the sub-chain of C between the two vertices can be computed in O (1) time. Usingthe data structure, | (cid:100) a a | can be computed in constant time. Since w ( a ) is already known, d ( α , q )can be computed in constant time. So is d ( α , q ). In the case where the two tangent rays do notintersect, either the tangent ray of A at a intersects the backward extension of the tangent ray of B at b or the tangent ray of B at b intersects the backward extension of the tangent ray of A at a . Inthe former case f ( a ) ≤ b holds while in the latter case f ( a ) ≥ b holds. Hence, whether f ( a ) ≥ b canbe determined in constant time.To determine whether g ( b ) ≥ a , we do the following. Find the intersection p a between (cid:96) and thetangent ray of A at a and the intersection p b between (cid:96) and the tangent ray of B at b . If p a is tothe left of p b , then g ( b ) ≥ a ; otherwise g ( b ) ≤ a . Note that by the definition of A , the tangent ray atany point of A intersects I ; the same is true for B . Hence, whether g ( b ) ≥ a can be determined inconstant time.The above algorithm returns a point q in O (log n ) time. If the intersection of (cid:96) and B ( α , α )exists, then q is the intersection. Because we do not know whether the intersection exists, we finallyvalidate q by computing d ( α , q ) and d ( α , q ) in O (log n ) time as well as checking whether q ∈ (cid:96) . Thepoint q is valid if and only if d ( α , q ) = d ( α , q ) and q ∈ (cid:96) . (cid:3) .5 The wavefront propagation procedure In this section, we discuss the wavefront propagation procedure, which is to compute the wavefront W ( e, g ) for all transparent edges g ∈ output ( e ) based on W ( e ). Consider a transparent edge g ∈ output ( e ). The wavefront W ( e, g ) refers to the portion of W ( e ) that arrives at g through the well-covering region U ( g ) of g if e ∈ input ( g ) and through U ( e ) otherwise (in the latter case g ∈ input ( e )).We will need to handle the bisector events, i.e., the intersections between bisectors and the intersectionsbetween bisectors and obstacle edges. The HS algorithm processes the bisector events in temporal order, i.e., in order of the simulation time τ . The HSY algorithm instead proposes a simpler approachthat processes the events in spatial order , i.e., in order of their geometric locations. We will adapt theHSY’s spacial-order method.Recall that each wavefront W ( e ) is represented by a list of generators, which are maintained inthe leaves of a fully-persistent balanced binary search tree T ( e ). We further assign each generator a“next bisector event”, which is the intersection of its two bounding bisectors (it is set to null if the twobisectors do not intersect). More specifically, for each bisector α , we assign it the intersection of thetwo bisectors B ( α l , α ) and B ( α, α r ), where α l and α r are α ’s left and right neighboring generators in W ( e ), respectively; we store the intersection at the leaf for α . Our algorithm maintains a variant thatthe next bisector event for each generator in W ( e ) has already been computed and stored in T ( e ).We further endow the tree T ( e ) with additional node-fields so that each internal node stores a valuethat is equal to the minimum (resp., maximum) x -coordinate (resp., y -coordinate) among all bisectorevents stored at the leaves of the subtree rooted at the node. Using these extra values, we can findfrom a query range of generators the generator whose bisector event has the minimum/maximum x -or y -coordinate in logarithmic time.The propagation from W ( e ) to g through U is done cell by cell, where U is either U ( e ) or U ( g ). Westart propagating W ( e ) to the adjacent cell c of e in U to compute the wavefront passing through alledges of c . Then by using the computed wavefronts on the edges of c , we recursively run the algorithmon cells of U adjacent to c . As U has O (1) cells, the propagation passes through O (1) cells. Hence,the essential ingredient of the algorithm is to propagate a single wavefront, say, W ( e ), across a singlecell c with e on its boundary. Depending on whether c is an empty rectangle, there are two cases. c is an empty rectangle We first consider the case where c is an empty rectangle, i.e., there is no island inside c and c does notintersect any obstacle. Without loss of generality, we assume that e is an edge on the bottom side of c , and thus all generators of W ( e ) are below e . Our goal is to compute W ( e, g ), i.e., the generators of W ( e ) claiming g , for all other edges g of c . Our algorithm is similar to the HSY algorithm in the highlevel but the low level implementations are quite different. The main difference is that each bisectorin the HSY algorithm is of constant size while this is not the case in our problem. Due to this, it takesconstant time to compute the intersection of two bisectors in the HSY algorithm while in our problemthis costs O (log n ) time.The technical crux of the algorithm is to process the intersections in c among the bisectors ofgenerators of W ( e ). Since all generators of W ( e ) are below e , their bisectors in c are y -monotone byCorollary 3. This is a critical property our algorithm relies on. Due to the property, we only need tocompute W ( e, g ) for all edges g on the left, right, and top sides of c . Another helpful property is thatsince we propagate W ( e ) through e inside c , if a generator of α of W ( e ) claims a point q ∈ c , thenthe tangent from q to α must cross e ; we refer it as the propagation property . Due to this property,the points of c claimed by α must be to the right of the tangent ray from the left endpoint of e to α (the direction of the ray is the from the tangent point to the left endpoint of e ), as well as to the leftof the tangent ray from the right endpoint of e to α (the direction of the ray is the from the tangentpoint to the right endpoint of e ). We call the former ray the left bounding ray of α and the latter the20 ight bounding ray of α . As such, for the leftmost generator of W ( e ), we consider its left boundingray as its left bounding bisector; similarly, for the rightmost generator of W ( e ), we consider its rightbounding ray as its right bounding bisector.Starting from e , we use a horizontal line segment (cid:96) to sweep c upwards until its top side. At anymoment during the algorithm, the algorithm maintains a subset W ( (cid:96) ) of generators of W ( e ) for (cid:96) by abalanced binary search tree T ( (cid:96) ); initially W ( (cid:96) ) = W ( e ) and T ( (cid:96) ) = T ( e ). Let [ x , x ] × [ y , y ] denotethe coordinates of c . Using the extra fields on the nodes of the tree T ( (cid:96) ), we compute a maximal prefix W ( (cid:96) ) (resp., W ( (cid:96) )) of generators of W ( (cid:96) ) such that the bisector events assigned to all generators in ithave x -coordinates less than x (resp., larger than x ). Let W m ( (cid:96) ) be the remaining elements of W ( (cid:96) ).By definition, the first and last generators of W m ( (cid:96) ) have their bisector events with x -coordinates in[ x , x ]. As all bisectors are y -monotone in c , the lowest bisector intersection in c above (cid:96) must be the“next bisector event” b associated with a generator in W m ( (cid:96) ), which can be found in O (log n ) timeusing the tree T ( (cid:96) ). We advance (cid:96) to the y -coordinate of b by removing the generator α associatedwith the event b . Finally, we recompute the next bisector events for the two neighbors of α in W ( (cid:96) ).Specifically, let α l and α r be the left and right neighboring generators of α in W ( (cid:96) ), respectively. Weneed to compute the intersection of the two bounding bisectors of α l , and update the bisector eventof α l to this intersection. Similarly, we need to compute the intersection of the bounding bisectors of α r , and update the bisector event of α r to this intersection. Lemma 8 below shows that each of thesebisector intersections can be computed in O (log n ) time by a bisector-intersection operation, using thetentative prune-and-search technique of Kirkpatrick and Snoeyink [32]. Note that if α is the leftmostgenerator, then α r becomes the leftmost after α is deleted, in which case we compute the left boundingray of α r as its left bounding generator. If α is the rightmost generator, the process is symmetric. Lemma 8 (Bisector-bisector intersection operation) Each bisector-bisector intersection opera-tion can be performed in O (log n ) time. Proof: We are given a horizontal line (cid:96) and three generators α , α , and α below (cid:96) such that theyclaim points on (cid:96) in this order. The goal is to compute the intersection of the two bisectors B ( α , α )and B ( α , α ) above (cid:96) , or determine that such an intersection does not exist. Using the tentativeprune-and-search technique of Kirkpatrick and Snoeyink [32], we present an O (log n ) time algorithm.To avoid the lengthy background explanation, we follow the notation in [32] without definition.We will rely on Theorem 3.9 in [32]. To this end, we need to define three continuous and decreasingfunctions f , g , and h . We define them in a way similar to Theorem 4.10 in [32] for finding a pointequidistant to three convex polygons. Indeed, our problem may be considered as a weighted case oftheir problem because each point in the underlying chains of the generators has a weight that is equalto its weighted distance from its generator.Let A , B , and C be the underlying chains of α , α , and α , respectively.We parameterize over [0 , 1] each of the three convex chains A , B , and C from one end to the otherin clockwise order, i.e., each value of [0 , 1] corresponds to a slope of a tangent at a point on the chains.For each point a of A , we define f ( a ) to be the parameter of the point b ∈ B such that the tangentray of A at a (following the designated direction) and the tangent ray of B at b intersect at a pointon the bisector B ( α , α ) (e.g., see Fig. 15 left); if the tangent ray at a does not intersect B ( α , α ),then define f ( a ) = 1. In a similar manner, we define g ( b ) for b ∈ B with respect to C and define h ( c )for c ∈ C with respect to A . One can verify that all three functions are continuous and decreasing(the tangent at an obstacle vertex of the chains is not unique but the issue can be handled [32]).The fixed-point of the composition of the three functions h · g · f corresponds to the intersection of B ( α , α ) and B ( α , α ), which can be computed by applying the tentative prune-and-search algorithmof Theorem 3.9 [32].To guarantee that the algorithm can be implemented in O (log n ) time, since each of the chains A , B , and C is represented by an array, we need to show that given any a ∈ A and any b ∈ B , we21an determine whether f ( a ) ≥ b in O (1) time. This can be done in the same way as in the proof ofLemma 7. Similarly, given any b ∈ B and c ∈ C , we can determine whether g ( b ) ≥ c in O (1) time,and given any c ∈ C and a ∈ A , we can determine whether h ( c ) ≥ a in O (1) time. Therefore, applyingTheorem 3.9 [32] can compute the intersection of B ( α , α ) and B ( α , α ) in O (log n ) time.The above algorithm is based on the assumption that the intersection of the two bisectors exists.However, we do not know whether that is true or not. To determine it, we finally validate the solutionas follows. Let q be the point returned by the algorithm. We compute the distances d ( α i , q ) for i = 1 , , 3. The point q is a true bisector intersection if and only if the three distances are equal.Finally, we return q if and only if q is above (cid:96) . (cid:3) The algorithm finishes once (cid:96) is at the top side of c . At this moment, no bisector events of W ( (cid:96) )are in c . Finally, we run the following wavefront splitting step to split W ( (cid:96) ) to obtain W ( e, g ) for alledges g on the left, right, and top sides c . Our algorithm relies on the following observation. Let ζ bethe union of the left, top, and right sides of c . Lemma 9 The list of generators of W ( (cid:96) ) are exactly those in W ( (cid:96) ) claiming ζ in order. Proof: It suffices to show that during the sweeping algorithm whenever a bisector B ( α , α ) of twogenerators α and α intersects ζ , it will never intersect ∂c again. Let q be such an intersection. Let ζ l , ζ t , and ζ r be the left, top, and right sides of c , respectively.If q is on ζ t , then since both α and α are below e , they are also below ζ t . By Lemma 5, B ( α , α )will not intersect the supporting line of ζ t again and thus will not intersect ∂c again.If q is on ζ l , then we claim that both generators α and α are to the right of the supporting line (cid:96) l of ζ l . Indeed, since both generators claim q , the bounding rays (i.e., the left bounding ray of theleftmost generator of W ( (cid:96) ) and the right bounding ray of the rightmost generator of W ( (cid:96) ) during thesweeping algorithm) guarantee the propagation property: the tangents from q to the two generatorsmust cross e . Therefore, both generators must be to the right of (cid:96) l . By Lemma 5, B ( α , α ) will notintersect the supporting line of ζ l again and thus will not intersect ∂c again.If q is on ζ r , the analysis is similar to the above second case. (cid:3) In light of the above lemma, our wavefront splitting step algorithm for computing W ( e, g ) of alledges g ∈ ζ works as follows. Consider an edge g ∈ ζ . Without loss of generality, we assume that thepoints of ζ are clockwise around c so that we can talk about their relative order.Let p l and p r be the front and rear endpoints of g , respectively. Let α l and α r be the generators of W ( (cid:96) ) claiming p l and p r , respectively. Then all generators of W ( (cid:96) ) to the left of α l including α l formthe wavefront for all edges of ζ in the front of g ; all generators of W ( (cid:96) ) to the right of α r including α r form the wavefront for all edges of ζ after g ; all generators of W ( (cid:96) ) between α l and α r including α l and α r form W ( e, g ). Hence, once α l and α r are known, W ( e, g ) can be obtained by splitting W ( (cid:96) ) in O (log n ) time. It remains to compute α l and α r . Below, we only discuss how to compute the generator α l since α r can be computed analogously.Starting from the root v of T ( (cid:96) ), we determine the intersection q between B ( α , α ) and ζ , where α is the rightmost generator in the left subtree of v and α is the leftmost generator of the rightsubtree of v . If q is in the front of g on ζ , then we proceed to the right subtree of v ; otherwise, weproceed to the left subtree of v .It is easy to see that the runtime of the algorithm is bounded by O ( η · log n ) time, where η is thetime for computing q . In the HSY algorithm, each bisector is of constant size and an oracle is assumedto exist that can compute q in O (1) time. In our problem, since a bisector may not be of constantsize, it is not clear how to bound η by O (1). But η can be bounded by O (log n ) using the bisector-lineintersection operation in Lemma 7. Thus, α l can be computed in O (log n ) time. However, this is not22ufficient for our purpose, as this would lead to an overall O ( n + h log h ) time algorithm. We insteaduse the following binary search plus bisector tracing approach.During the wavefront expansion algorithm, for each pair of neighboring generators α = ( A, a ) and α (cid:48) = ( A (cid:48) , a (cid:48) ) in a wavefront (e.g., W ( e )), we maintain a special point z ( α, α (cid:48) ) on the bisector B ( α, α (cid:48) ).For example, in the above sweeping algorithm, whenever a generator α is deleted from W ( (cid:96) ) at abisector event b = B ( α l , α ) ∩ B ( α, α r ), its two neighbors α l and α r now become neighboring in W ( (cid:96) ).Then, we initialize z ( α l , α r ) to b (the tangent points from b to α l and α r are also associated with b ).During the algorithm, the point z ( α l , α r ) will move on B ( α, α (cid:48) ) further away from the two defininggenerators α and α (cid:48) and the movement will trace out the hyperbolic-arcs of the bisector. We call z ( α, α (cid:48) ) the tracing-point of B ( α, α (cid:48) ). Our algorithm maintains a variant that the tracing point of eachbisector of W ( (cid:96) ) is below the sweeping line (cid:96) (initially, the tracing point of each bisector of W ( e ) isbelow e ).With the help of the above z -points, we compute the generator α l as follows. Like the abovealgorithm, starting from the root v of T ( (cid:96) ), let α and α be the two generators as defined above. Tocompute the intersection q between B ( α , α ) and ζ , we trace out the bisector B ( α , α ) by moving itstracing-point z ( α , α ) upwards (each time trace out a hyperbolic-arc of B ( α , α )) until the currenttracing hyperbolic-arc intersects ζ at q . If q is in the front of e on ζ , then we proceed to the rightsubtree of v ; otherwise, we proceed to the left subtree of v .After W ( e, g ) is obtained, we compute W ( e, g (cid:48) ) for other edges g (cid:48) on ζ using the same algorithmas above. For the time analysis, observe that each bisector hyperbolic-arc will be traced out at mostonce in the wavefront splitting step for all edges of ζ because the tracing point of each bisector willnever move “backwards”.This finishes the algorithm for propagating W ( e ) through the cell c . Except the final wavefrontsplitting step, the algorithm runs in O ((1 + h c ) log n ) time, where h c is the number of bisector eventsin c . Because c has O (1) edges, the wavefront splitting step takes O (log n + n c ) time, where n c is thenumber of hyperbolic-arcs of bisectors that are traced out. c is not an empty rectangle We now discuss the case in which the cell c is not an empty rectangle. In this case, c has a squarehole inside or/and the boundary of c contains obstacle edges. Without loss of generality, we assumethat e is on the bottom side of c .If c contains a square hole, then we partition c into four subcells by cutting c with two lines parallelto e , each passing through an edge of the hole. If c has obstacle edges on its boundary, recall that theseobstacles edges belong to O (1) convex chains (each of which is a fragment of an elementary chain);we further partition c by additional edges parallel to e , so that each resulting subcell contains at mosttwo convex chains, one the left side and the other on the right side. Since ∂c has O (1) convex chains, O (1) additional edges are sufficient to partition c into O (1) subcells as above. Then, we propagate e through the subcells of c , one by one. In the following, we describe the algorithm for one such subcell.By slightly abusing the notation, we still use c to denote the subcell with e on its bottom side.Since ∂c has obstacle edges, the propagation algorithm becomes more complicated. As in the HSYalgorithm, comparing with the algorithm for the previous case, there are two new bisector events. • First, a bisector may intersect a convex chain (and thus intersect an obstacle). The HSY algo-rithm does not explicitly compute these bisector events because such an oracle is not assumedto exist. In our algorithm, however, because the obstacles in our problem are polygonal, we canexplicitly determine these events without any special assumption. This is also a reason that thehigh-level idea of our algorithm is simpler than the HSY algorithm. • Second, new generators may be created at the convex chains. We still sweep a horizontal line (cid:96) from e upwards. Let W ( (cid:96) ) be the current wavefront at some moment during the algorithm.23 ( α , α ) ceζ l q (cid:48) α α q (cid:48)(cid:48) Figure 16: Illustrating the creation of a new generator at q (cid:48) . Consider two neighboring generators α and α in W ( (cid:96) ) with α on the left of α . We use ζ l to denote the convex chain on the left side of c . Let q (cid:48) be the tangent point on ζ l of thecommon tangent between ζ l and α and let q (cid:48)(cid:48) be the tangent point on α (e.g., see Fig. 16).If d ( α , q (cid:48) ) < d ( α , q (cid:48) ), then a new generator α on ζ l with initial vertex q (cid:48) and weight equal to d ( α , q (cid:48) ) is created (designated counterclockwise direction) and inserted into W ( (cid:96) ) right before α . The bisector B ( α, α ) is the ray emanating from q (cid:48) and extending away from q (cid:48)(cid:48) . The regionto the left of the ray has α as its predecessor. When the sweeping line (cid:96) is at q (cid:48) , all wavelets in W ( (cid:96) ) to the left of α have already collided with ζ l and thus the first three generators of W ( (cid:96) )are α , α , and α .In what follows, we describe our sweeping algorithm to propagate W ( e ) through c . We begin withan easier case where only the left side of c is a convex chain, denoted by ζ l (and the right side is avertical transparent edge, denoted by ζ r ). We use ζ t to denote the top side of c , which is a transparentedge. As in the previous case, we sweep a line (cid:96) from e upwards until the top side ζ t . During thealgorithm, we maintain a list W ( (cid:96) ) of generators by a balanced binary search tree T ( (cid:96) ). Initially, W ( (cid:96) ) = W ( e ) and T ( (cid:96) ) = T ( e ).We compute the intersection q of the convex chain ζ l and the bisector B ( α , α ), for the leftmostbisectors of α and α of W ( (cid:96) ). We call it the bisector-chain intersection operation . The followinglemma shows that this operation can be performed in O (log n ) time. Lemma 10 (Bisector-chain intersection operation) Each bisector-chain intersection operationcan be performed in O (log n ) time. Proof: We are given a convex chain ζ l above a horizontal line (cid:96) and two generators α = ( A , a )and α = ( A , a ) below (cid:96) such that they claim points on (cid:96) in this order. The goal is to compute theintersection of B ( α , α ) and ζ l , or determine that they do not intersect. In the following, using thetentative prune-and-search technique of Kirkpatrick and Snoeyink [32], we present an O (log n ) timealgorithm.To avoid the lengthy background explanation, we follow the notation in [32] without definition.We will rely on Theorem 3.9 in [32]. To this end, we need to define three continuous and decreasingfunctions f , g , and h .Suppose q is the intersection of B ( α , α ) and ζ l . Let p and p be the tangent points from q to A and A , respectively. Then, qp (resp., qp ) does not intersect ζ l other than q . We determine theportion C of ζ l such that for each point p ∈ C , its tangent to A does not intersect ζ l other than p .Hence, q ∈ C . C can be determined by computing the common tangents between ζ l and A , whichcan be done in O (log n ) time [20, 40]. Also, we determine the portion B of A such that the tangentray at any point of B must intersect C . This can be done by computing the common tangents between C and A in O (log n ) time [20, 40]. Let A = A . 24 Cg ( b ) A Bb (cid:96) CcA Bh ( c )0 00 0 001 11 1 11 Figure 17: Illustrating the definitions of g ( b ) and h ( c ). We parameterize over [0 , 1] each of the three convex chains A , B , and C from one end to the otherin clockwise order, i.e., each value of [0 , 1] corresponds to a slope of a tangent at a point on the chains A and B , while each value of [0 , 1] corresponds to a point of C . For each point a of A , we define f ( a )to be the parameter of the point b ∈ B such that the tangent ray of A at a (following the designateddirection of α ) and the tangent ray of B at b intersect at a point on the bisector B ( α , α ) (e.g., seeFig. 15 left); if the tangent ray at a does not intersect B ( α , α ), then define f ( a ) = 1. For each point b of B , define g ( b ) to be the parameter of the point c ∈ C such that cb is tangent to B at b (e.g., seeFig. 17 left); note that by the definition of B , the tangent ray from any point of B must intersect C and thus such a point c ∈ C must exist. For each point c ∈ C , define h ( c ) to be the parameter of thepoint of a ∈ A such that ac is tangent to A at a (e.g., see Fig. 17 right); note that by the definitionof C , such a point a ∈ A must exist and ac does not intersect C other than c . One can verify that allthree functions are continuous and decreasing (the tangent at an obstacle vertex of the chains is notunique but the issue can be handled [32]). The fixed-point of the composition of the three functions h · g · f corresponds to the intersection q of B ( α , α ) and ζ l , which can be computed by applying thetentative prune-and-search algorithm of Theorem 3.9 [32].To make sure that the algorithm can be implemented in O (log n ) time, since each convex chainis part of an elementary chain and thus is represented by an array, it suffices to show the following:(1) given any a ∈ A and any b ∈ B , whether f ( a ) ≥ b can be determined in O (1) time; (2) given any b ∈ B and any c ∈ C , whether g ( b ) ≥ c can be determined in O (1) time; (3) given any c ∈ C and any a ∈ A , whether h ( c ) ≥ a can be determined in O (1) time. We prove them below.Indeed, for (1), it can be done in the same way as in the proof of Lemma 7. For (2), g ( b ) ≥ c if andonly if c is to the right of the tangent ray of B at b , which can be easily determined in O (1) time. For(3), h ( c ) ≥ a if and only c is to the right of the tangent ray of A at a , which can be easily determinedin O (1) time.Therefore, applying the tentative prune-and-search technique in Theorem 3.9 [32] can compute q in O (log n ) time.Note that the above algorithm is based on the assumption that the intersection of B ( α , α ) and ζ l exists. However, we do not know whether this is true or not. To determine that, we finally validate thesolution as follows. Let q be the point returned by the algorithm. We first determine whether q ∈ ζ l .If not, then the intersection does not exist. Otherwise, we further compute the two distances d ( α i , q )for i = 1 , O (log n ) time. If the distances are equal, then q is the true intersection; otherwise, theintersection does not exist. (cid:3) If the intersection q of ζ l and B ( α , α ) does not exist, then we compute the tangent between ζ l and α , which can be done in O (log n ) time [40]; let q (cid:48) be the tangent point at ζ l . Regardless whether q exists or not, we compute the lowest bisector intersection b in c above (cid:96) in the same way as in thealgorithm for the previous case where c is an empty rectangle. Depending on whether q exists or not,we proceed as follows. For any point p in the plane, let y ( p ) denote the y -coordinate of p .1. If q exists, then depending on whether y ( q ) ≤ y ( b ), there are two subcases. If y ( q ) ≤ y ( b ),25hen we process the bisector event q : remove α from W ( (cid:96) ) and then recompute q , q (cid:48) , and b .Otherwise, we process the bisector event at b in the same way as in the previous case and thenrecompute q , q (cid:48) , and b .2. If q does not exist, then depending on whether y ( b ) ≤ y ( q (cid:48) ), there are two subcases. If y ( b ) ≤ y ( q (cid:48) ), then we process the bisector event at b in the same way as before and then recompute q , q (cid:48) , and b . Otherwise, we insert a new generator α = ( A, q (cid:48) ) to W ( (cid:96) ) as the leftmost generator,where A is the fragment of the elementary chain containing ζ l from q (cid:48) counterclockwise to theend of the chain, and α is designated the counterclockwise direction of A and the weight of q (cid:48) is d ( α , q (cid:48) ); e.g., see Fig. 16. The ray from q (cid:48) in the direction from q (cid:48)(cid:48) to q (cid:48) is the bisector of α and α , where q (cid:48)(cid:48) is the tangent point on α of the common tangent between α and ζ l . We initializethe tracing-point z ( α, α ) of B ( α, α ) to q (cid:48) . Finally, we recompute q , q (cid:48) , and b .Once the sweep line (cid:96) reaches the top side ζ t of c , the algorithm stops. Finally, as in the previouscase, we run a wavefront splitting step. Because the left side ζ l consists of obstacle edges, we split W ( (cid:96) ) to compute W ( e, g ) for all transparent edges g on the top side ζ t and the right side ζ r of c . Thealgorithm is the same as the previous case.The above discusses the case that only the left side ζ l of c is a convex chain. For the general casewhere both the left and right sides of c are convex chains, the algorithm is similar. The difference isthat we have to compute a point p corresponding to q and a point p (cid:48) corresponding to q (cid:48) on the rightside ζ r of c . More specifically, p is the intersection of B ( α (cid:48) , α (cid:48) ) with ζ r , where α (cid:48) and α (cid:48) are the tworightmost generators of W . If p does not exist, then we compute the common tangent between ζ r and α (cid:48) , and p (cid:48) is the tangent point on ζ r .In the following, if q does not exist, we let y ( q ) be ∞ ; otherwise, q (cid:48) is not needed and we let y ( q (cid:48) )be ∞ . We apply the same convention to p and p (cid:48) . We define b as the bisector event in the same wayas before. In each step, we process the lowest point r of { q, q (cid:48) , b, p, p (cid:48) } . If r is q or p , we process it inthe same way as before for q . If r is q (cid:48) or p (cid:48) , we process it in the same way as before for q (cid:48) . If r is b ,we process it in the same way as before. After processing r , we recompute the five points. Each steptakes O (log n ) time. After the sweep line (cid:96) reaches the top side ζ t of c , W ( (cid:96) ) is W ( e, ζ t ) for the topside ζ t of c because both the left and right sides of c are obstacle edges. Finally, we run the wavefrontsplitting step on W ( (cid:96) ) to compute the wavefronts W ( e, g ) for all transparent edges g on ζ t .In summary, propagating W ( e ) through c takes O ((1 + h c ) · log n + n c ) time, where h c is the numberof bisector events (including both the bisector-bisector intersection events and the bisector-obstacleintersection events) and n c is the number of hyperbolic-arcs of bisectors that are traced out in thewavefront splitting step.We use the following lemma to summarize the algorithm for both cases (i.e., regardless whether c is an empty rectangle or not). Lemma 11 Suppose W ( e ) is a wavefront on a transparent edge of a cell c of the subdivision S (cid:48) .Then, W ( e ) can be propagated through c to all other transparent edges of c in O ((1 + h c ) log n + n c ) time, where h c is the number of bisector events (including both the bisector-bisector intersection eventsand bisector-obstacle intersection events) and n c is the number of hyperbolic-arcs of bisectors that aretraced out in the wavefront splitting step. In this section, we show that the running time of our wavefront expansion algorithm described aboveis bounded by O ( n + h log h ). For this and also for the purpose of constructing the shortest pathmap SPM (cid:48) ( s ) later in Section 3.8, as in the HS algorithm, we mark generators in the way that if agenerator α is involved in a true bisector event of SPM ( s ) (either a bisector-bisector intersection or a26isector-obstacle intersection) in a cell c of the subdivision S (cid:48) , then α is guaranteed to be in a set ofmarked generators for c (but a marked generator for c may not actually participate in a true bisectorevent in c ). The generator marking rules are presented below, which are consistent with those in theHS algorithm. Generator marking rules: 1. For any generator α = ( A, a ), if its initial vertex a lies in a cell c , then mark α in c .2. Let e be a transparent edge and let W ( e ) be a wavefront coming from some generator α ’s sideof e .(a) If α claims an endpoint b of e in W ( e ), or if it would do so except for an artificial wavefront,then mark α in all cells c incident to b .(b) If α ’s claim in W ( e ) is shortened or eliminated by an artificial wavelet, then mark α for c ,where c is the cell having e as an edge and on α ’s side of e .3. Let e and g be two transparent edges with g ∈ output ( e ). Mark a generator α of W ( e ) in bothcells having e as an edge if one of the following cases happens:(a) α claims an endpoint of g in W ( e, g );(b) α participates in a bisector event either during the wavefront propagation procedure forcomputing W ( e, g ) from W ( e ), or during the wavefront merging procedure for computing W ( g ). Note that α is also considered to participate in a bisector event if its claim on g isshortened or eliminated by an artificial wavelet.4. If a generator α of W ( e ) claims part of an obstacle edge during the wavefront propagationprocedure for propagating W ( e ) toward output ( e ) (this includes the case in which α participatesin a bisector-obstacle intersection event), then mark α in both cells having e as an edge.Note that each generator may be marked multiple times and each mark applies to one instance ofthe generator. Lemma 12 The total number of marked generators during the algorithm is at most O ( h ) . Because the proof of Lemma 12 is quite technical and lengthy, we devote the entire Section 3.7 toit. In the rest of this subsection, we use Lemma 12 to show that the running time of our wavefrontexpansion algorithm is bounded by O ( n + h log h ). Our goal is to prove the following lemma. Lemma 13 The wavefront expansion algorithm runs in O ( n + h log h ) time and space. First of all, by Lemma 1, constructing the conforming subdivision S (cid:48) can be done in O ( n + h log h )time and O ( n ) space.The wavefront expansion algorithm has two main subroutines: the wavefront merging procedureand the wavefront propagation procedure.The wavefront merging procedure is to construct W ( e ) based on W ( f, e ) for the edges f ∈ input ( e ).By Lemma 6, this step takes O ((1 + k ) log n ) time, where k is the total number of generators in allwavefronts W ( f, e ) that are absent from W ( e ). According to the algorithm, if a generator α is absentfrom W ( e ), it must be deleted at a bisector event. Thus, α must be marked by Rule 3(b). Dueto Lemma 12, the total sum of k in the entire algorithm is O ( h ). As such, the wavefront mergingprocedure in the entire algorithm takes O ( h log n ) time in total.27 ( α, α (cid:48) ) eα α (cid:48) uv oq Figure 18: Illustrating the definitions of u , v , q , and o . The wavefront propagation procedure is to compute W ( e, g ) by propagating W ( e ) to all edges g ∈ output ( e ) either through U ( g ) or U ( e ). By Lemma 11, the running time of the procedure is O ((1 + h c ) log n + n c ) time, where h c is the number of bisector events (including both the bisector-bisectorintersection events and bisector-obstacle intersection events) and n c is the number of hyperbolic-arcsof bisectors that are traced out in the wavefront splitting step. For each bisector-bisector intersectionevent, at least one involved generator is marked by Rule 3(b). For each bisector-obstacle intersectionevent, at least one involved generator is marked by Rule 4. Hence, by Lemma 12, the total sum of h c in the entire algorithm is O ( h ). In addition, Lemma 14 below shows that the total sum of n c inthe entire algorithm is O ( n ). Therefore, the wavefront propagation procedure in the entire algorithmtakes O ( n + h log n ) time in total. Lemma 14 The total number of traced hyperbolic-arcs of the bisectors in the entire algorithm is O ( n ) . Proof: First of all, notice that each extension bisector consists of a single hyperbolic-arc, which is astraight line. As each generator is marked by Rule 1, by Lemma 12, the total number of generatorscreated in the algorithm is O ( h ). Since each generator can define at most one extension bisector,the number of hyperbolic-arcs on extension bisectors is at most O ( h ). In the following, we focus onhyperbolic-arcs of non-extension bisectors. Instead of counting the number of traced hyperbolic-arcs,we will count the number of their endpoints.Consider a hyperbolic-arc endpoint o that is traced out. According to our algorithm, o is tracedout during the wavefront propagation procedure for propagating W ( e ) to compute W ( e, g ) for sometransparent edge e and g ∈ output ( e ). Suppose o belongs to a non-extension bisector B ( α, α (cid:48) ) of twogenerators α and α (cid:48) in W ( e ). Then, o must be defined by an obstacle edge uv of either α or α (cid:48) ,i.e., the ray ρ ( u, v ) emanating from v along the direction from u to v (which is consistent with thedesignated direction of the generator that contains uv ) hits B ( α, α (cid:48) ) at o (e.g., see Fig 18). Withoutloss of generality, we assume that uv belongs to α . In the following, we argue that uv can define O (1)hyperbolic-arc endpoints that are traced out during the entire algorithm (a hyperbolic-arc endpointdefined by uv is counted twice is it is traced out twice), which will prove Lemma 14 as there are O ( n )obstacle edges in total.We first discuss some properties. Since o ∈ B ( α, α (cid:48) ), both α and α (cid:48) claim o . As both generators arein W ( e ), the ray ρ ( u, v ) must cross e , say, at a point q (e.g., see Fig 18). Because o is traced out whenwe propagate W ( e ) to compute W ( e, g ), oq must be in either U ( g ) or U ( e ), i.e., oq ⊆ U ( e ) ∪ U ( g ). Wecall ( e, g ) the defining pair of o for uv . According to our algorithm, during the propagation from W ( e )to g , uv defines only one hyperbolic-arc endpoint, because it is uniquely determined by the wavefront W ( e ). As such, to prove that uv can define O (1) hyperbolic-arc endpoints that are traced out duringthe entire algorithm, it suffices to show that there are at most O (1) defining pairs for uv . Let Π denotethe set of all such defining pairs. We prove | Π | = O (1) below.For each pair ( e (cid:48) , g (cid:48) ) ∈ Π, according to the above discussion, uv and ( e (cid:48) , g (cid:48) ) define a hyperbolic-arcendpoint o (cid:48) such that o (cid:48) is on the ray ρ ( u, v ) and vo (cid:48) crosses e (cid:48) at a point q (cid:48) . Without loss of generality,28 ( α, α (cid:48) ) eα α (cid:48) uv oq o (cid:48) q (cid:48) e (cid:48) Figure 19: Illustrating the definitions of u , v , q , and o . ev oq o (cid:48) q (cid:48) e (cid:48) s b π W ( e ) ( s, v ) Figure 20: The dashed (red) path is π W ( e ) ( s, v ), whichcrosses e (cid:48) at b . we assume that ( e, g ) is a pair that minimizes the length | vq (cid:48) | . Let W ( e (cid:48) ) refer to the wavefront at e (cid:48) whose propagation to g (cid:48) traces out o (cid:48) .We partition Π into two subsets Π and Π , where Π consists of all pairs ( e (cid:48) , g (cid:48) ) of Π such that qo intersects e (cid:48) and Π = Π \ Π . Since oq ⊆ U ( e ) ∪ U ( g ), each well-covering region contains O (1) cells,and | output ( e (cid:48) ) | = O (1) for each transparent edge e (cid:48) , the size of Π is O (1).For Π , we further partition it into two subsets Π and Π , where Π consists of all pairs ( e (cid:48) , g (cid:48) )of Π such that e is in the well-covering region U ( e (cid:48) ) of e (cid:48) or e (cid:48) ∈ U ( e ) ∪ U ( g ), and Π = Π \ Π .Since e is in U ( e (cid:48) ) for a constant number of transparent edges e (cid:48) , each of U ( e ) and U ( g ) contains O (1)cells, and | output ( e (cid:48) ) | = O (1) for each e (cid:48) , it holds that | Π | = O (1). In the following, we argue thatΠ = ∅ , which will prove | Π | = O (1).Assume to the contrary that | Π | (cid:54) = ∅ and let ( e (cid:48) , g (cid:48) ) be a pair of Π . Since ( e (cid:48) , g (cid:48) ) ∈ Π , by thedefinition of Π , e (cid:48) does not intersect oq . By the definition of e , vq \ { q } does not intersect e (cid:48) . Recallthat q (cid:48) is the intersection of e (cid:48) and ρ ( u, v ). Therefore, the points v , q , o , q (cid:48) , and o (cid:48) appear on the ray ρ ( u, v ) in this order (e.g., see Fig 19). Further, since ( e (cid:48) , g (cid:48) ) ∈ Π , by the definition of Π , e is not in U ( e (cid:48) ) and e (cid:48) is not in U ( e ) ∪ U ( g ).Without loss of generality, we assume that e (cid:48) is horizontal and the wavefront W ( e (cid:48) ) is from below e (cid:48) (thus v is below e (cid:48) while o (cid:48) is above e (cid:48) ). Let F (cid:48) be the modified free space by replacing e (cid:48) with anopaque edge of open endpoints. Since the generator in W ( e (cid:48) ) that contains v claims q (cid:48) , π (cid:48) ( s, q (cid:48) ) is ashortest path from s to q (cid:48) in F (cid:48) , where π (cid:48) ( s, q (cid:48) ) is the path following the wavefront W ( e (cid:48) ). Since v is in the generator of W ( e (cid:48) ) that claims q (cid:48) , v is the anchor of q (cid:48) in π (cid:48) ( s, q (cid:48) ), i.e., the edge of the pathincident to q (cid:48) is vq (cid:48) . Let π (cid:48) ( s, v ) be the sub-path of π (cid:48) ( s, q (cid:48) ) between s (cid:48) and v . Then, π (cid:48) ( s, v ) is ashortest path from s to v in F (cid:48) .Recall that v is in the generator α of W ( e ). Let π W ( e ) ( s, v ) be the path from s to v following W ( e ). We claim that | π (cid:48) ( s, v ) | = | π W ( e ) ( s, v ) | . Assume to the contrary this is not true. Then, either | π (cid:48) ( s, v ) | < | π W ( e ) ( s, v ) | or | π W ( e ) ( s, v ) | < | π (cid:48) ( s, v ) | . • If | π (cid:48) ( s, v ) | < | π W ( e ) ( s, v ) | , then since π W ( e ) ( s, v ) is a shortest path from s to v in the modifiedfree space by considering e as an opaque edge of open endpoints, π (cid:48) ( s, v ) must cross the interiorof e . This means that π (cid:48) ( s, q (cid:48) ) = π (cid:48) ( s, v ) ∪ vq (cid:48) crosses e twice. Because π (cid:48) ( s, q (cid:48) ) is a shortest pathin F (cid:48) and e ∈ F (cid:48) , it cannot cross e twice, a contradiction. • If | π W ( e ) ( s, v ) | < | π (cid:48) ( s, v ) | , then since π (cid:48) ( s, v ) is a shortest path from s to v in F (cid:48) , the path π W ( e ) ( s, v ) cannot be in F (cid:48) and thus must cross the interior of e (cid:48) , say, at a point b (e.g., seeFig 19). Let π W ( e ) ( s, b ) be the sub-path of π W ( e ) ( s, v ) between s and b . Let π W ( e ) ( b, q (cid:48) ) be ashortest path from b to q (cid:48) along e (cid:48) by considering it as an opaque edge with open endpoints. Note29 voq o (cid:48) q (cid:48) e (cid:48) sv (cid:48) pq (cid:48)(cid:48) b Figure 21: The dashed (red) path is π W ( e ) ( s, v (cid:48) ), which crosses e (cid:48) at b . that if b and q (cid:48) are on different sides of e (cid:48) , then π W ( e ) ( b, q (cid:48) ) must be through an open endpointof e (cid:48) . It is not difficult to see that | π W ( e ) ( b, q (cid:48) ) | ≤ | e (cid:48) | . Let π W ( e ) ( s, q (cid:48) ) be the concatenation of π W ( e ) ( s, b ) and π W ( e ) ( b, q (cid:48) ). Hence, π W ( e ) ( s, q (cid:48) ) is a path in F (cid:48) . Since π (cid:48) ( s, q (cid:48) ) is a shortest pathfrom s to q (cid:48) in F (cid:48) , it holds that | π (cid:48) ( s, q (cid:48) ) | ≤ | π W ( e ) ( s, q (cid:48) ) | .Notice that | π W ( e ) ( s, q (cid:48) ) | = | π W ( e ) ( s, b ) | + | π W ( e ) ( b, q (cid:48) ) | ≤ | π W ( e ) ( s, v ) | + | e (cid:48) | < | π (cid:48) ( s, v ) | + | e (cid:48) | .On the other hand, | π (cid:48) ( s, q (cid:48) ) | = | π (cid:48) ( s, v ) | + | vq (cid:48) | ≥ | π (cid:48) ( s, v ) | + | qq (cid:48) | . We claim that | qq (cid:48) | ≥ | e (cid:48) | .Indeed, since e is outside U ( e (cid:48) ), q ∈ e , and q (cid:48) ∈ e (cid:48) , qq (cid:48) must cross ∂ U ( e (cid:48) ) at a point b (cid:48) . By theproperty of well-covering regions of S (cid:48) , | b (cid:48) q (cid:48) | ≥ | e (cid:48) | . Since | qq (cid:48) | ≥ | b (cid:48) q (cid:48) | , we obtain | qq (cid:48) | ≥ | e (cid:48) | .In light of the claim, we have | π (cid:48) ( s, q (cid:48) ) | ≥ | π (cid:48) ( s, v ) | + 2 | e (cid:48) | > | π (cid:48) ( s, v ) | + | e (cid:48) | > | π W ( e ) ( s, q (cid:48) ) | . Butthis incurs contradiction since | π (cid:48) ( s, q (cid:48) ) | ≤ | π W ( e ) ( s, q (cid:48) ) | .Therefore, | π (cid:48) ( s, v ) | = | π W ( e ) ( s, v ) | holds.As q (cid:48) (cid:54)∈ qo , we define p as a point on oq (cid:48) \ { o } infinitely close to o (e.g., see Fig. 21). Hence, p ∈ vq (cid:48) .Since π (cid:48) ( s, q (cid:48) ) = π (cid:48) ( s, v ) ∪ vq (cid:48) is a shortest path from s to q (cid:48) in F (cid:48) , π (cid:48) ( s, v ) ∪ vp is a shortest path from s to p to F (cid:48) .Recall that o is on the non-extension bisector B ( α, α (cid:48) ) of two generators α and α (cid:48) in W ( e ) and v is on α . Let v (cid:48) be the anchor of o in α (cid:48) (e.g., see Fig. 21). Since | π (cid:48) ( s, v ) | = | π W ( e ) ( s, v ) | , we have | π (cid:48) ( s, v ) | + | vo | = | π W ( e ) ( s, v ) | + | vo | = | π W ( e ) ( s, v (cid:48) ) | + | v (cid:48) o | , where π W ( e ) ( s, v (cid:48) ) is the path from s to v (cid:48) following W ( e ). By the definition of p and because B ( α, α (cid:48) ) is a non-extension bisector, it holdsthat | π W ( e ) ( s, v ) | + | vp | > | π W ( e ) ( s, v (cid:48) ) | + | v (cid:48) p | . As | π (cid:48) ( s, v ) | = | π W ( e ) ( s, v ) | , we have | π (cid:48) ( s, v ) | + | vp | > | π W ( e ) ( s, v (cid:48) ) | + | v (cid:48) p | . Since π (cid:48) ( s, v ) ∪ vp is a shortest path from s to p in F (cid:48) , π W ( e ) ( s, v (cid:48) ) ∪ v (cid:48) p cannotbe a path from s to p in F (cid:48) . Therefore, π W ( e ) ( s, v (cid:48) ) ∪ v (cid:48) p must intersect the interior of e (cid:48) .We claim that v (cid:48) p cannot intersect e (cid:48) . Indeed, since v (cid:48) p is covered by the wavefront W ( e ) when W ( e ) propagates to g in either U ( e ) or U ( g ), v (cid:48) p is in U ( e ) ∪ U ( g ). Since ( e (cid:48) , g (cid:48) ) ∈ Π and by thedefinition of Π , e (cid:48) is not in U ( e ) ∪ U ( g ). Therefore, v (cid:48) p cannot intersect e (cid:48) .The above claim implies that π W ( e ) ( s, v (cid:48) ) must intersect the interior of e (cid:48) at a point, say, b (e.g.,see Fig. 21). Let π W ( e ) ( s, b ) be the sub-path of π W ( e ) ( s, v (cid:48) ) between s and b . Let π W ( e ) ( b, q (cid:48) ) bea path from b to q (cid:48) along e (cid:48) by considering e (cid:48) as an opaque edge of open endpoints. As discussedabove, | π W ( e ) ( b, q (cid:48) ) | ≤ | e (cid:48) | . Hence, π W ( e ) ( s, b ) ∪ π W ( e ) ( b, q (cid:48) ) is a path in F (cid:48) , whose length is at most | π W ( e ) ( s, b ) | + | e (cid:48) | ≤ | π W ( e ) ( s, v (cid:48) ) | + | e (cid:48) | . On the other hand, | π W ( e ) ( s, v (cid:48) ) | + | v (cid:48) o | + | oq (cid:48) | = | π (cid:48) ( s, v ) | + | vo | + | oq (cid:48) | = | π (cid:48) ( s, v ) | + | vq (cid:48) | . Since ( e (cid:48) , g (cid:48) ) ∈ Π and by the definition of Π , e is outside U ( e (cid:48) ). Because q (cid:48) ∈ e (cid:48) and ov (cid:48) crosses e at a point q (cid:48)(cid:48) , by the property of well-covering regions of S (cid:48) , | q (cid:48)(cid:48) o | + | oq (cid:48) | ≥ | e (cid:48) | .30herefore, we obtain | π (cid:48) ( s, v ) | + | vq (cid:48) | = | π W ( e ) ( s, v (cid:48) ) | + | v (cid:48) o | + | oq (cid:48) |≥ | π W ( e ) ( s, v (cid:48) ) | + | q (cid:48)(cid:48) o | + | oq (cid:48) |≥ | π W ( e ) ( s, v (cid:48) ) | + 2 | e (cid:48) | ≥ | π W ( e ) ( s, b ) | + 2 | e (cid:48) | > | π W ( e ) ( s, b ) | + | e (cid:48) | ≥ | π W ( e ) ( s, b ) | + | π W ( e ) ( b, q (cid:48) ) | . This means that π W ( e ) ( s, b ) ∪ π W ( e ) ( b, q (cid:48) ) is a path from s to q (cid:48) in F (cid:48) that is shorter than the path π (cid:48) ( s, v ) ∪ vq (cid:48) . But this incurs a contradiction as π (cid:48) ( s, v ) ∪ vq (cid:48) is a shortest path from s to q (cid:48) in F (cid:48) .This completes the proof of the lemma. (cid:3) In summary, the total time of the wavefront expansion algorithm is O ( n + h log n ), which is O ( n + h log h ).For the space complexity of the algorithm, since each wavefront W ( e ) is maintained by a persistentbinary tree, each bisector event costs additional O (log n ) space. As discussed above, the total sum of k + h c in the entire algorithm, which is the total number of bisector events, is O ( h ). Hence, the spacecost by persistent binary trees is O ( h log n ). The space used by other parts of the algorithm is O ( n ).Hence, the space complexity of the algorithm is O ( n + h log n ), which is O ( n + h log h ). This provesLemma 13. In this section, we prove Lemma 12. The proof follows the same strategy in the high level as that inthe HS algorithm, although many details are different. We start with the following lemma. Lemma 15 Suppose Π is the set of pairs ( e, B ) of transparent edges e and bisectors B such that B crosses e in the wavefront W ( e ) of e during our wavefront expansion algorithm but the same crossingdoes not occur in SPM (cid:48) ( s ) . Then | Π | = O ( h ) . Proof: Let α = ( A , a ) and α = ( A , a ) be the two generators of B . Recall that the well-coveringregion U ( e ) of e is the union of O (1) cells of S (cid:48) and each cell has O (1) elementary chain fragments on itsboundary. Hence, U ( e ) has O (1) convex chains on its boundary. Recall also that W ( e ) is computed bymerging all contributing wavefronts W ( f, e ) for f ∈ input ( e ) in the wavefront merging procedure, and W ( f, e ) is computed by propagating W ( f ) to e through U ( e ) in the wavefront propagation procedure.We first argue that the number of pairs ( e, B ) of Π in the case where at least one of α and α is in U ( e ) is O ( h ). Indeed, according to our wavefront propagation procedure, if a subcell c ∈ U ( e )is an empty rectangle, then no new generators will be created when W ( f ) is propagating through c ;otherwise, although multiple generators may be created, at most two (one on each side of c ) are inthe wavefront existing c . As U ( e ) has O (1) cells and each cell may be partitioned into O (1) subcellsduring the wavefront propagation procedure, only O (1) generators of W ( e ) are inside U ( e ). Since eachgenerator of W ( e ) may define at most two bisectors in W ( e ) and the total number of transparentedges of S (cid:48) is O ( h ), the number of pairs ( e, B ) of Π such that at least one generator of B is in U ( e ) isat most O ( h ).In the following, we assume that both α and α are outside U ( e ), i.e., their initial vertices a and a are outside U ( e ). Let q be the intersection of B and e . Let π (cid:48) ( a , q ) denote the path from q to a following the tangent from q to A and then to a along A ; define π (cid:48) ( a , q ) similarly. Clearly, both α and α claim q in W ( e ). Since α is outside U ( e ), α must be in W ( f ) for some transparent edge f of ∂ U ( e ) so that α of W ( e ) is from W ( f , e ). Since α is outside U ( e ), π (cid:48) ( a , q ) must intersects f , say, at point q (e.g., see Fig. 22), and π (cid:48) ( q , q ) is inside U ( e ), where π (cid:48) ( q , q ) is the sub-path of31 f qq α a B ( α , α ) Figure 22: Illustrating q and q . ef qq α a B ( α , α ) α a q Figure 23: Illustrating the psuedo-triangle (cid:52) ( q, q , q )(the shaded region). π (cid:48) ( a , q ) between q and q . Also, f is processed (for the wavefront propagation procedure) earlierthan e because W ( f ) contributes to W ( e ).We claim that q is claimed by α in SPM (cid:48) ( s ). Assume to the contrary this is not true. Then, let π ( s, q ) be a shortest path from s to q in the free space F . • If π ( s, q ) does not intersect e , then π ( s, q ) is in the modified free space F (cid:48) by replacing e withan opaque edge of open endpoints. Since α claims q on e , π = π W ( e ) ( s, a ) ∪ π (cid:48) ( a , q ) is ashortest path from s to q in the modified free space F (cid:48) , where π W ( e ) ( s, a ) is the path from s to a following the wavefront W ( e ). Since q is in π and q is not claimed by α in SPM (cid:48) ( s ), if wereplace the portion between s and q in π by π ( s, q ), we obtain a shorter path from s to q in F (cid:48) than π . But this incurs a contradiction since π is a shortest path from s to q in F (cid:48) . • If π ( s, q ) intersects e , say, at a point b , then since f ∈ ∂ U ( e ) and thus d ( f , e ) ≥ · max {| f | , | e |} by the properties of the well-covering regions of S (cid:48) , we show below that e must be processedearlier than f , which incurs a contradiction because f is processed earlier than e .Indeed, covertime ( e ) = ˜ d ( s, e ) + | e | ≤ d ( s, b ) + · | e | + | e | = d ( s, b ) + · | e | . On the other hand, covertime ( f ) = ˜ d ( s, f ) + | f | ≥ d ( s, q ) − · | f | + | f | = d ( s, b ) + d ( b, q ) + · | f | . Since f ∈ ∂ U ( e ), b ∈ e , and q ∈ f , d ( b, q ) ≥ d ( f , e ) ≥ · | e | . Hence, we obtain covertime ( f ) >covertime ( e ) and thus e must be processed earlier than f .Therefore, q must be claimed by α in SPM ( s ).We define f , q , and π (cid:48) ( q, q ) analogously with respect to α . Similarly, we can show that q isclaimed by α in SPM (cid:48) ( s ).For each f ∈ ∂ U ( e ), the generators of W ( f ) that are also in W ( e ) may not form a single subsequenceof the generator list of W ( e ), but they must form a constant number of (maximal) subsequences.Indeed, since U ( e ) is the union of O (1) cells of S (cid:48) , the number of islands in U ( e ) is O (1). Thus, thenumber of topologically different paths from f to e in U ( e ) is O (1); each such path will introduce atmost one subsequence of generators of W ( e ) that are also from W ( f ). Therefore, the generators of W ( f ) that are also in W ( e ) form O (1) subsequences of the generator list of W ( e ).Since ∂ U ( e ) has O (1) transparent edges, the generator list of W ( e ) can be partitioned into O (1)subsequences each of which is from W ( f ) following topologically the same path for a single transparentedge f of ∂ U ( e ). Due to this property, we have the following observation: the number of pairs ofadjacent generators of W ( e ) that are from different subsequences is O (1).Due to the above observation, the number of pairs of edges f and f on ∂ U ( e ) in the case where f (cid:54) = f , or f = f but π (cid:48) ( q, q ) and π (cid:48) ( q, q ) are topologically different in U ( e ) is only O (1). Therefore,32he total number of pairs ( e, B ) of Π in that case is O ( h ). In the following, it suffices to consider the casewhere f = f , and π (cid:48) ( q, q ) and π (cid:48) ( q, q ) are topologically the same in U ( e ); for reference purpose, weuse Π (cid:48) to denote the subset of pairs ( e, B ) of Π with the above property. Below we prove | Π (cid:48) | = O ( h ).Let f = f = f . Since π (cid:48) ( q, q ) and π (cid:48) ( q, q ) are topologically the same in U ( e ), the region boundedby π (cid:48) ( q, q ), π (cid:48) ( q, q ), and q q must be in U ( e ); we call the above region a pseudo-triangle , for both π (cid:48) ( q, q ) and π (cid:48) ( q, q ) are convex chains, and we use (cid:52) ( q, q , q ) to denote it (e.g., see Fig. 23). Because( e, B ) is not an incident pair in SPM (cid:48) ( s ), the point q must be claimed by a different generator α in SPM (cid:48) ( s ), which must be from the side of e different than α and α . We have proved above that q isclaimed by α and q is claimed by α in SPM (cid:48) ( s ). Hence, there must be at least one bisector eventin SPM (cid:48) ( s ) that lies in the interior of the pseudo-triangle (cid:52) ( q, q , q ) . We charge the early demise of B to any one of these bisector events in (cid:52) ( p, q , q ).Note that the path π (cid:48) ( q, q ) (resp., π (cid:48) ( q, q )) is in a shortest path from s to q following the wavefront W ( e ) in the modified free space by replacing e with an opaque edge of open endpoints. Hence, if Π (cid:48) has other pairs ( e, B (cid:48) ) whose first element is e , then the corresponding paths π (cid:48) ( q, q ) and π (cid:48) ( q, q ) ofall these pairs are disjoint, and thus the corresponding pseudo-triangles (cid:52) ( q, q , q ) are also disjointin U ( e ). Hence, each bisector event of SPM (cid:48) ( s ) inside (cid:52) ( q, q , q ) is charged at most once for all pairsof Π (cid:48) that have e as the first element. Since each cell of S (cid:48) belongs to the well-covering region U ( e )of O (1) transparent edges e , each bisector event of SPM (cid:48) ( s ) is charged O (1) times for all pairs of Π (cid:48) .Because SPM (cid:48) ( s ) has O ( h ) bisector events by Corollary 2, the number of pairs of Π (cid:48) is at most O ( h ).This completes the proof of the lemma. (cid:3) Armed with Lemma 15, we prove the subsequent five lemmas, which together lead to Lemma 12. Lemma 16 The total number of marked generators by Rule 2(a) and Rule 3 is O ( h ) . Proof: Suppose α is a generator marked by Rule 2(a) for a transparent edge e . Then, α must be thefirst or last non-artificial generator of W ( e ). Hence, at most two generators of W ( e ) can be markedby Rule 2(a). As S (cid:48) has O ( h ) transparent edges, the total number of generators marked by Rule 2(a)is O ( h ).Suppose α is a generator marked by Rule 3(a) for the two transparent edges ( e, g ) with g ∈ output ( e ). Since α claims an endpoint of g , α must be the first or last generator of W ( e, g ). Hence, atmost two generators of W ( e, g ) can be marked by Rule 3(a). Since | output ( e ) | = O (1), at most O (1)generators can be marked for the pairs of transparent edges with e as the first edge. Since S (cid:48) has O ( h )transparent edges, the total number of generators marked by Rule 3(a) is O ( h ).Suppose α is a generator marked by Rule 3(b) for the two transparent edges ( e, g ) with g ∈ output ( e ). Note that α is a generator in W ( e ). We assume that α is not the first or last non-artificialgenerators of W ( e ); the first and last non-artificial generators of W ( e ), countered separately, sum to O ( h ) for all transparent edges e of S (cid:48) . Let α and α be the two neighboring generators of α in W ( e ).Thus, both α and α are non-artificial. We assume that the bisectors B ( α , α ) and B ( α, α ) intersect e in SPM (cid:48) ( s ); by Lemma 15, there are only O ( h ) bisector and transparent edge intersections thatappear in some approximate wavefront but not in SPM (cid:48) ( s ).Since α is marked by Rule 3(b), at least one of B ( α , α ) and B ( α, α ) fails to reach the boundaryof D ( e ) during the algorithm, where D ( e ) is the union of cells through which W ( e ) is propagated toall edges g (cid:48) ∈ output ( e ). Note that D ( e ) ⊆ U ( e ) ∪ (cid:83) g (cid:48) ∈ output ( e ) U ( g (cid:48) ), which contains O (1) cells of S (cid:48) as | output ( e ) | = O (1) and the well-covering region of each transparent edge is the union of O (1) cellsof S (cid:48) ; also, each cell of D ( e ) is within a constant number of cells of e . Without loss of generality, weassume that B ( α , α ) does not reach the boundary of D ( e ) and a bisector event on B ( α , α ) is detected This is also due to that for any point p ∈ q q , the shortest path π ( s, p ) cannot intersect e ; this can be proved by asimilar argument to the above for proving that q is claimed by α in SPM (cid:48) ( s ). This observation implies that α cannotclaim any points on f in SPM (cid:48) ( s ). B ( α , α ) also implies that an actual bisectorevent in SPM (cid:48) ( s ) happens to B ( α , α ) no later than the detected event; we charge the marking of α tothat bisector actual endpoint (i.e., a vertex) in SPM (cid:48) ( s ), and the endpoint is in D ( e ) because B ( α , α )intersects e in SPM (cid:48) ( s ). Since each cell of D ( e ) is within a constant number of cells of e , each cellof S (cid:48) belongs to D ( e (cid:48) ) for a constant number of transparent edges e (cid:48) of S (cid:48) . Therefore, each vertex of SPM (cid:48) ( s ) is charged O (1) times. Since SPM (cid:48) ( s ) has O ( h ) vertices by Lemma 3, the total number ofgenerators marked by Rule 3(b) is O ( h ). (cid:3) Lemma 17 The total number of marked generators by Rule 1 is O ( h ) . Proof: According to our algorithm, generators are only created during the wavefront propagationprocedure, which is to propagate W ( e ) to compute W ( e, g ) for all edges g ∈ output ( e ) through U ,where U is U ( e ) or U ( g ). The procedure has two cases depending on whether a cell c of U is an emptyrectangle. If yes, then no generators will be created in c . Otherwise, c may be partitioned into O (1)subcells, each of which may have a convex chain on its left side and/or its right side. Let c be such asubcell. Without loss of generality, we assume that the current wavefront W is propagating in c frombottom to top. We consider the generators created on the left side ζ l of c .According to our algorithm, a generator on ζ l is created by the leftmost generator of the currentwavefront W . As such, if the leftmost wavelet of W does not create a generator, then no generator on ζ l will be created. Otherwise, let α be a generator created on ζ l , which becomes the leftmost generatorof W at the point of creation. Let α (cid:48) be the right neighbor of α in W . According to our algorithm, if α does not involve in any bisector event in c , i.e., the bisector B ( α, α (cid:48) ) does not intersect any otherbisector in c during the propagation of W in c , then no more new generators will be created on ζ l .Otherwise, we charge the creation of α to the bisector event involving B ( α, α (cid:48) ); each such event canbe charged at most twice (one for a generator on ζ l and the other for a generator created on the rightside of c ). Recall that for each bisector event a generator is marked by Rule 3(b).In light of the above discussion, since each cell of S (cid:48) belongs to the well-covering region U ( e (cid:48) ) of O (1) transparent edges e (cid:48) , the total number of generators marked by Rule 1 is O ( h ) + O ( k ), where k is the total number of generators marked by Rule 3(b), which is O ( h ) by Lemma 16. The lemma thusfollows. (cid:3) Lemma 18 The total number of marked generators by Rule 4 is O ( h ) . Proof: Let α be the generator and e be the transparent edge specified in the rule statement. Accordingto Rule 4, α is marked because it claims part of an obstacle edge during the wavefront propagationprocedure for propagating W ( e ) to compute W ( e, g ) for edges g ∈ output ( e ). As in the proof ofLemma 16, define D ( e ) as the union of cells through which W ( e ) is propagated to all edges g ∈ output ( e ). Recall that D ( e ) contains O (1) cells of S (cid:48) , and each cell of D ( e ) is within a constantnumber of cells of e , which implies that each cell of S (cid:48) belongs to D ( e (cid:48) ) for a constant number oftransparent edges e (cid:48) of S (cid:48) .First of all, for the case where α is a generator created during the wavefront propagation procedure,the total number of such marked generators is no more than the total number of marked generatorsby Rule 1, which is O ( h ) by Lemma 17.Second, for the case where α is the first or last generator in W ( e ), the total number of suchgenerators is clearly O ( h ) since S (cid:48) has O ( h ) transparent edges e .Third, for the case where α claims a rectilinear extreme vertex, the total number of such markedgenerators is O ( h ). To see this, since D ( e ) has O (1) cells of S (cid:48) and each cell contains at most onerectilinear extreme vertex, D ( e ) contains O (1) rectilinear extreme vertices. Therefore, at most O (1)34enerators will be marked in D ( e ). Since each cell of S (cid:48) belongs to D ( e (cid:48) ) for O (1) transparent edges e (cid:48) of S (cid:48) and S (cid:48) contains O ( h ) transparent edges e , the total number of marked generators in this caseis O ( h ).Finally, any Rule 4 marked generator α that does not belong to any of the above cases has thefollowing property: α does not claim a rectilinear extreme vertex, and is not the first or last non-artificial generator in W ( e ), and is not a generator created in D ( e ). Let α (cid:48) be a neighbor of α in W ( e ).Due to the above property, α (cid:48) is a non-artificial generator. We can assume that B ( α (cid:48) , α ) intersects e in SPM (cid:48) ( s ); by Lemma 15, there are only O ( h ) bisector and transparent edge intersections that appearin some approximate wavefront but not in SPM (cid:48) ( s ). Hence, in SPM (cid:48) ( s ), B ( α (cid:48) , α ) terminates in D ( e ),either on an obstacle edge or in a bisector event before reaching an obstacle edge. In either case, wecharge the mark of α at e to this endpoint of B ( α (cid:48) , α ) in SPM (cid:48) ( s ), which is vertex of SPM (cid:48) ( s ). Becauseeach cell of S (cid:48) belongs to D ( e (cid:48) ) for a constant number of transparent edges e (cid:48) of S (cid:48) , each vertex of SPM (cid:48) ( s ) is charged at most O (1) times. As SPM (cid:48) ( s ) has O ( h ) vertices by Lemma 3, the total numberof marked generators by Rule 4 is O ( h ). (cid:3) Lemma 19 The total number of marked generators by Rule 2(b) is O ( h ) . Proof: Let α be the generator and e be the transparent edge specified in the rule statement.First of all, the total number of marked generators in the case where α is in U ( e ) is O ( h ), becausethe number of Rule 1 marked generators is O ( h ) by Lemma 17 and U ( e ) has O (1) cells. In thefollowing, we consider the other case where α is outside U ( e ).Since α is outside U ( e ), there is a transparent edge f on ∂ U ( e ) (i.e., f ∈ input ( e )), such that α is in W ( f ) and it is in W ( e ) because W ( f ) is propagated to e through U ( e ). Since α ’s claim isshortened or eliminated by an artificial wavefront, there must be a bisector event involving α duringthe computation of W ( f, e ) from W ( f ). Hence, α is also marked by Rule 3(b). We charge α to thisRule 3(b) mark for f . Since there are O ( h ) generators marked by Rule 3(b) by Lemma 16, the totalnumber of marked generators of Rule 2(b) is O ( h ). (cid:3) SPM ( s ) In this section, we compute the shortest path map SPM ( s ). To this end, we need to mark generatorsfollowing our rules during our wavefront merging and propagation procedures. As the total numberof all marked generators is O ( h ), marking these generators do not change the running time of ouralgorithm asymptotically. In the following, we show that SPM ( s ) can be constructed in O ( n + h log h )time with the help of the marked generators. In light of Lemma 4, we will focus on constructing SPM (cid:48) ( s ).We first show in Section 3.8.1 that the marked generators are sufficient in the sense that if agenerator participates in a true bisector event of SPM (cid:48) ( s ) in a cell c of S (cid:48) , then α must be marked in c . Then, we present the algorithm for constructing SPM (cid:48) ( s ), which consists of two main steps. Thefirst main step is to identify all vertices of SPM (cid:48) ( s ) and the second one is to compute the edges of SPM (cid:48) ( s ) and assemble them to obtain SPM (cid:48) ( s ). The first main step is described in Section 3.8.2 whilethe second one is discussed in Section 3.8.3. In this section we show that if a generator participates in a bisector event of SPM (cid:48) ( s ) in a cell c of S (cid:48) , then α must be marked in c . Our algorithm in the next section for computing SPM (cid:48) ( s ) relies onthis property. Our proof strategy again follows the high-level structure of the HS algorithm. We firstprove the following lemma. 35 eaαα α pqα (cid:48) R B ( α , α ) B ( α, α Figure 24: Illustrating the case where q ∈ g (cid:48) . gea αα α pqα (cid:48) R B ( α , α ) B ( α, α s Figure 25: Illustrating the case where q ∈ B ( α , α ). Lemma 20 Let α be a generator in an approximate wavefront W ( e ) for some transparent edge e .Suppose there is a point p ∈ e that is claimed by α in W ( e ) but not in SPM (cid:48) ( s ) (because the approximatewavefront from the other side of e reaches p first). Then, α is marked in the cell c on α ’s side of e . Proof: Assume to the contrary that α is not marked in c . Then, α must have two neighboringnon-artificial generators α and α in W ( e ) since otherwise Rule 2 would apply. Also, because alltransparent edges of ∂ U ( e ) are in output ( e ), the two bisectors B ( α , α ) and B ( α, α ) must exist thewell-covering region U ( e ) of e through the same transparent edge g ∈ ∂ U ( e ), since otherwise Rule 3or 4 would apply. Further, the region R bounded by B ( α , α ), B ( α, α ), e , and g must be a subset of U ( e ). Indeed, if R contains an island not in U ( e ), then α would claim an endpoint of a boundary edgeof the island, in which case Rule 3 would apply.Let α (cid:48) = ( A, a ) be the true predecessor of p in SPM (cid:48) ( s ). Without loss of generality, we assumethat e is horizontal and α is below e and thus α (cid:48) is above e . Let π (cid:48) ( a, p ) be the path from p along itstangent to A and then following A to a .We first consider the case where α (cid:48) is not in the interior of the well-covering region U ( e ), i.e., theinitial vertex a is not in the interior of U ( e ). Since R is a subset of of U ( e ) without non- U ( e ) islands, a is not in the interior of R and thus π (cid:48) ( a, p ) intersects ∂R , say, at a point q . In the following, we arguethat α has been involved in a bisector event detected by our algorithm, and thus marked in c . Let g (cid:48) be the portion of g between its intersections with B ( α , α ) and B ( α, α ). Depending on whether q ∈ g (cid:48) , there are two subcases. • If q ∈ g (cid:48) (e.g., see Fig. 24), then ˜ d ( s, g ) ≤ d ( s, a ) + | π (cid:48) ( a, q ) | + | g | / 2, where π (cid:48) ( a, q ) is the sub-pathof π (cid:48) ( a, p ) between a and q . Hence, the time τ g when artificial wavefronts originating from theendpoints of g cover g is no later than ˜ d ( s, g )+ | g | ≤ d ( s, a )+ | π (cid:48) ( a, q ) | +3 | g | / 2. Because g ∈ ∂ U ( e ), | pq | ≥ · max {| g | , | e |} . Hence, τ g ≤ d ( s, a ) + | π (cid:48) ( a, q ) | + 3 | g | / < d ( s, a ) + | π (cid:48) ( a, q ) | + | qp | = d ( s, a ) + | π (cid:48) ( a, p ) | < d W ( e ) ( s, p ), where d W ( e ) ( s, p ) is the length of the path from s to p followingthe wavefront W ( e ); the last inequality holds because α (cid:48) is the true predecessor of p while α isnot.On the other hand, the time τ e when the wavefront W ( e ) reaches an endpoint of e cannot beearlier than d W ( e ) ( s, p ) − | e | . Hence, the wavelet from α cannot reach g earlier than d W ( e ) ( s, p ) −| e | + d ( e, g ) ≥ d W ( e ) ( s, p ) −| e | +2 | e | ≥ d W ( e ) ( s, p )+ | e | , which is larger than τ g since τ g < d W ( e ) ( s, p )as proved above. Therefore, by the time the wavelet from α reaches g , the artificial waveletsof g have already covered g , eliminating the wavelet from α from reaching g . Thus, α must bemarked by Rule 3(b). • If q (cid:54)∈ g (cid:48) , then q is on one of the bisectors B ( α , α ) and B ( α, α ). Let π ( s, a ) be a shortest pathfrom s to a and let π ( s, p ) = π ( s, a ) ∪ π (cid:48) ( a, p ), which is a shortest path from s to p . If π ( s, p )36ntersects g , then we can use the same analysis as above to show that the wavelet from α willbe eliminated from reaching g by the artificial wavelets of g and thus α must be marked byRule 3(b). In the following, we assume that π ( s, p ) does not intersect g .Without of loss of generality, we assume that q ∈ B ( α , α ) (e.g., see Fig. 25). Since π (cid:48) ( a, p )is a subpath of the shortest path π ( s, p ), every point of π (cid:48) ( a, p ) has α (cid:48) as its predecessor. As q ∈ π (cid:48) ( a, p ), the predecessor of q is α (cid:48) . Let π ( s, q ) be the subpath of π ( s, p ) between s and q ,which is a shortest path. Since π ( s, p ) does not intersect g , π ( s, q ) does not intersect g . Hence, π ( s, q ) is a shortest path in the modified free space F (cid:48) by replacing g with an opaque edge ofopen endpoints. Therefore, during the wavefront propagation procedure for computing W ( e, g )or during the wavefront merging procedure for computing W ( g ), the point q , which is on thebisector B ( α , α ), must be claimed by α (cid:48) . Hence, a bisector event must be detected for B ( α , α )during the computation of W ( e, g ) or W ( g ). In either case, α must be marked by Rule 3(b).We next consider the case where α (cid:48) lies inside U ( e ), i.e., its initial vertex a is inside U ( e ). If a is notbetween the two bisectors B ( α , α ) and B ( α, α ), then π (cid:48) ( a, p ) must intersect one of the bisectors andthus we can use a similar argument as above second case to show that α must be marked. Otherwise, a is in the region R . This implies that not all points in R are claimed by α when W ( e ) is propagating to g , and therefore, a bisector event involving α must happen during the wavefront propagation procedureto propagate W ( e ) to compute W ( e, g ) and thus α is marked by Rule 3(b). (cid:3) With the help of the preceding lemma, we prove the following lemma. Lemma 21 If a generator α participates in a bisector event of SPM (cid:48) ( s ) in a cell c of S (cid:48) , then α mustbe marked in c . Proof: If a bisector has an endpoint on an obstacle edge of c , it either emanates from an obstaclevertex a on the edge (i.e., the bisector is an extension bisector) or defined by two generators that claimpart of the opaque edge. In the first case, a is the initial vertex of a new created generator in c andthus the generator is marked by Rule 1. In the second case, both generators are marked by Rule 4.Let α be a generator that participates in a bisector event of SPM (cid:48) ( s ) in a cell c of S (cid:48) . Assume tothe contrary that α is not marked for c . Then by Rule 2(a) there must be transparent edges e and f on the boundary of c such that both W ( e ) and W ( f ) contain the three consecutive generators α , α ,and α . Without loss of generality, we assume that W ( e ) enters c and W ( f ) leaves c . Let R be theregion of c bounded by e , f , and the two bisectors B ( α , α ) and B ( α, α ); e.g., see Fig. 26. The region R must be a subset of c . Indeed, if R contains an island not in c , then α would claim an endpoint ofa boundary edge of the island, in which case Rule 3 would apply.Since α participates in a bisector event of SPM (cid:48) ( s ) in c , at least one point p in R is not claimedby α in SPM (cid:48) ( s ). Let α (cid:48) = ( A, a ) be the true predecessor of p . Note that the initial vertex a must beoutside R since otherwise a bisector event involving α must happen when W ( e ) is propagating through c and thus α would be marked by Rule 3(b). Let π (cid:48) ( a, p ) be the path from p along its tangent to A and then following A to a . Since a is outside R and p is inside R , π (cid:48) ( a, p ) must intersect the boundaryof R .Because α , α , and α are three consecutive generators of W ( e ), no generator other than α onthe same side of e as α claims any point of R . Thus, π (cid:48) ( a, p ) does not cross e . Let b and b bethe intersections of f with B ( α , α ) and B ( α, α ), respectively. Then, if α claims both b and b in SPM (cid:48) ( s ), then π (cid:48) ( a, p ) cannot cross either bisector on ∂R , and thus it must cross b b , say, at a point q (e.g., see Fig. 26). Note that q satisfies the hypothesis of Lemma 20 and thus α is marked for c . It α does not claim either b or b , then that point satisfies the hypothesis of Lemma 20 and thus α ismarked for c . (cid:3) eaαα α pqα (cid:48) Rb b c B ( α , α ) B ( α, α Figure 26: Illustrating the proof of Lemma 21. SPM (cid:48) ( s )In this section, we compute the vertices of SPM (cid:48) ( s ). The next section will compute the edges of SPM (cid:48) ( s ) and assemble them to obtain SPM (cid:48) ( s ). Computing active regions. Because the approximate wavefronts are represented by persistenttrees, after the above wavefront expansion algorithm finishes, the approximate wavefronts W ( e ) forall transparent edges e of S (cid:48) are still available. Also, for each cell c and each transparent edge e of c , a set of marked generators in W ( e ) are known. Using these marked generators, we first break c into active and inactive regions such that no vertices of SPM (cid:48) ( s ) lie in the inactive regions. Notethat each unmarked generator does not participate in a bisector event in SPM (cid:48) ( s ) by Lemma 21.If α and α are neighboring generators on ∂c such that one of them is marked while the other isunmarked, their bisector belongs to SPM (cid:48) ( s ). Therefore, all such generators are disjoint and theytogether partition c into regions such that each region is claimed only by marked generators or onlyby unmarked generators; the former regions are active while the latter are inactive. We can computethese active regions as follows.Since the wavefronts W ( e ) for all transparent edges e of c are available, the list of all generatorsordered along the boundary of c is also available. For each pair of adjacent generators α and α , ifone of them is marked and the other is unmarked, then we explicitly compute the portion of theirbisector B ( α , α ) in c . We describe the details of this step below.First of all, observe that α and α must be from W ( e ) for the same transparent edge e of c , sinceotherwise each generator must claim an endpoint of their own transparent edge thus must have beenmarked by Rule 2(a). Without loss of generality, we assume that e is horizontal and both α and α arebelow e . Note that we cannot afford computing the entire bisector B ( α , α ) and then determiningits portion in c because the running time would be proportional to the size of B ( α , α ), i.e., thenumber of hyperbolic-arcs of B ( α , α ). Instead, we first compute the intersection b of B ( α , α ) and e , which can be done in O (log n ) time by the bisection-line intersection operation in Lemma 7. Notethat b is one endpoint of the portion of B ( α , α ) in c . To determine the other endpoint b (cid:48) , we do thefollowing. Our goal is to compute b (cid:48) in O ( | B c | + log n ) time, with B c = B ( α , α ) ∩ c . To this end,we could trace B c in c from b , and for each hyperbolic-arc of B c , we determine whether it intersects ∂c . Recall that ∂c consists of O (1) transparent edges and at most O (1) convex chains. To achievethe desired runtime, we need to determine whether each hyperbolic-arc of B c intersects ∂c in O (1)time. However, it is not clear to us whether this is possible as the size of each convex chain maynot be of O (1) size. To circumvent the issue, we use the following strategy. Before tracing B c , wefirst compute the intersection between B ( α , α ) and each convex chain on ∂c , which can be done in38 (log n ) time by the bisector-chain intersection operation in Lemma 10. Among all intersections, let b (cid:48)(cid:48) be the closest one to α . Then, we start to trace B c from b , for each hyperbolic-arc e (cid:48) , we determinewhether e (cid:48) intersects each of the transparent edges of ∂c . If not, we further check whether e (cid:48) contains b (cid:48)(cid:48) . If b (cid:48)(cid:48) (cid:54)∈ e (cid:48) , then e (cid:48) is in c and we continue the tracing; otherwise, b (cid:48) is b (cid:48)(cid:48) and we stop the algorithm.If e (cid:48) intersects at least one of the transparent edges of ∂c , then among all such intersections as well as b (cid:48)(cid:48) , b (cid:48) is the one closest to α , which can be determined in O (log n ) time by computing their tangentsto α . In this way, computing the portion of B ( α , α ) in c can be done in O (log n + n c ( α , α )) time,where n c ( α , α ) is the number of hyperbolic-arcs of B ( α , α ) in c .In this way, all active regions of c can be computed in O ( h c log n + n c ) time, where h c is the totalnumber of marked generators in the wavefronts W ( e ) of all transparent edges e of c and n c is the totalnumber of hyperbolic-arcs of the bisector boundaries of these active regions. Computing the vertices of SPM (cid:48) ( s ) in each active region. In what follows, we compute thevertices of SPM (cid:48) ( s ) in each active region R of c .The boundary ∂R consists of O (1) pieces, each of which is a transparent edge fragment, an el-ementary chain fragment, or a bisector in SPM (cid:48) ( s ). Unlike the HS algorithm, where each of thesepieces is of O (1) size, in our problem both an elementary chain fragment and a bisector may not beof constant size. Let e be a transparent edge of R . Without loss of generality, we assume that e ishorizontal and R is locally above e . Without considering other wavefronts, we use W ( e ) to partition R into subregions, called Voronoi faces , such that each face has a unique predecessor in W ( e ). We use V or ( e ) to denote the partition of R . Note that if α is the predecessor of a Voronoi face, then for eachpoint p in the face, its tangent to α must cross e and we assign p a weight that is equal to d ( α, p ).Also, it is possible that these faces together may not cover the entire R ; for those points outside thesefaces, their weights are ∞ .We now discuss how to compute the partition V or ( e ). To this end, we can use our wavefrontpropagation procedure to propagate W ( e ) inside R . To apply the algorithm, one issue is that theboundary of R may contain bisectors, which consists of hyperbolic-arcs instead of polygonal segments.To circumvent the issue, an observation is that the bisectors of W ( e ) do not intersect the bisectorson ∂R . Indeed, for each bisector B of ∂R , one of its defining generators is unmarked and thus doesnot participate in any bisector event in c , and therefore, B does not intersect any bisector in c .In light of the observation, we can simply apply the wavefront propagation procedure to propagate W ( e ) to c instead of R to partition c into Voronoi faces, denoted by V or c ( e ). The above observationimplies that each bisector on the boundary of R must lie in the same face of V or c ( e ). Hence, we cansimply cut V or c ( e ) to obtain V or ( e ) using the bisectors on the boundary of R . When we apply thewavefront propagation procedure to propagate W ( e ) to c , here we add an initial step to compute theintersection of e and B ( α, α (cid:48) ) for each pair of neighboring generators α and α (cid:48) of W ( e ) and use it asthe initial tracing-point z ( α, α (cid:48) ). Since each such intersection can be computed in O (log n ) time bythe bisector-line intersection operation in Lemma 7, this initial step takes O ( h e log n ) time, where h e is the number of generators of W ( e ). Hence, the total time for propagating W ( e ) into c to obtain V or c ( e ) is O ( h e log n + n B ( e )), where n B ( e ) is the number of hyperbolic-arcs of the bisectors of W ( e )in c . After having V or c ( e ), cutting it to obtain V or ( e ) can be easily done in O ( h e + n B ( e ) + n R )time, where n R is the number of hyperbolic-arcs on the bisectors of ∂R . Therefore, the total time forcomputing the partition V or ( e ) is O ( h e log n + n B ( e ) + n R ) time. Note that since the bisectors of ∂R do not intersect any bisectors involving the generators of W ( e ), all bisector events of W ( e ) in c areactually in R .After having the partition V or ( e ) for all transparent edges e of R , we can now compute verticesof SPM (cid:48) ( s ) in R . Consider a transparent edge e of ∂R . We process it as follows. For each transparentedge f of ∂R other than e , we merge the two partitions V or ( e ) and V or ( f ) by using the merge stepfrom the standard divide-and-conquer Voronoi diagram algorithm to compute the sub-region of R W ( e ) than to W ( f ). This can be done in O ( h e + n B ( e ) + n R + h f + n B ( f )) time by findinga finding a curve γ in R that consists of points equal to W ( e ) and W ( f ), and the algorithm is similarto the HS algorithm.Intersecting the results for all such f produces the region R ( e ) claimed by W ( e ) in R . Intersecting R ( e ) with V or ( e ) gives the vertices of SPM (cid:48) ( s ) in R that W ( e ) contributes. Repeating the above forall transparent edges e of ∂R gives the vertices of SPM (cid:48) ( s ) in R . Since ∂R has O (1) transparent edges,the total time is O ( h R log n + n B ( R ) + n R ), where h R is the number of generators in all wavefrontsof all transparent edges of ∂R and n B ( R ) is the total number of hyperbolic-arcs of the bisectors allwavefronts of all transparent edges of ∂R .We do the above for all active regions R of c , and then all vertices of SPM (cid:48) ( s ) in c can be computed.The total time is O ( h c log n + n R ( c ) + n B ( c )), where h c is the total number of marked generators inthe wavefronts W ( e ) of all transparent edges e of c , n R ( c ) is the total number of hyperbolic-arcs ofthe bisector boundaries of the active regions of c , and n B ( c ) is the total number of hyperbolic-arcsof the bisectors of the wavefronts of all transparent edges of c . Processing all cells c of S (cid:48) as abovegives all vertices of SPM (cid:48) ( s ). For the running time, the total sum of n R ( c ) among all cells c ∈ S (cid:48) is O ( n ) because each hyperbolic-arc of a bisector on the boundary of any active region also appearsin SPM (cid:48) ( s ), whose total size is O ( n ) by Corollary 1. The total sum of h c among all cells c is O ( h )by Lemma 12. For n B ( c ), Lemma 14 is essentially a proof that the total sum of n B ( c ) over all cells c is O ( n ). To see this, since c is a cell incident to e , the wavefront W ( e ) will be propagated into c during the wavefront propagation procedure to propagate W ( e ) to edges of output ( e ), and thus thehyperbolic-arcs of the bisectors of W ( e ) are counted in the proof analysis of Lemma 14. In summary,the total time for computing all vertices of SPM (cid:48) ( s ) is O ( n + h log n ), which is O ( n + h log h ). SPM (cid:48) ( s )With all vertices of SPM (cid:48) ( s ) computed above, as in the HS algorithm, we next compute the edgesof SPM (cid:48) ( s ) separately and then assemble them to obtain SPM (cid:48) ( s ). One difference is that the HSalgorithm uses a standard plane sweep algorithm to assemble all edges to obtain SPM ( s ), which takes O ( n log n ) time as there are O ( n ) edges in SPM ( s ). In order to achieve the desired O ( n + h log h ) timebound, here instead we propose a different algorithm.We first discuss how to compute the edges of SPM (cid:48) ( s ), given the vertices of SPM (cid:48) ( s ). Recall thateach vertex v of SPM (cid:48) ( s ) is either an intersection between a bisector and an obstacle edge, or anintersection of bisectors. By our general position assumption, in the former case, v is in the interiorof an obstacle edge; in the latter case, v is the intersection of three bisectors, i.e., a triple point .During our above algorithm for computing vertices of SPM (cid:48) ( s ), we associate each vertex v with thetwo generators of the bisector that contains v (more specifically, we associate v with the initial verticesof the generators). We create a list of all vertices of SPM (cid:48) ( s ), each identified by a key consisting ofits two generators. If a vertex v is a triple point, then we put it in the list for three times (each timewith a different pair of generators); if v is a bisector-obstacle intersection, then it is put in the listonce. We now sort all vertices of SPM (cid:48) ( s ) with their keys; by traversing the sorted list, we can grouptogether all vertices belong to the same bisector. This sorting takes O ( h log n ) time as there are O ( h )vertices in SPM (cid:48) ( s ).We take all SPM (cid:48) ( s )-vertices in the same bisector and sort them along the bisector determined bytheir weighted distances from the generators of the bisector (each of these distances can be computedin additional O (log n ) time by computing the tangents from the vertex to the generators). Thesesorting takes O ( h log n ) time altogether. The above computes the bisector edges e of SPM (cid:48) ( s ) thatconnect two adjacent vertices. In fact, each edge e is only implicitly determined in the sense thatthe hyperbolic-arcs of e are not explicitly computed. In addition, for each obstacle P , we sort allvertices of SPM (cid:48) ( s ) on ∂P to compute the convex-chain edges of SPM (cid:48) ( s ) connecting adjacent verticesof SPM (cid:48) ( s ) on ∂P . This sorting can be easily done in O ( n + h log n ) time for all obstacles.40 ub Fu (cid:48) v (cid:48) b (cid:48) Figure 27: Illustrating v , u , u (cid:48) , v (cid:48) , b (cid:48) , and b . For each vertex v of SPM (cid:48) ( s ), the above computes its adjacent vertices. As discussed above, dueto the general position assumption, v has at most three adjacent vertices. We next define a set E ( v )of at most three points for v . For each adjacent vertex u of v in SPM (cid:48) ( s ), let e ( v, u ) be the edgeof SPM (cid:48) ( s ) connecting them. e ( v, u ) is either a bisector edge or a convex-chain edge. In the formercase, we refer to each hyperbolic-arc of e ( v, u ) as a piece of e ( v, u ); in the latter case, we refer to eachobstacle edge of e ( v, u ) as a piece . If we traverse e ( v, u ) from v to u , the endpoint u (cid:48) of the first piece(incident to v ) is added to E ( v ); we define h ( u (cid:48) ) to be u . Since v is adjacent to at most three verticesin SPM (cid:48) ( s ), | E ( v ) | ≤ 3. In fact, E ( v ) is exactly the set of vertices adjacent to v in the shortest pathmap SPM ( s ).With all edges of SPM (cid:48) ( s ) and the sets E ( v ) of all vertices v of SPM (cid:48) ( s ) computed above, weconstruct the doubly-connected-edge-list (DCEL) of SPM (cid:48) ( s ), as follows. For convenience, we assumethat there is a bounding box that contains all obstacles so that no bisector edge of SPM (cid:48) ( s ) will go toinfinity and thus each face of SPM (cid:48) ( s ) is bounded.For each vertex v of SPM (cid:48) ( s ), we store its coordinate in the DCEL data structure. As | E ( v ) | ≤ SPM (cid:48) ( s ) incident to v , we construct them one by one. For eachpoint u (cid:48) in E ( v ), we construct the face F clockwise (with respect to v ) incident to the edge e ( v, u )of SPM (cid:48) ( s ) connecting v and u = h ( u (cid:48) ), as follows (e.g., see Fig. 27). We first trace out the edge e ( v, u ), which is either a bisector edge or a convex chain edge. In the former case, we compute thehyperbolic-arcs of e ( v, u ), one at a time, until we reach u , and add them to the DCEL data structure.Each hyperbolic-arc can be computed in constant time using the two generators of the bisector. Inthe latter case, we trace out the obstacle edges of e ( v, u ), one at a time, and add them to the DCELdata structure. Let v (cid:48) be the last endpoint of the piece of e ( v, u ) containing u (e.g., see Fig. 27). Notethat v (cid:48) is in the set E ( u ). Let b (cid:48) be the first point of E ( u ) counterclockwise around u after v (cid:48) (e.g., seeFig. 27). Let b = h ( b (cid:48) ), which is a vertex of SPM (cid:48) ( s ) adjacent to u . Hence, the edge e ( u, b ) of SPM (cid:48) ( s )connecting u to b is incident to the face F . We trace out the edge e ( u, b ) in the same way as above.When we reach b , we continue to trace the next edge of F . Since each face of SPM (cid:48) ( s ) is bounded,eventually we will arrive back to the vertex v again . This finishes the construction of the face F .We do the same for all other faces of SPM (cid:48) ( s ), after which the DCEL data structure for SPM (cid:48) ( s ) isconstructed. For the running time, since each edge of SPM (cid:48) ( s ) is traced at most twice, by Corollary 1,the total time of the above procedure for constructing the DCEL data structure is O ( n ).In summary, SPM (cid:48) ( s ) can be computed in O ( n + h log h ) time. By Lemma 4, the shortest pathmap SPM ( s ) can be built in additional O ( n ) time. We could lift the assumption that each face of SPM (cid:48) ( s ) is bounded in the following way. During the above algorithmfor constructing F , if F is unbounded, then we will reach a bisector edge that extends to the infinity. If that happens,then we construct other edges of F from the other direction of v . More specifically, starting from the first point of E ( v )clockwise around v after u (cid:48) , we trace out the edges of F in the same way as above, until we reach a bisector edge thatextends to the infinity, at which moment all edges of F are constructed. .9 Reducing the space to O ( n ) The above provides an algorithm for computing SPM ( s ) in O ( n + h log h ) time and O ( n + h log h )space. In this subsection, we discuss how to reduce to the space to O ( n ), using the technique givenin [46].The reason that the above algorithm needs O ( n + h log h ) space is two-fold. First, it uses fullypersistent binary trees (with the path-copying method) to represent wavefronts W ( e ). Because thereare O ( h ) bisector events in the wavefront expansion algorithm and each event costs O (log n ) additionalspace on a persistent tree, the total space of the algorithm is O ( n + h log n ). Second, in order toconstruct SPM ( s ) after the wavefront expansion algorithm, the wavefronts W ( e ) of all transparentedges e of S (cid:48) are needed, which are maintained in those persistent trees. We resolve these two issuesin the following way. We still use persistent trees to represent wavefronts. However, as there are O ( h ) bisector events intotal in the algorithm, we divide the algorithm into O (log h ) phases so that each phase has no morethan h/ log h events. The total additional space for processing the events using persistent trees ineach phase is O ( h ). At the end of each phase, we “reset” the space of the algorithm by only storinga “snapshot” of the algorithm (and discarding all other used space) so that (1) the snapshot containssufficient information for the subsequent algorithm to proceed as usual, and (2) the total space of thesnapshot is O ( h ).Specifically, we make the following changes to the wavefront propagation procedure, which is tocompute the wavefronts W ( e, g ) for all edges g ∈ output ( e ) using the wavefront W ( e ). We nowmaintain a counter count to record the number of bisector events that have been processed so far sincethe last space reset; count = 0 initially. Consider a wavefront propagation procedure on the wavefront W ( e ) of a transparent edge e . The algorithm will compute W ( e, g ) for all edges g ∈ output ( e ), bypropagating W ( e ). We apply the same algorithm as before. For each bisector event, we first do thesame as before. Then, we increment count by one. If count < h/ log h , we proceed as before (i.e.,process the next event). Otherwise, we have reached the end of the current phase and will start anew phase. To do so, we first reset count = 0 and then reset the space by constructing and storing asnapshot of the algorithm (other space occupied by the algorithm is discarded), as follows.1. Let g refer to the edge of output ( e ) whose W ( e, g ) is currently being computed in the algorithm.We store the tree that is currently being used to compute W ( e, g ) right after the above event.To do so, we can make a new tree by copying the newest version of the current persistent treethe algorithm is operating on. The size of the tree is bounded by O ( h ). We will use this tree to“resume” computing W ( e, g ) in the subsequent algorithm.2. For each g (cid:48) ∈ output ( e ) \ { g } whose W ( e, g (cid:48) ) has been computed, we store the tree for W ( e, g (cid:48) ).We will use the tree to compute the wavefronts W ( g (cid:48) ) of g (cid:48) in the subsequent algorithm.3. We store the tree for the wavefront W ( e ). Note that the tree may have many versions due toprocessing the events and we only keep its original version for W ( e ). Hence, the size of the treeis O ( h ). This tree will be used in the subsequent algorithm to compute W ( e, g (cid:48) ) for those edges g (cid:48) ∈ output ( e ) \ { g } whose W ( e, g (cid:48) ) have not been computed yet.4. We check every transparent edge e (cid:48) of S (cid:48) with e (cid:48) (cid:54) = e . If e (cid:48) has been processed (i.e., the wavefrontpropagation procedure has been called on W ( e (cid:48) )) and there is an edge g (cid:48) ∈ output ( e (cid:48) ) that has not been processed, we know that W ( e (cid:48) , g (cid:48) ) has been computed and is available; we store thetree for W ( e (cid:48) , g (cid:48) ). We will use the tree to compute the wavefronts W ( g (cid:48) ) of g (cid:48) in the subsequentalgorithm. 42e refer to the wavefronts stored in the algorithm as the snapshot ; intuitively, the snapshot containsall wavelets in the forefront of the wavelet expansion. By the same analysis as in [46], we can showthat the snapshot contains sufficient information for the subsequent algorithm to proceed as usual andthe total space of the snapshot is O ( h ).The above discusses our changes to the wavefront propagation procedure. For the wavefrontmerging procedure, which is to construct W ( e ) from W ( f, e ) for the edges f ∈ input ( e ), notice thatwe do not need the old versions of W ( f, e ) anymore after W ( e ) is constructed. Therefore, it is notnecessary to use the path-copying method to process each event in the procedure. Hence, the totalspace needed in the wavefront merging procedure in the entire algorithm is O ( n ). SPM (cid:48) ( s )For the second issue of constructing SPM (cid:48) ( s ), our algorithm relies on the wavefronts W ( e ) for alltransparent edges e , which are maintained by persistent trees. Due to the space-reset, our algorithmdoes not maintain the wavefronts anymore, and thus we need to somehow restore these wavefrontsin order to construct SPM (cid:48) ( s ). To this end, a key observation is that by marking a total of O ( h )additional wavelet generators it is possible to restore all historical wavefronts that are needed forconstructing SPM (cid:48) ( s ). In this way, SPM (cid:48) ( s ) can be constructed in O ( n + h log h ) time and O ( n )space.More specifically, our algorithm considers each cell c of S (cid:48) individually. For each cell c , the algorithmhas two steps. First, compute the active regions of c . Second, for each active region R , compute thevertices of SPM (cid:48) ( s ) in R . For both steps, our algorithm utilizes the wavefronts W ( e ) of the transparentedges e on the boundary of c . Due to the space reset, the wavefronts W ( e ) are not available anymorein our new algorithm. We use the following strategy to resolve the issue.First, to compute the active regions in c , we need to know the bisectors defined by an unmarkedgenerator α and a marked generator α (cid:48) in the wavefronts W ( e ) of the transparent edges e of c . Weobserve that α is a generator adjacent to α (cid:48) in W ( e ). Based on this observation, we slightly modifyour wavefront expansion algorithm so that it also marks the neighbors of the generators that aremarked in our original algorithm and we call them newly-marked generators (the generators marked inour original algorithm are called originally-marked generators). If a generator is both newly-markedand originally-marked, we consider it as originally-marked. As each generator has two neighbors in awavefront, the total number of marked generators is still O ( h ). The newly-marked generators and theoriginally-marked generators are sufficient for computing all active regions of each cell c . Indeed, theactive regions are decomposition of c by the bisectors of adjacent generators with one originally-markedand the other newly-marked in W ( e ) of the transparent edges e of c .Second, to compute the vertices of SPM (cid:48) ( s ) in each active region R of c , we need to restore thewavefronts W ( e ) of the transparent edges e on the boundary of R . To this end, we observe that W ( e )consists of exactly the originally-marked generators that claim e . Consequently, the same algorithmas before can be applied to construct SPM (cid:48) ( s ). When computing the partition V or ( e ) of R , wepropagate W ( e ) using the wavefront propagation procedure, which uses persistent trees to represent W ( e ). Since here we do not need to keep the old versions of W ( e ) any more, we can use an ordinarytree (without the path-copying method) to represent W ( e ). In this way, processing each bisector eventonly introduces O (1) additional space. Hence, constructing SPM (cid:48) ( s ) takes O ( n + h log h ) time and O ( n ) space. In this section, we extend our algorithm for the convex case in Section 3 to the general case whereobstacles of P may not be convex. To this end, we resort to an extended corridor structure of P ,43igure 28: Illustrating a triangulation of P with twoobstacles. There are two junction triangles indicated bythe large dots inside them, connected by three solid (red)curves. Removing the two junction triangles results inthree corridors. cd bay ( cd ) a fb e afzycanal ( x, y ) xd eb Figure 29: Illustrating an open hourglass (left) and a closedone (right) with a corridor path connecting the apices x and y of the two funnels. The dashed segments are diagonals. Thepaths π C ( a, b ) and π C ( e, f ) are marked by thick solid curves.A bay bay ( cd ) with gate cd (left) and a canal canal ( x, y ) withgates xd and yz (right) are also shown. which has been used to solve various problems in polygonal domains [5, 8–10, 30, 38]. The structuredecomposes the free space F into three types of regions: an ocean M , O ( h ) canals , and O ( n ) bays .The details are given below. Let Tri ( P ) denote an arbitrary triangulation of P . Let G ( P ) be the (planar) dual graph of Tri ( P ),i.e., each node of G ( F ) corresponds to a triangle in Tri ( P ) and each edge connects two nodes of G ( P )corresponding to two triangles sharing a triangulation diagonal of Tri ( P ). We compute a corridorgraph G , as follows. First, repeatedly remove every degree-one node from G ( P ) until no such noderemains. Second, repeatedly remove every degree-two node from G ( P ) and replace its two incidentedges by a single edge until no such node remains. The resulting graph is G (e.g., see Fig. 28), whichhas O ( h ) faces, nodes, and edges [30]. Each node of G corresponds to a triangle in Tri ( P ), whichis called a junction triangle (e.g., see Fig. 28). The removal of all junction triangles results in O ( h ) corridors (defined below), and each corridor corresponds to an edge of G .The boundary of a corridor C consists of four parts (see Fig. 29): (1) A boundary portion of anobstacle, from an obstacle vertex a to an obstacle vertex b ; (2) a triangulation diagonal of a junctiontriangle from b to an obstacle vertex e ; (3) a boundary portion of an obstacle from e to an obstaclevertex f ; (4) a diagonal of a junction triangle from f to a . The two diagonals be and af are calledthe doors of C , and the other two parts of the boundary of C are the two sides of C . Note that C isa simple polygon. A point is in the interior of C if it is in C excluding the two doors. Let π C ( a, b )(resp., π C ( e, f )) be the shortest path from a to b (resp., e to f ) in C . The region H C bounded by π C ( a, b ) , π C ( e, f ), be , and f a is called an hourglass , which is open if π C ( a, b ) ∩ π C ( e, f ) = ∅ and closed otherwise (see Fig. 29). If H C is open, then both π C ( a, b ) and π C ( e, f ) are convex chains and calledthe sides of H C ; otherwise, H C consists of two “funnels” and a path π C = π C ( a, b ) ∩ π C ( e, f ) joiningthe two apices of the funnels, called the corridor path of C . Each side of a funnel is also convex.Let M be the union of all O ( h ) junction triangles, open hourglasses, and funnels. We call M the ocean , whose boundary ∂ M consists of O ( h ) convex chains that are sides of open hourglasses andfunnels. The other space of P , i.e., P \ M , is further partitioned into two types of regions: bays and canals , defined as follows. Consider the hourglass H C of a corridor C .If H C is open (see Fig. 29), then H C has two sides. Let S be a side of H C . The obstacle verticeson S all lie on the same side of C . Let c and d be any two consecutive vertices on S such that cd isnot an obstacle edge of P (e.g., see Fig. 29 left). The region enclosed by cd and the boundary portionof C between c and d is called a bay , denoted by bay ( cd ). We call cd the gate of bay ( cd ).If H C is closed, let x and y be the two apices of the two funnels. Consider two consecutive vertices44 and d on a side of a funnel such that cd is not an obstacle edge of P . If c and d are on the sameside of the corridor C , then cd also defines a bay. Otherwise, either c or d is a funnel apex, say, c = x ,and we call xd a canal gate at x = c (e.g., see Fig. 29 right). Similarly, there is also a canal gate atthe other funnel apex y , say yz . The region of C between the two canal gates xd and yz is the canal of H C , denoted by canal ( x, y ).Each bay or canal is a simple polygon. All bays and canals together constitute the space P \ M .Each vertex of ∂ M is a vertex of P and each edge of ∂ M is either an edge of P or a gate of a bay/canal.Gates are common boundaries between M and bays/canals. After P is triangulated, M and all baysand canals can be obtained in O ( n ) time [30].The reason that the extended corridor structure can help find a shortest path is the following.Suppose we want to find a shortest s - t path for two points s and t . We consider s and t as two specialobstacles and build the extended corridor structure. If a shortest s - t path π ( s, t ) contains a point inthe interior of a corridor C , then π ( s, t ) must cross both doors of C and stay in the hourglasses of C , and further, if the hourglass is closed, then its corridor path must be contained in π ( s, t ). In fact, π ( s, t ) must be in the union of the ocean M and all corridor paths [30].In light of the above properties, we propose the following algorithm. Let s be a given sourcepoint. By considering s as a special obstacle of P , we construct the extended corridor structure of P .Consider any query point t , which may be in the ocean M , a bay bay ( cd ), or a canal canal ( x, y ). • If t ∈ M , then the union of M and all corridor paths contains a shortest s - t path. To handlethis case, we will build a shortest path map SPM ( M ) in M with respect to the union of M and all corridor paths. In face, SPM ( M ) is exactly the portion of SPM ( s ) in M , i.e., SPM ( s ) ∩ M . To build SPM ( M ), a key observation is that the boundary ∂ M consists of O ( h ) convex chains. Therefore, we can easily adapt our previous algorithm for the convex case.However, the algorithm needs to be modified so that the corridor paths should be taken intoconsideration. Intuitively, corridor paths provide certain kind of “shortcuts” for wavefronts topropagate. • If t is in a bay bay ( cd ), then any shortest s - t path must cross its gate cd . To handle this case,we will extend SPM ( M ) into bay ( cd ) through the gate cd to construct the shortest path map in bay ( cd ), i.e., the portion of SPM ( s ) in bay ( cd ), SPM ( s ) ∩ bay ( cd ). • If t is in a canal canal ( x, y ), then any shortest s - t path must cross one of the two gates of thecanal. To handle this case, we will extend SPM ( M ) into canal ( x, y ) through the two gatesto construct the shortest path map in canal ( x, y ), i.e., the portion of SPM ( s ) in canal ( x, y ), SPM ( s ) ∩ canal ( x, y ).In the following, we first describe our algorithm for constructing SPM ( M ) in Section 4.2. We thenexpand SPM ( M ) into all bays in Section 4.3 and expand SPM ( M ) into all canals in Section 4.4. Thealgorithm for the canal case utilizes the bay case algorithm as a subroutine. SPM ( M ) As the boundary of M consists of O ( h ) convex chains, we can apply and slightly modify our algorithmfor the convex case. To do so, for each convex chain of ∂ M , we define its rectilinear extreme verticesin the same way as before. Let V be the set of the rectilinear extreme vertices of all convex chains.Hence, |V| = O ( h ). In addition, to incorporate the corridor paths into the algorithm, we include theendpoints of each corridor path in V . As there are O ( h ) corridor paths, the size of V is still boundedby O ( h ). Note that each corridor path endpoint is also an endpoint of a convex chain of ∂ M . Weconstruct the conforming subdivision S based on the points of V and then insert the convex chainsof M into S to obtain S (cid:48) . The algorithm is essentially the same as before. In addition, we make45he following changes to S (cid:48) , which is mainly for incorporating the corridor paths into our wavefrontexpansion algorithm, as will be clear later.Let v be an endpoint of a corridor path π . Since v is in V , v is incident to O (1) transparent edgesin S (cid:48) . For each such transparent edge e , if | π | < · | e | , then we divide e into two sub-edges such thatthe length of the one incident to v is equal to | π | / 2; for each sub-edge, we set its well-covering regionthe same as U ( e ). Note that this does not affect the properties of S (cid:48) . In particular, each transparentedge e is still well-covered. This change guarantees the following property: for each corridor path π , | π | ≥ · | e (cid:48) | holds, where e (cid:48) is any transparent edge of S (cid:48) incident to either endpoint of π . For referencepurpose, we refer to it as the corridor path length property .Next we apply the wavefront expansion algorithm. Here we need to incorporate the corridor pathsinto the algorithm. Intuitively, each corridor path provides a “shortcut” for the wavefront, i.e., if awavelet hits an endpoint of a corridor path, then the wavelet will come out of the corridor path fromits other endpoint but with a delay of distance equal to the length of the corridor path. More detailsare given below.Since all corridor path endpoints are in V , they are vertices of transparent edges of S (cid:48) . Consideran endpoint v of a corridor path π . Let u be the other endpoint of π . Recall that the wavefrontpropagation procedure for W ( e ) is to propagate W ( e ) to compute W ( e, g ) for all edges g ∈ output ( e ).In addition to the previous algorithm for the procedure, we also propagate W ( e ) through the corridorpath π to u . This is done as follows. Recall that when e is processed, since v is an endpoint of e , theweighted distance of v through W ( e ) is equal to d ( s, v ). Hence, the wavefront W ( e ) hits u through π at time d ( s, v ) + | π | . We then update covertime ( g ) = min { covertime ( g ) , d ( s, v ) + | π | + | g |} , for eachtransparent edge g incident to u . We also set the wavefront W ( e, g ) consisting of the only wavelet with u as the generator with weight equal to d ( s, v ) + | π | . Since there are O (1) transparent edges g incidentto u , the above additional step takes O (1) time, which does not change the time complexity of theoverall algorithm asymptotically. The corridor path length property assures that if W ( e ) contributesto a wavefront W ( g ) at g , then e must be processed earlier than g . This guarantees the correctness ofthe algorithm.In this way, we can first construct a decomposition SPM (cid:48) ( M ) of M in O ( n + h log h ) time and O ( n )space, where SPM (cid:48) ( M ) is defined similarly as SPM (cid:48) ( s ) in Section 3. Then, by a similar algorithm asthat for Lemma 4, SPM ( M ) can be obtained in additional O ( n ) time. SPM ( M ) into all bays We now expand SPM ( M ) into all bays in O ( n + h log h ) time and O ( n ) space. In fact, we expand SPM (cid:48) ( M ) to the bays. We process each bay individually. Consider a bay bay ( cd ) with gate cd .Without loss of generality, we assume that cd is horizontal, c is to the left of d , and bay ( cd ) is locallyabove cd .Let v , v , . . . , v m be the vertices of SPM (cid:48) ( M ) on cd ordered from left to right (e.g., see Fig. 30). Let c = v and d = v m +1 . Hence, each v i v i +1 is claimed by a generator α i = ( A i , a i ) for all i = 0 , , . . . , m .Let b i and c i be the tangent points on A i − and A i from v i , respectively, for each i = 1 , , . . . , m (e.g.,see Fig. 30). For v , only c is defined; for v m +1 , only b m +1 is defined. Observe that for any point p ∈ v i v i +1 , which is claimed by α i , its tangent point on A i must be on the portion of A i between c i and b i +1 and we use A (cid:48) i to denote that portion. So with respect to bay ( cd ), we use α (cid:48) i = ( A (cid:48) i , a (cid:48) i )to refer to the generator, where a (cid:48) i refers to the one of c i and b i +1 that is closer to a i . Hence, forany point t ∈ bay ( cd ), any shortest path π ( s, t ) from s to t must be via one of the generators α (cid:48) i for i = 0 , , . . . , m . Consider the region R bounded by A (cid:48) i for all i ∈ [0 , m ], the tangents from v i to theirgenerators for all i ∈ [0 , m + 1], and the boundary of the bay excluding its gate. Notice that R is asimple polygon. For any point t ∈ bay ( cd ), the above observation implies any shortest s - t path π ( s, t )is the concatenation of a shortest path from s to a generator initial vertex a (cid:48) i and the shortest pathfrom a (cid:48) i to t in R . 46 = v d = v m +1 v i b i c i v i +1 b i +1 b m +1 c B ( α i − , α i ) A (cid:48) i bay ( cd ) Figure 30: Illustrating bay ( cd ) and the generators. The thick segments on obstacles are A (cid:48) i , i = 0 , , . . . , m . According to the above discussion, expanding SPM (cid:48) ( M ) into bay ( cd ) is equivalent to the followingweighted geodesic Voronoi diagram problem in a simple polygon: Partition R into cells with respect tothe point sites a (cid:48) , a (cid:48) , . . . , a (cid:48) m (with weights equal to their geodesic distances to s ) such that all points inthe same cell have the same closest site. Let n b be the number of vertices in bay ( cd ) (the subscript “b”represents “bay”). Let n g be the total number of obstacle vertices in A (cid:48) i for all i ∈ [0 , m ] (the subscript“g” represents “generator”). Note that v i for all i = 1 , . . . , m are also vertices of R . Hence, the numberof vertices of R is n b + n g + m . The above problem can be solved in O ( m log m + n b + n g ) time by thetechniques of Oh [39]. Indeed, given m (cid:48) point sites in a simple polygon P (cid:48) of n (cid:48) vertices, Oh [39] gavean algorithm that can compute the geodesic Voronoi diagram of the sites in P (cid:48) in O ( n (cid:48) + m (cid:48) log m (cid:48) )time and O ( n (cid:48) + m (cid:48) ) space. Although the point sites in Oh’s problem do not have weights, ourproblem is essentially an intermediate step of Oh’s algorithm because all weighted point sites in ourproblem are on one side of cd . Therefore, we can run Oh’s algorithm from “the middle” and solve ourproblem in O ( m log m + n b + n g ) and O ( n b + n g ) space. In fact, our problem is a special case of Oh’sproblem because there are no sites in bay ( cd ). For this reason, we propose our own algorithm to solvethis special case and the algorithm is much simpler than Oh’s algorithm; our algorithm also runs in O ( m log m + n b + n g ) and O ( n b + n g + m ) space. This also makes our paper more self-contained.Before presenting the algorithm, we analyze the total time for processing all bays. Since SPM (cid:48) ( M )has O ( h ) vertices, the total sum of m for all bays is O ( h ). The total sum of n b for all bays is at most n . Notice that the obstacle edges on A (cid:48) i are disjoint for different bays, and thus the total sum of n g for all bays is O ( n ). Hence, expanding SPM (cid:48) ( M ) to all bays takes O ( n + h log h ) time and O ( n ) spacein total.The above actually only considers the case where the gate cd contains at least one vertex of SPM (cid:48) ( M ). It is possible that no vertex of SPM (cid:48) ( M ) is on cd , in which case the entire gate is claimedby one generator α of SPM (cid:48) ( M ). We can still define the region R in the same way. But now R hasonly one weighted site and thus the geodesic Voronoi diagram problem becomes computing a shortestpath map in the simple polygon R for a single source point. This problem can be solved in O ( n b + n g )time [22]; note that m = 0 in this case. Hence, the total time for processing all bays in this specialcase is O ( n ). We present an algorithm for the above special case of the weighted geodesic Voronoi diagram problemand the algorithm runs in O ( m log m + n b + n g ) and O ( n b + n g + m ) space.47f we consider all generators α i for i = 0 , , . . . , m as a wavefront at cd , denoted by W ( cd ), thenour algorithm is essentially to propagate W ( cd ) inside bay ( cd ). To this end, we first triangulate bay ( cd ) and will use the triangulation to guide the wavefront propagation. Each time we propagatethe wavefront through a triangle. In the following, we first discuss how to propagate W ( cd ) throughthe triangle (cid:52) cda with cd as an edge and a as the third vertex. This is a somewhat special case as (cid:52) cda is the first triangle the wavefront will propagate through; later we will discuss the general casebut the algorithm is only slightly different.Recall that each convex chain A (cid:48) i is represented by an array and the generator list W ( cd ) isrepresented by a balanced binary search tree T ( W ( cd )). We build a point location data structure onthe triangulation of bay ( cd ) in O ( n b ) time [17, 31], so that given any query point p , we can determinethe triangle that contains p in O (log n b ) time.We begin with computing the intersection of the adjacent bisectors B ( α i − , α i ) and B ( α i , α i +1 ) forall i = 1 , , . . . , m − 1. Each intersection can be computed in O (log n g ) time by the bisector-bisectorintersection operation in Lemma 8. Computing all intersections takes O ( m log n g ) time. For eachintersection q , called a bisector event , we use the point location data structure to find the triangle ofthe triangulation that contains q and store q in the triangle.Since all generators are outside bay ( cd ), by Corollary 3, all bisectors are monotone with respect tothe direction orthogonal to cd . Our algorithm for propagating the wavefront through (cid:52) cda is basedon this property. We sort all bisector events in (cid:52) cda according to their perpendicular distances to thesupporting line of cd . Then, we process these events in the same way as in our wavefront propagationprocedure. Specifically, for each bisector event q of two bisectors B ( α (cid:48) , α ) and B ( α, α (cid:48)(cid:48) ), we remove α from the generator list. Then, we compute the intersection q (cid:48) of B ( α (cid:48) , α (cid:48)(cid:48) ) with B ( α (cid:48) , α (cid:48) ), where α (cid:48) isthe other neighboring bisector of α (cid:48) than α , and we use the point location data structure to find thetriangle that contains q (cid:48) and store q (cid:48) in the triangle. If q (cid:48) ∈ (cid:52) cda , then we insert it to the bisectorevent sorted list of (cid:52) cda . We do the same for B ( α (cid:48) , α (cid:48)(cid:48) ) and B ( α (cid:48)(cid:48) , α (cid:48)(cid:48) ), where α (cid:48)(cid:48) is the other neighborof α (cid:48)(cid:48) than α . After all events in (cid:52) cda are processed, we split the current wavefront W at the vertex a . To this end, we first find the generator α ∗ of W that claims a . For this, we use the same algorithmas before, i.e., binary search plus bisector tracing. So we need to maintain a tracing-point for eachbisector as before (initially, we can set v i as the tracing-point for B ( α i − , α i ), i = 1 , , . . . , m ). Thecorrectness of the above algorithm for finding the generator α ∗ relies on the property of the followinglemma. Lemma 22 The bisector B ( α, α (cid:48) ) intersects ad ∪ ac at most once for any two bisectors of α and α (cid:48) of W ( cd ) . Proof: Note that we cannot apply the result of Lemma 5 since to do so we need to have ad and ac parallel to cd . But the proof is somewhat similar to that for Lemma 5.Assume to the contrary that B ( α, α (cid:48) ) intersects ad ∪ ac at two points, q and q . Let A and A (cid:48) bethe underlying arcs of α and α (cid:48) , respectively. Let v and u be the tangents points of q on A and A (cid:48) ,respectively (e.g., see Fig. 31). Let v and u be the tangents points of q on A and A (cid:48) , respectively.Since both A and A (cid:48) are on one side of cd while ad ∪ ac is on the other side, if we move a point q from q to q on ad ∪ ac , the tangent from q to A will continuously change from q v to q v and the tangentfrom q to A (cid:48) will continuously change from q u to q u . Therefore, either q u intersects q v in theirinteriors or q v intersects q u in their interiors; without loss of generality, we assume that it is theformer case. Let p be the intersection of q u and q v (e.g., see Fig. 31). Since q ∈ B ( α, α (cid:48) ), pointsof q u other than q have only one predecessor, which is α (cid:48) . As p ∈ q u and p (cid:54) = q , p has only onepredecessor α (cid:48) . Similarly, since q ∈ B ( α, α (cid:48) ) and p ∈ q v , α is also p ’s predecessor. We thus obtaina contradiction. (cid:3) After α ∗ is found, depending on whether α ∗ is the first or last generator of W , there are threecases. 48 A (cid:48) q v v u u p c daq Figure 31: Illustrating the proof of Lemma 22. A ∗ a p p p c d Figure 32: Splitting the generator α ∗ . Assume that A ∗ is the convex chain from p to p with p as the initial vertexof the generator. ap is tangent to A ∗ at p . After the split, the chain from p to p becomes a generator with p as theinitial vertex and the chain from p to p becomes another generator with p as the initial vertex. 1. If α ∗ is not the first or last generator of W , then we split W into two wavefronts, one for ca and the other for ad . To do so, we first split the binary tree T ( W ) that represents the currentwavefront W at α ∗ . Then, we do binary search on A ∗ to find the tangent point from a , where A ∗ is the underlying chain of α ∗ . We also split α ∗ into two at the tangent point of A ∗ , i.e., split A ∗ into two chains that form two generators, one for ac and the other for ad (e.g., see Fig. 32). As A ∗ is represented by an array, splitting α ∗ can be performed in O (1) time by resetting the endindices of the chains in the array. This finishes the propagation algorithm in (cid:52) acd . The abovesplits W into two wavefronts, one for ac and the other for ad ; we then propagate the wavefrontsthrough ac and ad recursively.2. If α ∗ is the first generator of W , then α ∗ must be α , i.e., the leftmost generator of W ( cd ). Inthis case, we do not need to split W . But we still need to split the generator α at the tangentpoint p of α from a . To find the tangent point p , however, this time we do not use binarysearch as it is possible that we will need to do this for Ω( n b ) vertices of bay ( cd ), which wouldtake Ω( n b log n g ) time in total. Instead, we use the following approach. Recall that the vertex c = v connects to α by a tangent with tangent point c (e.g., see Fig. 30), and c is an endpointof A (cid:48) . We traverse A (cid:48) from c to b , i.e., the other endpoint of A (cid:48) ; for each vertex, we checkwhether it is the tangent from a . In this way, p can be found in time linear in the number ofvertices of A (cid:48) between c and the tangent point p (e.g., see Fig. 33). After that, we still split α into two generators at p ; one is the only generator for ac and the other becomes the firstgenerator of the wavefront for ad .For ad , we propagate its wavefront through ad recursively. For ac , to propagate its wavefrontthrough ac , since the wavefront has only one generator, we can simply apply the linear timeshortest path map algorithm for simple polygons [22]; we refer to this as the one-generator case .49 (cid:48) ac dc b p R ( ac ) Figure 33: Illustrating the case for propagating the one-generator wavefront to R ( ac ). Indeed, ac partitions bay ( cd ) into two sub-polygons and let R ( ac ) denote the one that does notcontain (cid:52) cda (e.g., see Fig. 33). Hence, all points of R ( ac ) are claimed by the only generatorfor ac , whose underlying chain A (cid:48) is the sub-chain of A (cid:48) from c to p and whose initial vertexis p . Consider the region R (cid:48) ( ac ) bounded by cc , A (cid:48) , p a , and the boundary of R ( ac ) excluding ac . It is a simple polygon with a single weighted source p . Therefore, our problem is equivalentto computing the shortest path map in R (cid:48) ( ac ) with respect to the source point p , which can bedone in O ( | R ( ac ) | + | A (cid:48) | ) time [22], where | R ( ac ) | and | A (cid:48) | are the numbers of vertices of R ( ac )and | A (cid:48) | , respectively.3. If α ∗ is the last generator of W , the algorithm is symmetric to the above second case.The above describes our algorithm for propagating the wavefront W ( cd ) through the first triangle cda . Next, we discuss the general case where we propagate a wavefront W of more than one generatorthrough an arbitrary triangle of the triangulation of bay ( cd ). For the sake of notational convenience,we consider the problem of propagating the wavefront W ( ad ) at ad through ad into the region R ( ad ),where R ( ad ) is the one of the two sub-polygons of bay ( cd ) partitioned by ad that does not contain (cid:52) cda . Let (cid:52) adb be the triangle in R ( ad ) with ad with as an edge, i.e., b is the third vertex of thetriangle. We describe the algorithm for propagating W ( ad ) through (cid:52) adb . The algorithm is actuallyquite similar as before, with an additional event-validation step.We first sort all bisector events in (cid:52) adb , following their perpendicular distances to the supportingline of ad . Then we process these events as before. One difference is that we now need to check whethereach event is still valid. Specifically, for each event q , which is associated with three generators α , α (cid:48) ,and α (cid:48)(cid:48) , i.e., q is the intersection of the bisectors B ( α, α (cid:48) ) and B ( α (cid:48) , α (cid:48)(cid:48) ), we check whether all threegenerators are still in the current wavefront W . If not (this is possible if one of the three generatorswas deleted before, e.g., when the wavefront was propagated through (cid:52) cda ), then q is not valid andwe ignore this event; otherwise q is still valid and we process it in the same way as before.The above algorithm is based on the assumption that ad is a triangulation diagonal. If it is anobstacle edge, then the wavefront W ( ad ) stops at ad . Notice that each bisector of the wavefront W ( ad )must intersect ad . For each bisector of the wavefront, starting from its current tracing-point, we traceit out until the current traced hyperbolic-arc intersects ad .The algorithm stops once all triangles are processed as above. Time analysis. We now analyze the time complexity. As discussed before, the initial step forcomputing the intersections of adjacent bisectors of the wavefront W ( cd ) and locating the trianglescontaining them together takes O ( m log( n g + n b )) time. During the entire algorithm, each traced Note that the order of the generators in W must be consistent with their initial index order, α , α , . . . , α m . Hence,checking whether a generator is in W can be easily done in O (log m ) time by using the generator indices. bay ( cd ), i.e., the portion of SPM ( s ) in bay ( cd ), whose size is O ( n b + n g + m ). Hence, the total time on the bisector tracing in the entirealgorithm is O ( n b + n g + m ). For the one-generator case where only one generator is available for ac ,the time for processing the sub-polygon R ( ac ) is O ( | R ( ac ) | + | A (cid:48) | ). Notice that all such sub-polygons R ( ac ) in the one-generator case are interior disjoint. Hence, the total sum of their sizes is O ( n b ). Also,all such generator underlying chains A (cid:48) are also interior disjoint, and thus the total sum of their sizesis O ( n g ). Therefore, the overall time for processing the one-generator case sub-polygons is O ( n g + n b ).For the general case of processing a triangle (cid:52) of the triangulation, the total time for processingall events is O ( m (cid:48)(cid:48) log( n g + n b )) , where m (cid:48)(cid:48) is the number of events in (cid:52) , both valid and invalid.Each valid or invalid event is computed either in the initial step or after a generator is deleted. Thetotal number of bisector events in the former case in the entire algorithm is at most m − 1. Thetotal number in the latter case in the entire algorithm is no more than the number of generators thatare deleted in the algorithm, which is at most m because once a generator is deleted it will neverappear in any wavefront again. Hence, the total time for processing events in the entire algorithm is O ( m log( n g + n b )). Once all events in (cid:52) are processed, we need to find the generator α ∗ of the currentwavefront W that claims the third vertex b of the triangle, by binary search plus bisector tracing. Thetime is O (log m ) plus the time for tracing the bisector hyperbolic-arcs. Hence, excluding the time fortracing the bisector hyperbolic-arcs, which has been counted above, the total time for this operationin the entire algorithm is O ( m (cid:48) log m ), where m (cid:48) is the number of triangles the algorithm processed forthe case where α ∗ is not the first or last generator. We will show later in Lemma 23 that m (cid:48) = O ( m ).After α ∗ is found, depending on whether α ∗ is the first or last generator of W , there are threecases. If α ∗ is not the first or last generator of W , then we find the tangent from b to α ∗ by binarysearch in O (log n g ) time and split both W and α ; otherwise, we find the tangent by a linear scan on α ∗ and only need to split α ∗ . Splitting W takes O (log m ) time while splitting a generator only takes O (1)time as discussed before. Therefore, the total time for splitting generators in the entire algorithm is O ( n b ), as there are O ( n b ) triangles in the triangulation. If α ∗ is either the first or last generator of W ,then a one-generator case subproblem will be produced and the time of the linear scan for finding thetangent is O ( | A (cid:48) | ), where A (cid:48) is the sub-chain of α ∗ that belongs to the one-generator case subproblem.As discussed above, all such A (cid:48) in all one-generator case subproblems are interior disjoint, and thusthe total time on the linear scan in the entire algorithm is O ( n g ). Therefore, the total time for findingthe tangent point and splitting W is O ( m (cid:48) log n g + n b ) as m ≤ n g . Lemma 23 m (cid:48) ≤ m − . Proof: Initially W = W ( cd ), which consists of m generators. Each split operation splits a wavefront W at a generator α ∗ into two wavefronts each of which has at least two generators (such that bothwavefronts contain α ∗ ). More specifically, if a wavefront of size k is split, then one wavefront has k (cid:48) generators and the other has k − k (cid:48) + 1 generators with k (cid:48) ≥ k − k (cid:48) + 1 ≥ 2. The value m (cid:48) isequal to the number of all split operations in the algorithm.We use a tree structure T to characterize the split operations. The root of T corresponds to theinitial generator sequence W ( cd ). Each split on a wavefront W corresponds to an internal node of T with two children corresponding to the two subsequences of W after the split. Hence, m (cid:48) is equalto the number of internal nodes of T . In the worst case T has m leaves, each corresponding to twogenerators α i and α i +1 for i = 0 , , . . . , m . Since each internal node of T has two children and T hasat most m leaves, the number of internal nodes of T is at most m − 1. Therefore, m (cid:48) ≤ m − (cid:3) With the preceding lemma, the total time of the algorithm is O ( m log( n g + n b ) + n g + n b ), which More specifically, sorting all events takes O ( m (cid:48)(cid:48) log m (cid:48)(cid:48) ) time. For each event, removing a generator from W takes O (log m ), finding the intersections of new adjacent bisectors takes O (log n g ) time, and then locating the triangles con-taining the intersections takes O (log n b ) time. Note that m (cid:48)(cid:48) ≤ m ≤ n g . O ( m log m + n g + n b ) by similar analysis as Observation 1. The space is O ( n g + n b + m ) as nopersistent data structures are used. SPM (cid:48) ( M ) into all canals Consider a canal canal ( x, y ) with two gates xd and yz . The goal is to expand the map SPM (cid:48) ( M )into canal ( x, y ) through the two gates to obtain the shortest path map in the canal, denoted by SPM ( canal ( x, y )).The high-level scheme of the algorithm is similar in spirit to that for the L problem [10]. Thealgorithm has three main steps. First, we expand SPM (cid:48) ( M ) into canal ( x, y ) through the gate xd , byapplying our algorithm for bays. Let SPM ( canal ( x, y )) denote the map of canal ( x, y ) obtained by thealgorithm. Second, we expand SPM ( M ) into canal ( x, y ) through the gate yz by a similar algorithmas above; let SPM ( canal ( x, y )) denote the map of canal ( x, y ) obtained by the algorithm. Third, wemerge the two maps SPM ( canal ( x, y )) and SPM ( canal ( x, y )) to obtain SPM ( canal ( x, y )). This isdone by using the merge step from the standard divide-and-conquer Voronoi diagram algorithm tocompute the region closer to the generators at xd than those at yz . We will provide more details forthis step below and we will show that this step can be done in time linear in the total size of the twomaps SPM ( canal ( x, y )) and SPM ( canal ( x, y )). Before doing so, we analyze the complexities of thealgorithm.The first step takes O ( n + h log h ) time for all canals. So is the second step. The third step takeslinear time in the total size of SPM ( canal ( x, y )) and SPM ( canal ( x, y )). Since the total size of thetwo maps over all canals is O ( n ), the total time of the third step for all canals is O ( n ). In summary,the time for computing the shortest path maps in all canals is O ( n + h log h ) and the space is O ( n ).In the following, we provide more details for the third step of the algorithm.Recall that x and y are the two endpoints of the corridor path π in canal ( x, y ). It is possible thatthe shortest s - x path π ( s, x ) contains y or the shortest s - y path π ( s, y ) contains x . To determine that,we can simply check whether d ( s, x ) + | π | = d ( s, y ) and whether d ( s, y ) + | π | = d ( s, x ). Note that both d ( s, x ) and d ( s, y ) are available once SPM ( M ) is computed.We first consider the case where neither π ( s, x ) contains y nor π ( s, y ) contains x . In this case,there must be a point p ∗ in π such that d ( s, x ) + | π ( x, p ∗ ) | = d ( s, y ) + | π ( y, p ∗ ) | , where π ( x, p ∗ ) (resp., π ( y, p ∗ )) is the subpath of π between x (resp., y ) and p ∗ . We can easily find p ∗ in O ( | π | ) time. Tomerge the two maps SPM ( canal ( x, y )) and SPM ( canal ( x, y )) to obtain SPM ( canal ( x, y )), we find adividing curve Γ in canal ( x, y ) such that W ( xd ) claims all points of canal ( x, y ) on one side of Γ while W ( yz ) claims all points of canal ( x, y ) on the other side of Γ, where W ( xd ) is the set of generators of SPM (cid:48) ( M ) claiming xd (one may consider W ( xd ) is a wavefront) and W ( yz ) is the set of generatorsof SPM (cid:48) ( M ) claiming yz . The curve Γ consists of all points in canal ( x, y ) that have equal weighteddistances to W ( xd ) and W ( yz ). Therefore, the point p ∗ must be on Γ. Starting from p ∗ , we can trace Γout by walking simultaneously in the cells of the two maps SPM ( canal ( x, y )) and SPM ( canal ( x, y )).The running time is thus linear in the total size of the two maps.We then consider the case where either π ( s, x ) contains y or π ( s, y ) contains x . Without loss ofgenerality, we assume the latter case. We first check whether π ( s, p ) contains x for all points p on thegate yz . To do so, according to the definitions of corridor paths and gates of canals, for any point p ∈ yz , its shortest path to x in canal ( x, y ) is the concatenation of the corridor path π and yp . Hence,it suffices to check whether d ( s, z ) = d ( s, y ) + | yz | . • If yes, then all points of yz are claimed by the generators of W ( xd ) (in fact, they are claimedby x because for any point q ∈ xd and any point p ∈ yz , their shortest path in canal ( x, y ) isthe concatenation of qx , π , and yp ). Hence, all points in canal ( x, y ) are claimed by W ( xd ) and SPM ( canal ( x, y )) is SPM ( canal ( x, y )). 52 Otherwise, some points of yz are claimed by W ( xd ) while others are claimed by W ( yz ). As inthe above case, we need to find a dividing curve Γ consisting of all points with equal weighted dis-tances to W ( xd ) and W ( yz ). To this end, we again first find a point p ∗ in Γ. For this, since π ( s, y )contains x , y is claimed by W ( xd ). On the other hand, since d ( s, z ) (cid:54) = d ( s, y ) + | yz | , z is claimedby W ( yz ). Therefore, yz must contains a point p ∗ ∈ Γ. Such a point p ∗ can be found by travers-ing yz simultaneously in the cells of both SPM ( canal ( x, y )) and SPM ( canal ( x, y )). After p ∗ is found, we can again trace Γ out by walking simultaneously in the cells of SPM ( canal ( x, y ))and SPM ( canal ( x, y )). The running time is also linear in the total size of the two maps. Theorem 1 Suppose P is a set of h pairwise disjoint polygonal obstacles with a total of n vertices inthe plane and s is a source point. Assume that a triangulation of the free space is given. The shortestpath map SPM ( s ) with respect to s can be constructed in O ( n + h log h ) time and O ( n ) space. Proof: Using the triangulation, we decompose the free space F into an ocean M , canals, and baysin O ( n ) time [30]. Then, the shortest path map SPM ( M ) in the ocean M can be constructed in O ( n + h log h ) time and O ( n ) space. Next, SPM ( M ) can be expanded into all bays and canals inadditional O ( n + h log h ) time and O ( n ) space. The shortest path map SPM ( s ) is thus obtained. (cid:3) The current best algorithms can compute a triangulation of the free space in O ( n log n ) time or in O ( n + h log (cid:15) h ) time for any small (cid:15) > P are convex, then the triangulationcan be done in O ( n + h log h ) time [28].After SPM ( s ) is computed, by building a point location data structure [17, 31] on SPM ( s ) inadditional O ( n ) time, given a query point t , the shortest path length from s to t can be computed in O (log n ) time and a shortest s - t path can be produced in time linear in the number of edges of thepath. Corollary 4 Suppose P is a set of h pairwise disjoint polygonal obstacles with a total of n verticesin the plane and s is a source point. Assume that a triangulation of the free space is given. A datastructure of O ( n ) space can be constructed in O ( n + h log h ) time and O ( n ) space, so that given anyquery point t , the shortest path length from s to t can be computed in O (log n ) time and a shortest s - t path can be produced in time linear in the number of edges of the path. References [1] The open problem project. http://cs.smith.edu/ orourke/TOPP/.[2] S.W. Bae and H. Wang. L shortest path queries in simple polygons. Theoretical Computer Science , 790:105–116,2019.[3] R. Bar-Yehuda and B. Chazelle. Triangulating disjoint Jordan chains. International Journal of ComputationalGeometry and Applications , 4(4):475–481, 1994.[4] D.Z. Chen, J. Hershberger, and H. Wang. Computing shortest paths amid convex pseudodisks. SIAM Journal onComputing , 42(3):1158–1184, 2013.[5] D.Z. Chen, R. Inkulu, and H. Wang. Two-point L shortest path queries in the plane. Journal of ComputationalGeometry , 1:473–519, 2016.[6] D.Z. Chen, K.S. Klenk, and H.-Y.T. Tu. Shortest path queries among weighted obstacles in the rectilinear plane. SIAM Journal on Computing , 29(4):1223–1246, 2000.[7] D.Z. Chen and H. Wang. Computing shortest paths among curved obstacles in the plane. ACM Transactions onAlgorithms , 11, 2015. Article No. 26.[8] D.Z. Chen and H. Wang. A new algorithm for computing visibility graphs of polygonal obstacles in the plane. Journal of Computational Geometry , 6:316–345, 2015. 9] D.Z. Chen and H. Wang. Computing the visibility polygon of an island in a polygonal domain. Algorithmica ,77:40–64, 2017.[10] D.Z. Chen and H. Wang. Computing L shortest paths among polygonal obstacles in the plane. Algorithmica ,81:2430–2483, 2019.[11] Y.-J. Chiang and J.S.B. Mitchell. Two-point Euclidean shortest path queries in the plane. In Proceedings of theAnnual ACM-SIAM Symposium on Discrete Algorithms , pages 215–224, 1999.[12] K. Clarkson, S. Kapoor, and P. Vaidya. Rectilinear shortest paths through polygonal obstacles in O ( n log n ) time.In Proc. of the 3rd Annual Symposium on Computational Geometry , pages 251–257, 1987.[13] K. Clarkson, S. Kapoor, and P. Vaidya. Rectilinear shortest paths through polygonal obstacles in O ( n log / n )time. Manuscript, 1988.[14] K.L. Clarkson, R. Cole, and R.E. Tarjan. Randomized parallel algorithms for trapezoidal diagrams. InternationalJournal of Computational Geometry and Application , 2:117–133, 1992.[15] P.J. de Rezende, D.T. Lee, and Y.F. Wu. Rectilinear shortest paths in the presence of rectangular barriers. Discreteand Computational Geometry , 4:41–53, 1989.[16] J. Driscoll, N. Sarnak, D. Sleator, and R.E. Tarjan. Making data structures persistent. Journal of Computer andSystem Sciences , 38(1):86–124, 1989.[17] H. Edelsbrunner, L. Guibas, and J. Stolfi. Optimal point location in a monotone subdivision. SIAM Journal onComputing , 15(2):317–340, 1986.[18] S.D. Eriksson-Bique, J. Hershberger, V. Polishchuk, B. Speckmann, S. Suri, T. Talvitie, K. Verbeek, and H. Yıldız.Geometric k shortest paths. In Proceedings of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algo-rithms (SODA) , pages 1616–1625, 2015.[19] S.K. Ghosh and D.M. Mount. An output-sensitive algorithm for computing visibility graphs. SIAM Journal onComputing , 20(5):888–910, 1991.[20] L. Guibas, J. Hershberger, and J. Snoeyink. Compact interval trees: A data structure for convex hulls. InternationalJournal of Computational Geometry and Applications , 1(1):1–22, 1991.[21] L.J. Guibas and J. Hershberger. Optimal shortest path queries in a simple polygon. Journal of Computer andSystem Sciences , 39(2):126–152, 1989.[22] L.J. Guibas, J. Hershberger, D. Leven, M. Sharir, and R.E. Tarjan. Linear-time algorithms for visibility and shortestpath problems inside triangulated simple polygons. Algorithmica , 2(1-4):209–233, 1987.[23] J. Hershberger. A new data structure for shortest path queries in a simple polygon. Information Processing Letters ,38(5):231–235, 1991.[24] J. Hershberger, V. Polishchuk, B. Speckmann, and T. Talvitie. Geometric k th shortest paths. In Proceedings of the30th Annual Symposium on Computational Geometry (SoCG) ComputationalGeometry: Theory and Applications , 4(2):63–97, 1994.[26] J. Hershberger and S. Suri. An optimal algorithm for Euclidean shortest paths in the plane. SIAM Journal onComputing , 28(6):2215–2256, 1999.[27] J. Hershberger, S. Suri, and H. Yıldız. A near-optimal algorithm for shortest paths among curved obstacles in theplane. In Proceedings of the 29th Annual Symposium on Computational Geometry (SoCG) , pages 359–368, 2013.[28] S. Hertel and K. Mehlhorn. Fast triangulation of the plane with respect to simple polygons. Information andControl , 64:52–76, 1985.[29] R. Inkulu, S. Kapoor, and S.N. Maheshwari. A near optimal algorithm for finding Euclidean shortest path inpolygonal domain. In arXiv:1011.6481v1 , 2010.[30] S. Kapoor, S.N. Maheshwari, and J.S.B. Mitchell. An efficient algorithm for Euclidean shortest paths amongpolygonal obstacles in the plane. Discrete and Computational Geometry , 18(4):377–383, 1997.[31] D. Kirkpatrick. Optimal search in planar subdivisions. SIAM Journal on Computing , 12(1):28–35, 1983.[32] D. Kirkpatrick and J. Snoeyink. Tentative prune-and-search for computing fixed-points with applications to geo-metric computation. Fundamenta Informaticae , 22(4):353–370, 1995.[33] D.T. Lee and F.P. Preparata. Euclidean shortest paths in the presence of rectilinear barriers. Networks , 14(3):393–410, 1984.[34] J.S.B. Mitchell. An optimal algorithm for shortest rectilinear paths among obstacles. Abstracts of the , 1989. 35] J.S.B. Mitchell. A new algorithm for shortest paths among obstacles in the plane. Annals of Mathematics andArtificial Intelligence , 3(1):83–105, 1991.[36] J.S.B. Mitchell. L shortest paths among polygonal obstacles in the plane. Algorithmica , 8(1):55–88, 1992.[37] J.S.B. Mitchell. Shortest paths among obstacles in the plane. International Journal of Computational Geometryand Applications , 6(3):309–332, 1996.[38] J.S.B. Mitchell and S. Suri. Separation and approximation of polyhedral objects. Computational Geometry: Theoryand Applications , 5:95–114, 1995.[39] E. Oh. Optimal algorithm for geodesic nearest-point Voronoi diagrams in simple polygons. In Proceedings of the20th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA) , pages 391–409, 2019.[40] M. Overmars and J. van Leeuwen. Maintenance of configurations in the plane. Journal of Computer and SystemSciences , 23(2):166–204, 1981.[41] H. Rohnert. Shortest paths in the plane with convex polygonal obstacles. Information Processing Letters , 23(2):71–76, 1986.[42] M. Sharir and A. Schorr. On shortest paths in polyhedral spaces. SIAM Journal on Computing , 15(1):193–215,1986.[43] J.A. Storer and J.H. Reif. Shortest paths in the plane with polygonal obstacles. Journal of the ACM , 41(5):982–1012,1994.[44] H. Wang. Quickest visibility queries in polygonal domains. Discrete and Computational Geometry , 62:374–432,2019.[45] H. Wang. A divide-and-conquer algorithm for two-point L shortest path queries in polygonal domains. Journal ofComputational Geometry , 11:235–282, 2020.[46] H. Wang. Shortest paths among obstacles in the plane revisited. In Proceedings of the 32nd Annual ACM-SIAMSymposium on Discrete Algorithms (SODA) , pages 810–821, 2021., pages 810–821, 2021.