Abstract 1 Introduction 2 BFS in ball-intersection graphs 3 Reverse shortest paths in unit-ball graphs in three dimensions 4 Reverse shortest paths in ball-intersection graphs in 𝟑 5 Further extensions References

BFS and Reverse Shortest Paths for Ball Intersection Graphs in Three and Higher Dimensions

Matthew J. Katz ORCID Department of Computer Science, Ben-Gurion University, Beer Sheva, Israel Rachel Saban ORCID Department of Computer Science, Ben-Gurion University, Beer Sheva, Israel Micha Sharir ORCID School of Computer Science and AI, Tel Aviv University, Israel
Abstract

Let be a collection of n arbitrary balls in 3, and let G0() be their intersection graph. We provide an algorithm for performing BFS on G0(), which runs in O(n4/3) time, where the O() notation hides subpolynomial factors. For r0, let Gr() be the intersection graph of the set r={B+rB}, where B+r is the ball concentric with B whose radius is larger by r than the radius of B. We provide an efficient algorithm for the reverse shortest path (RSP) problem, where we are given two designated balls Bs, Bt of and a parameter 0<λ<n, and seek the smallest value r for which Gr() contains a path from Bs to Bt of at most λ edges. For the special case of congruent balls (equivalently, for points in 3), the algorithm runs in O(n29/21)O(n1.381) time. For the general case, the algorithm runs in O(n56/39)O(n1.436) time. We also extend the technique to handle other measures of expansion and higher dimensions.

Keywords and phrases:
Computational geometry, reverse shortest paths, breadth-first search, shrink-and-bifurcate, intersection graphs
Funding:
Matthew J. Katz: Work partially supported by Grant 2019715/CCF-20-08551 from the US-Israel Binational Science Foundation/US National Science Foundation, and by Grant 495/23 from the Israel Science Foundation.
Rachel Saban: Work partially supported by Grant 495/23 from the Israel Science Foundation.
Micha Sharir: Work partially supported by Grant 495/23 from the Israel Science Foundation.
Copyright and License:
[Uncaptioned image] © Matthew J. Katz, Rachel Saban, and Micha Sharir; licensed under Creative Commons License CC-BY 4.0
2012 ACM Subject Classification:
Theory of computation Computational geometry
Editors:
Ho-Lin Chen, Wing-Kai Hon, and Meng-Tsung Tsai

1 Introduction

The reverse shortest path (RSP) problem has received considerable attention recently. In general, we are given a set P of n geometric objects, and a parameter r0. We define a graph Gr(P) on P, whose edges are all the pairs (u,v) that satisfy some property, given by a predicate Π(u,v;r), which is monotone in r, meaning that if Π(u,v;r) holds then Π(u,v;r) holds for all r>r. That is, the graphs Gr(P) are monotone increasing in r. We are given two designated objects s,tP and a parameter 0<λ<n, and wish to find the smallest value r for which Gr(P) contains a path from s to t with at most λ edges.

Among the simplest examples is the case where P is a set of n points in the plane, and Gr(P)={(u,v)|uv|r}, under the Euclidean distance |uv|. This can also be interpreted as the intersection graph of the disks of radius r/2 centered at the points of P (the so-called unit-disk graph, with the unit being r/2). This case and many variants thereof are reviewed later in the introduction.

One generalization, studied here, is to the case of points in three dimensions. That is, Gr(P) in this case is the intersection graph of n congruent balls of radius r/2 centered at the points of P. This is a special instance of the general problem studied in this work, where the balls have arbitrary radii. That is, the input consists of a set of n balls in 3 of arbitrary radii, and r is a parameter that measures how the balls of are expanded. The simplest case is where we increase the radius of each ball by adding r to it. Let Gr() denote the intersection graph of the expanded balls, so G0() is the intersection graph of the original balls. Our technique also applies to any other well-behaved measure of expansion, e.g., multiplying each radius by r1, but for conciseness we will mainly stick to the additive measure of expansion, as just introduced, and discuss the general setup later in the paper.

A standard high-level approach to the RSP problem is to design a decision procedure which, with r as input, determines whether Gr() contains a path from s to t of length at most λ. If it does then rr, and if not then r>r. The decision procedure is typically implemented using Breadth-First Search (BFS) on Gr(). Once we have such a procedure, we can use binary search on r to home in on the correct value of r.

There are two major problems with this approach. First, the cost of a naïve implementation of BFS is linear in the number of edges of Gr(), which can be Θ(n2) in the worst case. Second, BFS appears not to be parallelizable, and thus it is not amenable to parametric search, which is the standard technique for turning the decision procedure into an efficient optimization procedure for solving the RSP problem. In this work we overcome these issues, for the special case of ball-intersection graphs in 3 (and in higher dimensions), and obtain efficient algorithms for solving the problem, albeit they are (not surprisingly) less efficient than in the planar case.

Related work.

Breadth-First Search (BFS) is recognized as one of the fundamental graph algorithms. On general graphs G=(V,E), its running time is O(|V|+|E|), which is inefficient when |E| is large. Its importance has led authors to seek families of graphs in which BFS can be implemented more efficiently. In the geometric context, one such example is the family of disk graphs in the plane. The vertices of a disk graph represent disks in the plane and there exists an edge between two vertices if and only if the corresponding disks intersect. Notice that the number of edges in a disk graph can be quadratic in n, the number of disks. Nevertheless, for unit-disk graphs (where all the disks are congruent), Cabello and Jejčič [8] presented an O(nlogn) implementation of BFS, and subsequently Chan and Skrepetos [10] presented an alternative O(n) implementation (after pre-sorting the points by their x- and y-coordinates). For arbitrary disks, Klost [14] described an O(nlog2n) implementation of BFS (see also [11, 12]).

