Curves And Surfaces For CAGD: A Practical Guide Curves+and+Surfaces+for+CAGD +A+Practical+Guide,+Fifth+Editi

User Manual:

Open the PDF directly: View PDF .
Page Count: 521

Download
Open PDF In Browser	View PDF

Curves and Surfaces
for CAGD
A Practical Guide
Fifth Edition

The Morgan Kaufmann Series in
Computer Grapiiics and Geometric Modeling
Series Editor: Brian A. Barsl^y, University of California, Berl^eley
Curves and Surfaces for CAGD: A Practical
Guide, Fifth Edition
Gerald Farin

Jim Blinn's Corner: A Trip Down the
Graphics Pipeline

Subdivision Methods for Geometric Design:
A Constructive Approach
Joe Warren and Henrik Weimer

Interactive Curves and Surfaces: A
Multimedia Tutorial on CAGD
Alyn Rockwood and Peter Chambers

The Computer Animator's Technical
Handbook
Lynn Pocock and Judson Rosebush

Wavelets for Computer Graphics: Theory
and Applications
Eric J. Stollnitz, Tony D. DeRose,
and David H. Salesin

Computer Animation: Algorithms
Techniques
Rick Parent

and

Principles of Digital Image Synthesis
Andrew S. Glassner

Advanced RenderMan: Creating CGI for
Motion Pictures
Anthony A. Apodaca and Larry Gritz
Curves and Surfaces in Geometric
Theory and Algorithms
Jean GaUier
Andrew Glassner's Notebook:
Computer Graphics
Andrew S. Glassner

Modeling:

Recreational

Warping and Morphing of Graphical Objects
Jonas Gomes, Lucia Darsa, Bruno Costa,
and Luis Velho
Jim Blinn's Corner: Dirty Pixels
Jim BUnn
Rendering with Radiance: The Art and
Science of Lighting Visualization
Greg Ward Larson and Rob Shakespeare
Introduction to Implicit Surfaces
Edited by Jules Bloomenthal

Jim Blinn

Radiosity & Global Illumination
Francois X. Sillion and Claude Puech
Knotty: A B-Spline Visualization Program
Jonathan Yen
User Interface Management
Models and Algorithms
Dan R. Olsen, Jr.

Systems:

Making Them Move: Mechanics, Control,
and Animation of Articulated Figures
Edited by Norman L Badler,
Brian A. Barsky, and David Zeltzer
Geometric and Solid Modeling:
An Introduction
Christoph M. Hoffmann
An Introduction to Splines for Use in
Computer Graphics and Geometric Modeling
Richard H. Bartels, John C. Beatty,
and Brian A. Barsky

Curves and Surfaces
for CAGD
A Practical Guide
Fifth Edition
Gerald Farin
Arizona State University

MORGAN
AN

U

KAUFMANN

IMPRINT

PUBLISHERS

OF A C A D E M I C

PRESS

A Harcourt Science and Technology Company
SAN FRANCISCO

SAN DIEGO

LONDON

SYDNEY

NEW YORK
TOKYO

BOSTON

Executive Director Diane D. Cerra
Publishing Services Manager Scott Norton
Senior Production Editor Cheri Palmer
Assistan Editor Belinda Breyer
Editorial Assistant Mona Buehler
Cover Design Yvo Riezbos Design
Cover Image © 1999 Jeff Greenberg / Rainbow / PictureQuest
Glass, curved rods of steel meet in the center to form a shape similar to the bow
of a boat at the Tampa Aquarium, Tampa, Florida.
Text Design Rebecca Evans & Associates
Composition
Windfall Software, using ZzT^X
Technical Illustration
LM Graphics
Copyeditor
Robert Fiske
Proofreader Erin Milnes
Printer Edwards Brothers
Designations used by companies to distinguish their products are often claimed as
trademarks or registered trademarks. In all instances in which Morgan Kaufmann
Publishers is aware of a claim, the product names appear in initial capital or all capital
letters. Readers, however, should contact the appropriate companies for more complete
information regarding trademarks and registration.
Morgan Kaufmann Publishers
340 Pine Street, Sixth Floor, San Francisco, CA 94104-3205, USA
http://www.mkp.com
ACADEMIC PRESS
A Harcourt Science and Technology Company
525 B Street, Suite 1900, San Diego, CA 92101-4495, USA
http://www.academicpress.com
Academic Press
Harcourt Place, 32 Jamestown Road, London, N W l 7BY, United Kingdom
http://www.academicpress. com
© 2002 by Academic Press
All rights reserved
Printed in the United States of America
06

05

04

03

02

5

4

3

2

1

No part of this publication may be reproduced, stored in a retrieval system, or
transmitted in any form or by any means—electronic, mechanical, photocopying,
recording, or otherwise—without the prior written permission of the publisher.
Library of Congress Control Number: 2001094373
ISBN: 1-55860-737-4
This book is printed on acid-free paper.

To My Parents

This Page Intentionally Left Blank

Contents

Preface

XV

Chapter 1

P. B e z i e r : H o w a Simple System Ws

Chapter 2

Introductory M a t e r i a l

13

2.1

Points and Vectors

13

2.2

Affine Maps

17

2.3

Constructing Affine Maps

20

2.4

Function Spaces

22

2.5

Problems

24

Chapter 3

Linear I n t e r p o l a t i o n
Linear Interpolation
3.1
3.2
Piecewise Linear Interpolation
3.3
3.4

Chapter 4

Menelaos' Theorem

3.5

Blossoms
Barycentric Coordinates in the Plane

3.6

Tessellations

3.7

Triangulations

3.8

Problems

1

25
25
29
30
31
34
37
39
41

The de Casteljau Algorithm

43

4.1

Parabolas

43

4.2

The de Casteljau Algorithm

45

4.3

Some Properties of Bezier Curves

47
vii

viii

Contents
4.4

The Blossom

4.5

Implementation

4.6

Problems

Chapter 5 The Bernstein Form of a Bezier Curve
5.1

Bernstein Polynomials

5.2

Properties of Bezier Curves

5.5

The Derivatives of a Bezier Curve

5.4

Domain Changes and Subdivision

5.5

Composite Bezier Curves

5.6

Blossom and Polar

5.7

The Matrix Form of a Bezier Curve

5.8

Implementation

5.9

Problems

Chapter 6 Bezier Curve Topics

Chapter 7

6.1

Degree Elevation

6.2

Repeated Degree Elevation

6.5

The Variation Diminishing Property

6.4
6.5
6.6
6.7
6.8
6.9
6.10
6.11
6.12

Degree Reduction
Nonparametric Curves
Cross Plots
Integrals
The Bezier Form of a Bezier Curve
The Weierstrass Approximation Theorem
Formulas for Bernstein Polynomials
Implementation
Problems

50
53
54
57
57
60
62
68
71
73
74
75
78

81
81
83
84
85
86
87
88
89
90
91
92
93

Polynomial Curve Constructions

95

7.1

Aitken's Algorithm

7.2

Lagrange Polynomials

7.5

The Vandermonde Approach

95
98
99
101
102
106
107

7.4

Limits of Lagrange Interpolation

7.5

Cubic Hermite Interpolation

7.6

Quintic Hermite Interpolation

7.7

Point-Normal Interpolation

Contents

ix

7.8

Least Squares Approximation

108

7.9

Smoothing Equations

110

7.10

Designing with Bezier Curves

112

7.11

The Newton Form and Forward Differencing

114

7.12

Implementation

116

7.13

Problems

117

Chapter 8 B-Spline Curves

119

8.1

Motivation

120

8.2

B-Spline Segments

122

8.3

B-Spline Curves

126

8.4

Knot Insertion

130

8.5

Degree Elevation

133

8.6

Greville Abscissae

134

8.7

Smoothness

135

8.8

B-Splines

140

8.9

B-Spline Basics

143

8.10

Implementation

144

8.11

Problems

146

Chapter 9 Constructing Spline Curves

147

9.1

Greville Interpolation

147

9.2
9.3
9.4

149
153
154
157

9.8

Least Squares Approximation
Modifying B-Spline Curves
C^ Cubic Spline Interpolation
More End Conditions
Finding a Knot Sequence
The Minimum Property
C^ Piecewise Cubic Interpolation

9.9

Implementation

174

9.10

Problems

176

9.5
9.6
9.7

Chapter 10 W. B o e i i m : D i f f e r e n t i a l G e o m e t r y 1

161
167
169

179

10.1

Parametric Curves and Arc Length

10.2

The Frenet Frame

181

10.3

Moving the Frame

182

10.4

The Osculating Circle

184

179

Contents
10.5
10.6

Nonparametric Curves
Composite Curves

Chapter 11 Geometric Continuity
1
11 2
11 3
11 4
11 5
6
7
11.8

Motivation
The Direct Formulation
The y, V, and p Formulations
C^ Cubic Splines
Interpolating C^ Cubic Splines
Higher-Order Geometric Continuity
Implementation
Problems

Chapter 12 Conic Sections
12.1
12.2
12.5
12.4
12.5
12.6
12.7
12.8
12.9
12.10

Projective Maps of the Real Line
Conies as Rational Quadratics
A de Casteljau Algorithm
Derivatives
The Implicit Form
Two Classic Problems
Classification
Control Vectors
Implementation
Problems

Chapter 13 Rational Bezier and B-Spline Curves
1 3.1 Rational Bezier Curves
1 3.2 The de Casteljau Algorithm
13.3 Derivatives
1 3.4 Osculatory Interpolation
1 3.5 Reparametrization and Degree Elevation
13.6 Control Vectors
13.7 Rational Cubic B-Spline Curves
13.8 Interpolation with Rational Cubics
13.9 Rational B-Splines of Arbitrary Degree
1 3.10 Implementation
13.11 Problems

187
188
191
191
192
194
195
199
200
203
203
205
205
209
214
215
216
219
220
223
224
225
227
227
230
233
234
235
238
238
240
241
242
243

Contents
Chapter 14 Tensor Product Patches
14.1

Bilinear Interpolation

14.2

The Direct de Casteljau Algorithm

1 4.3

The Tensor Product Approach

14.4

Properties

14.5

Degree Elevation

14.6

Derivatives

14.7

Blossoms

14.8

Curves on a Surface

14.9

Normal Vectors

14.10 Twists
14.11 The Matrix Form of a Bezier Patch
14.12 Nonparametric Patches
14.1 3 Problems

Chapter 15 Constructing Polynomial Patches
15.1

Ruled Surfaces

15.2

Coons Patches

15.3

Translational Surfaces

15.4

Tensor Product Interpolation

15.5

Bicubic Hermite Patches

15.6
1 5.7
1 5.8
15.9
15.10
15.11

Least Squares
Finding Parameter Values
Shape Equations
A Problem with Unstructured Data
Implementation
Problems

Chapter 16 Composite Surfaces
1 6.1

Smoothness and Subdivision

1 6.2

Tensor Product B-Spline Surfaces

1 6.3

Twist Estimation

1 6.4

Bicubic Spline Interpolation

1 6.5
16.6
1 6.7

Finding Knot Sequences
Rational Bezier and B-Spline Surfaces
Surfaces of Revolution

xi

245
245
247
250
253
254
255
258
260
261
264
265
266
267

269
269
270
272
274
276
278
281
282
282
283
284
285
285
288
290
293
295
297
299

xii

Contents
16.8
16.9
16.10
16.11

Volume Deformations
CONS and Trimmed Surfaces
Implementation
Problems

301
304
306
308

Chapter 17 Bezier IV-iangles
17.1 The de Casteljau Algorithm
17.2 Triangular Blossoms
17.3 Bernstein Polynomials
17.4 Derivatives
17.5 Subdivision
17.6 Differentiability
17.7 Degree Elevation
17.8 Nonparametric Patches
17.9 The Multivariate Case
17.10 S-Patches
17.11 Implementation
17.12 Problems

309

Chapter 18 Practical Aspects of Bezier Ik'iangies

335
335
337
341
341
343
345
346
347

18.1 Rational Bezier Triangles
18.2 Quadrics
18.3 Interpolation
18.4 Cubic and Quintic Interpolants
18.5 The Clough-Tocher Interpolant
18.6 The Poweli-Sabin Interpolant
18.7 Least Squares
18.8 Problems

Chapter 19 W. Boelim: Differentiai Geometry 11
19.1 Parametric Surfaces and Arc Element
19.2 The Local Frame
19.3 The Curvature of a Surface Curve
19.4 Meusnier's Theorem
19.5 Lines of Curvature
19.6 Gaussian and Mean Curvature
19.7 Euler's Theorem

309
313
315
316
320
323
326
326
328
330
331
332

349
349
352
352
354
355
357
358

Contents

xili

19.8

Dupin's Indicatrix

359

19.9

Asymptotic Lines and Conjugate Directions

360

19.10 Ruled Surfaces and Developabies

361

19.11 Nonparametric Surfaces

363

19.12 Composite Surfaces

364

Chapter 20 Geometric Continuity for Surfaces

367

20.1

Introduction

367

20.2

Triangle-Triangle

20.3

Rectangle-Rectangle

368
372

20.4

Rectangle-Triangle

373

20.5

"Filling in" Rectangular Patches

374

20.6

"Filling in" Triangular Patches

375

20.7

Theoretical Aspects

375

20.8

Problems

376

Chapter 21 Surfaces witii Arbitrary Topology

377

21.1

Recursive Subdivision Curves

377

21.2

Doo-Sabin Surfaces

380

21.5

Catmull-Clark Subdivision

383

21.4

Midpoint Subdivision

386

21.5

Loop Subdivision

387

21.6
21.7

\/3 Subdivision
Interpolating Subdivision Surfaces

21.8
21.9
21.10
21.11

Surface Splines
Triangular Meshes
Decimation
Problems

389
389
392
394
395
398

Chapter 22 Coons Patclies

399

22.1

Coons Patches: Bilinearly Blended

400

22.2

Coons Patches: Partially Bicubically Blended

402

22.3

Coons Patches: Bicubically Blended

404

22.4

Piecewise Coons Surfaces

406

22.5

Two Properties

407

22.6
22.7

Compatibility
Gordon Surfaces

408
410

xiv

Contents
22.8 Boolean Sums
22.9 Triangular Coons Patches
22.10 Problems

412
414
417

Chapter 25 Shape
25.1 Use of Curvature Plots
23.2 Curve and Surface Smoothing
23.3 Surface Interrogation
23.4 Implementation
23.5 Problems

419
419
421
425
427
429

Chapter 24 Evaluation of Some Methods
24.1 Bezier Curves or B-Spline Curves?
24.2 Spline Curves or B-Spline Curves?
24.3 The Monomial or the Bezier Form?
24.4 The B-Spline or the Hermite Form?
24.5 Triangular or Rectangular Patches?

431
431
431
432
435
435

Appendix A Quick Reference of Curve and Surface Terms

437

Appendix B List of Programs

445

Appendix C Notation

447

References

449

Index

491

Preface

Computer aided geometric design (CAGD) is a discipline dealing with computational aspects of geometric objects. It is best explained by a brief historical
sketch.^
Renaissance naval architects in Italy were the first to use drafting techniques
that involved conic sections. Prior to that, ships were built "hands on" without
any mathematics being involved. These design techniques were refined through
the centuries, culminating in the use of splines—wooden beams that were bent
into optimal shapes. In the beginning of the twentieth century, airplanes made
their first appearance. Their design (or rather, the design of the outside fuselage)
was streamlined by the use of conic sections, as pioneered by R. Liming [390]. He
devised methods that went beyond traditional drafting with conies—for the first
time, certain conic coefficients could be used to define a shape—thus numbers
could be used to replace blueprints!
The automobile, one of the defining cultural icons of the twentieth century,
also needed new design approaches as mass production started. In the late 1950s,
hardware became available that allowed the machining of 3D shapes out of
blocks of wood or steel.^ These shapes could then be used as stamps and dies for
products such as the hood of a car. The bottleneck in this production method was
soon found to be the lack of adequate software. In order to machine a shape using
a computer, it became necessary to produce a computer-compatible description
of that shape. The most promising description method was soon identified to
be in terms of parametric surfaces. An example of this approach is provided in

1 For more details, see [204].
2 A process that is now called CAM for computer aided manufacturing.
XV

xvi

Preface
Color Plates I and III: Color Plate I shows the actual hood of a car; Color Plate
III shows how it is represented internally as a collection of parametric surfaces.
The major breakthroughs in CAGD were the theory of Bezier curves and
surfaces, later combined with B-spline methods. Bezier curves and surfaces were
independently developed by P. de Casteljau at Citroen and by P. Bezier at Renault.
De Casteljau's development, slightly earlier than Bezier's, was never published,
and so the whole theory of polynomial curves and surfaces in Bernstein form
now bears Bezier's name. CAGD became a discipline in its own right after the
1974 conference at the University of Utah (see Barnhill and Riesenfeld [34]).
Several other disciplines have emerged and interacted with CAGD. Computational geometry is concerned with the analysis of geometric algorithms. An
example would be finding a bound on the time it takes to triangulate a set of
points. Knowledge of such bounds allows a comparison and evaluation of different algorithms. The literature includes Prepata and Shamos [497] and de Berg
et al. [135]. Ironically, another book with the term computational geometry in
it is the one by Faux and Pratt [228]. It was a very influential text, but today, it
would be classified as a CAGD text.
Another related discipline is solid modeling. It is concerned with the representation of objects that are enclosed by an assembly of surfaces, mostly very
elementary ones such as planes, cylinders, or tori. The literature includes Hoffmann [327] and Mantyla [416]. CAGD has also influenced fields such as medical
imaging, geographic information systems, computer gaming, and scientific visualization. It should go without saying that computer graphics is one of the earliest
and most important applications of CAGD; see [238] or [9].
For this fifth edition, the most notable addition is a chapter on subdivision
surfaces that were of academic interest at best when the first edition appeared in
1988. A recent special issue of the journal Computer Aided Geometric Design
highlights some of the new developments; see [393], [432], [470], [579], [598],
and [631]. Other new topics include triangle meshes, more in-depth treatment
of least squares techniques, and pervasive use of the blossoming principle.
Each chapter is concluded by a set of problems. They come in three categories:
simpler exercises at the beginning of each Problem section, harder problems
marked by asterisks, and programming problems marked by "P." Many of these
programming problems use data on the Web site. Students should thus get a better
feeling for "real" situations. In teaching this material, it is essential that students
have access to computing and graphics facilities; practical experience greatly
helps the understanding and appreciation of what might otherwise remain dry
theory.
The C programs on the Web site are my implementations of some (but not
all) of the most important methods described here. The programs were tested for
many examples, but they are not meant to be "industrial strength." In general.

Preface

^*'^K^

xvii

no checks are made for consistency or correctness of input data. Also, modularity
was valued higher than efficiency. The programs are in C but with non-C users in
mind—in particular, all modules should be easily translatable into FORTRAN.
The Web page for the book is www.mkp.com/cagd5e. This page includes C
programs, data sets, and errata.
As for all previous editions, sincere thanks go to Dianne Hansford for help
and advice with all aspects of the book.

This Page Intentionally Left Blank

p. Bezier

M

K

How a Simple System Was Born

I n order to solve CAD/CAM mathematical problems, many solutions have been
offered, each adapted to specific matters. Most of the systems have been invented
by mathematicians, but UNISURF v^as developed by mechanical engineers from
the automotive industry w^ho v^ere familiar v^ith parts mainly described by lines
and circles. Fillets and other blending auxiliary surfaces v^ere scantly defined;
their final shape w^as left to the skill and experience of patternmakers and diesetters.
Around 1960, designers of stamped parts, that is, car-body panels, used french
curves and sw^eeps, but in fact, the final standard was the "master model," the
shape of w^hich, for many valid reasons, could not coincide w^ith the curves
traced on the drawling board. This resulted in discussions, arguments, haggling,
retouches, expenses, and delay.
Obviously, no significant improvement could be expected so long as a method
w^as not devised that could prove an accurate, complete, and undisputable
definition of freeform shapes.
Computing and numerical control (NC), at that time, had made great progress,
and it w^as certain that only numbers, transmitted from drawing office to tool
drawing office, manufacturing, patternshop, and inspection could provide an
answer; of course, drawings would remain necessary, but they would only be
explanatory, their accuracy having no importance, and numbers being the only
and final definition.
Certainly, no system could be devised without the help of mathematics—yet
designers, who would be in charge of operating it, had a good knowledge of
geometry, especially descriptive geometry, but no basic training in algebra or
analysis.
1

Chapter 1 P. Bezier: How a Simple System Was Born

Figure 1.1 An arc of a hand-drawn curve is approximated by a part of a template.

In France, at that time, very httle vv^as know^n about the work performed in
the American aircraft industry; the papers from James Ferguson v^ere not much
displayed before 1964; Citroen was secretive about the results obtained by Paul
de Casteljau, and the famous technical report MAC-TR-41 (by S. A. Coons) did
not appear before 1967; The w^orks of W. Gordon and R. Riesenfeld w^ere printed
in 1974.
At the beginning, the idea of UNISURF w^as oriented tow^ard geometry rather
than analysis, but v^ith the idea that every datum should be exclusively expressed
by numbers.
For instance, an arc of a curve could be represented (Figure 1.1) by the
coordinates, cartesian, of course, of its limit points (A and B), together w^ith
their curvilinear abscissae, related w^ith a grid traced on the edge.
The shape of the middle line of a sw^eep is a cube, if its cross section is constant,
its matter is homogeneous, and neglecting the effect of friction on the tracing
cloth. How^ever, it is difficult to take into account the length betw^een endpoints;
moreover, the curves employed for softw^are for NC machine tools, that is, 2D
milling machines, v^ere lines and circles, and sometimes, parabolas. Hence, a
spline shape should be divided and subdivided into small arcs of circles put end
to end.
To transform an arc of circle into a portion of an ellipse, imagine (Figure 1.2) a
square frame containing tw^o sets of strings, v^hose intersections w^ould be located
on an arc of a circle; the frame sides being hinged, the square is transformed into
a diamond (Figure 1.3), and the circle becomes an arc of an ellipse, which would

p. Bezier: How a Simple System Was Born

