• Intersecting biquadratic B´ zier surface patches

  • FileName: biquad.pdf [read-online]
    • Abstract: Intersecting biquadratic B´ zier surface patcheseSt´ phane Chau , Margot Oberneder†, Andr´ Galligo and Bert J¨ ttler†e e uLaboratoire J.A. Dieudonn´ Universit´ de Nice - Sophia-Antipolis, France

Download the ebook

Intersecting biquadratic B´ zier surface patches
St´ phane Chau , Margot Oberneder†, Andr´ Galligo and Bert J¨ ttler†
e e u
Laboratoire J.A. Dieudonn´ Universit´ de Nice - Sophia-Antipolis, France
e, e

Institute of Applied Geometry, Johannes Kepler University, Austria
Abstract. We present three symbolic–numeric techniques for computing the in-
tersection and self–intersection curve(s) of two B´
ezier surface patches of bidegree
(2,2). In particular, we discuss algorithms, implementation, illustrative examples
and provide a comparison of the methods.
1 Introduction
The intersection of two surfaces is one of the fundamental operations in Computer
Aided Design (CAD) and solid modeling. Closely related to it, the elimination of self–
intersections (which may arise. e.g., from offsetting) is needed to maintain the correct-
ness of a CAD model.
Tensor–product B´ zier surface patches, which are parametric surfaces defined by
vector–valued polynomials x : [0, 1]2 → 3 of certain bidegree (m, n), are extensively
used to model surfaces in CAD and solid modeling. However, even for relatively small
bidegrees m, n ≤ 3, the intersection and self–intersection loci of such patches can be
fairly complicated. Consequently, standard algorithms for surface–surface intersections
[27, 31] generally do not take the properties of special classes of such tensor–product
surfaces into account.
In the case of two general surfaces, a brute–force approach to compute the inter-
section curve(s) consists in (step 1) approximating the surface by triangular meshes and
(step 2) intersecting the planar facets of these meshes. Clearly, in order to achieve high
accuracy, a very fine approximation with a mesh may be needed. Alternatively, one may
consider to choose another, more complicated representation, where the basic elements
are capable of capturing more of the geometric features. For instance, quadratic trian-
gular B´ zier patches were already studied in [2, 3, 13, 33, 34] and one might try to use
biquadratic patches. Clearly, this approach would need robust intersection algorithms
for the more complicated basic elements.
In this paper we address the computation of the intersection curve of two surface
patches of bidegree (2,2), i.e., biquadratic patches. Our aim is to compute the intersec-
tion by using – as far as possible – symbolic techniques, in order to avoid problems with
numerical robustness.
The remainder of the paper is organized as follows. After some preliminaries, Sec-
tions 3 to 5 present three different techniques for computing the intersection curves,
which are based on resultants, on approximate implicitization (which is an important
2 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
research topic in the GAIA II project), and on intersections of parameter lines, respec-
tively. Section 6 discusses the computation of self–intersections. We apply the three
techniques to three representative examples and report the results in Section 7. Finally,
we conclude this paper.
2 Intersection and self–intersection curves
We consider the intersection curves of two biquadratic B´ zier surfaces x(u, v) and
y(r, s), both with parameter domains [0, 1]2 . They are assumed to be given by their
parametric representations with rational coefficients (control points). More precisely,
these representations have the form
2 2
x(u, v) = ci,j Bi (u)Bj (v) (1)
i=0 j=0
with certain rational control points ci,j ∈ Q3 and the quadratic Bernstein polynomials
Bj (t) = 2 ti (1 − t)2−i (and similarly for the second patch y(r, s)).
The intersection curve is defined by the system of three non–linear equations
x(u, v) = y(r, s) (2)
which defines the intersection as a curve (in the generic case) in [0, 1] 4 . Similarly, self
intersections of one of the patches are characterized by
x(u, v) = x(¯, v ).
u ¯ (3)
In this case, the set of solutions contains the 2–plane u = u∗ , v = v ∗ as a trivial
While these equations could be solved by using numerical methods, we plan to ex-
plore how far it is possible to compute the intersections by using symbolic computations,
in order to avoid rounding errors and robustness problems.
The “generic” algorithm for computing the (self–) intersection curve(s), consists of
three steps:
1. Find at least one point on each component of the intersection,
2. trace the segments of the intersection curve, and
3. collect and convert the segments into a format that is suitable for further processing
(depending on the application).
We will focus on the first step, since the second step is a standard numerical problem,
and step 3 depends on the specific background of the problem. Several parts of the
intersection curve may exist. Some possible types are shown in Fig. 1 in the parameter
domain of a B´ zier surface x(u, v). Points with horizontal or vertical tangent are called
turning points, and intersections with the boundaries of the patches generate boundary
points. Note that also isolated points (where both surfaces touch each other) may exist.
Intersecting biquadratic patches 3
boundary points
turning points
PSfrag replacements
Fig. 1. Intersection curves in one of the parameter domains.
3 A resultant–based approach
In this section, we will use the resultant to compute the intersection locus between
x(u, v) and y(r, s). We consider the algebraic variety
C = {(u, v, r, s) | x(u, v) = y(r, s)} (4)
and we will suppose that C ∩ [0, 1]4 is a curve.
3.1 Resultant basics
Let f1 , f2 and f3 be three polynomials in two variables with given monomial supports
and N the number of terms of these 3 supports. For each i ∈ {1, 2, 3} we denote by
coeffs(fi ) the sequence of the coefficients of fi . The resultant of f1 , f2 and f3 is, by
definition, an irreducible polynomial R in N variables with the property, that
R(coeffs(f1 ), coeffs(f2 ), coeffs(f3 )) = 0 (5)
if and only if these 3 polynomials have a common root in a specified domain D. For a
more precise description of resultants, see e.g. [4, 10, 11].
In our application to surface–surface–intersections, the resultant can be used as a
projection operator. Indeed, if f1 , f2 and f3 are the three components of x(u, v)−y(r, s)
which are considered as polynomials in the two variables r and s, then the resultant
of f1 , f2 and f3 is a polynomial R(u, v) and it gives an implicit plane curve which
corresponds to the projection of C in the (u, v) parameters. More precisely, if f 1 , f2 and
f3 are generic, then the two sets
(u, v) ∈ [0, 1]2 | R(u, v) = 0 (6)
(u, v) ∈ [0, 1]2 | ∃(r, s) ∈ D : x(u, v) = y(r, s) (7)
are identical. Several families of multivariate resultants have been studied and some
implementations are available, see [7, 25].
4 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
3.2 Application to the intersection problem
A strategy to describe the intersection between x(u, v) and y(r, s) consists in projecting
C on a plane (by using the resultant). Many authors propose to project C on the (u, v)
(or (r, s)) plane and then the resulted plane curve is traced (see [19] and [23] for the
tracing method) and is lifted to the 3D space by the corresponding parameterization.
Note that this method can give some unwanted components (the so called ”phantom
components”) which are not in x([0, 1]2 ) ∩ y([0, 1]2 ). So, another step is needed to cut
off the extraneous branches. This last part can be done with a solver for multivariate
polynomial systems (see [28]) or an inversion of parameterization (see [5]).
As an alternative to these existing approaches, we propose to project the set C onto
the (u, r) space. Note that, in the equations x(u, v) = y(r, s), the two variables v and s
are separated, so they can be eliminated via a simple resultant computation. It turns out
that such a resultant can be computed via the determinant of a Bezoutian matrix (see
[18]). First, consider the (3, 3) determinant:
x(u, v) − x(u, v1 ) y(r, s) − y(r, s1 )
b = det x(u, v) − y(r, s), , . (8)
v − v1 s − s1
The determinant b is a polynomial and its monomial support with respect to (v, s) is
S = {1, v, s, vs} and similarly for (v1 , s1 ), where S1 = {1, v1 , s1 , v1 s1 }. So, a mono-
mial of b is a product of an element of S and of an element of S1 . Then, we form
the 4 × 4 matrix whose entries are the coefficients of b indexed by the product of the
two sets S and S1 . This matrix contains only the variables u and r and is called the Be-
zoutian matrix. In our case, its determinant is a polynomial in (u, r) equal to the desired
resultant R(u, r) (deg(R)=24 and degu (R)=degr (R)=16) and it gives an implicit curve
which corresponds to the projection of C in the (u, r) space.
Then, we analyse the topology of this curve (see [20] and [35]) and we trace it (see
[19] and [23]). Finally, for each (u0 , r0 ) ∈ [0, 1]2 such that R(u0 , r0 ) = 0, we can
determine if there exists a pair (v0 , s0 ) ∈ [0, 1]2 such that x(u0 , v0 ) = y(r0 , s0 ) (solve
a polynomial system of three equations with two separated unknowns of bidegree (2,2))
and thus we can avoid the problem of the phantom components (see Fig. 2). We lift the
obtained points in the 3D space to give the intersection locus. Note that this method can
also give the projection of C in the (v, s) space by the same kind of computation.
4 Approximate implicitization by a quartic surface
In this section, we apply the technique of approximate implicitization to compute the
intersection of two biquadratic patches.
4.1 Approximate implicitization
The implicitization problem – which consists in finding an implicit equation (an alge-
braic representation) for a given parameterized rational surface – can be adressed by
using several approaches, e.g., using resultants or Groebner bases [10, 11, 21]. How-
ever, the implicitization is very time consuming because of the degree of the implicit
Intersecting biquadratic patches 5
r r
u u
Fig. 2. Projection of C in the (u, r) space with (left) and without (right) phantom components.
This curve corresponds to the example of Figure 6, page 14.
equation: for a generic parameterized surface of bidegree (n1 ,n2 ), the implicit equation
has degree 2n1 n2 . Also, all rational parametric curves and surfaces have an algebraic
representation, but the reverse is not true; the relationship between the parametric and
the algebraic representations can be very complex (problem of ”phantom components”).
Thus, we can try to find an algebraic approximation of a given parameterized surface for
which the computation is more efficient and which contains less phantom components.
Consider a polynomial parameterized surface x(u, v) with the domain [0, 1] 2 , and
let d be a positive integer (the degree of the approximate implicit equation) and ≥ 0
(the tolerance). Following [15], the approximate implicitization problem consists in
finding a non–zero polynomial P ∈ R[x, y, z] of degree d such that
∀(u, v) ∈ [0, 1]2 , P (x(u, v) + α(u, v) g(u, v)) = 0 (9)
with |α(u, v)| ≤ and ||g(u, v)||2 = 1. Here, α is the error function and g is the
direction for error measurement, e.g., the unit normal direction of the surface patch.
4.2 Approximate implicitization of a biquadratic surface
The main question of the approximate implicitization problem is how to choose the de-
gree. A key ingredient for this choice seems to be the topology, especially if the initial
surface has self–intersections. The use of degree 4 was suggested by Tor Dokken; after
several experiments he concluded that the algebraic surfaces of degree 4 provide suffi-
ciently many degrees of freedom to approximate most cases encountered in practice. In
the case of a biquadratic surface, where the exact implicit equation has degree 8, using
degree 4 seems to be a reasonable trade-off.
We describe two methods for approximate implicitization by a quartic for a bi-
quadratic surface. The approximate implicit equation is
4 4−i 4−i−j
P (x, y, z) = bijk xi y j z k (10)
i=0 j=0 k=0
6 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
PSfrag replacements
An interior conic
A border conic
Fig. 3. Characterization of a conic in a biquadratic patch by 9 points
with the unknown coefficients b = (b000 , b100 , . . . , b004 ) ∈ R35 . Let β(u, v) be the
vector formed by the tensor–product Bernstein polynomials of bidegree (8,8).
Dokken’s method. This method, which is described in more detail in [15], proceeds
as follows:
1. Factorize P (x(u, v)) = (Db)T β(u, v) where D is a 81 × 35 matrix.
2. Generate a singular values decomposition (SVD) of D.
3. Choose b as the vector corresponding to the smallest singular value of D.
Note that this method is general and does not use the fact that we have a biquadratic
surface. Hereafter, we use an adapted method based on the geometry of the surface of
bidegree (2,2). Also, the computation of the singular value decomposition needs floating
point numbers.
Geometric method using evaluation: This approach consists in constructing some
pertinent geometrical constraints to give a linear system of equations (with the un-
knowns b000 , b100 , . . . , b004 ), and then solving the resulting system by a singular values
decomposition. In our method, we characterize some conics, especially the four border
conics and two interior conics:
C1 = x([0, 1] × {0}), C2 = x([0, 1] × {1})
C3 = x({0} × [0, 1]), C4 = x({1} × [0, 1]) (11)
1 1
C5 = x({ 2 } × [0, 1]), C6 = x([0, 1] × { 2 })
Lemma 1. If the quartic surface {P = 0} contains 9 points of any of the 6 conics C i ,
then Ci ⊂ {P = 0}, see Fig. 3.
Proof. Ci is of degree 2 and P is of degree 4, so by B´ zout’s theorem, if there are more
than 8 elements in Ci ∩ {P = 0}, then Ci ⊂ {P = 0}.
Using this geometric observation, we construct a linear system and solve it approx-
imately via SVD; this leads to an algebraic approximation of x(u, v) by a degree 4
Intersecting biquadratic patches 7
4.3 Application to the intersection problem
In order to compute the intersection curves, we apply the approximate implicitization
to one of the patches and compose it with the second one. This leads to an implicit
representation of the intersection curve in one of the parameter domains, which can
then be traced and analyzed using standard methods for planar algebraic curves.
These two approximate implicitization methods are very efficient and suitable for
general cases, but the results are not always satisfactory. When the given biquadratic
patch is simple (i.e. with a certain flatness and without singularity and self–intersection)
the approximation is very close to the initial surface. So, to use this method for a general
biquadratic surface, we combine it, if needed, with a subdivision method (Casteljau’s
algorithm). The advantage is twofold, we exclude domains without intersections (by
using bounding boxes) and avoid some unwanted configurations with a curve of self-
intersection (use Hohmeyer’s criterion [22]). For more complicated singularities, the
results are definitively not satisfactory.
Note that even if we have a good criterion in the subdivision step, we still may
have problems with phantom components (but in general fewer), so we have to cut
off the extraneous branches as in the resultant method. This has to be done carefully
in order to not discard points which do not correspond to phantom components. As
another drawback – because of the various approximations – it is rather difficult to
obtain certified points on the intersection locus. The use of approximate implicitization
is clearly a numerical method, and it can give only approximate answers, even in the
case of exact input.
5 Tracing intersections of parameter lines
In order to be able to trace the (self–) intersection curve(s), we have to find at least one
point for each segment. We generate these points by intersecting the parameter lines of
the first B´ zier surface with the second one (see also [22]).
5.1 Intersection of a parameter line
A parameter line of x(u, v) for a constant rational value u = u0 takes the form
p(v) = x(u0 , v) = a0 (u0 ) + a1 (u0 ) v + a2 (u0 ) v 2
with certain rational coefficient vectors ai ∈ Q3 . It is a quadratic B´ zier curve, hence
we can represent it as the intersection of a plane and a quadratic cylinder, see Fig. 4,
left. Since we are only interested in the intersection of these two surfaces in a certain
region, we introduce two additional bounding planes π1 and π2 . In the particular case
that the parameter line is a straight line, we represent it as an intersection curve of two
orthogonal planes.
In order to compute the intersection of the parameter line with the second surface
patch y(r, s), we use the following algorithm.
1. Describe the parameter line as the intersection of a plane and a cylinder.
8 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
PSfrag replacements
x(u, v) r1
kπ 2
PSfrag replacements
kπ 1 kπ 2
s0 s1 s
Fig. 4. Left: Representation of a parameter line as the intersection of a plane and a quadratic
cylinder. Right: Identifying the intervals with feasible values of s.
2. Intersect the plane with the second patch y(r, s) and compute the intersection I.
3. Restrict the intersection curve(s) I to the region of interest.
4. Intersect the cylinder with the restricted intersection curve(s).
The four steps of the algorithm will now be explained in some more detail.
Defining the plane, the cylinder and the two bounding planes. The parameter line
and its three control points are coplanar. For computing the normal vector n of the
plane, we have to evaluate the cross product of two difference vectors of the control
points. The plane is given by the zero set of a linear polynomial
π(u0 )(x, y, z) = e0 (u0 ) + n1 (u0 ) x + n2 (u0 ) y + n3 (u0 ) z. (12)
By extruding the parameter line in the direction of the normal vector of the plane, we
obtain the parametric form of the quadratic cylinder, which intersects the plane orthog-
w(u0 )(p, q) = x(u0 , p) + q · n. (13)
The implicitation of the cylinder is slightly more complicated. There exist two possibil-
ities: we can either use Sylvester resultants or the method of comparing coefficients. In
both cases we will get an equation of the form
ζ(u0 )(x, y, z) := a0 (u0 ) + a1 (u0 ) x + a2 (u0 ) y + a3 (u0 ) z
+ a4 (u0 ) x y + a5 (u0 ) x z + a6 (u0 ) y z
+ a7 (u0 ) x2 + a8 (u0 ) y 2 + a9 (u0 ) z 2 = 0. (14)
Now we have both the plane and the cylinder in their implicit representation. Note that
this is a semi-implicit representation in the sense of [8].
Intersecting biquadratic patches 9
If the parameter line degenerates into a straight line, then we choose two planes
through it which intersect orthogonally. Note that we use exact rational arithmetic, in
order to avoid any robustness problems.
Finally, we create the two planes π1 (x, y, z) and π2 (x, y, z) which bound the pa-
rameter line. For instance, one may choose the two normal planes of the parameter line
at its boundary points; this choice is always possible, provided that the curve segment
is not too long (which can be enforced by using subdivision). Alternatively one may
use the planes spanned by the boundary curves, but these planes may have an additional
intersection with the parameter line in the region of interest.
Intersection of the plane and the second patch y(r, s). Substituting the second
B´ zier surface y(r, s) into the equation (12) of the plane leads to a biquadratic equation
in r and s. We can treat it as a quadratic polynomial in r with coefficients depending
on s.
π(y(r, s)) = a(s) r 2 + b(s) r + c(s) = 0. (15)
For each value of s, we obtain two solutions r1 (s) and r2 (s) of the form
b(s) b(s)2 c(s)
r1,2 (s) = − ± d(s) with d(s) = − . (16)
2 a(s) 4 a(s)2 a(s)
These solutions parameterize the two branches of the intersection curve I in the rs–
parameter domain of the second patch. By solving several quadratic equations we de-
termine the intervals Si,j ⊂ [0, 1], where d(s) ≥ 0 and 0 ≤ ri (s) ≤ 1 holds; this leads
to a (list of) feasible domain(s) (i.e., intervals) for each branch of the intersection curve.
By composing (16) with y we obtain the two branches k1 (s) and k2 (s) of the
intersection curve I,
1 d(s)
k1,2 (s) = y(r1,2 (s), s) = h(s) ± l(s) + d(s)m(s) (17)
a(s)2 a(s)
where the components of h(s), l(s) and m(s) are polynomials of degree 6, 4, and 2,
Restriction to the region of interest. Since the region of interest is located between
the planes π1 (x, y, z) and π2 (x, y, z), the two inequalities
π1 (x, y, z) ≥ 0 and π2 (x, y, z) ≤ 0 (18)
have to be satisfied. By intersecting each bounding plane with the second B´ zier surface
y(r, s) in a similar way as described for π(u0 )(x, y, z), we obtain
kπ1 (s) := π1 (y(r(s), s)) ≥ 0 and kπ2 (s) := π2 (y(r(s), s)) ≤ 0 (19)
This leads to additional constraints for the feasible values of the parameter s. For each
branch of the intersection curve we create the (list of) feasible domain(s) and store
10 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
it. The bounds of the intervals can be computed by solving three systems of two bi-
quadratic equations or – equivalently – by solving a system of three polynomials of
degree 8, which are obtained after eliminating the parameter r. Here, we represent the
polynomials in Bernstein–B´ zier form and use a B´ zier–clipping–type technique see
e e
[17, 28, 29, 31], applied to floating point numbers.
Example 2. For a parameter line u = u0 of two biquadratic B´ zier surface patches
x(u, v) and y(r, s), Fig. 4, right, shows the rs–parameter domain of the second patch.
Only the first branch r1 (s) of the intersection curve is present. The bounds 0 ≤ r ≤ 1
do not impose additional bounds on s in this case. However, the intersection with the
bounding planes π1 and π2 produces two additional curves, which have to be intersected
with the curve s = r1 (s), leading to two bounds s0 and s1 of the feasible domain.
Intersection of the cylinder and the intersection curves. We substitute the parametric
representation of the intersection curve into the implicit equation (14) of the cylinder
and obtain
ζ(u0 )(s) = p1 (s) + p2 (s) d(s) + p3 (s) d(s) +
3 4
+ p4 (s) d(s) + p5 (s) d(s) = 0 (20)
where the polynomials pj (s) are of degree 12. In order to eliminate the square root, we
use the following trick. We split ζ(s, d(s)) = A − B, where A and B contain all even