The proximity graph of a set S of n segments in the plane, for a real parameter r0, is Gr(S)(S,E), where E={(e1,e2)dist(e1,e2)r} and dist(e1,e2) is the Euclidean distance between e1 and e2 (which is 0 if they intersect). Agarwal et al. [5] devised an O(n4/3) implementation of BFS in Gr(S). (As in the abstract, the O() notation hides subpolynomial factors.) For the special case where the segments in S are pairwise disjoint, Agarwal et al. [4] have later provided an O(nlog2n) implementation.

The bottleneck path problem (also known as the minimax path problem) and its complementary problem, i.e., the widest path problem (or the maximum capacity problem), are well-known problems in graph theory. In suitable contexts, the bottleneck path problem is strongly related to the RSP problem: Given s and t, if π is a path from s to t of at most λ edges with the minimum bottleneck r, then r is the solution to the RSP problem with s, t, λ as input, and vice versa.

Abu-Affash et al. [1] studied the problem of placing at most k Steiner points to minimize the bottleneck of a path between two designated points in the plane, for a given integer parameter k0, and developed an O(nlog2n)-time algorithm for the problem.

The RSP problem for unit disks (as mentioned above) was studied by Wang and Zhao [16], who proposed an O(n5/4)-time solution. An improved solution with running time O(n6/5), using the shrink-and-bifurcate technique (briefly reviewed later; see [7]), was subsequently presented by Kaplan et al. [11]. A very recent improvement of this technique by Chan and Huang [9] yields an implementation with running time O(n8/7), which can be further improved, for the case of unit-disk graphs, to O(n9/8) time. For arbitrary disks, Kaplan et al. [11] obtained a solution with randomized expected time O(n5/4), which can be improved to randomized expected time O(n6/5), using Chan and Huang’s technique [9].

The RSP problem was also studied for other objects, including wedges of some fixed angle, which are viewed as directional antennas, and unit-height towers placed on a 1.5-dimensional terrain. The former version was solved in O(n4/3) time by Agarwal et al. [5], and the latter version was studied in Katz et al. [13] (see also [9]).

Our results.

Our main contributions are efficient algorithms for the BFS and RSP problems in ball-intersection graphs in three and higher dimensions. We first consider the problem of efficient implementation of BFS, as a decision procedure for guiding the search for the optimum r in the RSP problem (and also an interesting problem in itself). In three dimensions we can perform BFS on G0(), for a set of n (congruent or arbitrary) balls in 3, in O(n4/3) time, and in d dimensions in O(n2e/(e+1)) time, where e=d/2+1. As these algorithms are not parallelizable (at least we do not know how to run them in small parallel depth),111This is a standard issue in BFS on any graph. we need to use a different approach. Such an approach, already mentioned above, known as the shrink-and-bifurcate technique, has originally been developed in Ben Avraham et al. [7] for computing a rather special case of Fréchet distances. But it has recently found applications to the RSP problem in the plane [11]. Very recently, as already mentioned, it has been dramatically improved by Chan and Huang [9].

Once a BFS procedure (serving as the decision procedure) is available, our technique follows closely the improved technique in [9], which however requires certain nontrivial modifications and enhancements to fit into higher-dimensional contexts. The technique of [9] consists of two parts, a general-purpose part, and an additional enhancement which improves the running time further for the special case of unit-disk graphs. We follow the general-purpose part for arbitrary balls, and both parts for congruent balls.

In three dimensions, for the special case of congruent balls (equivalently, of a set of n points in 3), we obtain an algorithm that runs in O(n29/21) randomized expected time. For the case of balls of arbitrary radii, we obtain an algorithm that runs in O(n56/39) randomized expected time. More generally, in d3 dimensions, the algorithm for arbitrary balls runs in

O(n2(d+1)(3e+1)(3d+4)(e+1))

randomized expected time, and, for the special case of congruent balls, the algorithm runs in randomized expected time

O(ng2e/(e+1)+g+2ee+1), where g=2(de)(e+1)(3d+1).

2 BFS in ball-intersection graphs

Let be a collection of n balls in 3, of arbitrary radii.222The case of congruent balls does not lead to an improved BFS implementation, although the algorithm becomes considerably simpler; see a discussion at the end of the section. We denote a ball with center c and radius ρ as B(c,ρ). Let G0() be the intersection graph of , i.e., the vertices of G0() are the balls of , and its edges consist of all the intersecting pairs of balls. The condition for two balls B(c1,ρ1), B(c2,ρ2) to intersect is |c1c2|ρ1+ρ2.

Let Bs be a start ball in . Our goal is to run BFS on G0() starting from Bs. We follow the standard approach, of constructing the layers of the BFS in order, starting from layer L0 that contains only Bs. Consider the step of passing from some layer Li to the next layer Li+1. Let 𝒰i denote the set of all balls that the BFS has not yet reached up to this step. We then need to find, for each ball B=B(c,ρ) in Li, all the balls of 𝒰i that it intersects. Each such ball is added to Li+1 and is immediately deleted from 𝒰i, to ensure that it is not detected again by other balls of Li (or by balls in future layers). We continue in this way until no more balls of 𝒰i intersect any ball in Li. At this point Li+1 has been computed, and we go on to construct the next layer. (For our decision procedure, executed on Gr(), we stop when either the target ball Bt has been reached or layer Lλ has been constructed, whichever happens first. If Bt has been reached, we report success (i.e., rr), otherwise we report failure (i.e., r>r).)

