Shortest Paths in Geodesic Unit-Disk Graphs

In this problem, we preprocess a set $S$ of $n$ points inside a simple polygon $P$ of $m$ vertices to support two operations: (1) given a query point $q$, determine whether there exists a point $p\in S$ whose geodesic distance to $q$ is at most $1$, and if so return such a point; and (2) delete a point from $S$. We refer to the first operation as a geodesic unit-disk range emptiness (GUDRE) query.

Using the dynamic randomized nearest neighbor data structure of Agarwal, Arge, and Staals [1], each deletion can be supported in $O({\log }^{7}n\log m+{\log }^{6}n{\log }^{3}m)$ expected amortized time, and each GUDRE query can be answered in $O({\log }^{2}n{\log }^{2}m)$ time. In the Euclidean case (i.e., when $P$ is the entire plane), one can preprocess $S$ in $O(n\log n)$ time and support both queries and deletions in $O(\log n)$ time [14, 30].

We obtain the following result. After an $O(m)$-time global preprocessing step for $P$, we can, for any given set $S$ of $n$ points in $P$, build a data structure in $O(n{\log }^{3}m+n{\log }^{2}m\log n)$ time such that deletions take $O(\log n{\log }^{3}m)$ amortized time and GUDRE queries take $O(\log n{\log }^{2}m)$ time. Our result is remarkable in three aspects. First, compared to applying [1], our result improves deletion time by a factor of $O({\log }^{5}n)$ and improves query time by a factor of $O(\log n)$, while being deterministic rather than randomized (note that [1] also supports insertions and more general queries). Second, when applied to the Euclidean case (with $m=O(1)$), our bounds match the optimal Euclidean results of [14, 30]. Third, our result has an interesting implicit nature: after the global preprocessing of $P$, the construction time for any given $S$ depends only logarithmically on $m$.

In this problem, we wish to construct a nearest neighbor data structure for a set $S$ of $n$ additively weighted points in a simple polygon $P$ with $m$ vertices, under the requirement that the data structure be implicitly constructed. That is, after a global preprocessing of $P$, the time to build the data structure for $S$ should depend only logarithmically on $m$. Such a data structure is referred to as an implicit Voronoi diagram.

An explicit diagram for $S$ can be built in $O((n+m)\log (n+m))$ time [19], after which weighted nearest neighbor queries can be answered in $O(\log (n+m))$ time via standard point location queries [21, 13]. However, this explicit approach cannot be used in the implicit setting, as even constructing a point location data structure requires $\mathrm{\Omega }(m)$ time.

For the unweighted version of this implicit problem, the method in [1] obtains the following: after $O(m)$-time global preprocessing, one can build a structure for $S$ in $O(n\log n\log m+n{\log }^{3}m)$ time that answers nearest neighbor queries in $O(\log n{\log }^{2}m)$ time.

The weighted case is more challenging, as several geometric properties used in [1] no longer hold in the presence of weights. We develop a new method and show that, after $O(m)$-time global preprocessing of $P$, a data structure for $S$ can be built in $O(n\log n{\log }^{2}m+n{\log }^{3}m)$ time, and each weighted nearest neighbor query can be answered in $O(\log n{\log }^{2}m)$ time.

In the rest of the paper, after introducing preliminary concepts in Section 2, we present our algorithm for the weighted case in simple polygons in Section 3. This algorithm relies on the solutions to the two subproblems. Section 4 discusses the first subproblem. For the unweighted case, we describe our algorithm for simple polygons in Section 5, and extend it to polygons with holes in Section 6. Thus, Section 6 is the only section dealing with polygons with holes. Due to space constraints, many details and proofs are omitted but can be found in the full paper.

Although the algorithm looks simple, its correctness is not immediate. See the full paper for the proof. In particular, we show that the algorithm maintains the following two invariants. (I) For every site $v\in A$, we have $\mathrm{dist}[v]=d_G(v)$. (II) For every site $v\in A_1$, every site $u$ adjacent to $v$ in $G(S)$ is already in $A$. Using these, we can show that after the algorithm terminates, $\mathrm{dist}[v]=d_G(v)$ holds for every site $v\in S$.