and odd powers of d, respectively. The equation A − B = 0 is then replaced with
A2 · d(s) − (B · d(s))2 = 0. This leads to a polynomial of degree 24 in one variable.
After factoring out the discriminant, we obtain a polynomial of degree 16 in s. Note that
this agrees with the theoretical number of intersections of a biquadratic surface, which
has algebraic order 8, with a quadratic curve.
Finally, we solve this polynomial within all the feasible intervals of s, which were
detected in the previous steps. Until this point we used symbolic computations. Now
– after generating the Bernstein–B´ zier representation – we change to floating-point
numbers and use a B´ zier–clipping–type method to find all roots within the feasible
domain(s). These roots correspond to intersection points of the parameter line of the
first patch with the second patch.
5.2 Global structure of the intersection curve
For each value u = u0 , the parameter line x(u0 , v) has a certain number of intersection
points with the second patch. If u0 varies continuously, then the number of intersection
points may change only if
(1) one of the intersection points is at the boundary of one of the patches (boundary
points) or
(2) the parameter line of the first patch touches the second patch (turning points).
Intersecting biquadratic patches 11
The algorithm for analyzing the global structure of the intersection curve proceeds in
two steps: First we detect those values of u0 where the number of intersection points
changes, and order them. This leads to a sequence of critical u 0 – values,
(0) (1) (K)
0 = u 0 < u0 < . . . < 1 = u 0 . (21)
In the second step, we analyze the intersection of the parameter lines u 0 = (u0 +
u0 )/2 with the second patch. Since the number of intersection points between any
two critical values remains constant, we can now either trace the segment using conven-
tional techniques for tracing surface–surface intersections (see [23]) or generate more
points by analyzing more intersections with parameter lines.
In the remainder of this section we address the computation of the critical u 0 values.
Boundary points. Such points correspond to intersections of the boundary parameter
lines of one surface with the other one. In order to compute them, we apply the algo-
rithm for intersecting parameter lines with a biquadratic patch to the 2 · 4 boundary
parameter lines of the two surfaces.
Turning points. We consider the turning points of x(u, v) in respect to u. Let y r and ys
denote the partial derivatives of y(r, s). Several possibilities for computing the turning
points exist.
1. The two surfaces x(u0 , v) and y(r, s) intersect, x(u0 , v) = y(r, s), and the tangent
vector of the parameter line lies in the tangent plane of the second patch,
xu · (yr × ys ) = 0. (22)
These conditions lead to a system of four polynomial equations for four unknowns,
which has to be solved for u.
2. By using the previous geometric result, we may eliminate the variable v, as follows.
First, the plane spanned by the parameter line has to contain the point y(r, s),
π(u0 )(y(r, s)) = 0, (23)
which gives an equation of degree (6, 2, 2) in (u0 , r, s). Second, the cylinder has to
contain the point,
ζ(u0 )(y(r, s)) = 0, (24)
which leads to an equation of degree (16, 4, 4). Finally, the tangent vector of the
parameter line has to be contained in the tangent plane of the second patch. Since
the tangent of the parameter line is parallel to the cross product of the gradient of
the plane and the gradient of the cylinder, the third condition gives an equation of
degree (18, 5, 5),
det [yr , ys , π(u0 )(y(r, s)) × ζ(u0 )(y(r, s))] = 0. (25)
For solving either of these two systems of polynomial equations, we use again a B´ zier–
clipping–type algorithm [17, 28, 31].
12 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
6 Self–intersections of biquadratic surface patches
In order to detect the self–intersection curves of any of the two patches, the methods
for surface–surface intersections have to be modified. Note that the computation of the
self–intersection locus by using the approximate implicitization was already treated in
6.1 Resultant-based method
In the parameter domain [0, 1]4 , the self–intersection curve of the first patch forms the
(u1 , v1 , u2 , v2 ) ∈ [0, 1]4 | (u1 , v1 ) = (u2 , v2 ) and x(u1 , v1 ) = x(u2 , v2 ) . (26)
This locus is the real trace of a complex curve. We assume that it is either empty or of
dimension 0 or 1. We do not consider degenerate cases, such as a plane which is covered
twice. In the examples presented below (see Section 7), the self–intersection locus is a
curve in R4 .
We use the following change of coordinates to discard the unwanted trivial compo-
nent (u1 , v1 ) = (u2 , v2 ). Let (u2 , v1 ) be a pair of parameters in [0, 1]2 , (l, k) ∈ R2 and
let u1 = u2 + l, v2 = v1 + lk. If we suppose that we have (u1 , v1 ) = (u2 , v2 ), then
l = 0. Hence x(u1 , v1 ) = x(u2 , v2 ) if and only if x(u2 + l, v1 ) = x(u2 , v1 + lk). We
suppose now that (u2 , v1 , l, k) verifies this last relation.
Let T (u2 , v1 , l, k) be the polynomial 1 [x(u2 + l, v1 ) − x(u2 , v2 + lk)], its degree
in (u2 , v1 , l, k) is (2, 2, 1, 2) and the monomial support with respect to (l, k) contains
only k 2 l, k, l and 1. We can decrease the degree by introducing
˜ 1
T (u2 , v1 , m, k) = mT (u2 , v1 , , k). (27)
Then in T (u2 , v1 , m, k), the monomial support in (m, k) consists only of 1, m, k 2 and
km. So, we can write T in a matrix form:
 
  1