• [ ^f

f

f

»[ 1 1 1
^ 1 1 1

•[II

1

f

1
1

JT

Figure 1.2

*

-f i^

c

1 ^^"^^^^^^^^^
Jkr
\
1 M

1

UK
•1 4*

-f

3

1

lit

1 II
*

*

*

1
* 1*

A circular arc is obtained by connecting the points in this rectangular grid.

Figure 1.3 If the frame from the previous figure is sheared, an arc of an ellipse is obtained.
be entirely defined as soon as the coordinates of points A, B, and C were known;
if the hinged sides of the frame were replaced by pantographs (Figure 1.4), the
diamond would become a parallelogram, and the definition of the arc of ellipse
still results from the coordinates of the three points A, B, and C (Figure 1.5).

4

Chapter 1 P. Bezier: How a Simple System Was Born

Figure 1 A Pantograph construction of an arc of an ellipse.
Of course, this idea was not realistic, but it was easily replaced by the computation of coordinates of successive points of the curve. Harmonic functions were
available with the help of analog computers, which were widely used at that time
and gave excellent results.
But employing only arcs of ellipses limited by conjugate diameters was far too
restrictive, and a more flexible definition was required.
Another idea came from the practice of a speaker projecting, with a flashlight,
a small sign, cross, or arrow, onto a screen displaying a figure printed on a
slide. Replacing the arrow with a curve and recording the exact location and
orientation of the flashlight (Figure 1.6) would define the image of the curve
projected on the wall of the drawing office. One could even imagine having a

p. Bezier: How a Simple System Was Born

5

Figure 1.5 A "control polygon" for an arc of an ellipse.

^c2^Z

Figure 1.6 A projector producing a "template curve" on the drawing of an object.

variety of slides, each of which would bear a specific curve: circle, parabola,
astroid, and so on.
Of course, this was not a realistic idea because the focal plane of the zoom
would seldom be square to the axis—an optician's nightmare! But the principle could be translated, via projective geometry and matrix computation, into
cartesian coordinates.

Chapter 1 P. Bezier: How a Simple System Was Born

Figure 1.7 Two imaginary projections of a car.

At that time, designers defined the shape of a car body by cross sections located
100 mm apart, and sometimes less. The advantage was that, from a drawing, one
could derive templates for adjusting a clay model, a master, or a stamping tool.
The drawback was that a stylist does not define a shape by cross sections but with
so-called character lines, which seldom are plane curves. Hence, a good system
should be able to manipulate and define directly "space curves" or "freeform
curves." Of course, one could imagine working alternately (Figure 1.7) on two
projections of a space curve, but it is very unlikely that a stylist would accept
such a solution.
Theoretically as least, a space curve could be expressed by a sweep having a
circular section, constrained by springs or counterweights (Figure 1.8), but this
would prove quite impractical.
Would it not be best to revert to the basic idea of a frame .^ But instead of a
curve inscribed in a square, it would be located in a cube (Figure 1.9) that could
become any parallelepiped (Figure 1.10) by a linear transformation that is easy to
compute. The first idea was to choose a basic curve that would be the intersection
of two circular cylinders; the parallelepiped would be defined by points O, X, Y,
and Z, but it is more practical to put the basic vectors end to end so as to obtain
a polygon OMNB (Figure 1.10), which defines directly the endpoint B and its

p. Bezier: How a Simple System Was Born

Figure 1.8

A curve held by springs.

Figure 1.9

A curve defined inside a cube.

7

8

Chapter 1 P. Bezier: How a Simple System Was Born

Y
Figure 1.10

A curve defined inside a parallelepiped.

tangent NB. Of course, points O, M, N, and B need not be coplanar and can
define a space curve.
Polygons v^ith three legs can define quite a variety of curves, but in order
to increase it, w^e can imagine making use of cubes and hypercubes of any order
(Figure 1.11) and the relevant polygons (Figure 1.13 and see Figure 4.4 in Section
4.3).
At that moment, it became necessary to do away with harmonic functions and
revert to polynomials. This was even more desirable since digital computers were
gradually replacing analog computers. The polynomial functions were chosen
according to the properties that were considered best: tangency, curvature, and
the like. Later, it was discovered that they could be regarded as sums of Bernstein's
functions.
When it was suggested that these curves could replace sweeps and french
curves, most stylists objected that they had invented their own templates and
would not change. It was solemnly promised that their "secret" curves should
be translated in secret listings, and buried in the most secret part of the memory
of the computer, and that nobody but they would keep the key of the vaulted
cellar. In fact, the standard curves were flexible enough, and secret curves were
soon forgotten. Designers and draftsmen easily understood the polygons and
their relation with the shape of the corresponding curves.
In the traditional process of body engineering, a set of curves is carved in a
3D model, between which interpolation is left to the experience of highly skilled
patternmakers. However, in order to obtain a satisfactory numerical definition,
the surface had to be totally expressed with numbers.

p. Bezier: How a Simple System Was Born

Figure 1.11

Higher-order curves can be defined inside higher-dimensional cubes.
At that time, around 1960, very Uttle, if anything, had been published about
biparametric patches. The basic idea of UNISURF came from a comparison v^ith
a process often used in foundries to obtain a core. Sand is compacted in a box
(Figure 1.12), and the shape of the upper surface of the core is obtained by
scraping off the surplus w^ith a timber plank cut as a template. Of course, a
shape obtained by such a method is relatively simple since the shape of the plank
is constant and that of the box edges is generally simple. To make the system
more flexible, one might w^ish to change the shape of the template as it moves.
In fact, this takes us back to a very old, and sometimes forgotten, definition of
a surface: it is the locus of a curve that is at the same time moved and distorted.
About 1970, a Dutch laboratory sculptured blocks of styrofoam w^ith a flexible
electrically heated strip of steel, the shape of w^hich w^as controlled by the flexion
torque imposed on its extremities.
This process could not produce a large variety of shapes, but the principle
could be translated into a mathematical solution. The guiding edges of the box
are similar to the curves AB and CD of Figure 1.13, v^hich can be considered
directrices of a surface, defined by their characteristic polygon. If a curve such as
EF is generatrix, defined by its own polygon, the ends of which run along lines AB
and CD, and the intermediate vertices of the polygon are on curves GH and JK,
the surface ABDC is known as soon as the four polygons are defined. Connecting
the corresponding vertices of the polygons defines the "characteristic net" of

10

Chapter 1 P. Bezier: How a Simple System Was Born

Figure 1.12 A surface is being obtained by scraping off excess material with wooden templates.

Figure 1.13 The characteristic net of a surface.

the patch, which plays, regarding the surface, the same part as a polygon of a
curve. Hence, the cartesian coordinates of the points of the patch are computed
according to the values of two parameters.
After expressing this basic idea, a good many problems remained to be solved:
choosing adequate functions, blending curves and patches, dealing with degener-

p. Bezier: How a Simple" System Was Born

11

ate patches, to name only a few. The solutions were a matter of relatively simple
mathematics, the basic principle remaining untouched.
So, a system has been progressively created. If we consider the way an initial
idea evolved, we observe that the first solution—parallelogram, pantograph—
is the result of an education oriented toward kinematics, the conception of
mechanisms. Next appeared geometry and optics, which very likely came from
some training in the army, when geometry, cosmography, and topography played
an important part. Then reflexion was oriented toward analysis, parametric
spaces, and finally, data processing, because a theory, as convenient as it may
look, must not impose too heavy a task to the computer and must be easily
understood, at least in its principle, by the operators.
The various steps of this conception have a point in common: each idea must
be related with the principle on a material system, however simple and primitive
it may look, on which a variable solution could be based.
Engineers define what is to be done and how it could be done; they not only
describe the goal, they lead the way toward it.
Before looking any deeper into this subject, it should be observed that elementary geometry played a major part, and it should not gradually disappear from
the courses of a mechanical engineer. Each idea, each hypothesis was expressed
by a figure, or a sketch, representing a mechanism. It would have been extremely
difficult to build a purely mental image of a somewhat elaborate system without
the help of pencil and paper. Let us consider, for instance, Figures 1.9 and 1.11;
they are equivalent to equations (5.6) and (14.6) in the subsequent chapters.
Evidently, these formulas, conveniently arranged, are best suited to express data
given to a computer, but most people would better understand a simple figure
than the equivalent algebraic expression.
Napoleon said: "A short sketch is better than a long report."
Which is the part played by experience, by theory, and by imagination in
the creation of a system? There is no definite answer to such a query. The
importance of experience and of theoretical knowledge is not always clearly
perceived. Imagination seems a gift, a godsend, or the result of a beneficial
heredity; but is not, in fact, imagination the result of the maturation of the
knowledge gained during education and professional practice? Is it not born from
facts apparently forgotten, stored in the dungeon of a distant part of memory, and
suddenly remembered when circumstances call them back? Is not imagination
based, partly, on the ability to connect notions that, at first sight, look quite
unrelated, such as mechanics, electronics, optics, foundry, data processing, to
catch barely seen analogies, as Alice in Wonderland, to go "through the mirror"?
Will, someday, psychologists be able to detect in man such a gift that would be
applicable to science and technology? Has it a relation with the sense of humor
that can detect unexpected relations between facts that look quite unconnected?
Shall we learn how to develop it? Will it forever remain a gift, devoted by

12

Chapter 1 P. Bezier: How a Simple System Was Born
pure chance to some people, whereas for others carefulness and cold blood
prevail?
It is important that, sometimes, "sensible" men give free rein to imaginative
people. "I succeeded," said Henry Ford, "because I let some fools try what wise
people had advised me not to let them try."

Introductory Material

2.1 Points and Vectors
When a designer or stylist works on an object, he or she does not think of that
object in very mathematical terms. A point on the object would not be thought of
as a triple of coordinates, but rather in functional terms: as a corner, the midpoint
between two other points, and so on. The objective of this book, however, is to
discuss objects that are defined in mathematical terms, the language that lends
itself best to computer implementations. As a first step toward a mathematical
description of an object, one therefore defines a coordinate system in which it
will be described analytically.
The space in which we describe our object does not possess a preferred
coordinate system—we have to define one ourselves. Many such systems could
be picked (and some will certainly be more practical than others). But whichever
one we choose, it should not affect any properties of the object itself. Our interest
is in the object and not in its relationship to some arbitrary coordinate system.
Therefore, the methods we develop must be independent of a particular choice
of a coordinate system. We say that those methods must be coordinate-free or
coordinate-independent}
We stress the concept of coordinate-free methods throughout this book. It
motivates the strict distinction between points and vectors as discussed next.
(For more details on this topic, see R. Goldman [262].)

1 More mathematically, the geometry of this book is affine geometry. The objects that we
will consider "live" in affine spaces, not in linear spaces.
13

14

Chapter 2 Introductory Material

Figure 2.1

Points and vectors: vectors are not affected by translations.
We shall denote points^ elements of three-dimensional (or 3D) euclidean (or
point) space E^, by lowercase boldface letters such as a, b, and so on. (The term
euclidean space is used here because it is a relatively familiar term to most people.
More correctly, v^e should have used the term affine space.) A point identifies a
location, often relative to other objects. Examples are the midpoint of a straight
line segment or the center of gravity of a physical object.
The same notation (low^ercase boldface) w^ill be used for vectors^ elements of
3D linear (or vector) space M^. If we represent points or vectors by coordinates
relative to some coordinate system, we shall adopt the convention of writing
them as coordinate columns.
Although both points and vectors are described by triples of real numbers,
we emphasize that there is a clear distinction between them: for any two points
a and b, there is a unique vector v that points from a to b. It is computed by
componentwise subtraction:
v= b-a;

a,bGE^v€Rl

On the other hand, given a vector v, there are infinitely many pairs of points a,
b such that v = b — a. For if a, b is one such pair and if w is an arbitrary vector,
then a + w, b + w is another such pair since v = (b + w) — (a + w) also. Figure
2.1 illustrates this fact.
Assigning the point a + w to every point a G E^ is called a translation^ and the
above asserts that vectors are invariant under translations while points are not.
Elements of point space E^ can be subtracted from each other—this operation
yields a vector. They cannot be added—this operation is not defined for points.
(It is defined for vectors.) Figure 2.2 gives an example.

2.1 Points and Vectors

15

*.'--.-,
Figure 2.2 Addition of points: this is not a well-defined operation, since different coordinate systems
would produce different "solutions."
However, additionlike operations are defined for points: they are barycentric
combinations? These are weighted sums of points where the weights sum to one:
b= ^

Qfyby;

by G E-^, QfQ H

h Qf^ = 1.

(2.1)

j=0

At first glance, this looks like an undefined summation of points, but we can
rewrite (2.1) as
n

b = bo + ^Qfy(b^—bo),
which is clearly the sum of a point and a vector.
An example of a barycentric combination is the centroid g of a triangle with
vertices a, b, c, given by
1
1,
1
g= -a + -b + -c.
The term barycentric combination is derived from "barycenter," meaning
"center of gravity." The origin of this formulation is in physics: if the by are
centers of gravity of objects with masses my, then their center of gravity b is
located at b = X! ^/by/ Yl ^j ^^^ has the combined mass ^ my. (If some of the
mj are negative, the notion of electric charges may provide a better analogy;
see Coxeter [130], p. 214.) Since a common factor in the mj is immaterial for
2 They are also called affine combinations.

16

Chapter 2 Introductory Material

Figure 2.3 Convex hulls: a point set (a polygon) and its convex hull, shown shaded.
the determination of the center of gravity, we may normalize them by requiring
An important special case of barycentric combinations are the convex combinations. These are barycentric combinations v^here the coefficients ofy, in addition
to summing to one, are also nonnegative. A convex combination of points is alw^ays "inside" those points, w^hich is an observation that leads to the definition
of the convex hull of a point set: this is the set that is formed by all convex combinations of a point set. Figure 2.3 gives an example (see also Problems). More
intuitively, the convex hull of a set is formed as follows: for a 2D set, imagine a
string that is loosely circumscribed around the set, with nails driven through the
points in the set. Now pull the string tight—it will become the boundary of the
convex hull.
The convex hull of a point set is a convex set. Such a set is characterized by
the following: for any two points in the set, the straight line connecting them is
also contained in the set. Examples are ellipses or parallelograms. It is an easy
exercise to verify that affine maps (see next section) preserve convexity.
Let us return to barycentric combinations, which generate points from points.
If we want to generate a vector from a set of points, we may write
n

7=0

where we have a new restriction on the coefficients: now we must demand that
the a, sum to zero.

2.2 AffineMaps

17

If we are given an equation of the form

and a is supposed to be a point, then we must be able to split the sum into three
groups:

a= E ^A+ E '^A+ E ^AS/6y=l

E^y=0

remaining ^s

Then the by in the first sum are points, and those in the second sum may be
interpreted as either points or vectors. The by in the third sum are vectors.
Whereas the second and third sums may be empty, the first one must contain
at least one term.
The interplay betv^een points and vectors is unusual at first. Later, it w^ill turn
out to be of invaluable theoretical and practical help. For example, w^e can perform quick type checking v^hen wt derive formulas. If the point coefficients fail
to add up to one or zero—depending on the context—w^e know^ that something
has gone w^rong. In a more formal w^ay, T. DeRose has developed the concept of
"geometric programming," a graphics language that automatically performs type
checks [1601, [161]. R. Goldman's article [262] treats the validity of point/vector
operations in more detail.

2.2 Affine Maps
Most of the transformations that are used to position or scale an object in a
computer graphics or CAD environment are affine maps. (More complicated,
so-called projective maps are discussed in Chapter 12.) The term affine map is
due to L. Euler; affine maps were first studied systematically by R Moebius [429].
The fundamental operation for points is the barycentric combination. We
will thus base the definition of an affine map on the notion of barycentric
combinations. A map O that maps E^ into itself is called an affine map if it
leaves barycentric combinations invariant. So if
X= ^

Qfyay;

Y^ aj = 1, x, ay G E^

and 0 is an affine map, then also
p^, we
simply compute it; then O^^ is the desired affine map.
In general, the two point sets are not related by an affine map; this procedure
will still produce an affine map that approximately maps the two point sets.

22

Chapter 2 Introductory Material

Qx

o

o

Figure 2.5 Norm ellipses: a point set and an associated ellipse.

This method works even if the number of q/ does not equal that of the p/.
An obvious generalization v^orks in 3D. Applications are in image registration,
where two images (typically represented by point sets) have to be mapped to each
other.

2.4 Function Spaces
This section contains material that will later simplify our work by allowing very
concise notation. Although we shall try to develop our material with an emphasis
on geometric concepts, it will sometimes simplify our work considerably if we
can resort to some elementary topics from functional analysis. Good references
are the books by Davis [133] and de Boor [138].
Let C[a^ h\ be the set of all real-valued continuous functions defined over the
interval \a^ h\ of the real axis. We can define addition and multiplication by a
constant for elements /^,g € C[a, h\ by setting (af + Pg)(t) = af(t) + Pg(t) for all
t 6 [a, b]. With these definitions, we can easily show that C[a^ b] forms a linear
space over the reals. The same is true for the sets C^[^,fe],the sets of all realvalued functions defined over [a, b] that are fe-times continuously differentiable.
Furthermore, for every k, C^^^ is a subspace of C^.
We say that n functions / ^ i , . . . , /^ G C[a^ b] are linearly independent if ^ Cifi =
0 for all t e [a^ b] implies ci = .. . = c^ = 0.
We mention some subspaces of C[a, b] that will be of interest later. The spaces
V^ of all polynomials of degree n are
p'^it) = ao-\- ait + ait^ + • • • + aj""-,

t e [a, b].

2.4 Function Spaces

23

For fixed n^ the dimension of V^\sn-\-1: each p^ e V^ is determined uniquely by
the w + 1 coefficients UQ^, ..., a^. These can be interpreted as a vector in {n + 1)dimensional finear space R " + \ which has dimension w + 1. We can also name
a basis for V^: the monomials \,t,t^^. .. ,t^ are w + 1 finearly independent
functions and thus form a basis.
Another interesting class of subspaces of C[^, b] is given by piecewise linear
functions: let ^ = ^Q < ^l < ' * * < ^w = ^ be a partition of the interval [a, b\ A
continuous function that is linear on each subinterval [^/, tij^^ is called a piecewise
linear function. Over a fixed partition of [a., b\ the piecev^ise linear functions form
a linear function space. A basis for this space is given by the hat functions: a hat
function Hi{t) is a piecev^ise linear function v^ith H^(^^) = 1 and H^(^y) = 0 if / /= /•
A piecew^ise linear function f v^ith f{tj) = ^ can alw^ays be v^ritten as
n

f{t) = Y.fi^i^^)'
Figure 2.6 gives an example.
We v^ill also consider linear operators that assign a function Af to a given
function f. An operator ^ : C[a, b] -^ C[a, b] is called linear if it leaves linear
combinations invariant:
A(af + Pg) = aAf + pAg;

a, ^ G R.

Figure 2.6 Hat functions: the piecewise linear function f can be written as /" = HQ + 3Hi + 2H2

24

Chapter 2 Introductory Material
An example is given by the derivative operator that assigns the derivative f^ to a
given function f: Af = f.

2.5 Problems
1 Of all affine maps, shears seem to be the least familiar to most people.^
Construct a matrix that maps the unit square with points (0,0), (1,0), (1,1),
(0,1) to the parallelogram v^ith image points (0,0), (1,0), (2,1), (1,1).
* 2 We have seen that affine maps leave the ratio of three collinear points
constant (i.e., they are ratio preserving). Show that the converse is also
true: every ratio-preserving map is affine.
* 3 Show that the n-\-l functions fi(t) = f; / = 0 , . . . , « are linearly independent.
P1 Fix two distinct points a, b on the x-axis. Let a third point x trace out all the
X-axis. For each location of x, plot the value of the function ratio(a, x, b),
thus obtaining a graph of the ratio function.

4

Recall that Figure 2.4 illustrates a shear.

Linear Interpolation

M ost of the computations that we use in CAGD may be broken down into
seemingly trivial steps—sequences of linear interpolations. It is therefore important to understand the properties of these basic building blocks. This chapter
explores those properties and introduces a related concept, called blossoms.

5.1 Linear Interpolation
Let a, b be two distinct points in E^. The set of all points x e E^ of the form
x = x(t) = (l-t)ei-^th;

t eR

(3.1)

is called the straight line through a and b. Any three (or more) points on a straight
line are said to be collinear.
For ^ = 0, the straight line passes through a and for ^ = 1, it passes through b.
For 0 < ^ < 1, the point x is between a and b, whereas for all other values of t it
is outside; see Figure 3.1.
Equation (3.1) represents x as a barycentric combination of two points in
E-^. The same barycentric combination holds for the three points 0, t, 1 in E^:
t = (l — t)-0-\-t'l.Sotis
related to 0 and 1 by the same barycentric combination
that relates x to a and b. Hence, by the definition of affine maps, the three points
a, X, b are an affine map of the three ID points 0, ^, 1! Thus linear interpolation
is an affine map of the real line onto a straight line in E^.^
1 Strictly speaking, we should therefore use the term affine interpolation instead of linear
interpolation. We use linear interpolation because its use is so widespread.

25

26

Chapter 3 Linear Interpolation

0
Figure 3.1

t

1

Linear interpolation: two points a,b define a straight hne through them. The point t in
the domain is mapped to the point x in the range.
It is now almost a tautology when we state: linear interpolation is affinely
invariant. Written as a formula: if 0 is an affine map
onto itself, and (3.1)
holds, then also
2i + t^h.

(3.2)

Since affine maps may be applied to vectors as well as to points, it makes sense
to ask what linear interpolation will do to vector arguments. These vectors "live"
in ID domain space, and will be denoted by v.
If c and d are two ID points in the domain, they define a vector v by setting
y — d — c. The corresponding vector \(v) in the range is then defined as
\{v)=\{d)-\{c),

(3.3)

Figure 3.2 illustrates. For the special case of v being the ID zero vector ? = 0,
we have
1(0) = 0.^

(3.4)

Closely related to linear interpolation is the concept of barycentric coordinates^ due to Moebius [429]. Let a, x, b be three coUinear points in E^:
x = a a + )Sb;
2

Here, 0 denotes the zero vector.

a + jS = l .

(3.5)

3.1 Linear Interpolation

0
Figure 3.2

1

27

-o«-

1

Linear interpolation: the vector 1 in the domain is mapped to the vector 1(1) in the range.
Then a and ^ are called harycentric coordinates of x with respect to a and b.
Note that by our previous definitions, x is a barycentric combination of a and b.
The connection between barycentric coordinates and linear interpolation is
obvious: we have a = \ — t and ^ — t. This shows, by the way, that barycentric
coordinates do not always have to be positive: for t ^ [0,1], either a or ^ is
negative. For any three collinear points a, b, c, the barycentric coordinates of b
with respect to a and c are given by
voli(b,c)
a = —-^
,
voli(a,c)
voli(a,b)
voli(a, c ) '
where vol^ denotes the one-dimensional volume, which is the signed distance
between two points. Barycentric coordinates are not only defined on a straight
line, but also on a plane. Section 3.5 has more details.
Another important concept in this context is that of ratios. The ratio of three
collinear points a, b, c is defined by
ratio(a, b, c) =

voli(a,b)
voli(b,c)

(3.6)

If a and ^ are barycentric coordinates of b with respect to a and c, it follows that
ratio(a, b, c) =

(3.7)

28

Chapter 3 Linear Interpolation
The barycentric coordinates of a point do not change under affine maps, and
neither does their quotient. Thus the ratio of three cohinear points is not affected
by affine transformations. So if (3.7) holds, then also
ratio(4)a, Ob, Oc) = - ,
a

(3.8)

where O is an affine map. This property may be used to compute ratios efficiently.
Instead of using square roots to compute the distances between points a, x,
and b, we would project them onto one of the coordinate axes and then use
simple differences of their x- or ^-coordinates.^ This shortcut works since parallel
projection is an affine map!
Equation (3.8) states that affine maps are ratio preserving. This property may
be used to define affine maps. Every map that takes straight lines to straight lines
and is ratio preserving is an affine map.
The concept of ratio preservation may be used to derive another useful property of linear interpolation. We have defined the straight line segment [a, b] to be
the affine image of the unit interval [0,1], but we can also view that straight line
segment as the affine image of any interval [a,, b]. The interval [^, b] may itself be
obtained by an affine map from the interval [0,1] or vice versa. With t € [0,1]
and u G [a, b\ that map is given by t = (u — a)/(b — a). The interpolated point
on the straight line is now given by both
x(t) = (1 - t)a + ^b
and
x(u) =

a+
b.
b—a
b—a

(3.9)

Since a^ w, b and 0, ^, 1 are in the same ratio as the triple a, x, b, we have shown
that linear interpolation is invariant under affine domain transformations. By
affine domain transformation, we simply mean an affine map of the real fine
onto itself. The parameter t is sometimes called a local parameter of the interval
A more general way to express this is by saying that any barycentric combination of three domain points r, s, t (not necessarily involving any interval
endpoints) carries over to the corresponding range points:
s = {l-ot)r
3

+ at=^ x(s) = (1 - Qf)x(r) + a x ( 0 .

(3.10)

But be sure to avoid projection onto the x-axis if the three points are parallel to the y-axis!

3.2 Piecewise Linear Interpolation

29

A concluding remark: we have demonstrated the interplay between the two
concepts of linear interpolation and ratios. In this book, we will often describe
methods by saying that points have to be collinear and must be in a given
ratio. This is the geometric (descriptive) equivalent of the algebraic (algorithmic)
statement that one of the three points may be obtained by linear interpolation
from the other two.

5.2 Piecewise Linear Interpolation
Let b o , . . . ,b„ G E^ form a polygon B. This polygon consists of a sequence of
straight line segments, each interpolating to a pair of points b^, b^_^i. It is therefore
also called the piecewise linear interpolant VC to the points b^. If the points b^
lie on a curve c, then B is said to be a piecewise linear interpolant to c, and we
write

B = H:C.

(3.11)

One of the important properties of piecewise linear interpolation is affine
invariance. If the curve c is mapped onto a curve Oc by an affine map O, then
the piecewise linear interpolant to Oc is the affine map of the original piecewise
linear interpolant:
P £ O c = ^ VLz,

(3.12)

Another property is the variation diminishing property. Consider a continuous
curve c, a piecewise linear interpolant VC c, and an arbitrary plane. Let cross c be
the number of crossings that the curve c has with this plane, and let cross(7X c)
be the number of crossings that the piecewise linear interpolant has with this
plane. (Special cases may arise; see Section 3.8.) Then we always have
cross(P£ c) < cross c.

(3.13)

This property follows from a simple observation: consider two points b/, hj^i.
The straight line segment through them can cross a given plane at one point at
most, whereas the curve segment from c that connects them may cross the same
plane in many arbitrary points. The variation diminishing property is illustrated
in Figure 3.3.

30

Chapter 3 Linear Interpolation

Figure 3.3 The variation diminishing property: a piecewise hnear interpolant to a curve has no more
intersections with any plane than the curve itself.

5.5 Menelaos' Theorem
We use the concept of piecewise linear interpolation to prove one of the most
important geometric theorems for the theory of CAGD: Menelaos' theorem.
This theorem can be used for the proof of many constructive algorithms, and
its importance v^as already realized by de Casteljau [146] and W. Boehm [67].
Referring to Figure 3.4, w^e define
b[0, ^] = (1 - ^)bo + thi,
b[s, 0] = (1 - s)bo + sbi,
b [ l , ^ ] = ( l - ^ ) b i + fb2,
b[s,l] = ( l - s ) b i + sb2.
Let us further define tv^o points
b[s, ^] = (1 - Ob[s, 0] + th[s, 1] and
b[^ s] = (1 - s)b[0, t] + sh[t, 1].

(3.14)

Menelaos' theorem now states that these points are identical:
b[5,^] = b [ ^ 4

(3.15)

For a proof, we simply verify that
b[s, t] = h[t, s\ = {l- t){l - s)bo + [(1 - t)s + t{l - s)]bi + sth2. (3.16)
Some interesting special cases are given by b[0, 0] = bo or by b[0,1] = b^.
Equation (3.15) is a "CAGD version" of the original Menelaos' theorem,
which may be stated as (see Coxeter [130]):

3.4 Blossoms

31

Figure 3.4 Menelaos' theorem: the point b[s, t] may be obtained from Hnear interpolation at t or
at s.
ratio(b[s, l ] , b [ l , 4 b i ) • ratio(bi,b[0,4b[s, 0])ratio(b[s, 0],b[s, 4 b [ s , 1]) = - 1 .

(3.17)

The proof of (3.17) is a direct consequence of (3.15). Note the ordering of points
in the second ratio! Menelaos' theorem is closely related to Ceva's, which is given
in Section 3.5.

5.4 Blossoms
The bivariate function b[^i, ti] from (3.16) will be very important for the remainder of this book. Functions of that type are called blossoms. Before we introduce
the general concept, we will further explore properties of (3.16).
The first property is called symmetry. It states that the order of the blossom
arguments does not matter—which is exactly Menelaos' theorem.
In Section 3.1, we saw that linear interpolation carries domain relationships
over to corresponding range relationships; see (3.10). Since blossoms are evaluated using linear interpolations, we now have: if the first argument ti of a blossom
is a barycentric combination of two (or more) ID points r and s, we may compute
the blossom values for each argument and then form their barycentric combination:
h[ar + ^s, t2] = ah[r, ti] + ^h{s, til

a + ^6 = 1.

(3.18)

Equation (3.18) states that the blossom b is affine with respect to its first argument, but it is affine for the second one as well because of the symmetry property.
This is the reason the blossom is called multiaffine—the second of its main properties.
For a third property, we study what happens if both blossom arguments are
equal: ti = t2 = t. Then the expression b[^, t\ denotes a point that depends on

32

Chapter 3 Linear Interpolation
one variable t—^thus it traces out a polynomial curve.^ This property is called the
diagonal property.
Our special blossom b[^i, ^2] has two arguments. Blossoms with an arbitrary
number n of arguments are easily defined by the preceding three properties. A
blossom is an w-variate function b[^i, • • •, ^w] from W into E^ or E^. It is defined
by three properties:
Symmetry:
b[^i,...,y = b[7rai,...,y]

(3.19)

where 7t{ti^..., ^«) denotes a permutation of the arguments ^ 1 , . . . , t^. Thus, for
example b[^i, ti, t^] = h[t2, t^, t^],

Multiafftnity:
h[{ar + Ps), *] = ah[r, *] + ph[s, *]; a + )6 = 1.

(3.20)

Here, the symbol * indicates that there are the same arguments on both sides of
the equation, but their exact meaning is not of interest. Because of symmetry,
this property holds for all arguments, not just the first one.
Diagonality:
If all arguments of the blossom are the same: ^ = ^ j , . . . , ^„, then we obtain a
polynomial curve (to be discussed later). We will use the notation
h[t,...,t]=h[t<">]
if the argument t is repeated n times.
We defined vector arguments for linear interpolation in Section 3.1. Blossoms
may also have vector arguments, resulting in expressions such as b[/7, r, s]. If we
assume (without loss of generality) that the first argument of a blossom is a vector
h = b — a, then the multiaffine property becomes
h[b - ^, *] = b[fo, *] - h[a, *].

(3.21)

Thus if (at least) one of the blossom arguments is a vector, then the blossom
value is a vector. For example, if we denote by 1 the ID unit vector, then
b[l, r, s] = b[l, r, s] - b[0, r, s] or b[l, r, s] = b[3, r, s] — b[2, r, s\
4 This kind of curve will later be called a Bezier curve.

3.4 Blossoms

33

As an application of the blossom properties, let us derive a formula that will
be used later. We consider the special case when a blossom argument is of the
form {ar + ^s)^^^. For this, we get
b[(ar + )Ss)<"^] = Y. (''\'P''~Hr^'^.
.=0

s^^"^^].

(3.22)

-'

We refer to this equation as the Leibniz formula,^
The proof is by induction. The case « = 1 is a trivial start. The inductive step
proceeds as follows (keeping in mind that (^^-^) = (^J = 0):

r=0

^'^

Now we transform the index of the first sum and let the second sum run to w + 1:

n+l

. .

i=0

The first sum may start with / = 0. Keeping in mind the recursion

n-\-t\ ( ^ \
i

f^\

)~\i-l)^\i)

we can combine the last two sums and get
h[(ar + ^s)<^+l^] - Y ] f"" ^ -^V^'^^+^-^bfr^^'^, s^^+l-^"^],
which concludes our proof. This result will be used several times later on.
5 It has the structure of Leibniz's rule for higher-order derivatives of a product of functions.

54

Chapter 3 Linear Interpolation
A different form of (3.22) is sometimes useful:
b[(ar + )Ss)<">]= Y.

(.''.)«'>b[r<^=^,s<^=^]

(3.23)

where
/ w \ _ w!
V,/7

^ '

5.5 Barycentric Coordinates in tiie Plane
Barycentric coordinates were discussed in Section 3.1, where they were used in
connection with straight lines. Now we will use them as coordinate systems when
dealing with the plane. Planar barycentric coordinates are at the origin of affine
geometry—they were introduced by F. Moebius in 1827; see his collected works
[429].
Consider a triangle with vertices a, b, c and a fourth point p, all in E^. It IS
always possible to write p as a barycentric combination of a, b, c:
(3.24)

p = u2i + vh + wc.

A reminder: if (3.24) is to be a barycentric combination (and hence geometrically
meaningful), we require that
(3.25)

u-\-v-]-w=l,

The coefficients u := (u, v^ w) are called barycentric coordinates of p with respect
to a, b, c. We will often drop the distinction between the barycentric coordinates
of a point and the point itself; we then speak of "the point u."
If the four points a, b, c, and p are given, we can always determine p's
barycentric coordinates u,v^w: Equations (3.24) and (3.25) can be viewed as
a linear system of three equations^ in three unknowns u, v, w. The solution is
obtained by an application of Cramer's rule:

^=

6

area(p, b, c)
7—r—:'
area(a, b, c)

^^

area(a, p, c)
TTT'
area(a, b, c)

Recall that (3.24) is shorthand for two scalar equations.

^"^

area(a, b, p)
TTT^'
area(a, b, c)

^^'^^^

3.5 Barycentric Coordinates in the Plane

35

Actually, Cramer's rule makes use of determinants; they are related to areas by
the identity
b.
area(a, b, c) = -

(3.27)
1

1

1

We note that in order for (3.26) to be well defined, we require area(a, b, c) ^ 0,
which means that a, b, c must not lie on a straight line.
Because of their connection with barycentric combinations, barycentric coordinates are affinely invariant: let p have barycentric coordinates w, z/, w with
respect to a, b, c. Now map all four points to another set of four points by an
affine map O. Then Op has the same barycentric coordinates w, v^ w with respect
to Oa, Ob, Oc.
Figure 3.5 illustrates more of the geometric properties of barycentric coordinates.

Figure 3.5 Barycentric coordinates: let p = wa -h i/b + wc. The two figures show some of the ratios
generated by certain straight lines through p.

36

Chapter 3 Linear Interpolation

f

^

Figure 3.6 Barycentric coordinates: a triangle defines a coordinate system in the plane. Points with
three positive barycentric coordinates: white. With one negative barycentric coordinate:
light gray. With two negative barycentric coordinates: dark gray.
An immediate consequence of Figure 3.5 is known as Ceva's theorem:
ratio(a, p^, b) • ratio(b, p^, c) • ratio(c, p^, a) = 1.
More details on this and related theorems can be found in most geometry books
(e.g., Cans [253] or Berger [52], or Boehm and Prautzsch [85]).
Any three noncoUinear points a, b, c define a barycentric coordinate system in
the plane. The points inside the triangle a, b, c have positive barycentric coordinates, w^hereas the remaining ones have (some) negative barycentric coordinates.
Figure 3.6 shoves more.
We may use barycentric coordinates to define bivariate linear interpolation.
Suppose v^e are given three points pj, p2, P3 G E^. Then any point of the form
p = p(u) = p(w, V, w) = upi -f z;p2 + wp2>

(3.28)

with u-{-v-\-w=l
lies in the plane spanned by pi, p2, P3. This map from E^ to
E^ is called linear interpolation. Since u-\-v-\-w=l,wt
may interpret w, v, w as
barycentric coordinates of p relative to pi, p2, P3. We may also interpret w, v, w as
barycentric coordinates of a point in E^ relative to some triangle a, b, c G E-^. Then
(3.28) may be interpreted as a map of the triangle a, b, c G E^ onto the triangle
Pl5 P25 P3 ^ ^^- We call the triangle a, b, c the domain triangle. Note that the actual
location or shape of the domain triangle is totally irrelevant to the definition
of linear interpolation. (Of course, we must demand that it be nondegenerate.)
Since we can interpret w, v, w as barycentric coordinates in both two and three
dimensions, it follows that linear interpolation (3.28) is an affine map.

3.6 Tessellations

37

Barycentric coordinates are not restricted to one and two dimensions; they
are defined for spaces of higher dimensions as well. For example, in 3D, any
nondegenerate tetrahedron with vertices p^, p2, P3, P4 may be used to write any
point p as p = u^pi -\- U2P2 + ^3P3 + ^4p4-

5.6 Tessellations
When dealing with sequences of straight line segments, we were in the context of
piecewise linear interpolation. We may also consider more than one triangle, thus
introducing bivariate piecewise linear interpolation. Although straight line segments are combined into polygons in a straightforward way, the corresponding
concepts for triangles are not so obvious; they are the subject of this section.
We will first introduce the concept of a Dirichlet tessellation; this will lead to
an efficient way to deal with triangles. So consider a collection of points p^ in
the plane. We are going to construct influence regions around each point in the
following way: suppose each point is a transmitter for a cellular phone network.
As a car moves through the points p^, its phone should always be using the closest
transmitter. We may think of each transmitter as having an area of influence
around it: whenever a car is in a given transmitter's area, its phone switches to
that transmitter. More technically speaking, we associate with each point p^ a
tile T^ consisting of all points p that are closer to p^ than to any other point
p^. The collection of all these tiles is called the Dirichlet tessellation of the given
point set.^ Two points are called neighbors if their tiles share a common edge.
See Figure 3.7.
It is intuitively clear that the tile edges should consist of segments taken from
perpendicular bisectors of neighboring points. This observation directly leads to
a recursive construction that is due to R. Sibson [576]: suppose that we already
constructed the Dirichlet tessellation for a set of points, and we now want to
add one more point p^. First, we determine which of the previously constructed
tiles is occupied by pi; referring to Figure 3.8, let us assume it is T^. We now
draw all perpendicular bisectors between p^ and its neighbors, thus forming T^.
Continuing in this manner, we can construct the tessellation for an arbitrary
number of points. Each point is thus in the "center" of a tile, most of them finite,
but some infinite. It is not hard to see that all points with infinite tiles determine
the convex hull of the data points; see Section 2.1 for a definition.
Although the preceding method may not be the most efficient one to construct
the Dirichlet tessellation for a set of points, it is very intuitive, and also forms
the basis of the following fundamental theorem. The tile T^^ is formed by cutting
7 This structure is also known as a Voronoi diagram or Thiessen regions.

38

Chapter 3 Linear Interpolation

Figure 3.7 Dirichlet tessellations: a point set and its tile edges.

Figure 3.8 Dirichlet tessellations: a new point is inserted into an existing tessellation; its tile is
outlined.

3.7 Triangulations

39

out parts of p^^'s neighboring tiles. Let Aj be the area cut of T^, and let A be the
area of T^. Then we can write pi as a barycentric combination of its neighbors
(note that ^ ^ • = ^ ) :

p^^E-jp-

(3.29)

This identity is also due to R. Sibson [576]; in case the summation is over only
three neighbors, it reduces to the barycentric coordinates of Section 3.5.

5.7 Triangulations
The Dirichlet tessellation of a set of points determines another fundamental
structure that is connected with the point set: its Delaunay triangulation. If we
connect all neighboring points, we have created a set of triangles that cover the
convex hull of the point set and that have the given points as their vertices.
Figure 3.9 was created in this way from the configuration of Figure 3.7. The
points with infinite tiles are now connected; they are called boundary points of
the triangulation.
We should mention one problem: although the Dirichlet tessellation is unique,
the Delaunay triangulation may not be. As an example, consider four points
forming a square: either diagonal produces a valid Delaunay triangulation. Four
points that have no unique Delaunay triangulation are called neutral sets; such
points are always cocircular.

Figure 3.9 Delaunay triangulations: a point set and its Delaunay triangulation.

40

Chapter 3 Linear Interpolation
Clearly, there are many valid triangulations of a given point set. For example,
every convex set of four points allov^s tv^o different triangulations. It is now^ time
to introduce the concept of a triangulation of a point set that is more general
than the Delaunay triangulation. A triangulation T of a set of 2D points {pj is
a collection of triangles such that
• The vertices of the triangles consist of the p^
• The interiors of any tv^o triangles do not intersect
• If two triangles are not disjoint, then they share either a vertex or an edge
An important implementation aspect is the type of data structure to be used
for triangulations. Data sets v^ith several million points are not unheard of, and
for those, an intelligent structure is crucial. Such a structure should have the
foUov^ing elements:
1. A point collection of x, y-coordinate pairs
2. A collection of triangles, each pointing to three elements in the point list
and also to three elements in the triangle collection, namely, those that
designate a triangle's three neighbors^
These collections are best realized in the form of linked lists, for ease of inserting
and deleting points. This data structure goes back to F. Little, v^ho implemented
it in 1978 at the University of Utah.
As it turns out, the Delaunay triangulation is one of the "nicer" triangulations.
Intuitively, we might say that a triangulation is "nice" if it consists of triangles
that are close to being equilateral. If w^e compare two different triangulations of
a point set, w^e might then compute the minimal angle of each triangle. The
triangulation that has the largest minimal angle w^ould be labeled the better
one. Of all possible triangulations, the Delaunay triangulation is the one that
is guaranteed to produce the largest minimal angle; for a proof, see Lawson
[375]. The Delaunay triangulation is thus said to satisfy the maxmin criterion.
One might also consider the triangulation that satisfies the minmax criterion:
the triangulation w^hose maximal angle is minimal. These triangulations are not
easy to compute; one reason is that their neutral point sets are fairly complex,
see Hansford [312].
A major use of triangulations is in piecewise linear interpolation: suppose that
at each data point p^ v^e are given a function value Zk- Then v^e may construct a
linear interpolant—using linear interpolation from Section 3.5—over each of the
8

Boundary triangles may have only one or two neighbors.

3.8 Problems

41

triangles. We obtain a faceted, continuous surface that interpolates to all given
data. This surface is not smooth, but it will give a decent idea of the shape of
the given data. One application is in cartography: here, the given data points
might be coordinates obtained from satellite readings and the function values
might be their elevations. Our piecewise linear surface is an approximation to
the landscape being surveyed.
Once function values are involved, it may be advantageous to construct a
triangulation that reflects this information. Such triangulations are called data
dependent^ see Dyn, Levin, and Rippa [180] or Brow^n [92]. Here, one does not
just consider triangles in the plane, but rather the 3D triangles generated by the
data points (x^, y^^, Zk).

5.8 Problems
1 In the definition of the variation diminishing property, we counted the
crossings of a polygon with a plane. Discuss the case when the plane
contains a whole polygon leg.
* 2 We defined the convex hull of a point set to be the set of all convex
combinations formed by the elements of that set. Another definition is the
following: the convex hull of a point set is the intersection of all convex
sets that contain the given set. Show that both definitions are equivalent.
* 3 Our definition of barycentric combinations gives the impression that it
needs the involved points expressed in terms of some coordinate system.
Show that this is not necessary: draw five points on a piece of paper, assign
a weight to each one, and construct the barycenter of your points using a
ruler (or compass and straightedge if you are more classically inclined).
Remark: For this construction, it is not necessary for the weights to
sum to one. This is so because the geometric construction remains the
same if we multiplied all weights by a common factor. In fact, one may
replace the concept of points (having mass one and requiring barycentric
combinations as the basic point operation) by that of mass points^ having
arbitrary weights and yielding their barycenter (with the combined mass of
all points) as the basic operation. In such a setting, vectors would also be
mass points, but with mass zero.^
* 4 Let a triangulation consist of b boundary points and of / interior points.
Show that the number of triangles is 2i -\- b — 2.
9

I was introduced to this concept by A. Swimmer. It was developed by H. Grassmann in
1844.

42

Chapter 3 Linear Interpolation
PI Let three points be given by
bn
For s = 0 , 0 . 0 5 , 0 . 1 , . . . , 1 and ^ = 0 , 0 . 0 5 , 0 . 1 , . . . , 1, plot the points b[s, t]
as defined by (3.16). Mark each point by a circle with radius 0.4.
P2 There is a 2D triangulation data set on this book's web site. Plot that
triangulation using gray shades or colors such that no two neighboring
triangles have the same color.
P3 Use the recursive algorithm from Section 3.6 to implement Dirichlet tessellations.

The de Casteljau
Algorithm

I he algorithm described in this chapter is probably the most fundamental
one in the field of curve and surface design, yet it is surprisingly simple. Its
main attraction is the beautiful interplay between geometry and algebra: a very
intuitive geometric construction leads to a pow^erful theory.
Historically, it is v^ith this algorithm that the w^ork of de Casteljau started
in 1959. The only w^ritten evidence is in [145] and [146], both technical reports
that are not easily accessible. De Casteljau's v^ork w^ent unnoticed until W. Boehm
obtained copies of the reports in 1975. Since then, de Casteljau's w^ork has gained
more popularity.

4.1

Parabolas
We give a simple construction for the generation of a parabola; the straightforv^ard generalization w^ill then lead to Bezier curves. Let bo, b^, hi be any three
points in E^, and let t eR. Construct
bj(0 = (1 - t)ho + ^bi,
h\it) = (1 - Obi + th2,

hl(t) = (l-t)hl(t)

+ th\(t).

Inserting the first tv^o equations into the third one, v^e obtain
hlit) = (1 - t)% + 2t(l - Obi + t%.

(4.1)
43

44

Chapter 4 The de Casteljau Algorithm

0
Figure 4.1

t

Parabolas: construction by repeated linear interpolation.

This is a quadratic expression in t (the superscript denotes the degree), and so
hi^it) traces out a parabola as t varies from —oo to +oo. We denote this parabola
by b^. This construction consists of repeated linear interpolation^ its geometry is
illustrated in Figure 4.1. For t between 0 and 1, b^(^) is inside the triangle formed
by bo, bi, b2; in particular, b^(0) = bo and b^(l) = b2.
Inspecting the ratios of points in Figure 4.1, we see that

ratio(bo, b j , bi) = ratio(bi, bj, hi) = ratio(bJ, b^, b^) = t/(l — t).

Thus our construction of a parabola is affinely invariant because piecewise linear
interpolation is affinely invariant; see Section 3.2.
We also note that a parabola is a plane curve, since h^{t) is always a barycentric
combination of three points, as is clear from inspecting (4.1). A parabola is a
special case of conic sections^ which will be discussed in Chapter 12.
Finally we state a theorem from analytic geometry, closely related to our
parabola construction. Let a, b, c be three distinct points on a parabola. Let the
tangent at b intersect the tangents at a and c in e and f, respectively. Let the
tangents at a and c intersect in d. Then ratio(a, e, d) = ratio(e, b, f) = ratio(d, f, c).
This three tangent theorem describes a property of parabolas; the de Casteljau
algorithm can be viewed as the constructive counterpart. Figure 4.1, although
using a different notation, may serve as an illustration of the theorem.

4.2 The de Casteljau Algorithm

45

4.2 The de Casteljau Algorithm
Parabolas are plane curves. However, many applications require true space
curves.^ For those purposes, the previous construction for a parabola can be
generalized to generate a polynomial curve of arbitrary degree n\

de Casteljau algorithm:
Given: bo, b j , . . . , b„ G E^ and t e M,
set

mt) ^ (1 - t)hi\t)+th^-iit)

{'lo;•••;;_,

(4.2)

and h^(t) = b/. Then h^it) is the point with parameter value t on the Bezier
curve b", hence \y^{t) = ^(t).
The polygon P formed by b o , . . . , b„ is called the Bezier polygon or control
polygon of the curve h^? Similarly, the polygon vertices b^ are called control
points or Bezier points. Figure 4.2 illustrates the cubic case.
Sometimes we also write b^(^) = B [ b o , . . . , b„; ^] = B[P; t] or, shorter, b^ =
B [ b o , . . . , b„] = BP. This notation-^ defines B to be the (linear) operator that
associates the Bezier curve with its control polygon. We say that the curve
B [ b o , . . . , b„] is the Bernstein-Bezier approximation to the control polygon, a
terminology borrowed from approximation theory; see also Section 6.9.
The intermediate coefficients b[(^) are conveniently written into a triangular
array of points, the de Casteljau scheme. We give the example of the cubic case:

bl K
hi b\ hi
b3 bl b2 hi.

(4.3)

This triangular array of points seems to suggest the use of a two-dimensional
array in writing code for the de Casteljau algorithm. That would be a waste of
1 Compare the comments by P. Bezier in Chapter 1!
2 In the cubic case, there are four control points; they form a tetrahedron in the 3D case. This
tetrahedron was already mentioned by W. Blaschke [65] in 1923; he called it "osculating
tetrahedron."
3 This notation should not be confused with the blossoming notation used later.

46

Chapter 4 The de Casteljau Algorithm

Figure 4.2 The de Casteljau algorithm: the point h^it) is obtained from repeated linear interpolation.
The cubic case w = 3 is shown for t = 1/3.

Example 4.1

Computing a point on a Bezier curve with the de Casteljau algorithm.
A de Casteljau scheme for a planar cubic and for t = j :

ro]
0

[2

8

[2
[4'
0

"0"
1
4
2
''6''
1

2

3
- 2-

[5]
3

r 7 1
2
3
- 7J

Storage, however: it is sufficient to use the left column only and to overwrite it
appropriately.
For a numerical example, see Example 4.1. Figure 4.3 shows 60 evaluations
of a Bezier curve. The intermediate points W- are also plotted and connected."^

4

Although the control polygon of thefigureis symmetric, the plot is not. This is due to the
organization of the plotting algorithm.

4.3 Some Properties of Bezier Curves

47

Figure 4.3 The de Casteljau algorithm: 60 points are computed on a degree six curve; all intermedediate points b^ are shown.

4.5 Some Properties of Bezier Curves
The de Casteljau algorithm allows us to infer several important properties of
Bezier curves. We will infer these properties from the geometry underlying the
algorithm. In the next chapter, we will show how they can also be derived
analytically.
Afflne invariance. Affine maps were discussed in Section 2.2. They are in the
tool kit of every CAD system: objects must be repositioned, scaled, and so
on. An important property of Bezier curves is that they are invariant under
affine maps, which means that the following two procedures yield the same
result: (1) first, compute the point h^(t) and then apply an affine map to it; (2)
first, apply an affine map to the control polygon and then evaluate the mapped
polygon at parameter value t.
Affine invariance is, of course, a direct consequence of the de Casteljau
algorithm: the algorithm is composed of a sequence of linear interpolations
(or, equivalently, of a sequence of affine maps). These are themselves affinely
invariant, and so is a finite sequence of them.
Let us discuss a practical aspect of affine invariance. Suppose we plot a cubic
curve b^ by evaluating at 100 points and then plotting the resulting point array.
Suppose now that we would like to plot the curve after a rotation has been
applied to it. We can take the 100 computed points, apply the rotation to each
of them, and plot. Or, we can apply the rotation to the 4 control points, then
evaluate 100 times and plot. The first method needs 100 applications of the
rotation, whereas the second needs only 4!

48

Chapter 4 The de Casteljau Algorithm
Affine invariance may not seem to be a very exceptional property for a
useful curve scheme; in fact, it is not straightforward to think of a curve
scheme that does not have it (exercise!). It is perhaps v^orth noting that Bezier
curves do not enjoy another, also very important, property: they are not
projectively invariant. Projective maps are used in computer graphics when
an object is to be rendered realistically. So if we try to make life easy and
simplify a perspective map of a Bezier curve by mapping the control polygon
and then computing the curve, we have actually cheated: that curve is not the
perspective image of the original curve! More details on perspective maps can
be found in Chapter 12.
Invariance under affine parameter transformations. Very often, one thinks
of a Bezier curve as being defined over the interval [0,1]. This is done because it
is convenient, not because it is necessary: the de Casteljau algorithm is "blind"
to the actual interval that the curve is defined over because it uses ratios only.
One may therefore think of the curve as being defined over any arbitrary
interval a, t^'^, r ' > ] = (1 - ^)b[0<"-^-'+i>, t^'-^^, r ' > ]
+ fb[o<«-'-'>, t-"-^^, i<'+i=^].

(4.10)

The point on the curve is given by b[f^^^].
We may also consider the blossom of a Bezier curve that is not defined over
[0,1] but over the more general interval [a, b]. Proceeding exactly as above—
but nov^ using (4.4)—v^e find that the Bezier points b/ are found as the blossom
values

b,=b[^<^-^>,6^^>].

(4.11)

Thus a cubic over u e [a, b] has Bezier points h[a^ a, a\ h[a^ a,fo],h[a,fc,fc],b[fc,fo,b].
If the original Bezier curve was defined over [0,1], the Bezier points of the one
corresponding to [a, b] are simply found by four calls to a blossom routine! See
also Figure 4.5.

b[l,l,0]
0

b[l,l,l]

^w-#

h[b,bM

h[a,b,b]

h[a,a,b]
h[a,a,a

b[0,0,0]

0

a

b[0,0,l]

h

i

Figure 4.5 Subdivision: the relevant blossom values.

4.5 Implementation

53

We may also find explicit formulas for blossoms; here is the case of a cubic:

= (1 - t{)[(\ - ^2)b[0,0, ^3] + ^2b[0,1, ^3]] + h[(X - ^2)b[0,1, ^3]

= b[0,0,0](l-^i)(l-?2)a-^3)
+ b[0,0, !][(! - ^i)(l - ^2)^3 + (1 - ^1)^2(1 - ^3) + ^i(l - t2){l - ^3)]
+ b[0,1, l][^i^2(l - ^3) + ^i(l - ^2)^3 + (1 - ^1)^2^3]
+ b[l,l,l]^1^2^3.
For each step, we have exploited the fact that blossoms are multiaffine, following
the inductive proof of the Leibniz equation (3.22).
We should add that every multivariate polynomial function may be interpreted
as the blossom of a Bezier curve—as long as it is both symmetric and multiaffine.

4.5 Implementation
The header of the de Casteljau algorithm program is:
float decas(degree,coeff,t)
/* uses de Casteljau to compute one coordinate
value of a Bezier curve. Has to be called
for each coordinate (x,y, and/or z) of a control polygon.
Input:
degree: degree of curve.
coeff:
array with coefficients of curve,
t:
parameter value.
Output: coordinate value.

V
This procedure invites several comments. First, we see that it requires the use of
an auxiliary array coeff a. Moreover, this auxiliary array has to be filled for each
function call! So on top of the already high computational cost of the de Casteljau
algorithm, we add another burden to the routine, keeping it from being very
efficient. A faster evaluation method is given at the end of the next chapter.

54

Chapter 4 The de Casteljau Algorithm
To plot a Bezier curve, we would then call the routine several times:
void bez_to_points(degree,npoints,coeff,points)
/*
Converts Bezier curve into point sequence. Works on
one coordinate only.
Input:
degree: degree of curve.
npoints: # of coordinates to be generated, (counting
from 0!)
coeff:
coordinates of control polygon.
Output:

points:

coordinates of points on curve.

Remark: For a 2D curve, this routine needs to be called twice,
once for the x-coordinates and once for y.

*/
The last subroutine has to be called once for each coordinate, that is, two or
three times. The main program decasmai n. c on the enclosed disk gives an example
of how to use it and how to generate postscript output.

4.6 Problems
1 Suppose a planar Bezier curve has a control polygon that is symmetric
with respect to the y-axis. Is the curve also symmetric with respect to
the y-axis? Be sure to consider the control polygon (—1, 0), (0,1), (1,1),
(0, 2), (0,1), ( - 1 , 1), (0, 2), (0, 1), (1, 0). Generalize to other symmetry
properties.
2 Use the de Casteljau algorithm to design a curve of degree four that has its
middle control point on the curve. More specifically, try to achieve

Five collinear control points are a solution; try to be more ambitious!
'' 3 The de Casteljau algorithm may be formulated as
B[bo, ...,K;t]

= a - OB[bo,..., K-i; t] + ^ [ b i , . . . , b^; t].

Show that the computation count is exponential (in terms of the degree) if
you implement such a recursive algorithm in a language like C.

4.6 Problems

55

• 4 Show that every nonplanar cubic in E"^ can be obtained as an affine map of
the standard cubic (see Boehm [70])

xW =
PI Write an experimental program that replaces {1 — t) and t in the recursion
(4.2) by [1 — f{t)] and f{t)^ where f is some "interesting" function. Change
the routine decas accordingly and comment on your results.
P2 Rewrite the routine decas to handle blossoms. Evaluate and plot for some
"interesting" arguments.
P3 Experiment with the data set outl i ne_2D. dat on the floppy: try to recapture
its shape using one, two, and four Bezier curves. These curves should have
decreasing degrees as you use more of them.
P4 Then repeat the previous problem with outline_3D.dat. This data set is
three dimensional, and you will have to use (at least) two views as you
approximate the data points. The points, by the way, are taken from the
outline of a high heel shoe sole.

This Page Intentionally Left Blank

The Bernstein Form
of a Bezier Curve

D e z i e r curves can be defined by a recursive algorithm, which is how de Casteljau
first developed them. It is also necessary, however, to have an explicit representation for them; this will facilitate further theoretical development considerably.

5.1 Bernstein Polynomials
We will express Bezier curves in terms of Bernstein polynomials^ defined explicitly
by

B"^it)^{^^t\l-t)"-\

(5.1)

where the binomial coefficients are given by

V/7

1

0

else.

There is a fair amount of literature on these polynomials. We cite just a few:
Bernstein [53], Lorentz [399], Davis [133], and Korovkin [364]. An extensive
bibliography is given in Gonska and Meier [269].
Before we explore the importance of Bernstein polynomials to Bezier curves,
let us first examine them more closely. One of their important properties is that
they satisfy the following recursion:
B'lit) = (1 - t)B'l-\t)

+ tB'lzlit)

(5.2)
57

58

Chapter 5 The Bernstein Form of a Bezier Curve
with
^Om ^_ 1
Blit)

(5.3)

and
BJ(O = 0 for

;^{0,...,n}.

(5.4)

The proof is simple:

=C7>'<'-'»-'-(::0''<'-'>

n—t

Another important property is that Bernstein polynomials form a partition of
unity:
1 ] B « ( 0 = 1.

(5.5)

This fact is proved with the help of the binomial theorem:

1=[^+(1 - t)r=J2 C'^^i^ - tr-^=E ^^^)Figure 5.1 shows the family of the four cubic Bernstein polynomials. Note that
the B" are nonnegative over the interval [0,1].
We are now ready to see why Bernstein polynomials are important for the
development of Bezier curves. Recall that a Bezier curve may be written as b[^'^"^]
in blossom form. Since t = (I — t) - 0 -{-1 • 1^ the blossom may be expressed as
b[(l -t)'0
+ t' l)^"""^], and now the Leibniz formula (3.22) directly yields
n

hit) = h[t<">] = J2^iBlit)
«=o
since b, = h[0<"-'>, !<'>] according to (4.9).

(5.6)

5.1 Bernstein Polynomials
Ik 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

y1

11 i
T-l^

U ?^n

M//

^"^ /

\

Figure 5.1

111

T^^

59

Nl

P^

M

P^

/

^ M

Bernstein polynomials: the cubic case.

Similarly, the intermediate de Casteljau points W- can be expressed in terms of
Bernstein polynomials of degree r:

(5.7)

This follows directly from

and the Leibniz formula.
Equation (5.7) shows exactly how the intermediate point W- depends on the
given Bezier points hj. Figure 5.2 shows how these intermediate points form
Bezier curves themselves.
With the intermediate points W- at hand, we can write a Bezier curve in the
form

b^(0 = ^b;(OBpw.

(5.8)

/=0

This is to be interpreted as follows: first, compute r levels of the de Casteljau
algorithm with respect to t. Then, interpret the resulting points b[(t) as control
points of a Bezier curve of degree n — r and evaluate it at t.

60

Chapter 5 The Bernstein Form of a Bezier Curve

Figure 5.2 The de Casteljau algorithm: 50 points are computed on a quartic curve, and the intermediate points b^ are connected.

5.2 Properties of Bezier Curves
Many of the properties in this section have aheady appeared in Chapter 4. They
were derived using geometric arguments. We shall now^ rederive several of them,
using algebraic arguments. If the same heading is used here as in Chapter 3, the
reader should look there for a complete description of the property in question.
Affine invariance. Barycentric combinations are invariant under affine maps.
Therefore, (5.5) gives the algebraic verification of this property. We note again
that this does not imply invariance under perspective maps!
Invariance under affine parameter transformations. Algebraically, this property reads

f:hfiiit)=j2hfi"A'^).
i=0

{5.9)

/=0

Convex hull property. This follows, since for t e [0,1], the Bernstein polynomials are nonnegative. They sum to one as show^n in (5.5). For values of t outside
[0,1], the convex hull property does not hold; Figure 5.3 illustrates.

5.2 Properties of Bezier Curves

61

Figure 5.3 Convex hull property: a quartic Bezier curve is plotted for parameter values t e [—1,2].
Endpoint interpolation. This is a consequence of the identities

and {S.5), Here, 5^y is the Kronecker delta function: it equals one when its
arguments agree, and zero otherwise.
Symmetry. Looking at the examples in Figure 4.4, it is clear that it does not
matter if the Bezier points are labeled bo, b ^ , . . . , b„ or b„, b„_i, . . . , bg. The
curves that correspond to the two different orderings look the same; they differ
only in the direction in which they are traversed. Written as a formula:
n

n

^b,B;a) = ^b,_^B;(i-o.

(5.11)

This follows from the identity
(5.12)

B1(t) = Bl_.(l-t\

which follows from inspection of (5.1). We say that Bernstein polynomials are
symmetric with respect to t and 1 — t,
Invariance under barycentric combinations. The process of forming the Bezier
curve from the Bezier polygon leaves barycentric combinations invariant. For
Of + jS = 1, we obtain
n

n

n

j=0

j=0

j=0

In words: we can construct the weighted average of two Bezier curves either
by taking the weighted average of corresponding points on the curves, or
by taking the weighted average of corresponding control vertices and then
computing the curve.
This linearity property is essential for many theoretical purposes, the most
important one being the definition of tensor product surfaces in Chapter 14.
It is illustrated in Figure 5.4.

62

Chapters The Bernstein Form of a Bezier Curve

Figure 5.4 Barycentric combinations: the middle curve (black) is the average of the two outer curves
(gray).
Linear precision. The following is a useful identity:

J2-B';(t) = t,

(5.14)

which has the following application: suppose the polygon vertices by are
uniformly distributed on a straight line joining two points p and q:

^i=\^--)p+-^'^
\
nj

n

7 = 0,...,^.

The curve that is generated by this polygon is the straight line between p and
q, that is, the initial straight line is reproduced. This property is called linear
precision}
Pseudolocal control. The Bernstein polynomial B^ has only one maximum and
attains it at ^ = i/n. This has a design application: if we move only one of
the control polygon vertices, say, b^, then the curve is mostly affected by this
change in the region of the curve around the parameter value i/n. This makes
the effect of the change reasonably predictable, although the change does
affect the whole curve. As a rule of thumb (mentioned to me by P. Bezier),
the maximum of each JB^ is roughly 1/3; thus a change of b/ by three units will
change the curve by one unit.

5.5 The Derivatives of a Bezier Curve
We start with an identity, closely resembling Leibniz's formula for derivatives.
Let ^ be a point on the real line, and let ? be a vector in the associated ID linear
1 If the points are not uniformly spaced, v^e will also recapture the straight line segment.
However, it will not be linearly parametrized.

5.3 The Derivatives of a Bezier Curve

63

space. Then

h[(t + ?)<"^] = J2 ('')b[^<^-^^, v^^n

(5.15)

i=0
This is an immediate consequence of the Leibniz formula (3.22).
The derivative of a curve x(^) is typically defined as

dt

h-^0 h

We will be a little more precise and observe that f is a ID point, whereas /? is a
ID vector. We thus denote it by h and obtain
^

dt

= lim4[xa + ^)-x(0].
h^o \h\

Invoking (5.15), we have
dx(t)
,.
1
—-— = lim -^
dt

h-^o \h\

J2 ('')b[^"''"'^, i^^'n - Ht^'^n
.1=0

(5.16)

^^^

For / = 0, two terms b[^^"^] cancel. We expand the rest and factor in the term
\h\:

dt

+ ( : > t<"-^>,i,h + .
1^1 J
1^1 J V2/ I

_^

— lim I nb f
h^o \

L

We observe that -§- == 1. Taking the limit annihilates all other terms containing
\h\

h, and we thus have
dx(t)

"dT

=

nh[t^''-'^,l].

(5.17)

Figure 5.5 illustrates the cubic case.
From now on, we use the expression x(^) for the first derivative.
This has two possible interpretations. For the first one, we perform a de
Casteljau step with respect to 1, and then n — 1 steps with respect to t; as an
equation:
n—l

Kit) = n l](b;+i - hi)Bl-\t).
;=0

64

Chapter 5 The Bernstein Form of a Bezier Curve

Figure 5.5 Blossoms and derivatives: the underlying geometry.
This can be simplified somewhat by the introduction of the forward difference
operator A:
Aby = by+i-by.

(5.18)

We now have for the derivative of a Bezier curve:
n-l

x(t) = nY^ AbyBp^O;

Aby € R \

(5.19)

The derivative of a Bezier curve is thus another Bezier curve, obtained by differencing the original control polygon. However, this derivative Bezier curve does
not "live" in E^ any more! Its coefficients are differences of points, that is, vectors,
which are elements of R^. To visualize the derivative curve and polygon in E^, we
can construct a polygon in E^ that consists of the points a + Abo, • • • 5 a + Ab„_i.
Here a is arbitrary; one reasonable choice is a = 0. Figure 5,6 illustrates a Bezier
curve and its derivative curve (with the choice a = 0). This derivative curve is
sometimes called a hodograph. For more information on hodographs, see Forrest
[244], Bezier [59], or Sederberg and Wang [559].
For a second interpretation of (5.17), we first perform n — 1 steps of the
de Casteljau algorithm, resulting in the two points h^'^it) and h^~^(t). Now
performing one step with respect to 1 yields (after multiplication by n):

x(t) = n{h"f\t)-h"f\t)).

(5.20)

Thus the first derivative vector is a "byproduct" of the de Casteljau algorithm;
see Figure 4.2. The de Casteljau algorithm is not the fastest way to evaluate a
Bezier curve, but this property makes it a desirable tool: very often, we not only
need a point on a curve, but the derivative vector as well. Using (5.20), we get
both in parallel. The two ways of computing the derivative are shown in Example
5.1.

5.3 The Derivatives of a Bezier Curve

65

Figure 5.6 Derivatives: a Bezier curve and its first derivative curve (scaled down by a factor of three).
Note that this derivative curve does not change if a translation is applied to the original

Example 5.1 Two ways to compute derivatives.
To compute the derivative of the Bezier curve from Example 4.1, we could
form the first differences of the control points and evaluate the corresponding
quadratic curve at ^ = 1/2:

ro]
[ij

[8]

[oj
-4
-2

r4]

[ij
•

2

'3'
0

•

-1

Alternatively, v^e could compute the difference b^ — bg:
'5"

"2"

3
3 -? -? .

=

"3"
0

In both cases, the resuh needs to be multiphed by a factor of 3.

Higher derivatives follow the same pattern:

dr

hit""-'^, !<'•>].

(n — r)\

(5.21)

66

Chapter 5 The Bernstein Form of a Bezier Curve
To compute these derivatives from the Bezier points, wt first generalize the
forward difference operator (5.18): the iterated forward difference operator A^
is defined by
A^by = A^-^by+i - A^-^by.

(5.22)

We Ust a fev^ examples:
A \ = b,A b, = b,+i - b,
A^by = b,-+2-2b,-+i + b,^ \ = b,+3 - 3b,+2 + 3b^+i - b,-.
The factors on the right-hand sides are binomial coefficients, forming a Pascallike triangle. This pattern holds in general:

;=0 ^^^

The r^ derivative of a Bezier curve is now^ given by
^h"(t)
dt^

= - ^
£ A^b^Bj-W.
(n — r)\ ^
' '

(5.24)

Two important special cases of (5.24) are given by ^ = 0 and f = 1. Because of
(5.10), we obtain
^b«(0) = — ^ A %
at^
{n — r)\

(5.25)

ilb^(l) = _!!!_A%_^.
dt'
{n-r)\
"" '

(5.26)

and

Thus the r^^ derivative of a Bezier curve at an endpoint depends only on the r -\-l
Bezier points near (and including) that endpoint. For r = 0, we get the already
established property of endpoint interpolation. The case r = 1 states that bg and

5.3 The Derivatives of a Bezier Curve

67

Figure 5.7 Endpoint derivatives: the first and second derivative vectors at ^ = 0 are mukiples of the
first and second difference vectors at bo-

bi define the tangent at ^ = 0, provided they are distinct.^ Similarly, b„_i and b„
determine the tangent at ^ = 1. The cases r = 1, r = 2 are illustrated in Figure 5.7.
If we knov^ all derivatives of a function at one point, corresponding to ^ = 0,
say, we can generate its Taylor series. The Taylor series of a polynomial is just
that polynomial itself, in the monomial form:

xit) = J2 -.x^'how.
,"''•

Using (5.25), we have

(5.27)
7=0 ^ ' ^

The monomial form should be avoided wherever possible; it is very unstable
for floating-point operations.
If x(^) is defined over an interval [a^fc],(5.17) becomes
dx(t)
b- a

2

b[^<"-^>, 1].

(5.28)

In general, the tangent at bo is determined by bo and the first b^ that is distinct from boThus the tangent may be defined even if the tangent vector is the zero vector.

68

Chapter 5 The Bernstein Form of a Bezier Curve

5.4 Domain Changes and Subdivision
A Bezier curve b" is usually defined over the interval (the domain) [0,1], but it can
also be defined over any interval [0, c]. The part of the curve that corresponds to
[0, c] can also be defined by a Bezier polygon, as illustrated in Figure 5.8. Finding
this Bezier polygon is referred to as subdivision of the Bezier curve.
The unknown Bezier points Cj are found v^ithout much w^ork if we use the
blossoming principle from Section 4.4. There, (4.11) gave us the Bezier points
of a polynomial curve that is defined over an arbitrary interval [a, b]. We are
currently interested in the interval [0, c], and so our Bezier points are:
q = b[0<^-^'>,c<^^].
Thus each q is obtained by carrying out / de Casteljau steps with respect to c, in
nonblossom notation:
Cy-b^W.

(5.29)

This formula is called the subdivision formula for Bezier curves.
Thus it turns out that the de Casteljau algorithm not only computes the point
b"(c), but also provides the control vertices of the Bezier curve corresponding to
the interval [0, c]. Because of the symmetry property (5.11), it follows that the
control vertices of the part corresponding to [c, 1] are given by the \y~\ Thus,
in Figures 4.1 and 4.2, we see the two subpolygons defining the arcs from b"(0)
to b"(0 and from b"(0 to b"(l).
Instead of subdividing a Bezier curve, we may also extrapolate it: in that case,
we might be interested in the Bezier points dj corresponding to an interval [l,d].
They are given by

It should be mentioned that extrapolation is not a numerically stable process,
and should be avoided for large values of d.
Subdivision for Bezier curves, although mentioned by de Casteljau [146],
was rigorously proved by E. Staerk [578]. Our blossom development is due to
Ramshaw [498] and de Casteljau [147].
Subdivision may be repeated: we may subdivide a curve at ^ = 1/2, then split
the two resulting curves at ^ = 1/2 of their respective parameters, and so on. After
k levels of subdivisions, we end up with 2^ Bezier polygons, each describing a
small arc of the original curve. These polygons converge to the curve if we keep
increasing k, as was shown by Lane and Riesenfeld [369].
Convergence of this repeated subdivision process is very fast (see Cohen and
Schumaker [123] and Dahmen [131]), and thus it has many practical applica-

5.4 Domain Changes and Subdivision

69

bi

0

1

Figure 5.8 Subdivision: two Bezier polygons describing the same curve: one (the b^) is associated v^ith
the parameter interval [0,1], the other (the Cj) with [0, c].
tions. We shall discuss here the process of intersecting a straight line with a Bezier
curve. Suppose we are given a planar Bezier curve and we wish to find intersection
points with a given straight line L, if they exist.
If the curve and L are far apart, we would like to be able to flag such
configurations as quickly as possible, and then abandon any further attempts
to find intersection points. To do this, we create the minmax box of the control
polygon: this is the smallest rectangle, with sides parallel to the coordinate axes,
that contains the polygon. It is found very quickly, and by the convex hull
property of Bezier curves, we know that it also contains the curve. Figure 5.9
gives an example.
Having found the minmax box, it is trivial to determine if it interferes with
L; if not, we know we will not have any intersections. This quick test is called
trivial reject.
Now suppose the minmax box does interfere with L. Then there may be an
intersection. We now subdivide the curve at ^ = 1/2 and carry out our trivial
reject test for both subpolygons.^ If the outcome is still inconclusive, we repeat.
Eventually the size of the involved minmax boxes will be so small that we can
simply take their centers as the desired intersection points.
The routine intersect employs this idea, and a little more: as we keep subdividing the curve, zooming in toward the intersection points, the generated
subpolygons become simpler and simpler in shape. If the control points of a
3 The choice t = 1/2 is arbitrary, but works well. We might try to find better places to
subdivide, but it is cheaper to just perform a few more subdivisions instead.

70

Chapter 5 The Bernstein Form of a Bezier Curve

Figure 5.9 The minmax box of a Bezier curve: the smallest rectangle that contains the curve's control
polygon.

Figure 5.10

Subdivision: finding the intersections of a curve with a line (dashed). Note the clustering
of minmax boxes near the intersection points.
polygon are almost collinear, we may replace them with a straight line. We could
then intersect this straight line with L in order to find an intersection point. The
extra work here lies in determining if a control polygon is "linear" or not. In our
case, this is done by the routine checkflat. Figure 5.10 gives an example. Note
how the subdivision process finds all intersection points. These points will not,
however, be recorded by increasing values of t.

5.5 Composite Bezier Curves

Figure 5.11

71

Font design: the characters in this book are stored as a sequence of cubic Bezier curves.

5.5 Composite Bezier Curves
Curves may be composed of several Bezier curves in order to generate shapes
that are too complex for a single Bezier curve to handle. For example, Figure
5.11 shows how composite Bezier curves may be used in font design^
In piecing Bezier curves together, we need to control the smoothness of the
resulting curve. Let bo,. . ., b3 and b3,. . ., b^ be the Bezier points of two cubic
curve segments x_ and x^. Since they both share the point b3, they clearly form
a continuous, or C^, curve. With this minimal continuity requirement, the two
curves may form a corner; for several examples, see Figure 5.11.
But if we want to ensure that the two pieces meet smoothly, more care is called
for. Based on our knowledge of endpoint derivatives from Section 5.3, the three
points b2, b3, h^ must be coUinear. That condition ensures that the tangent^ at
b3 is the same for both curves. Again, consult Figure 5.11 for examples. Curves
with a continuously changing tangent are called G^ or first-order geometrically
continuous; see Chapter 11.
A stronger condition is to require that the two curve segments form a C^, or
continuously differentiable curve. Since the derivative of a curve (more precisely,
the length of the derivative vector) depends on the domain of the curve, we need to
introduce domains for our two curve segments. We adopt the convention that x_
is defined over an interval [a, b] and that x+ is defined over [b, c]. The derivatives

4 This book was printed using the PostScript language. It represents all characters as
piecewise cubic Bezier curves in order to have a scalable font set. As an estimate, the
text in this book is made up using about 10 million cubic Bezier curves.
5

By "tangent," we refer to the tangent line, not to the derivative vector!

72

Chapter 5 The Bernstein Form of a Bezier Curve

?o
c
Figure 5.12

b

a

Composite curves: a C^ example.
y

1^
3 1

\\\M^

1

1

1i

|TK||
ttlM

W^

Figure 5.13 Composite curves: a C^ example.
of both segments at parameter value b are nov^ obtained using (5.28):
^

b-a

[b3-b2]=^[b4-b3].

-y

(5.30)

A geometric interpretation is that the ratio of the three points b2, b3, b4 is the
same as the ratio of the three parameter values a^fc,c. This is a much stronger
condition than that for G^ continuity above!
Figures 5.12 and 5.13 illustrate this difference. The composite parametric
curves—in the x, y-coordinate systems—are identical. The difference is their
domains: in Figure 5.12, w^e chose a^b,c = 0,1,2. Thus ratio(^, b^c) = 1 v^hile
the figure suggests that ratio(b2, b3, b4) = 1/3. Hence the composite curve is not
C \ despite the collinearity of the points b2, b3,b4. This is demonstrated by the
cross plot (y-part only): each component must form a C^ function for a curve to
be C^. Clearly, the }'-component is not C \
If v^e adjust the domain, hovs^ever, such that the range geometry is reflected by
the domain geometry, wt can achieve C^. This is show^n in Figure 5.13, v^here
now^ ratio(<3t,fc,c) = 1/3. This results in C^ components, and hence also in a C^
composite curve.

5.6 Blossom and Polar

73

Higher-order smoothness of composite curves is best dealt with in the context
of B-spline curves and blossoms; see Section 8.7.

5.6 Blossom and Polar
After the first de Casteljau step with respect to a parameter value ti, the resulting
b j c ^ l ) , . . . , h^_^(ti) may be interpreted as a control polygon of a curve pi(t) of
degree n — lAn the blossoming terminology from Section 4.4, we can write:
Pt{t)=b[ti,t<"-'>].
Invoking our knowledge about derivatives, we have:
n-l
i=0
n-l

n-l

= ^ [(1 - t,)h,+^A>i - h](t)]B'i-\t)+j2 W(^)^r'(^)
i=0

i=0
n-l

n-l

= {t, -1) ^[b,+i - h,w-\t)+Y.
i=0

b-w^r'(^)i=0

Therefore,
p^(t)=h(t)+^-^-^^h(t),
n at

(5.31)

The polynomial p^ is called first polar of h(t) with respect to t^. Figure 5.14
illustrates the geometric significance of (5.31): the tangent at any point h(t)
intersects the polar at pi(^). Keep in mind that this is not restricted to planar
curves, but is equally valid for space curves!
For the special case of a (nonplanar) cubic, we may then conclude the following: the polar pi lies in the osculating plane (see Section 11.2) of the cubic at b(^i).
If we intersect all tangents to the cubic with this osculating plane, we will trace
out the polar. We can also conclude that for three different parameters t^^ ti^ ^3,
the blossom value b[^i, ^2? ^3] is the intersection of the corresponding osculating
planes.
Another special case is given by b[0, f^"~^^]: this is the polynomial defined
by b o , . . . , ^n-l' Similarly, b[l, t^^"^^] is defined by b ^ , . . . , b^. This observation
may be used for a proof of (4.9).

74

Chapter 5 The Bernstein Form of a Bezier Curve

Figure 5.14 Polars: the polar pi(^) with respect to ti = 0.4 is intersected by the tangents of the given
curve h(t).
Returning to the general case, we may repeat the process of forming polars,
thus obtaining a second polar pi^2(0 = b[^i, ^25 ^^""^^L ^^^ so on. We finally
arrive at the « polar, which we have already encountered as the blossom
b[^i,.. • 5 ^«] oih(t). The relationship between blossoms and polars was observed
by Ramshaw in [499]. The preceding geometric arguments are due to S. JoUes,
who developed a geometric theory of blossoming as early as 1886 in [346].^

5.7 The Matrix Form of a Bezier Curve
Some authors (Faux and Pratt [228], Mortenson [433], Chang [106]) prefer to
write Bezier curves and other polynomial curves in matrix form. A curve of the
form
n

x(o-;^c,Q(f)
/=o

can be interpreted as a dot product:

Coit)
xit) = [ Co . . .

c„]

L C„{t) J
One can take this a step further and write
6 W. Boehm first noted the relevance of JoUes's w^ork to the theory of blossoming.

5.8 Implementation

Coit)

mm

C„{t) J

f"nO

mo„
m„

75

n^On

(5.32)
t".

The matrix M = {m^y} describes the basis transformation between the basis polynomials Q ( 0 and the monomial basis f.
If the Q are Bernstein polynomials, Q = B ^ the matrix M has elements

m/;-(-iy-^

(5.33)

a simple consequence of (5.27).
We list the cubic case explicitly:

M =

1 -3
3 -1
0
3 -6
3
0
0
3 -3
1
0
0
0

The matrix form (5.32) does not describe an actual Bezier curve; it is rather
the monomial form, which is numerically unstable and should be avoided where
accuracy in computation is of any importance. See the discussion in Section 24.3
for more details.

5.8 Implementation
First, we provide a routine that evaluates a Bezier curve more efficiently than
decas from the last chapter. It will have the flavor of Horner's scheme for the
evaluation of a polynomial in monomial form. To give an example of Horner's
scheme, also called nested multiplication^ we list the cubic case:
Co + tCi + t^C2 + t^C^ = Co + ^[Ci + t{C2 + ^€3)].
A similar nested form can be devised for Bezier curves; again, the cubic case:

h\t)
where s = 1 — ^. Recalling the identity

Q)A,)».QA3,

76

Chapter 5 The Bernstein Form of a Bezier Curve
we arrive at the following program (for the general case):
float hornbez(clegree,coeff,t)
/*
uses a Horner-like scheme to compute one coordinate
value of a Bezier curve. Has to be called
for each coordinate (x,y, and/or z) of a control polygon.
Input:
degree: degree of curve.
coeff: array with coefficients of curve,
t:
parameter value.
Output:
coordinate value.

V

To use this routine for plotting a Bezier curve, we would replace the call to decas
in bez_to_points by an identical call to hornbez. Replacing decas with hornbez
results in a significant savings of time: we do not have to save the control polygon
in an auxiliary array; also, hornbez is of order «, whereas decas is of order n^.
This is not to say, however, that we have produced superefficient code for
plotting points on a Bezier curve. For instance, we have to call hornbez once for
each coordinate, and thus have to generate the binomial coefficients n_choose_i
twice. This could be improved by writing a routine that combines the two calls. A
further improvement could be to compute the sequence of binomial coefficients
only once, and not over and over for each new value of t. All these (and possibly
more) improvements would speed up the program, but would be less modular
and thus less understandable. For the code in this book, modularity is placed
above efficiency (in most cases).
We also include the programs to convert from the Bezier form to the monomial
form:
voi d bezi er_to_power(degree,bez,coeff)
/*Converts Bezier form to power (monomial) form. Works on
one coordinate only.
Input:
Output:

degree:
bez:
coeff:

degree of curve.
coefficients of Bezier form
coefficients of power form.

Remark: For a 2D curve, this routine needs to be called twice,
once for the x-coordinates and once for y.

V
The conversion program internally calls iterated forward differences:

5.8 Implementation

77

void differences(degree,coeffjdiffs)

/*
Computes all forward differences Delta^i(b_0).
Has to be called for each coordinate (x,y, and/or z) of a control polygon.
Input:
degree: length (from 0) of coeff.
coeff: array of coefficients.
Output: diffs: diffs[i]= Delta'^i (coeff [0]).

V

Once the power form is found, it may be evaluated using Horner's scheme:
float

horner(degree,coeff,t)

/*
uses Horner's scheme to compute one coordinate
value of a curve in power form. Has to be called
for each coordinate (x,y, and/or z) of a control polygon.
Input:
degree: degree of curve.
coeff: array with coefficients of curve.
t:
parameter value.
Output: coordinate value.

V

The subdivision routine:
void

subdiV(degree,coeff,weight,t,bleft,bright,wleft,Wright)

/*
subdivides ratbez curve at parameter value t.
Input: degree:
degree of Bezier curve
coeff:
Bezier points (one coordinate only)
weight:
weights for rational case
t:
where to subdivide
Output:
b l e f t , b r i g h t : l e f t and right subpolygons
wleft,wright: their weights

Note:

1. For the polynomial case, set a l l entries in weight to 1.
2. Ordering of right polygon bright is reversed.

*/
Actually, this routine computes a more general case than is described in this
chapter; namely, it computes subdivison for a rational Bezier curve. This will be

78

Chapter 5 The Bernstein Form of a Bezier Curve
discussed later; if the entries in weight are all unity, then wleft and wright will
also be unity and can be safely ignored in the context of this chapter.
Now we present the routine to intersect a Bezier curve with a straight line (the
straight line is assumed to be the x-axis):
void intersect(bx,by,w,degree,tol)
/ * Intersects Bezier curve with x-axis by adaptive subdivision.
Subdivision is controlled by tolerance t o l . There is
no check for stack depth! Intersection points are not found in
'natural' order. Results are written into f i l e o u t f i l e .
Input: bx,by,w:
rational Bezier curve
degree:
i t s degree
tol:
accuracy for results
Output:
intersection points, written into a f i l e

*/
This routine (again covering the rational case as well) uses a routine to check
if a control polygon is flat:
int check_flat(bx,by,degree,tol)
/ * Checks i f a polygon is f l a t . I f all points
are closer than tol to the connection of the
two endpoints, then i t is f l a t . Crashes i f the endpoints
are identical.
Input:

V

Output:

bx,by, degree: the Bezier curve
tol:
tolerance
1 i f f l a t , 0 else.

5.9 Problems
1 Consider the cubic Bezier curve given by the planar control points

11 [-1] 1 ^1 1"-^
0 1L 1 J1111 0
At ^ = 1/2, this curve has a cusp\ its first derivative vanishes and it shows
a sharp corner. You should verify this by a sketch. Now perturb the
x-coordinates of b^ and hi by opposite amounts, thus maintaining a symmetric control polygon. Discuss what happens to the curve.

5.9 Problems

79

2 Show that a nonplanar cubic Bezier curve cannot have a cusp. Hint: use
the fact that b^" , b^~ , bQ are identical v^hen we evaluate at the cusp.
3 Show that the Bernstein polynomial B^ attains its maximum at ^ = i/n. Find
the maximum value. What happens for large n}
* 4 Show that the Bernstein polynomials B^ form a basis for the linear space
of all polynomials of degree n,
PI Compare the run times of decas and hornbez for curves of various degrees.
P2 Use subdivision to create smooth fractals. Start with a degree four Bezier
curve. Subdivide it into two curves and then perturb the middle control
point b2 for each of the two subpolygons. Continue for several levels. Try
to perturb the middle control point by a random displacement and then by
a controlled displacement. Literature on fractals: [35], [411].
P3 Use subdivision to approximate a high-order (n > 2) Bezier curve by a
collection of quadratic Bezier curves. You will have to write a routine
that determines if a given Bezier curve may be replaced by a quadratic one
within a given tolerance. Literature on approximating higher-order curves
by lower-order ones: [336], [341].

This Page Intentionally Left Blank

Bezier Curve Topics

6.1 Degree Elevation
i^uppose we were designing with Bezier curves as described in Section 4.3, trying
to use a Bezier curve of degree n. After modifying the polygon a few times, it may
turn out that a degree n curve does not possess sufficient flexibiUty to model the
desired shape. One way to proceed in such a situation is to increase the flexibility
of the polygon by adding another vertex to it. As a first step, we might want to
add another vertex yet leave the shape of the curve unchanged—this corresponds
to raising the degree of the Bezier curve by one. We are thus looking for a curve
with control vertices b^ ,. . . , b ^ ^ that describes the same curve as the original
polygon b o , . . . , b„.
Using the identities (6.24) to (6.26)—each easy to prove—we rewrite our given
curve as x ( 0 = (1 — t)iL{t) + ^x(0, or

The upper limit of the first sum may be extended to « + 1 since the corresponding
term is zero. The summation of the second sum may be shifted to the limits 1 and
w + 1, and then changed to the lower limit 0 since only a zero term is added. We
thus have

81

82

Chapter 6 Bezier Curve Topics

rb3
Figure 6.1

bI

^ ^ hi

Degree elevation: both polygons define the same (degree three) curve.

Combining both sums and comparing coefficients yields the desired result:

M+ 1

V

(6.1)

M + 1/

Thus the new vertices h- are obtained from the old polygon by piecewise linear
interpolation at the parameter values i/{n + 1). It follov^s that the nev^ polygon
£V lies in the convex hull of the old one. Figure 6.1 gives an example. Note how
SV is "closer" to the curve BV than the original polygon P.
Although our proof is based on straightforward algebraic manipulations, a
more elegant proof is provided through the use of blossoms. If we had the blossom
b(l)[ ^ 1 , . . . , ^„+i] of the degree elevated curve, then we could compute its control
polygon using (4.9). After some experimentation (try the case n = 2!), it is easy
to see that the blossom is given by

h^\h,...,

^„+i] = — - Y

b [ ^ i , . . . , t^^^\t^\

(6.2)

Here, the notation b [ ^ i , . . . , t^^i\tj\ indicates that the argument tj is omitted from
b [ ^ i , . . . , ^«+i]- The control points are now given by application of (4.9):
b|^^ = b^^^[o<"+^-^'^,r^"^].

6.2 Repeated Degree Elevation

83

Inspection of all terms that now arise in (6.2) reveals that the point b^_i appears
/ times and that the point b/ appears n-j-1 — i times, thus reproving our previous
result.^
Degree elevation has important applications in surface design: for several
algorithms that produce surfaces from curve input, it is necessary that these
curves be of the same degree. Using degree elevation, v^e may achieve this by
raising the degree of all input curves to the one of the highest degree. Another
application lies in the area of data transfer between different CAD/CAM or
graphics systems: suppose you have generated a parabola (i.e., a degree two Bezier
curve), and you want to feed it into a system that knows only about cubics. All
you have to do is degree elevate your parabola.

6.2 Repeated Degree Elevation
The process of degree elevation assigns a polygon £V to an original polygon P.
We may repeat this process and obtain a sequence of polygons P, f P, f ^P, and
so on. After r degree elevations, the polygon E^V has the vertices b^j , . . ., b|^^^^,
and each b • is explicitly given by

This formula is easily proved by induction.
Let us now investigate what happens if we repeat the process of degree
elevation again and again. As we shall see, the polygons E^V converge to the
curve that all of them define:
lim E'V = BV.

(6.4)

To prove this result, fix some parameter value t. For each r, find the index /
such that i/{n + r) is closest to t. We can think of i/{n + r) as a parameter on
the polygon f ^P, and as r -> oo, this ratio tends to t. We can now show (using
Stirling's formula) that

lim

:±j;-^t\\-tr-i,

1 Again, work out the example n = 2to build your confidence in this technique!

(6.5)

84

Chapter 6 Bezier Curve Topics

Figure 6.2 Degree elevation: a sequence of polygons approaching the curve that is defined by each
of them.
and therefore

i/{n+r)^t

^

c

'

Equation {6.S) will look familiar to readers with a background in probability: it
states that the hypergeometric distribution converges to the binomial distribution.
Figure 6.2 shows an example of the limit behavior of the polygons £^V,
The polygons £^V approach the curve very slowly; thus our convergence result
has no practical consequences. However, it helps in the investigation of some
theoretical properties, as is seen in the next section.
The convergence of the polygons £^V to the curve was conjectured by A. R. Forrest [244] and proved in Farin [187]. The preceding proof follows an approach
taken by J. Zhou [630]. Degree elevation may be generalized to "corner-cutting";
for a brief description, see Section 8.4.

6.5 The Variation Diminisiiing Property
We can now show that Bezier curves enjoy the variation diminishing property?
the curve BV has no more intersections with any plane than does the polygon P.
Degree elevation is an instance of piecewise linear interpolation, and we know
that operation is variation diminishing (see Section 3.2). Thus each £^V has fewer
intersections with a given plane than has its predecessor £^^~^^V. Since the curve is
The variation diminishing property was first investigated by I. Schoenberg [543] in the
context of B-spline approximation.

6.4 Degree Reduction

85

the limit of these polygons, we have proved our statement. For high-degree Bezier
curves, variation diminution may become so strong that the control polygon no
longer resembles the curve.
A special case is obtained for convex polygons: a planar polygon (or curve) is
said to be convex if it has no more than tw^o intersections with any plane. The
variation diminishing property thus asserts that a convex polygon generates a
convex curve. Note that the inverse statement is not true: convex curves exist
that have a nonconvex control polygon!
Though the variation diminishing property seems straightforward enough, it
is still not totally intuitive. Consider the following statement: two Bezier curves
with common endpoints do not intersect more often than their control polygons.
This appears to be true just after jotting down a few examples. Yet it is false, as
shown by Prautzsch [494].

6.4 Degree Reduction
Degree elevation can be viewed as a process that introduces redundancy: a curve
is described by more information than is actually necessary. The inverse process
might seem more interesting: can we reduce possible redundancy in a curve
representation? More specifically, can we write a given curve of degree n-\-1
as one of degree n} We shall call this process degree reduction.
In general, exact degree reduction is not possible. For example, a cubic with a
point of inflection cannot possibly be written as a quadratic. Degree reduction,
therefore, can be viewed only as a method to approximate a given curve by one
of lower degree. Our problem can now be stated as follows: given a Bezier curve
with control vertices b | ; / = 0 , . . ., w -h 1, can we find a Bezier curve with control
vertices b^; / = 0,. . ., w that approximates the first curve in a "reasonable" way?
The equations for degree elevation may be combined into one matrix equation:
1
•

•
•

bo'

*

r K(1) n
{6.6)

•

*

b„.

L "n+1

1
Abbreviated:
MB = B^i\
M being a matrix with n-\-l rows and n-\-l columns.

(6.7)

86

Chapter 6 Bezier Curve Topics
In degree reduction, we seek to approximate a degee n + 1 curve by one of
degree n. In terms of (6.7), this means that we would be given B^^^ and wish to
find B. Clearly this is not possible in terms of solving a linear system, since M is
not a square matrix.
A trick will help: simply multiply both sides of (6.7) by M^, thus getting
(6.8)

M^MB = AfB^^\

Now we have a linear system for the unknown B with a square coefficient matrix
M^M—and any linear system solver will do the job!^
The linear system (6.8) is called the system of normal equations. It guarantees
that B is optimal in a least squares sense; for more discussion on this technique,
see Section 7.8. It is also optimal in the sense that the original and the degree
reduced curves are as close as possible in the least squares sense; see [406].
In many cases, it might be desired that bg = bQ and b„ = b ^ ^ . Our least
squares solution will not meet these conditions in most cases—the simplest
solution is to enforce them after having found B.
Degree reduction has received a fair amount of attention in the literature; we
cite [97], [182], [406], [472], [612].

6.5 Nonparametric Curves
We have so far considered three-dimensional parametric curves b(^). Now we
shall restrict ourselves to functional curves of the form y = f(x), where f denotes
a polynomial. These (planar) curves can be written in parametric form:
x(t)

Ht) = yit)

t

fit)

We are interested in functions f that are expressed in terms of the Bernstein basis:

Note that now the coefficients bj are real numbers, not points. The bj therefore
do not form a polygon, yet functional curves are a subset of parametric curves
and therefore must possess a control polygon. To find it, we recall the linear
precision property of Bezier curves, as defined by (5.14). We can now write our
3 This linear system is bidiagonal, and may be solved much faster using simple backward
substitution.

6.6 Cross Plots

87

/T\

/

/ \

/

b0

r

-0

1

1/3-

i

1 1

s

1-

1 /I

Z/J

1 1 1

1

Figure 6.3 Functional curves: the control polygon of a cubic polynomial has abscissa values of 0,
1/3,2/3,1.
functional curve as

hit) = J2
7=0

j/n
b,

Blit).

(6.9)

Thus the control polygon of the function f(t) = ^ bjB^ is given by the points
(j/n, bj);j — 0,. . . , w. If we want to distinguish clearly between the parametric
and the nonparametric cases, we call f(t) a Bezier function. Figure 6.3 illustrates
the cubic case. We also emphasize that the bi art real numbers, not points; we
call the bj Bezier ordinates.
Because Bezier curves are invariant under affine reparametrizations, we may
consider any interval [a., b] instead of the special interval [0,1]. Then the abscissa
values are a + i{b — a)/n; / = 0 , . . . , ^.

6.6 Cross Plots
Parametric Bezier curves are composed of coordinate functions: each component
is a Bezier function. For two-dimensional curves, this can be used to construct
the cross plot of a curve. Figure 6.4 shows the decomposition of a Bezier curve
into its two coordinate functions.

88

Chapter 6 Bezier Curve Topics

?Oiorte
[0,1], it
follows that (6.14) is not a barycentric combination of the Cy. In fact, CQ is a point
whereas the other Cy are vectors. The following relations hold:
co = bo,
Cy = Aby_i;

(6.18)
;>0.

(6.19)

This undesirable distinction between points and vectors was abandoned soon
after Forrest's discovery that the Bezier form (6.14) of a Bezier curve could be
written in terms of Bernstein polynomials (see the appendix in [59]).
Comparing both forms, we notice that the Bernstein form is symmetric with
respect to t and 1 — t, whereas the Bezier form is not. Let us assume the defining
coefficients b^ or C/ are affected by some numerical error and then let us check
the effect on the point x(l). In the Bernstein form, x(l) changes its value only if
b„ is in error. In the Bezier form, the value of x(l) is the sum of all errors in the
c^. If those cancel out, no harm is done—but if they do not, we may see serious
error accumulation!

6.9 The Weierstrass Approximation Tiieorem
One of the most important results in approximation theory is the Weierstrass
approximation theorem. S. Bernstein invented the polynomials that now bear
his name in order to formulate a constructive proof of this theorem (see Davis
[133]orKorovkin[364]).
We will give a "customized" version of the theorem, namely, we state it in the
context of parametric curves. So let c be a continuous curve that is defined over
[0,1]. For some fixed n, we can sample c at parameter values i/n. The points
c{i/n) can now be interpreted as the Bezier polygon of a polynomial curve x„:

X„W-^C(-)B^(0.
We say that x„ is the n^"^ degree Bernstein-Bezier approximation to c.

6.10 Formulas for Bernstein Polynomials

91

We are next going to increase the density of our samples, that is, we increase
n. This generates a sequence of approximations x„, x ^ ^ ^ , . . . . The Weierstrass
approximation theorem states that this sequence of polynomials converges to
the curve c:
lim x„(0 = c(0.
n-^OQ

At first sight, this looks like a handy way to approximate a given curve by
polynomials: we just have to pick a degree n that is sufficiently large, and we
are as close to the curve as we like. This is only theoretically true, however. In
practice, we would have to choose values of n in the thousands or even millions
in order to obtain a reasonable closeness of fit (see Korovkin [364] for more
details).
The value of the theorem is therefore more of a theoretical nature. It shows
that every curve may be approximated arbitrarily closely by a polynomial curve.

6.10 Formulas for Bernstein Polynomials
This section is a collection of formulas; some appeared in the text, some did not.
Credit for some of these goes to R. Farouki and V. Rajan [225].
A Bernstein polynomial is defined by

^

I

0

else.

The power basis [f] and the Bernstein basis {B^} are related by

and
n

B"(t) = Y{-\i-'{

\{'\tL

Recursion:

B'i{t) =

(i-t)B';-\t)^tB';:i(t),

(6.21)

92

Chapter 6 Bezier Curve Topics
Subdivision:
B^(CO = ^ B ^ W B ; ( 0 .

(6.22)

Derivative:
±B^(t) =

n[B^:l(t)-B^-\t)l

Integral:
/ B^(x)dx=:
V B"+i(0,
Jo
n + 1 ^^^ ^
Jo

(x)dx ••

(6.23)

1
n-\-l

Three degree elevation formulas:
n +1
tBIit) = l±lBl+lit),
n+1 ^
Blit) = !l±l^B';+\t)
n+1

+ i±lBl+l(t).
n+1 ^

(6.24)
(6.25)
(6.26)

Product:
B-(«)B;(«) = ) ^ B - ) " ( « ) .

6.11

(6.27)

Implementation
A C routine for degree elevation follows. Note that we have to treat the cases
/ = 0 and / = « + 1 separately; the program would not like the corresponding nonexisting array elements. The program actually handles the rational
case, which will be covered later. For the polynomial case, fill wb with I's and
ignore wc.

6.12 Problems

93

void degree_elevdte(bx,by,wb,degree,ex,cy,wc)
/*
input: two-d Bezier polygon in bx, by and with weights
in wb. Degree is degree.
Output:degree elevated curve in cx,cy and with weights in wc.
Note: for nonrational (polynomial) case, fill wc with I ' s .

V

6.12 Problems
*1 Prove (6.17).
* 2 Prove the relationship betv^een the "Bezier" and the Bernstein form for a
Bezier curve (6.14).
3 Prove that
/.
* 4 With the result from the previous problem, prove
f;^(0 = n ( Bpi(x)dx.
PI The recursion formula for Bernstein polynomials is equivalent to the
de Casteljau algorithm. Devise a recursive curve evaluation algorithm for
curves in Chebychev form based on the recursion for Chebychev polynomials. Program it up and experiment!
P2 Program up degree reduction with some of the methods outlined in Section
6.4. Work with the Bezier polygon supplied in the file degred.dat.

This Page Intentionally Left Blank

Polynomial Curve
Constructions

1 olynomial interpolation is a fundamental concept for all of CAGD. Although
its uses are limited to low degrees, the basic concept still needs to be understood
in order to develop new algorithms. If the amount of data is too large for
interpolation to be successful, one uses approximation methods instead.

7.1 Aitken's Aigoritiim
A common problem in curve design is point data interpolation: from data points
p^ with corresponding parameter values f^, find a curve that passes through the
p^.^ One of the oldest techniques to solve this problem is to find an interpolating polynomial through the given points. That polynomial must satisfy the
interpolatory constraints
ip{ti) = p -

/• = 0 , . . . , ^.

Several algorithms exist for this problem—any textbook on numerical analysis
will discuss several of them. In this section, we shall present a recursive technique
that is due to A. Aitken.
We have already solved the linear case, n=\
in Section 3.1. The Aitken
recursion computes a point on the interpolating polynomial through a sequence
The shape of the curve depends heavily on the parameter values t^. Methods for their
determination will be discussed later in the context of spline interpolation; see Section
9.6.
95

96

Chapter 7 Polynomial Curve Constructions

Figure 7.1

Polynomial interpolation: a cubic interpolating polynomial may be obtained as a "blend"
of two quadratic interpolants.
of repeated linear interpolations^ starting with

Let us now suppose (as one does in recursive techniques) that we have already
solved the problem for the case « — 1. To be more precise, assume that we have
found a polynomial pQ~ that interpolates to the n first data points p o , . . . , pn-h
and also a polynomial Pj~ that interpolates to the n last data points p j , . . . , p„.
Under these assumptions, it is easy to write down the form of the final interpolant,
now called p^:
PoW = T ^ P r ' W + —TPT'it)-

(7.1)

Figure 7.1 illustrates this form for the cubic case.
Let us verify that (7.1) does in fact interpolate to all given data points p^: for
^ = ^05

Po(^o) = 1 * Po"^(^o) + 0 * p'l-\to)

= Po.

A similar result is derived for t = tn. Under our assumption, we have pQ~ {fi) =
p\~^{ti) = Pi for all other values of /.
Since the weights in (7.1) sum to one identically, we get the desired p^{ti) = p/.
We can now generalize (7.1) to solve the polynomial interpolation problem:
starting with the given parameter values tj and the data points p/ = p ^ we set

7.1 Aitken's Algorithm

o

•

o

97

o

Figure 7.2 Aitken's algorithm: a point on an interpolating polynomial may be found from repeated
linear interpolation.

It is clear from the preceding consideration that p^(t) is indeed a point on
the interpolating polynomial. The recursive evaluation (7.2) is called Aitken's
algorithm?
It has the follow^ing geometric interpretation: to find p[, map the interval
[tj, tj^j.] onto the straight line segment through p p , p^~^. That affine map takes
t to p^. The geometry of Aitken's algorithm is illustrated in Figure 7.2 for the
quadratic case.
It is convenient to write the intermediate p^ in a triangular array; the cubic
case would look like
Po

^' 1

P2 Pi
P3 pi

2

Po
Pi Po-

(7.3)

We can infer several properties of the interpolating polynomial from Aitken's
algorithm:
Affine invariance: This follows since Aitken's algorithm uses only barycentric combinations.

2 The particular organization of the algorithm as presented here is due to Neville.

98

Chapter 7 Polynomial Curve Constructions
Linear precision: If all p^ are uniformly distributed^ on a straight line
segment, all intermediate p[(^) are identical for r > 0. Thus the straight
line segment is reproduced.
No convex hull property: The parameter t in (7.2) does not have to lie
between ti and tij^y. Therefore, Aitken's algorithm does not use convex
combinations only: PQ(^) is not guaranteed to lie within the convex hull
of the p/. We should note, however, that no smooth curve interpolation
scheme exists that has the convex hull property.
No variation diminishing property: By the same reasoning, we do not
get the variation diminishing property. Again, no "decent" interpolation
scheme has this property. However, interpolating polynomials can augment
variation to an extent that renders them useless for practical problems.

7.2 Lagrange Polynomials
Aitken's algorithm allows us to compute a point p^(^) on the interpolating
polynomial through w + 1 data points. It does not provide an answer to the
following questions: (1) Is the interpolating polynomial unique? (2) What is
a closed form for it? Both questions are resolved by the use of the Lagrange
polynomials L".
The explicit form of the interpolating polynomial p is given by
n

p(0 = ^ P . L f ( 0 ,

(7.4)

where the L^ are Lagrange polynomials

Before we proceed further, we should note that the L" must sum to one in order
for (7.4) to be a barycentric combination and thus be geometrically meaningful;
we will return to this topic later.
3 If the points are on a straight line, but distributed unevenly, we will still recapture the
graph of the straight line, but it will not be parametrized linearly.

7.3 The Vandermonde Approach

99

We verify (7.4) by observing that the Lagrange polynomials are cardinal: they
satisfy
L-(?;) = 5,„,

(7.6)

vv^ith 5^ y being the Kronecker delta. In other words, the /th Lagrange polynomial
vanishes at all knots except at the /th one, w^here it assumes the value 1. Because
of this property of Lagrange polynomials, (7.4) is called the cardinal form of the
interpolating polynomial p. The polynomial p has many other representations,
of course (v^e can rewrite it in monomial form, for example), but (7.4) is the only
form in which the data points appear explicitly.
We have thus justified our use of the term the interpolating polynomial. In
fact, the polynomial interpolation problem always has a solution, and it always
has a unique solution. The reason is that, because of (7.6), the L^ form a basis
of all polynomials of degree n. Thus, (7.4) is the unique representation of the
polynomial p in this basis. This is why one sometimes refers to all polynomial
interpolation schemes as Lagrange interpolation.^
We can now be sure that Aitken's algorithm yields the same point as does (7.4).
Based on that knowlege, we can conclude a property of Lagrange polynomials
that was already mentioned right after {7.5)^ namely, that they sum to 1:

This is a simple consequence of the affine invariance of polynomial interpolation,
as shown for Aitken's algorithm.

7.5 The Vandermonde Approach
Suppose we want the interpolating polynomial p^ in the monomial basis:
n

/=0

4

More precisely, we refer to all those schemes that interpolate to a given set of data points.
Other forms of polynomial interpolation exist and are discussed later.

100

Chapter 7 Polynomial Curve Constructions
The standard approach to finding the unknown coefficients from the known data
is simply to write down everything one knows about the problem:
P^'C^o) = Po = ^0 + ^1^0 + • • • + a„^Q,
P'^ih) = Pi = ao + ai^i + . . . + a^t^,

P"(^«) = P« = ao + ^itn + . . . + ^n^n'

matrix form:
"Po"
Pi

-P«-

~1
1

^0

h

A

_1 tn tl

•

t^

ai

(7.8)

• . t""n -J -a„_

e can shorten this to
p = Ta.

(7.9)

We already know that a solution a to this linear system exists, but one can
show independently that the determinant det T is nonzero (for distinct parameter
values ti). This determinant is known as the Vandermonde of the interpolation
problem. The solution, that is, the vector a containing the coefficients a^, can be
found from

a = T-V

(7.10)

This should be taken only as a shorthand notation for the solution—not as an
algorithm! Note that the linear system (7.9) really consists of three linear systems
with the same coefficient matrix, one system for each coordinate. It is known
from numerical analysis that in such cases the LJJ decomposition of T is a more
economical way to obtain the solution a. This will be even more important when
we discuss tensor product surface interpolation in Section 15.4.
The interpolation problem can also be solved if we use basis functions other
than the monomials. Let {ff }^^o ^^ ^^^'^ ^ basis. We then seek an interpolating
polynomial of the form

p«(o-^cyp;w.
7=0

(7.11)

7.4 Limits of Lagrange Interpolation

101

This reasoning again leads to a linear system (three linear systems to be more
precise) for the coefficients c., this time with the generalized Vandermonde F:

F=
F"oit„)

F1(t„)

••

F"„(h)

••

F"„'it„)j

(7.12)

Since the F^ form a basis for all polynomials of degree n, it follows that the
generalized Vandermonde det F is nonzero.
Thus, for instance, we are able to find the Bezier curve that passes through a
given set of data points: the F^ would then be the Bernstein polynomials B^.

1A

Limits of Lagrange Interpolation
We have seen that polynomial interpolation is simple, unique, and has a nice
geometric interpretation. One might therefore expect this interpolation scheme
to be used frequently; yet it is virtually unknown in a design environment. The
main reason is illustrated in Figure 7.3: polynomial interpolants may oscillate.
The top curve in that figure is the Lagrange interpolant to 21 points read off
from a quarter of an ellipse. The data points were computed to a precision of six
digits. Slightly changing the input data points, namely, by reducing their accuracy
to four digits, produces the bottom interpolant. This is a disturbing phenomenon:
miniscule changes in the input data may result in serious changes of the result.
Processes with that behavior are called ill conditioned. From a more geometric
viewpoint, we may state that polynomial interpolation is not shape preserving.

Figure 7.3 Lagrange interpolation: The top and bottom input data differ only by the amount of
accuracy: six digits after the decimal point, top; four digits, bottom.

102

Chapter 7 Polynomial Curve Constructions
This phenomenon is not due to numerical effects; it is actually inherent in
the polynomial interpolation process. Suppose we are given a finite arc of a
smooth curve c. We can then sample the curve at parameter values t^ and pass
the interpolating polynomial through those points. If we increase the number of
points on the curve, thus producing interpolants of higher and higher degree, we
would expect the corresponding interpolants to converge to the sampled curve
c. But this is not generally true: smooth curves exist for which this sequence of
interpolants diverges. This fact is dealt with in numerical analysis, where it is
known by the name of its discoverer: it is called the Runge phenomenon [513].
Note, however, that the Runge phenomenon does not contradict the Weierstrass
approximation theorem!
As a second consideration, let us examine the cost of polynomial interpolation, that is, the number of operations necessary to construct and then evaluate
the interpolant. Solving the Vandermonde system (7.8) requires roughly n^ operations; subsequent computation of a point on the curve requires n operations.
The operation count for the construction of the interpolant is much smaller for
other schemes, as is the cost of evaluations (here piecewise schemes are far superior). This latter cost is the more important one, of course: construction of the
interpolant happens once, but it may have to be evaluated thousands of times!

7.5 Cubic Hermite Interpolation
Polynomial interpolation is not restricted to interpolation to point data; we can
also interpolate to other information, such as derivative data. This leads to an
interpolation scheme that is more useful than Lagrange interpolation: it is called
Hermite interpolation. We treat the cubic case first, in which one is given a
set of points p/, associated parameter values ^/, and associated tangent vectors
(i.e., derivatives) m^. We just consider the case of two points po, pi and two
tangent vectors mo, m^, setting tQ = 0 and t^ = 1. The objective is to find a cubic
polynomial curve p that interpolates to these data:
P(0) = Po,
p(0) = mo,
p(l) = mi,
P(l)=Pi.
where the dot denotes differentiation.
We will write p in cubic Bezier form, and therefore must determine four Bezier
points b o , . . . , b3. Two of them are quickly determined:

7.5 Cubic Hermite Interpolation

103

Figure 7.4 Cubic Hermite interpolation: the given data—points and tangent vectors—together v^ith
an interpolating cubic.

bo = Po?

b3 = pi.

For the remaining two, we recall (from Section 5.3) the endpoint derivative for
Bezier curves:
p(0) = 3Abo,

p(l) = 3Ab2.

We can easily solve for b^ and hi'u

1

oi = Po + THIO,

1.
b2 = pi -

1
-mi.

This situation—for the case of a general set of points and tangent vectors—is
shown in Figure 7.4.
Having solved the interpolation problem, we now attempt to write it in
cardinal form; we would like to have the given data appear explicitly in the
equation for the interpolant. So far, our interpolant is in Bezier form:

Pit) = poBlit) + L o + ^mo j B\it) + L i - ^mi j Bl(t) + p^B^W-

To obtain the cardinal form, we simply rearrange:
Pit) = PoH^it) + moHlit) + m^Hlit) + PiH^it),

(7.13)

104

Chapter 7 Polynomial Curve Constructions

\\H^\%i\\\\\\\\\\\\\]/^H^\\\
0 \ JX.
III

R3

K

^1 MJ4

MM
4 1

MM

MX

\Jr

3

ill

PV

n\J
l^il 1

nin4^44-U44fTr Hj
M M M M ITM M M ^ i

Figure 7.5 Cubic Hermite polynomials: the four Hf are shown over the interval [0,1].

where we have set^
WQ(t) =

Bl(t)+B\(t),

Hl(t) = h\(t),
(7.14)

Hl(t) = -lBl(t\
Hl(t) = Bl(t)-]-Bl{t).
The Hf are called cubic Hermite polynomials and are shown in Figure 7.5.
What are the properties necessary to make the Hf cardinal functions for the
cubic Hermite interpolation problem? They must be cardinal with respect to
evaluation and differentiation at t = 0 and ^ = 1, that is, each of the Hf equals 1
for one of these four operations and is 0 for the remaining three:

5 This is a deviation from standard notation. Standard notation groups by orders of
derivatives (i.e., first the two positions, then the two derivatives). The form of (7.13)
was chosen since it groups coefficients according to their geometry.

7.5 Cubic Hermite Interpolation

H 3 ( 0 ) = 1,

^ H 3 ( 0 ) = 0,

^ H 3 ( 1 ) = 0,

d^

d^

H | ( 0 ) = 0,

^ H | ( 0 ) = 0,

^^2^1) = !,

H | ( 1 ) = 0,

H|(0)

A H | ( 0 ) = 0,

A H 3 ( 1 ) = 0,

H | ( 1 ) = 1.

= 0,

105

H 3 ( 1 ) = 0,

Another important property of the Hf follows from the geometry of the
interpolation problem; (7.13) contains combinations of points and vectors. We
know that the point coefficients must sum to 1 if (7.13) is to be geometrically
meaningful:

This is, of course, also verified by inspection of (7.14).
Cubic Hermite interpolation has one annoying pecufiarity: it is not invariant
under affine domain transformations. Let a cubic Hermite interpolant be given
as in (7.13), that is, having the interval [0,1] as its domain. Now apply an affine
domain transformation to it by changing t to i = (1 — t)a + tb, thereby changing
[0,1] to some [a, b]. The interpolant (7.13) becomes
Pit) = poH^(?) + moHld) + m^Hld) + PiH|(?),

(7.15)

where the Hf(i) are defined through their cardinal properties:

H|(^)-0,

jMl{a)^0,

j.Hlib)

= l,

Hlib) = 0,

To satisfy these requirements, the new Hf must differ from the original Hf. We
obtain

106

Chapter 7 Polynomial Curve Constructions
j5ti\

_ u3/

Hl{t) =

{b-a)Hl{t),

Hl{i) =

{b-a)Hl(t\

(7.16)

where t e [0,1] is the local parameter of the interval [a^ b\
Evaluation of (7.15) at ? = ^ and i = b yields ^{d) = po, p(b) = p^. The derivatives have changed, hov^ever. Invoking the chain rule, we find that dp(^)/d^ =
(b — a)mQ and, similarly, dp(fc)/d^ = (b — a)mi.
Thus an affine domain transformation changes the curve unless the defining
tangent vectors are changed accordingly—a drawback that is not encountered
with the Bernstein-Bezier form.
To maintain the same curve after a domain transformation, we must change
the length of the tangent vectors: if the length of the domain interval is changed by
a factor a, we must replace mg and m^ by mo/a and m^/a, respectively. There is an
intuitive argument for this: interpreting the parameter as time, we assume we had
one time unit to traverse the curve. After changing the interval length by a factor
of 10, for example, we have 10 time units to traverse the same curve, resulting in
a much smaller speed of traversal. Since the magnitude of the derivative equals
that speed, it must also shrink by a factor of 10.
We also note that the Hermite form is not symmetric: if we replace ^ by 1 — ^
(assuming again the interval [0,1] as the domain), the curve coefficients cannot
simply be renumbered (as in the case of Bezier curves). Rather, the tangent vectors
must be reversed. This follows from the above by applying the affine map to the
[0,1] that maps that interval to [1, 0], thus reversing its direction.
The dependence of the cubic Hermite form on the domain interval is rather
unpleasant—it is often overlooked and can be blamed for countless programming
errors by both students and professionals. We will use the Bezier form whenever
possible.

7.6 Quintic Hermite Interpolation
Instead of prescribing only position and first derivative information at two points,
one might add information for second-order derivatives. Then our data are
Po, mo, SQ and p^, m^, s^, where So and s^ denote second derivatives. The lowestorder polynomial to interpolate to these data is of degree five. Its Bezier points
are easily obtained following the preceding approach. If we rearrange the Bezier
form to obtain a cardinal form of the interpolant p, we find

1.1 Point-Normal Interpolation

107

p(0 = PoH^W + moH^^W + soH|(0 + ^x^\if) + v^x^\(t) + P i H | ( 0 ,

(7.17)

where

^

20 ^

^5
3

1 R5
20 3

It is easy to verify the cardinal properties of the H^^: they are the straightforward
generalization of the cardinal properties for cubic Hermite polynomials. If used
in the context of piecewise curves, the quintic Hermite polynomials guarantee
C^ continuity since adjoining curve pieces interpolate to the same second-order
data. For most applications, one will have to estimate the second derivatives that
are needed as input. This estimation is a very sensitive procedure—so unless the
quintic form is mandated by a particular problem, the simpler C^ cubic splines
presented in Chapter 9 are recommended.

7.7 Point-Normal Interpolation
In a surface generation environment, one is often given a set of points p^ G E^
and a surface normal vector n^ at each data point, as illustrated in Figure 7.6.
Thus we know only the tangent plane of the desired surface at each data point,
not the actual endpoint derivatives of the patch boundary curves.
If we know that two points p/ and py have to be connected, then we must
construct a curve leading from p^ to py that is normal to n^ at p^ and to n. at p..
A cubic will suffice to solve this generalized Hermite interpolation problem. In
Bezier form, we already have bo = p/ and b3 = py. We still need to find b^ and h^.
There are infinitely many solutions, so we may try to pick one that is both
convenient to compute and of reasonable shape in most cases. Two approaches to
this problem appear in Piper [483] and Nielson [447]. Both approaches, although
formulated differently, yield the same result.

108

Chapter 7 Polynomial Curve Constructions

Figure 7.6 Finding cubic boundaries: although the endpoints of a boundary curve are fixed, its end
tangents only have to lie in specified planes.

In order to find b^, project b3 into the plane defined by bo = p/ and n^. This
defines a tangent at bg. As a rule of thumb, the distance ||bi — boll should be
roughly 0.4||b3 — boll; if our current b^ violates this rule, it may have to be
adjusted accordingly. The remaining point hi is then obtained analogously.

7.8 Least Squares Approximation
In many applications, v^e are given more data points than can be interpolated by
a polynomial curve. In such cases, an approximating curve w^ill be needed. Such
a curve does not pass through the data points exactly; rather, it passes near them,
still capturing the shape inherent to the given points. The technique best known
for finding such curves is knov^n as least squares approximation. An example is
given in Figure 7.7.
To make matters more precise, assume wt are given P + 1 data points
p o , . . . , pp, each pi being associated with a parameter value tj. We wish to find a
polynomial curve x(^) of a given degree n such that the distances ||p^ — x(^/) || are
small. Ideally, we would have p^ = x(^/); / = 0 , . . . , P. If our polynomial curve
x(t) is of the form

xw = coqw + ... + c^c^w
for some set of basis functions Cf(t), then we would have

7.8 Least Squares Approximation

Figure 7.7

109

Least squares approximation: data points are sampled from the cross section of an airplane
wing and a quinitic Bezier curve is fitted to them.

oqc^o) + •. • + c„qao)

=Po
Pp.

This may be condensed into matrix form:
Co

Clitp)

qcfp) J

Po

LPPJ

Or, even shorter:
MC = P.

(7.18)

Since we assume the number P of data points is larger than the degree n of
the curve, this hnear system is clearly overdetermined. We attack it by simply
multiplying both sides by M^:
M^MC = M^V.

(7.19)

This is a linear system with w + 1 equations in n-\-1 unknow^ns, with a square
and symmetric coefficient matrix M^M, Its solution is straightforward, provided
that M^M is always invertible. This is in fact the case: since the C" are assumed
to be linearly independent, the « + 1 columns of M are linearly independent, thus
ensuring full rank of M^M,

110

Chapter 7 Polynomial Curve Constructions
But is the solution meaningful? After all, we employed a rather crude trick in
going from (7.18) to (7.19). It turns out that our solution is not only meaningful—
in fact it is optimal.
In order to justify this claim (and also making it more precise), we give a
second derivation of (7.19).
Consider the following expression:
/"(Co,..., c„) = ^

(7.20)

lip,- - x(^,)

/=0

If the curve x(t) does not deviate from the data points p/ by much, then f will
attain a small value. Ideally, if x were to pass through all pj exactly, we would
have /" = 0.
If we substitute the full definition of x(^) into (7.20), we obtain

/ • ( C o , . . . , c^) =

(7.21)

^
/=0

j=0

We wish to find a set of C/ such that the value of f becomes minimal. Since f is
a multivariate function of all components of the C/, that minimum is achieved if
f's partials with respect to all these components vanish:
= 0;

^ = 0, . . . , w ; d = 1,2 or J = 1,2, 3,

where the superscript d labels the individual components of the c^^. Computing
the required derivatives yields
p

E
/=0

n

pf-J^cfCfit,)

Cl(t,) = 0;

k = 0,,..,n;

J = 1,2 or J = 1,2, 3.

'=0

Upon rearranging, we see that this is identical to (7.19)!

7.9 Smoothing Equations
The solution to the least squares problem aims only at minimizing the error
function f in (7.20). It does not "care" about the shape of the resulting curve.
It may wiggle more than we would like, yet it is as close to the data points as
possible. Sometimes wiggles are undesired, even if closeness of approximation is

13 Smoothing Equations

111

Figure 7.8 Least squares approximation: a degree 13 Bezier curve fitted to the airplane wing data
set.
lost. Typically, this is the case for noisy data: if the data points exhibit extraneous
wiggles, there is no point in trying to reproduce these. Figure 1.% shows a
"wiggly" control polygon. The next paragraph shows a way to improve the shape
of the approximating curve, albeit at the cost of deviating more from the given
data points.
We now assume that the basis functions C^^{t) are in fact the Bernstein polynomials W^(t) and the unknown coefficients are contained in a vector B. Instead of
addressing the shape of the curve^ we will simply look at the shape of the control
polygon. Clearly, if the polygon behaves nicely, then so does the curve.
How do we measure the shape of a polygon.^ Many such measures are conceivable; here, we pick one of the simplest ones. If all second differences A^hj are
small, then the polygon does not wiggle much. In addition to (7.18), we might
thus aim at also achieving all the following:
:0

bo — 2bi + b2
b„_2 - 2b„

b.

=0.

This we abbreviate as
5B = 0.

(7.22)

If we simply append these equations to the overdetermined system (7.18), we
obtain one that is even more overdetermined:
M
S

B=

(7.23)

It is solved in the same way as (7.18), that is, by forming the symmetric linear
system of normal equations. The system is solvable because the coefficient matrix
of (7.23) still has n-\-l linearly independent columns, as inherited from the initial
matrix M.

112

Chapter 7 Polynomial Curve Constructions

Figure 7.9 Least squares approximation: a degree 13 Bezier curve fitted to an incomplete airplane
wing data set. For this example, a = 0.1.
In (7.23), we have little control over the effect of the added equations deahng
with the shape of the curve. We gain such control by weighting the two different
components:

[";r]-["Tl-

(7.24)

For a e [0,1), this allows a weighting between the initial data fitting equations
(7.23) and the shape optimizing equations (7.22). For a = 0, we retrieve (7.23);
for values of a near 1, we give more weight to the fitting equations. As a rule
of thumb, the noisier the data, the larger a. If no information about noise is
available, a = 0.1 works well.
The effect of these shape equations is shown in Figure 7.9: the airplane wing
data set was artificially decimated, resulting in a least squares solution that looked
worse than Figure 7.8. Employing shape equations with a = 0.1 results in the
shown solution.

7.10 Designing witii Bezier Curves
According to Bezier, designers at Renault were quickly getting used to manipulating control points of a curve in order to create a particular shape. Other designers
may not like the concept of control points^ and would prefer to manipulate points
on the curve directly.
In that context, a designer would "grab" a point on a curve and designate a
new location for it, intending that the new curve should pass through the changed
location. We now describe one way of doing this, following ideas from Barrels
and Forsey [48].
Suppose we are given a Bezier curve

6

Such as the designers at Mercedes-Benz, where I worked in the 1980s.

7.10 Designing with Bezier Curves

113

n
i=0

We would like to change a point x{i) to a new location y. How do we have to
change the b/ such that the curve passes through y?
We need to find new control points b/ such that

y = X]W^^
i=0

which can be viewed as one equation inn— 1 unknowns, assuming that bo and
b„ should stay the same. This is an underdetermined linear system (for n > 2)
with infinitely many solutions. We would like the one that is the closest to the
original control points.
We write our linear system as
= y-boB^(?)-b,B^(?) = z

[B",ii)...B"„_,ii)]

(7.25)

Lb«-i J
and shorten this to
AB = z
with A = [B^(i)...

(7.26)

B^_-^(i)], We now use a little trick and write our unknowns

as
B = B + A^c

(7.27)

with an unknown vector c. The column vector B contains the current control
points bi through b„_i. This may be viewed as an overdetermined system for c,
which may be solved using a least squares approach:
A B = A B + AA'^C.

Since AB = z, this simplifies to
z - AB = AA'^C.

We observe that AA^ is a scalar and denote it by a. Thus'^
7

Keep in mind that AB = x - hoB^ii) - b„B^(?).

114

Chapter 7 Polynomial Curve Constructions

Figure 7.10 Modifying Bezier curves: the location of a point on the curve is changed.
c=-(y-x)
a
and finally

B = B+-A'^(y-x).
All Bezier points b^,. . . , b„_i move in the same direction y — x, but v^ith
different displacement magnitudes. Figure 7.10 gives an example.

7.11 The Newton Form and Forward Differencing
AH methods in this chapter—and in the Bezier curve chapters as well—^were
concerned with the construction of polynomial curves. We shall now introduce
a way to display or plot such curves. The underlying theory makes use of the
Newton form of a polynomial; the resulting display algorithm is called forward
differencing and is well established in the computer graphics community. For
this section, we deal only with the cubic case; the general case is then not hard
to work out.
So suppose that we are given a cubic polynomial curve p(^). Also suppose
that we are given four points pC^o)? P(^i)? p(h)'> p(h) ^^ i^ such that tj^i — tj = h,
that is, at equally spaced parameter intervals. Then it can be shown that this
polynomial may be written as
1

1

P(^) = Po + ^(* - ^o)Apo + ^ i t - to)it - ti)A^Po
1

(7.28)

The derivation of this Newton form is in any standard text on numerical analysis.

7.11 The Newton Form and Forward Differencing

115

The differences A^py are defined as
A'p,. = A ' - V / + i - A ' - i p y

(7.29)

and A^py = py.
The coefficients in (7.28) are conveniently written in a table such as the
following (setting g = 1/h):
Po
Pi gApo
P2 g^Pl
P3 g^Pl

g^^^Po
g^^^Pl

g^^W

The diagonal contains the coefficients of the Newton form. The computation of
this table is called the startup phase of the forward differencing scheme.
We could now evaluate p at any parameter value t by simply evaluating (7.28)
there. Since our evaluation points tj are equally spaced, a much faster way exists.
Suppose we had computed py = p(^y), and so on, from (7.28). Then we could
compute all entries in the following table:
Po
Pi gApo
P2 g^Pi

g^A^po

P3 g^Pi

g^^^Pi

g^^^PO'

P4

gAp3

g^A2p2

g^A^pi

P5

^Ap4

g^A2p3

g^A^p2

(7-30)

Now consider the last column of this table, containing terms of the form g^A^py.
All these terms are equal! This is so because the third derivative of a third degree
polynomial is constant, and because the third derivative of (7.28) is given by
g^A^P0 = g^A^Pl = . . . .
We thus have a new way of constructing the table (7.30) from right to left:
instead of computing the entry p4 from (7.28), first compute g'^ A^p2 from (7.29):
g^A^P2=g^A^Pl+g^A^Pl,
then compute gAp3 from
gAp3 = g^A^P2+gAp2,

116

Chapter 7 Polynomial Curve Constructions
and finally
P4=gAp3 + P3.
Then compute P5 in the same manner, and so on. The general formula is, with
qj = g^A^Py:
q; = q;:^l + q;_i;

/ = 2,1,0.

(7.31)

It yields the points py = q^^.^
This way of computing the py does not involve a single multiplication after
the startup phase! It is therefore extremely fast and has been implemented in
many graphics systems. Given four initial points po, pi, p2, P3 and a stepsize /?, it
generates a sequence of points on the cubic polynomial through the initial four
points. Typically, the polynomial will be given in Bezier form, so those four points
have to be computed as a startup operation.
In a graphics environment, it is desirable to adjust the stepsize h such that
each pixel along the curve is hit. One way of doing this is to adjust the stepsize
while marching along the curve. This is called adaptive forward differencing and
is described by Lien, Shantz, and Pratt [389] and by Chang, Shantz, and Rochetti
[109].
Although fast, forward differencing is not foolproof. As we compute more
and more points on the curve, they begin to be affected by roundoff. So while we
intend to march along our curve, we may instead leave its path, deviating from it
more and more as we continue. For more literature on this method, see Abi-Ezzi
[1], Bartels, Beatty, and Barsky [47], or Shantz and Chang [571].

7.12

Implementation
The code for Aitken's algorithm is very similar to that for the de Casteljau
algorithm. Here is its header:
float aitken(degree,coeff,t)
/*
uses Aitken to compute one coordinate
value of a Lagrange interpolating polynomial. Has to be called
for each coordinate (x,y, and/or z) of data points.
Input:
degree: degree of curve.

8

It holds for any degree n if we replace / = 2,1,0 by / = « — 1, w — 2 , . . . , 0.

7.13 Problems

Output:

117

coeff: array with coordinates to be interpolated,
t:
parameter value.
coordinate value.

Note: we assume a uniform knot sequence!

V

7.15 Problems
1 Show that the cubic and quintic Hermite polynomials are linearly independent.
2 Generahze Hermite interpolation to degrees 7, 9, and so on.
* 3 The de Casteljau algorithm for Bezier curves has as its "counterpart" the
recursion formula (5.2) for Bernstein polynomials. Deduce a recursion
formula for Lagrange polynomials from Aitken's algorithm.
* 4 The de Casteljau algorithm may be generalized to yield the concept of
blossoms. This is not possible for Aitken's algorithm. Why?
* 5 The Hermite form is not invariant under affine domain transformations,
v^hereas the Bezier form is. What about the Lagrange and monomial forms?
What are the general conditions for a curve scheme to be invariant under
affine domain transformations?
* 6 When P = n in least squares approximation, we should be back in an
interpolation context. Shov^ that this is indeed the case.
PI Aitken's algorithm looks very similar to the de Casteljau algorithm. Use
both to define a whole class of algorithms, of which each would be a special
case (see [205]). Write a program that uses as input a parameter specifying
if the output curve should be "more Bezier" or "more Lagrange."
P2 The function that was used by Runge to demonstrate the effect that now
bears his name is given by
fW = 7 ^ ;

xG[-l,l].

Use the routine ait ken to interpolate at equidistant parameter intervals.
Keep increasing the degree of the interpolating polynomial until you notice
"bad" behavior on the part of the interpolant.
P3 In Lagrange interpolation, each p^ is assigned a corresponding parameter
value tj. Experiment (graphically) by interchanging two parameter values
ti and tj without interchanging p^ and py. Explain your results.

This Page Intentionally Left Blank

B-Spline Curves

D-splines were investigated as early as the nineteenth century by N. Lobachevsky
(see Renyi [506], p. 165); they were constructed as convolutions of certain probability distributions.^ In 1946,1. J. Schoenberg [542] used B-splines for statistical
data smoothing, and his paper started the modern theory of spline approximation. For the purposes of this book, the discovery of the recurrence relations
for B-splines by C. de Boor [137], M. Cox [129], and L. Mansfield was one of
the most important developments in this theory. The recurrence relations were
first used by Gordon and Riesenfeld [284] in the context of parametric B-spline
curves.
This chapter presents a theory for arbitrary degree B-spline curves. The original development of these curves makes use of divided differences and is mathematically involved and numerically unstable; see de Boor [138] or Schumaker
[546]. A different approach to B-splines was taken by de Boor and Hollig [143];
they used the recurrence relations for B-splines as the starting point for the theory. In this chapter, the theory of B-splines is based on an even more fundamental
concept: the blossoming method proposed by L. Ramshaw [498] and, in a different form, by R de Casteljau [147]. More literature on blossoms: Gallier [252],
Boehm and Prautzsch [87].

1 However, those were only defined over a very special knot sequence.

119

120

Chapter 8 B-Spline Curves

8.1

Motivation
B-spline curves consist of many polynomial pieces, offering much more versatility
than do Bezier curves. Many B-spline curve properties can be understood by
considering just one polynomial piece—that is how w^e start this chapter.
The Bezier points of a quadratic Bezier curve may be written as blossom values
b[0,0],b[0,l],b[l,l].
Based on this, we could get the de Casteljau algorithm by repeated use of the
identity t = (l — t)'0-\-t'l.
The pairs [0,0], [0,1], [1,1] may be viewed as being
obtained from the sequence 0 , 0 , 1 , 1 by taking successive pairs.
Let us now generalize the sequence 0 , 0 , 1 , 1 to a sequence WQ, U^, ui-, Uy The
quadratic blossom b[w, u\ may be written as
1 r

n

^2 — ^ 1 r

b[w, u\ = -^

U2-U1
+

i

U —

b[^/i, u\ H

UA

^ ^

^

b[w, U2\

\U2-U0

U2-UQ

U — U^{u^
— U.^
-,
—
b[wi, U2] +
U2 — Ui \W3 — Ui

U — U^^^
\
^ b [ ^ 2 , ^3] •
U^ — Ui
)

This uses the identity
u=

uy — u
U2 — U\

ui H

u — u^

-U2

^2 ~ ^1

for the first step and the two identities
u—

U2 — U

^2 ~ ^0

^0 "^

U — Ur)

^2

^2 ~ ^0

and
U-i — U

U — U\

^3 — ^1

^3 — ^1

u = —^

Ux H

/

-u^

8.1 Motivation

121

for the second step. Note that we successively express u in terms of intervals of
growing size.
Starting with the b[w^, w^.^^] as input control points, we may rewrite this as:

b[w2? ^3] b[w, t/2] b[w, u\.
This is a first instance of the de Boor generalization of the de Casteljau
algorithm. See Example 8.1 for a detailed computation.
Figure 8.1 illustrates the algorithm, but using the knot sequence WQ? ^l? ^2? ^3 =
0 , 1 , 3,4 and ^ = 2.0.

Example 8.1

The de Boor algorithm for n = l.
Let t/Q, ^1, ^2? ^3 = 0,2,4, 6. Let the control points be given by
h[uQ, ui] =

b[wi, U2] =

b[w2, u^] =

0

Setting w = 3, we now compute the point b[3, 3]. At the first level, we compute
two points
b[2,3] =

4-3
4-0

b[3,4] =

6 - 3 [8] 3 - 2 [8] [8]
— _6_
6 - 2 _8_ + 6 3 1 _0_

b[3,3]:

4-3
4-2

3-0

+4-0

[t\

and

Finally,
3-2

+4-2

7
6

122

Chapter 8 B-Spline Curves

h[u2,u^]

1
UQ

1
Ui

____^___

Figure 8.1

•
u

1
«2

1
^3

ul

The de Boor algorithm: the quadratic case.

8.2 B-Spline Segments
B-spUne curves consist of a sequence of polynomial curve segments. In this
section, w^e focus on just one of them.
Let U be an interval [t/j, t^/+i] in a sequence {wj of knots. We define ordered
sets W- of successive knots, each containing ui or Ui^i. The set W- is defined such
that:
W- consists of r + 1 successive knots.
ui is the (r — /)th element of U[, with / = 0 denoting the first of I7['s elements.
We also observe

u; = u[+inu[+/.
When the context is unambiguous, we also refer to the U[ as intervals^ having
the first and last elements of W- as endpoints. In that context, l]\ = 17. We also
define U = [Ug, Uj] and use the term interval only if U^^U^,
A degree n curve segment corresponding to the interval U is given by « + 1
control points d^ which are defined by

8.2 B-Spline Segments
d,=b[Uf-l];

/ = 0,...,«.

123
(8.1)

The point x(u) = h[u^^^] on the curve is recursively computed as
d'.(u) = b[^<'^, Uf " ^ " 1 ;

r=l,,..,n;i

= 0,...,n-r

(8.2)

with x(u) = d-Qiu)? This is knov^n as the de Boor algorithm after Carl de Boor
see [137]. See Example 8.2 for the case n = 3 and Figure 8.2 for an illustration.
Equation (8.2) may alternatively be written as
d;(«) = (1 - C r ^ d p l ^ C r M ; ; l ;

r=l,...,n;i

= 0,...,n-r,

(8.3)

where t^7('^ is the local parameter in the interval Uf~-[~^ .
A geometric interpretation is as follows. Each intermediate control polygon
leg dp d^-^-^ may be viewed as an affine image of Uf_^/^ . The point d^^ is then
the image of u under that affine map.
For the special knot sequence 0^"^, 1^^^ and U = [0,1], the de Boor algorithm
becomes
d^(«)=b[«->,0<«—'>,l<-];

r =!,...,«;

i^O,...,n-r,

(8.4)

which is simply the de Casteljau algorithm.
If the parameter u happens to be one of the knots, the algorithm proceeds as
before, except that we do not need as many levels of the algorithm. For example,
if a quadratic curve segment is defined by h[uQ, u^], b[wi, ^2]? ^Wi^ ^3] ^^^ we
want to evaluate at u = ^2? then two of the intermediate points in the de Boor
algorithm are already known, namely, h[ui, U2] and b[t/2, ^3]. From these two,
we immediately calculate the desired point h[u2, ^2]? thus the de Boor algorithm
now needs only one level instead of two.
Derivatives of a B-spline curve segment are computed in analogy to the Bezier
curve case (5.17)
x(u) = nh[u^''-^^,ll

(8.5)

Expanding this expression and using the control point notation, this becomes

x(«) = ^ ( d r i - d ^ - i ) ,
2 This notation is different from the one used in previous editions of this book.

(8.6)

124

Chapter 8 B-Spline Curves

Example 8.2

The de Boor algorithm for « = 3.
Let part of a knot sequence be given by
. . . ^ 3 , ^ 4 , ^ 5 , ^ 5 , Wy, Wg? • • •

and let U = [^5, w^]. The standard blossom computation of b[w, w, u\ proceeds as
follows:
b[W4, ^5, U^] h[u, U4, Us]
h[us, u^, uj]

h[u. Us, u^] h[u, u, us\

h[u^,uj^ug]

h[u^u^^uj\ h[u,u,u^]

h[u,u,u].

We now write this as

h[Uf]

h[u,Ul]

b[Uf] h[u,Ul]

h[u,u,U^^]

HUJ]

h[u,u,U^]

h[u,Ul]

h[u,u,u].

In terms of control points:

dl 4
di

d\

d3

dl

dl
dl dl

The labeling in the first scheme depends on the subscripts of the knots, whereas
the last two employ a relative numbering.

where \U\ = U^ — UQ denotes the length of the interval 17. Thus the last two
intermediate points dQ~^ and d^~ span the curve's tangent, in complete analogy
to the de Casteljau algorithm.
Higher derivatives follow the same pattern:

/ - X ( « ) = -Jll—h[u<"-r>,
du^
{n — r)i

ll

(8.7)

i.2 B-Spline Segments

125

b[W2,W3,W4]

h[u^,U4,Us]

I

UQ

m-

h[UQ,Ui,U2\

U^

m

U2

U

W3

M4

W5

ui.

ul

Figure 8.2 The de Boor algorithm: a cubic example. The solid point is the result h[u, u, u]; it is on
the line through b[w, w, U2] and b[w, w, W3].

In the case of Bezier curves, we could use the de Casteljau algorithm for curve
evaluation, but we could also write a Bezier curve explicitly using Bernstein
polynomials. Since we changed the domain geometry, we will now obtain a
different explicit representation, using polynomials^ P":

x(^) = ^d,Pf(^).

(8.8)

/=o
The polynomials P" satisfy a recursion similar to the one for Bernstein polynomials, and the following derivation is very similar to that case:

3 These will later become building blocks of B-splines.

126

Chapter 8 B-Spline Curves
n-\

xiu) = J2djP^-\u)
i=0
n—\

n—\

= E d - ti^)d,pr\^)+E tfd^pth")
i=0

i=l

For the first step, we used the de Boor algorithm, letting tf be the local parameter
in the interval U^^^ For the second step, we used the convention P^~^ = 0 to
modify the first term. Using a similar argument: P^ = 0, we may extend the
second term to start with / = 0. We conclude
P^(u) = (1 - tl^)P'^-\u)

+ tfP'Izliu).

(8.10)

This recursion has to be anchored in order to be useful. This is straightforward
for the case n = l:

^,

phu) =
ov

1^1

^,

P](u) =
IV

1^1

where |U| = Uj - Ug. For the special knot sequence 0<"^, 1<"> and U = [0,1],
this is the Bernstein recursion.

8.5 B-Spline Curves
A B-spline curve consists of several polynomial curve segments, each of which
may be evaluated using a de Boor algorithm. A B-spline curve is defined by
the degree n of each curve segment,
the knot sequence UQ, ... ^ Uf^^ consisting of iC + 1 knots Uj < /^/+i,
the control polygon d g , . . . , d^ with L = K ~ n + 1.
Some comments: the numbering of the control points d^ in this definition is
global, whereas in Section 8.2 it was local relative to an interval U. Each dj
may be written as a blossom value with n subsequent knots as arguments. Hence
the number L + 1 of control points equals the number of w-tuples in the knot
sequence.

8.3 B-Spline Curves

127

Example 8.3 Some examples of B-spline curve definitions.
Let n= 1^ and let the knot sequence be 0,1,2, hence K = 2. There will be control
points do, d^, d3. The curve's domain is [UQ^ U^\ and there are two linear curve
segments.
Let n = l with the knot sequence 0, 0.2, 0.4, 0.5, 0.7,1, hence K = S. There will
be control points do, d^, d2, d3, A^ and three quadratic curve segments. If we now
change the knot sequence to 0, 0.2, 0.45, 0.45, 0.7,1, then the number of curve
segments will drop to two.

Each knot may be repeated in the knot sequence up to n times. In some cases it
is approriate to simply list those knots multiple times. For other applications, it is
better to list the knot only once and record its multiplicity in an integer array. For
example, the knot sequence 0.0,0.0,1.0,2.0,3.0,3.0,4.0,4.0 could be stored as
0.0,1.0,2.0, 3.0,4.0 and a multiplicity array 2 , 1 , 1 , 2 , 2 .
There is a different de Boor algorithm for each curve segment. Each is "started"
with a set of U^", that is, by sequences o{n-\-l knots. In order for local coordinates
to be defined in (8.3), no successive w + 1 knots may coincide.
In Section 8.2, we assumed that we could find the requisite V^ for each interval
U. This is possible only if U is "in the middle" of the knot sequence; more
precisely, the first possible de Boor algorithm is defined for U = [u^-i-, u^ and
the last one is defined for U = WK-n-> ^X-w+il = Wx-m ^LI- ^ ^ ^^^^ ^^11 Wn-h ^LI
the domain of the B-spline curve. A B-spline curve has as many curve segments
as there are nonzero intervals U in the domain. Example 8.3 illustrates these
comments.
For more examples of B-spline curves, see Figure 8.3.
Since a B-spline curve consists of a number of polynomial segments, one might
ask for the Bezier form of these segments. For a segment U = [uj, uj^i] of the
curve, we simply evaluate its blossom b ^ and obtain the Bezier points h^,. . ., b ^
as
j^ =b

[Uj

, Uj^^ J.

An example is shown in Figure 8.4. Several constituent curve pieces of the same
curve are shown in Figure 8.5.
When dealing with B-spline curves, it is convenient to treat it as one curve, not
just as a collection of polynomial segments. A point on such a curve is denoted
by d(w), with u e [u^-i-, uj^_^_^i]. In order to evaluate, we perform the following
steps:

128

Chapter 8 B-Spline Curves

\

^

00123456789

n=3

n=3

n=5

n=9
000000000

555555555

Figure 8.3 B-spline curves: several examples.

1. Find the interval U = [uj, ui_^i) that contains u,
2. Find the n-\-l control points that are relevant for the interval U, They are,
using the global numbering, given by dj_„_^i,..., d/.^!.
3. Renumber them as d o , . . . , d„ and evaluate using the de Boor algorithm
(8.3).
In terms of the global numbering of knots, w^e observe that the intervals U^
from the previous section are given by

8.3 B-Spline Curves

Figure 8.4

129

Conversion to Bezier form: the Bezier points of a segment of a cubic B-spline curve.

Figure 8.5 Individual curve segments: the first four cubic segments of a cubic B-sphne curve are
shown, ahernating between dashed and black.

The steps in the de Boor algorithm then become
d^iu) = (1 - af)d^-\u)

+ afd^-/(«)

with
k^

^ ^IM ~

^l-k+i
^l_k+i

for ^ = r + 1 , . . ., w, and / = 0 , . . ., w — ^. Here, r denotes the multiplicity of u,
(Normally, u is not already in the knot sequence; then, r = 0.)
The fact that each curve segment is only affected by w + 1 control points is
called the local control property.
We also use the notion of a B-spline blossom d [ t ' i , . . . , f „], keeping in mind
that each domain interval U has its ov^n blossom b ^ and that consequently
d [ ^ ' i , . . . , f „] is piecewise defined.
B-spline curves enjoy all properties of Bezier curves, such as affine invariance,
variation diminution, etc. Some of these properties are more pronounced now
because of the local control property. Take the example of the convex hull

1 30

Chapter 8 B-Spline Curves

Figure 8.6

The local convex hull property: top, quadratic; bottom, cubic.

Figure 8.7

The local control property: changing one control point of a cubic B-spline curve has only
a local effect.

property: the curve is in the convex hull of its control polygon, but also each
segment is in the convex hull of its defining n-\-l control points; see Figure 8.6.
A consequence of the local control property is that changing one control point
w^ill only affect the curve locally. This is illustrated in Figure 8.7.

8.4 Knot Insertion
Consider a B-spline curve segment defined over an interval U. It is defined by
all blossom values d[[7"~ ]; / = 0 , . . . , « w^here each «-tuple of successive knots
U^ contains at least one of the endpoints of U. If we now split U into two
segments by inserting a new knot ii, the curve will have two corresponding

8.4 Knot Insertion
Example 8.4

131

Knot insertion.
In Figure 8.8, just one de Boor step is carried out for the parameter value u. The
two new resulting points, in blossom notation, are b[wi, u] and b[ii, U2\. Let us
now consider the points

These are the B-spline control points of our curve b[w, u] for the interval [wj, u\\
Similarly, the B-spline control points for the interval [ii, Uj} are given by

h[u,U2]

bK,Wi
h[U2,U^]

UQ

Ui

U

U2

H

^3

Figure 8.8 Knot insertion: a quadratic example.

segments. What are the control points for these two segments.^ The answer is
surprisingly simple: all blossom values d[l7^"~^]; / = 0 , . . . , w + 1 where each ntuple of successive knots U^~^ contains at least one of the endpoints of U. This
result is due to W. Boehm [68], although it was not originally derived using
blossoms. See Example 8.4.
Knot insertion works since B-spline control points are nothing but blossom
values of successive knots—now they involve the new knot u. We may also view
the process of knot insertion as one level of the de Boor algorithm, as illustrated
in Example 8.4.

132

Chapter 8 B-Spline Curves

Figure 8.9

Chaikin's algorithm: starting with a (closed) control polygon B-spline curve, two levels of
the algorithm are shown.
An interesting application of repeated knot insertion is due to G. Chaikin
[105]. Consider a quadratic B-spline curve over a uniform knot sequence. Insert
a new^ knot at the midpoint of every interval of the knot sequence. If the "old"
curve had control vertices d/ and those of the nev^ one are d^- , it is easy to show^
that

4 L = ^d, + ld,_i and d(j> = ^d, + ld.+i.
If this procedure is repeated infinitely often, the resulting sequence of polygons
w^ill converge to the original curve, as follov^s from our previous considerations.
Figure 8.9 show^s the example of a closed quadratic B-spline curve; two levels of
the iteration are shov^n.
Chaikin's algorithm may be described as corner cutting: at each step, vs^e chisel
away the corners of the initial polygon. This process is, on a high level, similar
to that of degree elevation for Bezier curves, w^hich is also a convergent process.
One may ask if corner-cutting processes v^ill alw^ays converge to a smooth curve.
The answer is yes^ with some mild stipulations on the corner-cutting process, and
was first proved by de Boor [140]. One may thus use a corner-cutting procedure
to define a curve—and only very few of the curves thus generated are piecewise
polynomial! Recent work has been carried out by Prautzsch and Micchelli [495]
and [426], based on earlier results by de Rham [150], [151].
R. Riesenfeld [508] realized that Chaikin's algorithm actually generates uniform quadratic B-spline curves. A general algorithm for the simultaneous insertion of several knots into a B-spline curve has been developed by Cohen, Lyche,
and Riesenfeld [121]. This so-called Oslo algorithm needs a theory of discrete
B-splines for its development (see Barrels, Beatty, and Barsky [47]).

8.5 Degree Elevation

133

8.5 Degree Elevation
We may degree elevate in (almost) the same way we could degree elevate Bezier
curves using (6.2). The difference: a given B-spline is a piecewise degree n
curve over a given knot sequence. Its differentiability is determined by the knot
multiplicities. If we write it as a piecewise degree w + 1 curve, we need to
increase the multiplicity of every knot by one, thus maintaining the original
differentiability properties. For example, if we degree elevate a C^ piecewise
linear curve to piecewise quadratic, it is still C^. But for a piecewise quadratic
to be C^, it has to have double knots. Let us denote the knots in this augmented
knot sequence by ii^.
Let V" be a sequence of « + 1 real numbers i/^,..., t'^+i- Let y^\vi denote the
sequence V" with the value Vi removed. Then the degree n-\-l blossom b may be
expressed in terms of the degree n blossom b via

b[y(-+i)]= - ^

(b[V("+^Vi] + . . . +b[y^^+iV.+i]) .

W+ 1 V

/

(8.11)

The proof is identical to that for degree elevation of Bezier curves. The control
points are then recovered from the blossom as before (see Example 8.5).
The inverse process—degree reduction is more important for practical applications. Following the example of the analogous Bezier case, we write the
elevation process as a matrix product and invert it by a least squares technique
for the reduction process; see Section 6.4. This method is described in detail in
[617]. Other methods exist, see [481] and [624].

Example 8.5 B-spline degree elevation and blossoms.
Let a cubic B-spline curve be defined over {^Q = ui = u^^ ^ 3 , . . .}. Then the
interval \u^^ u^] corresponds to [iiy, u^]. We denote the corresponding blossoms
by d^la^fe,c] and d-jla,fe,c, d]. The new control point d^ is computed as follows:
d4 =

dy[u4,us,u^,uy]

\
= - (d4[^4. Us, u^] + d4[^4. Us, uj] + d4[^4, u^, uj] + d^ius, u^, uj^
\
= - (d4[^3. ^3. ^4] + d4[w3, ^4, U4]) .
For the last step, we have used U4 = us = W3 and u^ = uj = u^.

1 34

Chapter 8 B-Spline Curves

8.6 Greville Abscissae
Let l[u] = uht the blossom of the (nonparametric) hnear function u. If we want
to write this linear blossom as a quadratic one: /^[w, v] = l[u], we easily see that

i(2^[u,v]=h[u] + h[v]
gives the desired quadratic form of our linear blossom. If we asked for a cubic
form /^^^[w, V, w] of l[u\ we find that
P\u,

1
V, w] = -f[v,

1
1
w] + -f[u, w] + -fi[u,

v\

If we denote a degree n version of the linear blossom by /(")[V"] with V" =
1^1,..., !/„, it follows that
/(")[V«]=-(t;i+ . . . + £/„).
n
The proof is by induction and was anchored by the earlier examples. The
inductive step starts with the degree elevation formula (8.11):"^

n-\-\ ^

n

This is easily transformed to

n-\-l
thus finishing the proof.
If we are given a knot sequence UQ^. .. ,Uj^ and a degree w, then we know that
any B-spline function d(u) has control vertices d[uQ,..., w„_i],...,  ^h ^l? ^2? ^2? ^2They are
= Co, or, in terms of the associated blossoms, if b[^i, w^, Ui] =
c[ui^ Ui, ui]. Two such curves are shown in Figure 8.14.
The two curves are C^ if they may be written as a B-spline curve with a double,
not a triple knot u^. Then our triple knot at u^ is the result of knot insertion and
the three points b2, b3, c^ are collinear and in the ratio Ao : A | with AQ = UI — UQ
5 The osculant of order r of an n^^ degree polynomial curve x(u) at paramter value w is the
degree r polynomial that agrees w^ith x for all derivatives up to order r.

1 38

Chapter 8 B-Spline Curves

UQ

U2

UQ

Figure 8.14

Smoothness of Bezier curves: the C^ case.
and Ai = U2 — u^. In terms of blossoms:
Ci = h[ui, wj, ui]

and

hi =

C[UQ, W^,

ui\.

For C^ smoothness, the knot ui must have been the resuh of two knot
insertions. It follows that
b[wo, ^1, U2] = c[uo, ux, U2].

(8.12)

is our desired C^ condition. It is illustrated in Figure 8.15.
If we are to check if two given Bezier curves are C^ or not, all we have to do is
construct the two points appearing in (8.12). If they disagree, as in Figure 8.16,
we conclude that the given curve is not C^.
In most practical cases, a C^ check would have to check for approximate
satisfaction of (8.12), since reals or floats are rarely equal. In other words, a
tolerance has to be used. The practical value of (8.12) lies in the fact that it is
amenable to using a point tolerance that determines when two distinct points
are to be considered the same point. Checking for C^ smoothness by comparing
second derivatives would require a different, and less intuitive, tolerance.

8.7 Smoothness
b[Mo,Wi,W2]

Figure 8.15

1

1

1

UQ

Wi

U2

UQ

Ux

U2

Smoothness of Bezier curves: the C^ case.

€>0

h[Uo,Ux,U2]
C[UQ,UX,U2]

U2

Figure 8.16

UQ

U2

UQ

U2

Smoothness of Bezier curves: the C^ condition is violated.

1 39

140

Chapter 8 B-Spline Curves

8.8 B-Splines
Consider a knot sequence ^Q? • • • ? ^M ^^^ ^^e set of piecewise polynomials
of degree n defined over it, where each function in that set is n — TJ times
continuously differentiable at knot Uj. All these piecewise polynomials form a
linear space, with dimension
M-l

dim=(w + l ) + Y^Vi,

(8.13)

i=l

For a proof, suppose we want to construct an element of our piecewise polynomial linear space. The number of independent constraints that we can impose
on an arbitrary element, or its number of degrees of freedom, is equal to the
dimension of the considered linear space. We may start by completely specifying
the first polynomial segment, defined over [UQ, U^]; we can do this in w + 1 ways,
which is the number of coefficients that we can specify for a polynomial of degree
n. The next polynomial segment, defined over [wj, ui], must agree with the first
segment in position and n — r^ derivatives at u^, thus leaving only r^ coefficients
to be chosen for the second segment. Continuing further, we obtain (8.13).
We are interested in B-spline curves that are piecewise polynomials over the
special knot sequence [w„_i, ui]. The dimension of the linear space that they form
is L + 1, which also happens to be the number of B-spfine vertices for a curve in
this space. If we can define L + 1 linearly independent piecewise polynomials in
our linear function space, we have found a basis for this space. We proceed as
follows.
Define functions N^{u), called B-splines by defining their de Boor ordinates
to satisfy dj = 1 and dj = 0 for all / ^ i. The N^(u) are clearly elements of the
linear space formed by all piecewise polynomials over [w„_i, w^]. They have local
support:
Nf(u) 7^ 0 only Hue [w,_i, w,+„].
This follows because knot insertion, and hence the de Boor algorithm, is a local
operation; if a new knot is inserted, only those Greville abscissae that are "close"
will be affected.
B-splines also have minimal support: if a piecewise polynomial with the same
smoothness properties over the same knot vector has less support than N^", it
must be the zero function. All piecewise polynomials defined over [w^_i, w^_^„],
the support region of N^, are elements of a function space of dimension In + 1,
according to (8.13). A support region that is one interval "shorter" defines a
function space of dimension In, The requirement of vanishing n — rj_i derivatives
at Ui_i and of vanishing n — r/_^„ derivatives at Uj_^^ imposes 2n conditions on

8.8 B-Spiines

141

any element in the linear space of functions over [w/_i, Uj_^^_i\. The additional
requirement of assuming a nonzero value at some point in the support region
raises the number of independent constraints to 2n + 1, too many to be satisfied
by an element of the function space with dimension 2n,
Another important property of the N " is their linear independence. To demonstrate this independence, we must verify that
L

^CjN1{u)

(8.14)

=0

implies Cj = 0 for all;. It is sufficient to concentrate on one interval [wj, ui_^i] with
Ui < ui^i. Because of the local support property of B-splines, (8.14) reduces to
/+1

J2

^j^^i^) = 0

for w G K , wj+i].

j=I-n-\-l

We have completed our proof if we can show that the linear space of piecewise
polynomials defined over [w/_„, Uj^^^i] does not contain a nonzero element that
vanishes over [w/, w/_^i]. Such a piecewise polynomial cannot exist: it would have
to be a nonzero local support function over [w/_|_i, wj_^„+i]. The existence of such
a function would contradict the fact that B-splines are of minimal local support.
Because the B-splines N^ are linearly independent, every piecewise polynomial
s over [«„_!, ui] may be written uniquely in the form
L

(8.15)

s(u) = J2diNJ{u),
j=o

The B-splines thus form a basis for this space. This reveals the origin of their
name, which is short for Basis splines. Figure 8.17 gives examples of some cubic
B-splines.
If we set all 4 = 1 ii^ (8.15), the function s(u) will be identically equal to 1,
thus asserting that B-splines form a partition of unity.
Nl Nl

Figure 8.17

B-splines: some cubic examples.

Nl

Nl N^

142

Chapter 8 B-Spline Curves
B-spline curves are simply the parametric equivalent of (8.15):
L

x(u) = J2djNf{u),
;=0

Just as the de Casteljau algorithm for Bezier curves is related to the recursion
of Bernstein polynomials, the de Boor algorithm yields a recursion for B-splines.
It is given by
^"^^-^

Nf(u) =

+ ^l^±^ZiiN-i(^),

N^Hu)

(8.16)

w^ith the "anchor" for the recursion being given by
1

"?<«'

if Ui_i ^X
..„

+-

•

.

ffM

-

4.

-

...

y^.^

*•,
("•^

- - ^.- '

4-

-

Figure 8.18 The B-spline recursion: top, two linear B-splines yield a quadratic one; bottom, two
quadratic B-splines yield a cubic one.

8.9 B-Spline Basics

143

A comment on end knot multiplicities: the widespread data format IGES uses
two additional knots at the ends of the knot sequence; in our terms, it adds knots
u_i and ^L+2n-l- The reason is that formulas like (8.16) seemingly require the
presence of these knots. Since they are multiplied only by zero factors, their values
have no influence on any computation. There is no reason to store completely
inconsequential data, and hence the "leaner" notation of this chapter.

8.9 B-Spline Basics
Here, we present a collection of the most important formulas and definitions of
this chapter. As before, n is the (maximal) degree of each polynomial segment,
L + 1 is the number of control points, and K is the number of intervals.
Knot sequence: {WQ, . . . , Uf^},
Control points: d o , . . . , d^, with L = K -n + 1.
Domain: Curve is only defined over [u^-\^...,
Greville abscissae: ^i = ^(Uj-\

u^],

h ^/-f^-i).

Support: N^ is nonnegative over [w/_i, ^/+„].
Knot insertion: To insert uj ^1-2>^-^K-\)

= ^e

with
rs = P i - d i N ^ ( r i )

and

r^ = P K - 1 - dL-iN^_i(rK-i).

Each of the remaining unknowns d 2 5 . . . , dj^_2 is related to the data points by
one equation. Because of the local support property of cubic B-splines, it is of
the form
p, = d,N3(r,) + d,+iN^i(r,) + Ai^rn]^-^{x>); i = 2,...,K-2.

(9.11)

Together, these are iC — 1 equations for the K — 1 unknown control points. In
matrix form, we have

N|(T2)

Ni_3(rK_2)

Ni_2(TK_2)

H-l(-'K-2)

P2
(9.12)
L^L-2

PK-2

L r^

Schematically, the case K = 5 looks like this:
d2

*

•

•

•

•

•

•

•

•

•

ds
d4

Ids]

ts

P2
P3

.^e]

The entries in this tridiagonal matrix are easily computed from the definitions
of the cubic B-splines Nf,

9.5 More End Conditions

157

In the case of uniform knots Ui = /, the interpolation conditions (9.11) take
on a particularly simple form:
6p,- = d, + 4d,+i + d,+2;

/ = 2 , . . . , X - 2.

(9.13)

We conclude with a method for cubic spline interpolation that occasionally
appears in the literature (e.g., in Yamaguchi [621]). It is possible to solve the
interpolation problem without setting up a linear system! Just do the following:
define do, d^, dj^_i, d^^ as before, and set d/ = p/_i for all other control points.
Then, for / = 2 , . . . , L — 2, correct the location of d^ such that the curve passes
through p^_i. Repeat until the solution is found.
This method will always converge and will not need many steps in order to
do so. So why bother with linear systems? The reason is that tridiagonal systems
are most effectively solved by a direct method, whereas the earlier iterative
method amounts to solving the system via Gauss-Seidel iteration. So though
geometrically appealing, the iterative method needs more computation time than
the direct method.

9.5 More End Conditions
For cubic spline interpolation, the choice of end conditions is important for the
shape of the interpolant near the endpoints—they do not matter much in the
interior. We have seen a clamped end condition earlier. It works well in situations
where the end derivatives are actually known. But in most applications, we do
not have this knowledge. Still, two extra equations are needed in addition to the
basic interpolation conditions (9.10).
We list several possibilities:
The natural end conditions are derived from the physical analogy of a wooden
beam that is clamped at some positions; see Figure 9.22. Beyond the first and last
clamps, such a beam assumes the shape of a straight line. A line is characterized
by having a zero second derivative, and hence the end conditions
x(ro) = 0,

x(rx) = 0

are called natural end conditions.
The second derivative of a Bezier curve at bo is given by bo — 2bi + b2. In
terms of B-spline control vertices (using triple end knots), this becomes

do - 2di + - 4 ^ d i + - 4 ^ d 2 = 0,
AA+AI

An +

A,

158

Chapter 9 Constructing Spline Curves
where we set A^ = r^+i — r^ We rearrange and obtain
(Ao + Ai)do - (2Ao + Ai)di + Aod2 = 0.

(9.14)

A similar condition holds at the other endpoint. Unless required by a specific
application, this end condition should be avoided since it forces the curve to
behave linearly near the endpoints.
Typically another end condition, called Bessel end condition, yields better
results. It works as follows: the first three data points and their parameter values
determine an interpolating quadratic curve. Its first derivative at po is taken to
be the one for the spline curve. If we set
—

a=

a n d j6 = 1 — of

^2 - "^0

and
1

lap
then
2
1
d i = - ( a p o + ^a) + -po.
This value for dj is then used in the earlier clamped end condition. The control
point d£^_i is obtained in complete analogy; it is then also used for a clamped
end condition.
Other end conditions exist: requiring x(ro) = x(ri) is called a quadratic end
condition. If all data points and parameter values were read off from a quadratic
curve, then this condition would ensure that the spline interpolant reproduces
that quadratic. The same is true, of course, for the Bessel end conditions. If
the first and second cubic segment are parts of one cubic, then their third
derivatives at X\ would agree. The corresponding end condition is called nota-knot condition since the knot r^ does not act as a breakpoint between two
distinct cubic segments.
We finish this section with a few examples, courtesy T. Foley, using uniform parameter values in all examples.^ Figure 9.8 shows equally spaced data
5

Owing to the symmetry inherent in the data points, all parametrizations discussed later
yield the same knot spacing. All circle plots are scaled down in the }/-direction.

9.5 More End Conditions

159

Figure 9.8 Exact clamped end condition spline.

K=l

K=0

Figure 9.9 Curvature plot of exact clamped end condition spline.

points read off from a circle of radius 1 and the cubic spline interpolant
tained with clamped end conditions, using the exact end derivatives of
circle. Figure 9.9 show^s the curvature plot^ of the spline curve. Ideally,
curvature should be constant, and the spline curvature is quite close to
ideal.

6

The graph of curvature versus arc length; see also Chapter 23.

obthe
the
this

160

Chapter 9 Constructing Spline Curves

Figure 9.10 Bessel end condition spline.

K=l

K=0

Figure 9.11

Curvature plot of Bessel end condition spline.

Figure 9.10 shows the same data, but now using Bessel end conditions. Near
the endpoints, the curvature deviates from the ideal value, as shown in Figure
9.11.
Finally, Figure 9.12 shows the curve that is obtained using natural end conditions. The end curvatures are forced to be zero, causing considerable deviation
from the ideal value, as shown in Figure 9.13. This end condition should be
avoided unless the linear behavior near the ends is desired.

9.6 Finding a Knot Sequence

Figure 9.12

161

Natural end condition spline.

K=l

K=0

Figure 9.13 Curvature plot of natural end condition spline.

9.6 Finding a Knot Sequence
The spline interpolation problem is usually stated as "given data points p^
and parameter values w^, . . . . " Of course, this is the mathematician's way of
describing a problem. In practice, parameter values are rarely given and therefore
must be made up somehow. The easiest way to determine the ui is simply to set
Ui = i. This is called uniform or equidistant parametrization. This method is too

162

Chapter 9 Constructing Spline Curves
simplistic to cope with most practical situations. The reason for the overall poor^
performance of the uniform parametrization can be blamed on the fact that it
"ignores" the geometry of the data points.
The following is a heuristic explanation of this fact. We can interpret the
parameter u of the curve as time. As time passes from time UQ to time u^^
the point x(w) traces the curve from point x(wo) ^^ point x(wj^). With uniform
parametrization, x(w) spends the same amount of time between any two adjacent
data points, irrespective of their relative distances. A good analogy is a car driving
along the interpolating curve. We have to spend the same amount of time between
any two data points. If the distance between two data points is large, we must
move with a high speed. If the next two data points are close to each other, we
will overshoot since we cannot abruptly change our speed—we are moving with
continuous speed and acceleration, which are the physical counterparts of a C^
parametrization of a curve. It would clearly be more reasonable to adjust speed
to the distribution of the data points.
One way of achieving this is to have the knot spacing proportional to the
distances of the data points:
A. _ IIAp.ll
A,+i
||Ap,-+i||

(9^^^)

A knot sequence satisfying (9.15) is called chord length parametrization. Equation (9.15) does not uniquely define a knot sequence; rather, it defines a whole
family of parametrizations that are related to each other by affine parameter
transformations. In practice, the choices UQ = 0 and ui^ = lor:uQ^ = 0 and ui = L
are reasonable options.
Chord length usually produces better results than uniform knot spacing,
although not in all cases. It has been proven (Epstein [186]) that chord length
parametrization (in connection with natural end conditions) cannot produce
curves with corners^ at the data points, which gives it some theoretical advantage
over the uniform choice.
Another parametrization has been named "centripetal" by E. Lee [378]. It is
derived from the physical heuristics presented earlier. If we set

There are cases in which uniform parametrization fares better than other methods. An
interesting example is in Foley [239], p. 86.
A corner is a point on a curve where the tangent (not necessarily the tangent vector!)
changes in a discontinuous way. The special case of a change in 180 degrees is called a
cusp; it may occur even using the chord length parametrization.

9.6 Finding a Knot Sequence
lApil'
|Api+ill

A,+i

163

,1/2

(9.16)

the resulting motion of a point on the curve will "smooth out" variations in the
centripetal force acting on it.
Yet another parametrization v^as developed by G. Nielson and T. Foley [449].
It sets
A,- = di

^

3 Q/4-1
2 di_i + dj

3 0^+i4fl
2 dj -h di_^i

(9.17)

where dj = || Ap^|| and
0 ; = min

(.-e,.|).

and 0/ is the angle formed by p/-i, P/, p^+i- Thus 0^ is the "adjusted" exterior
angle formed by the vectors Ap/ and Ap^.j. As the exterior angle 0^ increases,
the interval A^ increases from the minimum of its chord length value up to a
maximum of four times its chord length value. This method was created to cope
with "wild" data sets.
We note one property that distinguishes the uniform parametrization from its
competitors: it is the only one that is invariant under affine transformations of the
data points. Chord length, centripetal, and the Foley methods all involve length
measurements, and lengths are not preserved under affine maps. One solution to
this dilemma is the introduction of a modified length measure, as described in
Nielson [446].^
For more literature on parametrizations, see Cohen and O'Dell [122], Hartley
and Judd [314], [315], McConalogue [421], and Foley [239].
Figures 9.14 to 9.21^^ show the performance of the discussed parametrization
methods for one sample data set. For each method, the interpolant is shown
together with its curvature plot. For all methods, Bessel end conditions were
chosen.
Although the figures are self-explanatory, some comments are in order. Note
the very uneven spacing of the data points at the marked area of the curves. Of
all methods, Foley's copes best with that situation (although we add that many
examples exist where the simpler centripetal method wins out). The uniform
9 The Foley parametrization was in fact first formulated in terms of that modified length
measure.
10 Kindly provided by T. Foley.

164

Chapter 9 Constructing Spline Curves

Figure 9.14

Chord length sphne.
K = 0.8

Figure 9.15 Curvature plot of chord length spUne.

Figure 9.16 Foley spline.

9.6 Finding a Knot Sequence

K = 1.6

K=0

Figure 9.17

Curvature plot of Foley spline.

Figure 9.18

Centripetal spline.

K = 2.0

K=0

Figure 9.19

Curvature plot of centripetal spline.

165

166

Chapter 9 Constructing Spline Curves

Figure 9.20 Uniform spline.

K = 70,000

• • • •

•

• •

•

•• • •

•

K = -80,000
Figure 9.21

Curvature plot of uniform spline.

Spline curve seems to have no problems there, if one just inspects the plot of the
curve itself. However, the curvature plot reveals a cusp in that region! The huge
curvature at the cusp causes a scaling in the curvature plot that annihilates all
other information. Also note hov^ the chord length parametrization yields the
"roundest" curve, having the smallest curvature values, but exhibiting the most
marked inflection points.
There is probably no "best" parametrization, since any method can be defeated by a suitably chosen data set. The foUov^ing is a (personal) recommendation. You may improve the shape of the curve, at an increase of computation
time, by the follov^ing hierarchy of methods: uniform, chord length, centripetal,
Foley. The best compromise betv^een cost and result is probably achieved by the
centripetal method.

9.7 The Minimum Property

167

9.7 The Minimum Property
In the early days of design, say, ship design in the 1800s, the problem had to be
handled of how to draw (manually) a smooth curve through a given set of points.
One way to obtain a solution was the following: place metal weights (called
ducks) at the data points, and then pass a thin, elastic wooden beam (called a
spline) between the ducks. The resulting curve is always very smooth and usually
aesthetically pleasing. The same principle is used today when an appropriate
design program is not available or for manual verification of a computer result;
see Figure 9.22.
The plastic or wooden beam assumes a position that minimizes its strain
energy. The mathematical model of the beam is a curve s, and its strain energy
£ is given by

h

(s)fds,

where /c denotes the curvature of the curve. The curvature of most curves involves
integrals and square roots and is cumbersome to handle; therefore, one often
approximates the preceding integral with a simpler one:

- /

s{u)

du.

(9.18)

Note that E is a vector; it is obtained by performing the integration on each
component of s.
Equation (9.18) is more directly motivated by the following example: when
an airplane is scheduled to fly from A to B, it will have to fly over a number of
intermediate "way points." The amount of fuel used by an airplane is mostly affected by its acceleration, which is essentially equivalent to the second derivative

j^
Figure 9.22

Spline interpolation: a plastic beam, the spline, is forced to pass through data points,
marked by metal weights, the ducks.

168

Chapter 9 Constructing Spline Curves
of its trajectory. Thus if the plane follows a cubic spline curve passing through
all the way points, it will be guaranteed to use the least amount offuell^i
In a more general setting, we may word this as: among all C^ curves interpolating the given data points at the given parameter values and satisfying the
same end conditions, the cubic spline yields the smallest value for each component of E. For a proof, let s{u) be the C^ cubic spline and let y{u) be another C^
interpolating curve. We can write y as
y{u) = s(u) + [y(u) - s(u)l
The preceding integrals are defined componentwise; we will show the minimum
property for one component only. Let s(u) and y(u) be the first component of s
and y, respectively. The "energy integral" E of y's first component becomes
pUi

E=

pUi

(sfdu

+ 2 /

JUQ

pUi

s(y - s)du + /

JUQ

(y -

sfdu.

JUQ

We may integrate the middle term by parts
I

Jua

s(y — s)du = s(y — s)

'"0

—/

Juo

s\y — s)du,

The first term vanishes because of the common end conditions. In the second
term, s is piecewise constant:
^ - 1

UT

"s{y - s)du = ^ Sj(y - s)
/
"0
/=o
*/
All terms in the sum vanish because both s and y interpolate. Since
/

{y - 'iydu > 0

JUQ

for continuous y ^ s,

f ^(yfdu > f ^(sfdu,
JUQ

JUQ

we have proved the claimed minimum property.
11 I am grateful to P. Crouch for bringing the airplane analogy to my attention.

(9.19)

9.8 C^ Piecewise Cubic Interpolation

169

The minimum property of splines has spurred substantial research activity.
The replacement of the actual strain energy measure £ by E is motivated by the
desire for mathematical simplicity. The curvature of a curve is given by
, ,

IIXAXII

llxlP
But we need ||x|| ^ 1 in order for ||x|| to be a good approximation to /c. This
means, how^ever, that the curve must be parametrized according to arc length;
see (10.7). This assumption is not very realistic for cubic splines in a design
environment; see Problems.
While the classical spline curve merely minimizes an approximation to (9.18),
methods have been developed that produce interpolants that minimize the true
energy (9.18), see [425], [99]. Moreton and Sequin have suggested to minimize
the functional f[K\t)f-dt instead; see [431].

9.8 C^ Piecewise Cubic Interpolation
Spline curves come in the B-spline form, but they may also be described as
piecew^ise Bezier curves. We now^ consider that approach, applied to piecewise
cubic interpolation. First, v^e try to solve the foUow^ing problem:
Given: Data points XQ, . . . , x^^ and tangent directions IQ, . . . , 1^ at those data
points; see Figure 9.23.
Find: A C^ piecewise cubic polynomial that passes through the given data
points and is tangent to the given tangent directions there.
It is important to note that we only have tangent directions^ that is, we have
no vectors with a prescribed length since our problem statement did not involve

Figure 9.23

C^ piecewise cubics: example data set.

170

Chapter 9 Constructing Spline Curves
a knot sequence. We can assume without loss of generahty that the tangent
directions 1^ have been normahzed to be of unit length;
111; II = 1 .
The easiest step in finding the desired piecewise cubic is the same as before: the
junction Bezier points hi^i are again given by b3/ = X/, / = 0 , . . . , L.
For each inner Bezier point, we have a one-parameter family of solutions: we
only have to ensure that each triple h^i-i-, h^^i^ ^^i^^i is coUinear on the tangent at
hi^i and ordered by increasing subscript in the direction of 1^. We can then find a
parametrization with respect to which the generated curve is C^ [see (5.30)].
In general, we must determine the inner Bezier points from
b3/+i = b3,+a,l„

(9.20)

b3.-i = b 3 , - ^ , _ i l „

(9.21)

so that the problem boils down to finding reasonable values for oti and )6^.
Although any nonnegative value for these numbers is a formally valid solution,
values for a^ and Pi that are too small cause the curve to have a corner at x^,
whereas values that are too large can create loops. There is probably no optimal
choice for of/ and Pi that holds up in every conceivable application—an optimal
choice must depend on the desired application.
A "quick and easy" solution that has performed decently many times (but also
failed sometimes) is simply to set
a, = ft = 0.4||Ax,||.

(9.22)

(The factor 0.4 is, of course, heuristic.)
The parametrization with respect to which this interpolant is C^ is the chord
length parametrization. It is characterized by

A, _ ft _
A,+i

a,+i

||Ax,||

^^^^^^

||Ax,+i||

A more sophisticated solution is the following: if we consider the planar curve
in Figure 9.24, we see that it can be interpreted as a function, where the parameter
t varies along the straight line through b3^ and b3^_^3. Then
l|Ab3,||-'l''^'>^-^^'l'
3 cos 0/
3 COS ^ij^i

9.8 C^ Piecewise Cubic Interpolation

Figure 9.24

171

Inner Bezier points: this planar curve can be interpreted as a function in an oblique
coordinate system with b3/, b3/+3 as the x-axis.

We are dealing with parametric curves, however, which are in general not planar
and for which the angles 0 and ^ could be close to 90 degrees, causing the
preceding expressions to be undefined. But for curves with 0^, ^^_^i smaller than,
say, 60 degrees, these expressions could be used to find reasonable values for of/
and ^i'.
1
3 cos 0 /

P^

lAX;

1
3 cos ^ij^\

I Ax, I

Since cos 60^ = 1/2, we can n o w make a case distinction:
iiAx,r
31/Ax,
|||Ax;||

if | 0 , | < 60^

(9.24)

Otherwise

and
II A x , I

^,= 1 3 1 ; ; ^
lAX;

if|*ml<60°

(9.25)

Otherwise.

This method has the advantage of having linear precision. It is C^ when the knot
sequence satisfies Aj/Aj^i = ^i/oti^\.
N o t e that neither of these t w o methods is affinely invariant: the first method—
(9.22)—does not preserve the ratios of the three points b3,_i, b3p b3/_^i since

172

Chapter 9 Constructing Spline Curves

Figure 9.25 FMILL tangents: the tangent at x^ is parallel to the chord through x^_i and Xj^i,
the ratios ||Ax^_i|| : ||Ax^|| are not generally invariant under affine maps.^^ The
second method uses angles, which are not preserved under affine transformations.
However, both methods are invariant under euclidean transformations.
We now address a different version of piecewise cubic interpolation:
Given: Data points XQ, . . . , xj^ together with corresponding parameter values
WQ,

. . . , Ui.

Find: A C^ piecewise cubic polynomial that passes through the given data
points.
One solution to this problem is provided by C^ (and hence also C^) cubic
splines, which are discussed in Section 9.4. Although those provide a global solution, a local one might be preferred for some applications where C^ smoothness
is not crucial. We are then faced with the problem of estimating tangents from
the given data points and parameter values.
The simplest method for tangent estimation is known under the name FMILL.
It constructs the tangent direction 1/ at x^ to be parallel to the chord through x^_i
and Xi^i:
Vi = Xi^i - Xi_i;

/• = 1 , . . . , L - 1.

(9.26)

Once the tangent direction V/ has been found,^^ the inner Bezier points are placed
on it according to Figure 9.25:

12 Recall that only the ratio of three collinear points is preserved under affine maps!
1 3 Note that here we do not have ||v^|| = 1!

9.8 C^ Piecewise Cubic Interpolation

173

b 3 . - i - b 3 , - ^ ^ / ^ - ; ^ V,,
3(A,_i + A,)

(9.27)

b 3 m = b3,+ .^^ \ ^ y r
3(A,_i + A,)

(9.28)

This interpolant is also known as a Catmull-Rom spline.
This construction of the inner Bezier points does not work at XQ and x^. The
next method, Bessel tangents, does not have that problem.
The idea behind Bessel tangents^^ is as follows: to find the tangent vector m^ at
x^, pass the interpolating parabola q/(w) through x^_i, x^, x^_^i with corresponding
parameter values w^_i, w^, u^j^i and let m^ be the derivative of q^. We differentiate
q, at Ui',
d
au

, ,

Written in terms of the given data, this gives
m, = ^^^^^x,_,
A,_i

+ ^AXi;
A,-

(9.29)

i=l,...,L-l,

where
«,=
The endpoints are treated in the same way: mo = d/duqi(uQ)^mi =
which gives
mo = 2—
^ Ax^_l
m^ = 2—

d/duqi_i(ui),

mi,
m^_i.

AL-1

Another interpolant that makes use of the parabolas q^ is known as an
Overhauser spline, after work by A. Overhauser [452] (see also [91] and [154]).
The i^^ segment s^ of such a spline (defined over [uj, W/+i]) is defined by
Si(u) = -^±i

q,(M) + —_^q,.+i(M);

A,14

They are also attributed to Ackland [2].

A,-

t=

l,...,L-2.

174

Chapter 9 Constructing Spline Curves
In other words, each S/ is a Hnear blend between q/ and q/+i. At the ends, one
sets So(w) = qo(w) and S£^_i(w) = q i - i W .
On closer inspection, it turns out that the last two interpolants are not different
at all: they both yield the same C^ piecewise cubic interpolant (see Problems).
A similar way of determining tangent vectors was developed by McConalogue
[421], [422].
Finally, we mention a method created by H. Akima [3]. It sets
m,- = ( 1 - Ci)2ii_i + Ci2ii,

where
AX;

A;
and
lAa;
|Aa,_2ll + ||Aa,-||
This interpolant appears fairly involved. It generates very good results, however,
in situations where curves are needed that oscillate only minimally.

9.9 Implementation
The following routines produce the cubic B-spline polygon of an interpolating
C^ cubic spline curve. First, we set up the tridiagonal linear system:
void set_up_system(knot,l,alpha,beta,gamma)
/*
given the knot sequence, the linear system for clamped end
condition B-spline interpolation is set up.
Input: knot:
knot sequence (all knots are simple; but,
in the terminology of Chapter 10, knot[0]
and knot[l] are of multiplicity three.)
points:
points to be interpolated
1:
number of intervals
Output:alpha, beta,gamma: 1-D arrays that constitute
the elements of the interpolation matrix.
Note: no data points needed so far!

V

9.9 Implementation

175

The next routine performs the LU decomposition of the matrix from the
previous routine. (Note that we do not generate a full matrix but rather three
linear arrays!)
void l_u_system(alpha,beta,gamma,],up,low)
/*
perform LU decomposition of tridiagonal system with
lower diagonal alpha, diagonal beta, upper diagonal gamma.
Input: alpha,beta,gamma:
1:
Output:low:
up:

the coefficient matrix entries
matrix size [ 0 , l ] x [ 0 , l ]
L-matrix entries
U-matrix entries

V
Finally, the routine that solves the system for the B-spline coefficients d^:
solve_system(up,1ow,gamma,1,rhs,d)
/* solve tridiagonal linear system
of size (1+1)(1+1) whose LU decomposition has entries up and low,
and whose right hand side is rhs, and whose original matrix
had gamma as its upper diagonal. Solution is d[0]
d[l+2].
Input: up,low,gamma: as above.
1:
size of system: 1+1 eqs in 1+1 unknowns,
rhs:
right hand side, i.e, data points with end
'tangent Bezier points' in rhs[l] and rhs[l+l].
Output:d:
solution vector.
Note shift in indexing from text! Both rhs and d are from 0 to 1+2.

V
In case Bessel ends are desired instead of clamped ends, this is the code:
void bessel_ends(data,knot,l)
/*
Computes B-spline points data[l] and data[l+l]
according to bessel end condition.
Input: data:

sequence of data coordinates data[0] to data[l+2].
Note that data[l] and data[l+l] are expected to
be empty, as they will be filled by this routine.
They correspond to the Bezier points bez[l] and bez[31-l].
knot: knot sequence
1:
number of intervals
Output: data: now including ''tangent Bezier points'' data[l], data[l+l].
*/

176

Chapter 9 Constructing Spline Curves
The centripetal parametrization is achieved by the following routine:
void parameters(data_x,data_y,l,knot)
/* Finds a centripetal parametrization for a given set
of 2D data points.
Input:
data_x, data_y: input points, numbered from 0 to 1+2.
1:
number of intervals.
Output:
Note:

knot:
knot sequence. Note: not (knot[l]=1.0)!
data_x[l], data_x[l+l] are not used! Same for data_y.

*/
A calling sequence that uses the preceding programs might look like this:
parameters(data_x,data_y,l,knot);
set_up_system(knot,l,alpha,beta,gamma);
1 _u__system (al pha, beta, gamma, 1, up, 1 ow);
bessel_ends(data_x,knot,l);
bessel_ends(data_y,knot,l);
sol ve_system(up,low,gamma,!,data_x,bspl_x);
solve_system(up,low,gamma,1,data_y,bspl_y);
Here, we solved the 2D interpolation problem with given data points in data_
X, data_y, a knot sequence knot, and the resulting B-spline polygon in bspl_x,
bspl__y. This calling sequence is realized in the routine c2_spl ine.c.

9.10 Problems
1 For the case of closed curves, C^ quadratic spline interpolation with uniform knots does not always have a solution. Why?^^
* 2 Show that interpolating splines reproduce cubic polynomial curves—that
they have cubic precision. This means that if all data points p/ are read off
from a cubic, p^ = c(w^), and the end tangent vectors are read off from the
cubic, then the interpolating spline equals the original cubic.
15 T. DeRose pointed this out to me.

9.10 Problems

177

3 Show that Akima's interpolant always passes a straight line segment
through three subsequent points if they happen to lie on a straight line.
* 4 Show that Overhauser splines are piecewise cubics with Bessel tangents at
the junction points.
* 5 One can generalize the quintic Hermite interpolants from Section 7.6 to
piecewise quintic Hermite interpolants. These curves need first and second
derivatives as input positions. Devise ways to generate second derivative
information from data points and parameter values.
PI Using piecewise cubic C^ interpolation, approximate the semicircle with
radius 1 to within a tolerance of 6 = 0.001. Use as few cubic segments as
possible. Literature: [172], [260].
P2 Program the following: instead of prescribing end conditions for interpolating C^ splines at both ends, prescribe first and second derivatives at TQ. The
interpolant can then be built segment by segment. Discuss the numerical
aspects of this method (they will not be wonderful).

This Page Intentionally Left Blank

W. Boehm
Differential Geometry I

Uifferential geometry is based largely on the pioneering work of L. Euler
(1707-1783), C. Monge (1746-1818), and C F. Gauss (1777-1855). One of
their concerns was the description of local curve and surface properties such
as curvature. These concepts are also of interest in modern computer-aided
geometric design. The main tool for the development of general results is the use
of local coordinate systems, in terms of which geometric properties are easily
described and studied. This introduction discusses local properties of curves
independent of a possible embedding into a surface.

10.1 Parametric Curves and Arc Length
A curve in E^ is given by the parametric representation

X = x(t) =

x(t)
y(t)
z(t)

t

e[a,b]c]

(10.1)

where its cartesian coordinates x, y, z are differentiable functions of t. (We have
encountered a variety of such curves already, among them Bezier and B-spline
curves.) To avoid potential problems concerning the parametrization of the curve,
we shall assume that
x{t)
x(t) =

m
m

^0,

te [a, bl

(10.2)
179

180

Chapter 10 W. Boehm: Differential Geometry I

Xi = X(ti)

-o—o—a—o-

a
Figure 10.1

Parametric curve in space.
where dots denote derivatives vs^ith respect to t. Such a parametrization is called
regular. An example is shown in Figure 10.1
A change x = x{t) of the parameter t^ where r is a differentiable function of
t^ will not change the shape of the curve. This reparametrization
will be regular
if f 7«^ 0 for all t e [a, fc], that is, we can find the inverse t = ^(r). Let

' = s(t)= f ||x\\dt
Ja

(10.3)

be such a parametrization. Because
. .
dxdr ,
dx .
xd^ = -——d^ = -—dr,
dr at
dr
s is independent of any regular reparametrization. It is an invariant parameter
and is called arc length parametrization of the curve. One also calls ds = \\x\\dt
the arc element of the curve.
Remark 7 Arc length may be introduced more intuitively as follows: let ti = a-\- iAt and
letA^ > 0 be an equidistant partition of the ^-axis. Let x^ = x{tj) be the corresponding sequence of points on the curve. Chord length is then defined by
(10.4)
where Ax^ = x^.^^ — x^. It is easy to check that for A^ -^ 0, chord length S
converges to arc length s, while Axj/At converges to the tangent vector x^
at Xj,

10.2 The Frenet Frame
Remark 2

181

Although arc length is an important concept, it is primarily used for theoretical
considerations and for the development of curve algorithms. If, for some application, computation of the arc length is unavoidable, it may be approximated by
the chord length (10.4).

10.2 The Frenet Frame
We w^ill nov^ introduce a special local coordinate system, linked to a point x(0 on
the curve, that w^ill significantly facilitate the description of local curve properties
at that point. Let us assume that all derivatives needed below do exist. The first
terms of the Taylor expansion of x(^ + A^) at t are given by
1
1
x(f + AO = X + xA^ + X-A^^ + X-A^^ +
2
6

^

Let us assume that the first three derivatives are linearly independent. Then
X, X, X form a local affine coordinate system w^ith origin x. In this system, x(t) is
represented by its canonical coordinates
A^H-...
lAf^ + . . .
v^here . . . denotes terms of degree four and higher in A^.
From this local affine coordinate system, one easily obtains a local cartesian
(orthonormal) system w^ith origin x and axes t, m, b by the Gram-Schmidt process
of orthonormalization, as shov^n in Figure 10.2:
t=^^,

llxll

m = bAt,

b=^^—^,

llxAxll

(10.5)

w^here A denotes the cross product.
The vector t is called tangent vector (see Remark 1), m is called main normal
vector^'^ and b is called binormal vector. The frame (or trihedron) t, m, b is called
the Frenet frame\ it varies its orientation as t traces out the curve.

1 We use the abbreviation A?-^ = (A?)^.
2

One often sees the notation n for this vector. We use m to avoid confusion with surface
normals, which are discussed later.

182

Chapter 10 W. Boehm: Differential Geometry I

Figure 10.2 Local affine system (left) and Frenet frame (right).
The plane spanned by the point x and the two vectors t, m is called the
osculating plane O. Its equation is
det

y

X

X

x"j _

1 1 0 0j~

det[y — X, x,x] = 0,

where y denotes any point on O. Its parametric form is
0(w, v) =x-\-ux-\-

vx.

Remark 5 The process of orthonormalization yields
m =

XX • X — XX • X
XX • X — XX • X

This equation may also be used for planar curves, where the binormal vector
b = t A m agrees with the normal vector of the plane.

10.5 Moving the Frame
Letting the Frenet frame vary with t provides a good idea of the curve's behavior in
space. It is a fundamental idea in differential geometry to express the local change
of the frame in terms of the frame itself. The resulting formulas are particularly
simple if one uses arc length parametrization. We denote differentiation with
respect to arc length by a prime. Since x' = t is a unit vector, we find the following
two identities:
x^x' = l

and

x^x''=:0.

The first identity states that the curve is traversed with unit speedy the second
one states that the tangent vector is perpendicular to the second derivative vector,
provided the curve is parametrized with respect to arc length.

10.3 Moving the Frame

183

Some simple calculations yield the so-called Frenet-Serret formulas:
t'
m'
b'

=
=
=

+/cm
+rb,

—Kt

(10.6)

—rm

where the terms K and r, called curvature and torsion^ may be defined both
in terms of arc length s and in terms of the actual parameter t. We give both
definitions:
K=K{S)=\W'\1
K=K(t)

=

IIXAXII

(10.7)

r = r(s) = ^ d e t [ x > ^ x n ,
T = r(t) =

det[x, X, x]
X A X |2 *

(10.8)

Figure 10.3 illustrates the formulas of (10.6).
Curvature and torsion have an intuitive geometric meaning: consider a point
x(s) on the curve and a "consecutive" point x(s + As). Let Aa denote the angle
between the two tangent vectors t and t(s + As) and let A^ denote the angle
between the two binormal vectors b and b(s + As), both angles measured in
radians. It is easy to verify that Aa = KAS -\-.. . and A^ = — r As + . . . , where
. . . denotes terms of higher degree in As. Thus, when As -> ds, we find that
K =

da

d7'

r =

ds'

— xm

—Kt

Figure 10.3 The geometric meaning of the Frenet-Serret formulas.

184

Chapter 10 W, Boehm: Differential Geometry I
In other words, K and — r are the angular velocities of t and b, respectively,
because the frame is moved according to the parameter 5.
Remark 4 Note that K and r are independent of the current parametrization of the curve.
They are euclidean invariants of the curve; that is, they are not changed by a rigid
body motion of the curve. Moreover, any two continuous functions K =K(S) > 0
and r = r(s) define uniquely (except for rigid body motions) a curve that has
curvature K and torsion r.
Remark 5 The curve may be written in canonical form in terms of the Frenet frame. Then
it has the form
As

x(s + As) =
|/crAs^ + . . .
where . . . again denotes terms of higher degree in As.

10.4 The Osculating Circle
The circle that has second-order contact with the curve at x is called the osculating
circle (Figure 10.4). Its center is c = x + pm, and its radius p = ^ is called the
radius of curvature. We shall provide a brief development of these facts. Using
the Frenet-Serret formulas of (10.6), the Taylor expansion of x(s -f As) can be
written as

1

2

x(s + As) = x(s) + 1 A s + -/cmAs + . . . ,

Figure 10.4 The osculating circle.

10.4 The Osculating Circle

m

185

p'^m

Figure 10.5 Construction of the osculating circle.
Let p* be the radius of the circle that is tangent to t at x and passes through the
point y = X + Ax, where Ax = tAs + ^/cmAs^ (see Figure 10.5). Note that y lies in
the osculating plane O. Inspection of the figure reveals that (^ Ax — p*m) Ax = 0;
that is, we obtain
. ^ 1 (Ax)^
2 mAx
From the definition of Ax, we obtain (Ax)^ = As^ + . . . and mAx = ^'^(As)^.
Thus p* = ^ + . . . . In particular, p = ^ as As ^- 0. Obviously, this circle lies in
the osculating plane.
Remark 6 Let x be a rational Bezier curve of degree n as defined in Chapter 13. Its curvature
and torsion at bo are given by
K

_n
=

~ 1 WQWI b

n

1

1'

n — 1 WQW^ C
n WiWi ah

(10.9)

where a is the distance between bo and b^, h is the distance of h^ to the tangent
spanned by bo and b^, and c is the distance of b3 from the osculating plane
spanned by bo, bj, and b2 (Figure 10.6). Note that these formulas can be used
to calculate curvature and torsion at arbitrary points x ( 0 of a Bezier curve after
subdividing it there (see Section 13.2).
Remark 7 An immediate application of (10.9) is the following: let x be a point on an integral
quadratic Bezier curve, that is, a parabola. Let 28 denote the length of a chord

186

Chapter 10 W. Boehm: Differential Geometry I

Figure 10.6 Frenet frame and geometric meaning of a, b, c.

Parabola

Figure 10.7 Curvature of a parabola.
parallel to the tangent at x, and let € be the distance between the chord and the
tangent. The radius of curvature at x is then P = j ^ (see Figure 10.7).
Remark 8 An equivalent way to formulate (10.9) is given by
K=l-

n — 1 WQWI area[bo, bj, b2]

n

w\

dist^[bo,bi]

(10.10)

and
_ 3 « — 2 WQW2, volume [bo, bj, hi^ h^]
2 n wiWi
area^[bo5bi,b2]

(10.11)

The advantage of this formulation is that it can be generalized to "higher-order
curvatures" of curves that span M^, 3 < d  0, while the ajj are arbitrary parameters.
It is interesting to note that curvature and torsion continuous curves exist that
are not K^ continuous^ (see Remark 4). Conversely,
x'^' = t'' = K^m + K(-Kt + rb)
3

Recall that K^ = d/c(s)/ds, v^here the prime denotes differentiation with respect to arc
length s of the (composite) curve. A formula for K' is provided by (23.2).

10.6 Composite Curves

189

implies that x.'^' is continuous if K' is and vice versa. To ensure x.'^ = x^', the
coefficients a and of^y must be the result of the application of the chain rule; that
is, with Qf2i = )6 and a^^x^y^ one finds that 0^32 = 30?^. Now, as before, the curve
is tangent continuous if
x_^ = ax_,

a > 0,

it is curvature and osculating plane continuous if in addition
..

9 ..

.

X_^ = Of X_ + j6x_,

but it is /c' continuous if in addition
x^ = a\_

+ 3a^x_ + yx_

and vice versa.
Remark / /

For planar curves, torsion continuity is a vacuous condition, but K' continuity is
meaningful.

Remark / 2

The preceding results may be used for the definition of higher-order geometric
continuity, A curve is said to be G^, or rth order geometrically continuous if a
regular reparametrization exists after which it is U. This definition is obviously
equivalent to the requirement of C^~^ continuity of K and C~^ continuity of r.
As a consequence, geometric continuity may be defined by using the chain rule,
as in the example for r = 3.

Remark 13 The geometric invariants curvature and torsion may be generalized for higherdimensional curves. Continuing the process mentioned in Remark 8, we find
that a J-dimensional curve has d — 1 geometric invariants. Continuity of these
invariants makes sense only in E^, as was demonstrated for d = 2'm Remark 11.
Remark / 4 Note that although curvature and torsion are euclidean invariants, curvature and
torsion continuity (as well as the generalizations discussed in Remarks 12 and 13)
are affinely invariant properties of a curve. Both are also projectively invariant
properties; see Boehm [76] and Goldman and Micchelli [267].
Remark 15 If two curve segments meet with a continuous tangent and have (possibly different) curvatures K_ and K_^ at the common point, then the ratio K_/K^ is also a
projectively invariant quantity. This is known as Memke's theorem; see Bol [88].

This Page Intentionally Left Blank

Geometric
Continuity

i4

uu

Ns

11.1 Motivation
Before we explain in detail the concept of geometric continuity, we will give
an example of a curve that is curvature continuous yet not twice differentiate.
Such curves (and, later, surfaces) are the objects that we will label geometrically
continuous.
Figure 11.1 shows three parabolas with junction points at the midpoints of
an equilateral triangle. According to (10.10), where we have to set all Wi equal
to 1, all three parabolas have the same curvature at the junction points. We thus
have a closed, curvature continuous curve. It is C^ over a uniform knot sequence.
But it is not C^ as is easily seen by sketching the second derivative vectors at the
junction points.
Differential geometry teaches us that our closed curve can be reparametrized
such that the new parameter is arc length. With that new parametrization, the
curve will actually be C^. Details are explained in Chapter 10. We shall adopt the
term G^ curves (second-order geometrically continuous) for curves that are twice
differentiable with respect to arc length but not necessarily twice differentiate
with respect to their current parametrization. Note that curves with a zero tangent
vector cannot be G^ under this definition. Planar G^ curves have continuously
varying signed curvature; G^ space curves have continuously varying binormal
vectors and continuously varying curvature.

191

192

Chapter 11 Geometric

Figure 11.1

Continuity

G^ continuity: a closed quadratic G^ spline curve.
The concept of geometric continuity is more appropriate when dealing with
shape; parametric continuity is appropriate when speed of traversal is an issue.^
Historically, several methods have been developed to deal with G^ continuity.
In the following pages, we present a unified treatment for most of these.

11.2 The Direct Formulation
Let b o , . . . , b3 and CQ, . . . , C3 be the control polygons of two cubic Bezier curves.-^
Since we are interested in G^ continuity here, we need only consider the control
points bi, b2, b3 = CQ, CJ, C2, all of which we assume to be coplanar. Referring to
Figure 11.2, let d be the intersection of the lines bib2 and c ^ .
We set
r_ = ratio(bi,b2,d),

(11.1)

r + = ratio(d, Ci, C2),

(11.2)

r = ratio(b2,b3,Ci).

(11.3)

Speed of traversal is important, for example, when the given curve is a vertical straight
line and we consider the motion of an elevator: higher orders of continuity of its path
ensure smoother rides.
The G^ conditions for general degrees will be identical, and so nothing is lost by concentrating on the cubic case.

11.2 The Direct Formulation

193

Figure 11.2 G^ continuity: using the direct formulation.
Letting A_,A^,B_,B^
denote the triangle areas in Figure 11.2, we can invoke
(10.10) in order to express the curvatures K_ and /c+ of the left and right segments
atbj:
4
A.
3||b3-b2lP'

K,=

+

4
A,
3||ci-coll3"

If these two curvatures agree, we have that
A.
A;

(11.4)

= r\

Referring to the figure again, we see that
A_

T_='-

B4+5
A+ = r.-'

B^ = r.

Inserting this into (11.4) yields our desired G^ condition:
'-'+'

(11.5)

With the notation of Figure 11.2, and setting /_ = ||b3 — b2l|, /+ = ||ci — b3||,
we have

Thus
r_r.

194

Chapter 11 Geometric Continuity
This is known as Memke's theorem and states that the ratio of left and right
curvatures on any point of a curve is invariant under affine maps. This follows
since we only use affinely invariant ratios in the formulation of K_/K_^.

11.3 The y, Vj and fi Formulations
Using the setting of Section 11.2, we observe that our composite curve could be
made C^ if we introduced a knot sequence with interval lengths A_, A_^ satisfying
A_/A_^ = r. Using (11.5), we define

We then have
ratio(bi, hi, d) =

yA+

and

ratio(d, Cj, C2) =

In the case that y = 1, we have the special case of a C^ piecewise cubic curve. See
also Figure 11.3.
W. Boehm used this framework for his development of G^ cubic splines; see
[71].
If a C^ curve x is curvature continuous at a point x(w), then we must have

Thus the three vectors involved must satisfy a relationship
(11.6)
The constant v depends on the parametrization of the curve; if we change u to
ku, we will have to replace v by ^y. The v formulation of G'^ continuity is due
to G. Nielson [442].

Figure 11.3 G^ continuity: using the y formulation.

11.4 G^ Cubic Splines

195

There is a one-to-one relationship between the constants y and v.
„= 2 ( ^ ^ ) 1 ^ .

,11.7)

first found by W. Boehm [71].
A similar approach was taken by B. Barsky [40], [47];he uses
Pi = —^
^+

and

^2 = 1^

(11-8)

as the descriptors of G^ continuity and calls them bias and tension, respectively.
Why three or four different formulations for G^ continuity of piecewise cubic
curves? The reason is partly historical, and partly depends on applications. In
fact, the preceding formulations are by no means the only ones—the discussion
of G^ continuity goes back as far as Baer [12], Bezier [59], Geise [256], and
Manning [414].
Applications that aim at constructing surfaces will be better served by ^, y, or y
splines. These involve a knot sequence and thus lend themselves to the framework
of tensor product surfaces; see Chapter 16.
Freeform curve design, on the other hand, will benefit more from the direct
formulation since it is linked the closest to the curve geometry. The direct
approach is the most geometric, followed by the y formulation, which needs
a knot sequence. The least geometric are the v and fi formulations; their defining
quantities are not invariant under scaling of the knot sequence.

11.4 C^ Cubic Splines
We start with a control polygon d o , . . . , d^^. In the context of C^ cubic B-splines,
we now needed a knot sequence in order to place the inner Bezier points on the
control polygon legs; the junction points then were fixed by the C^ conditions.
In our case, we have more freedom: we may place the inner Bezier points
anywhere on the control polygon legs; the junction points are then fixed by the
conditions.
To be more precise, consider Figure 11.4. Placing h^j-i ^^ ^^e polygon leg
dp dj_^i amounts to picking a number aj (between 0 and 1) and then setting
b3,_2 = ( l - a M + M . + l .
Similarly, we place h^i_i by picking a number coj and setting

(11.9)

196

Chapter 11 Geometric Continuity

Figure 11.4 G^ conditions; inner Bezier points may be placed on the control polygon legs. The junction
points then may be found using the G^ condition.
b3/-i = (1 - C0i)di + coidi^^,

(11.10)

In the same manner, by choosing numbers aj^i and o^^+i, we determine h^i-^^ and
We still have to determine the junction point b3/. Upon comparing Figures
11.4 and 11.2, we see that we need the quantities A/ = ratio(b3/_2, b3/_i, d^_^i)
and Pi = ratio(d,+i, ^3i-\-h ba^+i)- Since
b3,-l = i ^ b 3 , _ 2 + ^ d , . i
1 — aj

1 — aj

(11.11)

and
b3/+l =

^/+1 +

b3/+2?

(11.12)

we have
A, = ^
Setting Yi — ^kiPilil

,

p, = - ^ i ± ^ .

(11.13)

+ y/T~pi)^ we find the desired junction point to be
b3, = ( l - r , ) b 3 , _ i + r,b3,+i.

(11.14)

Continuing in this manner for all /, we have completed the definition of a G^
spline curve. We note that it is advisable to restrict all a^ and coi to be between 0
and 1. It is possible, however, to violate that condition: we only have to ensure
that Xi and pi have the same sign. As long as they do, b3^ is computable from
(11.14).
For an open polygon, we set Q^Q = 0 and (Oi_2 = 1. This ensures the usual
bi = di and h^L-4 = ^L-i-

11.4 G^ Cubic Splines

197

Figure 11.5 G^ splines: these two cubics are a G^ spline but do not possess a G^ control polygon.

Our development of G^ splines is solely based upon ratios, hence G^ spline
curves will be mapped to G^ spline curves by affine maps. We may also say that
continuity is affinely invariant.
There is one interesting difference between the construction for a G spline
and the corresponding construction for a C^ spline: every C^ piecewise cubic
possesses a B-spline control polygon—but not every G^ piecewise cubic curve
possesses a G^ control polygon. The two cubics in Figure 11.5 are curvature
continuous, yet they cannot be obtained with the preceding construction: the
control point d would have to be at infinity.
In interactive design, one would use G^ cubic splines in a two-step procedure.
The design of the G^ control polygon may be viewed as a rough sketch. The
program would estimate the inner Bezier points automatically, and the designer
could fine-tune the curve shape by readjusting them where necessary. For this
fine-tuning, it is important to observe that h^j_i,h^i^i is tangent to the curve.
Instead of prescribing numbers aj and coj—not very intuitive!—a designer may
thus specify tangents to the curve, and the a^, coj can be computed. Figure 11.6
gives examples.
We have just described G^ splines using the direct G^ formulation. Using the
y formulation, we arrive at y-splines, which use a set of y^ and a knot sequence,
employing the principles of Section 11.3. We then have

198

C h a p t e r 11 Geometric

Figure 11.6

Continuity

G^ splines: a shape may be varied by prescribing tangents in addition to the control
polygon.

A/_2

Figure 11.7

A,_i

A,-

A,„i

y-splines: the Bezier points are connected to the G^ control polygon by the shown ratios.

a,_l:

A,-i +
VAi-l

vAi

+ A,_i +

Yi^i^i

(11.15)

ind

^/-i = vAi-l

vAi
+ A,_i +

Yi^i^i

(11.16)

The geometry of a y-spline curve is shov^n in Figure 11.7. N o t e that for all
Yl -^ 0, the curve v^ill tend tov^ard its control polygon.

11.5 Interpolating G^ Cubic Splines

199

11.5 Interpolating C^ Cubic Splines
We may also use G^ cubics to interpolate to given data points x^; / = 0,. . . , L.
In the C^ case, we had to supply a knot sequence in addition to the data points.
Now, we have to specify a sequence of pairs a^, coi. How to do this effectively is
still an unsolved problem, so let us assume for now that a reasonable sequence
of Qfp coi is given. Setting b3/ = x^, (11.14) yields

d,+i + y/^iOti+A+2, /• = 1 , . . . , L - 1.

(11.17)

Together with two end conditions, we then have L — 1 equations for the L + 1
unknowns d^. A suitable end condition is to make d^ a linear combination of
the first three data points: d^ = WXQ + v^\ + wi^i. In our experience, (w, v^ w) =
(5/6,1/2, —1/3) has worked well. A similar equation then holds for dj^^^. For
the limiting case of aj -> 0 and coi -^ 1, the interpolating curve will approach the
polygon formed by the data points. In terms of the y formulation, this spline
type was investigated in [209].
Nielson [442] derived the G^ interpolating spline from the v formulation.
Assuming that the data points x^ have parameter values TJ assigned to them, and
using the piecewise cubic Hermite form, the interpolant becomes
xiu) = x^H^ir) + m,A,Hl{r) + A,m,+iH|(r) + x,_,^Hl(r),

(11.18)

where the Hj* are cubic Hermite polynomials from (7.14) and r = (u — r^)/A^ is
the local parameter of the interval (r/, r/+i). In (11.18), the x^ are the known data
points, whereas the m^ are as yet unknown tangent vectors. The interpolant is
supposed to be G^; it is therefore characterized by (11.6), more specifically,
x+(T,)-x_(T,) = y,m,

(11.19)

for some constants v^, where m^ = x(r^). The Vj are constants that can be used to
manipulate the shape of the interpolant; they will be discussed soon. We insert
(11.18) into (11.19) and obtain the linear system

(

A;AX; 1

^ - ^ i +
A,_i

A;_iAX;\

'I

A,

/

^

^

1

= A,m,_i + (2A,_i + 2 A, + - A,_iA,y,)m,
2
+ A,_im,+i;/=l,...,L-l
(11.20)

200

Chapter 11 Geometric Continuity

Figure 11.8 Infinite tension: an interpolating y-spline curve with the cross plot of the y-component.

Together with two end conditions, (11.20) can be used to compute the unknown tangent vectors m^. The simplest end condition is prescribing mo and m^,
but any other end condition from Chapter 9 may be used as well. Note that this
formulation of the y-spline interpolation problem depends on the scale of the r^;
it is not invariant under affine parameter transformations.
If the Vi are chosen to be nonnegative, the linear system (11.20) is solvable;
in the special case of all V/ = 0, it results in the standard C^ cubic spline. For the
case of all V/ -^ oo, the interpolant approaches the polygon formed by the data
points. Although the resulting curve will look piecewise linear, it is actually C^,
as shown in the cross plot (see Section 6,6) in Figure 11.8. Only the y-component
is shown; it has zero slopes at the r^.

11.6 Higher-Order Geometric Continuity
Just as we can define higher-order parametric continuity C^, we may also define
higher-order geometric continuity. We say that a curve is rth order geometrically
continuous, or G^, at a given point, if it can be reparametrized such that it will
become C^(see Remark 12 in Section 10.6). In particular, the new parameter
might be arc length.
To derive conditions for G^ continuity, we start with a composite O curve
x(w) with a global parameter u. At a given parameter value w, derivatives from
the left and from the right agree:

—-x_ = —-x.
au^
aw +'

:0,...,r.

(11.21)

11.6 Higher-Order Geometric Continuity

201

Figure 11.9 G^ continuity: a segment of a C^ curve may be reparametrized. The resuhing curve is not
C^ anymore, but still G^.
Now let us reparametrize the right segment by introducing a new parameter
t = t{u); see Figure 11.9. By our earlier definition, the resulting composite curve
will be G^, while it is clearly not C any more. We will now study the conditions
for G^ continuity using this composite G^ curve.
Modifying (11.21) so as to incorporate the new parametrization yields:
d'
dw

d'
du'

(11.22)

The terms on the right-hand side of this equation may be expanded using the
chain rule. For / = 1, we obtain
X

d^
dw'

= X

(11.23)

where a prime denotes differentiation with respect to w, and a dot denotes
differentiation with respect to t. For / = 2, we have to apply both the chain and
the product rule to the right-hand side of (11.23):

„

.. /dty

. d^t

(11.24)

For the case / = 3:
X ^'/-x

(^\^

\duj

dtdh . X. dft
+ dt/3dudu^

3x — —

(11.25)

Let us define a/ = dH/duK Then these equations may be written in matrix form:
x'_

x_"

x_'"

•

_.

0
_«3

0 •

0

«2

3o;ia2

4\

(11.26)
_x+_

202

Chapter 11 Geometric Continuity
The lower triangular matrix in (11.26) is called a connection matrix^ it
connects the derivatives of one segment to that of the other. For rth order
geometric continuity, the connection matrix is a lov^er triangular r y. r matrix;
for more details, see Gregory [291] or Goodman [270]. See also the related
discussion in Section 10.6. The connection matrix is a powerful theoretical tool,
and has been used to derive variation diminishing properties of geometrically
continuous curves (Dyn and Micchelli [181]), to show the projective in variance
of torsion continuity (Boehm [76]), and for other theoretical pursuits (Goldman
andMiccheni[267]).
The definition of geometric continuity has been used by Manning [414],
Barsky [40], Barsky and DeRose [43], Degen [152], Pottmann [488], [489], and
Farin [192]. In terms of classical differential geometry, the concept of G^ is called
order two of contact; see do Carmo [170]. It was used in a constructive context
by G. Geise [256] as early as 1962.
An interesting phenomenon arises if we consider geometric continuity of order
higher than two. Consider a G^ space curve. It is easy to verify that it possesses
continuous curvature and torsion. But the converse is not true: there are space
curves with continuous curvature and torsion that are not G^ (Farin [192]). This
more general class of curves, called Frenet frame continuous^ has been studied by
Boehm [74]; see also Section 10.6 and Hagen [298], [299], They are characterized
by a more general connection matrix than that for G^ continuity; it is given by
«!

0

0•

0^2

oi^

0

L «3

^

«i .

where fi is an arbitrary constant. For higher-order Frenet frame continuity, we
have to resort to higher-dimensional spaces; this has been carried out by Dyn and
MiccheUi [181], Goodman [270], Goldman and MiccheUi [267], and Pottmann
[487]; see also the survey by Gregory [291]. An even more general concept than
that of Frenet frame continuity has been discussed recently by H. Pottmann [488].
A condition for torsion continuity of two adjacent Bezier curves with polygons
b o , . . . , b^ and CQ, . . . , c„ is given by
volume[b„_3,... ,b^] _ volume[co,... ,03]
|Ab^_ll|6
IIAq•Oil116

( H 27)

See Boehm [73], Farin [192], or Hagen [298].
A nice geometric interpretation of the fact that torsion continuity is more
general than G^ continuity is due to W. Boehm [73]. If b „ _ 3 , . . . , b„ and CQ, . . . , C3
are given such that the two curves are G^, can we vary C3 and still maintain G^
continuity? The answer is yes, and C3 may be displaced by any vector parallel

11.8 Problems

203

to the tangent spanned by b„_i and c^. But we may displace C3 by any vector
parallel to the osculating plane spanned by b„_2, b^, Ci and still maintain torsion
continuity!

11.7 Implementation
We include a direct G^ spline program. It assumes that the piecewise Bezier
polygon has been determined except for the junction points h^i^ which will be
computed:
void direct_gspline(l,bez_x,bez_y)
/ * From given interior Bezier points,
the junction Bezier points b3i are found from the G2 conditions.
Input: 1:
no of cubic pieces.
bez_x,bez_y: interior Bezier points b_{3i+l}, b_{3i-l}.
Output:bez_x,bez_y: completed piecewise Bezier polygon.
Note:
b_0 and b_{31+3} should be provided, too!

V

11.8 Problems
1 Figure 11.1 shows a triangle and an inscribed piecewise quadratic curve.
Find the ratio of the areas enclosed by the curve and the triangle.
2 Show that the average of two G^ piecewise cubics is in general not G^.
3 Find an example of a G^ torsion continuous curve that is not
* 4 Let a G^ curve consist of two cubic Bezier curves. The derivatives of the
two curves at the junction point are related by a connection matrix. Work
out the corresponding connection matrix for the Bezier points.
* 5 Show that a nonplanar cubic cannot have zero curvature or torsion anywhere.
* 6 The G^ piecewise cubic from Figure 11.5 cannot be represented as a direct
G^ spline. Can it be obtained from a v-spline interpolation problem.^
PI Change the programs for interpolating C^ cubics so that they compute
interpolating G^ splines.

This Page Intentionally Left Blank

Conic Sections

V^onic sections (or, simply, conies) have received the most attention throughout
the centuries of any knov^n curve type. Today, they are an important design tool in
the aircraft industry; they are also used in areas such as font design. A great many
algorithms for the use of conies in design were developed in the 1940s; Liming
[390] and [391] are two books with detailed descriptions of those methods. A
thorough development of conies can also be found in [85] and [202].
The first person to consider conies in a CAD environment was S. Coons [124].
Later, Forrest [240] further investigated conies and also rational cubics. We shall
treat conies in the rational Bezier form; a good reference for this approach is Lee
[377]. We present conies partly as a subject in its own right, but also as a first
instance of rational Bezier and B-spline curves (NURBS), to be discussed later.

12.1 Projective Maps of tiie Real Line
Polynomial curves, as studied before, bear a close relationship to affine geometry.
Consequently, the de Casteljau algorithm makes use of ratios, which are the
fundamental invariant of affine maps. Thus the class of polynomial curves is
invariant under affine transformations: an affine map maps a polynomial curve
onto another polynomial curve.
Conic sections, and later rational polynomials, are invariant under a more
general type of map: the so-called projective maps. These maps are studied in
projective geometry. This is not the place to outline the ideas of that kind of
geometry; the interested reader is referred to the text by Penna and Patterson
[461] or to [85] and [202]. All we need here is the concept of a projective map.

205

206

Chapter 12 Conic Sections

Figure 12.1 Projections: a straight Hne L is mapped onto another straight Hne \J by a projection. Note
how ratios of corresponding triples of points are distorted.
We start with a map that is famihar to everybody with a background in
computer graphics: the projection. Consider a plane (called image plane) P and
a point o (called center or origin of projection) in E-^. A point p is projected onto
P through o by finding the intersection p between the straight line through o and
p with P. For a projection to be well defined, it is necessary that o is not in P.
Any object in E-^ can be projected into P in this manner.
In particular, we can project a straight line, L, say, onto P, as shown in
Figure 12.1. We clearly see that our projection is not an affine map: the ratios
of corresponding points on L and V are not the same. But a projection leaves
another geometric property unchanged: the cross ratio of four collinear points.
The cross ratio, cr, of four collinear points is defined as a ratio of ratios [ratios
are defined by (3.6)]:
cr(a, b, c, d) =

ratio(a,b, d)
ratio(a, c, d)

(12.1)

This particular definition is only one of several equivalent ones; any permutation
of the four points gives rise to a vahd definition. Our convention (12.1) has
the advantage of being symmetric: cr(a, b, c, d) = cr(d, c, b, a). Cross ratios were
first studied by C. Brianchon and F. Moebius, who proved their invariance under
projective maps in 1827; see [429].
Let us now prove this invariance claim. We have to show, with the notation
from Figure 12.2, that
cr(a, b, c, d) = cr(a, b, c, d).
This fact is called the cross ratio theorem.

(12.2)

12.1 Projective Maps of the Real Line

207

Figure 12.2 Cross ratios: the cross ratios of a, b, c, d and a, b, c, d only depend on the angles shown
and are thus equal.
For a proof, consider Figure 12.2.
Denote the area of a triangle with vertices p, q, r by A(p, q, r). We note that,
for instance,
ratio(a, b, c) = A(a, b, o)/A(b, c, o).
This gives
cr(a,b, c, d) =

A(a, b, o)/A(b, d, o)
A(a,c,o)/A(c,d,o)

_ lili sin a/l2l4 sin(j6 + y)
lil^ sin((y + P)/l^U sin y
_ sin Of/ sin(^ + y)
sin(a + y6)/ sin y
Thus the cross ratio of the four points a, b, c, d only depends on the angles at
o. The four rays emanating from o may therefore be intersected by any straight
line; the four points of intersection will have the same cross ratio, regardless of

208

Chapter 12 Conic Sections
the choice of the straight Hne. All such straight lines are related by projections,
and we can therefore say that projections leave the cross ratio of four collinear
points invariant. Since the cross ratio is the same for any straight line intersecting
the given four straight lines, one also calls it the cross ratio of the four given lines.
A concept that is slightly more abstract than that of projections is that of
projective maps. Going back to Figure 12.1, we can interpret both L and \J as
copies of the real line. Then the projection of L onto \J can be viewed as a map of
the real line onto itself. With this interpretation, a projection defines a projective
map of the real line onto itself. On the real line, a point is given by a real number,
so we can assume a correspondence between the point a and a real number a.
An important observation about projective maps of the real line to itself is that
they are defined by three preimage and three image points. To observe this, we
inspect Figure 12.2. The claim is that a, b, d and their images a, b, d determine
a projective map. It is true since if we pick an arbitrary fourth point c on L, its
image c on V is determined by the cross ratio theorem.
A projective map of the real line onto itself is thus determined by three
preimage numbers a^fo,c and three image numbers 5,fc,2. The projective image
? of a point t can then be computed from
cr(^,fo,^, c) = cr(5, &, ?, c).
Setting p = (b — a)/(c — b) and p = (b — a)/{c — b)^ this is equivalent to

(t-a)/(c-t)

(i-a)/(c-iy

Solving for i:
^^{t-

a)pc + (c- t)ap
pic-t) + p(t-a) *

A convenient choice for the image and preimage points is a = a = 0^ c = c = 1.
Equation (12.3) then takes on the simpler form
^^-

'^
pil-t)-\-pt

(12.4)

Thus a projective map of the real line onto itself corresponds to a rational
linear transformation. It is left for the reader to verify that the projective map
becomes an affine map in the special case that p — p.

12.2 Conies as Rational Quadratics

209

12.2 Conies as Rational Quadratics
We will use the following definition for conic sections: a conic section in E^ is
the projection of a parabola in ¥? into a plane. We take this plane to be the plane
z=l. Figure 12.3 gives an example of how to obtain a conic as the projection of
a 3D parabola. Since we will study planar curves in this section, we may think
of this plane as a copy of e, thus identifying points [x y]^ with [x y 1 ]^.
Our special projection is characterized by
X

y ->
_z _

x/z
y/z
1

Note that a point [x y ]^ is the projection of a whole family of points: every
point on the straight line [ wx wy w ]^ projects to [ x y ]^. In the following,

Figure 12.3 Conic sections: a parabolic arc in 3D space is projected into the plane z = V^ the result, in
this example, is part of a hyperbola.

210

Chapter 12 Conic Sections

Figure 12.A Projections: the special projection that is used to write objects in the plane ;2 = 1 as
projections of objects in E^.
we will use the shorthand notation [ wx w ]^ with x E E^ for [ wx wy w ]^}
An illustration of this special projection is given in Figure 12.4.
Let c{t) G E^ be a point on a conic. Then there exist real numbers WQ^ Wi, wi
and points bo, b^, hi G E^ such that
c{t) =

wohoBlit) + wihiBlit)

+

woBlit) + wiB\it) +

wihiBlit)
wiBlit)

(12.5)

Let us prove (12.5). We may identify c(t) e E^ with [ c(0 1 F ^ E^. This point
is the projection of a point [ w(t)c(t)
w(t) ]^, which lies on a 3D parabola. The
third component w(t) of this 3D point must be a quadratic function in t, and
1 The set of all points [ wx wy w ]^ is called the homogeneous form or homogeneous
coordinates oi{x y ]^.

12.2 Conies as Rational Quadratics

211

may be expressed in Bernstein form:
wit) = woBlit) + w^Blit) + W2B\{t),
Having determined w{t)^ we may now write
wit)

\c{t)-\\^it)Y.w,B]{t)l

L 1 J L E^.^fW J

Since the left-hand side of this equation denotes a parabola, we may write

/=0

with some points p^ e E^. Thus
2

2

J2PtBJit) = c{t)J2wiBJit),
i=0

(12.6)

/=0

and hence
c{t) =

PoBl(t)^p^B](t)+p2Bl(t)
woBl(t) + wiB^it) + W2Bl(t)'

Setting p^ = Wihi now proves (12.5).
We call the points b/ the control polygon of the conic c; the numbers Wi are
called weights of the corresponding control polygon vertices. Thus the conic control polygon is the projection of the control polygon with vertices [ w^^ w^ ]^,
which is the control polygon of the 3D parabola that we projected onto c.
The form (12.5) is called the rational quadratic form of a conic section. If
all weights are equal, we recover nonrational quadratics, that is, parabolas. The
influence of the weights on the shape of the conic is illustrated in Figure 12.5. In
that figure, we have chosen
bo =

0
1
0
,bi =
1
0 ,b2 = 0

Note that a common nonzero factor in the w^ does not affect the conic at all.
If WQ 7^ 0, we may therefore always achieve WQ — lhydi simple scaling of all Wi,

212

Chapter 12 Conic Sections

Figure 12.5 Conic sections: in the two examples shown, WQ = W2 = 1. As wi becomes larger, that is,
as [t^ibi, Wi] moves "up" on the z-^xis, the conic is "pulled" toward b^.

12.2 Conies as Rational Quadratics

213

There are other changes of the weights that leave the curve shape unchanged:
these correspond to rational linear parameter transformations. Let us set

p{l-t)-\-t

„ (1-0-

''"'-'^

p{l-t)

+t

[corresponding to the choice p = 1 in (12.4)]. We may insert this into (12.5) and
obtain:
c(?) ^ ^^^0^0-^0^^^) + P^ibiB^(f) + W2h2B\{t)
p^woBlii) + pwiB\{i) + W2B\{i)
Thus the curve shape is not changed if each weight Wi is replaced by Wi — p^~^Wi
(for an early reference, see Forrest [240]). If, for a given set of weights iv^^ we
select

we obtain WQ = wi-, and, after dividing all three weights through by wi-, we have
t2/o = t2^2 = 1- A conic that satisfies this condition is said to be in standard form.
All conies with U/Q, wii^O may be rewritten in standard form with the choice of
yo, provided, of course, that WI/WQ > 0.
If in standard form, that is, WQ = W2 = 1, the point s = c ( | ) is called the
shoulder point. The shoulder point tangent is parallel to bob2. If we set m =
(bo + b2)/2, then the ratio of the three collinear points m, s, b^ is given by
ratio(m, s, b^) = Wi.

(12.8)

We finish this section with a theorem that will be useful in the later development of rational curves: Any four tangents to a conic intersect each other in the
same cross ratio. The theorem is illustrated in Figure 12.6. The proof of this four
tangent theorem is simple: one shows that it is true for parabolas (see Problems).
It then follows for all conies by their definition as a projection of a parabola and
by the fact that cross ratios are invariant under projections. This theorem is due
to J. Steiner. It is a projective version of the three tangent theorem from Section
4.1.

214

Chapter 12 Conic Sections

Figure 12.6 The four tangent theorem: four points are marked on each of the four tangents to the
shown conic. The four cross ratios generated by them are all equal.

12.5 A de Casteljau Algorithm
We may evaluate (12.5) by evaluating the numerator and the denominator
separately and then dividing through. A more geometric algorithm is obtained
by projecting each intermediate de Casteljau point [ w^W- w^- ] into E^:

W':

W,

(12.9)

w^here

w\(t) = (\-t)w'r\t)

,.r-l/
+ tw'^[{t)

(12.10)

This algorithm has a strong connection to the four tangent theorem: if wt
introduce weight points

OL'iit) •

(12.11)

12.4 Derivatives

215

then
cr(b;:,q;,b^+\b^^i) = i ^

(12.12)

assumes the same value for all r, /. Though computationally more involved than
the straightforward algebraic approach, this generalized de Casteljau algorithm
has the advantage of being numerically stable: it uses only convex combinations,
provided the weights are positive and t G [0,1].

12.4 Derivatives
To find the derivative of a conic section, that is, the vector c{t) = dc/d^, we may
employ the quotient rule. For a simpler derivation, let us rewrite (12.6) as
p(t) = w(t)c(t).
We apply the product rule:
p(^) = w{t)c(t) + w(t)c(t)
and solve for c(t):
c(t) = -^\p(t)
w(t)

- w{t)cm

(12.13)

We may evaluate (12.13) at the endpoint ^ = 0:
2
c(0) = —[wihi - woho - (wi - WQ)ho].
After some simplifications, we obtain
c(0) = ^ ^ A b o .

(12.14)

c(l) = ^ A b i .

(12.15)

WQ

Similarly, we obtain

W2

Let us now consider two conies, one defined over the interval [UQ, U^ with
control polygon bo, b^, hi and weights W/Q, ^ i , wi and the other defined over

216

Chapter 12 Conic Sections
the interval [wj, ui] with control polygon b2, b3, b4 and weights W2, w^, W4, Both
segments form a C^ curve if
^ A b , = ^Ab2,
Ao
Ai

(12.16)

where the appearance of the interval lengths A^ is due to the application of the
chain rule, which is necessary since we now consider a composite curve with a
global parameter u.

12.5 The Implicit Form
Every conic c(^) has an implicit representation of the form
f{x,y) = 0,
where /^ is a quadratic polynomial in x and y. To find this representation, recall
that c{t) may be written in terms of barycentric coordinates of the polygon
vertices bo, bj, hi:
c(t) = robo + ribi + tihi;

(12.17)

see Section 3.5. Since c(t) may also be written as a rational Bezier curve (12.5),
and since both representations are unique, we may compare the coefficients of
the hi'.
To = [woa-tf]/D,

(12.18)

r^ = [2wit(l-t)]/D,

(12.19)

T2 = [W2t^]/D,

(12.20)

where D = ^ WjBJ, We may solve (12.18) and (12.20) for (1 - t) and if, respectively. Inserting both expressions into (12.19) yields
Tj — t
WQW2

This may be written more symmetrically as

ror2

W0W2

(12.21)

12.5 The Implicit Form

217

This is the desired impHcit form, since the barycentric coordinates TQ, r^, X2 of
c(^) are given by

^0 =

1 c^ b\ ^2
\cy bl bl
1 1 1
1 b^ fcJ ^2

1 1 1

bl

^1=-

bl

b^ cy bl
1 1 1

X2=-

LX

1

1

b\

(f\

bl b\ cy\
1 1 1
^Q b\ b^ 1

K b\ fed

1

1 1 1

The implicit form has an important application: suppose we are given a conic
section c and an arbitrary point x G E^. Does x lie on c? This question is hard to
ansv^er if c is given in the parametric form (12.5). Using the implicit form, this
question is ansv^ered easily. First, compute the barycentric coordinates TQ, r^, Xi
of X v^ith respect to bg, b^, b2. Then insert TQ, r^, Xi into (12.21). If (12.21) is
satisfied, x lies on the conic (but see Problems).
The implicit form is also important when dealing with the IGES data specification. In that data format, a conic is given by its implicit form f(x^ y) = 0 and
two points on it, implying a start point and endpoint bo and hi of a conic arc.
Many applications, however, need the rational quadratic form. To convert to
this form, we have to determine b^ and its weight w/^, assuming standard form.
First, we find tangents at bo and hi', we know that the gradient of /^ is a vector that is perpendicular to the conic. The gradient at bo is given by /"'s partials:
V/*(bo) = [/x(bo)9/y(bo)F- The tangent is perpendicular to the gradient and thus
has direction V-^/"(bo) = {—fy(ho)',fx0^o)V' Thus our tangents are given by
t o ( 0 = b o + ^V-^/"(bo)

and

t2(s)=b2 + sVV(b2).
Their intersection determines b^. Next, we compute the midpoint m of bo and b2.
Then the line mb^ will intersect our conic in the shoulder point s. This requires
the solution of a quadratic equation,^ but then, using (12.8), we have found our
desired weight Wi\
If the input is not well defined—imagine bo and b2 being on two different
branches of a hyperbola!—then the preceding quadratic equation may have
complex solutions. An error flag would be appropriate here. If the arc between bo
2 The quadratic equation will in general have two solutions. We take the one inside the
triangle bo, bi,b2.

218

Chapter 12 Conic Sections
bi

i

Figure 12.7 Pascal's theorem: the intersection points pi, p2, P3 of the indicated pairs of straight Hnes
are colUnear.

and b2 subtends an angle larger than, say, 120 degrees, it should be subdivided.
For more details, see [619].
Any conic section is uniquely determined by five distinct points in the plane.
If the points have coordinates (x^, y i ) , . . . , (^5,3/5), the implicit form of the
interpolating conic is given by

f(x, y):

x\

xy
xiyi

^2

XlJl

4-\
^]

y\

A

y\
x^y^ y\
^33'3

XsJs

yl

y 1
X]^ yi 1
Xi
yi 1
^3 y3 1
Xj^
y4 1
xs ys 1

The fact that five points are sufficient to determine a conic is a consequence of
the most fundamental theorem in the theory of conies, Pascal's theorem. Consider
six points on a conic, arranged as in Figure 12.7. If w^e connect points as shown,
we form six straight lines. Pascal's theorem states that the three intersection points
pl, p2, p3 are always coUinear.
It can be used to construct a conic through five points: referring to Figure
12.7 again, let a^, bj, Cj, a2, b2 be given (no three of them collinear). Let pj be the
intersection of the two straight lines through a^, b2 and a25 bj. We may now fix a
line 1 through pi, thus obtaining p2 and P3. The sixth point on the conic is then
determined as the intersection of the two straight lines through a^, p2 and b^, P3.
We may construct arbitrarily many points on the conic by letting the straight line
1 rotate around pj.

12.6 Two Classic Problems

219

1 2.6 IWo Classic Problems
A large number of methods exist to construct conic sections from given pieces
of information, most based on Pascal's theorem. A nice collection is given in the
book by R. Liming [391]. An in-depth discussion of those methods is beyond the
scope of this book; we restrict ourselves to the solution of two problems.
1. Conic from two points and tangents plus another point. The given data
amount to prescribing bo^bi, b2. The missing weight Wi must be determined
from the point p, which is assumed to be on the conic. We assume, without loss
of generality, that the conic is in standard form (WQ = tv2 = Vj.
For the solution, we make use of the implicit form (12.21). We can easily
determine the barycentric coordinates TQ, t j , Xi of p with respect to the triangle
formed by the three b/. We can then solve (12.21) for the unknown weight Wi:
wi=

^^

.

(12.22)

If p is inside the triangle formed by bo, b^, b2, then (12.22) always has a solution.
Otherwise, problems might occur (see Problems). If we do not insist on the conic
in standard form, the given point may be given the parameter value t — 1/2, in
which case it is referred to as a shoulder point,
2. Conic from two points and tangents plus a third tangent. Again, we are
given the Bezier polygon of the conic plus a tangent, which passes through two
points that we call bj and bj. We have to find the interior weight u/j, assuming
the conic will be in standard form. The unknown weight Wx determines the two
weight points qo and q^, with q^qi parallel to bob2; see Figure 12.8.
We compute the ratios ro = ratio(bo, b j , b^) and r^ = ratio(bi, b j , b2). From
the definition of the q^ in (12.11), it follows that ratio(bo, qo? b^) = Wi and
ratio(bi, qi, b2) = l/w^. The cross ratio property (12.12) now yields
^=rxwi,

(12.23)

from which we easily determine w^ = y/ro/rl. The number under the square root
must be nonnegative for this to be meaningful (see Problems). Again, if we do
not insist on standard form, we may associate the parameter value t = 1/2 with
the given tangent—it is then called a shoulder tangent.
Figure 12.8 also gives a strictly geometric construction: intersect lines bobj
and b2bo. Connect the intersection with bj and intersect with the given tangent:
the intersection is the desired point p.

220

Chapter 12 Conic Sections

bo
Figure 12.8 Conic constructions: bo, b^, hi-, and the tangent through bj and h\ are given.

1 2.7 Classification
In a projective environment, all conies are equivalent: projective maps map conies
to conies. In affine geometry, conies fall into three classes, namely, hyperbolas,
parabolas, and ellipses. Thus, ellipses are mapped to ellipses under affine maps,
parabolas to parabolas, and hyperbolas to hyperbolas. Hov^ can we determine
vv^hat type a given conic is?
Before we answer that question (following Lee [377]), let us consider the
complementary segment of a conic. If the conic is in standard form, it is obtained
by reversing the sign of wx. Note that the implicit form (12.21) is not affected by
this; hence we still have the same conic, but with a different representation. If c{t)
is a point on the original conic and c{t) is a point on the complementary segment,
one easily verifies that b^, c{t)^ and c{t) are coUinear, as shown in Figure 12.9. If
we assume that W\ > 0, then the behavior of c(^) determines what type the conic
is: if c(^) has no singularities in [0,1], it is an ellipse; if it has one singularity, it
is a parabola; and if it has two singularities, it is a hyperbola.
The singularities, corresponding to points at infinity of c(^), are determined
by the real roots of the denominator w{t) of c{t). There are at most two real
roots, and they are given by

l-\-Wi± Jw\ — 1
^ 2=

•

12.7 Classification

221

Figure 12.9 The complementary segment: the original conic segment and the complementary segment,
both evaluated for all parameter values t e [0,1], comprise the w^hole conic section.

Ellipse

Figure 12.10 Conic classification: the three types of conies are obtained by varying the center v^eight
Wi, assuming WQ = W2 = 1.
Thus, a conic is an ellipse if w/^ < 1, a parabola if wi = 1, and a hyperbola if
Wi > 1. The three types of conies are shown in Figure 12.10 (see also Figure
12.5).
The circle is one of the more important conic sections; let us now pay some
special attention to it. Let our rational quadratic (with Wi].

(14.15)

For the special case [a^ b] = [c, d] = [0,1], we recover the original Bezier points.
Though (14.15) may look complicated, it really is not: all we have to do is to
write a tensor product blossom routine—a matter of about ten lines of code!
Blossoms may also be used to find derivatives, in complete analogy to Section
5.3. Mixed partials take the form

dWdv^

{ni — r)\{n — s)\

Evaluations with respect to the vector 1 in the w-domain are equivalent to
taking differences in the /-direction, those with respect to the z/-vector 1 correspond to differences in the /-direction. Again, it does not matter in which order
we perform the evaluations.
We may use the blossom formulation of derivatives to approach a practical
problem. It is often the case^ that not only a point on a surface is needed, but
also its u- and ^'-partials. Standard tensor product evaluation will give us only
either a w-partial or a i/-partial as a byproduct. However, (14.14) may always be
used to compute both partials. Algorithmically, here is what to do: for a given
2

For applications such as rendering or numerical methods.

260

Chapter 14 Tensor Product Patches
(u^ v) (no blossom notation here), perform evaluation with respect to u for all
rows of control points, but stop all evaluations at level m — 1. This gives us two
points per row. Then perform evaluation with respect to v for the resulting two
columns of points, now stopping at level n — \. We have generated four control
points, corresponding to the bilinear osculant t of (14.14). They may now be
used for evaluation of position and partials. For example, we find the w-partial
as:
^^^^
ou

= m{[(l - v)^Q + vtn] - [(1 - i^)too + t^toi])

with tjj the control points of t(w, v). This approach was first discussed by Mann
and DeRose [413]. See also Sederberg [551].

14.8

Curves on a Surface
Let p = (pj^, pjj) and q = (q^, q^^) be w, z^-coordinates of two points in the domain
of a tensor product patch of degree («, n)? Let
u(t) = (1 - t)p + tq
be the parametric form of a straight line through p and q. This line is mapped to
a curve on the surface. What are its Bezier points c^.'*
Let b [ w i , . . . , w„ I i / i , . . . , f„] be the blossom of the surface. Then a point on
the curve is given by
b[((l - t)p, + ^ , ) < " > I ((1 - t)p, + ^q,)<""].
Applying the Leibniz formula (3.23) to the w-part of the blossom, we can write
it as

Y^ (.^)(l-OVb[p->,q I ((l-Op, + ^q,)<-].
Applying this technique to the v-part also, we get

3 This is not a restriction: we may always degree elevate to achieve equal degrees in u and
in V.

14.9 Normal Vectors

261

Equivalently,
n

n

^=0 /=0

or
In

(n\(n\

fe=0 /+7=f^

Vfe/

Thus

This development is due to T DeRose [159].

14.9 Normal Vectors
The normal vector n of a surface is a normalized vector that is normal to the
surface at a given point. It can be computed from the cross product of any two
vectors that are tangent to the surface at that point. Since the partials d/du and
d/dv 2iTt two such vectors, we may set
n(u, V) = -^^

'\

^^—,

(14.17)

where A denotes the cross product.
At the four corners of the patch, the involved partials are simply differences
of boundary points, for example,
A^'^ooAA^'lboo
||Ai'%,oAA04bo^oll
The normal at one of the corners (we take bg^o ^^ ^^ example) is undefined if
A^'%0,0 ^^^ ^^'^bo^O ^^^ linearly dependent: if that were the case, (14.18) would
degenerate into an expression of the form 0/0. The corresponding patch corner
is then called degenerate. Two cases of special interest are illustrated in Figures
14.11, 14.12, and 14.13.

262

Chapter 14 Tensor Product Patches
n(l,0),
n(0, 0)
c - boo - "10
= " 2 0 = D30

Figure 14.11 Degenerate patches: a "triangular" patch is created by collapsing a whole boundary curve
into a point. The normal at that point may be undefined. Normals are shown for M = 0
and for u=l.

In the first of these, a whole boundary curve is collapsed into a single point. As
an example, we could set boo = b^o = • • • = b^o = c- Then the boundary b(w, 0)
v^ould degenerate into a single point. In such cases, the normal vector at i/ = 0
may or may not be defined. To examine this in more detail, consider the tangents
of the isoparametric lines u = u^ evaluated at v = 0, These tangents must be
perpendicular to the normal vector, if it exists. So a condition for the existence
of the normal vector at c is that all z/-partials, evaluated at f = 0, are coplanar.
But that is equivalent to boi, b ^ , . . . , b^i and c being coplanar.
A second possibility in creating degenerate patches is to allov^ tw^o corner
partials to be coUinear, for example, d/du and d/dv at (0,0), as shov^n in Figure
14.13. In that case, b^o^boi, and boo ^^^ coUinear. Then the normal at boo is
defined, provided that b^^ is not coUinear v^^ith bio?boi, and boo- Recall that
boo?bio?boi, and b ^ form the osculating paraboloid at {u, v) = (0, 0). Then it
follow^s that the tangent plane at boo is the plane through the four coplanar points
boo? bio? boi? and b ^ . The normal at boo is perpendicular to it.
A warning: when we say "the normal is defined," it should be understood that
this is a purely mathematical statement. In any of the preceding degeneracies,
a program using (14.17) will crash. A case distinction is necessary, and then

14.9 Normal Vectors

263

• boo - b i o
• ^ 2 0 = "30

Figure 14.12

Degenerate patches: if all h^ and c are coplanar, then the normal vector at c is perpendicular to that plane.

Figure 14.13

Degenerate patches: the normals at all four corners of this patch are determined by the
triangles that are formed by the corner subquadrilaterals (one corner highlighted).

the program can branch into the special cases that we just described. M o r e
complex situations are encountered when we also w a n t to compute curvatures of
a degenerate patch. A solution is offered in [616]. An a priori check for degenerate
normals is described in [357].

264

Chapter 14 Tensor Product Patches

14.10 IWists
The twist of a surface^ is its mixed partial d^/dudv. According to (14.12), the
twist surface of b'^'^ is a Bezier surface of degree (m — 1, w — 1), and its (vector)
coefficients have the form mwA^'^b/y. These coefficients have a nice geometric
interpretation. The bihnear case is shown in Figure 14.14. In general, the point p/y
is the fourth point on the parallelogram defined by b^ y, b/+i y, b^ y_^i. It is defined
by
Pw-bm,/=k;+l-kr

(14.19)

A^'^b,-,- = (b,+i,,+i - b,+i,) - (b,-,,+i - b,,),

(14.20)

A ' ' \ , = b,+i,+i-p,-,.

(14.21)

Since

it follows that

Thus the terms A^'^b^y measure the deviation of each subquadrilateral of the
Bezier net from a parallelogram.
The twists at the four patch corners determine the deviation of the respective
corner subquadrilaterals of the control net from parallelograms. For example,

dudv

b^'"(0, 0) = mwA^Xo-

(14.22)

This twist vector is a measure for the deviation of b ^ from the tangent plane at
booAn interesting class of surfaces is obtained if all subquadrilaterals b^ y, b^,y_^i,
b/+i,/9 b^_^i y^i are parallelograms; in that case the twist vanishes everywhere. Such
surfaces are called translational surfaces and will be discussed in Section 15.3;
an example is shown in Figure 15.6. They have an interesting shape property: if
all control points of a translational surface lie on the boundary of their convex
hull, then the surface is convex; see Schelske [540]. A surface is convex if it does
not contain a pair of points such that their connection by a straight line intersects
the surface.
4

In this chapter, we are dealing only with polynomial surfaces. For these, the twist is
uniquely defined. For other surfaces, it may depend on the order in which derivatives
are taken; see Section 22.6.

14.11 The Matrix Form of a Bezier Patch

265

Figure 14.14 Twists: the twist is proportional to the deviation of a control net quadrilateral from a
parallelogram.

14.11 The M a t r i x Form of a Bezier Patcii
In Section S.7^ we formulated a matrix expression for Bezier curves. This approach carries over v^ell to tensor product patches. We can write:

boo
[B^{u)

...

•••

\n

(14.23)

BZiu)]
L D^Q

. . .

^mn -I

The matrix {b^;}, defining the control net, is sometimes called the geometry
matrix of the patch. If we perform a basis transformation and write the Bernstein
polynomials in monomial form, we obtain

rb,00
b'^'^CM, v)^[u^

...

bn„1

u'"]M^

W^^
(14.24)

N

L b^o •• • b„„J

L f"J

The square matrices M and N are given by

'-<-"'-'(:)(:

(14.25)

and

,Mru'

«,/ = (-1>

//Vy

(14.26)

266

Chapter 14 Tensor Product Patches
In the bicubic case, m = n = 3^we have

1 -3
3 -1
0
3 -6
3
0
0
3 -3
0
0
1
0

M=N =

For reasons of numerical stabiUty, the use of the monomial form (14.24) is
not advisable (see the discussion in Section 24.3). It is included here since it is
still in widespread use.

14.12 Nonparametric Patches
This section is the bivariate analog of Section 6.5. Having outlined the main ideas
there, we can be brief here. A nonparametric surface is of the form z = f(x^ y).
It has the parametric representation

X(M, V)

u

=

V

L f{u, v) J
and we restrict both u and v to between zero and one. We are interested in
functions f that are in Bernstein form:

A^'3^)=EEMrw^F(y)I

/

Using the identity (5.14) for both variables u and z/, we see that the Bezier points
of X are given by
i/m

jIn
The points (//m, j/n) in the x, y-plane are called Bezier abscissas of the function
f\ the bij are called its Bezier ordinates, A nonparametric Bezier function is not
constrained to be defined over the unit square; if a point p and two vectors v and
w define a parallelogram in the x, y-plane, then the Bezier abscissas a^y G E of a
nonparametric Bezier function over this domain are given by a/y = p + /v + /w.
Figure 14.15 gives an example.

14.13 Problems

267

Figure 14.15 Nonparametric patches: the Bezier points are located over a regular partition of the
domain rectangle.
Integrals also carry over from the univariate case. With a proof analogous to
the one in Section 6.7, we can show that
-1 »i m n

I I TThB'"(x)B"(y)=

yo Jo Y

/

y'"y"bii

^' ^' " .

(m + l)(« + l)

(14.27)

14.15 Problems
1 Draw the hyperbolic paraboloid from Figure 14.2 over the square (—1, —1),
(1, —1), (1,1), (—1,1). Try to do it manually, that is, without graphics
support.
2 Show that the direct de Casteljau algorithm generates surfaces of the form
(14.6). Hint: use blossoms.
3 If a Bezier surface is given by its control net, we can use the de Casteljau
algorithm to compute b'^'^^Cw, v) in three ways: by the direct form from
Section 14.2, or by the two possible tensor product approaches, computing
the coefficients of a w (or v) isoparametric line, and then evaluating that

268

Chapter 14 Tensor Product Patches
curve at 1/ (or w). Though theoretically equivalent, the computation counts
for these methods differ. Work out the details.
* 4 Show that Bezier surfaces have bilinear precision: if hjj = x ( ^ , ^) and x is
bilinear, then b^'"(w, v) = x(w, v) for all w, v and for arbitrary m, n,
* 5 Generalize (5.31) to the tensor product case.
PI Generalize the routine degree_e1 evate to the tensor product case.
P2 Generalize the routine ait ken to the tensor product case, that is, program
tensor product Lagrange interpolation.
P3 The data file car .dat contains data points (slightly modified) of four boundary curves of one of the surfaces shov^n in Color Plate III. Try to fit a Bezier
patch (your pick of the degrees!) so that you get close to the corresponding
surface in the color plates.

Constructing
Polynomial Patches

We' e have discussed the underlying principles of polynomial patches; now it is
time to study how they can be used. First applications of tensor product patches
go back to General Motors, and Boeing in the United States, and to Renault and
Citroen in France.

15.1 Ruled Surfaces
One of the most elementary surface-building schemes is this: let two Bezier curves
h(u) and c(u) be given by their control polygons bo, • . . , b„ and CQ, .. . , c„, and
find a surface r(^, v) that connects them. The simplest such surface will use a
linear type of connection; it is called a ruled surface or lofted surface^ It is given
by
r(w, v) = (1 — v)h(u) -\- vc{u).

(15.1)

For z/ = 0, it interpolates to b(w); for v—l^'\t interpolates to c{u).
The ruled surface r is linear in v and of degree n in u. In Bezier form, it becomes

r(w, v) = [1 — V

v]

bo
Co

B"{u)
lB"„{u)A

1 The word lofted has an interesting history. In the days of completely manual ship design,
full-scale drawings were difficult to handle in the design office. These drawings were stored
and dealt with in large attics, called lofts.
269

270

Chapter 15 Constructing Polynomial Patches

Figure 15.1 Ruled surfaces: corresponding points are connected linearly.

Figure 15.2 Linear interpolation: the average of two convex polygons may not be convex itself.
Figure 15.1 illustrates.
The straight lines on a ruled surface are called its rulings. If all rulings are
parallel, our ruled surface is cylindrical; if all intersect in a point, we have
a conical surface. Both cylindrical and conical surfaces have zero Gaussian
curvature and are special cases of developable surfaces; see Chapter 19 for more
details.
Linear interpolation betw^een curves may not be as intuitive as might be
expected. To see v^hy, consult Figure 15.2. It show^s that the average of tw^o
convex curves (polygons in the case of the figure) is not necessarily convex. Linear
interpolation betv^een curves is thus not shape preserving.

15.2 Coons Patches
Consider the follov^ing design situation: four boundary curves of a surface are
given, all four in Bezier form, that is, by their control polygons. Let us assume that
opposite boundary curves are of the same degrees. The problem: find the control

15.2 Coons Patches

271

Figure 15.3 Coons patches: an example for w = 11 and m = 9.
net of a Bezier surface that fits between the boundary curves. This situation is
illustrated in Figure 15.3.
A configuration of the given data for m = n = 3 looks like this:
boo ^01 bo2 bo3
bio
bi3
b20
b23
bso b3i b32 b33
The construction for the Coons surface mesh relies on the ruled surface form
Section 15.1 as well as the bilinear interpolant from Section 14.1. In order to
find the control point b/y, we first construct points on ruled surfaces:

and

bL = ( l - M b , o + ^,,«
We also need the point on the bilinear interpolant to the four corner points:
hi -•

•-

n

n

bo,o bo,^
bw,0

^m,n

1-zr.
L.

m

-J

The desired control point is found as a combination of these three points:

hi=h,,+h,rhf

(15-2)

Figure 15.4 illustrates the construction method. Figure 15.5 gives an example
of a completed control net.
Coons patches were originally defined in a more general setting, see Chapter
22. Since the Coons technique in this section only deals with piecewise linear

272

Chapter 15 Constructing Polynomial Patches

Figure 15.4 Coons patches: the construction. Gray points, from bottom: b"'2, b" 2? h^ 2* Above them,
solid black: b^ 2-

Figure 15.5 Coons patches: an example.
boundary polygons and control meshes, these are often referred to as discrete
Coons patches.

15.5 Tkranslational Surfaces
A translational surface has the simple structure of being generated by two
curves: let Ci(u) and C2(^') be tw^o such curves, intersecting at a common point
a = ci(0) = C2(0). A translational surface t(w, v) is then defined by

15.3 Translational Surfaces

273

Figure 15.6 Translational surfaces: the Bezier net of a translational tensor product surface. The control
polygons in each direction are translates of each other.
(15.3)

t(u, v) = Ciiu) + C2(i^) - a.

Why the name translational} It is justified by considering an arbitrary isoparametric line of the surface, say, u = u. We obtain t(w, v) = C2(i^) + [—a + Ci(ii)],
that is, all isoparametric lines are translates of one of the input curves; see also
Figure 15.6.
An interesting property of translational surfaces is that their twist is identically
zero everywhere:

dudv

t(u, v) = 0.

This property follows directly from the definition (15.3). Since both input curves
may be arbitrarily shaped, the resulting surface may well have high curvatures.
This dispels the myth that zero twists are identical to flat spots. In fact, twists
are not related to the shape of a surface—rather, they are a result of a particular
parametrization. See also Section 16.3 on twist generation.
A translational surface may be viewed as the solution to an interpolation
problem: given two intersecting curves, find a surface that contains them as
boundary curves. If four boundary control polygons are given, as in the problem
definition for a Coons patch, we can form four translational surfaces, one for
each corner. Let us denote by tjj the translational surface that interpolates to the
boundary curves meeting at the corner (/,/); /,; € {0,1}.
Now the Coons patch x(w, v) can be written as
x(^, v) = [1 — u

u]

too(^, v)

toi(w, v)

1-v
V

(15.4)

This form of the Coons patch is called a convex combination. It blends together four surfaces, weighting each with a weight function. The weight functions

274

Chapter 15 Constructing Polynomial Patches
sum to one (a necessity: nonbarycentric combinations are disallowed) and are
nonnegative for u,v e {0,1}. Note that the weight functions are zero where the
corresponding t^y is "wrong."
Translational and Coons surfaces are related in yet another way. If four
boundary control polygons are related such that opposite polygons are translates
of each other, then the resulting Coons control net will be a translational surface;
an example is provided by Figure 15.6.

1 5.4 Tensor Product Interpolation
We could use curves for both free-form design and for interpolation; the same is
true for tensor product patches. As a preparatory step, let us rewrite (14.6) in an
equivalent matrix form:
^00

x{u,v) = [B^(u)

^On

Bl{v)
(15.5)

BZiuy
.B>),

LD^O

Suppose now that we are given an (m + 1) x (w + 1) array of data points
X/y; 0 \ o ^
5Vi =
0.1
0.1
1.1
L 1-1

This surface does not go through any of the data points exactly, but it is reasonably close to them.

15.7 Finding Parameter Values

281

to solve a linear system with 16 equations in 16 unknowns. A note of caution: if
the number of data points is very large (10^ or more), then the normal equations
become ill-conditioned and the least squares problem may not be solvable.

15.7 Finding Parameter Values
In a practical setting, we would not typically be given the parameter values
u^ = i'^h ^k)- Finding good values for the u^ is not an easy problem.
But often, the data points can be projected into a plane; let us assume they
can be projected into the x, }/-plane for simplicity. Each p^ is projected by simply
dropping its ;^-coordinate, leaving a pair (x^, y^). We scale the (x^, y^) so that
they fit into the unit square. Then we can simply set uj^ = xj^ and Vk^JkAnother solution may be obtained by using a triangulation of the data points.
This scenario is realistic for data obtained using a laser digitizer. Assuming
that the triangulation is isomorphic to the unit square, we can construct a
triangulation in the unit square with the same connectivity as the given one
in 3D.
The following method is due to Floater [236]. First, a convex polygon is built
in the (^, v) unit square with as many vertices as the 3D mesh has boundary
vertices. This polygon is somewhat arbitray; a circle or the boundary of the unit
square is a good candidate for forming it. In this way, we assign 2D parameters
to the 3D mesh boundary points.
Next, consider any interior point u^ of the 2D mesh with Uj. neighbors. We call
these neighbors u^ i,. . ., u^ „ . For a "nice" triangulation, the following condition
should be satisfied: all interior u^. may be expressed as^

1 "^

We now observe that there are as many equations (15.19) as there are interior
points in the mesh. Some of these equations involve boundary points, others will
not. Thus we have a linear system for the interior u^, with as many equations as
there are interior points. It is always solvable.

2 Triangulations with this property are the result of Laplacian smoothing of a mesh.

282

Chapter 15 Constructing Polynomial Patches

1 5.8 Shape Equations
The solution to a least squares problem may be close to the data points, yet the
control net might "behave badly," similar to the curve case. A way to combat
such behavior is to invoke shape equations. These are conditions that a "good"
control net v^ould satisfy. Out of many possibilities, we choose the following: a
translational surface (see 15.3) is characterized by the fact that its twist vanishes
everywhere. For the control net, this means

K; = 0;

. 1,1]

/ = 0 , . . . , m - 1, / = 0 , . . . , w - 1.

When we add these equations to our overdetermined linear system (15.17), we
will be less faithful to the data points, but achieve a better shape of the control
net. In practice, we would weight the shape equations, maybe by a factor of 0.1,
just as we did for the curve case.

1 5.9 A Problem with Unstructured Data
The least squares approximation problem for Bezier patches leads to an interesting question: what happens if the number of data points equals the number of
unknowns.^ In this case, we do not have to use the normal equations; the problem
directly yields a linear system with a square coefficient matrix.
Although this interpolation scenario appears simpler than the least squares
problem, it has a potential for a serious pitfall.
Our linear system is given by

Po

B^{UQ)BI{VO)

K(^O)B:(VO)

rb,'0,0
, (15.20)
L '^rn,n -•

LPKJ

.B^(UK)B-^{VK)

BZ(UK)B-{VK)^

but now we have K = {m -\- l)(n -\-1).
The assignment of parameter values (^^, Vj) to the data points p/ is (theoretically) arbitrary. We now perform a "thought experiment."^ Take two data points
PI and py. They have parameters (w^, Vj) and (wy, Vj), Now, leaving the data points
untouched, start changing their parameter values. Do this as follows: find the
circle having (w/, Vj) and (wy, Vj) as diameter and move both (w/, Vi) and (wy, Vj)
3

Or "Gedankenexperiment" as coined by A. Einstein.

15.10 Implementation

283

Figure 15.9 Unstructured data: two parameter locations u^ = (w^, Vj) and Uy = (wy, Vj) are interchanged.
along this circle in a counterclockwise fashion until they have interchanged their
locations. Figure 15.9 illustrates.
The determinant of the linear system (15.20) has now changed its sign: two
of its rows are interchanged. It follows that somewhere during the interchanging
procedure, the determinant must have been zero, corresponding to an unsolvable
interpolation problem.
For an arbitrary collection of data point parameters, we cannot know if we
might be in or near an unsolvable situation. Hence this interpolation problem is
ill-posed.

15.10 Implementation
The following is the header for a program to plot a tensor product Bezier surface,
in fact, a rational one. If the polynomial case is desired, just set all weights to
unity.
void
/*

Input:

plot_ratsurf(bx,by,bw,degree_u,degree_v,u_points,v_points,
seale_x,seale_y)
plots v_points isoparametric
curves of the rat Bez surface, each with u_points
points on it.
bx, by:

arrays with x- and y- coordinates of
control net.
degree_u,degree_v: degrees in u- and v- direction
u_points,v_points: plot resolution
seale_x,seale_y:
scale factor for postscript.

Output: postscript file

284

Chapter 15 Constructing Polynomial Patches
Next is a routine that fits a biHnearly blended Coons patch in between four
boundary control polygons, as described in Section 15.2. The routine works on
one coordinate only, and will have to be called separately for the x-, y-, and
^-components of a control net.
void netcoons(net,rows,columns)
/ * Uses bilinear Coons blending to complete a control
net of which only the four boundary polygons are used as input.
Works for one coordinate only.
Input:
net:
control net.
rows, columns: dimensions of net.
Output:
net:
the completed net, with the old boundaries.

*/

15.11 Problems
1 Justify that in tensor product interpolation (Section 15.4), it does not
matter if we start with the row interpolation process or with the column
interpolation process. Give computation counts for both strategies. (In
general, they are not equal!)
2 Generalize Lagrange interpolation to the tensor product case.
* 3 Generalize quintic Hermite interpolation to the tensor product case.
* 4 Suppose we want to find a parametrization {ui) for a tensor product interpolant. We may parametrize all rows of data points and then form the
averages of the parametrizations thus obtained. Or we could average all
rows of data points, for example, by setting p^ = ^- ^x^ y and we could
then parametrize the p^. Do we get the same result? Discuss both methods.
PI Figure 7.3 shows how Lagrange interpolation behaves badly for curves.
Write a program to exhibit this effect for tensor product patches.

Composite
Surfaces

1 ensor product Bezier patches were under development in the early 1960s; at
about the same time, people started to think about piecewise surfaces. One of the
first publications was de Boor's work on bicubic splines [136] in 1962. Almost
simultaneously, and apparently unaware of de Boor's work, J. Ferguson [231]
implemented piecewise bicubics at Boeing. His method was used extensively,
although it had the serious flaw of using only zero corner twist vectors. An
excellent account of the early industrial use of piecewise bicubics is the article by
G. Peters [462].

16.1 Smoothness and Subdivision
Let x(w, v) and y(w, v) be two patches, defined over [w/_i, uj] x [vj, Vj^i] and
[w/, w/+i] X [vj, ^'/+i], respectively. They are r times continuously differentiable
across their common boundary curve x(wj, v) = y(w/, v) if all w-partials up to
order r agree there:

Z •x(u, v)
dw

(16.1)
u—ui

Now suppose both patches are given in Bezier form; let the control net of
the "left" patch be {b^y}; 0

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.6
Linearized                      : No
Page Mode                       : UseOutlines
XMP Toolkit                     : 3.1-701
Modify Date                     : 2009:01:12 11:53:41+01:00
Create Date                     : 2007:07:04 23:43:37Z
Metadata Date                   : 2009:01:12 11:53:41+01:00
Format                          : application/pdf
Title                           : Curves and Surfaces for CAGD: A Practical Guide
Creator                         : Gerald Farin
Description                     : The Morgan Kaufmann Series in Computer Graphics
Subject                         : ISBN-13: 9781558607378
Document ID                     : uuid:5f3c8bdc-e027-4baa-918d-e399c940fa70
Instance ID                     : uuid:5de36d88-8c1a-4b57-9d45-0869eab4e2cc
Producer                        : ABBYY FineReader 8.0 Professional Edition
Has XFA                         : No
Page Count                      : 521
Page Layout                     : SinglePage
Author                          : Gerald Farin
Keywords                        : ISBN-13:, 9781558607378
Warning                         : [Minor] Ignored duplicate Info dictionary