To compute $B_a$, we repeatedly locate a point $b\in B$ with $d(a,b)\le 1$ and delete it from $B$, stopping once no such point remains. Thus, the task reduces exactly to supporting deletion-only geodesic unit-disk range emptiness (GUDRE) queries, as introduced in Section 1.1. Given a set $B$ of points in $P$, we require a data structure that supports:

Delete:

Remove a point from $B$.

GUDRE query:

Given a point $a\in P$, determine whether there exists $b\in B$ with $d(a,b)\le 1$, and if so, return such a point.

One of our main results in this paper is the following theorem, proved in Section 4.

Using Theorem 1, we implement Subroutine I as follows. Initially, $B=S\setminus \{s\}$. We build the GUDRE data structure for $B$ in $O(m+n{\log }^{3}m+n\log n{\log }^{2}m)$ time, where the $O(m)$ term is for constructing the GH data structure. Since the algorithm performs $O(n)$ deletions and GUDRE queries overall, the total time spent in Subroutine I is $O(m+n\log n{\log }^{3}m)$.

For the second subroutine, observe that each ${\sigma }_b$ is essentially the (additively) weighted nearest neighbor of $b$ in $A_2$, where each site $p\in A_2$ has weight $\mathrm{dist}[p]$, and the weighted distance from $p$ to a query point $q$ is defined as $\mathrm{dist}[p]+d(p,q)$. Thus, we require a weighted nearest neighbor data structure for $A_2$.

Furthermore, since at the end of each iteration, we will move all points of $B_a$ from $B$ to $A_2$ and move $a$ from $A_2$ to $A_1$, we need a dynamic data structure for $A_2$ that supports both insertions and deletions. A direct dynamic data structure supporting both insertions and deletions appears difficult to obtain. However, we show that the problem can be reduced so that only insertions are needed. The key is the following Lemma 2.

We maintain a balanced binary search tree $T_A$ whose leaves store the sites of $A$ in increasing order of $\mathrm{dist}[\cdot ]$. By Lemma 2, the points of $A$ can always be arranged so that all sites of $A_1$ appear to the left of all sites of $A_2$; we refer to this as the ordering property of $T_A$. For each node $v$ of $T_A$, let $A_v$ denote the set of sites whose leaves lie in the subtree rooted at $v$. For each $v$, we construct a weighted nearest neighbor data structure ${𝒟}_v$ for $A_v$, where each $p\in A_v$ has weight $\mathrm{dist}[p]$.

Given a point $b\in B$, we compute ${\sigma }_b$ as follows. Because of the ordering property and the definition of $a$, all sites of $A_2\setminus \{a\}$ lie in leaves strictly to the right of the leaf of $a$. In $T_A$, we find a set $V_b$ of $O(\log n)$ nodes whose leaf ranges form a partition of the leaves to the right of $a$ (including the leaf of $a$). For each $v\in V_b$, using ${𝒟}_v$, we compute a weighted nearest neighbor $p_v$ of $b$ in $A_v$ and treat each $p_v$ as a candidate for ${\sigma }_b$. We then return the candidate minimizing $\mathrm{dist}[p_v]+d(p_v,b)$. Thus, computing ${\sigma }_b$ requires $O(\log n)$ weighted nearest neighbor queries plus $O(\log n)$ geodesic distance queries. Recall that a GH data structure for $P$ can be built in $O(m)$ time, after which each distance query $d(p_v,b)$ is computable in $O(\log m)$ time.

Next we address updates. At the end of each iteration, for every $b\in B_a$, after computing $\mathrm{dist}[b]$, we insert $b$ into $T_A$. This places $b$ at a leaf, and for each node $v$ on the path to the root, we insert $b$ into ${𝒟}_v$. By Lemma 2, the ordering property of $T_A$ is preserved automatically. Therefore, there is no need to delete $a$ from $T_A$, and each ${𝒟}_v$ need only support insertions.

To build an insertion-only data structure ${𝒟}_v$, one might attempt Bentley’s logarithmic method [3] together with a static geodesic Voronoi diagram for $A_v$. Since computing the geodesic Voronoi diagram of $A_v$ in $P$ requires $\mathrm{\Omega }(m)$ time (note that the best algorithm can compute the diagram in $O(m+|A_v|\log |A_v|)$ time [24]), rebuilding these diagrams $n$ times would lead to $\mathrm{\Omega }(mn)$ total time, which is too slow. Instead, we require an implicit construction whose dependence on $m$ is only polylogarithmic. This is exactly the subproblem introduced in Section 1.1. In the full paper, we establish the following result.

