Efficient Distance Computation for Quadratic ... - Semantic Scholar
Efficient Distance Computation for Quadratic Curves and Surfaces Christian Lennerz Elmar Sch¨omer Max-Planck-Institute for Computer Science 66123 Saarbr¨ucken, Germany {lennerz, schoemer}@mpi-sb.mpg.de
Abstract Virtual prototyping and assembly planning require physically based simulation techniques. In this setting the relevant objects are mostly mechanical parts, designed in CADprograms. When exported to the prototyping and planning systems, curved parts are approximated by large polygonal models, thus confronting the simulation algorithms with high complexity. Algorithms for collision detection in particular are a bottleneck of efficiency and suffer from accuracy and robustness problems. To overcome these problems, our algorithm directly operates on the original CAD-data. This approach reduces the input complexity and avoids accuracy problems due to approximation errors. We present an efficient algorithm for computing the distance between patches of quadratic surfaces trimmed by quadratic curves. The distance calculation problem is reduced to the problem of solving univariate polynomials of a degree of at most 24. Moreover, we will identify an important subclass for which the degree of the polynomials is bounded by 8.
1. Introduction During the past few years an increasing interest in virtual reality (VR) techniques could be observed. Virtual reality plays an increasingly important part in social and business life, ranging from entertainment and education to VR-based training methods and numerous industrial applications. However, the ‘classical’ challenges that come from concepts like virtual prototyping and assembly planning still represent strong drivers for the development of VR-techniques. To evaluate industrial prototypes and assembly processes, virtual environments have to simulate physical behavior. Particularly collisions between moving objects must be detected and resolved. Robustness and efficiency of such a simulation heavily depend on the quality of the collision detection algorithm. Dynamic collision detection routines are required in this context. They identify the earliest time
of collision in a given time interval and return the touching pair of objects. In contrast to this technique, static detection algorithms only decide whether two objects intersect at a given point in time. Since they check discrete time frames they, however, cannot guarantee to find the earliest contact. In [9] a framework is developed in which static distance computation algorithms can be used to design dynamic detection routines. The basic idea is that information about the distance and dynamic states of all objects provides lower bounds on the earliest time of collision. The problem of distance computation has been well studied in the past for polyhedral objects, e.g. [5, 8, 3], but only little effort has spent on objects with curved surfaces [11]. There are two reasons for this situation. First, a polygonal representation is not considered a true restriction since real objects can be approximated arbitrarily precisely by a polyhedron. The second reason is the fact that the basic algorithms and predicates can be implemented robustly and very efficiently on polygons. However, one has to consider that there is a trade-off between the accuracy of the approximation and the efficiency of the distance computation. On the one hand, precise approximations increase the complexity of the boundary representation and can affect the running time of the algorithm disproportionately. On the other hand, a small number of polygons compromises the quality or even the correctness of the simulation. Our considerations therefore raise the question, whether or not it is worth to handle curved objects directly. It means that one has to spend more effort on the basic routines computing the distance between two patches, but can profit from higher accuracy and a smaller number of faces. In this paper we present an efficient algorithm for computing the distance between so-called ‘quadratic complexes’. These are generalized polyhedra where the boundary consists of quadratic surface patches trimmed by quadratic curves. This class of surfaces plays an important part in designing mechanical parts and represents a set of modeling primitives provided by common CAD-systems. We reduce the distance calculation problem to the prob-
A
Central Surfaces: det( ) 6= 0 Ellipsoids / = 0 a0 6= 0 Hyperboloids Cone = 0 a0 = 0 Non-Central Surfaces: det( ) = 0 Paraboloids A3 = 0 a3 6= 0 a0 = 0 Ellipt. /Hyperbol. A3 = 0 = 0 a0 = 6 0 Cylinder Parabol. Cylinder A1 = A3 = 0 a1 6= 0 a0 = 0
a a
10 10 10
5 5
5
–5
A
x1
10
5
–5
–10
10
–5
5
–5
5
–10 –5
–5 5
10
x2 x1
5 5
–5
–5
–10
–10
x2 x1
5 10
10
10
10
–5
–10
5
5
–10
10
10
x1
x2
a
x2
x2 x1
10
–5
–10
–10
–5 5
–5
–5
–10
–10
–10
–10
–10
Sphere
Ellipsoid
Cone
One-SheetHyperboloid
Two-Sheet Hyperboloid
10
5
Table 1. Normal Forms of Quadratic Surfaces.
–10
–5
–5 5
x2 x1
–10
10
–10
–10
–5
5
–5
10
5
5 x2
10
10
10
–5
5
5
–5
–10
x2
10
–5
–5
5
–10
5 –5
–10 10
10 x1
–10 x1
–10
Elliptic Paraboloid
lem of solving univariate polynomial equations. For general quadrics we show that the degrees of the polynomials are bounded by 24. We have, however, identified an important subclass, the so-called natural quadrics, for which distance computation queries only require solving univariate polynomials of a degree of at most 8. We will also see that this upper bound is strict.
Parabolic Cylinder
Hyperbolic Cylinder
sphere, the cylinder and the circular cone. Analogously the line and the circle form the set of ‘natural conics’. Quadratic complexes with faces that are embedded on natural quadrics and trimmed by natural conics are consequently called ‘natural quadratic complexes’. Due to its importance in CAD it seems reasonable to extend the class of quadrics by the torus, which can be implicitly described by a polynomial of degree 4. However, the set of quadratic complexes is not closed under Boolean operations, since the intersection curves cannot always be represented by quadratic curves.
Rigid bodies are usually modeled as polyhedra and are represented by topological and geometrical descriptions of their boundaries. Thereby faces are embedded on planes and edges are given by line segments. Quadratic complexes can be considered as the most natural generalization of this polyhedral model with respect to a curved boundary representation. The faces of a quadratic complex are geometrically described by quadratic surfaces (quadrics) with edges represented by quadratic curves (conics). A quadratic surface is implicitly given by an algebraic equation of degree 2:
3. The Distance Computation Problem To get a precise mathematical formulation of the problem, we consider a quadratic complex as a non-empty, compact and connected set of points. The Euclidean distance δ between two points p and q can be extended to arbitrary point sets P and Q in a natural way:
(1)
for a vector a ∈ R3 and a symmetric matrix A ∈ R3×3 . Table 1 gives an overview over the so-called ‘normal forms’ (except from some degenerate cases), providing a classification that we will use in this paper. Typical examples of quadratic surfaces are shown in figure 1. A quadratic curve is explicitly given as the following point set: {p ∈ R3 | p = c + r(t)u + s(t)v},
Cylinder
Figure 1. Typical examples of Quadratic Surfaces.
2. Quadratic Complexes and Extensions
{x ∈ R3 | xT Ax + 2aT x + a0 = 0},
Hyperbolic Paraboloid
δ(P, Q) := inf{δ(p, q)|p ∈ P, q ∈ Q}.
(3)
From a practical point of view it makes sense to consider the quadratic distance function δ 2 , since in many cases square roots can be eliminated by this monotonous transformation. Apart from the (quadratic) distance value we request a pair of closest points as witness of the distance minimum between both quadratic complexes.
(2)
where (r, s) ∈ {(sin, cos), (sinh, cosh), (id, 0), (id, id2 )} and u, v ∈ R3 with uT v = 0. Hence, conics are ellipses, hyperbolas, lines or parabolas. We will also consider a subset of quadrics, the so-called ‘natural quadrics’. This class consists of the plane, the
Definition 1. Given two quadratic complexes C 1 , C 2 at a given point and orientation, the distance computation prob2
lem is to determine the global minimum of the distance function δ between the respective point sets, i.e.
Q2
(i) the value δ ∗ := δ(C 1 , C 2 ), (ii) a pair of points (p, q), s.t. δ ∗ = δ(p, q).
Q1
To compute the distance minimum, we can restrict ourselves to the boundaries that are described by the sets of faces. The next section consequently deals with closest points between quadratic surfaces.
Q1
f1 ∩ f2 = ∅:
case (i).
f1 ∩ f2 = ∅:
case (ii).
Q2 Q1
4. Closest Points Between Faces The basic problem is to determine the distance between two faces, that are – from a geometrical point of view – trimmed patches of quadratic surfaces Q1 and Q2 : Q1 Q2
Q2
f1 ∩ f2 6= ∅: Precondition violated.
Figure 2. Illustration of Theorem 1.
:= {x | xT Ax + 2aT x + a0 = 0}, := {y | y T By + 2bT y + b0 = 0}.
the problem min(x − y)2 , x ∈ Q1 y ∈ Q2 :
Let us first assume that the considered face pair is disjoint. If it is not, the distance value δ is zero and we are done. As long as we consider quadratic surfaces this condition can be easily checked by solving polynomials of a degree of at most 4 [7, 10]. Extending this class of surfaces by the torus, the former result still holds with the only exception that the intersection problem between a torus and a nonnatural quadratic surface requires the solution of a higher degree polynomial [6, 10]. The disjointness assumption allows us to characterize a pair of closest points between two faces:
L(x, y; α, β)
= (x − y)2 + α(xT Ax + 2aT x + a0 ) + β(y T By + 2bT y + b0 )
with the following L AGRANGE-conditions: ∂L(.)/∂x = 0 ⇔ α(Ax + a) = y − x
(4)
∂L(.)/∂y = 0 ⇔ β(By + b) = x − y ∂L(.)/∂α = 0 ⇔ xT Ax + 2aT x + a0 = 0
(5) (6)
∂L(.)/∂β = 0 ⇔ y T By + 2bT y + b0 = 0
(7)
where n(x) k Ax + a and n(y) k By + b. Figure 2 provides a 2-dimensional illustration of both cases and the necessary precondition. Closest points are marked by enclosing squares whereas the points of case (i), forming a distance extremum between Q1 and Q2 , are connected by their interpolating line segment.
Theorem 1. Let f1 and f2 be disjoint faces of quadratic complexes that are embedded on the quadratic surfaces Q1 and Q2 . If (p1 , p2 ) is a pair of closest points between f1 and f2 , then either (i) (p1 , p2 ) is an extremum of the distance function between Q1 and Q2 , i.e. there are λ, µ ∈ R, λ, µ 6= 0, s.t.
4.1. A Generic Algorithm n(p1 ) = λ(p2 − p1 )
n(p2 ) = µ(p1 − p2 ), Theorem 1 suggests a simple recursive procedure (see algorithm 3) for computing the distance between two faces. Arguments of the algorithm are the entities of the boundary representation, i.e. a face, an edge or a vertex. When called on two faces, the algorithm first checks the disjointness condition by the routine I NTERSECT. If the faces collide the distance value is set to zero and a pair of points witnessing the intersection is returned. Otherwise, all point pairs representing extrema of the distance function between both surfaces are computed. To do so, the subroutine E XTREMA enumerates all pairs (p1 , p2 ) that fulfill the conditions of case (i) for the untrimmed entities. In the variable δG we
where n(pi ) denotes the normal of Qi in pi , or (ii) p1 or p2 lies on the boundary of the face f1 or f2 , respectively. Proof. If neither p1 ∈ f1 nor p2 ∈ f2 are located at the boundary of the respective face, then both points form a local extremum of the distance function between Q1 and Q2 . Since both faces do not intersect, p1 and p2 are points for which the interpolating line is perpendicular to both surfaces in the respective points. This observation can be formally proven by setting up the L AGRANGE-function L for 3
Input: Entities E1 and E2 of type face, edge or vertex. Output: δ(E1 , E2 ) and a pair of closest points (p1 , p2 ). E NTITY D ISTANCE(E1 , E2 ) (1) isDisjoint, (p1 , p2 ) ← I NTERSECT(E1 , E2 ) (2) if isDisjoint = f alse (3) return 0, (p1 , p2 ) (4) δG ← ∞ (5) while δ, (q 1 , q 2 ) ← E XTREMA(E1 , E2 ) (6) if (q 1 ∈ E1 )and (q 2 ∈ E2 ) (7) if δ < δG (8) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (9) foreach subentity E of E1 (10) δ, (q1 , q 2 ) ← E NTITY D ISTANCE(E, E2 ) (11) if δ < δG (12) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (13) foreach subentity E of E2 (14) δ, (q1 , q 2 ) ← E NTITY D ISTANCE(E1 , E) (15) if δ < δG (16) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (17) return δG , (p1 , p2 )
4.2. Degree Complexity of the Polynomial Systems Procedure E XTREMA solves the central task of algorithm 3. It identifies the local extrema of the distance function between untrimmed entities of the boundary representation. In the next sections we will prove the following result which gives an upper bound on the degree complexity of the systems that have to be solved. Theorem 2. The distance between two faces of quadratic complexes can be computed by solving systems of univariate and bivariate polynomials in which the degree of every variable is 6 at most. These systems can be solved by finding the roots of univariate polynomials of a degree of at most 24.
4.3. The Surface-Surface Case In this section we will set up a system of bivariate polynomial equations that accomplishes this task in the case of two quadratic surfaces.
Figure 3. A generic algorithm for computing the distance between two entities of the boundary description.
Solving The Lagrange System In the proof of theorem 1 we have seen that the perpendicularity conditions of case (i) are a result of setting up the L AGRANGE formalism for the distance minimization problem between quadrics. By setting λ := 1/α and µ := 1/β we can derive from conditions (4) and (5):
compute the minimal distance of all these ’extremal’ points that lie inside the trimmed patches. However, the theorem also states that the closest points can also be found at the boundary of at least one of the patches (case(ii)). Therefore two recursive calls to the distance procedure must be executed, with the entity E1 in the first call and the entity E2 in the second call restricted to their boundaries ∂E1 and ∂E2 , respectively. In our algorithm the elements of ∂Ei are called ‘subentities’ of Ei , i = 1, 2. If the original arguments of the procedure were two faces, then the recursive calls would reciprocally determine the distance between the edges of one face and the other patch. The recursion ends if case (i) holds for a pair of points computed by procedure E XTREMA. This happens at the latest if the procedure is called on two vertices.
C λ,µ x
=
−(Ba + λb + µa),
(8)
C Tλ,µ y
=
−(Ab + λb + µa),
(9)
with C λ,µ := BA + λB + µA. If we solve these equations for x and y, we get: T
x = −C −1 λ,µ cB
y = −C −1 λ,µ cA ,
with cB := Ba + λb + µa and cA := Ab + λb + µa. Substituting the expressions for x and y in (6) and (7) gives the following system of polynomial equations: f (λ, µ) =
Algorithm 3 represents an elegant high-level formulation of the whole procedure. A practical implementation, however, must prevent repeated distance computation between identical pairs of entities. Moreover, one can avoid the recursive calls on the boundary elements when both the trimmed as well as the untrimmed entities do not intersect and the extremal points with minimal distance between the un-trimmend entities lie inside the trimmed patches.
T
cTB C λ,µ AC λ,µ cB − T
(10) 2
2|C λ,µ |a C λ,µ cB + a0 |C λ,µ | = 0, g(λ, µ) =
T
cTA C λ,µ BC λ,µ cA − T 2|C λ,µ |bT C λ,µ cA
(11)
+ b0 |C λ,µ |2 = 0,
where C λ,µ denotes the adjoint and |C λ,µ | the determinant of C λ,µ . 4
The Inverse of C λ,µ
These roots can be found with the help of elimination methods. Eliminating one variable and solving the univariate resultant polynomial Res(f, g) give the bivariate solutions projected onto the non-eliminated variable. Since the leading coefficients in λ as well as in µ are constants, the degree of Res(f, g) is exactly the number of bivariate roots and therefore bounded by 36. However, it is not necessary to solve this polynomial directly. The following lemma will show that, in our setting, the resultant is always the product of two polynomials of a lower degree. Solving the latter polynomials improves numerical stability as well as computational efficiency.
In order to solve the given system, we must express the inverse of C λ,µ as a bivariate matrix polynomial in λ and µ. The following proposition shows how the adjoint as well as the determinant can be written in terms of matrix coefficients in A and B. Proposition 1. The adjoint and determinant of C λ,µ = BA + λB + µA is given by C λ,µ
= Bλ2 + Aµ2 + T A Bλ + AT B µ + (T B T A − T AB )λµ + A B,
Lemma 1. Let f, g be the polynomials, given in (10) and (11), and the system h defined as follows:
|C λ,µ | = |B|λ3 + |A|µ3 + |A||B| + |B|tr(A)λ2 + |A|tr(B)µ2 +
h(λ, µ) := (h1 , h2 , h3 )T = C λ,µ cB = 0.
|B|tr(A)λ + |A|tr(B)µ + tr(BA)λ2 µ + tr(AB)λµ2 +
Then the roots of h = 0 are also solutions to f = g = 0.
(tr(A)tr(B) − tr(A B))λµ,
Proof. W.l.o.g. we can assume that the quadric given in normal form does not contain the origin, i.e. x 6= 0. Now let us consider values of λ and µ for which h vanishes. By multiplying equation (8) with C λ,µ from the lefthand side, we get |C λ,µ |x = −C λ,µ cB = 0 and therefore |C λ,µ | = 0. Applying the same transformation to (9) gives
where T M := tr(M )E − M for a matrix M ∈ R3×3 . Proof. We start with a simple problem, asking for the inverse of Gλ := D + λC. Since the characteristic polynomial of a matrix M ∈ R3×3 is given by 3
T
|C λ,µ |y = −C λ,µ cA = 0. Since the determinant vanishes,
2
|M − λE| = −λ + tr(M )λ − tr(M )λ + |M |.
T
the same must hold for C λ,µ cA . Now it is easy to see that every term in f and g becomes zero and we can conclude that the solutions of h solve the equations f and g.
Setting M := −DC −1 we get |Gλ | = −|C||M − λE| =
The system h = 0 consists of three polynomial equations hi = 0, i = 1, 2, 3, that have degree 2 in λ as well as µ. For every hi , hj with 1 ≤ i < j ≤ 3, it holds that the resultant rij = Res(hi , hj ) is a polynomial of degree 3 in the non-eliminated variable. The three polynomials rij are identical up to a constant factor and by Lemma 1 we know that every root of rij solves the resultant polynomial of the system f = g = 0. Moreover, one can show that the multiplicity of these roots is 4. Hence, we have found 4 a polynomial of degree 12, i.e. rij , that divides Res(f, g). These considerations prove that the closest points between two central surfaces can be determined by solving polynomials of a degree of at most 24. If Q1 is a cone and Q2 a general central surface then f simplifies to a polynomial of degree 4 in λ as well as µ. Hence, the resultant polynomial is of a degree of at most 24 and can be written as the product of two polynomials of degree 12. In the case of two cones we analogously observe that the roots of the resultant polynomial can be determined by solving polynomials of a maximum degree 4.
|C|λ3 + tr(DC)λ2 + tr(CD)λ + |D|. (12)
Now one can easily verify that the adjoint can be written as the following matrix polynomial: Gλ = Cλ2 + (T D T C − T CD )λ + D.
(13)
The determinant and the adjoint of C λ,µ are the result of setting D := BA + λB, C := A and recursively applying rule (12) and (13). The theorem shows that the adjoint of C λ,µ is a bivariate matrix polynomial of degree 2 in λ and µ. Since the degree of the determinant polynomial is 3, we can conclude that both polynomials have degree 6 in λ and µ. Central Surfaces For central surfaces the vector a vanishes in normal form. However, the degree of f and g does not decrease in this special case since it is still determined by |C λ,µ |2 . To analyze the complexity of solving this bivariate system, one first observes that both polynomials are sparse, having a total degree of 6. Hence, it follows from the B EZOUT bound that over C 2 the number of complex roots is at most 36.
Non-Central Surfaces In the case of non-central surfaces we have |A| = |B| = 0 and therefore |C λ,µ | simplifies to a quadratic polynomial 5
in both, λ and µ. This observation should give an intuition why the polynomial systems for non-central surfaces are less complex than in the case of central surfaces. If Q1 and Q2 are both paraboloids, then f and g become bivariate polynomials of degree 5 in λ and µ. Since the total degree is also 5 in both polynomials, the resultant is a degree 25 polynomial in the worst case. However, our considerations above have shown that Res(f, g) can be factorized into two lower degree polynomials. Since the equations hi = 0, i = 1, 2, 3, do not change in degree compared to the general case, one of the two factors is a polynomial of a degree of at most 12 and the other of a degree of at most 13. For elliptic and hyperbolic cylinders f and g collapse to degree 2 polynomials in λ and µ. The resultant is therefore a quartic polynomial that cannot be factorized as in previous cases. To see this, one should observe that in the special case of cylinders h = 0 iff λ or µ = 0. However, this is a degenerate situation since |C λ,µ | vanishes and therefore C −1 λ,µ does not exist. At this point we will omit the cases if Q1 is a central surface and Q2 not. As one would expect, the degrees of the polynomials are lower than in the case of central surfaces, but higher compared to the situation when only noncentral surfaces are involved. The results of this section are recorded in the following table:
a0 6= 0 a0 = 0 6= 0 =0
a a
Central Surfaces a0 6= 0 a0 = 0 24 12 4
Substituting x in (15) we get the univariate system: f (α)
T
with D α := E + αA. From equations (12) and (13), it follows that the adjoint of D α is a quadratic polynomial in α, whereas its determinant is of degree 3: Dα
=
Aα2 + T A α + E diag[d2 d3 , d1 d3 , d1 d2 ], |A|α3 + tr(A)α2 + tr(A)α + 1 d1 d2 d3 ,
where di := 1 + αAi and Ai denotes the ith diagonal element of A. Expanding both quantities in (17) gives a degree 6 polynomial of the form f (α) = A1 p2α1 d22 d23 + A2 p2α2 d21 d23 + A3 p2α3 d21 d22 + 2(a1 pα1 d1 d22 d23 a0 d21 d22 d23
+
a2 pα2 d21 d2 d23
+
(18)
a3 pα3 d21 d22 d3 )+
with pα := p − αa.
A
Central Surfaces If we consider central faces, the vector a is equal to 0 and therefore f simplifies to f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 + a0 d21 d22 d23 . The degree of this polynomial is still determined by the term d21 d22 d23 and therefore does not decrease compared with the general case. However, if we additionally assume that a0 vanishes, i.e. if we are dealing with a cone, then f becomes a polynomial of degree 4 in α:
As in the section before we will start by setting up the L AGRANGE-formalism for the respective distance minimization problem. W.l.o.g we assume that the surface is given in normal form. Let p denote the (transformed) query point. Then the L AGRANGE-function L is given by
f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 .
(19)
Non-Central Surfaces
L(x; α) = (x − p)2 + α(xT Ax + 2aT x + a0 ).
In the case of non-central surfaces we first analyze the situation if a = 0, i.e. if the given face is embedded on an elliptic or hyperbolic cylinder. Since A3 is equal to zero, d3 becomes 1 and in comparison to (18) the degree of the polynomial f decreases to 4:
The partial derivatives with respect to x and α yield the L AGRANGE-conditions: ⇐⇒ α(Ax + a) = p − x,
=
= |D α | =
4.4. The Point-Surface Case
∂L(.) =0 ∂x ∂L(.) =0 ∂α
(17) 2
2a Dα (p − αa)|D α | + a0 |D α | = 0,
Non-Central Surfaces rg = 1 18 12 8 8 4 2 13 8 5 4 2
a 6= 0 a = 0
= (p − αa)T D α ADα (p − αa) +
(14)
⇐⇒ xT Ax + 2aT x + a0 = 0. (15)
f (α) = A1 p21 d22 + A2 p22 d21 + a0 d21 d22 .
Again the first equation reflects the geometric intuition that the interpolating line through p and x is parallel to the surface normal in x. It can be transformed into
If the cylinder is parabolic, i.e. A1= A3= a0= 0, a16= 0, then d1 as well as d3 are constants and we obtain the following polynomial of degree 3 in α:
x = (E + αA)−1 (p − αa).
f (α) = a1 pα1 d22 .
(16) 6
(20)
Finally we have to consider the case of paraboloids for which a3 6= 0 and A3 as well as a0 vanish. Again d3 is equal to 1, reducing the degree of f to 5. f (α) = A1 p21 d22 + A2 p22 d21 + 2a3 pα3 d21 d22 .
we use the following substitutions to rationaly parameterize ellipses and hyperbolas: 1 − t2 2t , s(t) = , (Ellipse) (22) 2 1+t 1 + t2 1 + t2 2t r(t) = , s(t) = . (Hyperbola) (23) 1 − t2 1 − t2 To eliminate the denominators in pi , i = 1, 2, 3, we multiply f and g by (1 ± t2 ) and obtain a system of polynomial equations with degree 4 in t. Since the total degree of f and g is 10 and 6 respectively, the B EZOUT bound for the degree of the resultant is 60. In contrast to the surfacesurface case, this value drastically overestimates the actual degree. Therefore we computed the so-called ‘mixedvolume’ (M V (f, g)) [1], which – in the case of sparse polynomials – leads to a better bound on the degree of the resultant. For the given system f = g = 0 the mixed-volume bound is tight indicating a degree of 32. However, it is easy to observe that if the polynomial equation di = 0 has a solution αi , then the vector (αi , ti ) is a solution of the bivariate system for every ti solving the equation pi = 0, i = 1, 2, 3. Moreover, one can show that every αi is a root of multiplicity 4 in Res(f, g, t), whereas every ti has multiplicity 2 in Res(f, g, α). Hence Res(f, g) can be written as the following product: r(t) =
(21)
The following table summarizes our resultsin the pointsurface-case: Point - Central Surface Ellipsoid Hyperboloid Cone 6 6 4 Point - Non-Central Surface Paraboloids Elliptic / Hyperbolic Parabolic Cylinders Cylinder 5 4 3
4.5. The Curve-Surface Case In the former approach, the point p was considered a constant, not a variable. To derive the system of polynomial equations in the curve-surface case, we simply substitute p by the explicit representation of quadratic curves. P :
p(t) = c + r(t)u + s(t)v.
The new variable t introduces a third L AGRANGE-condition ∂L(.) = 0 ∂t
⇐⇒
(x − p)T
Res(f, g, t) = hα
∂p . ∂t
By substituting x according to (16), we get the condition: g(α, t) := D α (p − αa) − |D α |p
∂p ∂t
Res(f, g, α)
i=1 3 Y
d4i = hα
p2i = ht
i=1
= 0.
3 Y
(α − αi Ai )4 ,
i=1 3 Y
(t − ti1 )2 (t − ti2 )2 ,
i=1
where hα and ht are univariate polynomials of a degree of at most 20. If P is a parabola, then g becomes a cubic polynomial in t. Since r(t) = t and s(t) = t2 , pi is quadratic and p′i is linear in t, leading to a degree 4 polynomial f and a degree 3 polynomial g. The mixed-volume of f and g decreases to 26 and therefore the degrees of the polynomials hα , ht to 14. For the curve being a line, i.e. r(t) = t, s(t) = 0, the polynomial f is quadratic in t and g turns out to be a linear function. The resultant polynomial also simplifies such that the mixed-volume becomes 10. In this simple setting, however, the roots αi , i = 1, 2, 3 have multiplicity 2 and therefore hα and ht are polynomials of a degree of at most 4. In the special case of a cone surface, f also becomes a degree 4 polynomial in α, whereas the degree in t does not change compared to the setting discussed above:
Expanding the polynomials Dα and |Dα | yields the second equation of the system: g(α, t) = p˜1 p′1 d2 d3 + p˜2 p′2 d1 d3 + p˜3 p′3 d1 d2 = 0, ˜ := Ap + a and p′ := with p
= ht
3 Y
∂p ∂t .
Central Surfaces Let us first consider central surfaces for which a = 0. By interpreting (18) as a polynomial in α and t we get the following system of bivariate polynomial equations: f (α, t)=A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 +
f (α, t)=A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 = 0
a0 d21 d22 d23 = 0,
g(α, t)=A1 p1 p′1 d2 d3
g(α, t)=A1 p1 p′1 d2 d3 + A2 p2 p′2 d1 d3 + A3 p3 p′3 d1 d2 = 0.
+
A2 p2 p′2 d1 d3
+
A3 p3 p′3 d1 d2
(24) = 0.
Similar to the previous case, one can show that computing the roots of the resultant requires solving univariate polynomials of a degree of at most 12 (P: ellipse / hyperbola), 8 (P:
Whereas the degree in α is 6 for f and 2 for g, the degree of t obviously depends on the curve involved. First 7
parabola) and 2 (P: line) respectively. The following table gives an overview of the results described in this section: Curve - Central Surface Ellipsoid Hyperboloids Ellipse 20 20 Hyperbola 20 20 Parabola 14 14 Line 4 4
Our considerations in this section are summarized in the following table:
Cone 12 12 8 2
Ellipse Hyperbola Parabola Line
Non-Central Surfaces
To compute the distance between a point p and a quadratic curve Q, we assume w.l.o.g. that Q is embedded on the x1 x2 -plane and centered around the origin, i.e.
f (α, t) = A1 p21 d22 + A2 p22 d21 + a0 d21 d22 = 0, g(α, t) = A1 p1 p′1 d2 + A2 p2 p′2 d1 = 0.
Q:
Since d3 = 1, we obtain the following factorization of the resultant polynomial: Res(f, g, t) =
hα
2 Y
d4i = hα
i=1
Res(f, g, α) =
ht
2 Y
i=1
(α − αi Ai )4 ,
t
Setting the derivative of the distance function equal to zero gives
i=1
f (t) = rr′ u2 + ss′ v 2 − r′ pT u − s′ pT v = 0,
A similar result can be derived for the point-hyperbola problem. If Q is a parabola, then the degree decreases and f becomes a cubic polynomial in t:
g(α, t) = a1 p′1 d2 + A2 p2 p′2 = 0. To compute the roots of the resultant polynomial one has to solve polynomials of a degree of at most 8 (P: ellipse / hyperbola), 5 (P: parabola) and 1 (P: line) respectively. If the surface is an elliptic or hyperbolic paraboloid, we have already seen (cf. 21) that f is polynomial of degree 5 in α. It is easy to see that g also simplifies to a quadratic polynomial in α:
f (t) = 2v 2 t3 + (u2 − pT v)t − pT u = 0.
+
+
(26)
In the case of a line, condition (25) simplifies to the following linear equation: f (t) = u2 t − pT v = 0.
f (α, t) = A1 p21 d22 d+ A2 p22 d21 + 2a3 pα3 d21 d22 = 0, g(α, t) =
(25)
with r′ ≡ dd rt and s′ ≡ dd st . In the case of an ellipse, condition (25) is equivalent to the following polynomial equation of degree 4 in t: pT vt4 + 2 pT u − (v 2 − u2 ) t3 f (t) = −pT v + 2 pT u + (v 2 − u2 ) t2 = 0.
f (α, t) = a1 pα1 d22 = 0,
a3 p′3 d1 d2
uT v = 0.
min(p − r(t)u − s(t)v)2 .
i=1
2 Y p2i = ht (t − ti1 )2 (t − ti2 )2 .
A2 p2 p′2 d1
q(t) = r(t)u + s(t)v,
By considering the projection p of p onto the x1 -x2 -plane, the task is reduced to a 2-dimensional problem:
Moreover, we can see from the mixed-volume function that the degree of the resultant polynomial is at most 20 (P: ellipse / hyperbola), 16 (P: parabola) and 8 (P: line) respectively. Hence, the maximum degree of the univariate polynomials to solve is 12, 8 or 2, depending on the type of the curve. In the case of a parabolic cylinder a no longer vanishes, but A1 , A3 as well as a0 are all zero. Therefore f becomes a cubic and g a linear polynomial in α:
A1 p1 p′1 d2
Parab. Cylinder 8 8 5 1
4.6. The Point-Curve Case
For non-central sufaces with a = 0 and a0 6= 0, f remains a bivariate polynomial of degree 6 in α, whereas g simplifies to a linear function in α:
2 Y
Curve - Non-Central Surface Paraboloids Ellipt. / Hyperb. Cylinders 16 12 16 12 11 8 3 2
(27)
The following table shows the maximum degrees of the univariate polynomials to solve in the Point-Curve case:
= 0.
Again our factorization shows that the degrees of the univariate polynomials to solve are much lower than the degrees of the respective resultant polynomial.
Ellipse 4
8
Point - Curve Hyperbola Parabola 4 3
Line 1
4.7. The Curve-Curve Case
we have that Res(f, g, t2 ) also vanishes for these values and the quadratic polynomial (1 ± t21 ) must be a factor. Analogously one can show that in the ellipse-parabolaas well as in the hyperbola-parabola-case the two highest coefficients of f and g vanish for t1 = ±i or t1 = ±1, respectively. Therefore the degree 16 polynomial Res(f, g, t2 ) can be represented as the product of (1 ± t21 )2 and a univariate polynomial of degree 12 in t1 . The following table gives an overview of the degree bounds in this section:
The remaining case leads to the problem of computing the distance between two quadratic curves. To simplify matters, we still assume that one curve is located as described in 4.6, i.e. P :
p(t) = r1 (t1 )u1 + s1 (t1 )v 1 ,
uT1 v 1 = 0,
whereas the other curve is in arbitrary position and orientation:
Curve-Curve
Q:
q(t) = c2 + r2 (t2 )u2 + s2 (t2 )v 2 ,
uT2 v 2
= 0.
Ellipse / Hyperbola Parabola Line
The partial derivatives of the distance function δ 2 (t1 , t2 ) = (q(t2 ) − p(t1 ))2 yield the following system of bivariate equations: ∂ r1 ∂ s1 T [q(t2 ) − p(t1 )] − (28) u1 − v1 = 0, ∂ t1 ∂ t1 ∂ s2 ∂ r2 (29) u2 + v2 = 0. [q(t2 ) − p(t1 )]T ∂ t2 ∂ t2
g(t1 , t2 ) = (1 +
t21 )g1 (t2 )
+ (1 +
t22 )g2 (t1 , t2 ),
Parabola
Line
12 9
4 3 1
4.8. Natural Quadrics, Conics and the Torus Our statements so far hold for arbitrary quadratic curves and surfaces. This class, however, contains very simple geometric objects for which the general results can be substantially improved. With the natural quadrics and conics we have identified a set of special cases, which additionally is of practical importance.
If P and Q are both ellipses, then conditions (28) and (29) can be written in the following manner: f (t1 , t2 ) = (1 + t21 )f1 (t1 , t2 ) + (1 + t22 )f2 (t1 ),
Ellipse / Hyberbola 16
(30)
4.8.1 Natural Conics
(31)
Distance queries between points and ‘linear’ objects (lines, planes) are standard problems in computational geometry and computer graphics [2]. They involve simple linear systems that we will not discuss here. Hence, the circle remains the only interesting natural conic. Computing its distance to linear objects leads to polynomial equations of a degree of at most 4. In the pointcircle case we can see from section 4.6 that computing the closest points requires solving a quadratic polynomial. The circle-plane case shows the same degree complexity. We can rotate the coordinate system so that the plane is given by x3 = 0. Then the distance between a point of the circle p(φ) and the plane P is given as
where fi and gi , i = 1, 2, are polynomials of degrees at most 2 in t1 and t2 . This representation immediately shows that every (ξ1 , ξ2 ) ∈ {−i, i}2 solves the bivariate system. Moreover, it follows that (1 + t21 )2 is a factor of the resultant Res(f, g, t2 ). Since the mixed-volume of the sparse polynomials f and g is 20, we can conclude that computing the distance between two ellipses requires solving univariate polynomials of a degree of at most 16. The same result holds for the hyperbola-hyperbola- as well as for the ellipse-hyperbola-problem. In both cases one analogously observes that the degree 4 polynomial, (1 − t21 )2 and (1 + t21 )(1 − t21 )2 respectively, is a divisor of Res(f, g, t2 ). When considering two parabolas, one obtains polynomials f and g, each having total degree 3. Therefore Res(f, g) is a polynomial of degree 9 in the non-eliminated variable. The case of two lines or a line and a parabola is also straightforward, leading to a linear and cubic equation respectively. Setting up the bivariate system (28, 29) for P being an ellipse or hyperbola and Q a line yields to polynomials f and g that have a mixed-volume of 6. By considering f and g to be polynomials in t2 , it turns out that the leading coefficient vanishes for t1 = ±i (P : ellipse) or t1 = ±1 (P : hyperbola) in both polynomials. From the resultant theory
δ(p(φ), P ) = c3 + r(cos(φ)u3 + sin(φ)v3 ). Minimizing δ over φ and using the t-parameterization for the derivative yields to a quadratic equation in t. The circle-circle case is the most difficult one since it requires the solutions of a polynomial of degree 8 [2]. Moreover, N EFF has shown that this is a lower bound on the degree complexity if the arithmetical operations are restricted to field operations and square roots [4]. 4.8.2 Natural Quadrics and the Torus Concerning the natural quadrics we first observe that the cases involving spheres and cylinders are already solved 9
since they can be reduced to the case of lower dimensional objects. The distance to a sphere is just the distance to its center minus the radius. In the case of a circular cylinder we can analogously restrict ourselves to its axis and finally subtract the radius. The same argument holds for the torus – although it is not a quadratic surface. Here, the distance can be computed by considering its main circle. Finally, we have to treat the case of a circular cone. From (19) we know that the extrema of the distance function between a point and a general cone is given by:
have therefore considered algebraic techniques that reduce the multivariate root finding problem to the univariate case. The degrees of the resulting polynomials are usually high and hence elimination methods are often not appropriate in the context of floating-point computations. However, we could show that for all cases in which the degrees of the resulting polynomials become critical, a factorization into two lower degree polynomials can be found. Thereby one of the polynomial factors is easy to solve because its roots are of high multiplicity, whereas the other is of moderate degree and can be solved using numerically stable techniques. Since in most CAD-systems the intersection curve between two quadratic surfaces is represented by NURBS, we will extend the set of trimming curves by this more complex type. Then our object class can be considered to be (approximately) closed under Boolean operations. In a final step, we will speed up the overall algorithm by setting up a bounding volume hierarchy that reduces the number of distance computations between pairs of faces.
f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 = 0. (32) In the special case of a circular cone, we have A1 = A2 and d1 = d2 . Condition (32) therefore simplifies to a polynomial equation of degree 2 in α: f (α) = A1 p21 d23 + A1 p22 d23 + A3 p23 d21 = 0. For the curve-cone problem we consider system (24). Exploiting the fact that the cone is circular, we get: f (α, t) g(α, t)
= =
A1 p21 d23 + A1 p22 d23 + A3 p23 d21 = 0, A1 p1 p′1 d3 + A1 p2 p′2 d3 + A3 p3 p′3 d1
References [1] J. Canny and I. Emiris. An efficient algorithm for the sparse mixed resultant. In Proc. Intern. Symp.Applied Algebra, Algebraic Algor. and Error-Corr. Codes, volume 263 of Lecture Notes in Computer Science., pages 89–104. SpringerVerlag, 1993. [2] D. H. Eberly. 3D game engine design : a practical approach to real-time computer graphics. Morgan Kaufmann, 2001. [3] S. Ehmann and M. Lin. Accurate and fast proximity queries between polyhedra using surface decomposition. In Computer Graphics Forum (Proc. Eurographics), 2001. [4] R. Farouki, C. Neff, and M. O’Connor. Automatic parsing of degenerate quadric-surface intersection. ACM Transactions on Graphics, 8(3):174–203, 1989. [5] E. Gilbert, D. Johnson, and S. Keerthi. A fast procedure for computing the distance between complex objects in 3d space. IEEE Journ. of Robot. Autom., 4(2):193–203, 1988. [6] K.-J. Kim. Torus and Simple Surface Intersection Based on a Configuration Space Approach. Ph.D. thesis, Dep. of Computer Science and Engineering, POSTECH, Feb. 1998. [7] J. Levin. A parametric algorithm for drawing pictures of solid objects composed of quadric surfaces. Commun. ACM, 19(10):555–563, 1976. [8] M. C. Lin and J. F. Canny. Efficient algorithms for incremental distance computation. In Proc. IEEE Internat. Conf. Robot. Autom., volume 2, pages 1008–1014, 1991. [9] B. Mirtich. Impulse-based dynamic simulation of rigid body systems. PhD thesis, University of California, Berkeley, 1996. [10] J. Reichel, E. Sch¨omer, T. Warken, and C. Lennerz. Efficient collision detection for curved solid objects. In Proc. 7th ACM Symp. on Solid Modeling and Applications, SM’ 02, 2002. [11] C. Turnbull and S. Cameron. Computing distances between nurbs-defined convex objects. In IEEE Int. Conf. of Robot. Autom., pages 3686–3690, 1998.
= 0.
For p being a line, the resultant of f and g is a polynomial of degree 4 for which we know a divisor of degree 2 (cf. our considerations in 4.5). In the case of a circle, the mixedvolume of f and g is 12. Since we have already realized that d43 is a factor of Res(f, g, t), we can restrict ourselves to solving polynomials of a degree of at most 8. The remaining case is the one in which two cones are involved. For general cones we have seen in 4.3 that polynomials of a degree of at most 4 must be solved to compute the extrema of the distance function. In the special case of two circular cones this bound cannot be improved. Our considerations above prove the following theorem. Theorem 3. The distance between two faces of natural quadratic complexes can be computed by solving univariate polynomials of a degree of at most 8.
5. Conclusions In this paper we have presented an algorithm for computing the minimal distance between two patches of quadratic surfaces bounded by quadratic curves. To find the closest points on both patches, systems of univariate and bivariate polynomial equations have to be solved. Though there are numerical techniques for finding the solutions of multivariate systems directly, our experience has shown that these methods either do not guarantee to find all roots of the system (N EWTON method) or they are not efficient enough (interval methods) to fulfill our real-time requirements. We 10
Abstract Virtual prototyping and assembly planning require physically based simulation techniques. In this setting the relevant objects are mostly mechanical parts, designed in CADprograms. When exported to the prototyping and planning systems, curved parts are approximated by large polygonal models, thus confronting the simulation algorithms with high complexity. Algorithms for collision detection in particular are a bottleneck of efficiency and suffer from accuracy and robustness problems. To overcome these problems, our algorithm directly operates on the original CAD-data. This approach reduces the input complexity and avoids accuracy problems due to approximation errors. We present an efficient algorithm for computing the distance between patches of quadratic surfaces trimmed by quadratic curves. The distance calculation problem is reduced to the problem of solving univariate polynomials of a degree of at most 24. Moreover, we will identify an important subclass for which the degree of the polynomials is bounded by 8.
1. Introduction During the past few years an increasing interest in virtual reality (VR) techniques could be observed. Virtual reality plays an increasingly important part in social and business life, ranging from entertainment and education to VR-based training methods and numerous industrial applications. However, the ‘classical’ challenges that come from concepts like virtual prototyping and assembly planning still represent strong drivers for the development of VR-techniques. To evaluate industrial prototypes and assembly processes, virtual environments have to simulate physical behavior. Particularly collisions between moving objects must be detected and resolved. Robustness and efficiency of such a simulation heavily depend on the quality of the collision detection algorithm. Dynamic collision detection routines are required in this context. They identify the earliest time
of collision in a given time interval and return the touching pair of objects. In contrast to this technique, static detection algorithms only decide whether two objects intersect at a given point in time. Since they check discrete time frames they, however, cannot guarantee to find the earliest contact. In [9] a framework is developed in which static distance computation algorithms can be used to design dynamic detection routines. The basic idea is that information about the distance and dynamic states of all objects provides lower bounds on the earliest time of collision. The problem of distance computation has been well studied in the past for polyhedral objects, e.g. [5, 8, 3], but only little effort has spent on objects with curved surfaces [11]. There are two reasons for this situation. First, a polygonal representation is not considered a true restriction since real objects can be approximated arbitrarily precisely by a polyhedron. The second reason is the fact that the basic algorithms and predicates can be implemented robustly and very efficiently on polygons. However, one has to consider that there is a trade-off between the accuracy of the approximation and the efficiency of the distance computation. On the one hand, precise approximations increase the complexity of the boundary representation and can affect the running time of the algorithm disproportionately. On the other hand, a small number of polygons compromises the quality or even the correctness of the simulation. Our considerations therefore raise the question, whether or not it is worth to handle curved objects directly. It means that one has to spend more effort on the basic routines computing the distance between two patches, but can profit from higher accuracy and a smaller number of faces. In this paper we present an efficient algorithm for computing the distance between so-called ‘quadratic complexes’. These are generalized polyhedra where the boundary consists of quadratic surface patches trimmed by quadratic curves. This class of surfaces plays an important part in designing mechanical parts and represents a set of modeling primitives provided by common CAD-systems. We reduce the distance calculation problem to the prob-
A
Central Surfaces: det( ) 6= 0 Ellipsoids / = 0 a0 6= 0 Hyperboloids Cone = 0 a0 = 0 Non-Central Surfaces: det( ) = 0 Paraboloids A3 = 0 a3 6= 0 a0 = 0 Ellipt. /Hyperbol. A3 = 0 = 0 a0 = 6 0 Cylinder Parabol. Cylinder A1 = A3 = 0 a1 6= 0 a0 = 0
a a
10 10 10
5 5
5
–5
A
x1
10
5
–5
–10
10
–5
5
–5
5
–10 –5
–5 5
10
x2 x1
5 5
–5
–5
–10
–10
x2 x1
5 10
10
10
10
–5
–10
5
5
–10
10
10
x1
x2
a
x2
x2 x1
10
–5
–10
–10
–5 5
–5
–5
–10
–10
–10
–10
–10
Sphere
Ellipsoid
Cone
One-SheetHyperboloid
Two-Sheet Hyperboloid
10
5
Table 1. Normal Forms of Quadratic Surfaces.
–10
–5
–5 5
x2 x1
–10
10
–10
–10
–5
5
–5
10
5
5 x2
10
10
10
–5
5
5
–5
–10
x2
10
–5
–5
5
–10
5 –5
–10 10
10 x1
–10 x1
–10
Elliptic Paraboloid
lem of solving univariate polynomial equations. For general quadrics we show that the degrees of the polynomials are bounded by 24. We have, however, identified an important subclass, the so-called natural quadrics, for which distance computation queries only require solving univariate polynomials of a degree of at most 8. We will also see that this upper bound is strict.
Parabolic Cylinder
Hyperbolic Cylinder
sphere, the cylinder and the circular cone. Analogously the line and the circle form the set of ‘natural conics’. Quadratic complexes with faces that are embedded on natural quadrics and trimmed by natural conics are consequently called ‘natural quadratic complexes’. Due to its importance in CAD it seems reasonable to extend the class of quadrics by the torus, which can be implicitly described by a polynomial of degree 4. However, the set of quadratic complexes is not closed under Boolean operations, since the intersection curves cannot always be represented by quadratic curves.
Rigid bodies are usually modeled as polyhedra and are represented by topological and geometrical descriptions of their boundaries. Thereby faces are embedded on planes and edges are given by line segments. Quadratic complexes can be considered as the most natural generalization of this polyhedral model with respect to a curved boundary representation. The faces of a quadratic complex are geometrically described by quadratic surfaces (quadrics) with edges represented by quadratic curves (conics). A quadratic surface is implicitly given by an algebraic equation of degree 2:
3. The Distance Computation Problem To get a precise mathematical formulation of the problem, we consider a quadratic complex as a non-empty, compact and connected set of points. The Euclidean distance δ between two points p and q can be extended to arbitrary point sets P and Q in a natural way:
(1)
for a vector a ∈ R3 and a symmetric matrix A ∈ R3×3 . Table 1 gives an overview over the so-called ‘normal forms’ (except from some degenerate cases), providing a classification that we will use in this paper. Typical examples of quadratic surfaces are shown in figure 1. A quadratic curve is explicitly given as the following point set: {p ∈ R3 | p = c + r(t)u + s(t)v},
Cylinder
Figure 1. Typical examples of Quadratic Surfaces.
2. Quadratic Complexes and Extensions
{x ∈ R3 | xT Ax + 2aT x + a0 = 0},
Hyperbolic Paraboloid
δ(P, Q) := inf{δ(p, q)|p ∈ P, q ∈ Q}.
(3)
From a practical point of view it makes sense to consider the quadratic distance function δ 2 , since in many cases square roots can be eliminated by this monotonous transformation. Apart from the (quadratic) distance value we request a pair of closest points as witness of the distance minimum between both quadratic complexes.
(2)
where (r, s) ∈ {(sin, cos), (sinh, cosh), (id, 0), (id, id2 )} and u, v ∈ R3 with uT v = 0. Hence, conics are ellipses, hyperbolas, lines or parabolas. We will also consider a subset of quadrics, the so-called ‘natural quadrics’. This class consists of the plane, the
Definition 1. Given two quadratic complexes C 1 , C 2 at a given point and orientation, the distance computation prob2
lem is to determine the global minimum of the distance function δ between the respective point sets, i.e.
Q2
(i) the value δ ∗ := δ(C 1 , C 2 ), (ii) a pair of points (p, q), s.t. δ ∗ = δ(p, q).
Q1
To compute the distance minimum, we can restrict ourselves to the boundaries that are described by the sets of faces. The next section consequently deals with closest points between quadratic surfaces.
Q1
f1 ∩ f2 = ∅:
case (i).
f1 ∩ f2 = ∅:
case (ii).
Q2 Q1
4. Closest Points Between Faces The basic problem is to determine the distance between two faces, that are – from a geometrical point of view – trimmed patches of quadratic surfaces Q1 and Q2 : Q1 Q2
Q2
f1 ∩ f2 6= ∅: Precondition violated.
Figure 2. Illustration of Theorem 1.
:= {x | xT Ax + 2aT x + a0 = 0}, := {y | y T By + 2bT y + b0 = 0}.
the problem min(x − y)2 , x ∈ Q1 y ∈ Q2 :
Let us first assume that the considered face pair is disjoint. If it is not, the distance value δ is zero and we are done. As long as we consider quadratic surfaces this condition can be easily checked by solving polynomials of a degree of at most 4 [7, 10]. Extending this class of surfaces by the torus, the former result still holds with the only exception that the intersection problem between a torus and a nonnatural quadratic surface requires the solution of a higher degree polynomial [6, 10]. The disjointness assumption allows us to characterize a pair of closest points between two faces:
L(x, y; α, β)
= (x − y)2 + α(xT Ax + 2aT x + a0 ) + β(y T By + 2bT y + b0 )
with the following L AGRANGE-conditions: ∂L(.)/∂x = 0 ⇔ α(Ax + a) = y − x
(4)
∂L(.)/∂y = 0 ⇔ β(By + b) = x − y ∂L(.)/∂α = 0 ⇔ xT Ax + 2aT x + a0 = 0
(5) (6)
∂L(.)/∂β = 0 ⇔ y T By + 2bT y + b0 = 0
(7)
where n(x) k Ax + a and n(y) k By + b. Figure 2 provides a 2-dimensional illustration of both cases and the necessary precondition. Closest points are marked by enclosing squares whereas the points of case (i), forming a distance extremum between Q1 and Q2 , are connected by their interpolating line segment.
Theorem 1. Let f1 and f2 be disjoint faces of quadratic complexes that are embedded on the quadratic surfaces Q1 and Q2 . If (p1 , p2 ) is a pair of closest points between f1 and f2 , then either (i) (p1 , p2 ) is an extremum of the distance function between Q1 and Q2 , i.e. there are λ, µ ∈ R, λ, µ 6= 0, s.t.
4.1. A Generic Algorithm n(p1 ) = λ(p2 − p1 )
n(p2 ) = µ(p1 − p2 ), Theorem 1 suggests a simple recursive procedure (see algorithm 3) for computing the distance between two faces. Arguments of the algorithm are the entities of the boundary representation, i.e. a face, an edge or a vertex. When called on two faces, the algorithm first checks the disjointness condition by the routine I NTERSECT. If the faces collide the distance value is set to zero and a pair of points witnessing the intersection is returned. Otherwise, all point pairs representing extrema of the distance function between both surfaces are computed. To do so, the subroutine E XTREMA enumerates all pairs (p1 , p2 ) that fulfill the conditions of case (i) for the untrimmed entities. In the variable δG we
where n(pi ) denotes the normal of Qi in pi , or (ii) p1 or p2 lies on the boundary of the face f1 or f2 , respectively. Proof. If neither p1 ∈ f1 nor p2 ∈ f2 are located at the boundary of the respective face, then both points form a local extremum of the distance function between Q1 and Q2 . Since both faces do not intersect, p1 and p2 are points for which the interpolating line is perpendicular to both surfaces in the respective points. This observation can be formally proven by setting up the L AGRANGE-function L for 3
Input: Entities E1 and E2 of type face, edge or vertex. Output: δ(E1 , E2 ) and a pair of closest points (p1 , p2 ). E NTITY D ISTANCE(E1 , E2 ) (1) isDisjoint, (p1 , p2 ) ← I NTERSECT(E1 , E2 ) (2) if isDisjoint = f alse (3) return 0, (p1 , p2 ) (4) δG ← ∞ (5) while δ, (q 1 , q 2 ) ← E XTREMA(E1 , E2 ) (6) if (q 1 ∈ E1 )and (q 2 ∈ E2 ) (7) if δ < δG (8) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (9) foreach subentity E of E1 (10) δ, (q1 , q 2 ) ← E NTITY D ISTANCE(E, E2 ) (11) if δ < δG (12) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (13) foreach subentity E of E2 (14) δ, (q1 , q 2 ) ← E NTITY D ISTANCE(E1 , E) (15) if δ < δG (16) δG ← δ, (p1 , p2 ) ← (q 1 , q 2 ) (17) return δG , (p1 , p2 )
4.2. Degree Complexity of the Polynomial Systems Procedure E XTREMA solves the central task of algorithm 3. It identifies the local extrema of the distance function between untrimmed entities of the boundary representation. In the next sections we will prove the following result which gives an upper bound on the degree complexity of the systems that have to be solved. Theorem 2. The distance between two faces of quadratic complexes can be computed by solving systems of univariate and bivariate polynomials in which the degree of every variable is 6 at most. These systems can be solved by finding the roots of univariate polynomials of a degree of at most 24.
4.3. The Surface-Surface Case In this section we will set up a system of bivariate polynomial equations that accomplishes this task in the case of two quadratic surfaces.
Figure 3. A generic algorithm for computing the distance between two entities of the boundary description.
Solving The Lagrange System In the proof of theorem 1 we have seen that the perpendicularity conditions of case (i) are a result of setting up the L AGRANGE formalism for the distance minimization problem between quadrics. By setting λ := 1/α and µ := 1/β we can derive from conditions (4) and (5):
compute the minimal distance of all these ’extremal’ points that lie inside the trimmed patches. However, the theorem also states that the closest points can also be found at the boundary of at least one of the patches (case(ii)). Therefore two recursive calls to the distance procedure must be executed, with the entity E1 in the first call and the entity E2 in the second call restricted to their boundaries ∂E1 and ∂E2 , respectively. In our algorithm the elements of ∂Ei are called ‘subentities’ of Ei , i = 1, 2. If the original arguments of the procedure were two faces, then the recursive calls would reciprocally determine the distance between the edges of one face and the other patch. The recursion ends if case (i) holds for a pair of points computed by procedure E XTREMA. This happens at the latest if the procedure is called on two vertices.
C λ,µ x
=
−(Ba + λb + µa),
(8)
C Tλ,µ y
=
−(Ab + λb + µa),
(9)
with C λ,µ := BA + λB + µA. If we solve these equations for x and y, we get: T
x = −C −1 λ,µ cB
y = −C −1 λ,µ cA ,
with cB := Ba + λb + µa and cA := Ab + λb + µa. Substituting the expressions for x and y in (6) and (7) gives the following system of polynomial equations: f (λ, µ) =
Algorithm 3 represents an elegant high-level formulation of the whole procedure. A practical implementation, however, must prevent repeated distance computation between identical pairs of entities. Moreover, one can avoid the recursive calls on the boundary elements when both the trimmed as well as the untrimmed entities do not intersect and the extremal points with minimal distance between the un-trimmend entities lie inside the trimmed patches.
T
cTB C λ,µ AC λ,µ cB − T
(10) 2
2|C λ,µ |a C λ,µ cB + a0 |C λ,µ | = 0, g(λ, µ) =
T
cTA C λ,µ BC λ,µ cA − T 2|C λ,µ |bT C λ,µ cA
(11)
+ b0 |C λ,µ |2 = 0,
where C λ,µ denotes the adjoint and |C λ,µ | the determinant of C λ,µ . 4
The Inverse of C λ,µ
These roots can be found with the help of elimination methods. Eliminating one variable and solving the univariate resultant polynomial Res(f, g) give the bivariate solutions projected onto the non-eliminated variable. Since the leading coefficients in λ as well as in µ are constants, the degree of Res(f, g) is exactly the number of bivariate roots and therefore bounded by 36. However, it is not necessary to solve this polynomial directly. The following lemma will show that, in our setting, the resultant is always the product of two polynomials of a lower degree. Solving the latter polynomials improves numerical stability as well as computational efficiency.
In order to solve the given system, we must express the inverse of C λ,µ as a bivariate matrix polynomial in λ and µ. The following proposition shows how the adjoint as well as the determinant can be written in terms of matrix coefficients in A and B. Proposition 1. The adjoint and determinant of C λ,µ = BA + λB + µA is given by C λ,µ
= Bλ2 + Aµ2 + T A Bλ + AT B µ + (T B T A − T AB )λµ + A B,
Lemma 1. Let f, g be the polynomials, given in (10) and (11), and the system h defined as follows:
|C λ,µ | = |B|λ3 + |A|µ3 + |A||B| + |B|tr(A)λ2 + |A|tr(B)µ2 +
h(λ, µ) := (h1 , h2 , h3 )T = C λ,µ cB = 0.
|B|tr(A)λ + |A|tr(B)µ + tr(BA)λ2 µ + tr(AB)λµ2 +
Then the roots of h = 0 are also solutions to f = g = 0.
(tr(A)tr(B) − tr(A B))λµ,
Proof. W.l.o.g. we can assume that the quadric given in normal form does not contain the origin, i.e. x 6= 0. Now let us consider values of λ and µ for which h vanishes. By multiplying equation (8) with C λ,µ from the lefthand side, we get |C λ,µ |x = −C λ,µ cB = 0 and therefore |C λ,µ | = 0. Applying the same transformation to (9) gives
where T M := tr(M )E − M for a matrix M ∈ R3×3 . Proof. We start with a simple problem, asking for the inverse of Gλ := D + λC. Since the characteristic polynomial of a matrix M ∈ R3×3 is given by 3
T
|C λ,µ |y = −C λ,µ cA = 0. Since the determinant vanishes,
2
|M − λE| = −λ + tr(M )λ − tr(M )λ + |M |.
T
the same must hold for C λ,µ cA . Now it is easy to see that every term in f and g becomes zero and we can conclude that the solutions of h solve the equations f and g.
Setting M := −DC −1 we get |Gλ | = −|C||M − λE| =
The system h = 0 consists of three polynomial equations hi = 0, i = 1, 2, 3, that have degree 2 in λ as well as µ. For every hi , hj with 1 ≤ i < j ≤ 3, it holds that the resultant rij = Res(hi , hj ) is a polynomial of degree 3 in the non-eliminated variable. The three polynomials rij are identical up to a constant factor and by Lemma 1 we know that every root of rij solves the resultant polynomial of the system f = g = 0. Moreover, one can show that the multiplicity of these roots is 4. Hence, we have found 4 a polynomial of degree 12, i.e. rij , that divides Res(f, g). These considerations prove that the closest points between two central surfaces can be determined by solving polynomials of a degree of at most 24. If Q1 is a cone and Q2 a general central surface then f simplifies to a polynomial of degree 4 in λ as well as µ. Hence, the resultant polynomial is of a degree of at most 24 and can be written as the product of two polynomials of degree 12. In the case of two cones we analogously observe that the roots of the resultant polynomial can be determined by solving polynomials of a maximum degree 4.
|C|λ3 + tr(DC)λ2 + tr(CD)λ + |D|. (12)
Now one can easily verify that the adjoint can be written as the following matrix polynomial: Gλ = Cλ2 + (T D T C − T CD )λ + D.
(13)
The determinant and the adjoint of C λ,µ are the result of setting D := BA + λB, C := A and recursively applying rule (12) and (13). The theorem shows that the adjoint of C λ,µ is a bivariate matrix polynomial of degree 2 in λ and µ. Since the degree of the determinant polynomial is 3, we can conclude that both polynomials have degree 6 in λ and µ. Central Surfaces For central surfaces the vector a vanishes in normal form. However, the degree of f and g does not decrease in this special case since it is still determined by |C λ,µ |2 . To analyze the complexity of solving this bivariate system, one first observes that both polynomials are sparse, having a total degree of 6. Hence, it follows from the B EZOUT bound that over C 2 the number of complex roots is at most 36.
Non-Central Surfaces In the case of non-central surfaces we have |A| = |B| = 0 and therefore |C λ,µ | simplifies to a quadratic polynomial 5
in both, λ and µ. This observation should give an intuition why the polynomial systems for non-central surfaces are less complex than in the case of central surfaces. If Q1 and Q2 are both paraboloids, then f and g become bivariate polynomials of degree 5 in λ and µ. Since the total degree is also 5 in both polynomials, the resultant is a degree 25 polynomial in the worst case. However, our considerations above have shown that Res(f, g) can be factorized into two lower degree polynomials. Since the equations hi = 0, i = 1, 2, 3, do not change in degree compared to the general case, one of the two factors is a polynomial of a degree of at most 12 and the other of a degree of at most 13. For elliptic and hyperbolic cylinders f and g collapse to degree 2 polynomials in λ and µ. The resultant is therefore a quartic polynomial that cannot be factorized as in previous cases. To see this, one should observe that in the special case of cylinders h = 0 iff λ or µ = 0. However, this is a degenerate situation since |C λ,µ | vanishes and therefore C −1 λ,µ does not exist. At this point we will omit the cases if Q1 is a central surface and Q2 not. As one would expect, the degrees of the polynomials are lower than in the case of central surfaces, but higher compared to the situation when only noncentral surfaces are involved. The results of this section are recorded in the following table:
a0 6= 0 a0 = 0 6= 0 =0
a a
Central Surfaces a0 6= 0 a0 = 0 24 12 4
Substituting x in (15) we get the univariate system: f (α)
T
with D α := E + αA. From equations (12) and (13), it follows that the adjoint of D α is a quadratic polynomial in α, whereas its determinant is of degree 3: Dα
=
Aα2 + T A α + E diag[d2 d3 , d1 d3 , d1 d2 ], |A|α3 + tr(A)α2 + tr(A)α + 1 d1 d2 d3 ,
where di := 1 + αAi and Ai denotes the ith diagonal element of A. Expanding both quantities in (17) gives a degree 6 polynomial of the form f (α) = A1 p2α1 d22 d23 + A2 p2α2 d21 d23 + A3 p2α3 d21 d22 + 2(a1 pα1 d1 d22 d23 a0 d21 d22 d23
+
a2 pα2 d21 d2 d23
+
(18)
a3 pα3 d21 d22 d3 )+
with pα := p − αa.
A
Central Surfaces If we consider central faces, the vector a is equal to 0 and therefore f simplifies to f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 + a0 d21 d22 d23 . The degree of this polynomial is still determined by the term d21 d22 d23 and therefore does not decrease compared with the general case. However, if we additionally assume that a0 vanishes, i.e. if we are dealing with a cone, then f becomes a polynomial of degree 4 in α:
As in the section before we will start by setting up the L AGRANGE-formalism for the respective distance minimization problem. W.l.o.g we assume that the surface is given in normal form. Let p denote the (transformed) query point. Then the L AGRANGE-function L is given by
f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 .
(19)
Non-Central Surfaces
L(x; α) = (x − p)2 + α(xT Ax + 2aT x + a0 ).
In the case of non-central surfaces we first analyze the situation if a = 0, i.e. if the given face is embedded on an elliptic or hyperbolic cylinder. Since A3 is equal to zero, d3 becomes 1 and in comparison to (18) the degree of the polynomial f decreases to 4:
The partial derivatives with respect to x and α yield the L AGRANGE-conditions: ⇐⇒ α(Ax + a) = p − x,
=
= |D α | =
4.4. The Point-Surface Case
∂L(.) =0 ∂x ∂L(.) =0 ∂α
(17) 2
2a Dα (p − αa)|D α | + a0 |D α | = 0,
Non-Central Surfaces rg = 1 18 12 8 8 4 2 13 8 5 4 2
a 6= 0 a = 0
= (p − αa)T D α ADα (p − αa) +
(14)
⇐⇒ xT Ax + 2aT x + a0 = 0. (15)
f (α) = A1 p21 d22 + A2 p22 d21 + a0 d21 d22 .
Again the first equation reflects the geometric intuition that the interpolating line through p and x is parallel to the surface normal in x. It can be transformed into
If the cylinder is parabolic, i.e. A1= A3= a0= 0, a16= 0, then d1 as well as d3 are constants and we obtain the following polynomial of degree 3 in α:
x = (E + αA)−1 (p − αa).
f (α) = a1 pα1 d22 .
(16) 6
(20)
Finally we have to consider the case of paraboloids for which a3 6= 0 and A3 as well as a0 vanish. Again d3 is equal to 1, reducing the degree of f to 5. f (α) = A1 p21 d22 + A2 p22 d21 + 2a3 pα3 d21 d22 .
we use the following substitutions to rationaly parameterize ellipses and hyperbolas: 1 − t2 2t , s(t) = , (Ellipse) (22) 2 1+t 1 + t2 1 + t2 2t r(t) = , s(t) = . (Hyperbola) (23) 1 − t2 1 − t2 To eliminate the denominators in pi , i = 1, 2, 3, we multiply f and g by (1 ± t2 ) and obtain a system of polynomial equations with degree 4 in t. Since the total degree of f and g is 10 and 6 respectively, the B EZOUT bound for the degree of the resultant is 60. In contrast to the surfacesurface case, this value drastically overestimates the actual degree. Therefore we computed the so-called ‘mixedvolume’ (M V (f, g)) [1], which – in the case of sparse polynomials – leads to a better bound on the degree of the resultant. For the given system f = g = 0 the mixed-volume bound is tight indicating a degree of 32. However, it is easy to observe that if the polynomial equation di = 0 has a solution αi , then the vector (αi , ti ) is a solution of the bivariate system for every ti solving the equation pi = 0, i = 1, 2, 3. Moreover, one can show that every αi is a root of multiplicity 4 in Res(f, g, t), whereas every ti has multiplicity 2 in Res(f, g, α). Hence Res(f, g) can be written as the following product: r(t) =
(21)
The following table summarizes our resultsin the pointsurface-case: Point - Central Surface Ellipsoid Hyperboloid Cone 6 6 4 Point - Non-Central Surface Paraboloids Elliptic / Hyperbolic Parabolic Cylinders Cylinder 5 4 3
4.5. The Curve-Surface Case In the former approach, the point p was considered a constant, not a variable. To derive the system of polynomial equations in the curve-surface case, we simply substitute p by the explicit representation of quadratic curves. P :
p(t) = c + r(t)u + s(t)v.
The new variable t introduces a third L AGRANGE-condition ∂L(.) = 0 ∂t
⇐⇒
(x − p)T
Res(f, g, t) = hα
∂p . ∂t
By substituting x according to (16), we get the condition: g(α, t) := D α (p − αa) − |D α |p
∂p ∂t
Res(f, g, α)
i=1 3 Y
d4i = hα
p2i = ht
i=1
= 0.
3 Y
(α − αi Ai )4 ,
i=1 3 Y
(t − ti1 )2 (t − ti2 )2 ,
i=1
where hα and ht are univariate polynomials of a degree of at most 20. If P is a parabola, then g becomes a cubic polynomial in t. Since r(t) = t and s(t) = t2 , pi is quadratic and p′i is linear in t, leading to a degree 4 polynomial f and a degree 3 polynomial g. The mixed-volume of f and g decreases to 26 and therefore the degrees of the polynomials hα , ht to 14. For the curve being a line, i.e. r(t) = t, s(t) = 0, the polynomial f is quadratic in t and g turns out to be a linear function. The resultant polynomial also simplifies such that the mixed-volume becomes 10. In this simple setting, however, the roots αi , i = 1, 2, 3 have multiplicity 2 and therefore hα and ht are polynomials of a degree of at most 4. In the special case of a cone surface, f also becomes a degree 4 polynomial in α, whereas the degree in t does not change compared to the setting discussed above:
Expanding the polynomials Dα and |Dα | yields the second equation of the system: g(α, t) = p˜1 p′1 d2 d3 + p˜2 p′2 d1 d3 + p˜3 p′3 d1 d2 = 0, ˜ := Ap + a and p′ := with p
= ht
3 Y
∂p ∂t .
Central Surfaces Let us first consider central surfaces for which a = 0. By interpreting (18) as a polynomial in α and t we get the following system of bivariate polynomial equations: f (α, t)=A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 +
f (α, t)=A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 = 0
a0 d21 d22 d23 = 0,
g(α, t)=A1 p1 p′1 d2 d3
g(α, t)=A1 p1 p′1 d2 d3 + A2 p2 p′2 d1 d3 + A3 p3 p′3 d1 d2 = 0.
+
A2 p2 p′2 d1 d3
+
A3 p3 p′3 d1 d2
(24) = 0.
Similar to the previous case, one can show that computing the roots of the resultant requires solving univariate polynomials of a degree of at most 12 (P: ellipse / hyperbola), 8 (P:
Whereas the degree in α is 6 for f and 2 for g, the degree of t obviously depends on the curve involved. First 7
parabola) and 2 (P: line) respectively. The following table gives an overview of the results described in this section: Curve - Central Surface Ellipsoid Hyperboloids Ellipse 20 20 Hyperbola 20 20 Parabola 14 14 Line 4 4
Our considerations in this section are summarized in the following table:
Cone 12 12 8 2
Ellipse Hyperbola Parabola Line
Non-Central Surfaces
To compute the distance between a point p and a quadratic curve Q, we assume w.l.o.g. that Q is embedded on the x1 x2 -plane and centered around the origin, i.e.
f (α, t) = A1 p21 d22 + A2 p22 d21 + a0 d21 d22 = 0, g(α, t) = A1 p1 p′1 d2 + A2 p2 p′2 d1 = 0.
Q:
Since d3 = 1, we obtain the following factorization of the resultant polynomial: Res(f, g, t) =
hα
2 Y
d4i = hα
i=1
Res(f, g, α) =
ht
2 Y
i=1
(α − αi Ai )4 ,
t
Setting the derivative of the distance function equal to zero gives
i=1
f (t) = rr′ u2 + ss′ v 2 − r′ pT u − s′ pT v = 0,
A similar result can be derived for the point-hyperbola problem. If Q is a parabola, then the degree decreases and f becomes a cubic polynomial in t:
g(α, t) = a1 p′1 d2 + A2 p2 p′2 = 0. To compute the roots of the resultant polynomial one has to solve polynomials of a degree of at most 8 (P: ellipse / hyperbola), 5 (P: parabola) and 1 (P: line) respectively. If the surface is an elliptic or hyperbolic paraboloid, we have already seen (cf. 21) that f is polynomial of degree 5 in α. It is easy to see that g also simplifies to a quadratic polynomial in α:
f (t) = 2v 2 t3 + (u2 − pT v)t − pT u = 0.
+
+
(26)
In the case of a line, condition (25) simplifies to the following linear equation: f (t) = u2 t − pT v = 0.
f (α, t) = A1 p21 d22 d+ A2 p22 d21 + 2a3 pα3 d21 d22 = 0, g(α, t) =
(25)
with r′ ≡ dd rt and s′ ≡ dd st . In the case of an ellipse, condition (25) is equivalent to the following polynomial equation of degree 4 in t: pT vt4 + 2 pT u − (v 2 − u2 ) t3 f (t) = −pT v + 2 pT u + (v 2 − u2 ) t2 = 0.
f (α, t) = a1 pα1 d22 = 0,
a3 p′3 d1 d2
uT v = 0.
min(p − r(t)u − s(t)v)2 .
i=1
2 Y p2i = ht (t − ti1 )2 (t − ti2 )2 .
A2 p2 p′2 d1
q(t) = r(t)u + s(t)v,
By considering the projection p of p onto the x1 -x2 -plane, the task is reduced to a 2-dimensional problem:
Moreover, we can see from the mixed-volume function that the degree of the resultant polynomial is at most 20 (P: ellipse / hyperbola), 16 (P: parabola) and 8 (P: line) respectively. Hence, the maximum degree of the univariate polynomials to solve is 12, 8 or 2, depending on the type of the curve. In the case of a parabolic cylinder a no longer vanishes, but A1 , A3 as well as a0 are all zero. Therefore f becomes a cubic and g a linear polynomial in α:
A1 p1 p′1 d2
Parab. Cylinder 8 8 5 1
4.6. The Point-Curve Case
For non-central sufaces with a = 0 and a0 6= 0, f remains a bivariate polynomial of degree 6 in α, whereas g simplifies to a linear function in α:
2 Y
Curve - Non-Central Surface Paraboloids Ellipt. / Hyperb. Cylinders 16 12 16 12 11 8 3 2
(27)
The following table shows the maximum degrees of the univariate polynomials to solve in the Point-Curve case:
= 0.
Again our factorization shows that the degrees of the univariate polynomials to solve are much lower than the degrees of the respective resultant polynomial.
Ellipse 4
8
Point - Curve Hyperbola Parabola 4 3
Line 1
4.7. The Curve-Curve Case
we have that Res(f, g, t2 ) also vanishes for these values and the quadratic polynomial (1 ± t21 ) must be a factor. Analogously one can show that in the ellipse-parabolaas well as in the hyperbola-parabola-case the two highest coefficients of f and g vanish for t1 = ±i or t1 = ±1, respectively. Therefore the degree 16 polynomial Res(f, g, t2 ) can be represented as the product of (1 ± t21 )2 and a univariate polynomial of degree 12 in t1 . The following table gives an overview of the degree bounds in this section:
The remaining case leads to the problem of computing the distance between two quadratic curves. To simplify matters, we still assume that one curve is located as described in 4.6, i.e. P :
p(t) = r1 (t1 )u1 + s1 (t1 )v 1 ,
uT1 v 1 = 0,
whereas the other curve is in arbitrary position and orientation:
Curve-Curve
Q:
q(t) = c2 + r2 (t2 )u2 + s2 (t2 )v 2 ,
uT2 v 2
= 0.
Ellipse / Hyperbola Parabola Line
The partial derivatives of the distance function δ 2 (t1 , t2 ) = (q(t2 ) − p(t1 ))2 yield the following system of bivariate equations: ∂ r1 ∂ s1 T [q(t2 ) − p(t1 )] − (28) u1 − v1 = 0, ∂ t1 ∂ t1 ∂ s2 ∂ r2 (29) u2 + v2 = 0. [q(t2 ) − p(t1 )]T ∂ t2 ∂ t2
g(t1 , t2 ) = (1 +
t21 )g1 (t2 )
+ (1 +
t22 )g2 (t1 , t2 ),
Parabola
Line
12 9
4 3 1
4.8. Natural Quadrics, Conics and the Torus Our statements so far hold for arbitrary quadratic curves and surfaces. This class, however, contains very simple geometric objects for which the general results can be substantially improved. With the natural quadrics and conics we have identified a set of special cases, which additionally is of practical importance.
If P and Q are both ellipses, then conditions (28) and (29) can be written in the following manner: f (t1 , t2 ) = (1 + t21 )f1 (t1 , t2 ) + (1 + t22 )f2 (t1 ),
Ellipse / Hyberbola 16
(30)
4.8.1 Natural Conics
(31)
Distance queries between points and ‘linear’ objects (lines, planes) are standard problems in computational geometry and computer graphics [2]. They involve simple linear systems that we will not discuss here. Hence, the circle remains the only interesting natural conic. Computing its distance to linear objects leads to polynomial equations of a degree of at most 4. In the pointcircle case we can see from section 4.6 that computing the closest points requires solving a quadratic polynomial. The circle-plane case shows the same degree complexity. We can rotate the coordinate system so that the plane is given by x3 = 0. Then the distance between a point of the circle p(φ) and the plane P is given as
where fi and gi , i = 1, 2, are polynomials of degrees at most 2 in t1 and t2 . This representation immediately shows that every (ξ1 , ξ2 ) ∈ {−i, i}2 solves the bivariate system. Moreover, it follows that (1 + t21 )2 is a factor of the resultant Res(f, g, t2 ). Since the mixed-volume of the sparse polynomials f and g is 20, we can conclude that computing the distance between two ellipses requires solving univariate polynomials of a degree of at most 16. The same result holds for the hyperbola-hyperbola- as well as for the ellipse-hyperbola-problem. In both cases one analogously observes that the degree 4 polynomial, (1 − t21 )2 and (1 + t21 )(1 − t21 )2 respectively, is a divisor of Res(f, g, t2 ). When considering two parabolas, one obtains polynomials f and g, each having total degree 3. Therefore Res(f, g) is a polynomial of degree 9 in the non-eliminated variable. The case of two lines or a line and a parabola is also straightforward, leading to a linear and cubic equation respectively. Setting up the bivariate system (28, 29) for P being an ellipse or hyperbola and Q a line yields to polynomials f and g that have a mixed-volume of 6. By considering f and g to be polynomials in t2 , it turns out that the leading coefficient vanishes for t1 = ±i (P : ellipse) or t1 = ±1 (P : hyperbola) in both polynomials. From the resultant theory
δ(p(φ), P ) = c3 + r(cos(φ)u3 + sin(φ)v3 ). Minimizing δ over φ and using the t-parameterization for the derivative yields to a quadratic equation in t. The circle-circle case is the most difficult one since it requires the solutions of a polynomial of degree 8 [2]. Moreover, N EFF has shown that this is a lower bound on the degree complexity if the arithmetical operations are restricted to field operations and square roots [4]. 4.8.2 Natural Quadrics and the Torus Concerning the natural quadrics we first observe that the cases involving spheres and cylinders are already solved 9
since they can be reduced to the case of lower dimensional objects. The distance to a sphere is just the distance to its center minus the radius. In the case of a circular cylinder we can analogously restrict ourselves to its axis and finally subtract the radius. The same argument holds for the torus – although it is not a quadratic surface. Here, the distance can be computed by considering its main circle. Finally, we have to treat the case of a circular cone. From (19) we know that the extrema of the distance function between a point and a general cone is given by:
have therefore considered algebraic techniques that reduce the multivariate root finding problem to the univariate case. The degrees of the resulting polynomials are usually high and hence elimination methods are often not appropriate in the context of floating-point computations. However, we could show that for all cases in which the degrees of the resulting polynomials become critical, a factorization into two lower degree polynomials can be found. Thereby one of the polynomial factors is easy to solve because its roots are of high multiplicity, whereas the other is of moderate degree and can be solved using numerically stable techniques. Since in most CAD-systems the intersection curve between two quadratic surfaces is represented by NURBS, we will extend the set of trimming curves by this more complex type. Then our object class can be considered to be (approximately) closed under Boolean operations. In a final step, we will speed up the overall algorithm by setting up a bounding volume hierarchy that reduces the number of distance computations between pairs of faces.
f (α) = A1 p21 d22 d23 + A2 p22 d21 d23 + A3 p23 d21 d22 = 0. (32) In the special case of a circular cone, we have A1 = A2 and d1 = d2 . Condition (32) therefore simplifies to a polynomial equation of degree 2 in α: f (α) = A1 p21 d23 + A1 p22 d23 + A3 p23 d21 = 0. For the curve-cone problem we consider system (24). Exploiting the fact that the cone is circular, we get: f (α, t) g(α, t)
= =
A1 p21 d23 + A1 p22 d23 + A3 p23 d21 = 0, A1 p1 p′1 d3 + A1 p2 p′2 d3 + A3 p3 p′3 d1
References [1] J. Canny and I. Emiris. An efficient algorithm for the sparse mixed resultant. In Proc. Intern. Symp.Applied Algebra, Algebraic Algor. and Error-Corr. Codes, volume 263 of Lecture Notes in Computer Science., pages 89–104. SpringerVerlag, 1993. [2] D. H. Eberly. 3D game engine design : a practical approach to real-time computer graphics. Morgan Kaufmann, 2001. [3] S. Ehmann and M. Lin. Accurate and fast proximity queries between polyhedra using surface decomposition. In Computer Graphics Forum (Proc. Eurographics), 2001. [4] R. Farouki, C. Neff, and M. O’Connor. Automatic parsing of degenerate quadric-surface intersection. ACM Transactions on Graphics, 8(3):174–203, 1989. [5] E. Gilbert, D. Johnson, and S. Keerthi. A fast procedure for computing the distance between complex objects in 3d space. IEEE Journ. of Robot. Autom., 4(2):193–203, 1988. [6] K.-J. Kim. Torus and Simple Surface Intersection Based on a Configuration Space Approach. Ph.D. thesis, Dep. of Computer Science and Engineering, POSTECH, Feb. 1998. [7] J. Levin. A parametric algorithm for drawing pictures of solid objects composed of quadric surfaces. Commun. ACM, 19(10):555–563, 1976. [8] M. C. Lin and J. F. Canny. Efficient algorithms for incremental distance computation. In Proc. IEEE Internat. Conf. Robot. Autom., volume 2, pages 1008–1014, 1991. [9] B. Mirtich. Impulse-based dynamic simulation of rigid body systems. PhD thesis, University of California, Berkeley, 1996. [10] J. Reichel, E. Sch¨omer, T. Warken, and C. Lennerz. Efficient collision detection for curved solid objects. In Proc. 7th ACM Symp. on Solid Modeling and Applications, SM’ 02, 2002. [11] C. Turnbull and S. Cameron. Computing distances between nurbs-defined convex objects. In IEEE Int. Conf. of Robot. Autom., pages 3686–3690, 1998.
= 0.
For p being a line, the resultant of f and g is a polynomial of degree 4 for which we know a divisor of degree 2 (cf. our considerations in 4.5). In the case of a circle, the mixedvolume of f and g is 12. Since we have already realized that d43 is a factor of Res(f, g, t), we can restrict ourselves to solving polynomials of a degree of at most 8. The remaining case is the one in which two cones are involved. For general cones we have seen in 4.3 that polynomials of a degree of at most 4 must be solved to compute the extrema of the distance function. In the special case of two circular cones this bound cannot be improved. Our considerations above prove the following theorem. Theorem 3. The distance between two faces of natural quadratic complexes can be computed by solving univariate polynomials of a degree of at most 8.
5. Conclusions In this paper we have presented an algorithm for computing the minimal distance between two patches of quadratic surfaces bounded by quadratic curves. To find the closest points on both patches, systems of univariate and bivariate polynomial equations have to be solved. Though there are numerical techniques for finding the solutions of multivariate systems directly, our experience has shown that these methods either do not guarantee to find all roots of the system (N EWTON method) or they are not efficient enough (interval methods) to fulfill our real-time requirements. We 10
