Ray Shooting Amidst Spheres in Three Dimensions and Related Problems.
We consider the problem of ray shooting amidst spheres in 3-space: given n arbitrary (possibly intersecting) spheres in 3-space and any $\epsilon$ > 0, we show how to preprocess the spheres in time $O(n^{3+\epsilon})$ into a data structure of size $O(n^{3+\epsilon})$ so that any ray-shooting query can be answered in time $O(n^\epsilon)$. Our result improves previous techniques (see [P. K. Aggarwal, L. Guibas, M. Pellegrini, and M. Sharir, "Ray shooting amidst spheres," unpublished note] and [P. K. Aggarwal and J. Matousek, Discrete Comput. Geom., 11 (1994), pp. 393-418]), where roughly $O(n^4)$ storage was required to support fast queries. Our result shows that ray shooting amidst spheres has complexity comparable with that of ray shooting amidst planes in 3-space.  Our technique applies to more general (convex) objects in 3-space, and we also discuss those extensions.
Introduction
The ray shooting problem can be defined as follows:
Given a collection S of n objects in IR d , preprocess S into a data structure, so
that one can quickly determine the first object of S intersected by a query ray.
The ray shooting problem has received considerable attention in the past few years because
of its applications in computer graphics and other geometric problems [1, 4, 5, 6, 10, 11, 13,
16, 21]. Most of the work to date has studied the planar case, where S is a collection of line
segments in IR 2 . Chazelle and Guibas proposed an optimal algorithm for the special case
where S is the boundary of a simple polygon [16]. Their algorithm answers a ray shooting
query in O(log n) time using O(n) space; simpler algorithms, with the same asymptotic
performance bounds, were recently developed in [13, 24]. If S is a collection of arbitrary
segments in the plane, the best known algorithm answers a ray shooting query in time
Work on this paper has been supported by NSF Grant CCR-91-22103, and by grants from the U.S.-
Israeli Binational Science Foundation, and the G.I.F., the German-Israeli Foundation for Scientific Research
and Development, and the Fund for Basic Research administered by the Israeli Academy of Sciences. This
paper is part of the first author's M.Sc. thesis, prepared under the supervision of the second author.
y School of Mathematical Sciences, Tel Aviv University
z School of Mathematical Sciences, Tel Aviv University, and Courant Institute of Mathematical Sciences,
New York University
s log O(1) n) using O(s 1+" ) space and preprocessing 1 [1, 6, 10], where s is a parameter
that can vary between n and n 2 . Although no lower bound is known for this case, it is
conjectured that this bound is close to optimal.
In spite of some recent developments, the three-dimensional ray shooting problem seems
much harder and it is still far from being fully solved. The general three-dimensional ray
shooting problem is to preprocess a collection S of n convex objects, so that the first object
hit by a query ray can be computed efficiently. Most studies of this problem consider the case
where S is a collection of triangles. If these triangles are the faces of a convex polyhedron,
then an optimal algorithm, with O(logn) query time and linear space, can be obtained using
the hierarchical decomposition scheme of Dobkin and Kirkpatrick [19]. If the triangles form
a polyhedral terrain (a piecewise-linear surface intersecting every vertical line in exactly one
point), then the technique of Chazelle et al. [15] yields an algorithm that requires O(n 2+" )
space and answers ray shooting queries in O(log n) time. Non-trivial solutions to the general
problem (involving triangles) were obtained only recently; see [4, 6, 11] for some of these
results. The best known algorithm for ray shooting among triangles in three dimensions is
due to Agarwal and Matou-sek [5]; it answers a ray shooting query in time O( n 1+"
s 1=4 ), with
O(s 1+" ) space and preprocessing. The parameter s can range between n and n 4 . If s assumes
its maximum value, queries can be answered in O(log 2 n) time; see [5, 6] for more details.
A variant of this technique was given recently in [7] for the case of ray shooting amidst
a collection of convex polyhedra. We remark that no non-trivial lower bounds are known
for the 3-dimensional problem (for triangles) as well, although such bounds are known for
the related simplex range searching problem [12], which is used as a subprocedure in the
solutions just mentioned. These bounds are close to \Omega\Gamma n
s 1=4 ), and thus might suggest that
the known ray shooting algorithms for triangles in IR 3 are close to optimal.
On the other hand, there are certain special cases of the 3-dimensional ray shooting
problem which can be solved more efficiently. For example, if the objects of S are planes
or halfplanes, ray shooting amidst them can be performed in time O( n 1+"
space and preprocessing; see [4] for details. It is therefore of interest to identify additional
classes of objects for which ray shooting can also be performed more efficiently.
In this paper we consider the case of spheres. This problem has been recently studied
in [5], and also in a yet unpublished work [3]. The algorithm presented in [5] has query
time O( n 1+"
s 1=4 ), with O(s 1+" ) space and preprocessing; that is, the same performance as in
the case of triangles. The algorithm of [3] achieves similar performance for
We present an improved solution for the case where large storage is allowed. Specifically,
we achieve query time O(n " ) with only O(n 3+" ) space and preprocessing. Thus we show
that, in some sense, the case of spheres is essentially no harder than the case of planes,
which might sound somewhat surprising, given that spheres are more complex objects (and
require 4 real parameters to specify, as opposed to only 3 for planes). Our solution is fairly
general-the spheres may have arbitrary radii and may also intersect one another.
1 Throughout this paper, bounds of this kind mean that, given any arbitrarily small positive constant ",
the algorithm can be fine-tuned so that its performance satisfies the bound; the multiplicative constants in
such bounds usually depend on " and tend to 1 as " # 0.
Our technique adapts some ideas from [3]. Roughly speaking, in both approaches we
reduce the ray-shooting problem to the following problem (although certain parts of our
reduction are simpler than those of [3]): we are given a collection S of n arbitrary spheres
in 3-space, S we want to preprocess them so that we can answer the following
queries efficiently: we are given a line ', such that for each sphere S i either ' intersects
S i or else S i contains no point that lies vertically above some point on ', and we want
to determine whether ' intersects any of these spheres. The reduction of the ray shooting
problem to this problem is fairly routine, using mostly standard machinery.
We solve this main subproblem in a more careful manner than in [3]. This subproblem
can be mapped onto the 4-dimensional parametric space that represents lines in 3-space,
and reduces to the problem of point location in the region of this 4-space lying above
the upper envelope of certain low-degree algebraic surface patches, each representing an
'upper tangency' between a line and one of the spheres. Almost a year after the original
submission of this paper, it was shown by Aronov et al. [2] that point location above
such an envelope can be performed fast (in polylogarithmic time), using only O(n 3+" )
storage and O(n 3+" ) randomized expected time. This gives an alternative solution with
roughly the same asymptotic performance as our solution. However, (a) the algorithm of
[2] is randomized, whereas our solution is deterministic, and (b) our solution analyzes and
exploits the geometric structure of the problem much more explicitly, and we hope that
its analysis will find additional applications. Although our approach is more geometric in
nature, it does not exploit any special properties of spheres, and is thus much more general.
As a matter of fact, our technique can be viewed as an explicit solution of the problem of
point location above an envelope in 4-space, for the special type of surfaces arising in the
context of ray shooting. In Section 3 we indeed present a generalization of our technique to
ray shooting among more general convex objects in 3-space, and exemplify it for the case
of axis parallel ellipsoids.
2 The Algorithm

Overview

We use the general approach to ray shooting due to Agarwal and Matou-sek [4]. This
approach, which is based on the parametric searching technique of Megiddo [27], reduces
the ray shooting problem to the segment emptiness problem. Namely, we need to preprocess
the spheres, so that, given any segment e in 3-space, we can quickly determine whether e
intersects any of the spheres. As shown in [4], the time for the actual ray shooting query is
only within a logarithmic factor off the query time for the segment emptiness problem.
We first describe the preprocessing scheme. We construct a multi-level data structure,
where each level of the structure filters out spheres that satisfy a different geometric relationship
with respect to the query segment, so that a sphere satisfying the conjunction of
all these relationships must intersect the query segment; see [26] for a more detailed discussion
of multi-level data structures of this kind. Each level of the data structure consists
e

Figure

1: Illustration of condition (ii) for segment intersection detection
of a tree-like structure T , where each node of T corresponds to some 'canonical' subset of
spheres. The root of the top-level structure represents the entire collection S of spheres.
First, here is a brief overview of our approach. It is fairly easy to show that a segment
e intersects a sphere S if and only if one of the following two conditions hold:
(i) One endpoint of e lies inside S and the other lies outside S.
(ii) Both endpoints of e lie outside S, the center of S lies inside the slab \Sigma e , bounded by
the two planes passing through the endpoints of e and perpendicular to e, and the
line containing e intersects S. See Figure 1 for an illustration.
Consequently, we construct two data structures, one for testing whether condition (i) holds
for any sphere, and the other for testing condition (ii).
First data structure
We start by describing the first, and simpler, data structure. Its first level T 1 is used to find
all spheres of S containing a query point. The root of T 1 is associated with the entire set S.
We fix some sufficiently large constant parameter r, and construct a (1=r)-net E of O(r log r)
spheres from S, for certain ranges (subsets of S) that will be defined in a moment. (E is a
(1=r)-net if each range that contains more than jSj=r elements must contain an element of
see [23] for more details.) E is constructed by the deterministic technique of Matou-sek
[25] (or by simply choosing a random sample of O(r log r) spheres). We next construct the
arrangement A(E) of these spheres, and apply the vertical decomposition technique of [18]
(see also [14]) to that arrangement; it decomposes 3-space into O(r 3+" ) cells of 'constant
description complexity', so that the interior of each cell is not crossed by any of the spheres
of see [18] for more details. The ranges with respect to which E has to be a (1=r)-net are
sets of spheres of S, each consisting of those spheres intersecting a region having the shape
of a cell of the vertical decomposition in an arrangement of spheres, as above. As is easy
to see, the resulting range space has so-called finite VC-dimension, which thus implies the
existence of a (1=r)-net E with only O(r log r) spheres (again we refer the reader to [23] for
more details). For each cell C of the decomposed arrangement we create a child of the root
of T 1 , and associate with it the set SC of spheres intersecting C. Since E is a (1=r)-net, none
of the sets SC contains more than n=r spheres. We also store with each cell C the set S 0
C of
all spheres that completely contain C in their interior. We then continue the preprocessing
recursively at each child of the root, with its associated set SC .
The second-level structure is built for each of the nodes of the first level. The set of
spheres associated with the root of the second-level structure at a node corresponding to
some cell C in the first-level decomposition, is the set S 0
C defined above. The purpose of the
second-level structure is to test whether any of the given spheres does not contain a query
point. This is achieved using almost the same preprocessing as for the first level, except
that here we define the set S 0
C to be the set of all spheres whose enclosed ball is disjoint
from the cell C.
Given a query segment first search with p through the first-level structure,
obtaining the set of spheres containing p as a disjoint union of sets S 0
C . For each of these
sets, we search with q through the corresponding second-level structure; if we find there
a node whose set S 0
C is nonempty, we stop and conclude that e intersects a sphere of S.
Otherwise, we repeat this procedure with searching with q in the first-level structure and
with p in the second-level structures. Again, we either detect an intersection between e and
some sphere, or else conclude that no intersection of type (i) occurs, and then we move to
search in the second data structure.
Second data structure
The first two levels of the second data structure are very similar to those of the first data
structure. Their purpose is to filter out all spheres that do not contain any of the endpoints
of the query segment e. These levels are constructed in much the same way as above, except
that the sets S 0
C , in both levels, are the sets of spheres whose enclosed balls are disjoint
from the corresponding cell C. We search with e in these two levels as described above, and
the search results in a collection of canonical sets S 0
C of spheres, whose (disjoint) union is
the subset of all spheres whose enclosed balls do not contain any endpoint of e.
The next two levels of the structure aim at finding all spheres whose centers lie in the
slab \Sigma e as defined in condition (ii) above. In both levels we take the set of centers of the
relevant spheres, and preprocess them for halfspace range searching in 3 dimensions (using,
e.g., the method of [17]). Since \Sigma e is the intersection of two halfspaces, H 1 , H 2 , we search
with H 1 in the third-level structures and with H 2 in the fourth-level structures, and wind
up with a collection of canonical sets of spheres, whose (disjoint) union is the set of all
spheres that satisfy the first two parts of condition (ii). Hence, for each of the resulting sets

Figure

2: The planes V and - and their relationships to a sphere
of spheres, we only need to determine whether the line ' containing the query segment e
intersects any sphere in that set.
The fifth level of our structure is also a halfspace range searching structure on the centers
of the relevant spheres. We search in that level with the two halfspaces bounded by the
plane - passing through ' and orthogonal to the vertical plane V passing through '. Let
denote the halfspace lying above -, and let - \Gamma denote the halfspace lying below -. Let
S be a sphere whose center lies in - + and which intersects V in a disc D. Then the center
of D lies above ', so either ' intersects S or else it passes below S, in the sense that ' and S
are disjoint, and there is a point on ' that lies vertically below a point in S. See Figure 2 for
an illustration. A symmetric property holds if the center of S lies in
the description of our data structure so that it only handles the latter situation, since the
handling of the former case is fully symmetric.
To recap, we have reduced our problem to the following subproblem. We are given a
collection of arbitrary spheres in 3-space, which is a canonical set of spheres in the output
of the fifth level of our structure, and we want to preprocess it so that we can answer
the following queries efficiently: we are given a line ', such that for each sphere S in our
collection either ' intersects S or else there is no point of S that lies vertically above '. We
want to determine whether ' intersects any of the given spheres. For brevity of notation,
we call the given collection S and assume its size to be n.
Detecting intersection between a line and the upper hull of a vertical planar
cross section of S
We begin by describing a few constructs that will be used in our technique. For each point
denote the vertical plane passing through the line We

Figure

3: Two bridges connecting a disc D i to two larger discs
denote by (ae; z) the coordinates within V , where ae is the horizontal coordinate (along the
line z is the vertical coordinate. Let be a subset of
t spheres of S. Let D denote the (possibly empty) 2-dimensional discs formed by
intersecting the corresponding spheres S . We consider the upper convex
hull of these discs. This hull consists of bridges, namely common tangents between pairs of
discs which form edges of the hull, and of circular arcs which are portions of the boundaries
of the discs and which lie between two adjacent bridges; see Figure 3. Suppose the disc D i
is part of the upper convex hull of these t discs (within V ), and consider the set of bridges
that connect D i with other discs, whose radius is greater than or equal to that of D i . It is
easily checked that there can be at most two such bridges, at most one connecting D i to a
disc D j following it in the positive ae-direction, and at most one connecting D i to a disc D k
preceding it; see again Figure 3. (This claim is trivial if the discs are pairwise disjoint, but
also holds if they are allowed to intersect.) For each bridge b we denote by b + the halfplane
lying above the line containing it. For each circular arc fl of the hull we denote by fl + the
region above it (that is, we erect a vertical ray upwards from each point on fl and take the
union of all these rays); see Figure 4.
Here is a quick overview of our approach: we fix some sufficiently large constant parameter
r and construct a set E ' S of size log r), which will serve as a (1=r)-net with
properties that will be explained in a moment. Consider a query line ' and the vertical
passing through it. We will build a data structure on E , from which we will be able
to determine quickly the set of bridges of the upper convex hull of the t intersection discs of
the spheres of E with V . We will then check if ' intersects this hull. If it does we are done
intersects a sphere of S. Otherwise, we find two bridges b; b 0 whose slopes are nearest
to the slope of ', so that ' lies completely in the union of the regions b
the circular arc between b and b 0 ; see Figure 4. We then continue the query recursively in
the set of spheres intersecting each of these 3 regions. For this approach to be efficient we
need each of the regions b to be intersected by a small number of spheres (at most n=r
spheres). We will actually define here two range spaces, where each of the ranges of each
of these spaces is defined in terms of a constant number of spheres, and both range spaces
have finite VC-dimension. The first range space deals with ranges related to the regions
, and the second range space deals with ranges related to the regions fl + . We will build
a (1=r)-net E 1 of size O(r log r) for the first range space and another (1=r)-net E 2 of size
O(r log r) for the second range space, so the union is a (1=r)-net for both range
spaces, of size

Figure

4: A line ' passing above the upper convex hull, and the three regions whose union
contains '
The planar maps M \Sigma (S i )
We now describe our technique in detail. We build a data structure so that given any
vertical plane
represented by the point (- the bridges and circular
arcs of the upper convex hull of the discs D i within V can be found quickly. For each sphere
construct the following two planar maps M in the - 1 - 2 -coordinate
system, where - 1 is the slope and - 2 is the intercept of the xy-projection of a line in 3-space.
If D i contibutes to the upper convex hull and the unique larger disc following it along the
hull in the positive ae-direction is D j , we label in M by the index j;
similarly, if D i contibutes to the upper convex hull and the unique larger disc following it
along the hull in the negative ae-direction is D k , we label in M \Gamma (S i ) the point (- by the
index k; otherwise we label (- in the respective map. (this alternative labeling
applies also in the cases where D i is empty or does not appear at all along the upper hull).
The faces of are maximal connected regions, all of whose points have
the same label; edges and vertices of these maps are defined accordingly. The edge between
two adjacent cells of the map M are labeled by j and k respectively,
consists of points (- the disc following (preceding) D i along the hull is about
to change from D j to D k , or the radius of D j becomes equal to the radius of D i , or D i is
getting out of the hull boundary by another disc. Thus each edge e of any of these maps
consists of points (- at which the upper hull of the discs in the plane
has a bridge
tangent to three discs, or has two adjacent discs of equal radius, or has a disc degenerating
to a single point. See below for a more detailed analysis of these maps; we will show that
the overall complexity of all these maps is only roughly cubic in r. We next construct a
vertical decomposition of each of the planar maps M erecting a vertical
segment from each vertex and from the points with - 2 -vertical tangency on each edge, and
extend it upwards and downwards until it hits another edge of the map. Note that the
set of labels of a point (- these maps provides complete information about the
structure of the upper convex hull of the intersection discs D i in the vertical plane
Let C be a cell in, say, the vertical decomposition of the map M suppose that
ijk
ik
ijl

Figure

5: A cell C labeled by j ? 0 in the map M
its label is j ? 0 (cells with label 0 can be ignored in this process). For each point (-
C, we consider the plane defined above. We know that S i contributes to the
upper convex hull of the discs of intersection of the spheres of E with this plane, and that
the unique next larger disc along the hull in the positive ae-direction is D j . Let b denote
the bridge connecting the two discs in the upper convex hull; we denote by b
halfplane of V above the line containing b (that is, the halfplane supporting the discs D i ,
D j from above). Consider the 3-dimensional region
be the subset of all the original spheres intersecting RC . We define similar regions RC for
all the cells C in all the other maps M \Sigma (S k ), for t.
Note that each of the regions RC is the union of halfplanes b is a bridge
determined by a fixed pair of spheres of E , and where (- ranges over a cell C of one
of the maps M \Sigma (S i ), which is a region of constant complexity. Actually, suppose C has
label defined in terms of at most 6 spheres: at most three spheres,
and another sphere S k , define the top edge of C, at most three spheres, S i , S j and another
l , define its bottom edge, at most one more sphere is needed to define the left vertical
edge of C, and at most one more sphere is needed to define its right edge. See Figure 5 for
an illustration. It follows that C is defined by at most 6 spheres of E . Moreover, the two
spheres forming the bridge b belong to this collection of 6 spheres, so RC is also defined
in terms of at most 6 spheres, and is thus a region of constant complexity. The ranges
with respect to which E 1 (and hence E) has to be a (1=r)-net are subsets of the form SC ,
defined for regions RC which are defined in terms of at most 6 spheres, as just outlined.
Again it is easy to verify that the resulting range space has finite VC-dimension. Since,
by construction, none of the spheres of E intersects any of the regions RC , it follows by
definition that the maximum cardinality of the sets SC , over all cells C of all the 2t maps
M \Sigma (S i ), is at most n=r.
We construct the sets SC using any brute-force method, which takes O(n) time, since
Figure

The possible xy-projections of the portions S
(the projections have smooth
boundaries, but some are drawn with nonsmooth boundary to illustrate the fact that they
are composed of two different arcs)
we assume r to be constant. We then create, for each cell C in each map, a child of the root
node of T , associate the set SC of spheres with that child, and continue the construction
of the data structure recursively at each child with the corresponding set of spheres. There
is, however, a second set of children of the root, obtained through a second decomposition
scheme, which we now proceed to describe.
The spatial maps N(S i )
We now define our second range space, with respect to which E 2 (and hence E) has to be
a (1=r)-net. Let N denote the 3-dimensional parametric space with coordinates (-
where, as above, (- give the dual representation of a line -
in the xy-plane, and where ae measures the distance along that each point
actually represents (in a many-to-one manner) a point in the xy-plane,
namely the point oe(- along the line - The reason for this
3-dimensional representation of the xy-plane will become clearer later on.
For each of the t spheres S we form in N a spatial subdivision N(S i ), consisting
of roughly O(t 2 ) cells of constant complexity, as follows. For any other sphere S
j denote the portion of S j that lies above S i , and let S
j denote the xy-projection of S
.
We obtain a collection of at most t planar regions, each of which is the intersection of Q
the xy-projection of S i , with either a disc, or an ellipse, or a convex region with a smooth
boundary which is the union of a circular arc and an elliptic arc; see Figure 6 and Figure 7.
Let K i denote the complement (within Q i ) of the union S
. We map K i into the
following set in our parametric space
~
In other words, ~
K i is the pre-image of K i under the

Figure

7: The arrangement of the xy-projections of the portions S
mapping oe from N onto the xy-plane, as defined above. Note that a point (- ae) is in
~
only if the vertical line (in actual 3-space), passing through the point on
with coordinate ae, does not meet any sphere of E above its highest intersection with S i .
We next decompose ~
cells by applying a standard vertical
decomposition (with ae being the 'vertical' coordinate), as in [14, 18]. This yields the spatial
As observed in the papers just cited (and is easily seen), the complexity (number
of resulting cells) of N(S i ) is proportional to the number of vertically visible pairs
of edges of ~
(where we include among these edges also the loci of points with vertical
tangency along the boundary of ~
Such a pair of edges vertically visible if there
exists a point (- such that the vertical line in N passing through that point intersects
both e and e 0 , and the segment between these points of intersection lies fully in ~
edge e of ~
ae) such that the line -
of K i or is tangent to the boundary of K i , and ae is the coordinate of that vertex or point of
tangency along Thus each vertically-visible pair (e; e 0 ) of edges of ~
to a line - through two points, each being either a vertex of K i or a
point of tangency between and the boundary of K i , so that the interval between
these two points is fully contained in K i .
In other words, the number of vertically-visible pairs (e; e 0 ) as above is proportional to
f , where the sum extends over all faces f of K i , and where c f is the complexity of
f . All these faces are faces of a planar arrangement of O(t) circular and elliptic arcs, and
each pair of these arcs intersect in at most 4 points. Hence, as shown in [20], we have
is the (nearly-linear) maximum length of a
(t; sequence [8]. We have thus shown that the number of cells (each
having constant complexity) in all the subdivisions N(S i ) is O(t- 2
6 (t)), and is thus nearly
cubic in r (recall that
We next associate a set SC of spheres with each cell C of N(S i ). For each point
denote the vertical line (in actual 3-space) passing through
the point of coordinate ae on the line - denote the portion of
above its highest intersection point with S i . Let RC denote the union of all the
C. The set SC consists of all spheres of S that intersect
RC . We note that RC is a region of constant complexity, defined in terms of only a constant
number of spheres of E . The ranges with respect to which E 2 (and hence E) has to be a
(1=r)-net are subsets SC , each consisting of those spheres that intersect a region of the
form RC , defined for cells C in N that have the same structure as the cells of the vertical
decompositions constructed above. Again, it is easily seen that the resulting range space
has finite VC-dimension. (Intuitively, this is due to the fact that each of the regions RC
has constant description complexity; see, e.g., [29] for more details.)
We next observe that none of the spheres of E intersect RC . Indeed, if one of these
spheres, S j , did intersect RC , then there would exist a point (-
K i such that the
vertical line '(- at a point that lies above S i . But then, by definition, the
line the xy-plane at a point belonging to S
, so that, by definition of K i ,
ae does not belong to K i does not belong to ~
We thus conclude that the maximum cardinality of the sets SC , over all cells C of N(S i ),
and over all spheres S i of E , is at most n=r. We construct the sets SC , by any brute-force
method, in linear time (r is assumed to be a constant), create a new set of additional
children of the root of the tree T , one child for each cell C in each of the subdivisions N(S i ),
associate with that node the corresponding set SC of spheres, and continue the construction
of the data structure recursively at each child with the corresponding set of spheres. This
concludes the description of the sixth (and last) level of our data structure.
Answering a query
We next describe how a query is processed. Let e = pq be a query segment in 3-space.
We have already described how to search with e in the first data structure, and how to
search with e in the first five levels of the second data structure. Let S be a canonical set
of spheres in the output of the search in the first five levels of the second structure, and
let ' be the line containing e. We assume, with no loss of generality, that for each sphere
or else S has no point that lies vertically above '. Our goal is to
determine whether ' intersects any sphere in S.
Suppose ' is represented by the two equations We start
at the root of the sixth-level tree of S, and locate the point (- in each of the 2t maps
of the root. Let C denote the cells of these
maps which contain (- whose label is not 0. As noted, the labels of these cells
provide a complete description of the upper convex hull of the t intersection discs D i of our
net spheres with the plane V - 1 ;- 2 . We compute the slopes - of the k bridges of that
hull, and compare each of them with the slope of the query line within V , which is, as is
easily verified, -
1 . If one of the bridges, say b, has slope - simply
check whether our query line passes below b, in which case ' intersects a sphere of S, so we
can return a positive answer to the query. Otherwise, ' passes above (the line containing)

Figure

8: An original sphere intersecting V in a disc meeting the small region between D i
and the extension of its two bridges
b, and we recurse in the sub-tree corresponding to the cell C that defines b (only spheres in
SC can meet '). In general, though, we will obtain two adjacent bridges, b, b 0 , both tangent
to the same sphere of E , say S i , such that - b - b 0 . In this case, if an original sphere S
intersects ', then either S meets one of the two halfplanes b
meets the small region R enclosed between S i and the two lines containing b, b 0 ; see Figure 8
for an illustration. Let fl denote the (upper) arc of D i between its two points of tangency
with b and b 0 . It is clear that no sphere of E intersects the region consisting of all points
of
that lie above fl, which implies that the xy-projection fl   of fl is fully contained in
(a single connected component of) the region K i defined above. It thus follows
that the set ~ fl, consisting of all points (- ae is the coordinate along
a point of fl   , is fully contained in ~
in fact is fully contained within a single cell of
the subdivision N(S i ) of ~
is a vertical segment contained in ~
and such a segment
must be fully contained in a single cell of its vertical decomposition).
We thus continue the search recursively at the two nodes of T associated with the cells
corresponding to the bridges b, b 0 , and with the node corresponding to the cell C of N(S i )
containing the arc ~
fl. This concludes the description of the query processing.
Complexity analysis
We next analyze the complexity of our algorithm in terms of space, preprocessing time, and
query time. The first five levels of our second structure (as well as the two levels of the
first structure) involve halfspace range searching and structures for point location among
spheres in IR 3 . At each of these levels, the size of the structure for a set of m spheres,
excluding substructures at deeper levels, is O(m 3+" ), for any " ? 0. Indeed, for halfspace
range searching this follows from [17]. For point location among spheres, this follows by
noting that the maximum storage S(m) and preprocessing cost P (m) for this structure on
spheres satisfy the following recurrences:
and
for appropriate constants c 1 , c 2 and c 0
2 (some of which depend on "). The solution of these
recurrences is easily seen to be O(m 3+" 0
where the constant of proportionality depends on ".
As observed in [26] (and is easy to verify), the overall size and preprocessing cost of a
multi-level data structure of the type considered here can be deduced from the maximum
size and preprocessing at any fixed level. In particular, if we show that the sixth level of
our structure also requires O(m 3+" ) storage and preprocessing for a set of m spheres, it will
follow that the overall storage and preprocessing cost of the full multi-level structure will
also be O(n 3+" ) (with a larger constant of proportionality); see [26] for more details.
We thus turn to consider the cost of the sixth level of our structure, which involves
the upper convex hull structure of the intersection discs. We first claim that the total
complexity of all the 2t planar maps M roughly O(r 3 ). To see this,
ignore for now the added complexity caused by the vertical segments forming the vertical
decompositions of the maps (this increases the complexity by only a small constant factor).
Fix a sphere S i 2 E and let C be a cell of, say M label j. Each edge
e on the boundary of C consists of points (- at which either the sphere S j following
in the upper hull is about to change, or S i itself is about to disappear from the upper
hull, or the sizes of the two intersection discs D i , D j become equal. In the first two cases,
either D i shrinks to a point and then disappears from V , or the (line containing the) bridge
connecting D i and D j in the plane becomes tangent to a third disc D k (it is also
possible that D i and D j become tangent to each other); in the first case D k is between D i
and D j or after D j , and in the second case D k precedes D i . Each vertex v of C is formed by
the intersection of two such boundary edges e, e 0 . Suppose first that both e and e 0 represent
events of triple tangency. Then v must be of one of the following three types:
ffl Each of the edges e, e 0 represents a change in the sphere following S i . Denote the two
corresponding new spheres that are about to replace S j as the next sphere by S k and
S l . Since v lies on e, the discs D i , D j and D k have a common tangent that appears
as a bridge in the upper convex hull. Similarly, since v also lies on e 0 , the same is true
for D i , D j and D l . Hence at v the four discs D i , D j , D k , and D l have a common
tangent bridge. See Figure 9(a) for an illustration.
ffl One of the edges, say e, represents a change in the sphere following say from
to some S k , and the other edge e 0 represents a situation where S i is displaced from
the upper convex hull by another sphere S l preceding S i along the hull. As is easily
verified, in this case too the four corresponding discs D i , D j , D k , and D l must have
a common tangent bridge at v. See Figure 9(b) for an illustration.
ffl Each of the edges e, e 0 represents a situation where S i is displaced from the upper
convex hull, by two respective other spheres S k , S l preceding S i along the hull. In
this case too the four corresponding discs D i , D j , D k , and D l must have a common
tangent bridge at v. See Figure 9(c) for an illustration.
l
l
(a) (b) (c)

Figure

9: The three possible configurations that can occur at a vertex of the map M
We have thus shown that each vertex of M these types is a point (-
the plane
contains a line tangent to four intersection discs and passing above all other
spheres. Similar arguments imply that this also holds for all the other maps M \Sigma (S k ). Note
that, assuming general position of the spheres, each such vertex can be present in at most
eight of the 2t maps, twice for each of the at most four spheres that define it (it has no
effect on the labeling of cells in other maps). Similarly, each edge in one of these maps is
defined in terms of at most three spheres, and can thus appear only in the maps of those
three spheres.
Note that if two edges e, e 0 , both representing situations where the size of D i becomes
equal to that of the following disc on the upper hull, meet at a point v, then v also lies on
an edge of a triple tangency, as is easily seen.
The next case to consider is thus of vertices formed by the intersection of two edges e,
consists of points (- at which the upper hull of the intersection discs in the
contains a bridge tangent to D i and to two other discs, D j , D k , and e 0 consists
of points (- at which the size of D i becomes equal to the size of one of these discs, say
. It is easily seen that the number of such vertices is only O(t 3 ), over the entire collection
of maps M \Sigma (S i ), because each such vertex is defined by 3 spheres of E , and each triple of
spheres give rise to only a constant number of such vertices.
Finally, consider cases where a disc D i shrinks to a point, or become tangent to an
adjacent disc along the upper hull. Again, it is easily verified that a vertex v of M \Sigma (S i ) at
which this occurs is defined in terms of at most 3 spheres (S i plus two other spheres defining
an edge of triple tangency on which v lies, or S i , S j , and another sphere becoming tangent
to the common tangent line to D i and D j ). Thus the total number of such vertices, over
all maps, is O(t 3 ).
We thus need to bound the number of points (- at which the upper convex hull of
the intersection discs between the spheres and the vertical plane V - 1 ;- 2 has a bridge tangent
to four discs. To do so, we apply the following transformation to the problem. For each
sphere
the intersection disc
define F S (- to be the value of - 4 for which the line
tangent to D from above. Note that F S is (partially) well defined, and that its graph is an
algebraic surface patch of low constant degree.
It follows by definition that the upper envelope of the functions F S represents the locus
of all lines that are tangent to the upper convex hull of any vertical planar cross section of
the spheres of E . Moreover, each vertex of the upper envelope represents a line tangent to
such an upper convex hull and touching four of the corresponding intersection discs, and
each such 'critical' line is indeed a vertex of the envelope.
By the recent results of [22, 28], the complexity of the upper envelope of n algebraic
surface patches of constant maximum degree in d-space is O(n d\Gamma1+" ), for any " ? 0. Hence
it follows that the complexity of the upper envelope of the log r) functions F S is
O(r 3+" ), for any " ? 0. This implies that the number of vertices of the types considered
above in our 2t maps M \Sigma (S i ) is O(r 3+" ), which is thus also a bound on the overall complexity
of the vertical decompositions of these maps. As already argued, this bound also dominates
the overall complexity of the spatial subdivisions N(S i ), and implies that the number of
children of any node in any sixth-level tree structure is O(r 3+" ).
We can finally complete the analysis of the performance of our algorithm. Assuming r
to be a (sufficiently large) constant, the construction of the (1=r)-net E can be done in O(n)
time, as in [25]. The construction of the 2t maps M \Sigma (S i ) can be done in constant time.
More precisely, we first need to construct the upper envelope of the functions F S i
. As shown
in [2], this can be done in time O(r 3+" ), for any " ? 0. The features of these envelopes are
then distributed among the maps M \Sigma (S i ), the additional edges and vertices of these maps
are constructed in a straightforward manner, and the maps are then vertically decomposed,
within the same asymptotic running time.
The construction of the t spatial subdivisions N(S i ) can also be performed in constant
time. For each fixed sphere S first compute the 'top portions' S
j of the other
spheres of E , project them onto the xy-plane, and compute the complement K i of their
union. We then construct all vertices of N(S i ) by examining every pair of vertices and/or
edges of each face of K i , in the manner described above, from which the complete subdivision
is easy to construct.
For any cell C in each of these maps, we compute the set SC of all spheres intersecting its
associated region RC . This can be done in O(n) time, since each such region is of constant
complexity and the number of regions is also constant. More precisely, all the canonical
sets SC can be computed in O(nr 3+" ) time. As already noted, the maximum size of any set
SC is - n=r.
respectively, the expected space complexity, pre-processing
time, and query time for the sixth level of the data structure on a canonical set
of n spheres. Then we have the following reccurence relations:
and
for appropriate constants c 1 , c 2 and c 0
2 (some of which depend on "). As above, the solution
of these recurrences is easily seen to be O(n 3+" 0
the constant of proportionality depends on ".
As for the query time Q 6 (n), we have the following recurrence:
(the factor 3 appears since we search recursively in at most 3 children of any node that is
visited during the query processing). The solution of this recurrence is easily seen to be
again with a constant of proportionality that depends on ".
If we now combine the sixth level of the structure with the preceding five levels, and
apply the observations of [26] concerning multi-level structures, as mentioned above, and
also take into account the overhead of the parametric searching, we easily conclude the
following main result:
Theorem 2.1 Given a collection S of n arbitrary spheres in 3-space, and any " ? 0, one
can preprocess S, in time O(n 3+" ), into a data structure of size O(n 3+" ), which supports
ray shooting queries among the spheres of S in time O(n " ) per query.
Extensions
In this section we extend our algorithm so that it applies to more general classes of objects.
As a matter of fact, we present a general approach to ray shooting in 3 dimensions amidst
a collection S of n (possibly intersecting) convex objects, each of constant description com-
plexity. The approach is modular, and consists of several stages, where each stage filters
out objects that satisfy a certain geometric constraint with respect to the query ray, so
that the query ray intersects all objects that satisfy the conjunction of all these constraints;
using this property, the actual ray shooting is then easy to perform. Thus, to obtain a
ray shooting algorithm for a specific class of objects, we need to provide appropriate and
efficient filtering mechanisms for this class at each level separately. These mechanisms and
their efficiency depend on the class in question; we will exemplify them for two classes of
objects: triangles and ellipsoids.
As above, we only consider the case where we want to obtain fast queries (O(n "
per query), and are willing to use large storage. We follow the general approach of [4],
reducing the ray shooting problem to the segment emptiness problem. We then proceed
through the following stages:
(1) We begin by reducing the segment emptiness problem to the line emptiness problem,
and we accomplish this in the first two stages. Let e = pq be the query segment and let
' be the line containing e. In the first stage we want to partition our objects into three
categories:
(i) Objects that contain exactly one endpoint of e in their interior. If there is any such
object, we can stop the query immediately with a positive answer. We can compute
this subset of objects (and also the two other subsets, in (ii) and (iii) below), using
the same technique that we have used for spheres, since we did not exploit there any
special properties of spheres. The main tool that we used there was the existence
of a vertical decomposition of an arrangement of objects in 3-space into a slightly
super-cubic number of cells of constant complexity each. This can be done, e.g., for
any collection of objects whose boundary is defined in terms of a constant number
of algebraic equalities and inequalities of constant maximum degree [14]. Thus, in
general, this stage can be performed using O(n 3+" ) storage and preprocessing, and
O(log n) query time.
(ii) Objects that contain both endpoints of e in their interior. Since the objects are assumed
to be convex, we can ignore such objects because e cannot intersect any of them.
(iii) Objects that do not contain any endpoint of e. These objects require further treatment,
and are passed to the subsequent stages.
(2) Let S be an object of the third type, with respect to our query segment e, and suppose
that e intersects S. Since p lies outside S, there exists a plane through p that avoids S,
and it is clear that q and S lie on the same side of that plane (as a matter of fact, this
holds for any plane that passes through p and avoids S). A symmetric property holds at
q. Conversely, if there exist a pair of such planes, then, as is easily seen, e intersects S if
and only if ' intersects S. Thus at this stage we want to sift out all objects for which there
exist a pair of planes with these properties.
Here is a way (perhaps not always the most efficient) of doing this. Suppose that each
object S 2 S is defined by a polynomial inequality Q(x; lies outside S,
there is - ? 1 such that at p. The plane - p through p tangent to the surface
is a plane that contains S fully on one side. The equation of the halfspace bounded by - p
and containing S is
where are the coordinates of p. Since q also has to lie in this halfspace, we obtain
the constraint
where are the coordinates of q. A symmetric condition has to hold at q, namely
Both inequalities are linear in the coefficients of Q. Thus, if we represent each S 2 S as a
point in IR k , where k is the number of distinct monomials appearing in the partial derivatives
of the corresponding polynomials Q, the testing for the present condition amounts to two
halfspace range searching queries in IR k .
Figure

10: A disc appearing more than once along the upper convex hull
Of course, this method is not ideal if k is large, but is nevertheless a fairly general
technique. It works well for spheres: indeed, if the equation of a sphere S i 2 S is
then the above two constraints become
and
which, interestingly enough, are identical to the testing of the condition that the center of
lies inside the slab \Sigma e , as used in our previous algorithm.
(3) At this stage we have reduced our problem to that of detecting an intersection between
a query line ' and the objects of S (or, more precisely, of some canonical subset of S). We
next need to further reduce it to the problem of detecting an intersection between ' and
the planar upper convex hull of the cross sections of the objects of S in the vertical plane
passing through ', as we did for spheres. To do that we wish to partition S into two
subsets; none of the objects in the first subset should intersect V in a cross section that lies
fully above ', and none of the objects in the second subset should intersect V in a cross
section that lies fully below '. We have shown that this step can be performed for spheres
using an appropriate halfspace range searching on the set of centers of the spheres.
(4) In the last stage we need to decide whether our query line ' passes above the upper
convex hull of the cross sections of the objects of (a canonical subset of) S in the vertical
passing through '. However, our solution for this stage is fairly general, and can be
adapted to apply to objects other than spheres. Indeed, a close inspection of our technique
shows that the only place where we have used specific geometric properties of spheres is in
the definition of the maps M \Sigma (S i ), where each map records the spheres that are adjacent to
a sphere S i along the upper hull, provided that the intersection of such a sphere with V is
larger than the intersection of S i with V . This was done in order to overcome the technical
difficulty that the sphere following (or preceding) S i along the hull need not be unique, since
one intersection disc may appear more than once along the hull; see Figure 10. Here is a
way to overcome this difficulty if we assume that the objects of S are pairwise disjoint. In
this case, their cross sections in any vertical plane are pairwise disjoint convex regions. We
say that one such region D i is wider than another region D j if the length of the projection
of D i on the x-axis is larger than that of D j . We now define the maps M \Sigma (S i ) in complete
analogy to their definition in the case of spheres, with the proviso that a bridge connecting a
cross section D i with a cross section D j is recorded in the appropriate map of the narrower
of the two regions. It is easily verified that there can be at most one bridge connecting D i
with a wider D j which follows D i in the upper hull, and at most one bridge connecting D i
with a wider region preceding D j along the hull. Hence the maps are all well defined, and
we can proceed in exactly the same manner as for spheres. If the objects in S can intersect,
we need to devise an alternative method for recording all bridges on the upper hull in a
well-defined manner. Assuming that this is indeed possible, the remainder of our algorithm
applies to more general classes of objects, with only some obvious and simple modifications.
We only need to require that each object in S has constant description complexity; this
will suffice to show that the overall complexity of the maps M \Sigma (S i ) and (the vertically
This implies that this stage of the algorithm can be
performed on a set of n objects with O(n 3+" ) preprocessing and storage, and with O(n "
query time. It is interesting that the most complex stage of our algorithm in the case
of spheres is the most general stage, which applies, with comparable efficiency, to general
classes of objects.
We now combine all these stages into a multi-level data structure, as we did for spheres.
The preprocessing time and storage of the overall structure are roughly dominated by those
of the most 'expensive' stage, as already mentioned above and will be illustrated below.
We next illustrate how the various stages of the algorithm can be performed for two
specific classes of objects: triangles and ellipsoids.
3.1 The Case of Triangles
Most of the stages of the algorithm are easy to perform in the case of triangles. The first
stage is void, assuming non-degenerate position of the query segment. In the second stage
it suffices to find the set of all triangles such that the two endpoints of the query segment
are separated by the plane containing the triangle. This is easy to do using halfspace range
searching in 3 dimensions. The fourth stage applies in full generality, as we have noted, and
the technical difficulty concerning relative sizes of the intersection objects (line segments in
this case) is easy to handle; we leave details of this to the reader.
The most expensive stage, in the case of triangles, is the third one. Here we want to find
all triangles passing below (resp. above) a query line or intersecting it. For each triangle \Delta
we can find an edge s so that it suffices to determine whether the query line ' passes above
(resp. below) the line containing s. This in turn can be done by mapping the problem into
a halfspace range-searching problem in Pl-ucker 5-space, which can be done using O(n 4+" )
preprocessing and storage; see [15] for details. Thus the third stage is the bottleneck in
obtaining an efficient ray shooting method for triangles, and our approach does not seem
to yield in this case a better solution than those already known [5].
3.2 The Case of Ellipsoids
We next consider the case where our objects are ellipsoids whose axes are parallel to the
coordinate axes. Each ellipsoid S i 2 S can be represented by the equation A i
are positive.
As noted above, the first stage can be applied in a general setting, which includes the
case of ellipsoids. The second stage calls for the testing of the inequalities
and
These are homogeneous linear inequalities in the six parameters
and so can be performed in projective 5-space, using O(n 5+" ) storage and preprocessing,
and logarithmic query time.
Remark: Since this stage is the bottleneck in the case of ellipsoids, it is useful to observe
that we can make this stage more efficient if our ellipsoids satisfy additional properties. For
example, if they are all ellipsoids of revolution (about their vertical axis) then A
the above inequalities become homogeneous linear inequalities in projective 4-space, so they
can be done in O(log n) time, using only O(n 4+" ) storage and preprocessing.
The third stage can be carried out as follows. Suppose that the query line ' is represented
by the equations Take an ellipsoid S i and project it onto the
xy-plane, obtaining the ellipse S
1. The intersection of S
with the projection '   leads to the equation
or
The midpoint between the two intersection points (if they exist) thus satisfies
2:
It suffices to compare z then the cross-section S
cannot lie fully above ', and if z cross-section cannot lie fully below '; this
yields the classification of the ellipsoids that we seek.
The required comparison is equivalent to the testing of the sign of the expression
which is a linear homogeneous expression in the six coefficients
Thus the testing in this stage reduces to halfspace range searching in projective 5-space.
Moreover, since we are looking for a solution that uses large storage and small query time,
the problem becomes, in a dual setting, that of point location among an arrangement of
hyperplanes in 5-space. Since the query points lie on a 4-dimensional low-degree algebraic
surface (they are parametrized by the 4 parameters (- )), we can exploit the Zone
Theorem of Aronov et al. [9] to solve this problem using O(n 4+" ) storage and preprocessing,
with logarithmic query time (see [5] for more details).
The fourth stage is applicable in general setting, and the technical difficulty concerning
multiple bridges from the same cross-section can be overcome for pairwise disjoint ellipsoids,
as described above. Hence, putting it all together, we obtain:
Theorem 3.1 Ray shooting amidst n axis-parallel pairwise disjoint ellipsoids in 3-space
can be performed in O(n " ) time, using O(n 5+" ) preprocessing time and storage. This can
be reduced to O(n 4+" ) for ellipsoids of revolution.
In this paper we have developed a new and improved technique for ray shooting amidst
spheres in 3 dimensions. Our method requires O(n 3+" ) preprocessing time and storage,
and performs a ray shooting query in time O(n " ), for any " ? 0, improving previous
solutions by roughly an order of magnitude (in terms of preprocessing and storage costs).
We have also shown that our technique can be extended to obtain a general algorithm for
ray shooting amidst convex objects (of constant description complexity) in 3-space, and
have demonstrated this extension for the case of triangles (where the resulting algorithm is
no worse than known solutions), and for the case of axis-parallel ellipsoids, where we obtain
an algorithm that appears to be new, and yields fast queries, using O(n 5+" ) preprocessing
and storage.
The weakness of our algorithm is that it does not seem to yield a good trade-off between
storage and query time (in contrast with the less efficient solution of [4], which does have
such a trade-off). Of course, we can combine our solution with that of [4] to obtain some
tradeoff. Readers familiar with this technique can easily check that the resulting algorithm
yields query time of O(n 9=8+" =s 3=8 ) using O(s 1+" ) storage and preprocessing, for
It would be interesting (and we pose it as an open problem) to design an algorithm for ray-shooting
amidst spheres, which takes close to linear storage and can answer ray shooting
queries in time close to O(n 2=3 ).
Another set of open problems is to apply our general technique for various specific classes
of objects, attempting to find solutions to each stage that are as efficient as possible. In
particular, can the algorithm be improved for the case of ellipsoids to require only O(n 4+" )
storage?
Finally, no lower bounds are known for the ray shooting problem. Is our solution for
the case of spheres close to optimal in the worst case?


--R


Computing envelopes in four dimensions with applications

Computing 22
Range searching with semi-algebraic sets
Applications of a new space partitioning technique

Sharp upper and lower bounds for the length of general Davenport-Schinzel sequences
On the zone of a surface in a hyperplane arrangement
Good splitters with applications to ray shooting
Algorithmica 12
Lower bounds on the complexity of polytope range searching

A singly exponential stratification scheme for real semi-algebraic varieties and its applications
Lines in space: Combinatorics and algorithms
Visibility and intersection problems in plane geometry

Combinatorial complexity bounds for arrangements of curves and spheres
Determining the separation of preprocessed polyhedra: a unified approach
Arrangements of curves in the plane - Topology
related queries in arrangements of segments
New bounds for lower envelopes in three dimensions with applications to visibility of terrains

A pedestrian approach to ray shooting: Shoot a ray
Approximations and optimal geometric divide-and-conquer
Range searching with efficient hierarchical cuttings
Applying parallel computation algorithms in the design of serial algo- rithms
Almost tight upper bounds for lower envelopes in higher dimensions

--TR

--CTR
Boris Aronov , Mark de Berg , Chris Gray, Ray shooting and intersection searching amidst fat convex polyhedra in 3-space, Proceedings of the twenty-second annual symposium on Computational geometry, June 05-07, 2006, Sedona, Arizona, USA
Micha Sharir , Hayim Shaul, Ray shooting amid balls, farthest point from a line, and range emptiness searching, Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms, January 23-25, 2005, Vancouver, British Columbia
Micha Sharir , Hayim Shaul, Ray shooting and stone throwing with near-linear storage, Computational Geometry: Theory and Applications, v.30 n.3, p.239-252, March 2005
Pankaj K. Agarwal , Boris Aronov , Vladlen Koltun , Micha Sharir, On lines avoiding unit balls in three dimensions, Proceedings of the twentieth annual symposium on Computational geometry, June 08-11, 2004, Brooklyn, New York, USA
Vladlen Koltun, Almost tight upper bounds for vertical decompositions in four dimensions, Journal of the ACM (JACM), v.51 n.5, p.699-730, September 2004
Pankaj K. Agarwal , Micha Sharir, Efficient algorithms for geometric optimization, ACM Computing Surveys (CSUR), v.30 n.4, p.412-458, Dec. 1998
ray shooting;computational geometry