To implement this algorithm efficiently, we therefore need a dynamic data structure for storing 𝒰=𝒰i, the set of unreached balls, that supports (efficiently) the following two kinds of operations.

Intersection detection queries: Given a query ball B, detect whether it intersects any ball of 𝒰, and, if so, report such a ball.

Deletions: Delete a ball B from 𝒰.

An intersection query with a ball B(c0,ρ0) can be rephrased as a range emptiness detection query. For this, we map each ball B=B(c,ρ)=B((xc,yc,zc),ρ) of 𝒰 to the point B^=(xc,yc,zc,ρ) in 4, and map each ball B0=B(c0,ρ0)=B((xc0,yc0,zc0),ρ0) of Li to the range

σB0={(x,y,z,ρ)4(xxc0)2+(yyc0)2+(zzc0)2(ρ+ρ0)2}.

Then a ball B intersects B0 iff B^σB0.

Using a standard lifting transform, we further map each point B^=(xc,yc,zc,ρ) to the point B=(xc,yc,zc,ρ,xc2+yc2+zc2ρ2) in 5. Each range σB0 is lifted to the halfspace

hB0:=x52xc0x+2yc0y+2zc0z+2ρ0ρ+(ρ02xc02yc02zc02).

As is easily checked, B^σB0 iff B lies in hB0.

In other words, in the transformed problem we have a set ={BB} of n points in 5, which we want to process into a dynamic data structure (under deletions) for answering halfspace range emptiness queries. The algorithm of Agarwal and Matoušek [6] provides such a structure. Specifically, for any given parameter nsn2, the structure can be constructed to have size O(s), an initial version of it can be constructed in O(s) time, each query takes O(n/s1/2) time, and each deletion takes amortized O(s/n) time. Choosing s=n4/3, we obtain a data structure of size O(n4/3), so that each operation (query or deletion) on the structure takes O(n1/3) time.

The number of operations is O(n), since each query either detects a new intersecting ball, which is promptly removed from 𝒰, or is the last query performed with a ball, and each ball becomes a query only at the layer it belongs to. We conclude that the overall running time of the BFS is O(n4/3), because the cost of maintaining and searching in the data structure dominates the cost of the other steps of the BFS. We thus obtain:

Theorem 1.

BFS on the ball intersection graph of a set of n balls in 3 can be performed in O(n4/3) time.

The case of congruent balls.

We first state the following theorem, whose results will be useful for handling congruent balls. Its proof is an easy consequence of the analysis in [6, 15], applied to the lifted points and halfspaces in 5, as above.

Theorem 2.

(a) Given two sets P, Q of m and n points, respectively, in three dimensions, and a prameter r>0, we can determine, for each pP whether there exists a point qQ such that |pq|r, in overall time

O(m2/3n2/3+m+n).

(b) Given two sets P, Q, as above, we can compute, for each pP, the nearest neighbor of p in Q, in overall time

O(m2/3n2/3+m+n).

(Note that (b) subsumes (a), but the algorithm in (a) is simpler, so we state it separately.)

When the balls of are congruent, say all of radius r, we do not need the above dynamic (and rather involved) structure. Instead, we follow the approach in the planar setup of unit-disk graphs [8, 10]. That is, we form a grid of cell size r/2, and distribute the ball centers among its cells. When we perform the step of passing from Li to Li+1, we process each active grid cell, i.e., a cell that contains centers of balls in Li. For each such cell τ, all unreached balls with centers in τ are immediately placed in Li+1. Then, for each cell τ in a suitable punctured constant-size neighborhood of τ (only such cells, and τ itself, can contain centers of balls in 𝒰i that can be reached, in one step, from balls with centers in τ), we take the set 𝒰i(τ) of unreached balls with centers in τ and test each such ball B whether its center lies within distance r from the center of some ball in Li(τ), the set of balls of Li with centers in τ. Theorem 2(a) implies that this reverse search can be implemented in time

O(|𝒰i(τ)|2/3|Li(τ)|2/3+|𝒰i(τ)|+|Li(τ)|).

The remaining straightforward implementation details, and the analysis of correctness and performance of the structure, as given in [8, 10, 11], carry over to the three-dimensional case, and imply that the overall cost of this implementation of the BFS is O(n4/3) time.

3 Reverse shortest paths in unit-ball graphs in three dimensions

We can now solve the corresponding optimization problem, which is the reverse shortest path (RSP) problem, or, as it is also called, the bounded-hop bottleneck path problem, for unit-ball graphs in 3: Given a set P of n points in 3, two designated points s,tP, and a parameter 0<λ<n, find the minimum value r of r for which the associated graph Gr(P) contains a path from s to t of at most λ edges.

We can solve the RSP problem using the aforementioned shrink-and-bifurcate technique [7, 9, 11]. We begin with a brief overview of the technique. It applies in more general setups, but let us restrict it to our setup of unit-ball graphs (and later also to intersection graphs of arbitrary balls, although the description below assumes that the balls are congruent). The shrinking stage constructs an interval I that contains r and at most L other critical values of r, each of which is the distance between two centers, where L is a prespecified parameter. It does so, in the improved approach of [9], by running a distance selection algorithm on various pairs of random samples of (that is, on the sets of centers of the balls in the samples).

