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

© Copyright 2018 AnyForm