Computing the Minimum Enclosing Sphere of Free

Computing the Minimum Enclosing Sphere of
Free-form Hypersurfaces in Arbitrary Dimensions
Ramanathan Muthuganapathy∗,
†
Gershon Elber†
Gill Barequet†
Myung-Soo Kim‡
Abstract
The problem of computing the minimum enclosing sphere (MES) of a point set is a classical
problem in Computational Geometry. As an LP-type problem, its expected running time on the
average is linear in the number of points. In this paper, we generalize this approach to compute
the minimum enclosing sphere of free-form hypersurfaces, in arbitrary dimensions. This paper
makes the bridge between discrete point sets (for which indeed the results are well-known) and
continuous curves and surfaces, showing that the general solution for the former can be adapted
for the latter.
To compute the MES of a pair of hypersurfaces, each one having a contact point (a point at
which the sphere touches the hypersurface), antipodal constraints are employed. For more than
a pair, equidistance constraints along with tangency constraints are applied. These constraints
yield a finite set of solution points which is used to identify the minimum enclosing sphere. The
algorithm uses the LP-characteristic of the problem to process the input set. Furthermore, an
optimization procedure that uses the convex hull of sampled points from the hypersurfaces is
also described. Finally, results from our implementation are presented.
Keywords: minimum enclosing circle, minimum enclosing sphere, minimum enclosing sphere,
LP-type problems, free-form hypersurfaces.
1
Introduction
Computing the minimum spanning/enclosing circle of a given set of points in the plane is considered
a classic problem in computational geometry [4]. (The original reference goes back to [14].) The
optimal solution (a minimum-radius circle) can be found in Θ(n) expected time, where n is the
number of points. In fact, the minimum enclosing sphere of n points in any dimension can be found
in expected Θ(n) time by using the same method. Sharir and Welzl [13] further generalized the
technique and introduced the notion of LP-type (linear-programming-like) problems. In a nutshell,
one can apply the LP-technique to an optimization problem in which adding a new constraint C does
not change the solution P, or otherwise the new solution P 0 is partially defined by C, thereby reducing
∗ Corresponding
author
for Graphics and Geometric Computing, Dept. of Computer Science, The Technion—Israel Institute
of Technology, Haifa 32000, Israel, Tel: +972-4-829-4338. The first author is currently in Department of Engineering Design, Indian Institute of Technology Madras, Chennai - 600036, India, Tel: +91-44-22574734. E-mail:
[email protected], [gershon|barequet]@cs.technion.ac.il
‡ School of Computer Science and Engineering, Seoul National University, Seoul 151-744, South Korea. E-mail:
[email protected]
† Center
1
the number of possible configurations of constraints defining P 0 . Amenta [2] later showed that special
properties of these problems imply so-called Helly-type theorems. Fischer et al. [8, 9] developed an
algorithm for computing the minimum enclosing ball of a set of N balls in a d-dimensional space
using the LP technique, while starting by explicitly addressing the small cases (N ≤ d + 1). They
appear to require rational arithmetic in their implementation. A deterministic algorithm for the
smallest enclosing ball of balls has been provided in [11, 12].
It has been shown that the computation of the minimum enclosing circle of a set of smooth planar
curves is possible to compute in expected linear time as it is also an LP-type problem [3] (by smooth
curves we mean piecewise C 1 polynomial or rational curves of bounded degree.). This type of problem
usually arises when there is a need to design efficient hierarchical spatial data structures for storing a
large amount of geometric objects. These data structures can then support fast clipping, searching,
collision detection, ray tracing, and other similar query-operations. For example, the sphere tree
in [10] supports a hierarchical modeling of complex 3-dimensional scenes using bounding spheres
of various 3-dimensional objects, and it has also been extensively used in 3-dimensional computer
games [1, 5].
In this paper, we show that the generalized problem, in which one aims to find the minimum
enclosing sphere in d dimensions of an arbitrary set of free-form hypersurfaces, is possible. Amenta’s
thesis [2] talks about relating Helly theorems with genearlized LP and convex programming methods
and gave some application algorithms on halfspaces, points (which can possibly be extended to compact sets), polygons/polytopes etc (refer [2, §9.1 to 11.3, pp. 64–94]). Though our work can be considered as based on the Helly theorems, we have also provided a methodology for computing the MES
of a set of closed curves/surfaces. Typically, these problems are solved for point-sets/cirles/ellipses.
The purpose of this paper is to show that the discrete approach (for point sets) can be adapted to
the continuous case, with the appropriate modifications and additional tools. The main contribution
of the paper is, we have shown that this problem can be solved (through the multivariate constraint
solver of Elber and Kim [7]) even if the input is of freeform curves/surfaces defined as Bezier/BSpline through an algebraic technique. Moreover, the tangent and equidistance conditions can be
generalized for any dimension.
1.1
Basic Idea
Given a set of N closed C 1 continuous hypersurfaces H1 , H2 , ...., HN , where a hypersurface can be
a point, curve, surface, or a freeform entity in more than two parameters in Rd , the problem is
to compute the minimum enclosing sphere (MES) of the set. Each Hi = Hi (ui1 , . . . , uin ), where
ui1 , . . . , uin are the parameters representing the hypersurface, with n <= d − 1. The MES is a
hypersphere that encloses the given set and is of minimum radius. For example, in 2D, the input set
2
Figure 1: A MES in two dimension. The input consists of points and free-form planar C 1 curves.
can consist of points and free-form curves (represented as H1 (u11 ), H2 (u21 ), . . .), and the problem
is to compute the minimum enclosing circle (Figure 1). For practical purposes, our implementation
assumes that the MES is tangent to any hypersurface in at most one point. We explain in Section 6
how to avoid this assumption. Moreover, the hypsersurfaces are assumed to be in general position
i.e., they overlap neither partially nor completely. In addition, it is assumed that no “spherical”
portions of any hypersurface lie completely on the MES.
The algorithm always starts with two input hypersurfaces. It is shown that, for the MES of two
hypersurfaces, the points on the hypersurfaces touching the MES have to be antipodal, i.e., they
lie in exactly opposite direction to each other with respect to the center of the MES. Based on the
fact that the points have to be antipodal, constraint equations are then formulated. The finite set
of solutions from these constraints are then used to compute the MES of the two hypersurfaces.
The rest of the input set is considered one by one (which is typical in an LP-type algorithm). Each
hypersurface is tested for its containment in the current MES. When three or more hypersurfaces
are involved, the algorithm uses the fact that their touching points have to be equidistant from the
center of the sphere, and the sphere has to be tangent to all of them. Hence, the constraints are
termed as equidistance and tangency constraints, respectively. Again, examining the finite set of
solutions from these constraints, the MES can then be computed.
It is a well-known fact that, when the input set is a point set, it is sufficient to use the points on
the boundary of the convex hull (hereafter, convex hull implies the boundary of the convex hull) of the
input set to compute its MES. An optimization procedure uses this fact to eliminate hypersurfaces
3
that lie in the interior of the convex hull of points sampled from the input data. It is also shown
that such an optimization procedure is indeed correct and speeds up the algorithm.
Our implementation presented here uses IRIT [6], and its multivariate piecewise-rational solver [7],
to solve sets of algebraic constraints. This paper is an extension of [3], a work that deals with the
computation of the minimum enclosing circle for inputs consisting only of planar curves.
The rest of the paper is organized as follows. The equidistance and tangency constraints for
d + 1 hypersurfaces in Rd are provided in Section 2. A description of constraints for less than
d + 1 hypersurfaces in Rd is described next. In particular, Section 3.1 describes the constraints
for two hypersurfaces and Section 3.2 describes the constraints for triplets in three dimensions.
The computation of the MES using the solution of the constraints is detailed in Section 4. Section 5
proves the LP characteristics of the problem. Details of the algorithm are described next in Section 6.
An optimization procedure which uses the convex hull of the input is described in Section 7. The
complexity of the algorithm is analyzed in Section 8. Experimental results are provided in Section 9,
and we terminate with some concluding remarks in Section 10.
2
Equidistance and Tangency Constraints for d + 1 Hypersurfaces in Rd
Let H1 , H2 , . . . , Hd+1 be d + 1 closed C 1 -continuous hypersurfaces in Rd (please note that Hi =
Hi (ui1 , . . . , uin ), where ui1 , . . . , uin are the parameters and n <= d − 1). Since the center of the
MES, C = (C1 , C2 , ...., Cd ), has to be equidistant from the touching points of all the hypersurfaces
(we assume, for now, that the MES is indeed tangent to all Hi , i = 1, ..., d + 1), the equidistance
constraints can be formulated as follows:
hC − H1 , C − H1 i =
hC − H2 , C − H2 i ,
hC − H1 , C − H1 i =
..
.
hC − H3 , C − H3 i ,
hC − H1 , C − H1 i =
hC − Hd+1 , C − Hd+1 i .
(1)
Equations (1) are equivalent to the following linear equations in C = (C1 , C2 , ...., Cd ):
2 hC, H1 − H2 i
= ||H1 ||2 − ||H2 ||2 ,
2 hC, H1 − H3 i
..
.
= ||H1 ||2 − ||H3 ||2 ,
2 hC, H1 − Hd+1 i
= ||H1 ||2 − ||Hd+1 ||2 .
Noting that the MES also has to be tangent to the C 1 -continuous closed hypersurfaces it touches,
4
Degree of Freedom (DOF)
d coefficientsP
of center of MES
Di
Number of Constraints
d
Pequidistance constraints
Di tangency constraints
Table 1: DOF vs. the number of constraints for the case of d + 1 hypersurfaces in Rd , where Di is
the dimension of the hypersurface Hi , i = 1, ..., d + 1
.
the tangency constraints can be written as follows:
hC − H1 , dH1 i
= 0,
hC − H2 , dH2 i
..
.
= 0,
hC − Hd+1 , dHd+1 i
= 0.
(2)
where dHi denotes a basis of vectors that spans the tangent vectors of Hi . Solving the equidistance
and tangency constraints results in a finite set of solutions for C which can be represented using
the parameters of the hypersurfaces. It is easy to see that this set of constraints is well-defined and
results in a zero dimensional solution space. It should be noted that, though the closed C 1 -continuous
hypersurfaces might have singular locations because of parameterization, we do assume regularity of
the constraints when solving for their zeros. Nonetheless, the probability of a tangency/contact at a
singularity, in a general arrangement, is close to zero. Practically, these cases could be dealt with as
special cases.
Let Di = dim(Hi ) be the dimension of each hypersurface. First, the equidistance constraints (1)
P
result in d equations. Second, the number of tangency constraints (2) is
Di , the total sum of the
dimensions of the hypersurfaces. On the other hand, we have d unknowns of the center C of the
P
sphere and
Di unknowns in all Hi hypersurfaces. Table 1 shows the degrees of freedom available
P
versus the number of constraints, d and
Di .
The following example shows how the equidistance and tangency constraints can be employed in
practice.
Example 2.1 Equidistance and Tangency Constraints for Three Curves in R2
Given three planar curves H1 (u11 ) (= C1 (t)), H2 (u21 ) (= C2 (r)), and H3 (u31 ) (= C3 (s)), their
equidistance constraints from the center of the MES, C = (Cx , Cy ) (see Figure 2 and Equations (1))
can be written as:
hC − C1 (t), C − C1 (t)i =
hC − C2 (r), C − C2 (r)i ,
hC − C1 (t), C − C1 (t)i =
hC − C3 (s), C − C3 (s)i .
5
(3)
C1 (t)
C3 (s)
C
C2 (r)
Figure 2: Constraints for three free-form planar curves in R2 .
Equations (3) are equivalent to the following linear equations in C = (Cx , Cy ):
2 hC, C1 (t) − C2 (r)i
= ||C1 (t)||2 − ||C2 (r)||2 ,
2 hC, C1 (t) − C3 (s)i
= ||C1 (t)||2 − ||C3 (s)||2 .
(4)
The tangency constraints (see Equation (2)) for the three curves can be written as:
hC − C1 (t), C10 (t)i =
0,
hC − C2 (r), C20 (r)i = 0,
hC − C3 (s), C30 (s)i =
0.
(5)
The equidistance constraints (Equations (4)) and tangency constraints (Equations (5)) for three
curves in R2 result in five equations, in as many unknowns, t, r, s, and C = (Cx , Cy ).
3
Constraints for Less than d + 1 hypersurfaces (having less
than d+1 contact points) in d Dimensions
When we have less than d+1 hypersurfaces in d dimensions, the equidistance and tangency constraints
toward the MES result in an under-constrained problem. For example, for two curves in 2D, there
is one equidistance constraint and two tangency constraints, resulting in three equations in four
unknowns. Yet, clearly there is a unique MES.
6
3.1
Antipodality Constraints for Two hypersurfaces
Lemma 1 A sphere, centered at C in Rd enclosing two different hypersurfaces (each having a contact
point) is minimum in radius if and only if the normals of the corresponding contact points are opposite
to each other with respect to C (i.e., antipodal.)
Proof: By contradiction. Assume that the unit vectors from the center C to the contact points on
the hypersurfaces are V1 and V2 and are not antipodal. There exists a small distance ε > 0 so that
the sphere can be translated in the direction of V1 + V2 , then the sphere has no contact point with
the two hypersurfaces though it still bounds them. We can then shrink the sphere slightly so as to
find an enclosing sphere which is smaller in radius. Hence, V1 and V2 have to be antipodal.
2
The above antipodality constraint can be formulated in terms of equations. Let H1 and H2 be
the two hypersurfaces. The antipodality constraint can be written as (noting that:
¿
À
H1 + H2
dH1 , H1 −
= 0,
2
¿
À
H1 + H2
dH2 , H2 −
= 0,
2
H1 +H2
2
= C):
or, equivalently,
hdH1 , H1 − H2 i = 0,
hdH2 , H2 − H1 i = 0.
(6)
Equations (6) represent the antipodality constraint in any dimension d. Again, dHi denotes
a basis of vectors that spans the tangent vectors of Hi and is not necessarily prescribing a single
constraint. The following two examples show the formulation of antipodal equations for two curves
in R2 and two surfaces in R3 , respectively.
Example 3.1.1 Antipodal constraint of two curves in R2
Consider two planar closed C 1 -continuous curves H1 (u11 ) (= C1 (t)) and H2 (u21 ) (= C2 (r)). The
antipodality constraint (Equations (6)) is:
¿
À
C1 (t) + C2 (r)
0
C1 (t), C1 (t) −
2
¿
À
C
(t)
+
C2 (r)
1
0
C2 (r), C2 (r) −
2
=
0,
=
0.
(7)
Clearly, the two equations have two unknowns (t and r). Figure 3 illustrates the antipodality
constraint for two planar curves.
Example 3.1.2 Antipodality constraint of two surfaces in R3
7
C1 (t)
C
C2 (r)
Figure 3: Antipodality constraint for two free-form planar curves in R2 .
S1 (u1 , v1 )
C
S2 (u2 , v2 )
Figure 4: Antipodality constraint for two free-form surfaces in R3 .
8
H3
C
H1
H2
Figure 5: Constraints for three free-form surfaces in R3 .
Consider two closed C 1 -continuous surfaces H1 (u11 , u12 ) (= S1 (u1 , v1 )) and H2 (u21 , u22 ) (=
S2 (u2 , v2 ) (Figure 4). Then, each of the antipodality constraint is represented by two equations.
À
¿
S1 (u1 , v1 ) + S2 (u2 , v2 )
∂S1 (u1 , v1 )
, S1 (u1 , v1 ) −
= 0,
∂u1
2
¿
À
∂S1 (u1 , v1 )
S1 (u1 , v1 ) + S2 (u2 , v2 )
, S1 (u1 , v1 ) −
= 0,
∂v1
2
¿
À
∂S2 (u2 , v2 )
S1 (u1 , v1 ) + S2 (u2 , v2 )
, S2 (u2 , v2 ) −
= 0,
∂u2
2
¿
À
∂S2 (u2 , v2 )
S1 (u1 , v1 ) + S2 (u2 , v2 )
(8)
, S2 (u2 , v2 ) −
= 0,
∂v2
2
Thus, in (8) we have four equations in four unknowns (u1 , v1 , u2 and v2 ) in all. Note that in
the antipodal cases, C is easily deduced as (S1 + S2 )/2. For the antipodal condition, the major
degeneracy appears to be when the hypersurfaces having hyperplane portions that are parallel to
each other. For example, when the input involves curves, this translates to having straight line
portions in the curves. The MES though will have antipodal points that do not lie on the straight
line portions.
3.2
Constraints for Triplets in Three and Higher Dimensions
Consider H1 , H2 , and H3 (Figure 5). The equidistance constraints:
hC − H1 , C − H1 i =
hC − H2 , C − H2 i ,
hC − H1 , C − H1 i =
hC − H3 , C − H3 i .
9
(9)
The corresponding tangency constraints are:
hC − H1 , dH1 i
= 0,
hC − H2 , dH2 i
= 0,
hC − H3 , dH3 i
= 0.
(10)
Let us examine the set of constraints assuming that the given hypersurfaces are surfaces (in two
parameters). Clearly, the three tangency Equations (10) become six equations, taking into account
the two parameters of each surface. This results in a total of eight equations (two equidistance
constraints and six tangency constraints) in nine unknowns (the center point for C and the parameters
of each surface), which is an under-constrained problem. The remaining equation comes from the
following lemma.
Lemma 2 A sphere enclosing three different hypersurfaces (each having a contact point) in R3 is
minimum in radius if and only if the vectors from the center of the sphere to the contact points on
the hypersurfaces lie in a plane, i.e, the vectors are coplanar.
Proof: By contradiction. Let the vectors from the center C of the MES to the contact points on
the hypersurfaces be V1 , V2 and V3 , ||Vi || > 0, and assume they are not coplanar and let L denote
the plane determined by the three contact points. (Note that the center C is not contained in the
plane L.) There exists a small distance ε > 0 so that the sphere can be translated toward the plane
L along the normal direction of the plane. The sphere will have no contact point with the three
hypersurfaces though it still bounds them. We can then shrink the sphere so as to find a smaller
enclosing sphere. Hence, the vectors have to lie on a plane.
2
The above is the generalization of the two hypersurfaces’ antipodality constraint (Lemma 1).
Corollary 1 The three contact points of the MES of the three hypersurfaces in R3 define a great
circle on the MES.
As the center C and the contact points of the MES (which is a sphere in R3 ) are on a plane (by
Lemma 2, Corollary 1 follows.
Writing the coplanarity constraint in Lemma 2 as a ninth constraint results in
­
®
(C − H1 ) × (C − H2 ), (C − H3 ) = 0.
(11)
Equations (9), (10) and (11) completely determine C in R3 and the parameters of the hypersurfaces, having nine equations in nine unknowns.
10
For m hypersurfaces in Rd and m < d + 1, the generalization is obtained from the following
lemma (Lemma 3).
Lemma 3 The contact points of the MES of m different hypersurfaces in Rd , m < d + 1, are all on
a sphere of dimension m − 1.
Proof: Consider Equations (9) and (10). The tangency Equations (10) determine the parameters of
the hypersurfaces (see Table 1). The dimension of the sphere or the number of coordinates of C (the
sphere defined by the center C and the radius defined by the center C and a contact point) depend
on Equations (9). If there are m hypersurfaces, then there are m − 1 equidistance constraints. This
implies that there are m − 1 coordinates for the sphere and hence it is of dimension m − 1. This
sphere is also minimum in radius, otherwise, it violates Corollary 1 (though Corollary 1 is for R3 ,
the generalization to Rd is rather straightforward).
4
2
Computing the MES
The solution of the constraint equations (detailed in Sections 2 and 3) results in a finite set of
candidates. Each solution defines a sphere which is tangent to all the hypersurfaces, although some
of the tangencies may be to the outside of the sphere as is clearly undesired. For example, in Figure 6,
the solutions B1 , B2 , and B3 are incorrect, while we seek the solution B4 . To identify the correct
MES from the set of candidate solution spheres, the hypersurfaces in consideration are subjected to
a containment test – of whether or not the hypersurface is contained in the sphere, i.e., all the points
on the hypersurface lie within (or on the boundary of) the sphere. Obviously, only the MES will
contain all the hypersurfaces.
Consider a hypersurface H and a sphere (C, R), centered at C and of radius R (see Figure 7). H
is inside the sphere if and only if
||Dm − C|| ≤ R,
where Dm is the maximum distance point of H from the center point C. While one can indeed
compute all radially-extreme locations from H to C, the following would be examined.
F = ||H − C||2 − R2 ≤ 0.
(12)
Given a piecewise rational hypersurface, H, Equation (12) is also rational, and hence one can
examine all the coefficients of F . If all the coefficients are nonpositive, then, due to the convex-hull
property of the NURBS representation, F must be nonpositive as well. Otherwise, one needs to
compute the zero set of F . If the zero-set is empty, then H is either completely inside or completely
outside the sphere, in which case a single evaluation of a point in H could reveal its inclusion status.
11
B3
B1
B4
B2
Figure 6: Constraints yield multiple sets of solution points.
Dm
H
C
R
Figure 7: Containment test.
12
Q
B0
S
C1
C0
B1
Figure 8: LP-type characteristics of the problem (see Theorem 1). B0 =MES(S) and B1 =MES(S ∪
{S}). Q is an intersection point of B0 and B1
While this test does not affect the complexity of the algorithm, in practice it can greatly save a lot
of running time.
5
LP Characteristics of the Problem
We now consider the computation of the MES of an arbitrary set of C 1 smooth surfaces as an LP-type
problem. This is merely a generalization of Theorem 1 in [3].
Toward this end, we need to prove that given the MES of a set of surfaces, if a newly introduced
surface is entirely or partially outside the current MES, then the MES of all the surfaces (including
the new one) is partially defined by the new surface. In other words, the new MES must be in contact
with the new surface.
More formally, we prove the following:
Theorem 1 Let S be a set of C 1 smooth surfaces, and let MES(S) be the MES of S. Also, let S 6∈ S
be a newly introduced surface. Then, either
1. The surface S is completely inside (possibly touching) MES(S), in which case MES(S ∪ {S}) =
13
M ES(S); or
2. Surface S is partially or completely outside MES(S), in which case S must be in contact with
MES(S ∪ {S}).
Proof: The proof follows step by step the proof of Lemma 4.14 in [4], which makes the same assertion
for points. Consider two possible cases:
1. MES(S) fully contains S. Assume to the contrary that some sphere B whose radius is less than
that of MES(S) contains S ∪{S}. In particular, B contains S, which contradicts the minimality
of MES(S).
2. S is partially (or fully) outside MES(S). Clearly, MES(S ∪ {S}) 6=MES(S) and the two spheres
must intersect; see Figure 8. Let B0 =MES(S) and B1 =MES(S ∪ {S}) with the respective
center points C0 and C1 . Denote by Q some point on the intersection of the boundaries of B0
and B1 . Define a continuous deformation between B0 and B1 as follows: the center of Bt (for
0 ≤ t ≤ 1) is Ct = (1 − t)C0 + t(C1 ), and its radius is D(Ct , Q), where D(·, ·) is the Euclidean
distance. Since S is partially or fully outside B0 , and S is fully contained in B1 , by continuity
there must be some 0 < t∗ ≤ 1 such that Bt∗ contains S and touches it. The set S is contained
entirely in both spheres B0 and B1 , therefore, S is contained in B0 ∩ B1 , which, by definition,
in turn is contained in Bt (for any 0 ≤ t ≤ 1). Also, the radius of any sphere Bt (for any
0 < t < 1) is strictly less than that of B1 . Therefore, we must have t∗ = 1; otherwise, we have
a contradiction to the minimality of B1 .
2
Remark: The proof of Theorem 1 holds for any dimension. The uniqueness of the MES of the
set S is shown by a similar argument of a continuous family of spheres. The reader is referred to [4,
§4.7, Lemma 4.14] for the exact details.
6
The Algorithm
After establishing the LP characteristics of the problem, we apply the well-known machinery to solve
it in an expected time that is linear in n, the number of hypersurfaces, on the average. We use the
randomized algorithm of [4, §4.7, pp. 85–88] that originally finds the minimum enclosing circle of a
planar point set. Assume for the moment that the MES of the hypersurfaces is defined by points that
belong to different hypersurfaces. The algorithm (described in Figure 9) outlines the computation
of the MES in three dimensions. In order to generalize it to d dimensions, one needs to extend the
algorithm to have d levels (possibly by implementing all levels in one function).
14
Algorithm MinSphere (S)
Input: A set S of n hypersurfaces H1 , H2 , ... , Hn in R3 .
Output: The minimum-radius sphere that fully contains S.
begin
Compute a random permutation H1 , . . . , Hn of the hypersurfaces in S.
Let B0 = B1 = B2 be the smallest sphere enclosing H1 and H2 .
for i = 3, . . . , n do
if Hi ⊂ Bi−1
then Bi := Bi−1 ;
else Bi := MinSphereWithOneHS ({H1 , . . . , Hi−1 }, Hi ).
end if
end for
return Bn .
end MinSphere
Function MinSphereWithOneHS (S, Q)
Input: A set S of N hypersurfaces and a hypersurface Q, such that there exists an
enclosing sphere of S that touches Q.
Output: The minimum-radius sphere that fully contains
S and that touches Q.
begin
Compute a random permutation H1 , . . . , HN of the hypersurfaces in S.
Let B0 = B1 be the smallest sphere enclosing H1 and Q.
for i = 2, . . . , N do
if Hi ⊂ Bi−1
then Bi := Bi−1 ;
else Bi := MinSphereWithTwoHSs ({H1 , . . . , Hi−1 }, Q, Hi ).
end if
end for
return Bn .
end MinSphereWithOneHS
Function MinSphereWithTwoHSs (S, Q1 , Q2 )
Input: A set S of n hypersurfaces and two hypersurfaces Q1 , Q2 , such that there exists an
enclosing sphere of S that touches Q1 and Q2 .
Output: The minimum-radius sphere that fully contains
S and that touches Q1 and Q2 .
begin
Compute a random permutation H1 , . . . , Hn of the hypersurfaces in S.
Let B0 be the smallest sphere enclosing {Q1 , Q2 }.
for i = 1, . . . , n do
if Hi ⊂ Bi−1
then Bi := Bi−1 ;
else Bi := MinSphereWithThreeHSs ({H1 , . . . , Hi−1 }, Q1 , Q2 , Hi ).
end if
end for
return Bn .
end MinSphereWithTwoHSs
Function MinSphereWithThreeHSs (S, Q1 , Q2 , Q3 )
Input: A set S of n hypersurfaces and three hypersurfaces Q1 , Q2 , Q3 , such that there exists an
enclosing sphere of S that touches Q1 and Q2 , Q3 .
Output: The minimum-radius sphere that fully contains
S and that touches Q1 , Q2 and Q3 .
begin
Compute a random permutation H1 , . . . , Hn of the hypersurfaces in S.
Let B0 be the smallest sphere enclosing {Q1 , Q2 , Q3 }.
for i = 1, . . . , n do
if Hi ⊂ Bi−1
then Bi := Bi−1 ;
else Bi := the minimum enclosing sphere of Q1 , Q2 , Q3 , and Hi .
end if
end for
return Bn .
end MinSphereWithThreeHSs
15
Figure 9: Computing the minimum-radius enclosing sphere of a set of hypersurfaces (following closely
the algorithm of [4] for enclosed points) in R3
In the upper level of the algorithm, we iteratively compute the minimum enclosing sphere of the
first i entities, where i goes from 2 to n (see function MinSphere in Figure 9). In the ith step, we
check whether the hypersurface Hi lies inside Bi−1 , the minimum enclosing sphere of the first i − 1
hypersurfaces. If so, then Bi = Bi−1 . Otherwise, if Hi lies partially or entirely outside Bi−1 , we need
to re-compute Bi , but now it is guaranteed by Theorem 1 that Bi touches Hi .
Accordingly, we now invoke a secondary function that performs the same task (namely, finding the
minimum enclosing sphere of a set of hypersurfaces), with the only restriction that the sought-after
sphere touches one specific hypersurface Q (see function MinSphereWithOneHS in Figure 9). Again,
we add one hypersurface at a time, and check whether the newly-added hypersurface is contained in
the previously-computed sphere. If this is not the case, we need to re-compute the sphere, but this
time it is guaranteed that the new MES touches both Q and the newly-added hypersurface.
The algorithm then invokes the function MinSphereWithTwoHSs and possibly the function MinSphereWithThreeHSs (Figure 9), with the respective restrictions that the sought-after sphere touches
three or four hypersurfaces. These two functions have the same logic as above. Note that we add one
hypersurface at a time while performing these functions. At the bottom level we need the services
of a function that computes the MES of four hypersurfaces.
We do not need to alter the algorithm to accommodate cases in which the MES is defined by two
or more points that lie on the same entity. The invariant of Theorem 1 remains true. That is, if the
ith entity lies (partially or entirely) outside the MES of the first through the (i − 1)st entity, then the
new MES must touch the ith entity. The algorithm only needs to consider all possible multi-touching
cases. For example, in the case of two points lying on the same entity, the constraint equations can be
solved considering H1 = H2 . To accommodate all possible multi-touching situations, the algorithm
must check all cases in increasing order of complexity. The correctness of the algorithm follows
directly from the LP characteristics of the problem.
7
Optimization Using the Convex Hull
It is well known that, for a given point set in the plane, the MES is entirely determined by the convex
hull (CH) of the point set. This is because the points defining the MES are always on the convex
hull of the set.
In higher dimensions, when free-form curves, surfaces, hypersurfaces etc., are involved, an optimization procedure can be applied. Such a procedure will eliminate those hypersurfaces that do
not contribute to the convex hull of the set. Since the convex hull of points in Rd is straightforward
to compute, points lying on all hypersurfaces are sampled and their convex hull is computed. The
optimization procedure, is the following:
16
C
Figure 10: The optimization procedure using the convex hull.
• Compute the convex hull of the sampled points (along with any input points).
• Eliminate all the hypersurfaces that are inside the convex hull. A simple test for this purpose
is to examine whether or not all the control points of the hypersurface are inside the hull.
• The remaining hypersurfaces are used to compute the MES using the algorithm described in
Section 6.
Figure 10 illustrates the optimization procedure. In the figure, the points sampled on the curves
are shown in grey, whereas the input points are shows in black. Also shown is, the convex hull of the
sampled and input points. Applying the test for containment of hypersurfaces in the convex hull, all
the input points and the curve C are eliminated in this example. The remaining hypersurfaces are
used for computing the MES.
Again, the correctness of the optimization procedure follows from the fact that CH(sampled +
input points) ⊂ CH(the input set) ⊂ MES(the input set).
17
8
Complexity Analysis
Denote by N the number of hypersurfaces in d dimensions, d fixed, and by p the maximum degree
of the polynomials describing the hypersurfaces. Let I(p) be the time needed to decide whether or
not a hypersurface is fully contained in a given sphere, and M (p) be the time needed to compute the
MES of d + 1 such hypersurfaces.
The main part of the algorithm, MinSphere, performs N iterations, in each of which it either
decides in I(p) time whether or not the MES (computed so far) has to change. Occasionally, it also
calls the function MinSphereWithOneHS. In the worst case, the main algorithm can call the latter
Θ(N ) times, in case all of the second through the N th hypersurfaces require an update of the MES.
Similarly, the function MinSphereWithOneHS performs O(N ) iterations. In each iteration it decides
in I(p) time whether or not to re-compute the MES by calling the function MinSphereWithTwoHSs.
This also can occur in O(N ) iterations. Its running time is, therefore, O(N (I(p) + M (p))). A similar
argument holds for the functions MinSphereWithTwoHSs and MinSphereWithThreeHSs. Overall, in
the worst case, the entire algorithm requires Θ(N 4 (I(p) + M (p))) time. Note that I(p) and M (p),
the practical running times of numerical methods, do not depend at all on N , but solely on p and on
the convergence parameters. When the latter are fixed or assumed constant, we have that both I(p)
and M (p) are low-degree polynomials in p. If we further fix p, then I(p) and M (p) can be regarded
as constants. Even if we take into consideration endpoints and C 1 -discontinuity areas of the entities,
the asymptotic running time of the algorithm does not change. Hence, the worst-case running time
of a practical implementation of our algorithm will be Θ(N 4 ) in three dimensions. It is easy to see
that the worst-case running time in d dimensions will be Θ(N d+1 ).
The average case is much more favorable. A straightforward analysis, that follows closely the
analysis in [4, §4.7, pp. 85–89], shows that the expected number of invocations of all functions is only
Θ(N ). Assuming that both I(p) and M (p) are constant, one obtains that the expected running time
of the entire algorithm is also Θ(N ), which is optimal.
The amount of space used by the algorithm is clearly Θ(N p), which is the space needed to store
all the hypersurfaces. At all times the algorithm holds no more than a constant number of spheres
needed to decide whether the MES of four hypersurfaces touches one, two, three, or all the four of
them.
9
Experimental Results
The algorithm described in Figure 9 has been implemented for computing the MES in two and three
dimensions. Optimization using the convex hull has also been implemented in two dimensions. The
18
(a)
(b)
Figure 11: The MES of points and free-form rational C 1 closed planar curves.
implementation was done with the aid of the IRIT [6] geometric kernel.
Here are two representation examples. Figure 11 shows two test results for the MES computation
in two dimensions. The input consists of points and free-form planar curves. Figure 12 shows some
test results for computing the MES in three dimensions, where the input consists of points, and
free-form curves and surfaces.
9.1
C 1 -discontinuous Hypersurfaces
The discussion so far has been for hypersurfaces which are assumed to be C 1 continuous. If C 1
discontinuous hypersurfaces are present, then their support is straightforward. It involves the introduction of low-dimensional hypersurfaces where the discontinuity occurs in the input set and
splitting the original hypersurface along the discontinuity into two pieces. For example, consider C 1
discontinuity points in the case of planar free-form curves. Figure 13 shows the MES computed for
an input set having C 1 discontinuous planar curves.
9.2
Optimization Procedure
Figure 14 shows some test cases after applying the convex hull optimization procedure described in
Section 7. Thicker points and curves are retained for computing the MES whereas the thinner are
the eliminated ones. Table 2 shows the optimization results in terms of statistical data. Note that
the number of points and curves that are eliminated are quite large and hence there is a significant
19
(a)
(b)
(c)
(d)
Figure 12: The MES in three dimensions. The surfaces and the open curves in the figure are of third
order.
20
Figure 13: The MES of points and planar free-form curves having C 1 discontinuities.
Curves
(a)
(b)
(c)
(d)
(e)
(f)
16
26
56
96
46
50
Degree
of
Curves
3
3
4
3
3
3
Control
Points
Input
points
5
5
7
5
3
3
15
100
50
1
11
11
Curves
eliminated
in opt.
9
17
32
61
24
30
Points
eliminated
in opt.
11
94
48
1
11
9
Running
time (sec.)
with opt.
2.0
2.53
16.7
0.84
6.42
4.9
Running
time (sec.)
w/o opt.
2.1
10.5
16.9
6.9
20.6
12.9
Table 2: Statistics of the test cases shown in Figure 14
change in the running times (the last two columns in the table give the running times with and
without optimization). Figure 14(a) also shows the convex hull of the input points and the points
sampled on the curves.
9.3
Large Data Sets
We have also experimented with the algorithm on several large data sets (of up to 40,000 curves),
indicating the robustness of the implementation. Figure 15(a) and 15(b) show the computed MESs
of 1,000 and 2,000 curves, respectively. The large data sets contain only curves. Table 3 shows the
average time taken to compute MES for all data sets (Each experiment was repeated five times.).
Figure 16 plots the running time as a function of the number of curves. It appears experimentally
that the running time is roughly linear in the number of curves when the convex-hull optimization is
not employed. However, when the optimization is performed, the average running time may not be
linear as the procedure will eliminate an arbitrary number of curves (which depends on the specific
21
(a)
(c)
(e)
(b)
(d)
(f)
Figure 14: Test cases showing the effect of using the convex-hull optimization procedure. Convex
22
hull is also shown in Figure 14(a).
(a) MES of 1,000 curves
(b) MES of 2,000 curves
Figure 15: Test cases for large data set without using optimization.
Number of curves
Avg Running Time (w/o opt.) (sec.)
Avg Running Time (with opt.) (sec.)
1,000
48.2
20.3
5,000
73.2
14.96
10,000
105.9
39.37
15,000
130.2
50.64
20,000
171.6
71.4
30,000
223.5
79.18
40,000
362.4
183.9
Table 3: Average running time for large data sets (curves).
set of curves).
10
Conclusion
In this paper, we presented an optimal-time (expected) solution for the problem of computing the
minimum spanning (enclosing) sphere of a set of points, free-form curves, surfaces, etc. It has been
shown that the theory holds for any dimensions, though our implementation was for two and three
dimensions only. The solution is inherently the same as for computing the minimum enclosing circle
(or sphere) of a set of points; however, the implementation details of our solution are significantly
more complicated. The expected running time of the algorithm is linear in the number of input
hypersurfaces, but also contingent on the other details, such as the maximum degree of the hypersurfaces. It has been shown that the well-known results for MES of point sets can be adapted to
continuous curves, surfaces, and hypersurfaces. In future, we also seek to overcome the assumption
that the MES is tangent to a hypersurface at most once. Finally, a similar algorithm could be sought
for hierarchical minimum enclosing computations, such as minimum enclosing ellipsoids, a problem
we are now investigating.
23
400
With Optimization
No Optimization
Average running time (sec)
350
300
250
200
150
100
50
0
0
0.5
1
1.5
2
2.5
Number of curves
3
3.5
4
4
x 10
Figure 16: Plot of the number of curves vs ṫhe average running time
Acknowledgements
Work on this paper has been supported in part by the Israel Science Foundation Grant 346/07, and
in part by a grant for Korean-Israeli Research Cooperation. Work on this paper by the first author
has been supported in part by a Lady Davis Fellowship, Israel. This work was also supported in
part by KICOS through the Korean-Israeli Binational Research Grant (K20717000006) provided by
MEST in 2007.
References
[1] Tomas Akenine-Möller, Eric Haines, and Natty Hoffman. Real-Time Rendering 3rd Edition. A.
K. Peters, Ltd., Natick, MA, USA, 2008.
[2] Annamaria Beatrice Amenta. Helly-type theorems and generalized linear programming. Discrete
Comput. Geom, 12:241–261, 1994.
[3] Gill Barequet, Gershon Elber, and Myung soo Kim. Computing the minimum enclosing circle
of a set of planar curves. Computer Aided Design and Applications, 2:301–308, 2005.
24
[4] Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars. Computational Geometry: Algorithms and Applications. Springer, 3rd edition, April 2008.
[5] David H. Eberly. 3D Game Engine Design, Second Edition: A Practical Approach to Real-Time
Computer Graphics (The Morgan Kaufmann Series in Interactive 3D Technology). Morgan
Kaufmann Publishers Inc., San Francisco, CA, USA, 2006.
[6] Gershon Elber. IRIT 10.0 User’s Manual. http://www.cs.technion.ac.il/~irit, 2009.
[7] Gershon Elber and Myung-Soo Kim. Geometric constraint solver using multivariate rational
spline functions. In SMA ’01: Proceedings of the sixth ACM symposium on Solid modeling and
applications, pages 1–10, New York, NY, USA, 2001. ACM.
[8] Kaspar Fischer and Bernd Gartner. The smallest enclosing ball of balls: combinatorial structure
and algorithms. In SCG ’03: Proceedings of the nineteenth annual symposium on Computational
geometry, pages 292–301, New York, NY, USA, 2003. ACM.
[9] Kaspar Fischer, Bernd Grtner, and Martin Kutz. Fast smallest-enclosing-ball computation
in high dimensions. In Proc. 11th European Symposium on Algorithms (ESA, pages 630–641.
Springer-Verlag, 2003.
[10] Philip M. Hubbard. Approximating polyhedra with spheres for time-critical collision detection.
ACM Trans. Graph., 15(3):179–210, 1996.
[11] Nimrod Megiddo. Linear-time algorithms for linear programming in r3 and related problems.
SIAM Journal on Computing, 12(4):759–776, 1983.
[12] Nimrod Megiddo. Linear programming in linear time when the dimension is fixed. J. ACM,
31(1):114–127, 1984.
[13] Micha Sharir and Emo Welzl. A combinatorial bound for linear programming and related
problems. In STACS ’92: Proceedings of the 9th Annual Symposium on Theoretical Aspects of
Computer Science, pages 569–579, London, UK, 1992. Springer-Verlag.
[14] Emo Welzl. Smallest enclosing disks (balls and ellipsoids). In Results and New Trends in
Computer Science, pages 359–370. Springer-Verlag, 1991.
25