The distance selection algorithm, on which the shrinking procedure is based, takes, ignoring the shrinking aspect, O(n3/2) time. More generally, it takes O(m3/4n3/4+m+n) time for selecting distances between two sets of m and n points, respectively, which is the setup that arises when applying the technique of [9]. Using (a suitable adaptation to three dimensions of) the shrinking mechanism of Chan and Huang [9], the construction of the desired interval I can be performed in expected O(n3/2/L3/4) time.

This is followed by a bifurcation stage, in which we simulate the decision procedure (in our setup, the BFS procedure of the preceding section). When we reach a comparison, we resolve it right away if the critical value r that it induces lies outside I. If however r lies in I, we bifurcate, exploring both outcomes rr and r>r. When we have either accumulated enough unresolved comparisons, or moved forward at least some number of steps in the simulation, along each path in the bifurcation tree, we stop and resolve all comparisons using binary search, guided by the (unsimulated) decision procedure. This shrinks I further, and we start a new phase of bifurcating simulation. As shown in [7, 11], with a suitable choice of parameters, the overall cost of the bifurcation is O(L1/2D(n)), where D(n)=O(n4/3) is the cost of the decision procedure (i.e., of the BFS).

Hence the overall cost of the RSP algorithm is O(n3/2/L3/4+L1/2n4/3). We choose L to balance these terms, that is, L=n2/15, and then the running time is O(n7/5). That is, we have our initial, albeit weaker, result:

Theorem 3.

For a set P of n points in 3, two designated points s,tP, and an integer parameter 0<λ<n, we can find the minimum radius r for which the associated graph Gr(P) contains a path from s to t of at most λ edges, in O(n7/5) time.

3.1 An improved solution

To obtain a faster algorithm, we follow the high-level approach of the second improvement in [9], with many nontrivial adjustments that cater to the three-dimensional nature of the problem. We follow the same setup as in the second implementation of the BFS algorithm. However, for technical reasons that will become clear later on, we first perform the following steps, having to do with the cell size of the grid. We first note that the precise size of the cells (which we have set to r/2 in the BFS algorithm) is not very important, as long as it is smaller than r/3 (but not much smaller), so that the propagation within a cell is immediate. We can therefore proceed as follows. Note that, by the nature of the RSP problem, r must be between |st|/2 and |st|/(2n) (or, more precisely, |st|/(2λ)). We therefore run binary search through this interval, using any of the BFS decision procedures of Section 2 to guide the search, and in O(logn) steps we locate r within a range [r1,(1+ε)r1], for some suitable small ε>0. We then set r to be333Note that this choice might force us to change the definition of the neighborhood of a cell to potentially contain more cells, but still only a constant number thereof. r1/2, and keep this size throughout the simulation. See below how this is used in the algorithm. We emphasize that this step precedes the actual simulation, as well as the step of distributing the points of P in the grid cells.

We introduce a new threshold parameter Δn, whose concrete value will be set later, and classify each nonempty grid cell as either heavy, if it contains at least Δ points of P, or light, otherwise. The number of heavy cells is at most k:=n/Δ. We now compress each heavy cell τ into a single (symbolic) “heavy” point pτ, and obtain a new compressed version of Gr(P), which we denote as Hr(P^), where P^ consists of all the points in the light cells and of all the representative heavy points pτ. The distance between any pair p,q of points in light cells remains the same, namely |pq|. The distance between a heavy point pτ and a light point q (a point in a light cell) is dist(q,Pτ)=minpτ|qp|. Finally, the distance between two heavy points pτ, pτ is the distance determined by the bichromatic closest pair (BCP) in Pτ×Pτ.

Running BFS on Hr(P^) is performed similarly to the previous implementation. A single propagation step, at some iteration i, involving two adjacent cells τ, τ, is executed depending on whether one of τ, τ, or both, are heavy. If both cells are light, we proceed as before, incurring the cost

O(|Ui(τ)|2/3|Li(τ)|2/3+|Ui(τ)|+|Li(τ)|); (1)

see Theorem 2(a). When τ is light and τ is heavy, Ui(τ) (if nonempty) is the singleton representative point pτ, and by running the all-nearest-neighbor procedure (Theorem 2(b)) on Li(τ) and Pτ, we can determine whether this step should place pτ in Li+1. The cost is then

O(|Pτ|2/3|Li(τ)|2/3+|Pτ|+|Li(τ)|). (2)

A similar procedure is applied when τ is heavy and τ is light. In this case each point of Ui(τ) searches for a near neighbor in Pτ, and the cost is

O(|Ui(τ)|2/3|Pτ|2/3+|Ui(τ)|+|Pτ|). (3)

Finally, when both cells are heavy we need to compute the distance of the BCP in Pτ×Pτ. This step is handled differently because its naïve implementation is too expensive, because it does not exploit the smaller parameter Δ; see shortly below.

We sum the bounds in (1), (2) and (3), over all iterations i and pairs τ, τ of adjacent cells that are active at the corresponding iteration, when both τ and τ are light (eq. (1)), when τ is light and τ is heavy (3), and when τ is heavy and τ is light (2). Note that in each of the bounds (1)–(3), at least one of the sets (Li(τ) and / or Ui(τ)) is from a light cell, so its size is at most Δ.

Applying Hölder’s inequality, one can show that the sums of the bounds in (1), (2) and (3) are all O(nΔ1/3). For example, summing the bound in (1) over all iterations i and active pairs τ,τ of neighboring cells, and using the facts that (i) |Li(τ)|Δ, and (ii) each cell is active, or a neighbor of an active cell, at only O(1) iterations, in each of these bounds, we get