With Theorem 3, we construct ${𝒟}_v$ for every node $v$ of $T_A$. Using Bentley’s method, each insertion into ${𝒟}_v$ takes amortized $O(\log n{\log }^{3}m+{\log }^{2}n{\log }^{2}m)$ time. Each weighted nearest neighbor query takes $O({\log }^{2}n{\log }^{2}m)$ time. Since $T_A$ has height $O(\log n)$, computing ${\sigma }_b$ takes $O({\log }^{3}n{\log }^{2}m)$ time, and inserting a site $b$ into $T_A$ takes amortized $O({\log }^{2}n{\log }^{3}m+{\log }^{3}n{\log }^{2}m)$ time. Thus, the total running time of Subroutine II over the entire algorithm is $O(m+n{\log }^{2}n{\log }^{3}m+n{\log }^{3}n{\log }^{2}m)$, including the $O(m)$ time required to build the GH data structure for $P$ in Theorem 3. It is actually bounded by $O(m+n{\log }^{3}n{\log }^{2}m)$.

Roughly speaking, the right envelope $\mathrm{\Xi }(S)$ of ${𝒟}_R(S)$ is the boundary of ${𝒟}_R(S)$ excluding $d$. Specifically, $\mathrm{\Xi }(S)$ is union of ${\mathrm{\Xi }}^{\circ }(S)$ and ${\mathrm{\Xi }}^{\diamond }(S)$, with ${\mathrm{\Xi }}^{\circ }(S)=\{p\in P_R:d(S,p)=1\}$ and . In general, ${\mathrm{\Xi }}^{\circ }(S)$ is a collection of circular arcs each belonging to $\partial D_s$ for some point $s\in S$, while ${\mathrm{\Xi }}^{\diamond }(S)$ is a collection of segments on $\partial P_R$; see Figure 3(a).

Note that ${𝒟}_R(S)$ may have multiple connected components. We can show that for each connected component $𝒟$ of ${𝒟}_R(S)$, $𝒟\cap d$ consists of a single segment of $d$ and $𝒟$ is simply connected (i.e., there are no holes), meaning that $\partial 𝒟$ is composed of a segment of $d$ and a simple curve connecting the two endpoints of $𝒟\cap d$.

For any point $s\in S$, define ${\xi }_s(S)$ as the set of points $p\in \mathrm{\Xi }(S)$ such that $d(s,p)\le d(s^{\prime },p)$ for any other point $s^{\prime }\in S$. In other words, if we consider the geodesic Voronoi diagram $𝒱𝒟(S)$ of $S$, then ${\xi }_s(S)$ consists of all points $p\in \mathrm{\Xi }(S)$ in the Voronoi region of $s$; see Figure 3(b). A useful property is that if ${\xi }_s(S)\ne \mathrm{\varnothing }$, then ${\xi }_s(S)$ is a single connected portion of $\mathrm{\Xi }(S)$.

Intuitively, for any point $q\in P_R$, $q\in {𝒟}_R(S)$ if and only if $q$ is “to the left of” $\mathrm{\Xi }(S)$. However, the question is how to interpret “left.”

In the Euclidean case where $P$ is the entire plane and $d$ becomes a vertical line, $\mathrm{\Xi }(S)$ is a $y$-monotone curve. Due to this property, there is a trivial way to determine whether $q$ is to the left of $\mathrm{\Xi }(S)$: First find by binary search the arc of $\mathrm{\Xi }(S)$ intersecting the horizontal line through $q$, and then determine whether $q$ is to the left of the arc [14, 30].

