13.472J/1.128J/2.158J/16.940J
COMPUTATIONAL GEOMETRY
Lectures 10 - 12
N. M. Patrikalakis
Massachusetts Institute of Technology
Cambridge, MA 02139-4307, USA
Copyright c 2003 Massachusetts Institute of Technology
Contents
10.1 Overview of intersection problems . . . . . . . . . . . . . . . . . . . . . . . . . 3
10.2 Intersection problem classi cation . . . . . . . . . . . . . . . . . . . . . . . . . 5
10.2.1 Classi cation by dimension . . . . . . . . . . . . . . . . . . . . . . . . 5
10.2.2 Classi cation by type of geometry . . . . . . . . . . . . . . . . . . . . . 5
10.2.3 Classi cation by number system . . . . . . . . . . . . . . . . . . . . . . 6
10.3 Point/point \intersection" . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
10.4 Point/curve intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
10.4.1 Point/Implicit curve intersection . . . . . . . . . . . . . . . . . . . . . 8
10.4.2 Point/Parametric curve intersection . . . . . . . . . . . . . . . . . . . . 10
10.4.3 Point/Procedural parametric (o set, evolute, etc.) curve intersection . 12
10.5 Point/surface intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
10.5.1 Point/Implicit (usually algebraic) surface intersection . . . . . . . . . . 13
10.5.2 Point/Rational polynomial surface intersection . . . . . . . . . . . . . . 13
10.5.3 Point/Procedural surface intersection . . . . . . . . . . . . . . . . . . . 19
10.6 Curve/curve intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
10.6.1 Case D3: RPP/IA curve intersection . . . . . . . . . . . . . . . . . . . 20
10.6.2 Case D1: RPP/RPP Curve Intersection . . . . . . . . . . . . . . . . . 27
10.6.3 Case D2/D5: RPP/PP and PP/PP Curve Intersections . . . . . . . . . 28
10.6.4 Case D6: PP/IA Curve Intersection . . . . . . . . . . . . . . . . . . . . 28
10.6.5 Case D8: IA/IA Curve Intersection . . . . . . . . . . . . . . . . . . . . 29
10.7 Curve/surface intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
10.7.1 Case E3: RPP Curve/IA Surface Intersection . . . . . . . . . . . . . . 30
10.7.2 Case E1: RPP Curve/RPP Surface Intersection . . . . . . . . . . . . . 31
10.7.3 Case E2/E6: RPP/PP, PP/PP Curve/Surface Intersection . . . . . . . 31
10.7.4 Case E7: PP Curve/IA Surface Intersection . . . . . . . . . . . . . . . 31
1
10.7.5 Case E11: IA Curve/IA Surface Intersection . . . . . . . . . . . . . . . 32
10.7.6 IA Curve/RPP Surface Intersection . . . . . . . . . . . . . . . . . . . . 33
10.8 Surface/Surface Intersections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
10.8.1 Case F3: RPP/IA Surface Intersection . . . . . . . . . . . . . . . . . . 34
10.8.2 Case F1: RPP/RPP Surface Intersection . . . . . . . . . . . . . . . . . 42
10.8.3 Case F8: IA/IA Surface Intersection . . . . . . . . . . . . . . . . . . . 46
10.9 Nonlinear Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
10.9.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
10.10Local Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
10.11Classi cation of Global Solution Methods . . . . . . . . . . . . . . . . . . . . . 54
10.12Subdivision Method (Projected Polyhedron Method) . . . . . . . . . . . . . . 55
10.13Interval Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
10.13.1Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
10.13.2De nitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
10.13.3Interval Arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
10.13.4Algebraic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
10.13.5Rounded Interval Arithmetic and its Implementation . . . . . . . . . . 62
10.14Interval Projected Polyhedron Algorithm . . . . . . . . . . . . . . . . . . . . . 65
10.15Interval Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Bibliography 72
Reading in the Textbook
Chapters 4 and 5, pp. 73 - 160
2
Lectures 10 - 12
Intersection Problems, Nonlinear
Solvers and Robustness Issues
10.1 Overview of intersection problems
Intersections are fundamental in computational geometry, geometric modeling and design,
analysis and manufacturing applications. Examples of intersection problems include:
Shape interrogation (eg. visualization) through contouring (intersection with a series of
parallel planes, coaxial cylinders, and cones etc.)
Numerical control machining (milling) involving intersection of o set surfaces with a
series of parallel planes, to create machining paths for ball (spherical) cutters.
Representation of complex geometries in the \Boundary Representation" scheme; for
example, the description of the internal geometry and of structural members of cars,
airplanes, ships, etc, involves
{ Intersections of free-form parametric surfaces with low order algebraic surfaces
(planes, quadrics, torii).
{ Intersections of low order algebraic surfaces.
in a process called boundary evaluation, in which the Boundary Representation is cre-
ated by \evaluating" the Constructive Solid Geometry model of the object. During this
process, intersections of the surfaces of primitives (see Figure 10.1) must be found during
Boolean operations.
Boolean opertations on point sets A, B include (see Figure 10.2)
Union: A[B,
Intersection: A\B, and
Di erence: A B.
3
Figure 10.1: Primitive solids.
)))
Figure 10.2: Example of a Boolean operation: union.
All such operations involve intersections of surfaces to surfaces. In order to solve general
surface to surface intersection problems, the following auxiliary intersection problems (similar
to distance computation problems used in CAM for inspection of manufactured objects) need
to be considered
point/point (P/P)
point/curve (P/C)
point/surface (P/S)
curve/curve (C/C)
curve/surface (C/S)
All above types of intersection problems are also useful in geometric modeling, robotics,
collision avoidance, manufacturing simulation, scienti c visualization, etc.
When the geometric elements involved in intersections are nonlinear (curved), intersection
problems typically reduce to solving complex systems of nonlinear equations, which may be
either polynomial or more general in character. Solution of nonlinear systems is a very com-
plex process in general in numerical analysis and there are specialized textbooks on the topic.
However, geometric modeling applications pose severe robustness, accuracy, automation, and
e ciency requirements on solvers of a nonlinear systems as we will see later. Therefore, geo-
metric modeling researchers have developed specialized solvers to address these requirements
explicitly using geometric formulations.
4
10.2 Intersection problem classi cation
Intersection problems can be classi ed according to the dimension of the problems and ac-
cording to the type of geometric equations involved in de ning the various geometric elements
(points, curves and surfaces). The solution of intersection problems can also vary according to
the number system in which the input is expressed and the solution algorithm is implemented.
10.2.1 Classi cation by dimension
P/P, P/C, P/S
C/C, C/S
S/S
10.2.2 Classi cation by type of geometry
polynomial Procedural
Polynomial
(algebraic)
(eg. f(x;y) = 0; z = 0)
Implicit
Rational
Parametric
(eg. R = R(t) )
Figure 10.3: Curve geometry classi cation
1. Points
Explicit: R = R0; R = [x;y;z]
Procedural: Intersection of two procedural curves, procedural curve and surface, or
three procedure surfaces, eg. o set or blending surfaces.
Algebraic: f(R) = g(R) = h(R) = 0; where f, g, and h are polynomials.
2. Curves
A classi cation of curves is illustrated in Figure 10.3.
Parametric: R = R(t) A t B
(a) Rational Polynomials (eg: NURBS, rational B ezier).
(b) Procedural, eg: o sets, evolutes, ie. the locus of the centers of curvature of a
curve.
Implicit: These require solution of (usually nonlinear) equations
5
(a) Algebraics (polynomial)
f(R) = g(R) = 0 space curves
z = 0; f(x;y) = 0 planar curves
(b) Procedural o sets (eg. non-constant distance o sets involving convolution, see
Pottmann 1997)
3. Surfaces
Parametric R = R(u;v) where u;v vary in some nite domain, the parametric
space.
(a) (Rational) Polynomial (eg: NURBS, B ezier, rational B ezier etc.)
(b) Procedural
{ o sets
{ blends
{ generalized cylinders
Implicit: Algebraics
f(R) = 0
where f is a polynomial.
10.2.3 Classi cation by number system
In our discussion of intersection problems, we will refer to various classes of numbers:
integer numbers;
rational numbers, m=n, n6= 0, where m;n are integers;
oating point numbers in a computer (which are a subset of the rational numbers);
radicals of rational numbers, eg.
q
m=n, n6= 0, where m;n are integers;
algebraic numbers (roots of polynomials with integer coe cients);
transcendental (e, , trigonometric, etc.) numbers.
real numbers;
interval numbers, [a;b], where a;b are real numbers;
rounded interval numbers, [c;d], where c;d are oating point numbers.
Issues relating to oating point and interval numbers a ecting the robustness of intersection
algorithms are addressed in the next section on nonlinear solvers.
6
10.3 Point/point \intersection"
Check if jR1 R2j< , where represents the maximum allowable tolerance.
Choice of \tolerances" in a geometric modeller is di cult{an open question.
Lack of transitivity, see Figure 10.4:
R1 = R2
R2 = R3 ) R1 6= R3
R2
R3
.
.
.
R1
Figure 10.4: Intersections of points within a tolerance is intransitive.
What should re ect?
7
10.4 Point/curve intersection
10.4.1 Point/Implicit curve intersection
R0 \fz = 0; f(x;y) = 0g
where f(x;y) is usually a polynomial (and f(x;y) = 0 represents an algebraic curve). In an
exact arithmetic context, we can substitute R0 in fz; f(x;y) = 0g and verify if the results are
zero. Similarly, we could handle:
R0 \ff(R) = g(R) = 0g
where f(R) = g(R) = 0 represents an implicit 3D space curve.
What does verify mean in \ oating point" arithmetic?
Example A:
Let z0 = 0 and x0, y0 satisfy
jf(x0;y0)j< 1 (10.1)
where is a small constant and jf(x;y)j 1 in the domain of interest including (x0;y0),
then a \distance" check can be performed by:
jf(x0;y0)j
j5f(x0;y0)j < 1 (10.2)
provided j5f(x0;y0)j 6= 0. Equation 10.1 is called the \algebraic distance" and Equa-
tion 10.2 is called the \non-algebraic distance". The true geometric distance is given
by:
d = minjR R0j; where R = (x;y); f(R) = 0 (10.3)
The true geometric distance is di cult and expensive to compute (particularly for implicit
f(R) = 0 and involves computing the global minimum of jR R0j. Equation 10.2 results
from the rst order approximation of Equation 10.3 as derived by Taylor expansion and
is exact when f(R) is represents a plane.
8
.R0 g(R) = 0
f(R) = 0
Figure 10.5: Curves meet at small angle.
Example B:
R0 \ff(R) = g(R) = 0g
When curves f = 0, g = 0 meet at a small angle ( 5fj5fj 5gj5gj = 1), then the condition
jfj< and 1 = jfjj5fj <
jgj< and 2 = jgjj5gj <
(where jfj;jgj; 1; 2 are evaluated with R = R0 and ; 1) are not enough to guarantee
proximity of R0 to the intersection of f, g, see Figure 10.5.
δ
δ
δ
1
2
3
φ
g
f
Figure 10.6: Approximate curves with straight lines.
Using a linear approximation, and letting
= cos 1 j 5fj5fj 5gj5gjj
be the angle of intersection as in Figure 10.6 near the intersection point, a better criterion
for evaluating if R0 is near the intersection of f and g is
3 = 1f jfjj5fj + jgjj5gjg< 1
9
10.4.2 Point/Parametric curve intersection
1. Rational polynomial curves
R0 \R = R(t) A t B
Brute force elementary method:
We solve each of the following three nonlinear polynomial equations separately and
we search for common real roots in A t B.
x(t) x0 = 0 ! t01; ;t0n
y(t) y0 = 0 ! t001; ;t00n
z(t) z0 = 0 ! t0001 ; ;t000n
In principle, this elementary approach is \easy" for polynomials. However, in prac-
tice, this process is complex and ine cient and prone to numerical inaccuracies.
Preprocessing and subdivision method
{ Use bounding box of R(t) to eliminate easily resolvable cases, with some level
of subdivision (splitting) to reduce box size.
{ Concept of subdivision in rational arithmetic: To eliminate numerical error
in the subdivision process (which can lead to erroneous decisions), rational
arithmetic may be employed (if the input coe cients of R(t) are rational or
oating point numbers). This can be easily done in object-oriented languages
such as C++ using operator overloading.
{ Continue subdivision until box is small.
{ Then, we could use a numerical technique , such as:
F(t) = minfjR0 R(t)j2g t2D1 [A;B]
and use some t from the interval D1 as the initial approximation. Use of the
square of the distance function is necessary to avoid possible divergence of the
derivative of the distance function, if it approaches zero.
{ If the minimization process converges to t0 and
q
F(t0) < , t = t0 is the desired
solution.
Implicitization (perhaps with box preprocessing) such as
(x0;y0) \ fx = x(t); y = y(t)g
Let us consider an example where x(t);y(t) are quadratic polynomials (the curve is
a parabola). We will attempt to eliminate t to get a polynomial f(x;y) = 0 which
the x;y coordinates of all points on the curve satisfy. We start with the system
x = a0t2 +b0t+c0
y = a00t2 +b00t+c00 )
a0t2 +b0t+c0 x = 0
a00t2 +b00t+c00 y = 0
10
which can be rewritten in matrix form as follows:
)
2
66
64
c0 x b0 a0 0
0 c0 x b0 a0
c00 y b00 a00 0
0 c00 y b00 a00
3
77
75
2
66
64
1
t
t2
t3
3
77
75 =
2
66
64
0
0
0
0
3
77
75 (10.4)
The maximum degree of t in the above vector is determined by the degree m of the
x polynomial and the degree n of the y polynomial, and is given by m+n 1. In
this case m+n 1 = 3.
A necessary and su cient condition for the above system to be solvable is
c0 x b0 a0 0
0 c0 x b0 a0
c00 y b00 a00 0
0 c00 y b00 a00
= f(x;y) = 0
The equation f(x;y) = 0 is the implicit equation of the curve. Consequently in an
exact arithmetic context, we need to check if f(x0;y0) = 0, to verify if (x0;y0) is on
the initial curve.
In general, if
x = x(tn); y = y(tn) )F(xn;yn) = 0
where n is the total degree.
Inversion: If f(x;y) = 0 then we could use the rst 3 equations 10.4:
2
64 b0 a0 0c
0 x0 b0 a0
b00 a00 0
3
75
2
64 tt2
t3
3
75 =
2
64 c0 x00
c00 y0
3
75
)t = (x0;y0) (x
0;y0)
Where and are polynomials in x0 and y0, and x0, y0 satisfy
f(x0;y0) = 0
{ The method is e cient and (usually) accurate for n 3 (but no real guaran-
tees on accuracy and robustness exist if the method is implemented in oating
point).
{ Subdivision methods are preferable for higher n, and as we will see later when
coupled with rounded interval arithmetic are robust, accurate and e cient.
Intersection of points (x0;y0;z0) and 3D polynomial curves R = R(t) via implicit-
ization of such curves involves a process of projection on x;y plane and nding t0
by inversion and veri cation of z0 = z(t0).
11
10.4.3 Point/Procedural parametric (o set, evolute, etc.) curve
intersection
R0 \R = R(t) A t B
In general there is no known and easily computable convex box decreasing arbitrarily
with subdivision!
An approximate solution method may involve minimization of
F(t) = jR(t) Rj2
where t2 [A;B]. This would involve
{ Checking end points, ie. if F(A), F(B) are very small.
{ Initial estimate for the possible minima, perhaps using linear approximation of R(t)
to start the process.
However,
{ Convergence of the above minimization processes is not guaranteed in general.
{ There may exist more than one minima.
{ Convergence to local and not global minimum (where F(t) 6= 0 ) is possible.
For certain classes of procedural curves such as o sets and evolutes of rational curves
involving radicals of polynomials, it is possible to use the \auxiliary variable method" to reduce
the point to curve intersection (or minimum distance) problem to a set of (a larger number
of) nonlinear polynomial equations. Such systems can be solved robustly and e ciently using
the nonlinear solver describe in the next section.
12
10.5 Point/surface intersection
10.5.1 Point/Implicit (usually algebraic) surface intersection
The condition for R0 \ff(R) = 0g, where f(R) = 0 is an implicit surface, is:
jf(R0)j< ; jf(R0)jj5f(R
0)j
<
where ; are small constants.
10.5.2 Point/Rational polynomial surface intersection
1. Implicitization is possible for all such surfaces but computationally expensive and possi-
bly inaccurate. For a tensor product rational polynomial surface with maximum degrees
in u and v equal to m and n, of the form
R = R(um;vn);
the implicit equation is
f(xq;yq;zq) = 0
where q 2mn
Therefore;for m = n = 3 ! q 18; m = n = 2 ! q 8
The above method is useful for special surfaces such as cylindrical and conical ruled
surfaces, surfaces of revolution, etc.
Examples:
(a) Implicitization of a surface of revolution.
y
x
z
r
r
R(t)
Figure 10.7: Surface of revolution.
13
Let us consider a pro le curve to be a rational polynomial of degree n, see Figure 10.7
R(t) = [r(t);z(t)]
By simple implicitization of R = R(t), we get:
fn(r;z) = 0 (10.5)
where n is the maximum total degree of f. Also,
r2 = x2 +y2 (10.6)
Next we eliminate r from equations 10.5 - 10.6 by rewriting equations 10.5 - 10.6 as
follows:
fn(r;z) = a0(z)rn +a1(z)rn 1 + +an(z) = 0
) r2 + (x2 +y2) = 0
The resultant of eliminating r from these two equations is
D =
a0(z) a1(z) an 1(z) an(z) 0
0 a0(z) an 2(z) an 1(z) an(z)
1 0 x2 +y2
1 0 x2 +y2
...
1 0 x2 +y2
= 0
and the degree of D f(x;y;z) = 0 is 2n. An example is a torus (degree 4 algebraic
surface).
(b) Implicitization of a cylindrical ruled surface
x
y
^n
R(t)
^t
a
z
Figure 10.8: Cylindrical ruled surface.
14
Let
R(t) = [X(t);Y(t)]
be a degree n planar rational polynomial curve in the x, y plane. The resulting
implicit equation of the curve
f(X;Y) = 0
is a polynomial of degree n. Let
a = [a1; a2; a3]
be a direction vector. Then the three equations
x = X +ua1
y = Y +ua2
z = ua3
describe a cylindrical ruled surface. Hence, the implicit surface equation becomes:
f(x za
3
a1;y za
3
a2) = 0
This equation can be transformed to the standard form using a symbolic manipu-
lation program such as Macsyma.
(c) Implicitization of a conical ruled surface
x
y
z
R(t)
u
R0
^t
Figure 10.9: Conical ruled surface.
Let
R0 = [x0; y0; z0]
be the apex of the conical ruled surface and
R(t) = [X(t); Y(t)]
15
a degree n planar rational polynomial curve on the x;y plane. Its implicit equation
f(x;y) = 0
is a degree n polynomial. The equation of the resulting conical ruled surface is
x = x0(1 u) +Xu
y = y0(1 u) +Yu
z = z0(1 u)
Eliminating u = 1 zz0 and solving for X; Y yields:
f
z
0
z0 zx
x0
z0 zz;
z0
z0 zy
y0
z0 zz
= 0
This equation can be transformed to the standard form using a symbolic manipu-
lation program such as Macsyma.
2. Newton’s method:
Solve x0 = x0(u;v), y0 = y0(u;v) and verify the third equation. Use a linear approxi-
mation to start the process. Preprocessing using convex bounding box should always be
used, coupled with some level of subdivision.
3. Convex box and possibly subdivision followed by minimization in 0 u;v 1 or
within a rectangular subdomain of the following function (see Figure 10.10).
F(u;v) = jR(u;v) R0j2 0
A point u0; v0 where F(u0;v0) = 0 yields the solution.
u
v
F(u,v)
O
Figure 10.10: Distance function squared
In order to solve this minimization problem, we need to compute
Minimum of all local minima of F(u;v) (Fu = Fv = 0) in domain;
16
Minimum of all local minima of boundary \curves" eg. F(0;v) (i.e. Fv = 0);
Values of F(u;v) at corners, ie. F(0;0), F(0;1), F(1;0), F(1;1);
and then choose u0;v0 from above solutions where F(u0;v0) = 0.
The disadvantages of a minimization method are:
(a) Initial approximation is required;
(b) Possibility of divergence;
(c) No guarantee that all minima are located. (We need to enhance con dence by
subdivision.)
(d) First and second derivativers of F(u;v) are required.
Note: When R = R(u;v) is a polynomial parametric surface patch, it is helpful to
reformulate F(u;v) to
F(u;v) =
kX
i=0
mX
j=0
wijBi;k(u)Bj;m(n) (10.7)
If wij > 0 for all i; j then there is no solution. We could use wij to construct initial
approximations for the various local minima to be computed by usual descent numerical
methods. These initial approximation may be obtained by discrete sampling or subdivi-
sion.
Let F(u;v) be expressed in the Bernstein basis, as in equation (10.7). Then, let us also
express Fu; Fv in the Bernstein basis:
Fu(u;v) =
k 1X
i=0
mX
j=0
AijBi;k 1(u)Bj;m(v) = 0
Fv(u;v) =
kX
i=0
m 1X
j=0
BijBi;k(u)Bj;m 1(v) = 0 (10.8)
The equations Fu(u;v) = 0 and Fv(u;v) = 0 represent planar algebraic curves illustrated
in Figure 10.11. Their intersection are the required extrema from which the minima can be
selected using elementary calculus.
A geometrically motivated solution of the system 10.8 is possible using the convex hull
property and subdivision to isolate an area where convex hulls intersect.
Taking G(u;v) = Fu(u;v) and n = k 1 for example, we can write
w = G(u;v) =
nX
i=0
mX
j=0
AijBi;n(u)Bj;m(v)
We can reformulate this \height" function into a parametric surface as follows:
w = [u;v;G] = XXLijBi;m(u)Bj;n(v)
Lij = [ im; in; Aij]
17
v
Fv = 0
Fu = 0
u 1
1
0
0
Figure 10.11: Intersection of algebraic curves.
w=0 curve
w=1?u ?v2 2 = 0
w
u
v
1
Figure 10.12: Control net
To solve Fu = Fv = 0, we nd the projection of convex hulls of w(u;v) and of the corre-
sponding surface for Fv on the coordinate planes u = 0, v = 0 , then intersect them with the
lines w = 0, see Figure 10.12 for an example. A more detailed description of this procedure is
given in the next section where a nonlinear solver is described.
18
10.5.3 Point/Procedural surface intersection
A procedural surface may be an o set surface or a generalized cylinder surface or a blending
surface. The typical solution method is minimization. In this case, no convex box assistance
is possible in general, and we need a dense sampling for an initial approximation (which may
be expensive) and no rigorous guarantees for the solution’s reliability are generally available.
For certain classes of procedural surfaces such as o sets and evolutes of rational surfaces
involving radicals of polynomials, it is possible to use the \auxiliary variable method" to reduce
the point to surface intersection (or minimum distance) problem to a set of (a larger number
of) nonlinear polynomial equations. Such systems can be solved robustly and e ciently using
the nonlinear solver describe in the next section.
19
10.6 Curve/curve intersection
The curve types we will consider can be classi ed as follows:
Rational polynomial parametric (RPP)
Procedural parametric (PP)
Implicit algebraic (IA)
Implicit procedural (IP)
Procedural curves may be general o sets, evolutes, etc.
Curve to curve intersection cases are identi ed in table 10.1.
RPP PP IA IP
RPP D1 D2 D3 D4
PP D5 D6 D7
IA D8 D9
IP D10
Table 10.1: Curve to curve intersection cases
Conceptually, D3 (RPP/IA curve intersection) is the \simplest" of the above cases of inter-
section to describe and use for illustrating various general di culties of intersection problems.
10.6.1 Case D3: RPP/IA curve intersection
We start with a planar RPP curve (which is conceptually the simplest case).
R(t) = [x(t);y(t)] = [X(t)w(t);Y(t)w(t)]
=
Pn
i=0 RihiBi;n(t)P
n
i=0hiBi;n(t)
0 t 1
where hi > 0 are weights and Bi;n(t) are Bernstein basis functions.
The implicit algebraic curve fm(x;y) = 0 of total degree m is described by
fm(x;y) =
mX
i=0
m iX
j=0
cijxiyj = 0
For convenience, we convert it to homogeneous form, by setting x = Xw , y = Yw and multiplying
by wm, which leads to
fm(X;Y;w) =
mX
i=0
m iX
j=0
cijXiYjwm i j = 0
20
Every term of the above sum has a total degree m. Substituting R(t) into fm(X;Y;w) leads
to a polynomial of degree up to mn
F(t) = 0
Therefore, now the problem of intersection is equivalent to nding the real roots of F(t) in
0 t 1. The most usual form of F(t) is the power basis. The coe cients ai can be evaluated
symbolically by substitution and collection of terms. This can be readily done in a standard
symbolic manipulation program (such as Macsyma, Reduce, Maple etc.). Such programs are
oriented to processing rational numbers exactly.
Example:
Let the algebraic curve be an ellipse x24 +y2 1 = 0, as illustrated in Figure 10.13. Let the
parametric curve be a cubic B ezier curve with control points:
" 0
1
#
;
" 1
4
#
;
" 2
1
#
; and
" 2
0
#
x
y
1
2
Figure 10.13: Ellipse and cubic B ezier curve polygon
Using Macsyma (or any similar symbolic manipulation program) and simplifying, we get
(in exact arithmetic mode):
F(t) = 1025t6 3840t5 + 5514t4 3728t3 + 1149t2 120t = 0
Next we nd the (real) roots of F(t) in t 2 [0;1] using Macsyma’s factoring capability over
integers, which leads to
F(t) t(t 1)2G(t)
where
G(t) = [1025t3 1790t2 + 909t 120]
21
Using a standard numerical solver for polynomials in oating point (such as NAG C02AEF
routine), we obtain the following numbers as solutions of G(t) = 0 (reported with four decimal
digits)
t = 0:9228; 0:61843; 0:2051
Alternately solving F(t) = 0 using the same routine leads to the following roots t = tR +itI
tR 1 1 0:9228 0:6183 0:2051 0
tI 0:22 10 6 0:22 10 6 0 0 0 0
Notice the sensitivity to errors for the 6th degree polynomial, especially for multiple roots
as t = 1. In oating point arithmetic, such roots split into a number of roots (complex or real).
Obviously, complex roots are not usable (as we require only the real intersection points). The
consequence is lost roots, which implies an erroneous solution of the intersection problem.
An alternate basis for the representation of F(t) = 0 is the Bernstein basis, which has better
stability of its real roots under perturbations of its coe cients than the power form. We will
introduce the concept of condition numbers for polynomial roots later in this subsection.Such
representation can be obtained without rst converting to a power basis and using a symbolic
manipulation program. It rather requires polynomial arithmetic involving products such as
Bi;k(t)Bj;l(t) =
k
i
! l
j
! k +l
i+j
! 1
Bi+j;k+l(t)
In this way we obtain the representation:
F(t) =
mnX
i=0
ciBi;mn(t) = 0
In the above example of the ellipse and cubic B ezier curve,
c0 = 0; c1 = 20; c2 = 36:6; c3 = 16:6; c4 = 1:6; c5 = 0; andc6 = 0
Using the linear precision property
t =
mnX
i=0
i
mnBi;mn(t)
we can construct
f(t) = [t;F(t)] =
mnX
i=0
i
mnc
i
!
Bi;mn(t)
which is now a standard degree mn B ezier curve as in Figure 10.14
Notice that in our example c0 = 0 which implies that t = 0 is a root. Also c5 = c6 = 0
implies that t = 1 is a double root. Note that
If all ci > 0 or ci < 0 then there are not roots in [0;1].
We can use subdivision (split in half) to identify subintervals of [0;1] where coe cients
change sign once. Variation diminishing property implies one root in such areas. The
Newton method can be used for fast convergence within such subintervals.
22
15
30
45
?15
?30
?45
1/6
2/6 3/6 4/6 5/6 1 t0
0.21
0.61
0.92
Figure 10.14: Intersection of a B ezier curve/straight line.
Another solution method is illustrated in Figure 10.15 for a quadratic curve (parabola). In this
method the curve’s convex hull is intersected with the axis w = 0 to give the [A;B] interval
as in Figure 10.15. The part of the [0;1] interval outside [A;B] does not contain roots. By
curve subdivision at A and B a smaller curve segment can be obtained in the B ezier form and
the process can be continued. In this case, we use binary subdivision of AB when the rate of
decrease of AB slows. See paper by Sherbrooke and Patrikalakis for details and also the next
section on the nonlinear solver.
0 1A B
w
w = 0
Figure 10.15: Subdivision at A; B.
Numerical condition of polynomials in Bernstein form
1. Condition numbers for polynomial roots
Consider the displacement x of a real root xo of a polynomial in the basis f k(x)g:
P(x) =
nX
k=0
k k(x);
23
due to a perturbation r in a single coe cient r in this basis. Since xo + x is a root
of the perturbed polynomial, it satis es:
P(xo + x) = r r(xo + x):
Performing a Taylor series expansion about xo on both sides of the above equation and
noting that P(xo) = 0, we obtain:
nX
k=1
( x)k
k! P
(k)(xo) = r
nX
k=0
( x)k
k!
(k)
r (xo):
If xo is a simple root of P(x), then P0(xo) 6= 0, and in the limit of in nitesimal pertur-
bations the above equation gives:
lim
r!0
x
r= r =
r r(xo)
P0(xo) :
We de ne
C = j r r(xo)=P0(xo)j
the condition number of the root xo with respect to the single coe cient r.
If xo is an m-fold root, m 2, then we de ne a multiple-root condition number C(m) in
the form
C(m) = [ m!jP(m)(x
o)j
nX
r=0
j r r(xo)j]1=m:
Theorem For an arbitrary polynomial P(x) with a simple root xo 2 [0;1], let Cp(xo) and
Cb(xo) denote the condition numbers of the root in the power and Bernstein bases on [0;1],
respectively. Then Cb(xo) Cp(xo) for all xo 2 [0;1]. In particular Cb(0) = Cp(0) = 0,
while for xo 2 (0;1] we have the strict inequality Cb(xo) <Cp(xo).
2. Example { Wilkinson polynomial
Consider the polynomial with the linear distribution of real rootsxo = k=n;k = 1;2;:::;n,
on the unit interval [0, 1] for n = 20:
P(x) =
20Y
k=1
(x k=20):
The condition numbers for each root with respect to a perturbation in the single co-
e cient a19 are shown in Table 10.2. It is evident from Table 10.2 that the Bernstein
form a ords a dramatic inprovement in root condition numbers compared to the power
form in this example, the condition number of the most unstable root being reduced by
a factor of about 107.
We now perturb the single power coe cient a19 = 212 by an amount 2 23=20. This
corresponds to a fractional perturbation = 2 23=210 = 5:7 10 10. The roots of the
24
Table 10.2: Condition numbers for Wilkinson polynomial
k Cp(xo) Cb(xo)
1 2:100 101 3:413 100
2 4:389 103 1:453 102
3 3:028 105 2:335 103
4 1:030 107 2:030 104
5 2:059 108 1:111 105
6 2:667 109 4:153 105
7 2:409 1010 1:115 106
8 1:566 1011 2:215 106
9 7:570 1011 3:321 106
10 2:775 1012 3:797 106
11 7:822 1012 3:321 106
12 1:707 1013 2:215 106
13 2:888 1013 1:115 106
14 3:777 1013 4:153 105
15 3:777 1013 1:111 105
16 2:833 1013 2:030 104
17 1:541 1013 2:335 103
18 5:742 1012 1:453 102
19 1:310 1012 3:413 100
20 1:378 1011 0
25
pertured polynomials are shown in Table 10.3. For comparison, we now perturb the
coe cient
c19 = 1484925542112800000000000000000
of the Bernstein form by the same fractional amount, and obtain approximations to the
roots of this perturbed polynomial, which are shown in Table 10.3.
Table 10.3: Roots of perturbed polynomials
k power Bernstein
1 0.05000000 0.0500000000
2 0.10000000 0.1000000000
3 0.15000000 0.1500000000
4 0.20000000 0.2000000000
5 0.25000000 0.2500000000
6 0.30000035 0.3000000000
7 0.34998486 0.3500000000
8 0.40036338 0.4000000000
9 0.44586251 0.4500000000
10 0:50476331 0.5000000000
11 0:03217504i 0.5499999997
12 0:58968169 0.6000000010
13 0:08261649i 0.6499999972
14 0:69961791 0.7000000053
15 0:12594150i 0.7499999930
16 0:83653687 0.8000000063
17 0:14063124i 0.8499999962
18 0:97512197 0.9000000013
19 0:09701652i 0.9499999998
20 1.04234541 1.0000000000
3D-Space geometry: Case D3: RPP/IA (continued)
R(t) = [x(t);y(t);z(t)] \ f(R) = g(R) = 0
1. Substitute
f(R(t)) F1(t) = 0
g(R(t)) F2(t) = 0
2. Compute the resultant of F1(t);F2(t) by eliminating t.
3. If R(F1(t);F2(t)) 0, then there is a common root between the two equations.
4. Use the inversion algorithm to nd t.
26
10.6.2 Case D1: RPP/RPP Curve Intersection
( R
1 = R1(t); 0 t 1
R2 = R2(v); 0 v 1
Setting R1(t) = R2(v) leads to 3 nonlinear equations with 2 unknowns (overdetermined
system). A possible approach is to choose 2 equations to solve for t;v, and then substitute the
results into the third equation for veri cation. Alternatively the subdivision-based nonlinear
solver described in the next section can directly solve such systems without the above 2 step
approach.
Preprocessing idea in 3D:
Check bounding boxes for intersection. If there is such intersection, examine x;y projection.
Method 1 .
1. Implicitize, e.g., x2 = x2(v), y2 = y2(v), to get f(x;y) = 0.
2. Substitute x1(t);y1(t) into f(x;y) to get F(t) = 0 and solve it for real roots in [0;1].
3. Use the inversion algorithm: v = (x;y)q(x;y)
Method 2 Recursive subdivision and use of bounding boxes.
Some issues:
Method 1 is more e cient for n 3 but its robustness is unclear in the presence of ill-
conditioned (tangent) intersections.
Method 2 is more e cient for n 4 and its robustness with respect to missing roots can
be guaranteed if the method is inplemented in rounded interval arithmetic (see next
section). If two bounding boxes intersect, and they are of nite size, we can nd roots
using linear approximation.
In Figure 10.16(a), the boxes intersect, the linear approximation do not, and the curves
intersect. Similar behavior is observed in (b) where polygon is used as the curve approx-
imation.
(a) (b)
Figure 10.16: Ill-conditioned curve intersections.
27
Use hodograph (as in Figure 10.17) to nd the range of tangent variation. Construct
bounding angular sectors of hodographs of two curves and make vertices coincident. If
the sectors do not intersect, then there is at most one root; otherwise, subdivide the two
curves. For a precisely \tangent" root, this method would lead to in nite subdivision
steps.
(1)
(2) (1)
(2)
angular sector
Figure 10.17: Hodograph concept.
10.6.3 Case D2/D5: RPP/PP and PP/PP Curve Intersections
3 equations with 2 unknowns
( R
1 = R1(t); 0 t 1
R2 = R2(v); 0 v 1
Possible Approach
Minimize
F(t;v) = jR(t) R(v)j2; 0 t;v 1
See comments on Point/RPP surface intersection.
More di cult to compute derivatives of F(t;v) exactly. May need numerical techniques
(slower and usually more inaccurate).
10.6.4 Case D6: PP/IA Curve Intersection
( R
1 = R1(t); 0 t 1
g(R) = f(R) = 0
It could be reduced to PP Curve/IA Surface intersection and comparison of solutions for g = 0
and f = 0.
28
10.6.5 Case D8: IA/IA Curve Intersection
The planar case is of interest in processing trimmed patches.
f(u;v) = 0
g(u;v) = 0
)
(u;v) 2 Parametric domain
Method 1
Eliminate v to form the resultant F(u), then solve F(u) = 0 for u and use the inversion
algorithm to get v.
Example: Let us consider an ellipse and a circle
f = x
2
4 +y
2 1 = 0
g = (x 1)2 +y2 1 = 0
as in Figure 10.18.
x
y
Figure 10.18: Ellipse and circle intersection
Let us eliminate y from these two equations. This leads to
3x2 8x+ 4 = 0
which has as rootsx = 2 andx = 23. Each of these leads toy2 = 0 andy2 = 89 respectively.
However there are possible numerical problems at the \tangent" solution x = 2; y = 0.
Let us assume that due to error
x = 2 +"
hence
y2 = "(1 + "4) < 0
This implies that y is imaginary and that no real roots exist. This would have as a
consequence missing an intersection solution, leading to a robustness problem.
Method 2
After tracing of f(u;v) = 0 and g(u;v) = 0, linear approximation is available. Find
intersections of linear approximations and minimum distance points between them. Use
this to drive a Newton Method on f = g = 0 or a minimization of F = f2 +g2.
29
10.7 Curve/surface intersection
Such intersections are classi ed as in Table 10.4. We will start with case E3 which is quite
representative.
Surface Type
Curve Type RPP PP IA IP
RPP E1 E2 E3 E4
PP E5 E6 E7 E8
IA E9 E10 E11 E12
IP E13 E14 E15 E16
Table 10.4: Curve/Surface Intersections
Curve to surface intersection involves the intersection of straight line to surface of all cases.
Such intersection is useful for
ray tracing
point classi cation
for procedural surface interrogation
10.7.1 Case E3: RPP Curve/IA Surface Intersection
Let us consider an (implicit) algebraic surface of total degree m
fm(x;y;z) =
mX
i=0
m iX
j=0
m i jX
k=0
cijkxiyjzk = 0
We rst convert it to a degree homogeneous form by setting x = X=W;y = Y=W;z = Z=W
and multiply by Wm leading to
fm(X;Y;Z;W) =
mX
i;j;k;q=0
i+j+k+q=m
cijkXiYjZkWq = 0
Let us also consider a rational curve of degree n
R(t) = [x(t);y(t);z(t)] = [X(t)W(t); Y(t)W(t); Z(t)W(t)]; 0 t 1
We can easily substitute R(t) in fm(X;Y;Z;W) to obtain a polynomial equation F(t) = 0 of
degree mn. We then nd its (real) roots in [0;1], see comments under section 10.6.1.
30
10.7.2 Case E1: RPP Curve/RPP Surface Intersection
3 nonlinear equations in 3 unknowns t;u;v;R(t) = Q(u;v) where
R = R(t); 0 t 1
Q = Q(u;v); 0 u;v 1
A Preprocessing step is to check bounding boxes for absence of intersection to eliminate
easy cases.
Method 1 Implicitization of Q(u;v) for simple surface forms and reduction to Case E3.
Method 2 Recursive subdivision and use of bounding boxes. See also next section for nonlin-
ear solver. Eventual use of a linear approximation technique for R and Q to obtain ap-
proximate solution, which is to be used to initiate a Newton method on R(t) Q(u;v) = 0
or a minimization method on F(t;u;v) = jR Qj2. See section 10.6.2 for problem areas
for Ft = Fu = Fv = 0 (3 equations.)
10.7.3 Case E2/E6: RPP/PP, PP/PP Curve/Surface Intersection
3 nonlinear equations in 3 unknowns t;u;v;R(t) = Q(u;v) where
( R = R(t); 0 t 1 Q = Q(u;v); 0 u;v 1
Possible Approach
Minimize
F(t;u;v) = jR(t) Q(u;v)j2
in a cube 0 t;u;v 1. See comments under point-surface intersection (e.g. examine vertices,
edges, etc.)
10.7.4 Case E7: PP Curve/IA Surface Intersection
4 non-linear equations in 4 unknowns t;R
( R = R(t); 0 t 1
f(R) = 0
Possible approach
Could use a Newton method initiated by a linear approximation of R = R(t), which can be
intersected more easily with f(R) = 0 using the method of Case E3.
Problems
1. All roots?
2. Convergence?
31
3. E ciency?
In case of convergence problems, embedding or continuation methods are normally helpful.
For example,
1. Let
Q(t) = Q0 + (Q1 Q0)t ab a; t2 [a;b]
where [a;b] [0;1], be an approximation of R(t) within [a;b] interval.
2. Compute Q(t) \f(R) = 0 in t2 [a;b] using the method of Case E3.
3. De ne sequence of problems
R(t;") = Q(t) +"(R(t) Q(t))
where " = "1;"2; ;"N such that 0 = "1 <"2 < <"N = 1.
4. Solve R(t;"i)\f(R) = 0 using as initial approximation the solution for " = "i 1.
Such a method does not by itself provide initial approximations of all possible solutions; rather,
it assists in the computation of a particular solution.
10.7.5 Case E11: IA Curve/IA Surface Intersection
3 equations in 3 unknowns
f(R) = g(R)| {z }
curve
= h(R)| {z }
surface
= 0
Possible approaches
1. Elimination Methods
2. Newton Methods
3. Minimization Methods
F(R) = f2 +g2 +h2
Possible reformulation in a box
F(R) =
MX
i=0
NX
j=0
QX
k=0
wijkBi;M( x)Bj;N( y)Bk;Q( z)
(x;y;z) 2 = [a1;a2] [b1;b2] [c1;c2]
and x = (x a1)=(a2 a1) and similarly for y and z. If wijk > 0 for all i;j;k, there is no
solution is .
Initial approximation: nd mini;j;kwijk < 0 and start at ( iM; jN; kQ).
4. Approximate f(R) = g(R) = 0 curve with a linear spline, reduce to E3 and re ne using
minimization.
32
10.7.6 IA Curve/RPP Surface Intersection
( curve f(R) = g(R) = 0
surface R = R(u;v); 0 u;v 1
Substitute R = R(u;v) in f(R) = 0; g(R) = 0 to obtain two algebraic curves f(u;v) =
0;g(u;v) = 0, as in Figure 10.19. This formulation reduces to Case D8 in Section 10.6.5:
IA/IA curve intersection. Algebraic curves are treated under intersections of algebraic and
RPP surfaces.
f=0 g=0
u
v
Figure 10.19: Intersection of two algebraic curves
33
10.8 Surface/Surface Intersections
The surface types we will consider can be classi ed as follows:
Rational polynomial parametric (RPP)
Procedural parametric (PP)
Implicit algebraic (IA)
Implicit procedural (IP)
Procedural surfaces may include general o sets, focal surfaces, etc.
Surface to surface intersection cases are identi ed in table 10.5.
Surface Type
Surface Type RPP PP IA IP
RPP F1 F2 F3 F4
PP F5 F6 F7
IA F8 F9
IP F10
Table 10.5: Classi cation of Surface/Surface Intersections
The solution of a surface/surface intersection problem may be empty, or include a curve
(possibly made of several branches), a surface patch, or a point. Conceptually, F3 (RPP/IA
surface intersection) is the \simplest" of the above cases of intersection to describe and use for
illustrating general di culties of surface intersection problems.
10.8.1 Case F3: RPP/IA Surface Intersection
We start with a RPP surface patch
R(u;v) = [X(u;v)W(u;v); Y(u;v)W(u;v); Z(u;v)W(u;v)] (10.9)
where X;Y;Z are all of degree p in u, q in v, (u;v) 2 [0;1] [0;1].
Next we consider an implicit algebraic surface fm(x;y;z) = 0 of total degree m described
by
fm(x;y;z) =
mX
i=0
m iX
j=0
m i jX
k=0
cijkxiyjzk (10.10)
Examples of such surfaces in practical use are low order surfaces such as planes (degree 1), the
natural quadrics (cylinder, sphere, cone) (degree 2), and torii (degree 4). In fact in a survey of
mechanical parts (mechanical elements), over 90% of all surfaces involved are of these types. It
is also well known that all these surfaces also have a low degree rational polynomial parametric
34
representation, so that when two surfaces of the above types are intersected, the methods of
this section may be used.
For convenience, we convert fm(x;y;z) = 0 to its homogeneous form by setting
x = XW; y = YW; z = ZW (10.11)
and multiplying by Wm, which leads to
fm(X;Y;Z;W) = X
i;j;k;q 0
i+j +k +q = m
cijkXiYjZkWq = 0 (10.12)
Consequently, the intersection problem
fm(X;Y;Z;W) = 0 (10.13)
X = X(u;v); Y = Y(u;v); Z = Z(u;v); W = W(u;v) (10.14)
may be thought of as a nonlinear system of 5 equations in 6 unknowns. A reduction of the di-
mensionality of the system may be obtained by substituting equations 10.14 in equation 10.13,
which since all functions involved are polynomial leads to an algebraic curve
f(u;v) = 0 (10.15)
of degree M = mp and N = mq in u;v, respectively. Consequently, the problem of intersection
reduces to the problem of tracing f(u;v) = 0 without omitting any special features of the curve,
e.g., small loops, singularities, and accurately computing all its branches. This is a fundamental
problem in algebraic geometry and much work has been done to understand its solution. In the
context of algebraic geometry the coe cients of f(u;v) = 0 are integers. In the context of CAD
and computer implementation, the coe cients of fm = 0, and R = R(u;v) are oating point
numbers. Consequently, if the above substitution is performed in oating point arithmetic the
coe cients of f(u;v) = 0 involve error. To avoid such error, rational arithmetic may be used
for robustness. These issues will also be discussed in the next section.
The algebraic curve
f(u;v) =
MX
i=0
NX
j=0
aijuivj = 0 (10.16)
can be reformulated as
f(u;v) =
MX
i=0
NX
j=0
wijBi;M(u)Bj;N(v) = 0 (10.17)
where (u;v) 2 [0;1]2.
As an example consider a plane in homogeneous form
AX +BY +CZ +DW = 0 (10.18)
35
u
v
1
1
border
points singular points
loop
Figure 10.20: Parameter space of R(u;v) and resulting algebraic curve f(u;v) = 0
and a rational B ezier patch of degree p in u, q in v
R(u;v) =
Pp
i=0
Pq
j=0hijRijBi;p(u)Bj;q(v)P
p
i=0
Pq
j=0hijBi;p(u)Bj;q(v)
(10.19)
where Rij = [xij;yij;zij] and weights hij 0.
The resulting algebraic curve is of the form of equation 10.17 with
wij = (Axij +Byij +Czij +D)hij (10.20)
In fact the power basis form off(u;v) = 0 need not be computed at all, if polynomial arithmetic
for Bernstein polynomials is used.
The advantage of the Bernstein form is its numerical stability and convex hull property.
If wij > 0 or < 0 for all i;j, there is no solution and the two surfaces do not intersect.
More precisely, the (entire) algebraic surface fm(R) = 0 does not intersect the surface patch
R = R(u;v) for (u;v) 2 [0;1]2.
What happens when all wij = 0 (or ij = 0)? Obviously, the two surfaces coincide in
their entirety. A somewhat complex algebraic curve f(u;v) = 0 is shown in Figure 10.20
involving various branches (from border to border), internal loops, and singularities, see also
Figure 10.24.
Given a point on every branch (connected component) of an algebraic curve, tracing the
curve using di erential curve properties is an e ective procedure (marching method).
Marching Method:
Let us expand f(u;v) as a Taylor series
36
δvL
δv
δu
(u,v)
Newton
Figure 10.21: A zoomed view of an algebraic curve near a point (u;v)
f(u+ u;v + v) = f(u;v) +fu u+fv v
+12(fuu u2 + 2fuv u v +fvv v2) + (10.21)
- To rst order and if f2u +f2v > 0, in order to have f(u;v) = 0 and f(u+ u;v + v) = 0,
fu u+fv vL = 0 ) vL = fuf
v
u (10.22)
see also Figure 10.21.
- The Newton method on f(u + u;v) = 0 with initial approximation vI = v + vL may be
used to compute vF = v + v with high accuracy and in an e cient manner.
- For \vertical" branches, ie. when jfvj is very small, we may use uL = fvfu v.
To avoid these special stepping procedures the equation fu u+fv v = 0 may be converted
to
fu _u+fv _v = 0 (10.23)
where u;v are considered functions of a parameter t, u = u(t); v = v(t). This equation is
satis ed if
_u = fv(u;v) (10.24)
_v = fu(u;v) (10.25)
where is an arbitrary constant. For example, can be chosen to be equal to [f2u +f2v ] 12 , in
order that t is an arc length parameter in the [u;v] 2 [0;1]2 parameter space. This is a system
37
of two rst order nonlinear di erential equations which can be solved by the Runge-Kutta or
other methods with adaptive step size.
Problems of Marching Methods
1. Starting points on all branches need to be provided in advance.
2. Marching through singularities (f2u +f2v ’ 0) is problematic.
3. Step size selection is complex and too large a step size may lead to straying or looping,
as in Figure 10.22.
straying looping
Figure 10.22: Step size problems in marching method
Computation of starting points
Starting points for tracing algebraic curves are certain \characteristic" points de ned below:
1. Border points, ie. the intersections of f(u;v) = 0 with the boundary of the parameter
space [0;1]2, e.g., f(0;v) = 0; 0 v 1.
2. Turning points f = fu = 0 or f = fv = 0, which are illustrated in Figure 10.23. If f is of
degree (M;N) in (u;v), then fu is of degree (M 1;N) and fv is of degree (M;N 1). It
can be seen that the total number of roots of two simultaneous bivariate polynomials of
degree (m;n) and (p;q), respectively, is mq +np.Thus, there can be at most 2MN M
u-turning points and 2MN N v-turning points.
3. Singular points f = fu = fv = 0. Notice fu = rf Ru;fv = rf Rv, and therefore
fu = fv = 0 means that rf k Ru Rv or that the normals of two surfaces are parallel
and since f(u;v) = 0 at these points the two surfaces intersect tangentially. If f is of
degree (M;N) in (u;v), fu is of degree (M 1;N) and fv is of degree (M;N 1), there
can be at most 2MN M N + 1 singular points.
38
f = fu = 0 f = f = 0v
Figure 10.23: Turning points
From the above discussions we can get upper bounds for the maximum number ofu-turning,
v-turning and singular points, see table 10.6. These bounds refer to the maximum possible
number of solutions (u;v) in the entire complex plane. It turns out that the number of such
points in the real square [0;1]2 is much smaller, but still quite large. Consequently methods
which focus on the possible solutions only in the [0;1]2 are advantageous. The subdivision
method of Section 10 is one such method. Interval Newton methods are also potentially useful
in this context.
algebraic curve max number max number max number
S1 S2 f(u;v) degree u-turning pts v-turning pts singular points
M;N 2MN M 2MN N 2MN M N + 1
plane biquadratic 2, 2 6 6 5
plane bicubic 3, 3 15 15 13
quadric biquadratic 4, 4 28 28 25
quadric bicubic 6, 6 66 66 61
torus biquadratic 8, 8 120 120 113
torus bicubic 12, 12 276 276 265
biquadratic biquadratic 16, 16 496 496 481
bicubic biquadratic 36, 36 2556 2556 2521
bicubic bicubic 54, 54 5778 5778 5725
Table 10.6: Number of turning and singular points in various cases
Analysis of singular points:
Let (u0;v0) satisfy f(u0;v0) = 0. We construct a straight line L
u = u0 + t; v = v0 + t (10.26)
and we nd its intersections with the algebraic curve f(u;v) = 0 by substitution as follows:
0 = f(u0 + t;v0 + t)
= f(u0;v0) + tfu(u0;v0) + tfv(u0;v0)
+12( 2t2fuu + 2 t2fuv + 2t2fvv)j(u=u0;v=v0)
39
+h.o.t. (up to t ; is nite)
= t( fu + fv) + t
2
2 (
2fuu + 2 fuv + 2fvv) + h.o.t. (10.27)
(1) Case A: f2u +f2v > 0, f = 0 at (u0;v0).
L is tangent to f = 0 at (u0;v0) if t = 0 is a double root. Thus
fu + fv = 0
= fv
= fu
This is also a proof of fact that rf is ? to curve.
(2) Case B: fu = fv = f = 0 at (u0;v0).
t = 0 is a triple root if f2uu + f2uv + f2vv > 0 and at least one of the 3rd derivatives is
nonzero and
2fuu + 2 fuv + 2fvvju=u0v=v
0
= 0 (10.28)
We can solve this quadratic equation for or and there are three possiblities:
(1) 2 real distinct roots ) 2 distinct tangents (self-intersection)
(2) 1 real double root ) 1 tangent (cusp)
(3) 2 complex roots ) no real tangents (isolated point)
See also Figure 10.24 for illustration of the three cases.
Example 1: Let f(u;v) = u3 +u2 +v2 = 0 so that
fu = u(3u+ 2); fv = 2v; fuu = 6u+ 2; fvv = 2; fuv = 0
Turning points:
(1) f = fu = 0;fv 6= 0 )u = 23
f = 0 )v2 = u2(1 +u) < 0 ) no real solution
(2) f = fv = 0;fu 6= 0 )v = 0;u = 1
Singular points:
f = fu = fv = 0 )u = v = 0
Tangents at u = v = 0 can be obtained from 2 2 + 2 2 = 0; ( )2 + 1 = 0, which has no real
solution, so u = v = 0 is an isolated point.
40
(1)self?intersection (2)cusp (2)cusp
(3)isolated point
(2)cusp
L
L
L L
L
Figure 10.24: Singularities of planar algebraic curves
u
v
1
1
?1
?1
?2
Figure 10.25: Example 1 algebraic curve with an isolated point
If the domain of interest is [ 2;1] [ 1;1], border points are ( 1:465; 1).
Example 2: Let f(u;v) = v2 u3 = 0. This curve has a cusp at u = v = 0 with tangent
v = 0, see Figure 10.26.
Example 3: Let us consider the equation
f(u;v) = (u+ 1)u(u 1)(v + 1)v(v 1) + 120 = 0 (10.29)
within the domain [ 2;2]2. This is a degree 6 algebraic curve illustrated in Figure 10.27.
On every border line segment, there are three border points. The curve has no singular
points, but involves two (internal) loops and six border-to-border branches. The algebraic
curve f(u;v) = 0 in this example has degrees M = 3;N = 3 in u and v. Consequently,
using the previous formulas the number of u turning points, v turning points and singular
41
?1 1
1
?1
u
v
Figure 10.26: Example 2 algebraic curve with a cusp at (u;v) = (0;0)
points (in the entire complex plane) is bounded by 2MN M = 15, 2MN N = 15, and
2MN M N + 1 = 13. However, as we can be seen in Figure 10.27, these numbers
overestimate the actual number of such points in the real square [ 2;2]2.
Computing starting points for all branches
1. Border points: This involves solution of a polynomial, eg.
f(0;v) =
NX
j=0
w0jBj;N(v) = 0
A robust and e cient solution of this type of equation is addressed in Section 10.
2. Turning and singular point computation: Here we use the fact that
fu(u;v) = M
M 1X
i=0
NX
j=0
(wi+1;j wij)Bi;M 1(u)Bj;N(v) (10.30)
fv(u;v) = N
MX
i=0
N 1X
j=0
(wi;j+1 wij)Bi;M(u)Bj;N 1(v) (10.31)
Consequently, computing turning points (f = fu = 0 and f = fv = 0) is equivalent to
solving a system of two nonlinear polynomial equations in two variables, and computing
singularities f = fu = fv = 0 is equivalent to solving a system of three nonlinear
polynomial equations in two variables. Robust and e cient solution of these systems of
nonlinear polynomial equations is addressed in Section 10.
10.8.2 Case F1: RPP/RPP Surface Intersection
In this case we have two rational polynomial surface patches
( R
1 = R1(u;v)
R2 = R2(s;t) (10.32)
42
2
2
?2
?2
Figure 10.27: Example 3 algebraic curve
eg. two rational B ezier patches and by setting R1 = R2 we obtain three nonlinear polynomial
equations for four unknowns u;v;s;t. This system can be solved by the nonlinear solver of
Section 10. However, as the solutions are typically not isolated points but curves, such approach
is ine cient. There are three major techniques for solving RPR/RPP surface intersections.
Method 1: Lattice method
In this method, one of the two surface patches is discretized to a grid of isoparameter
curves at some xed resolution. Each of the resulting curves is intersected with the other
patch (using a technique as in Section 9.7). The resulting solution points are connected
to form various curve branches based on empirical distance-based criteria. The method
is typically ine cient (because it does not use the convex hull properties of RPP surface
patches to their full extent) and generally not robust, leading to missing of small features
of the intersection such as small loops, singularities. Also the connection of the points
to form curves near constrictions and singularities is not robust as it is empirical.
Method 2: Subdivision method
Typically, subdivision methods involve the follwing steps, see Figure 10.28,
Preprocessing by bounding boxes to eliminate subpatches that do not intersect.
Subdivision (typically in four subpatches by multiple knot insertion at the mid point
of parameter axes in NURBS patches)
43
Approximation of the surface with triangular facets (either from the polyhedron or
a grid on the subdivided surface)
Intersections of triangular facets.
Quadtree
Figure 10.28: Subdivision method
Issues:
- Robust: resolution of loops and isolated points (under nite subdivision) is not guar-
anteed.
- E ciency is typically better than lattice methods but usually inferior to marching
methods.
Method 3: Marching along tangent to intersection curve
t = (R1u R1v) (R2s R2t) (10.33)
Issues:
- Finding starting points on all branches.
- Straying, looping, singularities
44
u
v
s
t
Figure 10.29: Parameter spaces of R1(u;v) and R2(s;t)
Let us consider the intersection, R1(u;v) = R2(s;t), eg. as illustrated in the parameter
spaces of R1, R2 in Figure 10.29. A resulting intersection curve branch can be expressed
as a parameter curve in terms of the parameter , as follows
u = u( ); v = v( ); s = s( ); t = t( ) (10.34)
The intersection curve tangent can be obtained as the tangent along the curve on the
surfaces R1(u( );v( )) and R2(s( );t( )) using the chain rule of di erentiation as follows
_R1 = R1u _u+ R1v _v (10.35)
_R2 = R2s _s+ R2t _t (10.36)
where ( _ ) denotes derivatives with respect to . However, _R1 = _R2 and this leads to
R1u _u+ R1v _v = R2s _s+ R2t _t (10.37)
This is a system of three linear equations with four unknowns _u; _v; _s; _t, which can be
solved to provide the following system of rst order nonlinear ordinary di erential equa-
tions (ODE):
ds
d =
x(2)t x(1)u x(1)v
y(2)t y(1)u y(1)v
z(2)t z(1)u z(1)v
= jA1j (10.38)
dt
d =
x(2)s x(1)u x(1)v
y(2)s y(1)u y(1)v
z(2)s z(1)u z(1)v
= jA2j (10.39)
du
d =
x(2)s x(2)t x(1)v
y(2)s y(2)t y(1)v
z(2)s z(2)t z(1)v
= jA3j (10.40)
dv
d =
x(2)s x(2)t x(1)u
y(2)s y(2)t y(1)u
z(2)s z(2)t z(1)u
= jA4j (10.41)
where
R1(u;v) =
h
x(1)(u;v);y(1)(u;v);z(1)(u;v)
i
; (10.42)
R2(s;t) =
h
x(2)(s;t);y(2)(s;t);z(2)(s;t)
i
: (10.43)
45
Here is an arbitrary non-zero factor that can be chosen to provide arc-length parametriza-
tion in the s;t domain as follows:
d = pds2 +dt2 =
q
2(jA1j2 +jA2j2)d
hence
= 1q
jA1j2 +jA2j2
:
This ODE system (10.38) to (10.41) can be solved using the Runge-Kutta method or a
multistep method.
In order to compute approximate starting points for the above marching method we need
to identify rst the possible presence of (internal) loops. This can be done using the
concept of collinear normal points.
Sederberg et al rst recognized the importance of collinear normals in detecting the
existence of closed intersection loops in intersection problems of two distinct parametric
surface patches. Two points on two surfaces are said to be collinear normal points if
their associated normal vectors lie on the same line.
TheoremIf two regular tangent plane continuous surface patches, R1 and R2, intersect
in a closed loop, then there exists a line that is perpendicular to both R1 and R2 if
the dot product of any two normal vectors (on the same patch or on di erent patches)
is never zero. This means that the total range of normal directions for both patches
considered simultaneously can not deviate more than 90 .
In other words, if the two surfaces do not contain a collinear normal (and do not turn
more than 90 ), then there are no closed loops of surface intersection. Denoting the two
surfaces by R1(u;v) and R2(s;t), collinear normal points satisfy the following equations
(R2 R1) R2s = 0; (R2 R1) R2t = 0;
(R2s R2t) R1u = 0; (R2s R2t) R1v = 0: (10.44)
If R1;R2 are RPP surface patches, equation (10.44) form a system of four nonlinear
polynomial equations that can be solved using the method of Section 10.
Now we split the patches in (at least) one parametric direction at these collinear normal
points. Consequently, starting points are only border points on the boundaries of all
subdomains created. Border points are intersections of the border curves of each patch
with the other patch. Their computation involves a system of three nonlinear polynomial
equations in three unknowns, which can be solved with the method of Section 10.
10.8.3 Case F8: IA/IA Surface Intersection
f(R) = g(R) = 0 (10.45)
46
where f;g are polynomials. Here we have two equations and three unknowns R.
One intersection method for low order f;g is to eliminate one variable (e.g. z) to nd
projection of intersection curves on plane of other two variables (e.g. x, y), then trace the
algebraic curve and use the inversion algorithm to nd z. A more complete analysis of this
problem is beyond the scope of these notes.
Example 1: See Figure 10.30
f = x2 +y2 +z2 1 = 0
g = x2 + (y 12)2 14 = 0
x
y
z
x
y1
g=0
y
z
1
1
?1
1
?1
x
z
z=0 plane
y=0 plane
x + z ? z =02 4 2
x=0 plane
y=1?z2
Figure 10.30: Intersection of two implicit quadrics (sphere and cylinder) from example 1
47
Appendix A
Tracing tangent intersection
Let an algebraic curve be such that
f(u;v) = fu(u;v) = fv(u;v) = 0 (10.46)
on all points, and that
f2uu +f2vv +f2uv = 0 (10.47)
and at least one of the 3rd derivative is nonzero.
Then at a point (u0;v0) on f(u0;v0) = 0, the tangent u = u0 + t;v = v0 + t is de ned by
2fuu + 2 fuv + 2fvv = 0 (10.48)
Now assume there is a single real tangent direction on all points of the curve. This occurs
when
f2uv fuufvv = 0 (10.49)
The concept of turning points generalizes to = 0 or = 0, because
= 0 )f = fv = fvv = 0
= 0 )f = fu = fuu = 0
)
de nition of turning points (10.50)
u
vα=0
β=0
Figure 10.31: Turning points
From (10.49), fuufvv 0, so that
fuv =
q
fuufvv
Hence, (10.48) becomes
2fuu + 2
q
fuufvv + 2fvv = 0
48
Case 1 fuu 6= 0
) fuu( )2 2
q
fuufvv( ) +fvv = 0
=
pf
uufvv
fuu =
sf
vv
fuusgn(fuu)
Case 2 fvv 6= 0
) fvv( )2 2
q
fuufvv( ) +fuu = 0
=
pf
uufvv
fvv =
sf
uu
fvv sgn(fvv)
We may choose
=
q
jfvvj = _u
=
q
jfuuj = _v
Normalize so that 2 + 2 = 1
_u = K
q
jfvvj
_v = K
q
jfuuj
9=
;) _u
2 + _v2 = K2(jfuuj+jfvvj) = 1
K = 1q
jfuuj+jfvvj
_u(t) =
q
jfvvj
q
jfuuj+jfvvj
_v(t) =
q
jfuuj
q
jfuuj+jfvvj
Example: f(u;v) = (u2 +v2 1)2
fu = 4u(u2 +v2 1)
fv = 4v(u2 +v2 1)
)
) if f = 0 )fu = fv = 0
fuu = 4(3u2 +v2 1)
fvv = 4(3v2 +u2 1)
fuv = 8uv
9
>=
>; ) if f = 0 )
8
><
>:
fuu = 8u2
fvv = 8v2
f2uv = fuufvv
) _u = jvj_v = juj
)
on f = 0
49
fuu=0
fvv=0
fuu=0
fvv=0
u
v
1
0 1uA
v
Figure 10.32: Tangent intersections
Case of in nite turning points
f(u;v) = (u A)kg(u;v) = 0
1. f(u;0) = 0f(u;1) = 0
)
) common solution u = A,
2. On which fv = (u A)kgv(u;v) = 0.
Try factoring (u A) sequentially until the v derivative of ratio is nonzero.
Condition (a) and (b) are necessary but not su cient. However, subdivision would work
in the limit, in determining f = fv = 0.
50
Nonlinear solvers and robustness issues
10.9 Nonlinear Solvers
10.9.1 Motivation
As we have seen in earlier chapters, the geometric shape of curves and surfaces is usually
represented by polynomial equations of various types (eg., implicit or parametric). As we
have seen in Chapter 9, intersection problems reduce to solving systems of nonlinear equations
which are usually polynomial if the geometries involved are polynomial.
Occassionally, the governing equations for general interrogation problems (intersections,
curvature extrema, etc.) reduce to systems of nonlinear equations, involving also square roots
of polynomials, which arise from normalization of the normal vector and analytical expressions
of the curvatures.
Example 1 : Intersection between two planar implicit polynomial (algebraic) curves, eg.
two circles, is shown in Figure 10.33. Their equations are
x2 +y2 916 = 0
(x 1)2 +y2 14 = 0
which form a nonlinear polynomial system. The roots can be obtained by eliminating y, solving
for x and then backsubstituting and solving for y, leading to
(x;y) = (2132;
p135
32 ) ’ (0:65625; 0:36309)
Example 2 : Self-intersection of o set of a parabola is shown in Figure 10.34. Such
intersections are needed in planning NC machining. Self-intersections of a parametric o set
curve can be obtained by seeking pairs of distinct parameter values s6= t such that
r(s) +dn(s) = r(t) +dn(t) (10.51)
where r(s) is a planar (progenitor) parametric curve, d is an o set distance and n(s) is a unit
normal vector to r(s) given by
n = t ez = ( _y(s); _x(s))q
_x2(s) + _y2(s)
(10.52)
51
y
x
1
1
?1
?1
Figure 10.33: Intersection between two circles.
where t is the unit tangent vector and ez the unit normal to the plane of the progenitor curve.
Figure 10.34: Self-intersection of an o set of a parabola.
Substituting equation (10.52) into equation (10.51) yields
x(s) + _y(s)dq
_x2(s) + _y2(s)
= x(t) + _y(t)dq
_x2(t) + _y2(t)
y(s) _x(s)dq
_x2(s) + _y2(s)
= y(t) _x(t)dq
_x2(t) + _y2(t)
(10.53)
If r(s) is a polynomial function, equations 10.53 form a system of two nonlinear equations
involving polynomials and radicals of polynomials. But if we set 2(s) = _x2(s) + _y2(s), we
obtain 3 polynomial equations in which is the auxiliary variable.
52
10.10 Local Solution Methods
They are designed to compute roots based on initial approximations.
Newton’s Method in one variable [6]
We want to nd roots for f(x) = 0. Iteration formula is given by
xn+1 = xn f(xn)f0(x
n)
; n = 0;1;2::::: (10.54)
This is illustrated in Figure 10.35(a). A modi ed Newton’s method for one variable is
shown in Figure 10.35(b), where we take a fractional step as follows in order to reduce the
possibility of divergence in Figure 10.35(b) of the full step method given by equation 10.54
xn+1 = xn f(xn)f0(x
n)
(10.55)
where = max[1; 12; 14;:::] such that jf(xn+1)j<jf(xn)j.
f(x)
xx
x
n
n+1
f(x ) n
f(x )n+1
(a) ( b)
xnxn+1
Figure 10.35: Newton’s method for f(x) = 0
Newton’s Method for m equations in m unknowns
We want to nd the roots of the system
f(r) = 0 (10.56)
where r = [r1;r2;:::;rm]T , and f = [f1;f2;:::;fm]T .
Iteration formula is given by
r(n+1) = r(n) + r(n) (10.57)
where
J(r(n)) r(n) = f(r(n)) (10.58)
and J = [@fi@rj (r(n))] is the Jacobian matrix [6].
53
Advantages
Quadratic convergence.
Straightforward to program.
Disadvantages
Needs good initial approximations for each root, otherwise it may diverge.
Cannot provide full assurance that all roots have been found.
Example Two intersecting circles x2 +y2 = 916, (x 1)2 +y2 = 14. Let
f(x;y) = x2 +y2 916 = 0 (10.59)
g(x;y) = (x 1)2 +y2 14 = 0 (10.60)
Then, the Jacobian matrix is evaluated as follows:
[J] = 2
" x y
x 1 y
#
(10.61)
10.11 Classi cation of Global Solution Methods
Global solution methods are designed to compute all roots in some area of interest. In recent
computational geometry related research, three classes of methods for the computation of
solutions of nonlinear polynomial systems have been favored:
Algebraic techniques [4] [5]
Advantages: Theoretically elegant, and well-suited for implementation in symbolic math-
ematical systems. They determine all roots with arbitrary precision if the coe cients of
the polynomials involved are integers or rational numbers.
Disadvantages: Ine cient in memory and processing time requirements and therefore
unattractive except for very low degree or dimensionality systems.
Homotopy methods [9] [28]
They tend to be numerically ill-conditioned for high degree polynomials and similarly
ine cient because they are exhaustive.
Subdivision methods [18] [10] [23] [27]
Advantages: Subdivision methods are generally e cient and stable. Therefore, they are
likely to be the most successful methods in practice. They can be combined with interval
methods to numerically guarantee that certain subdomains do not contain solutions.
Disadvantages: They are not as general as algebraic methods, since they are only capable
of isolating zero-dimensional solutions.
Furthermore, although the chances, that all roots have been found, increase as the res-
olution tolerance is lowered, there is no certainty that each root has been extracted.
Subdivision methods typically do not provide a guarantee as to how many roots there
may be in the remaining subdomains. However, if these subdomains are very small the
54
existence of a (single) root within these subdomains is a typical assumption.
Lastly, subdivision techniques provide no explicit information about root multiplicities
without additional computation.
10.12 Subdivision Method (Projected Polyhedron Method)
We want to nd roots of degree m polynomial equation f(x) = co +c1x+c2x2 + +cmxm = 0
over the region a x b.
The procedure is to rst make the a ne parameter transformation x = a+t(b a) such
that 0 t 1. The transition from the interval a x b to the interval 0 x 1 is
an a ne map, and the polynomials are invariant under a ne parameter transformation.
Now the polynomial equation is given by the monomial form
f(t) =
mX
i=0
cMi ti (10.62)
Change the basis from monomial to Bernstein.
f(t) =
mX
i=0
cBi Bi;m(t) (10.63)
with cBi =
iX
j=0
(ij)
(mj )c
M
j
where Bi;m(t) is the ith Bernstein polynomial given by
Bi;m(t) = (mi )ti(1 t)m i (10.64)
Use linear precision property of the Bernstein polynomial
t =
mX
i=0
i
mBi;m(t) (10.65)
In other words, the monomial t can be expressed as the weighted sum of Bernstein
polynomials with coe cients evenly spaced in the interval [0;1].
Create a graph of function f(t). Then the graph will become a B ezier curve
f(t) =
t
f(t)
!
=
MX
i=0
i
mcB
i
!
Bi;m(t) (10.66)
where ( im;cBi )T are control points.
Now the problem of nding roots of the univariate polynomial has been transformed into
a problem of nding the intersection of the B ezier curve with the parameter axis.
55
Using the convex hull property, we discard the regions which do not contain the roots
by applying the de Casteljau subdivision algorithm and nd a sub-region of [0,1] which
contain the root.
If the sub-region is su ciently small, we conclude that there is a root inside and return
it.
But when there are more than one root in the sub-region, the sub-region will not be
reduced. In such case we split the region evenly.
Scale the sub-region so that it will become [0,1] using a ne parameter transformation
and go back to *.
Example 1: Projected Polyhedron Method in one variable
Find the roots of f(x) = 1:1x2 + 1:4x 0:2 = 0 where 0 x 2. The roots are
approximately, 0.164 and 1.108.
A ne parameter transformation
Plug x = 0 + (2 0)t = 2t into f(x)
f(t) = 4:4t2 + 2:8t 0:2 = 0 where 0 t 1
Monomial to Bernstein basis
We have cM0 = 0:2, cM1 = 2:8 and cM2 = 4:4, thus using
cBi =
iX
j=0
(ij)
(mj )c
M
j (10.67)
we obtain cB0 = 0:2, cB1 = 1:2 and cB2 = 1:8
Use linear precision property of the Bernstein polynomial
t =
2X
i=0
i
2Bi;2(t) (10.68)
Create a graph of function f(t). Then the graph will become a B ezier curve
f(t) =
t
f(t)
!
=
2X
i=0
i
2cB
i
!
Bi;2(t) (10.69)
Now the control points of the B ezier curve is as follows:
(0, -0.2), (0.5, 1.2) and (1, -1.8).
* Construct a convex hull of the B ezier curve. In this example it is a triangle whose
vertices are the control points (0, -0.2), (0.5, 1.2) and (1, -1.8).
56
The convex hull intersects the t-axis with t = 0:0714 and t = 0:7.
Discard the regions 0 t 0:0714 and 0:7 t 1, which do not contain roots by
appling the de Casteljau algorithm.
Now we have a smaller convex hull which contains the roots. (Shaded triangular in
Figure 10.36)
Since there are two roots in the convex hull the sub-region will not reduce much. There-
fore we split the region evenly and scale the two boxes so that it will become [0,1] using
the a ne parameter transformation and repeat the process until the two sub-regions
become small enough.
t
f(t)
t=0.0714 t=0.7
(0.5, 1.2)
(0, ?0.2)
(1, ?1.8)
Figure 10.36: de Casteljau Algorithm applied to the quadratic B ezier curve
Example 2: Projected Polyhedron Method in two variables
Two intersecting circles x2 +y2 = 916, (x 1)2 +y2 = 14, see Figure 10.33. Let
f(x;y) = x2 +y2 916 = 0; 1 x; y 1 (10.70)
g(x;y) = (x 1)2 +y2 14 = 0; 1 x; y 1 (10.71)
A ne parameter transformation
s = x ( 1)1 ( 1) = x+ 12 (10.72)
t = y ( 1)1 ( 1) = y + 12 (10.73)
(10.74)
Plug x = 2s 1 and y = 2t 1 into equations (10.70) (10.71), then we have
f(s;t) = 4s2 4s+ 4t2 4t+ 2316 = 0 (10.75)
g(s;t) = 4s2 8s+ 4t2 4t+ 194 = 0 (10.76)
57
Monomial to Bernstein basis
cBij =
iX
k=0
jX
l=0
(ik)(jl )
(mk )(nl )c
M
kl (10.77)
We obtain
f(s;t) =
2X
i=0
2X
j=0
fBijBi;2(s)Bi;2(t) (10.78)
g(s;t) =
2X
i=0
2X
j=0
gBijBi;2(s)Bi;2(t) (10.79)
where
fB00 = 1:4375; fB01 = 0:5625; fB02 = 1:4375
fB10 = 0:5625; fB11 = 2:5625; fB12 = 0:5625
fB20 = 1:4375; fB21 = 0:5625; fB22 = 1:4375
and
gB00 = 4:75; gB01 = 2:75; gB02 = 4:75
gB10 = 0:75; gB11 = 1:25; gB12 = 0:75
gB20 = 0:75; gB21 = 1:25; gB22 = 0:75
Use linear precision property of the Bernstein polynomial and create a graph.
f(s;t) =
0
B@ st
f(s;t)
1
CA = 2X
i=0
0
B@
i
2j
2fB
ij
1
CAB
i;2(s)Bi;2(t) (10.80)
g(s;t) =
0
B@ st
g(s;t)
1
CA = 2X
i=0
0
B@
i
2j
2gB
ij
1
CAB
i;2(s)Bi;2(t) (10.81)
The two B ezier surfaces are shown in Figure 10.37.
Now we will nd the intersections of three surfaces, f(s;t), g(s;t) and xy-plane. Figure
10.38 shows the intersection between the plane and both B ezier surfaces. We can easily
observe that the two intersection curves are the circles in Figure 10.33.
Project the control points of f(s;t) and g(s;t) into xz and yz planes.
For each xz and yz plane, construct the 2D convex hulls. Solid line corresponds to f(s;t)
and the dotted line corresponds to g(s;t).
58
x
y
z
x
yz
Figure 10.37: B ezier surfaces and their control points
Figure 10.38: B ezier surfaces intersecting with xy-plane
Figure 10.39: Projections of Control Points
59
Intersect the convex hull with horizontal axis. Because the polygon is convex, the inter-
section may be either a closed interval (which may degenerate to a point) or empty. If
it is empty, then no root of the system exists within the given search box.
Intersect the intervals with one another. Again, if the result is empty, no root exists
within the given search box.
10.13 Interval Methods
10.13.1 Motivation
Nonlinear solvers operating in rational arithmetic are robust, but are generally time-consuming.
On the other hand, nonlinear solvers operating in oat point arithmetic are faster, but not
robust. Interval methods solve the two problems, namely, nonlinear solvers operating in interval
arithmetic are inexpensive compared to rational arithmetic, and they are robust.
x
y
o
Figure 10.40: Curves y = x4 and y = 0 contact tangentially at the origin.
Example
Suppose we have a degree four planar B ezier curve whose control points are given by
( 0:5;0:0625); ( 0:25; 0:0625); (0;0:0625); (0:25; 0:0625); (0:5;0:0625) (10.82)
as shown in Figure 10.40. This B ezier curve is equivalent to the explicit curve y = x4 ( 0:5
x 0:5). Apparently the curve intersects with x-axis tangentially at (x;y) = (0;0). However,
if the curve has been translated by +1 in the y direction and translated back to the original
position by moving by 13 three times during a geometric processing session, the curve will
generally not be the same as the original curve if oating point arithmetic (FPA) is used for the
computation. For illustration, let us assume a decimal computer with a four-digit normalized
mantissa, and the computer rounds o intelligently rather than truncating. Then the rational
number 13 will be stored in the decimal computer as 0:3333 100 and after the processing
the new control points will be
( 0:5;0:0631); ( 0:25; 0:0624); (0;0:0631); (0:25; 0:0624); (0:5;0:0631) (10.83)
60
If we evaluate the curve at parameter value t = 0:5, we obtain (0, 0.00035) instead of (0,0).
Therefore there exists a numerical gap which could later lead to inconsistency between topo-
logical structures and geometric equations. For example, if these new control points are used
for computing intersections with the x-axis, the computer will return no solutions when the
tolerance is smaller than 0.00035. The above problem illustrates the case when the error is
created during the formulation of the governing equations by various algebraic transformations.
10.13.2 De nitions
An interval is a set of real numbers de ned below [21]:
[a;b] = fxja x bg (10.84)
Two intervals [a;b] and [c;d] are said to be equal if
a = c and b = d (10.85)
The intersection of two intervals is empty or [a;b] \[c;d] = ;, if either
a>d or c>b (10.86)
Otherwise,
[a;b]\ [c;d] = [max(a;c);min(b;d)] (10.87)
The union of the two intersecting intervals is
[a;b] [[c;d] = [min(a;c);max(b;d)]: (10.88)
An order of intervals is de ned by
[a;b] < [c;d] if and only if b<c: (10.89)
The width of an interval [a;b] is b a.
The absolute value is
j[a;b]j = max(jaj;jbj): (10.90)
Examples
[2;4]\ [3;5] = [max(2;3);min(4;5)] = [3;4]
[2;4][ [3;5] = [min(2;3);max(4;5)] = [2;5]
j[ 7; 2]j = max(j 7j;j 2j) = 7
61
10.13.3 Interval Arithmetic
[a;b] [c;d] = fx y jx2 [a;b] andy 2 [c;d]g: (10.91)
where represents an arithmetic operation 2 f+; ; ;=g. Using the end points of the two
intervals, we can rewrite equation (10.91) as follows
[a;b] + [c;d] = [a+c;b+d]
[a;b] [c;d] = [a d;b c]
[a;b] [c;d] = [min(ac;ad;bc;bd);max(ac;ad;bc;bd)]
[a;b]=[c;d] = [min(a=c;a=d;b=c;b=d);max(a=c;a=d;b=c;b=d)] (10.92)
provided 0 62 [c;d] in the division relation.
Examples
[2;4] + [3;5] = [2 + 3;4 + 5] = [5;9]
[2;4] [3;5] = [2 5;4 3] = [ 3;1]
[2;4] [3;5] = [min(2 3;2 5;4 3;4 5);max(2 3;2 5;4 3;4 5)] = [6;20]
[2;4]=[3;5] = [min(2=3;2=5;4=3;4=5);max(2=3;2=5;4=3;4=5)] = [2=5;4=3]
10.13.4 Algebraic Properties
Interval arithmetic is commutative and associative.
[a;b] + [c;d] = [c;d] + [a;b]
[a;b] [c;d] = [c;d] [a;b]
[a;b] + ([c;d] + [e;f]) = ([a;b] + [c;d]) + [e;f]
[a;b] ([c;d] [e;f]) = ([a;b] [c;d]) [e;f]
But it is not distributive, however, it is subdistributive.
Examples
[1;2]([1;2] [1;2]) = [1;2]([ 1;1]) = [ 2;2]
[1;2][1;2] [1;2][1;2] = [1;4] [1;4] = [ 3;3]
10.13.5 Rounded Interval Arithmetic and its Implementation
One has to keep in mind that any sequence of operations on a digital computer is essentially
equivalent to a nite sequence of manipulations on a discrete grid of points. For example, a
oating point number in the general form is given by [11]
( ):d1d2 dp 2exp,
62
where d1d2 dp is mantissa with d1 6= 0, and p is the number of signi cant digits, and exp
is the integer exponent. If p = 2 and 2 exp 3 then a list of positive numbers in this
system is
:10 2 2 = 18 :11 2 2 = 316
:10 2 1 = 14 :11 2 1 = 38
:10 20 = 14 :11 20 = 34
:10 21 = 1 :11 21 = 32
:10 22 = 2 :11 22 = 3
:10 23 = 4 :11 23 = 6
0 1?
8
1
?
4
3
?
8
1
?
2
3
?
4
3
?
2
3
?
16
1 2 3 4 6
Figure 10.41: Nonnegative oating-point numbers on the interval [0,6], adapted from [11]
If oating{point arithmetic is used to evaluate these interval arithmetic equations there is
no guarantee that the roundings of the bounds are performed conservatively. 1.
Most commercial processors implement oating{point arithmetic using the representation
de ned by ANSI/IEEE Std 754{1985, Standard for Binary Floating{Point Arithmetic [2]. This
standard de nes the binary representation of the oating{point number X in terms of a sign
bit s, an integer exponent E, for Emin E Emax, and a p{bit signi cand (or mantissa) B,
where:
X = ( 1)s2EB (10.93)
The signi cand B is a sequence of p bits b0b1 bp 1, where bi = 0 or 1, with an implied
binary point (analogous to a decimal point) between bits b0 and b1. Thus, the value of B is
calculated as:
1This statement is true only for the default IEEE-754 rounding mode of round towards nearest
63
B = b0:b1b2 bp 1 = b020 +
p 1X
1
bi2 i (10.94)
For double precision arithmetic, the standard de nes p = 53, Emin = 1022, and Emax =
1023. The number X is represented as a 64{bit quantity with sign bit s, an 11{bit biased ex-
ponent e = E+1023, and a 52{bit fractional mantissa m composed of the bit string b1b2 b52.
Since the exponent can always be selected such that b0 = 1 (and thus, 1 B < 2), the value
of b0 is constant and it does not need to be stored in the binary representation.
63 62 52 51 0
s e m
The integer value of the 11-bit biased exponent e is calculated as:
e = e0e1 e10 =
10X
0
ei210 i (10.95)
ULP (One Unit in the Last Place)
If x and x0 are consecutive positive double-precision numbers, they di er by an amount
called ulp. To calculate ulp it is necessary to extract the integer value of the exponent from the
binary representation. Recall that the value of the signi cand B of a double precision number
X is:
B = 1 +b12 1 +b22 2 + +b522 52 (10.96)
and that the double precision value X = ( 1)s2EB. The value of the least signi cant bit b52
is 2 52. Thus, the value of ulp is 2E2 52 = 2E 52.
Rounded interval arithmetic [19, 20] ensures that the computed end points always contain
the exact interval as follows:
[a;b] + [c;d] [a+c ‘;b+d+ u]
[a;b] [c;d] [a d ‘;b c+ u]
[a;b] [c;d] [min(a c;a d;b c;b d) ‘;max(a c;a d;b c;b d) + u]
[a;b] = [c;d] [min(a=c;a=d;b=c;b=d) ‘;max(a=c;a=d;b=c;b=d) + u]
(10.97)
where ‘ and u are the units{in{last{place and are denoted by ulp‘ and ulpu. When per-
forming standard operations for interval numbers using RIA, the lower bound is extended to
include its previous consecutive oating{point number, which is smaller than the lower bound
by ulp‘. Similarly, the upper bound is extended by ulpu to include its next consecutive num-
ber. Thus, the width of the result is enlarged by ulp‘ +ulpu and the result will be reliable in
subsequent operations.
Example
64
main()
f
double a = 1.5;
Interval b = 1.5;
double dresult = pow(a, 20.);
Interval iresult = pow(b, 20.);
g
dresult 3325.2567300796509
iresult [3325.2567300796404,3325.2567300796613]
10.14 Interval Projected Polyhedron Algorithm
We extend the de Casteljau subdivision method to operate in rounded interval arithmetic in
order to nd all the roots of a polynomial system robustly. We illustrate the concept for a
single polynomial equation.
(x - 0.1)(x - 0.6)(x - 0.7) = 0
65
Floating Point Arithmetic
Iter Bounding Box (FPA) Message
1 [0,1]
2 [0.0763636363636364, 0.856]
3 [0.098187732239346, 0.770083868323999]
4 [0.0999880766853688, 0.72387404781026] Binary Sub.
5 [0.402239977003124, 0.704479954527487]
6 [0.550441290533288, 0.700214508664293]
7 [0.591018492648952, 0.700000534482207]
8 [0.599458794784619, 0.700000000003332] Binary Sub.
9 [0.649998841568898, 0.699999999999999] No Root
10 [0.599997683137796, 0.649998841568898] Root Found
11 [0.099999999478761, 0.402239977003124] Root Found
t
f(t)
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a2a0a1a0
[P1]
[P2][P0]
t=tAL t=tBU
Figure 10.42: Interval Projected Polyhedron Method
Rounded Interval Arithmetic
Iter Bounding Box (RIA) Message
1 [0, 1]
2 [0.076363636363635, 0.856000000000001]
3 [0.0981877322393447, 0.770083868324001]
4 [0.0999880766853675, 0.723874047810262] Binary Sub.
5 [0.402239977003124, 0.704479954527489]
6 [0.550441290533286, 0.700214508664294]
7 [0.591018492648947, 0.700000534482208]
8 [0.599458794784611, 0.700000000003333] Binary Sub.
9 [0.649998841568894, 0.7] Root Found
10 [0.599997683137788, 0.649998841568895] Root Found
11 [0.0999999994787598, 0.402239977003124] Root Found
Most applications mentioned above result in n nonlinear polynomial equations with n
unknowns, referred to as balanced systems. However, there exist some problems such as
66
tangential and overlapping intersection or implicit curve/surface rendering consisting of n
nonlinear polynomial equations with m unknowns, where n could be larger or smaller than m.
When n>m the system is called overconstrained and when n<m it is called underconstrained.
Here is an algorithm (the interval n-dimensional Projected-Polyhedron algorithm) such that
it can e ectively handle overconstrained problems [15].
Suppose we solve a system of nonlinear interval polynomial equations [f] = ([f1];[f2];:::;[fn]) =
0 over the box S 2 Rm (n>m, n = m, n<m) where S is de ned by
S = [a1;b1] [a2;b2] ::: [am;bm]: (10.98)
and each [fi] is an interval polynomial in m variables. That is, we wish to nd all u 2S such
that
[f1](u) = [f2](u) = ::: = [fn](u) = 0: (10.99)
By making the a ne parameter transformation [13] ui = ai + xi(bi ai) for i = 1; m, we
simplify the problem to the problem of determining all x 2 [0;1]m such that
[f1](x) = [f2](x) = ::: = [fn](x) = 0: (10.100)
Now furthermore suppose that each of the [fk] is polynomial in the independent parameters
x1;x2;:::;xm. Let d(k)i denote the degree of [fk] in the variable xi; then [fk] can be written in
the multivariate Bernstein basis:
[fk](x) =
d(k)1X
i1=0
d(k)2X
i2=0
:::
d(k)mX
im=0
[w(k)i1i2:::im]Bi
1;d(k)1
(x1)Bi
2;d(k)2
(x2):::Bi
n;d(k)m
(xm): (10.101)
where Bi;m is the ith Bernstein polynomial. The notation in (10.101) may be simpli ed by
letting I = (i1;i2;:::im), D(k) = (d(k)1 ,d(k)2 ,:::,d(k)m ) and writing (10.101) in the equivalent form
[27]
[fk](x) =
D(k)X
I
[w(k)I ]BI;D(k)(x): (10.102)
Here we have merely rewritten the product of Bernstein polynomials as a single Bernstein
multinomial BI;D(k)(x). Bernstein polynomials have a useful identity called linear precision
property [13], in other words, the monomial t can be expressed as the weighted sum of Bernstein
polynomials with coe cients evenly spaced in the interval [0;1]. Using this property, we can
rewrite equations (10.102) as follows:
[Fk](x) =
D(k)X
I
[v(k)I ]BI;D(k)(x) (10.103)
where
[v(k)I ] = ([ i1d(k)
1
];[ i2d(k)
2
];:::;[ imd(k)
m
];[w(k)I ]): (10.104)
These [v(k)I ] are called the interval control points of [Fk]. Now the algebraic problem of nding
roots of systems of polynomials has been transformed to the geometric problem involving the
67
intersection of the hypersurfaces. Because the problem is now phrased geometrically, we can
use the convex hull property of the multivariate Bernstein basis to bound the set of roots. We
can structure a root- nding algorithm as follows:
1. Start with an initial box of search.
2. Scale the box and, as we did in converting between equations (10.99) and (10.100),
perform an appropriate a ne parameter transformation to the functions fk, so that the
box becomes [0;1]m. However, keep track of the scaling relationship between this box
and the initial box of search. This transformation can be performed with multivariate
De Casteljau subdivision.
3. Using the convex hull property, nd a sub-box of [0;1]m which contains all the roots.
The essential idea behind the box generation scheme in this algorithm is to transform
a complicated m + 1-dimensional problem into a series of m two-dimensional problems.
Suppose Rm+1 can be coordinatized with the x1;x2;:::;xm+1 axes; we can then employ
these steps:
(a) Project the [v(k)I ] of all of the [Fk] into m di erent coordinate planes; speci cally,
the (x1;xm+1)-plane, the (x2;xm+1)-plane, and so on, up to the (xm;xm+1) plane.
(b) In each one of these planes,
i. Construct n two-dimensional convex hulls. The rst is the convex hull of the
projected control points of [F1], the second is from [F2] and so on.
ii. Intersect each convex hull with the horizontal axis (that is, xm+1 = 0). Because
the polygon is convex, the intersection may be either a closed interval (which
may degenerate to a point) or empty. If it is empty, then no root of the system
exists within the given search box.
iii. Intersect the intervals with one another. Again, if the result is empty, no root
exists within the given search box.
(c) Construct an m-dimensional box by taking the Cartesian product of each one of
these intervals in order. In other words, the x1 side of the box is the interval
resulting from the intersection in the (x1;xm+1)-plane, and so forth.
4. Using the scaling relationship between our current box and the initial box of search, see
if the new sub-box represents a su ciently small box in Rm. If it does not, then go to
step 5. If it does, then check the convex hulls of the hypersurface in the new box. If
the convex hulls cross each variable axis, conclude that there is a root or at least an
approximate root in the new box, and put the new box into a root list. Otherwise the
new box is discarded (see Appendix for elaboration of this point).
5. If any dimension of this sub-box are not much smaller than 1 unit in length (i.e., the
box has not decreased much in size along one or more sides), split the box evenly along
each dimension which is causing trouble. Continue on to the next iteration with several
independent sub-problems.
68
6. If none of the boxes is left, then the root- nding process is over. Otherwise, go back to
step 2, and perform it once for each new box.
If we assume that each equation 10.99 is degree m in each variable and the system is n-
dimensional, then the total asymptotic time per step is of O(n2mn+1). The number of steps
depends primarily on the accuracy required [27].
The Projected Polyhedron Algorithm achieves quadratic convergence in one dimension,
while for higher dimensions, it exhibits at best linear convergence.
10.15 Interval Newton method
Interval Newton methods 2 have been the focus of signi cant attention in numerical analysis.
A thorough review of various types of interval Newton methods is presented in [3]. In the
sequel we brie y review the interval Newton method. The interval Newton method solves a
system of nonlinear equations in a numerically veri able manner.
fi(x1;x2; ;xn) = 0; 1 i n (10.105)
within boxes
ai xi bi; 1 i n (10.106)
If we denote then-vector whoseith component isxi by X and then-vector whoseith component
is fi by F(X), the interval Newton methods can be described as nding a box Xk that contains
all the solutions of the interval linear system
J(Xk)( Xk Xck) = F(Xck) (10.107)
where the subscript k denotes the kth iteration, J(Xk) is the Jacobian matrix of F over the
box Xk, and Xck is some point in Xk. There is a theorem about unique solution in the box
given by [16]:
Theorem
If Xk is strictly contained in Xk, then the system of equations in (10.105) has a unique
solution in Xk, and Newton’s method starting from any point in Xk will converge to that
solution. Conversely, if Xk \ Xk is empty, then there are no solutions of the system in
(10.105).
The next iteration Xk+1 is evaluated by
Xk+1 = Xk \ Xk (10.108)
According to the mean value theorem, the solutions in the Xk must be in Xk+1. If the
coordinate intervals of Xk+1 are smaller than those of Xk, equations (10.107) and (10.108)
2When we use the term \interval Newton method", we assume that bisection is included in the process
when interval reduction is not substantial by the pure interval Newton step.
69
are iterated until the bounding boxes are smaller than a speci ed tolerance. If the coordinate
intervals of Xk+1 are not smaller than those of Xk, then one of these intervals is bisected to
form two new boxes. The boxes are pushed into a stack and iteration is continued until the
stack becomes empty. The rst interval Newton method introduced by Moore [21] involves
computing the inverse of the interval matrix J(Xk). Hansen [12] introduced the Gaussian
elimination procedure to solve the linear equation system in an interval Newton method.
Krawczyk [12] introduced a variation of the interval Newton method which avoids the Gaussian
elimination of an interval matrix by not attempting to obtain a sharp solution of (10.107). He
computes the new box Xk as follows
Xk = K(Xk) = Xck YkF(Xck) + (I YkJ(Xk))(Xk Xck) (10.109)
where Yk is a preconditioned matrix of midpoints of the elements of the interval Jocobian
matrix. Hansen and Sengupta [12] introduced a box which is generally smaller than K(Xk).
They simply solve the ith equation for the ith variable and replace the others by bounding
intervals, which is the non-linear version of the Gauss-Seidel operator for linear systems. Let xci
be the ith component of Xck and ki be the ith component of YkF(Xk) and Gij be the entry in
the ith row and jth column of YkJ(Xk) then, the step for the ith row of the Hansen-Sengupta
operator becomes
xi = xi [Gii] 1[ki +
i 1X
j=1
Gij(^xj xcj) +
nX
j=i+1
Gij(xj xcj)] (10.110)
^xi = xi \ xi
for i = 1; ;n.
70
Bibliography
[1] S. L. Abrams, W. Cho, C.-Y. Hu, T. Maekawa, N. M. Patrikalakis, E. C. Sherbrooke, and
X. Ye. E cient and reliable methods for rounded-interval arithmetic. Computer-Aided
Design, 30(8):657{665, July 1998.
[2] ANSI/IEEE Std 754{1985. IEEE Standard for Binary Floating{Point Arithmetic. IEEE,
New York, 1985. Reprinted in ACM SIGPLAN Notices, 22(2):9-25, February 1987.
[3] C. Bliek. Computer Methods for Design Automation. PhD thesis, Massachusetts Institute
of Technology, Cambridge, MA, July 1992.
[4] B. Buchberger. Gr obner bases: An algorithmic method in polynomial ideal theory. In
N. K. Bose, editor, Multidimensional Systems Theory: Progress, Directions and Open
Problems in Multidimensional Systems, pages 184{232. Dordrecht, Holland: D. Reidel
Publishing Company, 1985.
[5] J. F. Canny. Generalized characteristic polynomials. Journal of Symbolic Computation,
9:241{250, 1990.
[6] G. Dahlquist and A. Bj orck. Numerical Methods. Prentice-Hall, Inc., Englewood Cli s,
NJ, 1974.
[7] R. T. Farouki and V. T. Rajan. Algorithms for polynomials in Bernstein form. Computer
Aided Geometric Design, 5(1):1{26, June 1988.
[8] I. D. Faux and M. J. Pratt. Computational Geometry for Design and Manufacture. Ellis
Horwood, Chichester, England, 1987.
[9] C. B. Garcia and W. I. Zangwill. Global continuation methods for nding all solutions
to polynomial systems of equations in n variables. In A. V. Fiacco and K. O. Kortanek,
editors, Extremal Methods and Systems Analysis, pages 481{497. Springer-Verlag, New
York, NY, 1980.
[10] A. Geisow. Surface Interrogations. PhD thesis, School of Computing Studies and Accoun-
tancy, University of East Anglia, Norwich NR47TJ, U. K., July 1983.
[11] C. F. Gerald and P. O. Wheatley. Applied Numerical Analysis. Addison-Wesley, Reading,
MA, 4th edition, 1990.
71
[12] E. Hansen and S. Sengupta. Bounding solutions of systems of equations using interval
analysis. BIT, 21:201{211, 1981.
[13] J. Hoschek and D. Lasser. Fundamentals of Computer Aided Geometric Design. A.
K. Peters, Wellesley, MA, 1993. Translated by L. L. Schumaker.
[14] C. Y. Hu, T. Maekawa, N. M. Patrikalakis, and X. Ye. Robust interval algorithm for
surface intersections. Computer-Aided Design, 29(9):617{627, September 1997.
[15] C. Y. Hu, T. Maekawa, E. C. Sherbrooke, and N. M. Patrikalakis. Robust interval algo-
rithm for curve intersections. Computer-Aided Design, 28(6/7):495{506, June/July 1996.
[16] R. B. Kearfott and M. Novoa. INTBIS, a portable interval Newton/bisection package
(algorithm 681). ACM Transactions on Mathematical Software, 16(2):152{157, June 1990.
[17] G. A. Kriezis, P. V. Prakash, and N. M. Patrikalakis. Method for intersecting algebraic
surfaces with rational polynomial patches. Computer-Aided Design, 22(10):645{654, De-
cember 1990.
[18] J. M. Lane and R. F. Riesenfeld. Bounds on a polynomial. BIT: Nordisk Tidskrift for
Informations-Behandling, 21(1):112{117, 1981.
[19] T. Maekawa and N. M. Patrikalakis. Computation of singularities and intersections of
o sets of planar curves. Computer Aided Geometric Design, 10(5):407{429, October 1993.
[20] T. Maekawa and N. M. Patrikalakis. Interrogation of di erential geometry properties for
design and manufacture. The Visual Computer, 10(4):216{237, March 1994.
[21] R. E. Moore. Interval Analysis. Prentice-Hall, Englewood Cli s, NJ, 1966.
[22] M. E. Mortenson. Geometric Modeling. John Wiley and Sons, New York, 1985.
[23] T. Nishita, T. W. Sederberg, and M. Kakimoto. Ray tracing trimmed rational surface
patches. ACM Computer Graphics, 24(4):337{345, August 1990.
[24] N. M. Patrikalakis. Surface-to-surface intersections. IEEE Computer Graphics and Ap-
plications, 13(1):89{95, January 1993.
[25] N. M. Patrikalakis and P. V. Prakash. Surface intersections for geometric modeling.
Journal of Mechanical Design, Transactions of the ASME, 112(1):100{107, March 1990.
[26] H. Pottmann. General o set surfaces. Neural, Parallel and Scienti c Computations,
5:55{80, 1997.
[27] E. C. Sherbrooke and N. M. Patrikalakis. Computation of the solutions of nonlinear
polynomial systems. Computer Aided Geometric Design, 10(5):379{405, October 1993.
[28] W. I. Zangwill and C. B. Garcia. Pathways to solutions, xed points, and equilibria.
Prentice-Hall, Englewood Cli s, NJ, 1981.
72
| 课件名称： | Computational Geometry |
| 课件分类： | 航空航天 |
| 课件类型： | 教学课件 |
| 文件大小： | 2.7MB |
| 下载次数： | 0 |
| 评论次数： | 0 |
| 用户评分： | 0 |
- 1. Computational Geometry：lecnotes10
- 2. Computational Geometry：lecnotes13
- 3. Computational Geometry：lecnotes14_fixed
- 4. Computational Geometry：lecnotes19
- 5. Computational Geometry：lecnotes1_fixed
- 6. Computational Geometry：lecnotes2
- 7. Computational Geometry：lecnotes20_fixed
- 8. Computational Geometry：lecnotes21
- 9. Computational Geometry：lecnotes23
- 10. Computational Geometry：lecnotes3