O(iτ,τ(|Ui(τ)|2/3|Li(τ)|2/3+|Ui(τ)|+|Li(τ)|))
=O(Δ1/3iτ,τ|Ui(τ)|2/3|Li(τ)|1/3)+O(n)
=O(Δ1/3(iτ,τ|Ui(τ)|)2/3(iτ,τ|Li(τ)|)1/3)+O(n)=O(nΔ1/3),

and similarly for the sums of the bounds in (2) and (3).

It remains to handle the case where both τ and τ are heavy. Here we need to compute the distance of the BCP in Pτ×Pτ, over all adjacent heavy pairs (i.e., pairs where both τ and τ are heavy). As above, this can be done in a total of O(n4/3) time. Since it does not depend on Δ, this cost turns out to be too expensive to pay during the simulation of the procedure. To bypass this issue, we note that, once we have fixed the cell-size r1/2 of the grid (as we did, ahead of the simulation), this step is independent of r or r, and can too be performed ahead of the simulation. We thus compute the BCPs for all pairs of adjacent heavy cells, once and for all. The overall cost, of both the binary search to fix the cell-size of the grid and the BCP computations, is easily seen to be O(n4/3). In addition, after having computed all the heavy BCP distances, we run an additional binary search through them, to ensure that the interval with which we start the shrink-and-bifurcate procedure does not contain any heavy BCP distance in its interior.

All these considerations imply that we can run BFS on the compressed graph Hr(P^) in time O(nΔ1/3), after an initial phase that takes O(n4/3) time, which we will execute only once, outside (and preceding) the actual simulation during the search for r.

We now run this simulation of the BFS on Hr(P^), using the shrink-and-bifurcate paradigm, as in [9] and as briefly reviewed above. As before, the shrinking procedure is based on distance selection in 3, but we apply it separately to each pair of adjacent grid cells, at least one of which is light. (This is justified as in [9], because r must arise as a distance within such a pair, unless it is the BCP distance of some heavy pair, which will have been detected during the preliminary stage.) Arguing as in [9], and adapting the setup to 3, the cost of the shrinking stage is O(nΔ1/2/L3/4).

The subsequent bifurcation procedure runs, as before, in time O(L1/2D(n)), where D(n) is the cost of the decision procedure. Ignoring the cost of the preliminary part, which stays outside the simulation, we have D(n)=O(nΔ1/3). Hence the overall cost of the shrink-and-bifurcate procedure, adding back the preliminary cost as a one-time cost, is

O(nΔ1/2/L3/4+L1/2nΔ1/3+n4/3). (4)

Of course, so far we have only simulated the BFS on the compressed graph Hr(P^), so the output value r is still not the correct value. Denoting it by r1, it is the smallest value of r for which Hr(P^) has a path of at most λ edges between s and t. As argued in [9], denoting by dH(s,p) the graph distance between s and a point p in Hr1(P^), and by dG(s,p) the graph distance between s and p in Gr(P), we have, for each p,

dG(s,p)kdH(s,p)dG(s,p), (5)

where k=n/Δ. With some care. this also covers the cases where s and / or t lie in heavy cells.

Recovering 𝒓

As in [9], to find the true value r, we use dynamic programming, following the high-level approach of [9]. For the sake of completeness, we repeat here some details of this technique, adapted to our setup. We first comment that it is easy to solve the RSP problem using dynamic programming, without using the preceding machinery at all. To do so, let Di(p), for pP, denote the smallest bottleneck value of a path from s to p consisting of at most i edges. Initially, D0(s)=0 and D0(p)=+ for any ps. Then we have, for i1 and pP,

Di(p)=minqPmax{Di1(q),|qp|}. (6)

If Dλ(t)<+ then Dλ(t)=r, which follows easily from (6). The problem with this approach is that the number of times the step (6) is applied is Θ(n2) (it is actually O(nλ), for the number of points n and the number of iterations λ to reach, or not to reach, t). Moreover, without any efficient implementation of the steps (6), the overall cost can be Θ(n3).

In our setup, though, we can reduce the number of dynamic programming steps to O(nk). Indeed, write Pi={pPikdH(p)i}. As each point participates in k+1 sets Pi, we have i|Pi|=O(nk).

As in [9], we exploit this property by replacing the function D by another function D^, defined by the following modified iterative process. For each pPi put

D^i(p)=minqPi1max{D^i1(q),|qp|}, (7)

and D^0(s)=0 and D^0(p)=+ for all ps. As shown in [9], we also have D^λ(t)=r. (This argument is fully general and discrete, and does not depend on any geometric feature of the problem, except for the inequalities (5). In particular, it does not depend on the dimension.)

The number of steps in this modified dynamic programming process is thus O(nk). To make these steps (relatively) efficient, we follow a variant of the technique in [9, Lemma 8], adapted to three dimensions. Briefly, the input is a set Q of n points in 3, so that each point qQ carries a positive weight wq (the weight wq will be set to D^i1(q) at the i-th iteration), and we have a set P of m queries, where each query point p seeks the point qQ that minimizes max{wq,|qp|}. The algorithm in [9, Lemma 8] is based on Voronoi diagrams, which works well in the plane, but is not efficient in three dimensions. We use the following modified approach. As in [9], we store the points of Q at the leaves of a balanced binary tree T, sorted in increasing order of their weights. Each node v of T processes the set Qv of those points of Q that are stored at the subtree rooted at v. Each query point p follows a path in T. When p reaches a node v, with a left child vL and a right child vR, it finds its (standard) nearest neighbor qL in QvL. If |pqL|wqL we recurse in vL, and if |pqL|>wqL we recurse in vR. See [9] for details and justification of this rule.