In our problem, however, $\mathrm{\Xi }(S)$ is no longer $y$-monotone. Intuitively, this is because $P$ is not “large enough” and $\partial P$ changes the “shape” of $\mathrm{\Xi }(S)$. Thus, we need a new approach. At a very high level, our strategy is to develop some concepts analogous to those in the Euclidean case. Specifically, we introduce an order ${\prec }_R$ for any two points of $P_R$; this order corresponds to the $y$-coordinate order in the Euclidean case. We also define a traversal order ${\prec }_{\mathrm{\Xi }}$ for any two points of $\mathrm{\Xi }(S)$; in the Euclidean case, this would be the order of traversing on $\mathrm{\Xi }(S)$ from bottom to top as $\mathrm{\Xi }(S)$ is $y$-monotone. In addition, we define an order ${\prec }_L$ for all points of $S$ based on their geodesic distances to the lower endpoint $z_0^{\ast }$ of $d$ (i.e., for any two points $s,t\in S$, define $s{\prec }_{L}t$ if ). A useful property is that these three orders are consistent with each other, which is crucial to the success of our approach.

Our data structure $\mathrm{\Psi }(S)$ is a complete binary tree whose leaves from left to right store the points of $S$ under the order ${\prec }_L$. For each node $v$, let $S_v$ be the subset of points $S$ stored in the leaves of the subtree rooted at $v$. We wish to store an implicit representation of $\mathrm{\Xi }(S_v)$ at $v$. Let $I(S_v)$ be the set of points $s\in S_v$ such that ${\xi }_s(S_v)\ne \mathrm{\varnothing }$. We consider $I(S_v)$ a list of its points sorted under ${\prec }_L$.

If $v$ is the root, then we store $I(S_v)$ at $v$ explicitly. If we store $I(S_v)$ at $v$ explicitly for all nodes $v\in \mathrm{\Psi }(S)$, then $\mathrm{\Psi }(S)$ would take super-linear space. We instead employ a standard space saving technique [18, 26], where we store a point $s\in S$ in the highest node $v$ such that $s\in I(S_v)$. Using this strategy, the total space of $\mathrm{\Psi }(S)$ is $O(n)$.

For a node $v$ with left child $u$ and right child $w$, define $r_u={\max }_{{\prec }_L}I(S_v)\cap I(S_u)$, i.e., among all points in $I(S_v)\cap I(S_u)$, $r_u$ is the largest one following the order ${\prec }_L$. Symmetrically, define $r_w={\min }_{{\prec }_L}I(S_v)\cap I(S_w)$. We can show that $I(S_u)\cap I(S_v)$ is a prefix of $I(S_u)$ and $I(S_w)\cap I(S_v)$ is a suffix of $I(S_w)$. Hence, $I(S_v)$ is the concatenation of $I{(S_u)}_{\le r_u}$ and $I{(S_w)}_{\ge r_w}$, where the notation $I{(S_u)}_{\le r_u}$ represents the prefix of $I(S_u)$ up to and including $r_u$ and $I{(S_w)}_{\ge r_w}$ represents the suffix of $I(S_w)$ from and including $r_w$. We store both $r_u$ and $r_w$ at $v$. As $I(S_v)$ is the concatenation of $I{(S_u)}_{\le r_u}$ and $I{(S_w)}_{\ge r_w}$, we only store $I{(S_u)}_{>r_u}$ in a list $arcs(u)$ at $u$ and similarly store in $arcs(w)$ at $w$. This means that in a descent from the root of $\mathrm{\Psi }(S)$, we can do the following. Let $v$ be the current node in our descent and assume $arcs(v)=I(S_v)$, which is true initially when $v$ is the root. Suppose $v$ has left child $u$ and right child $w$. We can cut $arcs(v)$ at $r_u$ to obtain $I{(S_u)}_{\le r_u}$ and $I{(S_w)}_{\ge r_w}$. As $I{(S_u)}_{>r_u}$ is stored in $arcs(u)$, we can concatenate $I{(S_u)}_{\le r_u}$ and $I{(S_u)}_{>r_u}$ to obtain $I(S_u)$. We can do likewise to obtain $I(S_w)$. In an ascent, we can simply reverse this process to obtain $I(S_v)$ from $I(S_u)$ and $I(S_w)$ using $r_u$ and $r_w$. To facilitate the queries, additional information will also be maintained such as the endpoints of ${\xi }_s(S_v)$ for each point $s\in arcs(v)$.

It is not obvious why the algorithm is correct. See the full paper for the proof.