a1 (u2 , v1 ) b1 (u2 , v1 ) c1 (u2 , v1 ) d1 (u2 , v1 ) 
m 
T (u2 , v1 , m, k) =  a2 (u2 , v1 ) b2 (u2 , v1 ) c2 (u2 , v1 ) d2 (u2 , v1 )   2  (28)
k 
a3 (u2 , v1 ) b3 (u2 , v1 ) c3 (u2 , v1 ) d3 (u2 , v1 )
By Cramer’s rule, we get
D2 D3 D4
m= , k2 = , and km = (29)
D1 D1 D1
b1 c1 d 1 −a1 c1 d1 b1 −a1 d1 b1 c1 −a1
D1 = b2 c2 d2 , D2 = −a2 c2 d2 , D3 = b2 −a2 d2 , D4 = b2 c2 −a2 .
b3 c3 d 3 −a3 c3 d3 b3 −a3 d3 b3 c3 −a3
2 2
Let Q(u2 , v1 ) be the polynomial Q = D4 D1 − D2 D3 .
Intersecting biquadratic patches 13
Fig. 5. A self–intersection of a surface with a cuspidal point
Lemma 3. The implicitly defined curve (u2 , v1 ) ∈ [0, 1]2 | Q(u2 , v1 ) = 0 is the pro-
jection of the self–intersection locus (given by the set (26) but in C4 ) into the parameters
domain (u2 , v1 ) ∈ [0, 1]2 .
Proof. Q(u2 , v1 ) = 0 is the only algebraic relation (of minimal degree) between u2 and
v1 such that ∀(u2 , v1 ) ∈ [0, 1]2 , Q(u2 , v1 ) = 0 ⇒ ∃(m, k) ∈ 2 , T (u2 , v1 , m, k) = 0.
This lemma provides a method to compute the self–intersection locus, we just have
to trace the implicit curve Q(u2 , v1 ) = 0 and for every point (u2 , v1 ) on this curve,
we obtain by continuation the corresponding point (u1 , v2 ) ∈ [0, 1]2 if it exists (see the
results on Fig. 9). So it suffices to characterize the bounds of these segments of curves.
6.2 Parameter-line-based method
For computing the self–intersection curves, we use the same algorithm as described
in Section 5. We intersect the surface x(u0 , v) with itself x(r, s). In this case, both the
”plane” equation (23) and the ”cylinder” equation (24) contain the linear factor (r−u 0 ),
which has to be factored out. The computation of turning points as in section 5.2 leads
us to two different types: the usual ones and cuspidal points (see Fig. 5).
7 Examples
The three methods presented in this paper (using resultants, via approximate implici-
tization, and by analyzing the intersections with parameter lines) work well for most
standard situations usually encountered in practice. In this section, we present three
representative examples. Additional ones are available at [24].
14 S. Chau, M. Oberneder, A. Galligo and B. J¨ ttler
Fig. 6. First example. Left and center: Result of the resultant method after and before eliminating
phantom branches. Right: result of the approach using approximate implicitization.
x(u, v) y(r, s)
v s
PSfrag replacements PSfrag replacements
u r
Fig. 7. First example. The intersection curves in the parameter domains of both surface patches,
generated by the parameter–line based technique. Boundary points and turning points have been
marked by grey circles.
7.1 First example
We consider two biquadratic surfaces with an open and a closed component of the
intersection curve. The two surfaces have the control points
 1 3 3 1 3 7
  2 1 2 3 1 2 
7 , 0, 5 5, 5, 4 1, 0, 10 7, 7, 5 5 , 10 , 3 1, 0, 4
 3 4 2 2 3 1 6 3 5   3 4 2 1 1 5 3 2 
 8, 9, 3 3, 4, 3 7 , 8 , 7  and  8 , 9 , 3 3, 2, 1 7, 8, 7  .
1 6 4 3 7 3 7 7 5 1 6 3 3 7 5 7 4 1
5, 7, 7 4, 8, 4 8, 9, 8 5, 7, 7 4, 8, 8 8, 7, 2
x(u,v) y(r,s)
By using the resultant method, a phantom component appears (see Fig. 6, center). It can
be cut off as described in Section 3.2 (see Fig. 6, left).
Similar to the resultant method, the approximate implicitization produces a phantom
component (see Fig. 6, right). However, when we cut it off, we obtain only very few
certified points on the intersection locus as described in section 4.3.
The parameter-line-based approach finds both parts of the intersection curve, but
no phantom components. One segment is closed and has two turning points with respect

Use: 0.1783