Instead of using Voronoi diagrams, as in [9], we use the algorithm in Theorem 2(b). To do so, we run all the queries in parallel, proceeding level by level. When we reach a node v, we have available the set Pv of all the queries whose paths through T reach v. Using the algorithm in Theorem 2(b), we find, for each pPv, its nearest neighbor in QvL, in overall time

O(|Pv|2/3|QvL|2/3+|Pv|+|QvL|).

This allows us, by our search rule, to assign each pPv either to vL or to vR, as appropriate. Carrying out this process at each node of the current level, we obtain, for each node w of the next level, the set Pw of query points that want to access w (actually, wL). We continue in this recursive manner until we reach the leaves. The cost at a leaf v is O(1+|Pv|). Summing over all nodes of T, bearing in mind that

v|Pv|=O(mlogn)andv|QvL|=O(nlogn),

the overall cost is easily seen to be O(|P|2/3|Q|2/3+|P|+|Q|).

We now apply this procedure at each iteration i of the dynamic programming, where the input set is Pi1, the query set is Pi, and the weight of each input point q is D^i1(q). Clearly, the output of this procedure, with this input, implements the iterative process in (7). The cost of iteration i is therefore

O(|Pi|2/3|Pi1|2/3+|Pi|+|Pi1|).

Since each |Pi| is at most n, and i|Pi|=O(nk), the overall cost of the dynamic program is (using Hölder’s inequality)

O(i(|Pi|2/3|Pi1|2/3+|Pi|+|Pi1|))
= n1/3O(i(|Pi|1/3|Pi1|2/3))+O(i(|Pi|+|Pi1|))
= O(n1/3(nk)1/3(nk)2/3+nk)=O(n4/3k)=O(n7/3Δ). (8)

Altogether, the cost of the full algorithm is therefore

O(nΔ1/2L3/4+L1/2nΔ1/3+n4/3+n7/3Δ)=O(nΔ1/2L3/4+L1/2nΔ1/3+n7/3Δ). (9)

We first balance the first two terms, choosing L5/4=Δ1/6, or L=Δ2/15, and then the bound becomes

O(nΔ2/5+n7/3Δ).

We now choose Δ to balance these terms, that is Δ=n20/21, and the cost of the algorithm is therefore O(n29/21). That is, we have

Theorem 4.

Let P be a set of n points in 3, let s and t be two designated points, and let 0<λ<n be a given parameter. Then we can find the smallest value r for which Gr(P) contains a path between s and t that consists of at most λ edges, in O(n29/21)=O(n1.381) time.

3.2 Higher dimensions

The algorithm given above can be generalized to points in d, for any dimension d4. The geometric components in the algorithm are (i) constructing the grid, (ii) distance selection and the corresponding interval shrinking procedure, (iii) an efficient algorithm for the all-near-neighbor problem, (iv) an efficient algorithm for the all-nearest-neighbor problem, and, as a special case, (v) an efficient algorithm for the BCP problem.

The grid construction can trivially be extended to any fixed dimension. Distance selection in d takes O(n2d/(d+1)) time (see, e.g., [3, Appendix]). The bipartite version, where we want to select distances between two given sets of m and n points, respectively, takes O(md/(d+1)nd/(d+1)+m+n) time, and the corresponding interval shrinking procedure, implemented using the technique of [9], takes O(md/(d+1)nd/(d+1)/Ld/(d+1)) expected time, as long as m and n do not deviate significantly from one another. Agarwal and Matoušek [6] present an algorithm that solves problems (iii)–(v) for a pair of sets of respective sizes m and n, in d in time O(me/(e+1)ne/(e+1)+m+n), where e=d/2+1. (Recall that these problems are mapped, via the lifting transform given in Section 2, to dynamic halfspace range emptiness queries in d+2 dimensions, which explains the value of e; see [6].)

Substituting these components into the algorithm, with parameters L and Δ, the running time becomes

O(nΔ(d1)/(d+1)Ld/(d+1)+L1/2nΔ(e1)/(e+1)+n2e/(e+1)+n(3e+1)/(e+1)Δ),

where, as above, the third term is dominated by the fourth. We balance the terms in this bound, first by choosing L as a function of Δ, and then by choosing Δ as a function of n. As the exponents become rather messy, we simplify them somewhat by introducing another parameter g=2(de)(e+1)(3d+1). Straightforward, albeit rather tedious, calculations then yield

L=Δ2gandΔ=n2e/(e+1)2e/(e+1)+g,

and the bound then becomes

O(ng2e/(e+1)+g+2ee+1).

As a sanity check, we verify that in 3, where d=3, e=2, and g=1/15, the bound is indeed O(n29/21). In conclusion, we have

Theorem 5.

Let P be a set of n points in d, for d3, let s and t be two designated points, and let 0<λ<n be a given parameter. Then we can find the smallest value r for which Gr(P) contains a path between s and t that consists of at most λ edges, in O(ng2e/(e+1)+g+2ee+1) randomized expected time, where e=d/2+1 and g=2(de)(e+1)(3d+1).

As an additional example, in the case d=4, with e=3 and g=1/26, the algorithm runs in O(n61/40)=O(n1.525) randomized expected time.

4 Reverse shortest paths in ball-intersection graphs in 𝟑

As before, we use the shrink-and-bifurcate paradigm (see [7, 9, 11]). As we recall, the shrinking stage receives a number Ln as input, and constructs an interval I that contains r and at most L other critical values of r. For arbitrary balls, the problem is represented in 4, where a ball with center c and radius ρ is mapped to the point (c,ρ). A value r is critical if two expanded balls Br(c1,ρ1)=B(c1,ρ1+r) and Br(c2,ρ2)=B(c2,ρ2+r) become externally tangent to each other, namely when

|c1c2|=(ρ1+r)+(ρ2+r),orr=12(|c1c2|ρ1ρ2).

The improved shrinking procedure of [9] performs “distance selection” on various random samples from , where the selection has to find the k-th smallest critical value between pairs of balls in the cartesian product of the samples, for some given k. This selection process is obtained in turn from a procedure that counts the number of critical values that are smaller than or equal to some given r0. To do so, we rephrase the problem as an offline range searching problem, in which each ball B=B(c1,ρ1) is mapped to the point (c1,ρ1) in 4, and also to the range444Here the lifting transform to a halfspace in 5 does not seem to lead to a more efficient procedure. Technically, this is because the lifted problem to 5 now has to deal with halfspace range counting, whose performance is actually worse than the corresponding problem (involving semi-algebraic ranges instead of halfspaces) in 4.

σr(B)={(c,ρ)4|cc1|ρ2r+ρ1}.

Each point has four degrees of freedom, and, for r fixed, so does each range. Using known results on semi-algebraic range searching (see [2] and [3, Appendix]), the cost of this procedure is O(n8/5). In fact, running this “distance selection” procedure in a bipartite setting (which is the setting that arises when applying the technique of [9]), where we want to select critical values determined between pairs of balls in a pair of samples, consisting of m and n balls, respectively, the cost is O(m4/5n4/5+m+n). Plugging this into the machinery of [9], the cost of the shrinking stage is O(n8/5/L4/5).

The cost of the bifurcation stage, as in all the previous applications, is O(L1/2D(n)), where D(n) is a bound on the running time of the decision procedure, which is our BFS algorithm on Gr()=G0(r), so D(n)=O(n4/3). Hence the overall cost of the procedure is

O(n8/5L4/5+L1/2n4/3).

We optimize this bound by choosing L13/10=n4/15, or L=n8/39, and then the running time becomes O(n56/39). That is, we have:

Theorem 6.

The reverse shortest path problem on the ball-intersection graph of a set of n balls in 3 can be solved in O(n56/39)O(n1.436) time.

4.1 Higher dimensions

Here too, the machinery developed in Section 2 and the preceding part of Section 4 can be extended to any higher dimension d.

Breadth-First Search.

We proceed as in the previous algorithm. The lifting transform used in the BFS implementation maps each ball B0=B(c0,ρ0) to the point B=(c0,ρ0,|c0|2ρ02) in d+2, and also to the halfspace

hB0:xd+22c0c+2ρ0ρ+ρ02|c0|2

in d+2. The dynamic halfspace range reporting structure of [6] yields, for any parameter nsne, where e=(d+2)/2=d/2+1, a dynamic data structure of size O(s), with initial construction cost O(s), so that each query takes O(n/s1/e) time, and each deletion takes O(s/n) amortized time. Choosing s=n2e/(e+1), the size (and construction cost) becomes O(s)=O(n2e/(e+1)), and each operation takes O(n(e1)/(e+1)) time, for a total of O(n2e/(e+1)) time. Since this dominates the running time of the BFS, we obtain:

Theorem 7.

BFS on the ball intersection graph of a set of n balls in d can be performed in O(n2e/(e+1)) randomized expected time, where e=d/2+1.

Reverse Shortest Paths.

Here the selection procedure, rephrased as an offline semi-algebraic range searching problem, involves points and ranges, each point and range with d+1 degrees of freedom. Hence the running time of this procedure, in the bipartite setting of sets of m and n (ball-representing) points, respectively, is (again, see [3, Appendix])

O(m(d+1)/(d+2)n(d+1)/(d+2)+m+n).

This implies, as in [9], that the cost of the interval shrinking procedure is

O(n2(d+1)/(d+2)/L(d+1)/(d+2)).

As in all its implementations, the bifurcation procedure runs in time

O(L1/2D(n))=O(L1/2n2e/(e+1)).

The overall running time is thus

O(n2(d+1)/(d+2)/L(d+1)/(d+2)+L1/2n2e/(e+1)).

Optimizing the value of L, we choose

L12+d+1d+2=n2(d+1)d+22ee+1,orL=n4(de+1)(3d+4)(e+1),

and the overall running time of the algorithm is

O(n2(d+1)(3e+1)(3d+4)(e+1)).

(We verify that, for d=3 and e=2, this bound is indeed O(n56/39).) That is, we have:

Theorem 8.

The reverse shortest path problem on the ball intersection graph of a set of n balls in d can be solved in O(n2(d+1)(3e+1)(3d+4)(e+1)) time, where e=d/2+1.

(We remind the reader that the second improvement in [9] is not applicable here, when the balls do not have the same radius.)

For example, in d=4 dimensions, with e=3, the algorithm runs in O(n25/16)=O(n1.5625) time.

5 Further extensions

Other measures of expansion.

The RSP technique can easily be applied to other well-behaved measures of expansion of the balls. For example, consider, in three dimensions, the case where, for r1, each ball B(c,ρ) is expanded to the ball B(c,rρ) (all radii are multiplied by r). Now r is a critical value iff a pair B(c1,rρ1), B(c2,rρ2) of balls become externally tangent to each other, that is,

|c1c2|=r(ρ1+ρ2),orr=|c1c2|ρ1+ρ2

Hence, in the selection procedure, we map each ball B0=B(c0,ρ0) to the range

σB0={(c,ρ)4|cc0|=r(ρ+ρ0)}.

Since these ranges are also semi-algebraic regions of constant complexity with four degrees of freedom, the same previously used offline semi-algebraic range searching machinery (on a different kind of ranges) can be applied here too, with the same asymptotic running time bound. The remainder of the procedure is carried out as before. Any other expansion rule can also be handled in this manner, as long as the corresponding ranges σB0 are semi-algebraic of constant complexity (with four degrees of freedom).

The same observations hold in any larger dimension d.

In conclusion, Theorems 4, 5, 6 and 8 continue to hold for any well-behaved measure of expansion.

References

  • [1] A. Karim Abu-Affash, Paz Carmi, Matthew J. Katz, and Michael Segal. The Euclidean bottleneck Steiner path problem and other applications of (α, β)-pair decomposition. Discret. Comput. Geom., 51(1):1–23, 2014. doi:10.1007/S00454-013-9550-9.
  • [2] Pankaj K. Agarwal. Simplex range searching and its variants: A review. In Journey through Discrete Mathematics: A Tribute to Jiří Matoušek, pages 1–30. Springer Verlag, Berlin-Heidelberg, 2017.
  • [3] Pankaj K. Agarwal, Boris Aronov, Esther Ezra, Matthew J. Katz, and Micha Sharir. Intersection queries for flat semi-algebraic objects in three dimensions and related problems. ACM Trans. Algorithms, 21(3):25:1–25:59, 2025. doi:10.1145/3721290.
  • [4] Pankaj K. Agarwal, Haim Kaplan, Matthew J. Katz, and Micha Sharir. Segment proximity graphs and nearest neighbor queries amid disjoint segments. In 32nd Annual European Symposium on Algorithms, ESA 2024, volume 308 of LIPIcs, pages 7:1–7:20. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2024. doi:10.4230/LIPICS.ESA.2024.7.
  • [5] Pankaj K. Agarwal, Matthew J. Katz, and Micha Sharir. On reverse shortest paths in geometric proximity graphs. Comput. Geom., 117:102053, 2024. doi:10.1016/J.COMGEO.2023.102053.
  • [6] Pankaj K. Agarwal and Jirí Matousek. Dynamic half-space range reporting and its applications. Algorithmica, 13(4):325–345, 1995. doi:10.1007/BF01293483.
  • [7] Rinat Ben Avraham, Omrit Filtser, Haim Kaplan, Matthew J. Katz, and Micha Sharir. The discrete and semicontinuous Fréchet distance with shortcuts via approximate distance counting and selection. ACM Trans. Algorithms, 11(4):29:1–29:29, 2015. doi:10.1145/2700222.
  • [8] Sergio Cabello and Miha Jejcic. Shortest paths in intersection graphs of unit disks. Comput. Geom., 48(4):360–367, 2015. doi:10.1016/J.COMGEO.2014.12.003.
  • [9] Timothy M. Chan and Zhengcheng Huang. Faster algorithms for reverse shortest path in unit-disk graphs and related geometric optimization problems: Improving the shrink-and-bifurcate technique. In 41st International Symposium on Computational Geometry, SoCG 2025, volume 332 of LIPIcs, pages 32:1–32:14. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2025. doi:10.4230/LIPICS.SOCG.2025.32.
  • [10] Timothy M. Chan and Dimitrios Skrepetos. All-pairs shortest paths in unit-disk graphs in slightly subquadratic time. In 27th International Symposium on Algorithms and Computation, ISAAC 2016, volume 64 of LIPIcs, pages 24:1–24:13. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2016. doi:10.4230/LIPICS.ISAAC.2016.24.
  • [11] Haim Kaplan, Matthew J. Katz, Rachel Saban, and Micha Sharir. The unweighted and weighted reverse shortest path problem for disk graphs. In 31st Annual European Symposium on Algorithms, ESA 2023, volume 274 of LIPIcs, pages 67:1–67:14. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2023. doi:10.4230/LIPICS.ESA.2023.67.
  • [12] Haim Kaplan, Wolfgang Mulzer, Liam Roditty, Paul Seiferth, and Micha Sharir. Dynamic planar Voronoi diagrams for general distance functions and their algorithmic applications. Discret. Comput. Geom., 64(3):838–904, 2020. doi:10.1007/S00454-020-00243-7.
  • [13] Matthew J. Katz, Rachel Saban, and Micha Sharir. Near-linear algorithms for visibility graphs over a 1.5-dimensional terrain. In 32nd Annual European Symposium on Algorithms, ESA 2024, volume 308 of LIPIcs, pages 77:1–77:17. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2024. doi:10.4230/LIPICS.ESA.2024.77.
  • [14] Katharina Klost. An algorithmic framework for the single source shortest path problem with applications to disk graphs. Comput. Geom., 111:101979, 2023. doi:10.1016/J.COMGEO.2022.101979.
  • [15] Jirí Matousek. Reporting points in halfspaces. Comput. Geom., 2:169–186, 1992. doi:10.1016/0925-7721(92)90006-E.
  • [16] Haitao Wang and Yiming Zhao. Reverse shortest path problem for unit-disk graphs. J. Comput. Geom., 14(1):14–47, 2023. doi:10.20382/JOCG.V14I1A2.