From Spiral to Spline: Optimal Techniques in Interactive Curve DesignbyRaphael Linus Levien

A dissertation submitted in partial satisfactionof the requirements for the degree ofDoctor of Philosophy
inEngineering–Electrical Engineering and Computer Sciencesin theGRADUATE DIVISIONof the
UNIVERSITY OF CALIFORNIA, BERKELEY

Committee in charge:
Professor Carlo Séquin, Chair
Professor Jonathan Shewchuk
Professor Jasper Rine

Fall 2009

From Spiral to Spline: Optimal Techniques in Interactive Curve Design

Copyright© 2009byRaphael Linus Levien

Abstract

From Spiral to Spline: Optimal Techniques in Interactive Curve DesignbyRaphael Linus Levien
Doctor of Philosophy in Engineering–Electrical Engineering and Computer Sciences
University of California, Berkeley
Professor Carlo Séquin, Chair
A basic technique for designing curved shapes in the plane is interpolating splines. The designer
inputs a sequence of control points, and the computer fits a smooth curve that goes throughthese points. The literature of interpolating splines is rich, much of it based on the mathematical
idealization of a thin elastic strip constrained to pass through the points. Until now there is little
consensus on which, if any, of these splines is ideal. This thesis explores the properties of an idealinterpolating spline. The most important property is fairness, a property often in tension withlocality, meaning that perturbations to the input points do not affect sections of the curve at
a distance. The idealized elastic strip has two serious problems. A sequence of co-circular input
points results in a curve deviating from a circular arc. For some other inputs, no solution (withfinite extent) exists at all.
The idealized elastic strip has two properties worth preserving. First, any ideal spline must be
extensional, meaning that the insertion of a new point on the curve shouldn’t change its shape.
Second, curve segments between any two adjacent control points are drawn from a two-parameterfamily (and this propety is closely related to good locality properties). A central result of thisthesis is that any spline sharing these properties also has the property that all segments between
two control points are cut from a single, fixed generating curve. Thus, the problem of choosing an
ideal spline is reduced to that of choosing the ideal generating curve. The Euler spiral has excellent
all-around properties, and, for some applications, a log-aesthetic curve may be even better.
Shapes in applications such as font outlines contain extra features such as corners and transitions between straight lines and smooth curves. Attaching additional constraints to control points
expresses these features, and, carefully applied, give the designer a richer palette of curve types.
The splines presented in this thesis are entirely practical as well, especially for designing fonts.
Sophisticated new numerical techniques compute the splines at interactive speeds, as well as convert
to optimized cubic Bézier representation.

Professor Carlo Séquin
Dissertation Committee Chair

1

To my father.

iContents
Contents

ii

List of Figuresvii

List of Tablesxi

Acknowledgements

xii

1 Introduction

1

2 Properties of splines

6

Fairness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.1.1

Empirical testing of fairness . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Roundness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.4

Extensionality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.5

Existence and uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.6

Locality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.7

Transformational invariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.8Counting parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.9

Monotone curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.1

2.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 The Elastica

16

3.1

Statement of the elastica problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2

Variational study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.1

The Euler–Lagrange equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2.2

Application of the Euler–Lagrange equation to the elastica . . . . . . . . . . 18
ii

3.3

A family of solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.4

Equilibrium of moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.5

Pendulum analogy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.6

Equilibrium of forces: Finite element approach . . . . . . . . . . . . . . . . . . . . . 25

3.7

Solutions with length and endpoint constraints . . . . . . . . . . . . . . . . . . . . . 29

3.8

Elastica with general constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.9

Angle constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.10 Elastica as an interpolating spline. . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.11 SI-MEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.12 Real splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4 Two-parameter splines

42

4.1

Extensional two-parameter splines have a generating curve . . . . . . . . . . . . . . . 43

4.2

Properties of the MEC spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.3

SI-MEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.4

Cubic Lagrange (with length parametrization) . . . . . . . . . . . . . . . . . . . . . . 47

4.5

Ikarus (biarc) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.6

Euler’s spiral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.6.1

Euler’s spiral as a spline primitive . . . . . . . . . . . . . . . . . . . . . . . . 51

4.6.2

Variational formulation of Euler’s spiral . . . . . . . . . . . . . . . . . . . . . 51

4.6.3

Euler spiral spline in use . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.6.4

Existence and uniqueness of Euler spiral spline . . . . . . . . . . . . . . . . . 53

4.7

Hobby’s spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.8

Splines from fixed generating curves . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.8.1

4.9

Log-aesthetic curve splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Piecewise SI-MEC splines, or Kurgla 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.10 Circle splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5 History of the elastica

67

5.1

Jordanus de Nemore – 13th century . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2

Galileo sets the stage – 1638 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.3

Hooke’s law of the spring – 1678 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.4James Bernoulli poses the elastica problem – 1691 . . . . . . . . . . . . . . . . . . . 70

iii

5.5

James Bernoulli partially solves it – 1692 . . . . . . . . . . . . . . . . . . . . . . . . 71

5.6James Bernoulli publishes the first solution – 1694 . . . . . . . . . . . . . . . . . . . 73

5.7

Daniel Bernoulli proposes variational techniques – 1742 . . . . . . . . . . . . . . . . 77

5.8

Euler’s analysis – 1744 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.8.1

5.9

Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Elliptic integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

5.10 Laplace on the capillary – 1807 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.11 Kirchhoff’s kinetic analogy – 1859 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.12 Closed form solution: Jacobi elliptic functions – 1880 . . . . . . . . . . . . . . . . . . 89
5.12.1 Inflectional solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.12.2 Non-inflectional solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.13 Max Born – 1906 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.14 Influence on modern spline theory – 1946 to 1965 . . . . . . . . . . . . . . . . . . . . 92
5.15 Numerical techniques – 1958 through today . . . . . . . . . . . . . . . . . . . . . . . 93
6 History of the Euler spiral

95

6.1

James Bernoulli poses a problem of elasticity – 1694 . . . . . . . . . . . . . . . . . . 95

6.2

Euler characterizes the curve – 1744 . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

6.3

Euler finds the limits – 1781 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6.4

Relation to the elastica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.5

Fresnel on diffraction problems – 1818 . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.6

Talbot’s railway transition spiral – 1890 . . . . . . . . . . . . . . . . . . . . . . . . . 106

6.7

Mathematical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

6.8

Use as an interpolating spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

6.9

Efficient computation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7 Four-parameter curves

110

7.1

Why four parameters? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

7.2

Minimum Variation Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
7.2.1

Euler–Poisson equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

7.2.2

Euler–Poisson solution of MVC . . . . . . . . . . . . . . . . . . . . . . . . . . 115

7.3

Polynomial spiral spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

7.4

Generalized constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119iv

8 Numerical toolbox
8.1

122

Numerical integration of polynomial spirals . . . . . . . . . . . . . . . . . . . . . . . 123
8.1.1

Polynomial approximation of the integral . . . . . . . . . . . . . . . . . . . . 124

8.1.2

Higher-order polynomial approximations . . . . . . . . . . . . . . . . . . . . . 125

8.1.3

Hybrid approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

8.2

Determining Euler spiral parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

8.3

Global constraint solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

8.4

8.3.1

Two-parameter solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

8.3.2

Four-parameter solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

2-D Lookup table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

9 Converting to Bézier curves

135

9.1

Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

9.2

Error metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

9.3

Simple conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

9.4

Equal arc length, equal arm length solution . . . . . . . . . . . . . . . . . . . . . . . 139

9.5

Fully optimized solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

9.6

Global optimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
9.6.1

Hybrid error metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

9.6.2

Simple subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

9.6.3

Smarter subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

9.6.4

Alternative optimum algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 147

10 Applications in font design

149

11 Space curves

160

11.1 Curves in space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
11.2 Minimum energy curve in space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
11.2.1 MEC as spline in 3-space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
11.3 Mehlum’s spiral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
11.4 Log-aesthetic space curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
11.5 Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
12 Conclusion

167v

Bibliography

171

vi

List of Figures
1.1

A mechanical spline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2

Rectangular elastica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3

Euler spiral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4

Conversion to optimized Bézier curves. . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1

Test instrument for fairness user study. . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

Survey results for aesthetic curve preference. . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Continuity 6= fairness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.4

MEC roundness failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.5

Circle spline extensionality failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.6

Nonexistence of solution to MEC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.7

Uniqueness failure in Euler spiral spline. . . . . . . . . . . . . . . . . . . . . . . . . 12

2.8

G2 spline has better locality than G4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.9

Attneave’s curvature experiment [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1

The family of elastica solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.2

Curvature plots for the family of elastica solutions. . . . . . . . . . . . . . . . . . . 22

3.3

Equilibrium of moments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.4

An oscillating pendulum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.5

Examples of pendulum analogy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.6

Finite element approximation of elastica as chain of struts and pivots. . . . . . . . . 27

3.7

Forces on a single strut. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.8

Forces on a single pivot.

3.9

Forces on a single pivot point in a chain. . . . . . . . . . . . . . . . . . . . . . . . . 28

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.10 Values of λ for varying chord lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.11 Multiplicity of solutions for same tangent constraints. . . . . . . . . . . . . . . . . . 31

vii

3.12 Forces acting on rectangular elastica. . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.13 Angle-constrained minimum energy curves. . . . . . . . . . . . . . . . . . . . . . . . 33
3.14 Curvature plots, angle-constrained minimum energy curves. . . . . . . . . . . . . . . 34
3.15 MEC and SI-MEC energies for varying lengths (π/2). . . . . . . . . . . . . . . . . . 36
3.16 MEC and SI-MEC energies for varying lengths (π/3). . . . . . . . . . . . . . . . . . 37
3.17 SI-MEC and MEC curve energy for (π/2, −π/2) angle constraints. . . . . . . . . . . 38
3.18 A real spline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.19 High-friction spline constraints on a real spline. . . . . . . . . . . . . . . . . . . . . 40
3.20 Bricks as high friction constraints on a real spline. . . . . . . . . . . . . . . . . . . . 41
4.1

Incremental construction of the generator curve of a two-parameter spline. . . . . . 43

4.2

MEC roundness failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.3

Typical letter-shape using Ikarus. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.4

Euler’s spiral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.5

Euler’s spiral spline primitive over its parameter space. . . . . . . . . . . . . . . . . 50

4.6

Closed Euler example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.7

Complex Euler example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.8

Difficulty with straight-to-curved joins. . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.9

Multiple windings for Euler spiral segment connecting two points (from [82]). . . . . 56

4.10 Hobby spline primitive over its parameter space. . . . . . . . . . . . . . . . . . . . . 58
4.11 Curvature plot of Hobby spline over its parameter space. . . . . . . . . . . . . . . . 59
4.12 Test instrument for fairness user study. . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.13 Perturbation fall-off rate as a function of aesthetic curve exponent. . . . . . . . . . 62
4.14 MEC, SI-MEC, and Euler spiral compared. . . . . . . . . . . . . . . . . . . . . . . . 64
4.15 Circle spline extensionality failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.1

Galileo’s 1638 problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2

Hooke’s figure on compound elasticity. . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.3

Bernoulli poses the elastica problem. . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5.4

Drawing from James Bernoulli’s 1692 Med. CLXX, with modern reconstruction. . . 72

5.5

Bernoulli’s 1694 publication of the elastica. . . . . . . . . . . . . . . . . . . . . . . . 73

5.6

Bernoulli’s justification for his formula for curvature. . . . . . . . . . . . . . . . . . 74

5.7

A simplified diagram of moments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.8

Huygens’s 1694 objection to Bernoulli’s solution. . . . . . . . . . . . . . . . . . . . . 77

viii

5.9

Euler’s elastica figures, Tabula III. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.10 Euler’s elastica figures, Tabula IV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.11 The family of elastica solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.12 Bernoulli’s lemniscate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.13 Maxwell’s figures for capillary action. . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.14 An oscillating pendulum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.15 Rectangular elastica as roulette of hyperbola. . . . . . . . . . . . . . . . . . . . . . . 90
5.16 Born’s experimental apparatus for measuring the elastica.

. . . . . . . . . . . . . . 92

6.1

Euler’s spiral as an elasticity problem. . . . . . . . . . . . . . . . . . . . . . . . . . . 96

6.2

Bernoulli’s construction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

6.3

Euler’s drawing of his spiral, from Tabula V of the Additamentum. . . . . . . . . . 98

6.4

Reconstruction of Euler’s Fig. 17, with complete spiral superimposed. . . . . . . . . 99

6.5

The figure from Euler’s “De valoribus integralium”, with modern reconstruction. . . 101

6.6

Euler’s spiral, rectangular elastica, and cubic parabola. . . . . . . . . . . . . . . . . 102

6.7

Diffraction through a slit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.8

Cornu’s plot of the Fresnel integrals. . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

6.9

Talbot’s Railway Transition Spiral. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

7.1

Eccentricity 0.56199 oval with G2 and G4 continuity. . . . . . . . . . . . . . . . . . 111

7.2

The suitcase corners problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

7.3

Bracketed serifs and straight-to-curve joins. . . . . . . . . . . . . . . . . . . . . . . . 113

7.4

One-way constraints enforce locality.

7.5

G2 spline has better locality than G4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

7.6

Matching nonzero curvatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

7.7

A one-way constraint creates a cyclic dependency. . . . . . . . . . . . . . . . . . . . 120

7.8

Examples of different kinds of constraint points. . . . . . . . . . . . . . . . . . . . . 120

8.1

Precision vs. timing for different integration strategies. . . . . . . . . . . . . . . . . 126

8.2

Determining Euler spiral parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . 128

8.3

Python code for computing Euler spiral parameters. . . . . . . . . . . . . . . . . . . 130

8.4

Pseudocode for computing four-parameter splines. . . . . . . . . . . . . . . . . . . . 132

8.5

Example of four-parameter solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

8.6

Contour plot of curvature lookup table for Euler spiral. . . . . . . . . . . . . . . . . 134. . . . . . . . . . . . . . . . . . . . . . . . . . 113

ix

9.1

A cubic Bézier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

9.2

Two circle approximations with equal distance error. . . . . . . . . . . . . . . . . . 137

9.3

Euler spiral approximation drawn using .03 and .001 error threshold, respectively. . 138

9.4

Three Bézier solutions: fast, equal arc length, and fully optimized. Curvature plotis the top graph, the Bézier is on bottom. . . . . . . . . . . . . . . . . . . . . . . . . 139

9.5

Approximation errors over parameter space for cubic Bézier fit. . . . . . . . . . . . 141

9.6

Local minima of cubic Bézier fitting, t=[0.5, 1.0]. Curvature plots are the top twographs, the Béziers are on the bottom. . . . . . . . . . . . . . . . . . . . . . . . . . 142

9.7

Local minima of cubic Bézier fitting, t=[0.5, 0.85]. Curvature plots are the top
graphs, the Béziers are on the bottom. . . . . . . . . . . . . . . . . . . . . . . . . . 143

9.8

Conversion into optimized cubic Béziers. . . . . . . . . . . . . . . . . . . . . . . . . 144

9.9Four different levels of optimization: subdivision in half with equal arm length,
subdivision in half with fully optimized Bézier segments, greedy subdivision withfully optimized segments, and fully optimized. . . . . . . . . . . . . . . . . . . . . . 146

10.1 Selected glyphs from Inconsolata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
10.2 Selected glyphs from Inconsolata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
10.3 Selected glyphs from Baskerville Italic. . . . . . . . . . . . . . . . . . . . . . . . . . 153
10.4 Selected glyphs from Baskerville Italic. . . . . . . . . . . . . . . . . . . . . . . . . . 154
10.5 Lowercase ‘w’ from Baskerville Italic, scan with outline. . . . . . . . . . . . . . . . . 155
10.6 Selected glyphs from Cecco. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
10.7 Selected glyphs from Cecco. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
10.8 Typical letter-shape using IKARUS, with Spiro comparison. . . . . . . . . . . . . . 158
10.9 Interpolation masters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
10.10 Interpolation between extremes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
10.11 Interpolation between extremes, superimposed. . . . . . . . . . . . . . . . . . . . . . 159
10.12 Interpolation between cubic Béziers fails to preserve G1 -continuity. . . . . . . . . . 159
11.1 Space curves rendered into ironwork. . . . . . . . . . . . . . . . . . . . . . . . . . . 160
11.2 Curvature as a function of arc length in the Mehlum spiral. . . . . . . . . . . . . . . 164
12.1 Contour plot of curvature lookup table for α = 2 log-aesthetic curve. . . . . . . . . 168

x

List of Tables
4.1

Comparison of properties of two-parameter splines. . . . . . . . . . . . . . . . . . . 66

5.1

Euler’s characterization of the space of all elastica. . . . . . . . . . . . . . . . . . . . 80

7.1

Constraints for curve points, and their corresponding formulae. . . . . . . . . . . . . 121

10.1 Constraints used to draw font outlines. . . . . . . . . . . . . . . . . . . . . . . . . . 152

xi

Acknowledgements
Thanks to Derek Waters for the image of the spline and weights in Figure 1.1. The photograph ofarchitectural ironwork shown in Chapter 11 is by Terry Ross of drawmetal.com, as is the original
work.
This historical sketch in Chapter 5 draws heavily on Truesdell’s history in an introduction to Ser
2, Volume X of Euler’s Opera Omnia [86]. Indeed, a careful reading of that work reveals virtually
all that needs to be known about the problem of the elastica and its historical development throughthe time of Euler.
I am deeply grateful to Alex Stepanov, Seth Schoen, Martin Meijering and Ben Fortson for help
with the translations from the Latin in Chapters 5 and 6. Special mention is also due the excellent“Latin words” program by William Whitaker, which lets anyone with a rudimentary knowledge of
Latin and a good deal of patience and determination puzzle out the meaning of a text. Of course,
all errors in translation are my own responsibility.
Patricia Radelet has been most helpful in supplying high-quality images for James Bernoulli’s
original figures, reproduced here as Figures 5.4 and 5.5, and for providing scans of the original textof his 1694 publication.
Thanks also to the Stanford Special Collections Library for giving me access to Prof. Forsythe’snotes and correspondence, which were invaluable in tracing the influence of the elastica throughearly work on non-linear splines.
Thanks to Karl Berry for supporting the development of the Inconsolata font, and for building
the TEX configuration files so it could easily be used in this thesis.
I gratefully appreciate Prof. Norimasa Yoshida for a careful reading and helpful comments onan earlier draft of this thesis.
And of course, many thanks are due to my thesis committee for providing feedback invaluable
for refining earlier drafts, and especially my advisor, Prof. Carlo Séquin.

xii

xiii

Chapter 1

Introduction
A basic design task is drawing smooth curves in the plane. Curves are needed in purely aesthetic applications such as making vector drawings, purely industrial applications such as the
cross-sections of airplane wings or railway track layouts, and also applications both functional andartistic, such as font design.
There are several approaches to curve design. The most popular is drawing with Bézier curves,usually cubic. These curves are versatile and very easy to compute and manipulate, but they are
somewhat difficult to learn, and it is especially difficult to get smooth curves. In expert hands, theresults are good, but for less experienced designers, adding more control points to tweak a curve
often leads to more lumpiness.
Another approach is to sample real-world data (for example, from a sculpture), then fit a curve
to it. Since these points are noisy, it’s most common to fit an approximating spline through them
– the spline goes near the points but not necessarily through them.
The main focus of this thesis is interpolating splines, or curves that go through the points given,in sequence. With a good interpolating spline, the designer can enter a relatively sparse sequenceof data points, and the result is a smooth curve.
Interpolating splines have their roots in the mechanical spline, a thin strip of flexible material
constrained to go through the chosen points with weighted implements (called “ducks” or “whales”
due to their shape), but otherwise free. An example is shown in Figure 1.1. Such a strip will
naturally find its minimal energy configuration. Indeed, one approach to interpolating splines is
simply to simulate the mechanical spline. As we shall see, such an approach is severely flawed.
A part of the problem with interpolating splines is that there are many plausible constructions,with quite a few in actual use, but no real consensus on which one is best, or even how to adequately
measure or compare. Many of the interpolating splines in use seem to have been chosen based on
expedience, simple representation, or cheap computational cost. However, these criteria are notreally justified, not only because Moore’s law has made even intricate computations feasible atinteractive speeds, but also because excellent practical splines exist with a simple representation
and an efficient implementation.
This thesis will closely examine the question of what interpolating spline is ideal out of all1

Figure 1.1. A mechanical spline.

possibilities. Such a focus reduces the space of alternatives to consider, excluding all the manyrough and crude approximations.
By far, the most important property of an interpolating spline is fairness, (also known, lesstechnically, as smoothness). Indeed, one way to define an interpolating spline is to choose a fairness metric (such as the total bending energy of a flexible strip) and optimize it, given that the
curve passes through the control points. However, this approach may compromise other important
properties, such as robust convergence, locality (intuitively, lack of ripples), stability, and others.
Chapter 2 lists all these desirable properties of interpolating splines and describes them in detail.
The mechanical spline is the basis for essentially all following work in splines. Its mathematical
idealization is the Minimum Energy Curve (MEC), and the curve describing segments between
adjacent control points is the elastica. This curve has surprising and deep properties, and there are
a number of ways to fully understand it – as the result of an energy minimization, as the limit ofa finite element system consisting of rigid beams connected by flexible springs, and as the physicalanalogue of a pendulum (which explains why the curve is periodic). Chapter 3 explores all thesecharacterizations and the connections between them.
Yet, while the MEC is the optimum of a particular functional (bending energy), it is clearly
not the best possible spline. Experiment shows that bending energy does not accurately predict
perceived smoothness. This finding is also supported by theoretical analysis – in particular, the
MEC spline lacks roundness, the property that control points arranged on a circular arc yields
that arc. Further, optimizing this particular metric yields a spline with poor robustness – for many
inputs, no optimum solution exists at all. Even so, the MEC has many desirable properties, some ofwhich are worth preserving. Two in particular stand out. An optimum spline must be extensional,
meaning that adding an additional control point on the curve of an existing spline preserves theshape of the curve. Further, the MEC has the property that the family of curve segments between
any two adjacent control points is a manifold with a two-dimensional parameter space. Chapter
2

4 argues that, for two-dimensional drawing applications, any “best” spline should also have thisproperty.

Figure 1.2. Rectangular elastica.

Figure 1.3. Euler spiral.

Chapter 4 also presents one of the central technical results of this thesis: any spline that is both
extensional and two-parameter is derived from a single, rigid generating curve, with each segment
cut from the curve, then rotated, scaled, and translated into place to align with the endpoints.
In the MEC, this curve is the rectangular elastica (shown in Figure 1.2), but other choices yield
better splines. In particular, the Euler spiral (shown in Figure 1.3) produces results similar to the
MEC, but fixing the roundness and robustness problems. Chapter 4 then concludes by exploring
the potential for the family of log-aesthetic curves to produce splines even slightly “better” than the
Euler spiral spline, improving simultaneously fairness (as determined by empirical measurement)
and locality.
The two basic underlying curves for two-parameter splines, the elastica and the Euler spiral,
both have rich histories, dating back hundreds of years. But they have not been fully explainedin one concise source until now. For both of these curves, Jacob Bernoulli wrote down the precise
mathematical description and partially explored the nature of the curves, while Euler characterized
them fully. Chapter 5 deeply explores the history of the elastica, including its subsequent role infounding the field of study of elliptical integrals, which are now also useful for giving a closed-form,
analytical solution to the shape of the curve. Similarly, Chapter 6 explores the history of the Euler
spiral, including its repeated rediscovery and applications in other diverse fields, such as solving
diffraction problems and layout of railroad tracks.
3

thresh=0.01 optim=3 tot segs=18

Figure 1.4. Conversion to optimized Bézier curves.

Two-parameter splines are ideal as basic interpolating splines used for planar curve design, butthere are applications where a higher degree of geometric continuity is desired. Chapter 7 presents
the logical extension to a four-parameter space, providing G4 -continuity, but at the cost of reduced
locality. The underlying curve family is also useful for satisfying a richer family of constraints thanmerely that the curve passes through the control points. In particular, such constraints can specify
smooth transitions from straight sections to curves, a frequent pattern in font designs. The chapter
also explores the connection of this spline to the Minimum Variation Curve, a known spline basedon minimizing a metric defined in terms of the variation of curvature rather than curvature itself,
as in the MEC.
None of these theoretical results about optimal splines would be much use without an efficient,practical implementation. Chapter 8 provides a toolbox of such techniques, including: efficiently
computing the integral defining the Euler spiral and its generalization to a four-parameter spline;
solving the parameters of this curve to meet given tangent constraints; and providing an efficient
Newton-style solver for simultaneously solving all the constraints required for a given spline. Together, these numerical techniques add up to a very fast solver for splines, entirely suitable forinteractive work even in limited computational environments.
Even after solving the global spline and determining an analytical curve that meets the constraints, the resulting curve is still not in a form that can be widely used. Euler spirals (and theircubic polynomial generalization, even more so) are considered “special functions” and are very
unlikely to be present as primitives in graphics files formats, CAD applications, and so on. Thus,
another important component for actually using these splines in practice is conversion to a more
widely-supported format, of which by far the most popular is cubic Bézier curves. Chapter 9 describes how to convert these curves into Bézier form, with a tradeoff between the computational
effort expended on optimization and the conciseness of the representation. For contexts when thesize of the representation matters (for example, when downloading fonts over the network), thechapter presents a fully optimized technique, yielding the best possible Bézier representation with4

respect to the given error threshold. Figure 1.4 shows an example of a letterform converted intosuch an optimized Bézier representation.
A major motivation for this work was to improve tools for designing fonts. I designed several
complete fonts entirely using the splines described in this thesis, and one, Inconsolata, is released andis widely used. Chapter 10 presents some of these font designs, and also discusses particular aspectsof splines relevant to the specific application of font design, including prospects for more easily
creating fonts with variation, for example, a smooth blending between light and heavy weights.
While the primary focus of this thesis is planar curves, the principles are also useful for splines inspace, as explored in Chapter 11. An intrinsic formulation of planar curves is a relationship between
arc length and curvature. The corresponding intrinsic form for a space curve is a relationshipbetween arc length and both curvature and torsion. Three parameters are sufficient to attain G2 continuity of space curves, as opposed to two for planar curves. The MEC generalizes easily fromplanar to space curves, but of course has the same limitations as the planar MEC. A beautiful
curve primitive due to Mehlum [60] generalizes the Euler spiral to space curves. In addition toits excellent properties as a spline primitive, it has a beautiful geometric construction in terms ofrolling up an Euler spiral on the surface of the sphere.
Chapter 12 concludes and summarizes the main results of this thesis. The spline curves presented
in the following chapters are practical and well-suited to design tasks such as fonts. Theoreticalanalysis characterizes large regions of the space of all possible splines, suggesting that the idealinterpolating spline is, in fact, drawn from an easy-to-understand family. In addition, the discovery
that all extensional, two-parameter curves are derived from a characteristic generating curve unifies
previously disparate areas of the literature, in particular the Euler spiral spline and variational
techniques. It is hoped that these beautiful curves will become more popular in the future, nowthat their properties are better known, and their implementation need not be much more complex
or computationally expensive than various crude polynomial approximations have been in the past.

5

Chapter 2

Properties of splines
The space of all possible interpolating splines is infinite, and there is no consensus in the
spline literature of which is “best.” However, a number of properties have emerged as good criteria
for evaluating interpolating splines. Foremost among these is fairness, which is a technical term
for smoothness. There is no agreed-upon formal definition for fairness, but it is associated withsmoothly varying curvature.
Although fairness is important, other properties matter too. For example, some splines (including those modeling the mechanical spline) do not exist for arbitrary sequences of points. This
chapter explores fairness deeply, and lists other desirable properties of splines. In subsequent chapters, the evaluation of splines will be based on the properties described below.

2.1

Fairness

In almost all applications, the most important property of a spline is fairness, also known assmoothness. Fairness is ultimately a property of the way humans perceive curves, not an abstract
mathematical property. As in Justice Potter Stewart’s definition of pornography, “I know it when
I see it.”1
Even so, there are a number of metrics that correlate reasonably well with perception of fairness.
Knowledge of the human visual system and empirical evidence show that fairness is intimately
tied to curvature, which of course is an objective mathematical property. Curves with regions of
extremely high curvature are clearly not fair. In addition, while more subtle, discontinuities incurvature and regions with large curvature variation are also perceived as lack of fairness.
Curvature is central in many other application domains as well. There is compelling evidence
that the human visual system perceives curves in terms of curvature features, motivating a substantial body of literature on shape completion, or inferring missing segments of curves when shapes
1
Justice Potter Stewart, concurring opinion in Jacobellis v. Ohio 378 U.S. 184 (1964), regarding possible obscenity
in The Lovers.

6

are occluded in the visual field. A recent example, advocating the Euler spiral spline, is Kimia etal. [44].

2.1.1

Empirical testing of fairness

Bending energy has been proposed as a fairness metric, and minimizing that metric is thedefinition of the MEC spline (which is discussed in considerable detail throughout this thesis).
However, does it accurately measure human perception of smoothness?

0.5

1

1.5

2

2.5

3

3.5

4

Figure 2.1. Test instrument for fairness user study.

To address this question, I did a very simple user study, asking respondents to choose the fairest
curve from a range of possibilities, shown in Figure 2.1.2 The individual curves are all generated
with the formula κ = csα , with α the varying parameter and c chosen to make the curvature
2

The discussion thread is on the web at http://typophile.com/node/52821

7

continuous at the endpoints. These are part of a family of so-called “log-aesthetic curves” [63],discussed in more detail in Section 4.8.1.
18

survey results
MEC energy

16

14

number of votes

12

10

8

6

4

2

0
0

0.5

1

1.5

2
exponent

2.5

3

3.5

4

Figure 2.2. Survey results for aesthetic curve preference.

Among the curves in this family, minimal bending energy is obtained with an exponent α near
0.9; thus the Euler spiral, with the exponent α = 1, also yields very low bending energy. The on-line
study then asked respondents to select the “fairest” curve from the choices shown in Figure 2.1. The
results are plotted in Figure 2.2, along with the MEC bending energy. The mean preferred exponentis 1.78; this is a curve for which the bending energy is 6% higher than the Euler spiral. These results
confirm the belief that the concept of fairness is rather fuzzy, and that it varies significantly among
different designers and graphic artists. Nonetheless, the preference for an exponent in the range 1.5
to 2 is marked. In particular, this study establishes that MEC energy does not in fact accurately
predict the perception of fairness.

2.2

Continuity

A concrete property related to fairness is the degree of continuity, essentially the number ofhigher derivatives that exist. A parametrized polygon, for example, has a first derivative everywhere,
but that first derivative is not continuous. Thus, a polygon has only G0 -continuity (the property of
being a connected curve). All splines worthy of the name have at least G1 -continuity, meaning that8the second derivative exists everywhere, but relatively fewer have G2 -continuity, corresponding
to continuous curvature. It is worth remembering, though, that wild but continuous swings incurvature are no fairer in practice than tiny jumps at the control points.
Intuitively, curvature is the position of the steering wheel when driving a car along the curve.
Therefore, G2 continuity is equivalent to the lack of jerks of the steering wheel.

Figure 2.3. Continuity 6= fairness.

Figure 2.3 illustrates clearly that continuity is not the same as fairness. The spline on the left
is an Euler spiral spline, with G2 continuity, and the spline on the right is a circle spline, with G3 .
It seems clear that the one on the left, the G2 spline, is fairer, because it exhibits less variation incurvature.

2.3

Roundness

One specific aspect of fairness is roundness. A round spline always yields a circular arc when
the control points are co-circular. It seems intuitively obvious that a circular arc is the fairest (orsmoothest) curve.
It makes sense for circular arcs to be the result of splines, as the circle is often found in natureand is also essential in mechanics. However, many splines merely produce an approximation, somebetter than others. In particular, any spline based on polynomials cannot be exactly round.
Roundness is not an all-or-nothing property. Even the splines that do not generate perfectcircular arcs produce curves that are fairly close. Exactly how close can be quantified, as we will
see in the discussions of individual splines.
Splines other than polynomials can also fail the roundness property. For example, the curve
minimizing the bending energy of a thin strip (the Minimum Energy Curve) fails to describe aperfect circular arc. The deviation from roundness, interpolating three control points placed on an
equilateral triangle, is shown in Figure 2.4, and will be discussed in more detail in Section 4.2. In
this figure, the circle is on the inside, and the MEC describes a curve with somewhat larger total
arc length.
9

Figure 2.4. MEC roundness failure.

2.4

Extensionality

An extensional interpolating spline is one in which adding a new point to the control polygon,
lying on the curve generated from the original control polygon, does not change the resulting curve.
Unlike fairness, extensionality is not so much a property of the curve segments used to generate
the spline as of the effect of changes to the control polygon (adding and deleting points) on the
resulting curve.
Extensionality is closely related to the concept of choosing an optimal curve with respect tosome fairness metric. If adding an on-curve point changes the curve, either the new or the old curve
must not have been optimal. Therefore, splines defined in terms of minimizing a functional tend
to be extensional. Extensionality is a powerful filter for weeding out crude approximations, but it
does admit a variety of interesting splines in addition to those defined in terms of optimization.

Figure 2.5. Circle spline extensionality failure.

10Figure 2.5 shows a spline with particularly poor extensionality, the circle spline (described in
more detail in Section 4.10). Adding the point in the center changes the tangent radically at thatpoint, introducing new inflection points, and clearly degrading the overall quality. In this case,the problem arises because the tangents are derived from a local window; when three consecutive
control points are co-linear, that line determines the tangent of the middle one even when it is not
ideal given the global shape.
Roundness is a very weak extensionality property; for the special case of control points on acircle, adding an additional point on the circle does not change the resulting curve. However, eventhough these two properties are related, there are round splines that are not extensional (Ikarus[42]), and extensional splines that are not round (MEC).

2.5

Existence and uniquenessOne highly desirable property of a spline is that all control polygons yield interpolating curves.
Especially in interactive applications, where control polygons may well pass through ill-defined ordegenerate states on the way, it is inconvenient at best for there to be no solution.

Figure 2.6. Nonexistence of solution to MEC.

Figure 2.6 shows a configuration in which the MEC spline has no viable solution. In thisconfiguration, the tangent angles at the endpoints are constrained to be at greater than a rightangle with respect to the chord, but the curve is otherwise unconstrained. In particular, the lengthof the curve is allowed to assume the value that minimizes overall bending energy. Intuition suggests
that, for any given functional to optimize, some curve is better than none at all, so the lack of asolution is surprising.
Better intuition comes from considering the MEC as a model of a mechanical spline. Bending
energy tends to decrease as the curve gets longer. In particular, the longer curves in Figure 2.6
have a lower bending energy than the shorter ones. Thus, the mechanical springiness applies a force
to pull in more length. In a real spline, at some point this force becomes so small that it is easily

11defeated by the friction at the control point, but in the mathematical idealization, it would continue
indefinitely until the length became infinite and the bending energy zero.
Ironically, even though intuition would suggest that some curve is better than no curve atall, the splines with the worst existence properties are those defined in terms of minimizing afairness functional. Apparently, fairness is not the only measure of quality. In general, splines
defined in terms of meeting continuity constraints at the control points tend to have better existence
properties. For example, an interpolating Euler spiral spline always exists, as discussed in Section
4.6.4.

Figure 2.7. Uniqueness failure in Euler spiral spline.

Another desirable property is uniqueness, the principle that only one curve exists that satisfies
the definition of the spline, for a fixed input. For splines that are defined by optimization, theconcern is that there may be multiple local minima, where a local minimum is a solution that is
locally optimum under smooth deformation of the curve while it still interpolates the control points.
Correspondingly, for splines that are defined in terms of continuity constraints on parametrized
curve segments, there may be multiple assignments of parameters that satisfy the constraints andstill interpolate the points. Figure 2.7 shows an example, based on the Euler spiral spline. Bothcurves interpolate the same five input points, in the same order, both are composed of Euler
spiral segments between each pair of adjacent control points, and both satisfy the G2 -continuity
constraints.
It is possible to fix the uniqueness problems of a spline by adding an additional (global) criterion,
to prefer one solution to the other. Then, the revised spline can be defined as one that bothsatisfies the local constraints and also optimizes the global criterion. It may be difficult to properly
implement a spline revised in this way, as search for a global optimum is often expensive, but atleast such a spline is well defined.
As will be discussed in Section 4.7, Hobby used the criteria of existence and uniqueness tomotivate the design of his spline. Using “mock curvature” guarantees existence of a unique solution,even though it has less desirable fairness properties.

12

2.6

Locality

When a control point moves, how much of the resulting curve changes? A small section in theneighborhood of the control point, perhaps, or does the entire spline change? This is the propertyof locality.
Some splines have finite support, a particularly strong form of locality. In these splines, moving
a control point changes the curve only for a section bracketed by a finite number of control points
on either side of the one moving. For example, in circle-splines this section consists of four segments
bracketed by stepping out two control points in either direction.
However, finite support conflicts directly with extensionality, so in general we will measure
locality in terms of how quickly the effect of moving a control point dies out as we move away fromit. In splines with good locality, this effect falls off exponentially in the distance from the control
point (measured by counting control points). Then, locality can be quantified√ by the ratio of effectfrom one control point to the next. Section 4.8.1 shows, this ratio is 2 + 3 ≈ 3.73 for the Euler
spiral spline (as well as the MEC).

Figure 2.8. G2 spline has better locality than G4 .
Figure 2.8 shows two splines with different locality. The one on the left is an Euler spiral splinewith G2 -continuity, and the one on the right is its G4 counterpart. In general, there is a tradeoff
between the degree of continuity and the locality, so the latter spline has more wiggling farther
from the perturbed point.

2.7

Transformational invariance

Another basic property is invariance under various transformations and symmetries. All functions worthy of being considered splines are invariant under rigid body transformations (rotationand translation) as well as uniform scaling. If the control polygon undergoes a rigid body transformation, then the resulting curve undergoes exactly the same transformation.
Further, all of the splines we consider in this thesis are invariant under mirror symmetry andfront-to-back reversal of the order of the points in the control polygon. These are very common
properties, but it is worth noting that quaternion splines computed on the surface of the spherelack front-to-back symmetry.
A sometimes desirable property of splines is affine invariance, meaning that the shape of the
spline is preserved through affine transformation of the control points. This property is observed

13

in many polynomial-based splines, but it conflicts with the roundness property. Consider that any
round spline fits a circular arc to an input consisting of three control points. Yet, affine transformation of the resulting curve yields an ellipse, which in general would not match the spline resulting
from the three transformed control points – which would also need to be a circular arc due to theroundness property. For the types of designs considered in this thesis, roundness is more important
than affine invariance, so the latter will not be considered in detail.

2.8

Counting parameters

In an arbitrary spline, the family of curve segments between any two control points is potentiallydrawn from an infinite-dimensional space. In practice, most splines use curve segments chosenfrom a finite-dimensional manifold. Generally, there is a vector of real parameters that uniquelydetermines the shape of the curve between any two endpoints. In counting these parameters, wehold the endpoints fixed, and apply the rotation, scaling, and translation to make the curve coincidewith these endpoints.
Thus, the degenerate spline consisting of straight lines is effectively a zero-parameter spline.
A fair number of interesting splines fall into the category of two-parameter splines, and these are
discussed further in Chapter 4. Well constructed splines in this class can achieve G2 continuity.

2.9

Monotone curvature

The human visual system is particularly sensitive to minima and maxima of curvature. This
property suggests that all curvature extrema should coincide with the given control points, andthat the spline segments in between should exhibit monotone curvature. A fair amount of literature is concerned with constraining existing curves (such as cubic polynomials) to have monotone
curvature. One such citation is “Measures of Fairness for Curves and Surfaces” [75, p. 93], which
also discusses fairness metrics in considerable detail. The Euler spiral, in which curvature varies
linearly with arc length, obviously has this property, as do other related curves. However, the MEC
does not have this property; for symmetrical end conditions it has a tendency to increase arc lengthand to add a curvature maximum in the middle of the segment. The MVC also in general does notplace curvature extrema to coincide with the control points.
Evidence that the visual system is sensitive to minima and maxima of curvature goes all the way
back to Attneave’s seminal experiments in the mid 1950’s [5]. Attneave presented test subjects withsomewhat random, blobby drawings and asked them to annotate which points along the outline
were the most salient. An example is shown in Figure 2.9. The “porcupine quills” emanating fromthe shape represent the relative frequency that subjects identified the points. These strongly tend
to coincide with curvature extrema. As such, in interactive curve design it is reasonable to expect
designers to place these points first, then refine the curves with more control points only as needed.
A good spline should accommodate such a designer by preserving the curvature extrema as marked.

14

INFORMATIONAL ASPECTS OF VISUAL PERCEPTION

185

which the dots represented. A good
sample of the results is shown in Fig. 2:
radial bars indicate the relative frequency with which dots were placed oneach of the segments into which the contour was divided for scoring purposes.
It is clear that Ss show a great deal ofagreement in their abstractions of points
best representing the shape, and mostof these points are taken from regions
where the contour is most different froma straight line. This conclusion is verified by detailed comparisons of dot frequencies with measured curvatures onboth the figure shown and others.
Fie. 2. Subjects attempted to approximate
Common objects may be represented
the dosed figure shown above with a patternwith great economy, and fairly strikingof2.9.
10 dots.
Radiatingcurvaturebars indicateexperiment
the relaFiguretive
Attneave’s[5]. by copying the points at which
fidelity,frequency with which various portions of
the outline were represented by dots chosen. their contours change direction maximally, and then connecting these points
Evidence from other and entirely dif- appropriately with a straightedge. Fig2.10 Summary
ferent situations supports both of these ure 3 was drawn by applying this techinferences. The concentration of infor- nique, as mechanically as possible, to asleeping cat. The informational
mationin contours is illustrated
the real
The most important property
for interpolating
splines is byfairness.
Yet, othersuchcontentof aproperties
drawing likethisasmay beremarkably similar appearance of obconsideredtoconsistof
two comporobustness and locality are alsojects
important.
alike in contour and different
nents:
one
describing
thepositions ofotherwise.
Thea"same"for exSeveral previous authors haveproposed
set of triangle,
properties
desirablein splines.
Knuth’senu- whichthepoints,theotherindicating
ample, may be either white on black orpointsareconnected
with
which
others.
meration of such properties, forgreen
example,comprisedinvariance,
symmetry,
extensionality,
locality,on white. Even more impressive
The
first
of
these
components
will
alsmoothness, and roundness [47,ispp.
the39–42].
familiar fact that an artist's
mostalwayscontainmore
information
sketch, in which lines are substituted
than thesecond, butlist
its ofexact
share will
Similarly, Henry Moreton inforhissharp
Ph.D.
thesis[64, ch.may2] gives
exhaustive
propcolor
gradients,
consti-a rather
erties, including most of the ones
above.
tutedescribeda readily identifiable
representation depend upon the precision with which
positions are designated, and will furof a person or thing.
ther vary from object to object.
An experiment relevant to the second
Let us now return to the hypotheticalprinciple, i.e., that information is fursubject
whom we left between the corner
ther concentrated at points where acontour changes direction most rapidly,
may be summarized briefly.8 Eighty 5s
were instructed to draw, for each of 16
outline shapes, a pattern of 10 dots
which would resemble the shape as
closely as possible, and then to indicate
on the original outline the exact places
3

This study has been previously publishedonly in the form of a mimeographed note:"The Relative Importance of Parts of a Contour," Research Note P&MS Sl-8, Human Resources Research Center, November 1951.

15

FIG. 3. Drawing made by abstracting 38
points of maximum curvature from the contours of a sleeping cat, and connecting thesepoints appropriately with a straightedge.

Chapter 3

The Elastica
The word spline originally refers to a thin strip of wood used to smoothly interpolate control
points. A full understanding of splines requires an examination of the theory of such elastic strips.
This theory, that of the elastica, has a long and rich history, much of which will be explored in thenext chapter, but the literature is missing a concise, accessible treatment of the topic. This chapter
fills that gap.
One difficulty of the literature is that there are many different ways to approach the elastica:
as an equilibrium of moments in a static physical system; as an equilibrium of forces in the same
physical system (analyzed in terms of a finite element approximation composed of rigid struts andpivot springs); as the curve that minimizes the total bending energy of that system; as the solution
to an abstract differential equation; and through the analogue with another physical system, thependulum, which happens to share the same differential equation.
Each approach sheds light on a different facet of the elastica. For instance, treating it as astatics problem is the easiest way to understand the effect of different constraints operating on
the curve, and of course the physics interpretation is the only way to think about the differences
between a real elastic lamina and the mathematical ideal – especially, the effect of friction at theconstraint points. The variational formulation is most revealing when examining the minimizationof some other metric than total bending energy. And the pendulum analogy is helpful for developing
an intuitive understanding of the underlying differential equation, particularly the fact that it isperiodic.
My approach is to pull together all of these approaches, letting each tell its own story. For
a given problem involving the elastica, any one of these approaches may offer the simplest, mostintuitive solution. Yet, the underlying mathematics is the same. In the equilibrium of moments
formulation, the equation takes on a remarkably simple form: κ = λy, where κ is the curvature, yis the Cartesian ordinate, with x the abscissa, and λ is a constant.

16

3.1

Statement of the elastica problem

The elastica models a thin strip of flexible material (traditionally wood, but spring metal willdo). In the spirit of idealizing the problem, “thin” means that the thickness is negligible comparedwith the amount of bending. The mathematical idealization further assumes that the centerline ofthis strip lies entirely on a plane, and that there is no twisting, stretching, or compression.
Thus, the curve traced by the strip is characterized by its curvature at each point along its
length. A convenient definition of curvature is the derivative of tangent angle θ with respect to arc
length s. Unless otherwise stated, in this chapter the θ0 notation will represent dθ/ds, so curvatureκ is written simply θ0 .
Bending such a strip induces potential energy, much like that of a spring. According to thetheory of elasticity, the bending energy is proportional to the square of the curvature. Thus, thetotal bending energy of a strip of length l is
Z lE[κ(s)] =κ(s)2 ds .

(3.1)

0

When completely unconstrained, the elastica will assume the shape of a straight line, in which
the curvature is everywhere zero, and thus the total bending energy is also zero. When constrained,
the bending energy will tend to the minimum possible under the constraints.

3.2

Variational study

The elastica equation 3.1 is particularly well suited to study by the calculus of variations[19], a technique for transforming problems specified in terms of infinite-dimensional functions into
tractable ordinary differential equations. The calculus of variations was founded at the end of the
17th century, and was developed by Johann and Jacob Bernoulli, Leonhard Euler, and others.
In particular, Equation 3.1 yields readily to the Euler–Lagrange equation.

3.2.1

The Euler–Lagrange equation

If there exists a minimum value of the functional
Z x1
v[y(x)] =F (x, y(x), y 0 (x)) dx ,

(3.2)

x0

subject to the n constraints, each of the form 1 ≤ i ≤ n,
Z x1
ci (y(x)) dx = ki ,x0then this minimum is characterized by a tuple of values λi such that

17

(3.3)

d ∂F X dci
∂F
−
+
λi
=0.
∂y
dx ∂y 0
dy

(3.4)

i

Therefore, the problem of finding a function minimizing the cost functional v in an infinitedimensional space of functions is reduced to that of finding a tuple of values minimizing thatfunctional. In general, such an approach yields easily to numerical techniques for both integrating
the differential equation 3.4 and finding a set of values λi minimizing the cost functional.

3.2.2

Application of the Euler–Lagrange equation to the elastica

There are several ways to apply the Euler–Lagrange equation to find curves satisfying the elastica problem as defined by Equation 3.1. The main choice is that of independent and dependent
variables (x and y in the formulation of Equation 3.2). At least three such pairs appear in the literature of the elastica: arc length and tangent angle, arc length and signed curvature, and Cartesian
coordinates.
Of these, the choice of arc length s as the independent variable and tangent angle θ as thedependent yields the simplest variational solution. The most general endpoint constraints are specification of curve length, endpoint positions, and endpoint tangents. Of these, the endpoint position
constraint requires Lagrange multipliers in the form of Equation 3.3. Assume without loss of generality that the starting point is (0, 0) and the endpoint is (x, y) (simple translation will allow
arbitrary starting point). Then, the endpoint constraint is expressed as the integral of the unit
speed vector of direction θ(s):Rlcos θ(s) ds = x ,R0l0 sin θ(s) ds = y .

(3.5)

In terms of arc length and tangent, these constraints are easily expressible in the form ofequation 3.3.
It’s now easy to see why representing the curve in terms of tangent angle as a function ofarc length yields most readily to the variational approach. In terms of curvature κ, the endpoint
constraint is a double integral, which does not fit neatly into a formulation of constraints as a single
integral of some function. Alternatively, Cartesian coordinates would make the endpoint constraints
even simpler, but the variational function becomes more complex, requiring second derivatives aswell as the first derivatives of the Euler–Lagrange equation. The first solution was indeed doneusing Cartesian coordinates, but this approach is definitely simpler.
Applying the variational technique to minimizing the bending energy of a strip, Equation 3.1,
we arrive at the solution,
θ00 + λ1 sin θ + λ2 cos θ = 0 .

(3.6)

Equation 3.6 is the simplest form of the solution, but other forms can be readily derived, andoften appear in the literature. Taking the derivative with respect to s,
θ000 + λ1 θ0 cos θ − λ2 θ0 sin θ = 0 .
18

(3.7)

Conversely, multiplying both sides of 3.6 by dθ/ds and then taking the integral,
1 0 2(θ ) − λ1 cos θ + λ2 sin θ + C = 0 .
2

(3.8)

Combining 3.7 and 3.8, we can eliminate the Lagrange multipliers, retaining the single constantof integration C.
1
θ000 + (θ0 )3 + Cθ0 = 0 .
2

(3.9)

Since κ = θ0 , we can reformulate this differential equation entirely in terms of curvature:1
κ00 + κ3 + Cκ = 0 .
2(3.10)

Another formulation in the literature (for example, Eq. 12 of Lee and Forsythe [51]) multiplies
this last equation by 1/κ, yielding
1 00 1 2
κ + κ +C =0.
κ
2

3.3

(3.11)

A family of solutions

The equations derived in the previous section for the shape of the elastica don’t have a single
solution. Rather, they describe a family of solutions, characterized by a single scalar parameter. Allsuch solutions are realizable by a physical elastica; none are purely an artifact of the mathematical
approach.
For mathematical convenience, assign κ = 1 at the midpoint of the curve, and use λ1 as theparameter in Equation 3.6, setting λ2 = 0. If we write λ1 = λ cos φ and λ2 = λ sin φ, then all otherinstances of the elastica can be transformed to meet these conditions by scaling and rotating theentire system so that φ = 0. Since the assumption of φ = 0 obviates the need to distinguish multiple
λ parameters, it is convenient to drop the subscript 1 and simply write λ, resulting in a simplified
equation.
θ00 + λ sin θ = 0 .

(3.12)

Figure 3.1 shows a representative selection from this family, for various values of λ.
The remainder of this chapter is devoted to exploring these solutions and their physical meaning,
but a few highlights are immediately apparent. All solutions have a mirror (even) symmetry at eachcurvature maximum. All solutions save λ = 0.25 are periodic. The solutions for values larger than
0.25 have inflections, and have odd symmetry around each inflection point. For values smaller than
0.25, there is no inflection; the curvature always remains the same sign. Note that the curve has amirror symmetry for all values of λ.
As λ approaches 0, the elastica tends to a circle. As λ approaches ∞, the elastica tends to asine wave. The solution at λ = 0.5 is special, and is known as the rectangular elastica. The tangents19

λ = 0.1

λ = 0.2

λ = 0.249

λ = 0.25

λ = 0.251

λ = 0.28

λ = 0.3027

λ = 0.35

λ = 0.4

λ = 0.5λ=1
λ=2Figure 3.1. The family of elastica solutions.

20

at the inflection points of this curve are perpendicular to the horizontal axis; the curve is tangent
to the rectangle formed by inflections and curvature extrema.
Figure 3.2 shows the curvature plots of the instances of the elastica family shown in Figure 3.1,
plotted with curvature on the vertical axis and arc length on the horizontal. As is usual, inflection
points on the curve correspond to the curvature crossing (or touching) the horizontal axis.

3.4

Equilibrium of moments

Analyzing the elastica using the theory of elasticity yields three physical quantities which must
be in equilibrium along its length: tension (longitudinal force), shearing force, and bending moment.
The tension and shearing forces are useful when working out the effect of constraints, but it ispossible to derive the equation for the shape of the elastica from moments alone. Indeed, computing
the moments is by far the easiest derivation of the elastica equation.
This section uses elementary principles of statics to carry out the derivation of the elastica fromequilibrium of moments. A classic textbook introducing these concepts is Statics, by Goodman andWarner [28].
Consider the elastica in Figure 3.3. Two forces of magnitude F and opposite directions push theends of the elastica together. Because these forces are equal in magnitude, opposite in direction,and have parallel lines of action, they form a couple [28, p. 65] (the line of action is defined as avector in the same direction as the force, emerging from the point where the force is applied).
At any point along the elastica, then, the resultant moment is simply the force of the couple
times the distance from the line of action: M = F y, where F is the magnitude of the force vector
F.
The simple theory of elasticity says that the relation between bending moment and curvature
is linear. Thus, collecting all the constants together into a single λ yields this remarkably simple
formulation of the shape of a general elastica,
κ = λy .

(3.13)

This equation is readily shown to be equivalent to the variational solution above. Take thederivative with respect to arc length s, and note that y 0 = sin θ. Then,κ0 = λ sin θ,(3.14)and this equation is equivalent to Equation 3.12, assuming without loss of generality that the lineof action is parallel to the horizontal axis.
All elastica solutions are equivalent to the system shown in Figure 3.3, but in non-inflectional
cases the line of action of the couple is virtual, in that it doesn’t touch the curve.

21

λ = 0.1

λ = 0.2

λ = 0.249

λ = 0.25

λ = 0.251

λ = 0.28

λ = 0.3027

λ = 0.35

λ = 0.4

λ = 0.5λ=1

λ=2Figure 3.2. Curvature plots for the family of elastica solutions.

22M = FyyF

F

Figure 3.3. Equilibrium of moments.

3.5

Pendulum analogy

The differential equation 3.12 for the shape of the elastica is mathematically equivalent tothat of the dynamics of a simple swinging pendulum. This kinetic analog is useful for developing
intuition about the solutions of the elastica equation. In particular, it’s easy to see that all solutions
are periodic, and that the family of solutions is characterized by a single parameter (modulo scaling,translation, and rotation of the system).
Consider, as shown in Figure 3.4, a weight of mass m at the end of a pendulum of length r,with the angle from vertical specified as a function of time θ(t).
Then the velocity of the mass is rθ0 (where, in this section, the 0 notation represents derivativewith respect to time), and its acceleration is rθ00 . The net force of gravity, acting with the constraintof the pendulum’s angle is mg sin θ, and thus we have
Fnet = ma = mrθ00 = mg sin θ .

(3.15)

With the substitution λ1 = −g/r, and the replacement of arc length s for time t, this equation
becomes equivalent to Equation 3.12, the equation of the elastica. Note that angle (as a functionof arc length) in the elastica corresponds to angle (as a function of time) in the pendulum, and theelastica’s curvature corresponds to the pendulum’s angular momentum.
The swinging of a pendulum is perhaps the best-known example of a periodic system. Further,
it is intuitively easy to grasp the parameter space of the system. Transformations of scaling (varying
23

θ

r
F = mg

Fnet = mg sin θ

Figure 3.4. An oscillating pendulum.

the pendulum length) and translation (assigning different phases of the pendulum swing to t = 0)
leave the solutions essentially unchanged. Only one parameter, how high the pendulum swings,
changes the fundamental nature of the solution.
The solutions of the motion of the pendulum, as do the solutions of the elastica equation 3.12,
form a family characterized by a single scalar parameter, once the trivial transforms of rotationand scaling are factored out. Without loss of generality, let t = 0 at the bottom of the swing (i.e.
maximum velocity) and let the pendulum have unit length. The remaining parameter is, then,the ratio of the total energy of the system (which remains unchanged over time) to the potential
energy of the pendulum at the top of the swing (the maximum possible), in both cases counting
the potential energy at the bottom of the swing as zero. Thus, the total system energy is also equal
to the kinetic energy at the bottom of the swing.
For mathematical convenience, define the parameter λ as one fourth the top-of-swing potential
energy divided by the total energy. This convention is equivalent to varying the force of gravity
while keeping the maximum kinetic energy fixed, which may be justified by noting that the zerogravity case (which would require infinite kinetic energy if we were assuming unit gravity) is far
more relevant in the context of splines (it corresponds, after all, to a circular arc, so any elasticabased spline meeting the criterion of roundness must certainly have this solution) than the opposite
of zero kinetic energy.
The potential energy at the top of the pendulum is 2mgr. The kinetic energy at the bottom of
the swing is 21 mv 2 , where, as above, v = rθ0 . Thus, according to the definition above,
λ=1 2mgrg=.
10 )2
4 2m(rθ
r(θ
0 )2

24

In other words, assuming (without loss of generality) unit pendulum length and unit velocity
the bottom of the swing, λ simply represents the force of gravity.
Figure 3.2, the curvature plots for solutions of the elastica equation, thus also shows the velocityof the pendulum as a function of time for various solutions of the motion of the pendulum. For
each value of λ, velocity (θ0 ) is plotted against time. The bottom of the swing (maximum velocity)is the center of the plot, and 25 units of time are plotted to either side.
Note that the period has a singularity at λ = 0.25. For values of λ greater than 0.25, thependulum swings back and forth periodically, velocity reaching zero at the highest point in theswing. For values less than 0.25, the pendulum winds round and round, slowing at the top of thecycle, but never actually reaching zero. At exactly 0.25, the kinetic energy at the bottom is equal
to the potential energy at the top, so the velocity is zero at the top.
In the kinetic analogy with the elastica, velocity of the pendulum corresponds to curvature.
Thus, values of λ above 41 have inflection points. Figure 3.5 shows the relationship graphically,
including the right angle made by the pendulum at the top of its swing for the λ = 0.5 rectangularelastica case.
Note that while λ = 0.3027 is a special value for the elastica (it loops over itself and occupies
finite space), there is no special physical meaning in the pendulum analogy.

3.6

Equilibrium of forces: Finite element approach

When an elastic strip is in a minimal energy configuration, it is in mechanical equilibrium,
meaning that at all points along its length, the sum of forces on that point is zero.
The standard derivation of the equations of mechanical equilibrium from the theory of elasticity[52] expresses relations between shear forces, tension (longitudinal) forces, and flexural couples. The
equations are straightforward, but the relevant physics requires appeal to results in the theory ofelasticity. This section, therefore, presents a simplified finite-element model with only bending andtension forces.
The finite element model is a chain consisting of alternating rigid struts and pivots with torsion
springs. Such a chain, with spring torque indicated schematically at each pivot, is shown in Figure
3.6.
Each strut generates a tension force at both ends, in the same direction as the strut’s length.
In Figure 3.7, this tension force is indicated by the quantity T . By convention, compression force
is denoted with a negative sign.
At each pivot is a spring generating forces tending to straighten out the chain, i.e. towards aturning angle of zero. In Figure 3.8, the turning angle is denoted by θ and the the force generated
by the moment (or spring torque) by M . At the endpoints, the force generated by the moment isnormal to the strut joining the center to the endpoint. According to Newton’s third law, the sumof all forces arising from the element must be zero, so there is a balancing force at the center point
equal to the vector sum of the forces at the endpoints, pointing in the opposite direction.
The spring at the pivot point is assumed to be purely elastic, which is to say that the moment(torque) is proportional to the turning angle. Using c for this constant of proportionality (which,

25

5
1/λ
λ = 0.2432
10

λ = 0.26

4
1/λ
3
2
1
0

λ = 0.5

4
1/λ
3
2
1
0

λ=2

4
1/λ
3
2
1
0

Figure 3.5. Examples of pendulum analogy.

26

Figure 3.6. Finite element approximation of elastica as chain of struts and pivots.
T

T

Figure 3.7. Forces on a single strut.

in a real physical model, would be Young’s modulus of elasticity times the second moment of thesection of the beam about the axis of bending),
Mi = c∆θi .

(3.16)As shown in Figure 3.9, when these primitive elements are linked together in a chain, there are
five force vectors acting on a single pivot point: the tension force from each adjoining strut (T0 and−T1 ), the moment forces from the previous and next pivots in the chain (M0 and M2 ), and thebalancing force −2M1 cos(∆θ1 /2).
In flatland, the equilibrium condition that these vector forces sum to zero induces two scalar
equations. When the axes are chosen to be parallel and perpendicular to one of the struts, theseequations have a particularly simple form. Arbitrarily choosing the strut labeled 1 in Figure 3.9,
the parallel component is−T1 + cos ∆θ1 T0 − sin ∆θ1 (M0 − M1 ) = 0 ,(3.17)and the perpendicular component issin ∆θ1 T0 + cos ∆θ1 (M0 − M1 ) + M2 − M1 = 0 .

(3.18)In the continuous limit of this finite element model, the individual turning angles at the pivots

2Mcos(∆θ/2)
M

M

∆θ

Figure 3.8. Forces on a single pivot.

27

2M1 cos(∆θ1/2)

∆θ0

0

2
1

T1

∆θ2

∆θ1
T0
M0
M2

Figure 3.9. Forces on a single pivot point in a chain.

∆θi tend to zero, justifying the usual linear approximations sin ∆θ ≈ ∆θ and cos ∆θ ≈ 1. The resultof using these approximations, and substituting the elastic pivot assumption of Equation 3.16 is−T1 + T0 − c∆θ1 (∆θ0 − ∆θ1 ) = 0 ,(3.19)∆θ1 T0 + c(∆θ0 + ∆θ2 − 2∆θ1 ) = 0 .

(3.20)In the limit as ∆θ approaches zero (as the number of links in the chain increases), these difference
equations become the differential equations
−

dκ
dT
+ cκ
=0,
ds
ds

(3.21)

d2 κ
=0.
ds2

(3.22)c
T = κ2 + C .
2(3.23)

κT + cEquation 3.21 is readily integrated, yielding

Substituting this into Equation 3.22, and dividing by the common constant c, yields a restatement of Equation 3.10.
1
κ00 + κ3 + Cκ = 0 .
2(3.24)

The fact that the constant c, related to the coefficient of elasticity, factors out has the physicalinterpretation that the elastica takes the same shape regardless of the value of this coefficient, aslong as it is consistent across the length of the elastica.
28

It is satisfying but not surprising that the variational and the mechanical equilibrium treatmentsof the elastica problem result in the same equations. The theory of conservative systems states thatat a stable equilibrium, the energy is at a local minimum. In addition to confirming the variational
derivation, the mechanical equilibrium equations prove a useful physical interpretation for theconstant C in terms of curvature and tension (longitudinal force). In particular, if C is zero, then
the tension vanishes at inflection points.

3.7

Solutions with length and endpoint constraints

Now that we have the mathematical solution to the elastica equation, including plots for avariety of values of the parameter λ, we turn our attention to particular solutions corresponding toreal elastic strips under given constraints. Can real elastic strips take on all values of λ, or is someof the solution space a mathematical artifact?
Without loss of generality, assume that the length of the strip is 1, and that the endpoints lie
on the horizontal axis (all other cases can be recovered through scaling, rotation, and translation).
Three parameters remain: the tangent angles at the endpoints with respect to the horizontal axis,
and the distance between the endpoints (chord length).
Figure 3.10 shows one linear section of this three dimensional parameter space: the endpoint
tangents are fixed at π/2 from horizontal, and the chord length varies from 0 to 1.
Several properties are apparent: λ varies smoothly from about 0.34 to a maximum of 0.5, thenfalls quickly to 0, where the solution is a precise half-circle, then increases asymptotically to 0.25
as the strip is pulled tighter and tighter towards its full length. For these endpoint angles, the λmaximum at 0.5 also corresponds to zero curvature (inflection points) at the endpoints.
Note that the above discussion only reflects the principal solution of the elastica with respect tothe constraints, i.e. the one with the minimum total arc length from the original curve. The principal
solution always has at most one inflection point. The periodic nature of the elastica implies that
additional solutions may also exist, with more than one inflection point. In general, these additional
solutions are not interesting for the purpose of generating spline curves, but they are valid solutions.
Figure 3.11 shows a number of these solutions for the end constraints θ0 = π/4, θ1 = π/4, λ = 1/2.
It is interesting that some of these solutions lack symmetry even though the constraints themselves
are symmetrical.

3.8

Elastica with general constraints

Working towards the application of an elastica to splines, this section considers the effect ofgeneral constraints on an elastica. A goal is to develop physical intuition, so the constraints are
expressed in terms of forces.
There are four basic types of constraints. Each constraint, can either rotate freely or fix thetangent angle of the elastica. Further, it can either allow the elastica to slide freely, or can fix therelative arc length of the elastica. A segment between two constraints that fix relative arc length has

29

λ
0.5

0.45

0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0

0

0.1

0.2

0.3
0.4
0.5
chord length / arc length

0.6 2/π0.7Figure 3.10. Values of λ for varying chord lengths.

30

0.8

0.9

1

Figure 3.11. Multiplicity of solutions for same tangent constraints.

fixed arc length, even if there are sliding constraints in between. In all these types of constraints,
the position of the constraint is fixed, and the elastica is constrained to go through that point.
A position constraint, freely rotating, and with the elastica free to slide through the point,
exerts only a normal force on the elastica. If the tangent angle is fixed, the constraint exerts anormal force on the elastica. Similarly, if the relative arc length is fixed, the constraint exerts alongitudinal force.

Result 1 A normal force acting on an elastica preserves G2 -continuity across the point where theforce is applied, but creates a discontinuity in the first derivative of curvature.

Result 2 A shear force (torque) acting on an elastica preserves G1 -continuity of the elastica but
creates a discontinuity in curvature.

Result 3 A segment of an elastica with zero longitudinal force at the endpoints takes the form ofa rectangular elastica. Nonzero longitudinal force induces values of λ other than 0.5.
Proofs of these results are not included here, as they are known in the literature of the elastica.
An accessible source for them, in the context of interpolating splines, is Lee and Forsythe [51].
31

F

F

Figure 3.12. Forces acting on rectangular elastica.

Considering a full cycle of the rectangular elastica, as shown in Figure 3.12, it is clear thatthe longitudinal forces are zero at the endpoints, because the curve is perpendicular to the force F
acting on it.

3.9

Angle constraints

The general constraint approach of the previous section yields a specific result relevant to using
the elastica as a spline.

Result 4 An elastica segment with an angle constraint at each endpoint, but allowed to slide freely
through those constraints, is a section cut from the rectangular elastica.
Equivalently, the minimum energy curve (of arbitrary length) with given endpoints and given
angles at the endpoints is a section cut from the rectangular elastica.
Figure 3.13 shows the full parameter space of angle-constrained minimum energy curves. For
sufficiently large angles, no such solution exists; the envelope for which the MEC exists is shown as
well. This envelope was determined by implementing a very robust numerical solver and measuring
the boundary at which the solver did not converge. As far as is known, there is no analytical formulafor the envelope’s shape. Only half of the complete parameter space is shown, to avoid clutter. The
envelope is symmetrical around the θleft = θright diagonal.
Similarly, Figure 3.14 shows the curvature plots of all angle-constrained MECs within theirparameter space. It is worth noting that, for small angles, the curvature is very nearly linear witharc length.
Note that the curvatures at the endpoints are fully determined by the angles at the endpoints.

32θright

π/2

π/2

θleft

−π/2

Figure 3.13. Angle-constrained minimum energy curves.

33θright

π/2

π/2

θleft

−π/2

Figure 3.14. Curvature plots, angle-constrained minimum energy curves.

34

3.10

Elastica as an interpolating spline

The classical formulation of an interpolating spline is the minimum energy curve passing througha sequence of control points. Equivalently, it is an elastica passing through a sequence of constraints,
each of which rotates freely and allows the elastica to slide freely through it.
Because each of these constraints exerts only a normal force, the curve is G2 continuous. Since
the elastica is unbent past either endpoint, the curvature at the endpoints is zero (this is also known
as the “natural” endpoint condition).
Each segment of the spline, bounded by two control points, is an instance of an angle-constrained
MEC as described in the previous section.
To find an interpolating spline, it is necessary and sufficient to find a global solution assigning atangent angle to each point so that the resulting curvatures are continuous across each point (andzero at the endpoints).

3.11

SI-MECTo address some of the shortcomings of the MEC as a spline, Moreton and Séquin proposed
the Scale Invariant Minimum Energy Curve (SI-MEC) [65]. It is defined as the curve minimizing
the functional
Z l
SIMEC[κ(s)] = lκ(s)2 ds .

(3.25)

0

The SI-MEC derives its name from the fact that the functional is invariant with respect toscaling of the curve. The MEC, by contrast, is inversely proportional to the scale factor. However,this scaling property has little, if any, physical meaning, as the resulting shape is invariant to scalingof the entire physical system in either case.
The SI-MEC functional of Equation 3.25 is equivalent to the bending energy of a lamina withlength normalized to unity. Thus, a SI-MEC with angle constraints at the endpoints is equivalent
to this physical system: the endpoint constraints slide on a rail, but have their angles fixed. The
lamina is of fixed length and is not allowed to slide through the endpoint constraints.
The SI-MEC will, for any given constraints, be shorter than the corresponding MEC. Otherwise,
the SI-MEC would have even lower MEC energy than the MEC curve, which is impossible bydefinition.
Any elastica meeting the given angle constraints can be evaluated in its MEC energy and its
SI-MEC energy. In short, the MEC energy is the bending energy when the chord length is held
fixed and the arc length varies, and the SI-MEC energy is the bending energy when the arc lengthis held fixed.
Consider first the case of symmetrical angle constraints. Figure 3.15 shows the MEC and SIMEC energy for all elastica meeting the angle constraints (π/2, π/2), of chord length 0 to 1 (assuming unit arc length). The MEC energy is minimized, as expected, for λ = 0.5, the rectangularelastica. The SI-MEC energy is minimized for the circle, λ = 0. Figure 3.16 shows a similar plot

35

λ
0.5

0.45
λ

0.4

0.35

0.3

0.25

0.2

SIMEC energy

0.15

0.1
MEC energy
0.05

0

0

0.1

0.2

0.3
0.4
0.5
chord length / arc length

0.6 2/π0.7

0.8Figure 3.15. MEC and SI-MEC energies for varying lengths (π/2).

36

0.9

1

λ
1

0.9
λ
0.8

0.7

0.6

0.5

0.4

0.3

SIMEC energy

0.2

0.1
MEC energy
0

0

0.1

0.2

0.3
0.4
0.5
chord length / arc length

0.6

0.7

0.8

Figure 3.16. MEC and SI-MEC energies for varying lengths (π/3).

37

0.9

1

energy (relative units)

SIMEC energy

MEC energy

chord length

Figure 3.17. SI-MEC and MEC curve energy for (π/2, −π/2) angle constraints.

for angle constraints (π/3, π/3), showing that as the angles become less extreme, the difference inchord length between MEC and SI-MEC becomes less pronounced.
Figure 3.17 plots the MEC and SI-MEC energy as a function of the chord length / arc length
ratio, given angle constraints of (π/2, −π/2).
Along with each representative elastica curve is shown the continuation of the elastica curve(with dashed lines), and the line of action (of which the length is proportional to the force of action
applied by the constraints).
Note that the line of action for the curve minimizing the SI-MEC energy is vertical. In fact,
the line of action is always perpendicular to the chord for a SI-MEC with angle constraints at theendpoints, as can be deduced from the equilibrium of forces in the physical model. If the action had
nonzero component along the chord, the constraints would slide closer together or farther apartalong the rail.
Several interesting properties follow immediately from the vertical line of action. Based on
equilibrium of moments, the curvature varies linearly along the chord length. For symmetricalangle constraints, the curvature is also symmetrical, and thus must be equal at the two endpoints.
Thus, the curvature is constant, and the curve is a circle.
For the angle constraints (π/2, −π/2), the SI-MEC is a segment of the lemnoid elastica. Thisis also true for π/2 and the lemnoid elastica angle (≈ .71 radians).

38

Figure 3.18. A real spline.

Further, all such angle-constrained SI-MEC curves have λ ≤ 0.5, for the λ > 0.5 cases are
single-valued in x for any given y (where the line of action is drawn as vertical), and of course both
endpoints must lie on the x axis.
Because curvature is linear along the chord length, for angle constraints with absolute value≤ π/2, curvature is monotonic with arc length. Contrast with the MEC, which has a curvaturepeak for all symmetric and nearly-symmetric angle constraints, even as the angles approach zero.
As described in Section 2.9, monotonic curvature is a desirable property for curve segments used
in interpolating splines.
The domain where the angle-constrained SI-MEC has a solution is about double the size of the
MEC’s. In the symmetrical case, it reaches a singularity for the full circle (π, π) (exactly doublethat of the MEC), and in the anti-symmetrical case the singularity is a full figure-eight around thelemnoid elastica, approximately (5.57, −5.57), as opposed to the MEC’s (2.22, −2.22).
The SI-MEC has some problems when generalizing from angle constraints at the endpoints tomultiple position constraints. In particular, it does not form an extensional spline.

3.12

Real splines

As mentioned, the MEC is a mathematical idealization of a real wooden spline, the kind shownin Figure 3.18 that has been used in shipbuilding for hundreds of years.

39

Figure 3.19. High-friction spline constraints on a real spline.

For small angles, the two agree closely. But for large angles, how accurately does the MEC modela real spline? The fact that the MEC does not exist for certain angle constraints is particularly
suspicious.
The key to this question is, I believe, tension forces and friction. In a real spline, small tension
forces are countered by friction at the control points.
Indeed, shipbuilders use traditional spline weights in two different ways. In one such way, thespline weight pushes against the side of the spline, and the spline is free to slide with fairly low
friction.
Another common way to use spline weights is to place the point of the pin on the spline curve,
as shown in Figure 3.19. In this way, the weights effectively constrain the arc length as well as the
position, leaving only rotation entirely free.
It is understood that applying such constraints indiscriminately yields curves of poor fairness(applying tension at the constraints yields curvature discontinuity). Thus, designers use an iterativetechnique to improve fairing. Once the curve is laid out roughly, each weight is lifted, the spline
makes small, fairly localized adjustments when it is free from the weight, and then the weight isreplaced, in turn. Of course, if the resulting curve diverges too much from the desired shape, thespline is pushed when the weight is replaced, or a new weight is added.
Iterative local adjustment is reminiscent of successive over-relaxation, a common numerical
technique for solving global systems of equations. Note that, as a weight is lifted, the total arc
length bounded by the two neighboring weights remains constant. Thus, this process does notinevitably converge to the MEC solution, which has longer arc length than a spline under sometension. Solutions similar to the SI-MEC, or to circular arcs, remain stable.
Friction, then, seems an appealing technique for improving the robustness of the pure MEC,
and for making it conform to circles more closely. Shipbuilders have built up a tradition for using
40

Figure 3.20. Bricks as high friction constraints on a real spline.

splines effectively. For the spline, they often prefer poplar or white pine for the material, about aquarter inch thick (but thinner for more detailed work involving smaller radii, or thicker for large
boat lines). For the weights, they usually prefer a pound or two so it holds the spline firmly, arounded shape resembling a duck or whale for easier handling, and a configuration of the pin sothat the weight itself is separated from the spline. Even so, the traditional weight has no different
effect on the spline than would a brick (and simple bricks are pressed into service by do-it-yourself
boatbuilders, as well, as in Figure 3.20). Aside from shape of the weight, though, all these factors
affect the degree to which the real spline models a MEC, especially the friction at the constraints.
A number of researchers have been tempted by the possibility of using tension to make elasticabased splines better than the pure MEC. However, there is no clear way to accomplish this. Real
splines use spline weight and wood thickness to set the free parameters, but on a computer theseparameters would either need to be set explicitly by the user, or in some other way – at worst, a
computer representation of physical quantity such as “one pound weight” or “quarter inch thickwhite pine” are “voodoo constants,” which are well known for limiting the scope of usefulness of asystem.
All these problems are solved, though, when the elastica is left behind, and the Euler spiral isused in its stead. As described in Chapters 5 and 6, dedicated to the history of these curves, theirfates are closely intertwined. In fact, the complete mathematical characterization of both curves
can be traced to the same 1744 treatise by Leonhard Euler. But the details will wait for Chapter
5. The next chapter explores the properties of G2 -continuous interpolating splines in general, andquantifies exactly why the Euler spiral is optimal as a curve primitive for interpolation.

41

Chapter 4

Two-parameter splines
The broad goal of this work is to find the best spline, or at least provide a menu of best choices,each suited to a specific requirement. In attempting to characterize the space of all splines, including
making a taxonomy of existing splines in the literature, a remarkable number share this property:
A two-parameter spline is defined as one for which each curve segment between
two knots is uniquely determined by the two angles between tangent and chord at theendpoints of the segment.
Examples of two-parameter splines include the Minimum Energy Curve, the Euler spiral spline,
Mehlum’s kurgla 1, Hobby’s spline, Ikarus, and Séquin’s circle spline. Note that the count ofparameters excludes the scaling, rotation, and translation required to get the curve segment tocoincide with the two endpoints – those additional parameters are uniquely determined by theconstraints that the curve segment fits the endpoints, providing no additional flexibility for the
shape of the curve.
Any given two-parameter spline can be characterized by two properties: the family of curves
parametrized by the endpoint tangent angles, and the algorithm used to determine the tangent
angles. For example, the Euler spiral spline uses a curve family with linearly varying curvature,and uses a global solver to ensure G2 -continuity to determine the tangents.
Of the two-parameter splines listed above, the MEC and the Euler spiral are extensional. In
the case of the MEC, this property follows directly from the variational statement – if adding a
new on-curve point would produce a new curve with a lower cost functional, that curve would havesatisfied the old constraints as well. In addition, both splines are based on cutting segments froma fixed generating curve, the rectangular elastica or the Euler spiral.
These properties are not simply coincidental. In fact, all two-parameter, extensional splines
correspond to a single generating curve, and any generating curve can be used as the basis of a
two-parameter extensional spline. The two families, then, are essentially equivalent.

42

4.1

Extensional two-parameter splines have a generating curve

A central result of this thesis is that for any extensional two-parameter spline there exists
a generator curve so that all curve segments between adjacent control points can be cut fromthat generator curve, subject to scaling, rotation, translation, and mirror image transformations tomake the endpoints of the cut segments align with the endpoints in the spline. Thus, any splinein this family can be defined simply in terms of the generator curve. In some ways, this definition
circles back to one of the earliest formulations of nonlinear splines, “After looking carefully into therelevant equations, one realizes that the schemes of [53] and [23] do not approximate mechanicalsplines more closely than other curves – nor does it seem particularly desirable to have them doso. For example, they approximate equally well to Hermite interpolation by segments of Euler’s
spirals, joined together with continuous curvature.” [10, p. 171]. The Euler spiral is but one of afamily of desirable curves for constructing extensional interpolating splines.
Theorem 1 In an extensional, G2 -continuous, two-parameter spline, there exists a curve such thatfor any spline segment (parametrized by angles θ0 and θ1 ) there exist two points on the curve (s0and s1 ) such that the segment of the curve, when transformed by rotation, scaling and translation
so that the endpoints coincide, also coincides along the length of the segment.

κ, κ′, κ′′
s0

θ0

θ1

s1
s2

Figure 4.1. Incremental construction of the generator curve of a two-parameter spline.
A typical section of curve is shown in Figure 4.1. For some given endpoint angles, the generated
curve segment is the section from s0 to s2 . Because of the extensionality, for any s1 along this curve,
the subdivided curve from s0 to s1 (which is the result of applying θ0 and θ1 endpoint constraints)
must be consistent with the original segment. Modulo scaling, the curvature κ and its derivatives
κ0 and κ00 must of course also match between these two segments.
The proof of this theorem revolves around the quantity κ0 /κ2 , which is also the limit of 6(θ1 −θ0 )/(θ0 + θ1 )2 as θ0 and θ1 both approach zero, i.e. as the length of the curve segment becomes
infinitesimal (this limit follows from Equation 8.31, which is derived in detail in Section 8.2). Thisquantity represents the amount of curvature variation for a segment of a set amount of curvature –it is invariant to similarity transformations including uniform scaling. As such, for a given generator
curve, it uniquely identifies a point on the curve (modulo periodicity, if the curve is periodic as
opposed to monotone in curvature).
We propose that for any two-parameter spline that is also extensional, the second derivative ofcurvature, suitably normalized for scale invariance, κ00 /κ3 , is uniquely determined given a particular43

value of κ0 /κ2 . If there are two curve segments with different values for this second derivative, thenit is possible to cut subsegments from each with identical θ0 and θ1 values. However, the value ofκ00 /κ3 remains invariant even after such a cutting, so there are two different values for this. Yet,
the two parameter property states that there exists only a single curve for a given θ0 , θ1 pair – acontradiction!
Thus, there is a process for constructing the generator curve. Start with arbitrary κ and κ0values sampled from some point on the curve, then continually derive κ00 based on the lookup
process suggested by the above relation: find a curve segment that has the correct value for κ0 /κ2 ,then sample κ00 /κ3 at that point. Continually integrate (in the sense of an ordinary differential
equation), and the result is a curve.
The final component of the argument is that for any given θ0 , θ1 pair, there exists a segment
cut from this curve such that it is equal to the spline segment. Determine κ0 /κ2 for one side (s0 )of the spline segment, and find that point on the generator curve. Continue extending the curve
until the other side is reached. Since at each point κ00 /κ3 must be equal between the two curves,the curves themselves must also be equal.
This argument depends on the assumption that the curve is a continuous function of the parameters θ0 and θ1 . We conjecture that this assumption follows from the assumption of G2 -continuity
and extensionality, but do not have a rigorous proof. It is also possible the continuity assumptions
could be weakened and the theorem would still be valid.
Both the two-parameter and extensionality conditions are necessary for the generator curve
property to hold, and there are counterexamples when either is not present. The Ikarus and Metafont
splines are two-parameter curves, but not extensional, and they do not have a single generator curve.
The Minimum Variation Curve (discussed in much more detail in Section 7.2) is extensional but
does not fit into a two-parameter space, and also does not have a single generator curve.
The Minimum Energy Curve and the Euler spiral spline have a number of features in common,
among them that all segments between two adjacent control points are fully determined by two
parameters. This count of parameters does not include the scaling, rotation, and translation needed
to make the two endpoints coincide – just the shape of the curve once the positions of the endpoints
are fixed.
In fact, the parameter space for each such segment is described by the tangent angles at endpoints. Since our definition of “interpolating spline” requires G1 -continuity, determining a single
tangent angle at each control point suffices. Thus, combining a two-parameter curve family with atechnique for determining tangent angles at the control points yields an interpolating spline.
Note that there are also several splines in the literature designed to approximate the MEC(quite accurately in the case of small angles), but do not fit into the two-parameter framework.
In general, cubic polynomial curves have four parameters. In a standard cubic spline, the curve
segments can be any cubic polynomial, so the parameter space is four dimensional. Parametrizingby chord length approximates the MEC more closely but does not reduce the dimensionality of theparameter space.
Similarly, the scale invariant MEC (dubbed SI-MEC in [65]) has segments that are piecewise
elastica, and there is an additional parameter global to the spline (and associated with tension
forces at the endpoints, in a physical model of the spline) that increases the dimension space ofindividual segments from the two of the pure MEC to the three admitted by the elastica.

44

While the previous chapter examined the elastica with constrained endpoints, a more general
problem is the minimum energy curve that interpolates a sequence of control points. More formally,
we seek, over the space of all curves that pass through the control points (xi , yi ) in order, the onethat minimizes the MEC energy functional, restated here:
Z l
EMEC [κ(s)] =

κ2 ds

(4.1)

0

This problem is the mathematical idealization of a mechanical spline through control points
fixed in position but free to rotate, and with an infinitely thin beam free to slide through themwith zero friction.
The earliest discussion of this more general problem is probably Birkhoff and de Boor in 1964.
Over the next few years, several authors published algorithms for computing these “nonlinear
splines.” Probably the most illuminating and accessible discussion of the early nonlinear spline
theory is the 1971 report of Forsythe and Lee [51]. A brief summary of their results follows:• Each segment between knots is a MEC with λ = 0.5.
• Curvature is continuous across knots (G2 -continuity).
• Curvature is zero at endpoints (in the open case).
• The spline is extensional.
For each segment between knots, the elastica with λ fixed has two remaining parameters,
corresponding to the tangent angles at the endpoints. Recall that Figure 3.13 shows instances overthe entire two-dimensional space in which such solutions exist, as well as the boundary of thatsolution space. In particular, in the case of symmetrical endpoint tangents, the curve is not a
circular arc, so the MEC spline does not have the roundness property. Forsythe and Lee found thisunappealing, and propose adding an arc length penalty term to the functional,
Z lE[κ(s)] =κ(s)2 + c ds .

(4.2)

0

As we now know, this change admits the solutions to the elastica equation other than λ = 0.5.
Forsythe and Lee derive the value for c so that the solution is a circular arc (corresponding to λ = 0),
and speculate “whether such an approach could be generalized is an open question.” Answering
their question is a primary focus of this chapter.

4.2

Properties of the MEC spline

The most important property the MEC spline is lacking is roundness. For three control points
placed in an equilateral triangle (shown in Figure 4.2, with the spline drawn in a thick line, and thecircle going through the control points thin), the containing radius for the spline is approximately
1.035 times the radius of the circle going through the control points. With four points, this radius

45

Figure 4.2. MEC roundness failure.

deviation decreases to about 1.009, and asymptotically converges as O(n4 ) with the number ofpoints.
Another unpleasant property of the MEC spline is the failure to converge to a stable solution.
This property is intimately related to the relatively small envelope enclosing the solution space –for symmetric solutions the maximal endpoint tangent for λ = 0.5 is π/2, and it is not much betterfor any of the non-symmetric solutions.

4.3

SI-MEC

Several researchers have tried to address these shortcomings in the MEC spline. One is the scale
invariant minimum energy curve (or SI-MEC) proposed by Moreton and Séquin [65]. This spline isdefined as the curve minimizing the SI-MEC functional
Z lE[κ(s)] = lκ(s)2 ds .

(4.3)

0

The functional is termed “scale invariant” because its value is invariant as the entire curve isscaled up. By contrast, the MEC functional is proportional to the reciprocal of the scaling constant.
The MEC functional can be seen as rewarding increases in the length of the curve, which accounts
for the roundness failure – the actual curve obtained is slightly longer than the circle, which accounts
for the lower MEC functional. In the SI-MEC, the bias toward longer curves is removed, and, infact, the spline is round when the control points are co-circular.
The SI-MEC minimizes the MEC functional for a given length (a curve with a lower MEC
energy at the same length would also obviously have a lower SI-MEC cost functional). Thus, anappealing physical interpretation is that of a mechanical spline with a bit of tension pulling at theendpoints, making the overall curve slightly shorter. Such a curve is generally not two-parameter46

in the formulation above, and, in fact, all values of the elastica are possible (thus, it can be seen
as a three-parameter spline). The SI-MEC has an unpleasant global coupling (equivalent to finding
the tension to apply to the mechanical spline) that makes a global solver considerably slower thanthe corresponding pure MEC [65].
The SI-MEC is extensional, which follows directly from its variational formulation as the curve
minimizing an energy functional, such that it passes through all the control points. It shares G2 continuity and similar locality properties with the MEC.

4.4

Cubic Lagrange (with length parametrization)

The cubic Lagrange spline is defined as follows (see [49, p. 23] for more detail), given a sequenceof input data (ti , zi ), where zi can be interpreted as the vector (xi , yi ):Li,3 (t) =

i+2
X
j=i−1

zj

ωi−1,3 (t),0(t − tj )ωi−1,3(tj )ωi−1,3 (t) =i+2
Y(t − tj )

(4.4)

j=i−1

To obtain the Lagrange polynomial normalized by chord lengths, set successive ti+1 so thatti+1 − ti = |zi+1 − zi |. In other words, the parametrization is the same as the chord length. (Sett0 = 0 for convenience, but the shape of the spline doesn’t depend on it.)This spline is a reasonable approximation to MEC when the angles of the control polygon aresmall. It is also not round, as it is based on a polynomial approximation. It shares locality properties
with the MEC as well.
Further, the spline is not extensional, as adding additional points perturbs the chord lengthparametrization.
Without chord length parametrization, the cubic Lagrange spline is extensional, but prone towild excursions when the control points are not spaced evenly.
The cubic Lagrange spline is sometimes used on its own, but its main interest here is as thebasis of the Ikarus spline, discussed in the next section.

4.5

Ikarus (biarc)

The spline used by Peter Karow’s Ikarus system [42] also falls into the category of two-parameter
splines. It was intended for font design, and was particularly well suited for interpolating between“masters” of different weights.
The Ikarus spline uses two completely different mechanisms to achieve the two basic elementsof the two parameter spline: the determination of tangents and the shapes of the curve segments
between control points. To determine the tangents, first a cubic Lagrange polynomial (with normalization based on the chord lengths, as described in Section 4.4) is fit through the control points.
Then the polynomial curve is discarded (but the tangents preserved), and, for each segment, acurve defined by biarcs is substituted. An arbitrary biarc has one additional degree of freedom, cor47

responding to the split point where the two arcs join, and Ikarus defines this somewhat arbitrarily
to provide reasonable results.
Even though the Lagrange polynomial fit can easily incorporate inflection points, the biarc
primitive cannot, so inflection
pointsmust beinexplicitly(b) Changes
are marked
black marked in the input.

(c) Final IK-points
Figure 4.3. Typical letter-shape using Ikarus.
Figure 6. Comparison of the hand-digitized result with the original
The spline is not extensional, as adding points perturbs the chord length parametrization.
Further, roundness is fairly decent, as the biarc becomes a circular arc when the endpoint tangentsare symmetrical. Continuity is only G1 , but the biarc segments keep the spline from exhibiting
extreme variations in curvature typical of purely polynomial-based splines. The spline is reported
to work well when the control points are dense. Karow [43] recommends one control point for every
30 degrees of arc. A typical example is shown in Figure 4.3 (reproduced from [43]). See also Figure
10.8 for a comparison with the splines of this thesis.

4.6

Euler’s spiral

One particularly well-studied two-parameter spline is Euler’s spiral, also known as the spiralof Cornu or the clothoid. More precisely, it is the spline formed by joining piecewise Euler’s spiralsegments with curvature (G2 ) continuity at the knots.
Euler’s spiral, plotted in Figure 4.4 over the range [−5, 5], is most easily characterized by its
intrinsic definition: curvature is linear in arc length.
κ(s) =

s
.
2

(4.5)

Integrating the intrinsic equation yields this equation, known as the Fresnel integrals (historically, one each for the real and complex part):

48

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
-1

-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-0.9

Figure 4.4. Euler’s spiral

Z sx(s) + iy(s) =2

eit dt .

(4.6)

0This integral is also commonly expressedwith a scale factor applied to the parameter. Theq2
resulting curve is similar; its size is
π times that of the first formulation, with the velocity
similarly scaled so that it also has unit velocity. The traditional notation for these integrals is C(s)and S(s), which evoke cosine and sine, respectively.

Z sC(s) + iS(s) =
0

49

2

eiπt /2 dt .

(4.7)

θright

π/2

π/2

θleft

−π/2

Figure 4.5. Euler’s spiral spline primitive over its parameter space.

50

4.6.1

Euler’s spiral as a spline primitive

Euler’s spiral is readily adapted for use as a two-parameter primitive for a spline. Given tangent
angles at the endpoints, find parameters t0 and t1 from which to “cut” a segment from the spiral, sothat the resulting segment meets the endpoint constraints (for now, assume that such parameters
can be found; an efficient algorithm will be given in Section 8.2).
The behavior of this primitive is shown in Figure 4.5. Compare to Figure 3.13, showing thecorresponding behavior for the MEC spline. Note, among other things, that the Euler’s spiral spline
primitive is defined and well-behaved over a much wider range of endpoint tangent angles.

4.6.2

Variational formulation of Euler’s spiral

The spline literature has generally been split into two camps: variational techniques, and the useof explicitly stated curves such as the Euler spiral. Ironically, Euler’s spiral was born of the studyof elastic springs just as was the elastica, but for the most part the literature does not examine
Euler’s spiral from a variational point of view. This section ties these two threads together, showingthat the Euler spiral does indeed satisfy some fundamental variational properties.
First, as noted by Kimia et al [44], the Euler spiral is among the solutions of the equation for
the MVC functional, i.e. the curve that minimizes the MVC energy functional subject to various
endpoint constraints:

EM V C [κ(s)] =

Z l
0

dκ 2
ds
ds

(4.8)

This result can be readily derived by treating the problem as an Euler variational problem
expressed in terms of κ, i.e. F (x, y, y 0 ) = y 02 , with s for x and κ for y in the terms of Equation 3.2,
as stated in Chapter 3. The corresponding Euler–Lagrange equation is particularly simple:
κ00 = 0

(4.9)

The solution family is thus κ = k0 + k1 s, which are simple similarity transformations of thebasic Euler spiral.
While the simplicity of this result is pleasing, it is unsatisfying for a variety of reasons. First,
because the variational equation is written in terms of κ rather than θ, the boundary conditions
are written in terms of curvature rather than tangent angle. Similarly, the additional terms corresponding to Equation 3.5 representing the positional endpoint constraints are missing. Recall thatwithout these constraints, the only solution to the MEC equation is a straight line.
So this result only tells us that the Euler spiral is the MVC-optimum curve satisfying curvature
constraints on the endpoints, but with unconstrained endpoint positions and tangent angles. As
we shall explore more deeply in Section 7.2.2, the solution to the MVC under other constraint
configurations will in general not be an Euler spiral, a point not clearly made in the presentation
by Kimia et al.
A more fundamental problem is that the degree of continuity doesn’t match. A spline made of51

piecewise Euler spiral segments has G2 -continuity, the same as a MEC spline, while the solution to
the MVC functional with positional constraints for internal knots has G4 -continuity.
So the Euler spiral is a reasonably good first-order approximation to the MEC functional, andis a special case in the solution space of the MVC functional, but might there be a functional forwhich it is actually the optimum? In a sense, this is the reverse of the usual variational problem,in which a function is sought that minimizes a given function. Here, we seek a functional whichhappens to be minimized by the given function, in this case the Euler spiral.
We conjecture that there is indeed a functional: the maximum rate of curvature change overthe entire arc length,
EMAXVC [κ(s)] = max

dκ
.
ds

(4.10)

This functional is related to the MVC; while the MVC is the 2-norm of curvature variation,
EMAXVC corresponds to the ∞-norm.

4.6.3

Euler spiral spline in use

How does the Euler spiral perform in the intended application domain of drawing glyph outlines?
Experimentation yields encouraging results. One such example, a lowercase ‘c’ traced from 14 point
metal Linotype Estienne, is shown in Figure 4.6. Estienne is a distinguished book face designed by
the eminent English printer George W. Jones, and released by Linotype in 1930 [58, p. 133]. It is
not currently available in digital form, making it an attractive candidate for digitization.
The ‘c’ is a simple closed curve, reasonably smooth everywhere, but with regions of high curvature at both terminals. Below the glyph itself is a curvature plot, beginning at the knot marked
by an arrow and continuing clockwise. The long flat section on the left is the inner contour, thesharp spike is the lower terminal, and the other long flat section (which is both longer and lower
curvature) is the outer contour. Note that the curve has two inflection points, corresponding tozero crossings on the curvature plot and marked on the outline with short lines.
A more complex example is shown in Figure 4.7, in this case a lowercase ‘g’ from the same
font. In this case, there are three closed curves; the outer shape and two “holes.” While the top
hole is a simple smooth curve, the other curves contain both G2 -continuous smooth knots and“corner points”, with only G0 -continuity. Such curves can be factored into sections consisting of
the corner points at the endpoints and the smooth knots in the interior, then joined together. Notethat sections with no smooth knots are straight lines (there is one at the top of the “ear”), andsections with a single smooth knot are circular arcs (two of these make up the rest of the ear). Moregenerally, each spline segment bordered by a corner point is a circular arc, with zero variation ofcurvature.
Experience with drawing a significant number of glyphs is encouraging. The Euler spiral spline
can represent a wide range of curves expressively and with an economy of points. Further, in aninteractive editor, the correspondence between the input (the placement of knots and control points)
seems considerably more intuitive than the standard Bézier curve editing techniques.
The trickiest curves to draw using pure Euler’s spiral splines are transitions between straightlines and curves, as typical in the shape of a “U.” The difficulty is illustrated in Figure 4.8. Es52

Figure 4.6. Closed Euler example.

sentially, the problem is that the segment adjoining the corner point tends to be a circular arc
segment, rather than being constrained to be a straight line. It’s possible, by careful interactive
placement, to decrease this curvature so that it’s almost a straight line, but then any changes to thecurved segment will undo the effect of this careful placement. Similar problems occur if the points
are determined by interpolation between different masters, or some other parametric variation. We
shall return to the problem of straight-curve joins in Section 7.4.

4.6.4

Existence and uniqueness of Euler spiral spline

One of the serious shortcomings of the MEC spline is lack of existence for all configurations ofcontrol points (see Section 2.5 for a general discussion). An appealing feature of the Euler spiral
spline is that it solutions exist in many cases when no such solution exists for the MEC spline.

53

Figure 4.7. Complex Euler example.

54

Figure 4.8. Difficulty with straight-to-curved joins.

55

The obvious question is then, does a solution to the Euler spiral spline always exist, no matter
what sequence of points are given? Also, is the solution unique?
A partial answer is given in a 1982 paper by Stoer [82].1 Stoer considered a slightly weaker
formulation of the Euler spiral spline: each segment is piecewise an Euler spiral segment, and there
is G2 -continuity, but there are no additional constraints on the endpoints. Given these relaxed
constraints, there is always a solution, in fact infinitely many. In other words, for any given sequenceof control points, there exist infinitely many assignments of Euler spiral parameters such that theresulting segments have G2 -continuity.
The proof is dependent on the following theorem:
Stoer’s theorem 2.25 (restated): For all points P0 and P1 , tangent angles θ0 , andcurvatures κ0 , there exist countably many different Euler spiral segments connecting P0
to P1 such that the tangent at P0 is θ0 and the curvature at P0 is κ0 .

Cs(s, ))FIGURE6

Figure 4.9. Multiple windings for Euler spiral segment connecting two points (from [82]).
Because of (2.10) and using the same reasoning as with (1.8), every Cornu spiral CA(s), 0 < A A
Stoer’s
proof
is basedon the0 fact
that
Euler
with gently
curvaturecuts dA
infinitely
often
at abscissae
S sI(A)<s 2an(A)< . . .,spiral
whichsegment
satisfy estimatesof thevarying
form [cf.
(1.8)].
may wind around arbitrarily many times before reaching P1 . In the limit, the total winding numberof the spiral segment increases by approximately 2π for each successive solution. Such multiple+ 3n/2 4from*A(s,+1
windings are shown in Figure +A(sl(A))
4.9, reproduced[82, p.(Al)
329].
1 --

1

(2.14)

2
The result that
there
piecewise-Euler
spiral,2(n-1)nC are
2(n -infinitelyn C manyt+sn(A))C 2(n + - )r C
2(n+GI)u.-continuous interpolating
splines follows readily from Stoer’s theorem 2.25. Start at the first point, choosing arbitrary tangent
forn >curvature,
1 , sothat then choose one of the solutions guaranteed by the theorem, and repeat until theandfinal point is reached.
>
+1 3n
1+
+
3n'
s.+A(A)- sn(A)
With arbitrarychoices,such a procedure will tend to produce splines with excessive winding.

X+ASn(k)\[x+Asn(A)12(2.15)

Stoer proposed using the MEC energy to select a single “best” spline from the myriad alternatives.
Stoer also gave an algorithm based
standard
band-diagonal
Jacobian matrix, and claimed
S,-l0 on <thesnQd
6 sC +i(A),that it gave good results in general, but did not include a proof that the algorithm achieves a global
whereoptimum in all cases.
In addition, using the global minimum of MEC energy will also have an effect on the curvature
1+
4mnA
sm =MEC
4mrnover(i1 +all curves
at the endpoints. Since minimizing
assigns zero curvature to endpoints (see

VTx2

1

)

Although, the reader is warned, this paper contains numerous errors, including in the statement of the crucial
2.25.of the quadratic equation,
isTheorem
the solution=
hA(s)=-Xs+ A 2 56
2mn

As a consequence,Section 3.10), and the Euler spline primitive is a fairly good approximation of the MEC primitive,at least in the small-angle assumption the curvature at the endpoints would be approximately zero.
It is not currently known whether a solution always exists in the presence of additional endpoint
constraints, such as zero curvature variation for the terminal segments. Obviously, with the methodsof Stoer it is possible to apply such a constraint to one side, but it’s not clear whether there, for
all input sequences of points, an assignment of tangent for the left side such that the curvature
variation is zero on the right side.
The numerical techniques proposed later in this thesis (Section 8.3 converge on a solution (anda solution with no excess windings, as well) for reasonable input sequences of control points, but
can fail to converge when turning angles become excessive.

4.7

Hobby’s spline

A significant interpolating spline, especially in font design, is Hobby’s spline as implemented
in Metafont [37]. It is perhaps easiest to understand this spline as a very computationally efficientapproximation to the Euler spiral spline of the previous section.
The definition of the primitive segment, with endpoint tangent angles θ0 and θ1 , is√a = 2,1, √b = 16c = (3 − 5)/2,
α = a(sin θ0 − b sin θ1 )(sin θ1 − b sin θ0 )(cos θ0 − cos θ1 ),2+α
ρ = 1+(1−c) cos
θ0 +c cos θ1 ,2−α
σ = 1+(1−c) cos θ1 +c cos θ0 ,x0 , y0 = 0, 0,
x3 , y3 = 1, 0,ρρcos θ0 , 3τ
sin θ0 ,x1 , y1 = 3τσσx2 , y2 = 1 − 3τcos θ1 , 3τ
sin θ1 .

(4.11)The (xi , yi ) are the coordinates for a standard cubic Bézier curve. The adjustable parameterτ corresponds closely to the “tension” parameter of B-splines, and its default value is 1. In the
Metafont sources for Computer Modern, the most common other value found is 0.9, and that is
only for a few glyphs, including lowercase ‘c’ [46, p. 311].
The behavior of the Hobby spline primitive over its parameter space can be seen in Figure 4.10
(see Figure 3.13 for comparison to the MEC). It is defined for all values of θ0 and θ1 , but the figure
only shows values up to 2.156 radians (123.5 degrees).
Similarly, the curvature plot is shown in Figure 4.11. For relatively small angles, its approximation to linear curvature is good, but sharp curvature peaks are apparent at large angles.
While the Euler spiral spline is defined as the one which is a piecewise Euler spiral and G2 continuous at the control points, the Hobby spline uses a slightly different criterion. Hobby introduces the concept of mock curvature, which is the linear approximation (taking only the first-orderterm of the exact equation) of endpoint curvature as a function of the endpoint angles θ0 and θ1 .

57

θright

π/2

π/2

θleft

−π/2

Figure 4.10. Hobby spline primitive over its parameter space.

58θright

π/2

π/2

θleft

−π/2

Figure 4.11. Curvature plot of Hobby spline over its parameter space.

59

The global solution to the spline is the assignment to the vector of tangent angles so that mock
curvature is equal on either side of each knot.
A very important consequence is that solving the global spline is a linear problem. In fact, sinceeach tangent angle affects two neighboring segments, and each segment in turn contributes to themock curvature constraint at its two endpoints, the entire system of constraints can be represented
as a tridiagonal linear equation, which is readily solvable in O(n) using textbook techniques. Hobby
also proves that the matrix is nonsingular over a reasonable range of values for τ , including τ = 1.
Thus, Hobby’s is one of the few “advanced” splines for which there is a solution for all configurationsof input points.
Evaluating the spline properties, we find that the Hobby spline does fairly well. For small
angles, it closely approximates the Euler spiral spline. It is not round, but converges to roundness
as O(n6 ) as opposed to the MEC, which is only O(n4 ). Similarly, the extensionality and continuityof curvature are only approximate, and only good at small angles. However, as noted above, it is anexceedingly robust spline. It is thus a good choice if computational resources are limited, robustness
is key, and the input control points are closely enough spaced that angles are small.
Another two-parameter curve family closely related to the Hobby spline is the primitive usedin Potrace [79]. While Potrace is not an interpolating spline (rather, it is a tool for automaticallyconverting bitmaps into vector curve outlines), its use of a primitive curve family is conceptually similar. Like the Hobby spline, the potrace primitive curve family is based on cubic Bézier
curves (inherently four parameters), with mathematical equations to determine the additional two
parameters.

4.8

Splines from fixed generating curves

The MEC has the property that all segments between control points can be cut from a single
generating curve, the rectangular elastica (here, “cut” means selecting points on the generating
curve, then performing the unique rotation, scaling and translation to make these endpoints matchthe corresponding control points). This property can be derived from the equilibrium-of-forces
formulation of the elastica, and the fact that the constraints at the control points exert zero tension,
but these facts are not intuitively obvious.
Similarly, in the Euler spiral spline, segments between control points are cut from the Euler
spiral (“cut” may also include a mirror flip, because this curve lacks a mirror symmetry).
What do splines cut from a single generating curve have in common? Why do so few other
splines have this property? Might there be some additional flexibility that can be gained by choosing
segments from a curve family rather than a fixed generating curve?
The key to answering these questions is extensionality. In fact, as demonstrated in Section 4.1,
the class of two-parameter extensional splines is exactly equal to the class of two-parameter splines
cut from a generating curve.
Other properties, such as roundness, can then be derived from properties of the generating
curve. For example, for the spline to be round, the curve must tend to a circular shape in the limit.
More precisely, lims→∞ κ0 (s)/κ(s)2 = 0. The Euler spiral meets this condition because κ tends toinfinity, while κ0 is fixed. Another possibility would be a curve tending to fixed curvature in the60

asymptote, in which case κ0 would tend to 0. The MEC, by contrast, clearly does not have thisproperty.
The envelope of endpoint chord angles (as plotted in Figure 3.13 for the MEC) also has aneffect. A small envelope yields failure to exist. A generating curve without an inflection (such asa parabola) excludes all antisymmetric endpoint angle configurations from this envelope, even fortiny angles.

0.5

1

1.5

2

2.5

3

3.5

4

Figure 4.12. Test instrument for fairness user study.

4.8.1

Log-aesthetic curve splines

A very promising family of candidates for generator curves is the aesthetic curve, also known as
the log-aesthetic curve [63]. These curves are defined as having a straight line when κ0 /κ is plottedagainst κ on a log-log scale. In addition to a number of special cases (including the equiangularspiral and Nielsen’s spiral), the general formula for aesthetic curves is κ = sα , where the exponent61

α is the parameter that distinguishes the members of this family. An intuitive motivation for theaesthetic curve is that the eye is not as sensitive to absolute variation of curvature as it is to relative
variation. In areas of high curvature, high variation of curvature is more easily tolerated; but in
low-curvature areas, even small variations in curvature are visible.
The literature on aesthetic curves doesn’t tend to include curves with inflection points – on alog-log plot, these are off the scale. However, since it is well known that the aesthetic curve with
α = 1 is the Euler spiral, curves with inflection points are known to be part of the family, but it’s
possible to generalize the formula to κ = sign(s)|s|α to admit a continuous range of exponents. A
family of such curves (used to interpolate a zigzag sequence of control points) is shown in Figure
4.12.
As described in Section 2.1.1, an empirical user study establishes that MEC energy is not an
accurate predictor of perceived fairness. MEC energy is minimized at α = 0.9, so if MEC energy
were an accurate predictor of fairness, that would have been the preferred value. Since the mean
preferred value was α = 1.78, that hypothesis is refuted.
1

falloff

0.9
0.8
0.7

falloff factor

0.6
0.5
0.4
0.3
0.2
0.1
0
0

1

2

3

4

5

exponent

Figure 4.13. Perturbation fall-off rate as a function of aesthetic curve exponent.
The exponent of the aesthetic curve used as the generator curve in a spline also has a profoundeffect on the locality of the resulting spline. In general, the weaker the curvature variation near the
inflection point of the generator curve, the better the locality. The aesthetic curve family varies
from very steep to very shallow curvature variation across its parameter space. The falloff factor fora single perturbed point, a quantitative measurement of locality, is plotted against the exponentof the aesthetic curve in Figure 4.13. For α = 1, which represents the Euler spiral, the amplitude62

of the ripple diminishes by 0.268 from one crest to the next one. This ratio can be determined
analytically by using the small-angle approximationand solving for κ0 θ1 = κ1 θ0 for the endpoints√of a segment, which yields a value of 1/(2 + 3) for the Euler spiral. Note also that the MEC spline
has the same locality properties as the Euler spiral spline. For the user-preferred exponent of 1.78,
the falloff factor is 0.173, noticeably smaller and thus better. It is possible to achieve even better
locality at the expense of higher curvature variation, if that is desired.

4.9

Piecewise SI-MEC splines, or Kurgla 1

The Autokon system of Mehlum [59] was inspired by the goal of simulating a real, mechanicalspline, but its kurgla 1 algorithm actually implemented something a bit different, with bothadvantages and disadvantages with respect to the MEC it was intended to mimic.
Mehlum based his design on the fact that the curvature of an elastica is linear in distance fromthe line of action. His solver approximated the integration of this intrinsic equation as a sequenceof circular arcs, the curvature of each determined by the distance from the line of action.
As we’ve seen, the general elastica has three degrees of freedom; two for the position and angleof the line of action with respect to the control point chord, and one for the amplitude of thecurvature variation. Note that inflections occur at the line of action. In this formulation, choosing
all these parameters to minimize the total bending energy is a tricky global optimization problem,
and kurgla 1 did not attempt it.
In particular, it nailed down the angle of the line of action, choosing it to be perpendicular
to the chord. In this way, the parameter space was reduced from three to two, and the problem
became much more tractable. In particular, the curve segments were determined by the tangentsat the endpoints, and the remaining global problem of choosing tangents for G2 -continuity yields
to simple Newton-style solution.
As shown in Section 3.11, an elastica segment with the line of action perpendicular to the chordis the one minimizing the SI-MEC energy for the angle constraints at those endpoints. Thus, thekurgla 1 algorithm effectively computes a piecewise SI-MEC.
One consequence is that kurgla 1 is round, unlike the pure MEC, since the SI-MEC describes a circular arc in the case of symmetrical endpoint angle constraints. Another consequenceis considerably better robustness than the pure MEC, because the envelope of parameter space isapproximately twice as large.
However, the price paid for these improvements is loss of extensionality; when a new on-curve
point is added, it generally changes the chords, and thus the lines of action for the resulting elastica
sections.
Thus, though it is based on the elastica, and has some desirable properties, kurgla 1 is neither
a true energy minimizing spline like the MEC, nor does its extensionality and robustness compare
well to the Euler spiral spline. Indeed, Mehlum was aware of these limitations, and his subsequent
kurgla 2 algorithm is the earliest known implementation of the Euler spiral spline.2
The MEC, SI-MEC, and Euler spiral are closely related. One way to visualize this relation is
2

Mehlum wrote to George Forsythe on July 9, 1971, mentioning “Fresnel integral splines” as recent work. The
work was published in 1974 [59].

63Euler spiral

SI-MEC

MEC

Figure 4.14. MEC, SI-MEC, and Euler spiral compared.

to consider the gradient of curvature. In an Euler spiral, curvature increases linearly along the arc
length of the curves. In a SI-MEC defined by two endpoint tangents, curvature increases linearly
along the chord between the endpoints. And in a pure MEC, curvature increases along the line
tangent to the inflection point. Figure 4.14 shows these gradients of curvature visually – along eachof the curves, tick marks meter the rate of change of curvature, and the tick mark is angled sothat the rate of curvature change is constant between parallel lines. This figure also shows that the
SI-MEC and Euler spiral are visually very similar.

4.10

Circle splines

Interpolating splines have a tradeoff between extensionality and locality. The two-parameter
splines described above favor extensionality. It is possible to make different choices, and favor
locality instead. An excellent example is the “circle spline” of Séquin, Lee, and Yin [80].
In the circle spline, the tangent angle at each control point is chosen as the tangent of the circle
going through the control point and its two neighbors. Each segment is then constructed using a
blending of these circles from the tangents at the endpoints. Thus, moving a control point only
affects four segments – two on either side.
The construction of the primitive curves is robust and smooth. In fact, the construction guarantees zero κ0 at the endpoints, and that the curvature is identical to that of the circle used todetermine tangents, so the resulting spline has G3 -continuity.
Unfortunately, even though the spline has good locality and G3 -continuity, it is not particularlysmooth for all inputs, and adding additional on-curve control points can cause serious deviations

64

Figure 4.15. Circle spline extensionality failure.

from extensionality, in many cases adding inflections. Figure 4.15 clearly shows the lack of extensionality; the two paths differ only by the addition of the central point in the right-hand figure.

4.11

SummaryThis chapter described an important class of interpolating splines, those with a primitive curve
described by a two-dimensional parameter space. A great many splines in the literature fall into
this category, notably the Minimum Energy Curve.
A central result of this work is that any two-parameter, extensional spline is defined in terms ofa fixed generating curve. These properties hold across a large family of desirable splines, and makes
the choice of an optimum spline a much simpler problem. Choosing a most preferred generating
curve can be done empirically. By contrast, the variational approach would require developing an
accurate model of the perception of curve fairness. Empirical testing establishes that the mostcommon such metric, bending energy, does not coincide accurately with user preference.
Table 4.11 summarizes the properties of several splines in this two-parameter family. All ofthese are known and appear in the literature.
• MEC is the interpolating spline minimizing the total bending energy. Each segment between
control points is a piece of a rectangular elastica, and the tangents are chosen so that G2
continuity is preserved across each control point. (Section 3.10)• Euler spiral is the interpolating spline for which each segment between control points is asection of an Euler spiral, i.e. curvature is linear in arc length. Like the MEC, the tangentsare chosen for G2 -continuity. (Section 6.8)
65

Curve type

Tangents

Continuity

Extensionality

Roundness

Robustness

Locality

MEC

G2

G2

+++

−−

−−

+

Euler spiral

2

G

2

G

+++

+++

++

+

log-aesthetic

G2

G2

+++

+++

−

++

Hobby

≈ G2

≈ G2

−

−

+++

+

Biarc (Ikarus)

Lagrange

G1

−−

+++

++

+

kurgla 1

G2

G2

−

+++

+

+

local

3

−−−

+++

+++

+++

CircleG

Table 4.1. Comparison of properties of two-parameter splines.

• log-aesthetic is the two-parameter extensible spline based on log-aesthetic curves with an
exponent near α = 2. (Section 4.8.1)• Ikarus (Section 4.5).
• Hobby is a piecewise cubic spline developed for the Metafont system. It is best understoodas a polynomial approximation to the Euler spiral spline, with better robustness but lower
performance on other metrics such as smoothness and extensionality [37] (Section 4.7).
• kurgla 1 (Section 4.9).
• Circle spline (Section 4.10).
For the basic interpolating spline problem, two-parameter extensible splines seem likely to be
the best possible choice. From this family, while the MEC (based on the rectangular elastica) is
the original inspiration for virtually the entire field of splines, the Euler spiral spline is superior
and has excellent all-around properties, and is especially well-suited for interactive use due to its
stability. The elastica and Euler spiral are particularly important and beautiful curves. In addition,
they have a fascinating and intertwined history, explored in detail in the next two chapters.

66

Chapter 5

History of the elasticaThis chapter tells the story of a remarkable family of curves, known as the elastica, Latin for athin strip of elastic material. The elastica caught the attention of many of the brightest minds inthe history of mathematics, including Galileo, the Bernoullis, Euler, and others. It was present at
the birth of many important fields, most notably the theory of elasticity, the calculus of variations,and the theory of elliptic integrals. The path traced by this curve illuminates a wide range ofmathematical style, from the mechanics-based intuition of the early work, through a period oftechnical virtuosity in mathematical technique, to the present day where computational techniques
dominate.
There are many approaches to the elastica. The earliest (and most mathematically tractable),is as an equilibrium of moments, drawing on a fundamental principle of statics. Another approach,ultimately yielding the same equation for the curve, is as a minimum of bending energy in theelastic curve. A force-based approach finds that normal, compression, and shear forces are also inequilibrium; this approach is useful when considering specific constraints on the endpoints, which
are often expressed in terms of these forces.
Later, the fundamental differential equation for the elastica was found to be equivalent to thatfor the motion of the simple pendulum. This formulation is most useful for appreciating the curve’s
periodicity, and also helps understand special values in the parameter space.
In more recent times, the focus has been on efficient numerical computation of the elastica(especially for the application of fitting smooth spline curves through a sequence of points), andalso determining the range of endpoint conditions for which a stable solution exists. Many practicalapplications continue to use brute-force numerical techniques, but in many cases the insights ofmathematicians (many working hundreds of years earlier) still have power to inspire a more refinedcomputational approach.

67

5.1

Jordanus de Nemore – 13th century

In the recorded literature, the problem of the elastica was first posed in De Ratione Ponderis
by Jordanus de Nemore (Jordan of the Forest), a thirteenth century mathematician. Proposition 13
of book 4 states that “when the middle is held fast, the end parts are more easily curved.” He thenposes an incorrect solution: “And so it comes about that since the ends yield most easily, while theother parts follow more easily to the extent that they are nearer the ends, the whole body becomes
curved in a circle.” [17]. In fact, the circle is one possible solution to the elastica, but not to thespecific problem posed. Even so, this is a clear statement of the problem, and the solution (though
not correct) is given in the form of a specific mathematical curve. It would be several centuries
until the mathematical concepts needed to answer the question came into existence.

5.2

Galileo sets the stage – 1638

Three basic concepts are required for the formulation of the elastica as an equilibrium of moments: moment (a fundamental principle of statics), the curvature of a curve, and the relationshipbetween these two concepts. In the case of an idealized elastic strip, these quantities are linearly
related.

Figure 5.1. Galileo’s 1638 problem.

Galileo, in 1638, posed a fundamental problem, founding the mathematical study of elasticity.
Given a prismatic beam set into a wall at one end, and loaded by a weight at the other, howmuch weight is required to break the beam? The delightful Figure 5.1 illustrates the setup. Galileo

68considers the beam to be a compound lever with a fulcrum at the bottom of the beam meeting thewall at B. The weight E acts on one arm BC, and the thickness of the beam at the wall, AB, isthe other arm, “in which resides the resistance.” His first proposition is, “The moment of force at
C to the moment of the resistance... has the same proportion as the length CB to the half of BA,
and therefore the absolute resistance to breaking... is to the resistance in the same proportion as
the length of BC to the half of AB...”
From these basic principles, Galileo derives a number of results, primarily a scaling relationship. He does not consider displacements of the beam; for these types of structural beams, thedisplacement is negligible. Even so, this represents the first mathematical treatment of a problem
in elasticity, and firmly establishes the concept of moment to determine the force on an elastic
material. Many researchers elaborated on Galileo’s results in coming decades, as described in detailin Todhunter’s history [85].
One such researcher is Ignace-Gaston Pardies, who in 1673 posed one form of the elastica
problem and also attempted a solution: for a beam held fixed at one end and loaded by a weight at
the other, “it is easy to prove” that it is a parabola. However, this solution isn’t even approximately
correct, and would later be dismissed by James Bernoulli as one of several “pure fallacies.” Truesdell
gives Pardies credit for introducing the elasticity of a beam into the calculation of its resistance,
and traces out his influence on subsequent researchers, particularly Leibniz and James Bernoulli[86, p. 50–53].

5.3

Hooke’s law of the spring – 1678

Hooke published a treatise on elasticity in 1678, containing his famous law. In a short Latin
phrase posed in a cryptic anagram three years earlier to establish priority, it reads, “ut tensio sicvis; that is, The Power of any Spring is in the same proportion with the Tension thereof... Now as
the Theory is very short, so the way of trying it is very easie.” (Spelling and capitalization as in
the original). In modern formulation, it is understood as F ∝ ∆l; the applied force is proportional
to the change in length.

Figure 5.2. Hooke’s figure on compound elasticity.

Hooke touches on the problem of elastic strips, providing an evocative illustration (Figure 5.2),
but according to Truesdell, “This ‘compound way of springing’ is the main problem of elasticity

69

for the century following, but Hooke gives no idea how to relate the curvature of one fibre to thebending moment, not to mention the reaction of the two fibres on one another.” [86, p. 55]
Indeed, mathematical understanding of curvature was still in development at the time, andmastery over it would have to wait for the calculus. Newton published results on curvature in
his “Methods of Series and Fluxions,” written 1670 to 1671 [35, p. 232], but not published forseveral more years. Leibniz similarly used his competing version of the calculus to derive similar
results. Even so, some results were possible with pre-calculus methods, and in 1673, Christiaan
Huygens published the “Horologium oscillatorium sive de motu pendulorum ad horologia aptato
demonstrationes geometrica,” which used purely geometric constructions, particularly the involutes
and evolutes of curves, to establish results involving curvature. The flavor of geometric construction
pervades much of the early work on elasticity, particularly James Bernoulli’s, as we shall see.
Newton also used his version of the calculus to more deeply understand curvature, and provided
the formulation for radius of curvature in terms of Cartesian coordinates most familiar to us today[35],3ρ = (1 + y 02 ) 2 /y 00 .

(5.1)

An accessible introduction to the history of curvature is the aptly named “History of Curvature”
by Dan Margalit [56].
Given a solid mathematical understanding of curvature, and assuming a linear relation between
force and change of length, working out the relationship between moment and curvature is indeed“easie,” but even Bernoulli had to struggle a bit with both concepts. For one, Bernoulli didn’t
simply accept the linear law of the spring, but felt the need to test it for himself. And, since hehad the misfortune to use catgut, rather than a more ideally elastic material such as metal, hefound significant nonlinearities. Truesdell [88] recounts a letter from James Bernoulli to Leibniz
on 15 December 1687, and the reply of Leibniz almost three years later. According to Truesdell,this exchange is the birth of the theory of the “curva elastica” and its ramifications. Leibniz had
proposed: “From the hypothesis elsewhere substantiated, that the extensions are proportional to thestretching forces...” which is today attributed as Hooke’s law of the spring. Bernoulli questioned
this relationship, and, as we shall see, his solution to the elastica generalized even to nonlinear
displacement.

5.4

James Bernoulli poses the elastica problem – 1691

James Bernoulli posed the precise problem of the elastica in 1691:“Si lamina elastica gravitatis espers AB, uniformis ubique crassitiei & latitudinis, inferiore
extremitate A alicubi firmetur, & superiori B pondus appendatur, quantum sufficit ad laminam
eousque incurvandam, ut linea directionis ponderis BC curvatae laminae in B sit perpendicularis,
erit curvatura laminae sequentis naturae:”
And then in cipher form:

70Figure 5.3. Bernoulli poses the elastica problem.

“Portio axis applicatam inter et tangentem est ad ipsam tangentem sicut quadratum applicatae
ad constans quoddam spatium.”1
Assuming a lamina AB of uniform thickness and width and negligible weight of its own,
supported on its lower perimeter at A, and with a weight hung from its top at B, theforce from the weight along the line BC sufficient to bend the lamina perpendicular, thecurve of the lamina follows this nature:
The rectangle formed by the tangent between the axis and its own tangent is a constantarea.
This poses one specific instance of the general elastica problem, now generally known as therectangular elastica, because the force applied to one end of the curve bends it to a right angle withthe other end held fixed.
The deciphered form of the anagram is hardly less cryptic than the original, but digging throughhis 1694 explanation, it is possible to extract the fundamental idea: at every point along the curve,
the product of the radius of curvature and the distance from the line BC is a constant, i.e. the twoquantities are inversely proportional. And, indeed, that is the key to unlocking the elastica; theequation for the shape of the curve follows readily, given sufficient mathematical skill.

5.5

James Bernoulli partially solves it – 1692

By 1692, James Bernoulli had completely solved the rectangular case of the elastica posed
earlier. In his Meditatione CLXX dated that year, titled “Quadratura Curvae, e cujus evolutione
describitur inflexae laminae curvatura” [7], or, “Quadrature of a curve, by the the evolution ofwhich is traced out the curve of a bent lamina” he draws a figure of that curve (reproduced here
as Figure 5.4), and gives an equation for its quadrature2 .
1

The cipher reads “Qrzumu bapt dxqopddbbp...” and the key was published in the 1694 Curvatura Laminae withthe detailed solution to the problem. Such techniques for establishing priority may seem alien to academics today,
but are refreshingly straightforward by comparison to the workings of the modern patent system.
2
Today, we would say simply “integral” rather than quadrature.

71

B

C

A
F

E
Q

D

G

H
L

I

M

N

Figure 5.4. Drawing from James Bernoulli’s 1692 Med. CLXX, with modern reconstruction.

Radius circuli AB = a, AED lamina elastica ab appenso pondere in A curvata, GLM
illa curva, ex cujus evolutione AED describitur: AF = y, F E = x, AI = p, IL = z.
Aequatio differentialis naturam curvae AED exprimens,
xx dx
dy = √,
a4 − x4
ut suo tempore ostendam:
A rough translation:
Let the radius of the circle AB = a, AED be an elastic lamina curved by a suspended
weight at A, and GLM that curve, the involute of which describes AED. AF = y,
F E = x, AI = p, IL = z. The differential equation for the curve AED is expressed,
x2 dx
dy = √,
a4 − x4(5.2)

as I will show in due time.
Equation 5.2 is a clear, simple statement of the differential equation for the rectangular instanceof the elastica family, presented in readily computable form; y can be obtained as the integral of astraightforward function of x.
Based on the figure, AED is the elastica, and GLM is its evolute. Thus, the phrase “evolutione
describitur” appears to denote taking the involute of the curve GLM to achieve the elastica AED;72

L is the point on the evolute corresponding to the point E on the elastica itself. L has Cartesian
coordinates (p, z), likewise E has coordinates (x, y). It is clear that Bernoulli was working with theevolute as the geometric construction for curvature, which is of course central to the theory of theelastica. In particular, by the definition of the evolute, the length EL is the radius of curvature atpoint E on the curve AED.

5.6

James Bernoulli publishes the first solution – 1694

Figure 5.5. Bernoulli’s 1694 publication of the elastica.
Bernoulli held on to this solution for a couple of years, and finally published in his landmark

73

1694 Curvatura Laminae Elasticae 3 . See Truesdell [86, pp. 88–96] for a detailed description of thiswork, which we will only outline here.

Figure 5.6. Bernoulli’s justification for his formula for curvature.

Bernoulli begins the paper by giving a general equation for curvature (illustrated in Figure 5.6),
which he introduces thus, “in simplest and purely differential terms the relation of the evolventof radius of the osculating circle of the curve... Meanwhile, since the immense usefulness of this
discovery in solving the velaria, the problem of the curvature of springs we here consider, and otherrecondite matters makes itself daily more and more manifest to me, the matter stands that I cannot
longer deny to the public the golden theorem.” Bernoulli’s formulation is not entirely familiar to
the modern reader, as it mixes infinitesimals somewhat promiscuously. He set out this formula forthe radius of curvature z:
z=

dx ds
dy ds
=
ddy
ddx

(5.3)

From Bernoulli’s tone, it is clear he thought this was an original result, but Huygens and Leibniz
were both aware of similar formulas; Huygens had already published a statement and proof in termsof pure Cartesian coordinates (which would be the form most familiar today). However, in spiteof this knowledge, both considered the problem of the elastica impossibly difficult. Huygens, in a
letter to Leibniz dated 16 November 1691, wrote, “I cannot wait to see what Mr. Bernoulli theelder will produce regarding the curvature of the spring. I have not dared to hope that one would
come out with anything clear or elegant here, and therefore I have never tried.” [86, p. 88, footnote
4]
Bernoulli’s treatment of the elastica is fairly difficult going (as evidenced by the skeptical
reaction and mistaken conclusions from Leibniz and others), but again, it is possible to tease out
the central ideas. First, the idea that the moment at any point along the curve is proportional to the3

Curvatura Laminae Elasticae. Ejus Identitas cum Curvatura Lintei a pondere inclusi fluidi expansi. Radii Circulorum Osculantium in terminis simplicissimis exhibiti, una cum novis quibusdam Theorematis huc pertinentibus,&c.”, or “The curvature of an elastic band. Its identity with the curvature of a cloth filled out by the weight of the
included fluid. The radii of osculating circles exhibited in the most simple terms; along with certain new theorems
thereto pertaining, etc.”. Originally published in the June 1694 Acta Eruditorum (pp. 262–276), it is collected in the1744 edition of his Opera [8, p. 576–600], now readily accessible online.

74

distance from the line of force. Bernoulli writes in a 1695 paper4 explaining the 1694 publication (andreferring to Figure 5.5 for the legend): “I consider a lever with fulcrum Q, in which the thickness Qyof the band forms the shorter arm, the part of the curve AQ the longer. Since Qy and the attachedweight Z remain the same, it is clear that the force stretching the filament y is proportional
to the segment QP .” Next, Bernoulli carefully separated out the force and the elongation, and, infact, allows for a completely arbitrary functional relationship, not necessarily linear (this function isrepresented by the curve AF C in Figure 5.5, labeled “Linea Tensionum”). “And since the elongationis reciprocally proportional to Qn, which is plainly the radius of curvature, it follows that Qn is also
reciprocally proportional to x.” A central argument here is that the moment (and hence curvature,assuming a linear relationship) is proportional solely to the amount of force and the distance fromthe line of that force; the shape of the curve doesn’t matter.
F

M = Fxx

Figure 5.7. A simplified diagram of moments.

The moments can be seen in the simplified diagram in Figure 5.7, which shows the curve of theelastica itself (corresponding to AQR in Figure 5.5), the force F on the elastica (corresponding to
AZ in Figure 5.5), and the distance x from the line of force (corresponding to P Q in Figure 5.5).
According to the simple lever principle of statics, the moment is equal to the force F applied on
the elastica times the distance x from the line of force.
Here we present a slightly simplified version of Bernoulli’s argument (see Truesdell [86, p. 92]
for a more complete version). Write the curvature as a function of x, understanding that in theidealized case it is linear, i.e., f (x) = cx. Then, using Equation 5.3, the “golden theorem” forcurvature:d2 y= f (x) .
dx ds

(5.4)

Integrating with respect to dx and assuming dy/ds = 0 at x = s = 0 gives
4
Explicationes, Annotationes et Additiones ad eas quae in actis superiorum annorum de curva elastica, isochrona
paracentrua, et velaria, hinc inde memorata, et partim controversa leguntur; ubi de linea mediarum directionum,
aliisque novis. Originally published in Acta. Eruditorum, Dec. 1695, pp. 537–553, and reprinted in the 1744 Complete
Works [8, pp. 639–663].

75

dy
=
ds

Z x
f (ξ) dξ = S(x) .

(5.5)

0

Substituting the standard identity (which follows algebraically from ds2 = dx2 + dy 2 )
dy
dy/ds
=p
dx
1 − (dy/ds)2

(5.6)

S(x)
dy
=p
.
dx
1 − S(x)2(5.7)

we getAnd in the ideal case where f (x) = 2cx, we have S(x) = cx2 , and thus
cx2
dy
=√
.
dx
1 − c2 x4(5.8)

Which is the same as the Equation 5.2, with straightforward change of constants. As mentioned
before, the use of both ds and dx as the infinitesimal is strange by modern standards, even though
in this case it leads to a solution quite directly. A parallel derivation using similar principles but
only Cartesian coordinates can be found in Whewell’s 1833 Analytical Statics [89, p. 128]; there,
Bernoulli’s use of ds rather than dx can be seen as the substitution of variable that makes theintegral tractable.
It is worth noting that, while James Bernoulli’s investigation of the mathematical curve describing the elastica is complete and sound, as a more general work on the theory of elasticity
there are some problems. In particular, Bernoulli assumes (incorrectly) that the inner curve of theelastica (AQR in Figure 5.5) is the “neutral fiber,” preserving its length as the elastica bends. In
fact, determining the neutral fiber is a rather tricky problem whose general solution can still be
determined only with numerical techniques (see Truesdell [86, pp. 96–109] for more detail). Fortunately, in the ideal case where the thickness approaches zero, the exact location of the neutral fiber
is immaterial, and all that matters is the relation between the moment and curvature.
The limitations of Bernoulli’s approach were noted at the time [24]. Huygens published a shortnote in the Acta eruditorum in 1694, very shortly after Bernoulli’s publication in the same forum, illustrating several of the possible shapes the elastica might take on, and pointing out that
Bernoulli’s quadrature only expressed the rectangular elastica. His accompanying figure is reproduced here as Figure 5.8. The shapes are shown from left to right in order of increasing force atthe endpoints, and shape A is clearly the rectangular elastica.
Bernoulli acknowledged this criticism (while pointing out that he described these other cases
explicitly in his paper, something that Huygens apparently overlooked) and indicated that his
technique could be extended to handle these other cases (by using a non-zero constant for the
integration of Equation 5.5), and went on to give an equation with the general solution, reproducedby Truesdell [86, p. 101] as:(x2 ± ab)dx±dy = pa4 − (x2 ± ab)2
76

(5.9)

Figure 5.8. Huygens’s 1694 objection to Bernoulli’s solution.

It is not clear that the publication of this more general solution had much impact. Bernoulli
did not graph the other cases, nor did he seem to be aware of other important properties of thesolution, notably its periodicity, nor the fact that it includes solutions with and without inflections.
(Another indication that these results were not generally known is that Daniel Bernoulli, in a
November 8, 1738 letter to Leonhard Euler during a period of intense correspondence regarding
the elastica, wrote, “Apart from this I have since noticed that the idea of my uncle Mr. James
Bernoulli includes all elasticas.” [86, p. 109])
Without question, the Bernoulli family set the stage for Euler’s definitive analysis of the elastica.
The next breakthrough came in 1742, when Daniel Bernoulli (the nephew of James) proposed
to Euler solving the general elastica problem with the the technique that would finally crack it:
variational analysis. This general problem concerns the family of curves arising from an elastic bandof arbitrary given length, and arbitrary given tangent constraints at the endpoints.

5.7

Daniel Bernoulli proposes variational techniques – 1742

Daniel Bernoulli, in an October 1742 letter to Euler [6], discussed the general problem of theelastica, but had not yet managed to solve it himself:
Ich möchte wissen ob Ew. die curvaturam laminae elasticae nicht könnten sub hac faciesolviren, dass eine lamina datae longitudinis in duobus punctis positione datis fixirt sey,also dass die tangentes in istis punctis such positione datae seyen. ... Dieses ist die idea
generalissima elasticarum; hab aber sub hac facie noch keine Solution gefunden, als per
methodum isoperimetricorum, da ich annehme, dass die vis viva potentialis laminae
elasticae insita müsse minima seyn, wie ich Ew. schon einmal gemeldet. Auf diese Weise
bekomme ich eine aequationem differentialem 4ti ordinis, welche ich nicht hab genugsam
reduciren können, um zu zeigen, dass die aequatio ordinaria elastica general sey.
A rough translation into English reads:
I’d like to know whether you might not solve the curvature of the elastic lamina under
this condition, that on the length of the lamina on two points the position is fixed, and77

that the tangents at these points are given. ... This is the idea of the general elastica; I
have however not yet found a solution under this condition by the isoperimetric method,
given my assumption that the potential energy of the elastic lamina must be minimal,
as I’ve mentioned to you before. In this way I get a 4th order differential equation,
which I have not been able to reduce enough to show a regular equation for the generalelastica.
This previous mention is likely his 7 March 1739 letter to Euler, where he gave a somewhat lesselegant formulation of the potential energy of an elastic lamina, and suggested the “isoperimetricmethod,” an early name for the calculus of variations. Many founding problems in the calculus ofvariations concerned finding curves of fixed length (hence isoperimetric), minimizing or maximizing some quantity such as area enclosed. Usually additional constraints are imposed to make theproblem more challenging, but, even in the unconstrained case, though the answer (a circle) was
known as early as Pappus of Alexandria around 300 A.D, rigorous proof was a long time coming.
In any case, Bernoulli ends the letter with what is likely the first clear mathematical statementof the elastica as a variational problem in terms of the stored energy5 :Ew. reflectiren ein wenig darauf, ob man nicht könne, sine interventu vectis, die curvaturam ABC immediate ex principiis mechanicis deduciren. Sonsten exprimire ich
die
R ds
vim vivam potentialem laminae elasticae naturaliter rectae et incurvatae durch RR,
sumendo elementum ds pro constante et indicando radium osculi per R. Da Niemand die
methodum isoperimetricorum
R ds so weit perfectionniret, als Sie, werden Sie dieses problema, quo requiritur ut RR
faciat minimum, gar leicht solviren.
You reflect a bit on whether one cannot, without the intervention of some lever, immediately deduce the curvature of ABC from the principles of mechanics. Otherwise,
I’d express the potential energy
R dsof a curved elastic lamina (which is straight when in
its natural position) through RR , assuming the element ds is constant and indicating
the radius of curvature Rby R. There is nobody as perfect as you for easily solving theds
problem of minimizing RR
using the isoperimetric method.
Daniel was right. Armed with this insight, Euler was indeed able to definitively solve the general
problem within the year (in a letter dated 4 September 1743, Daniel Bernoulli thanks Euler for
mentioning his energy-minimizing principle in the “supplemento”), and this solution was publishedin book form shortly thereafter.6

5.8

Euler’s analysis – 1744

Euler, building on (and crediting) the work of the Bernoullis, was the first to completely characterize the family of curves known as the elastica, and published this work as an Additamentum
5

That said, James Bernoulli in 1694 expressed the equivalence of the elastica with the lintearia, which has avariational formulation in terms of minimizing the center of gravity of a volume of water contained in a cloth sheet,as pointed out by Truesdell [86, p. 201].
6
Truesdell also points out that this formulation wasn’tRentirely novel; Daniel Bernoulli and Euler had corresponded
in 1738 about the more general problem of minimizing rm ds, and they seemed to be aware that the special casem = −2 corresponded to the elastica [86, p. 202]. The progress of knowledge, seen as a grand sweep from far away,
often moves in starts and fits when seen up close.

78

(appendix) [20] in his landmark book on variational techniques. His treatment was quite definitive,
and holds up well even by modern standards. This Additamentum is a truly remarkable work,
deserving of deep study. Fortunately for the reader, it accessible both in facsimile due to the Euler
archive, and also in a 1933 English translation by Oldfather [68].
Closely following Daniel Bernoulli’s suggestion, Euler expressed the problem of the elastica veryclearly in variational form. He wrote [20, p. 247],
ut, inter omnes curvas ejusdem longitudinis, quæ non solum per puncta A & B
transeant, sed etiam in his punctis
R dsa rectis positione datis tangantur, definiatur ea inminimus.
qua sit valor hujus expressionis RR
In English:
That among all curves of the same length that not only pass through points A andB but are also tangent to given straightat these points, it is defined as the one
R linesds
minimizing the value of the expression RR .
Here, s refers to arc length, exactly as is common usage today, and R is the radius of curvature(“radius osculi curvæ” in Euler’s words), or κ−1 in modern notation. Thus, today we are more
likely to write that the elastica is the curve minimizing the energy E[κ] over the length of the curve
0 ≤ s ≤ l, whereZ lE[κ(s)] =κ(s)2 ds .

(5.10)

0

While today we would find it more convenient to work in terms of curvature (intrinsic
p equations),
Euler quickly moved to Cartesian coordinates, using the standard definitions ds = 1 + y 02 dx and
y 00
1
0
00
2
2
R = (1+y 02 )3/2 , where y and y represent dy/dx and d y/dx , respectively. Thus, the variational
problem becomes finding a minimum for
Z

y 002
dx .
(1 + y 02 )5/2

(5.11)

This equation is written in terms of the first and second derivatives of y(x), and so the simple Euler–Lagrange equation does not suffice. Daniel Bernoulli had run into this difficulty, as heexpressed in his 1742 letter. In hindsight, we now know that expressing the problem in terms oftangent angle as a function of arc length yields to a first-derivative variational approach, but thisapparently was not clear to either Bernoulli or Euler at the time.
In any case, by 1744, Euler had discovered what is now known as the Euler–Poisson equation, capable of solving variational problems in terms of second derivatives, and he could apply
it straightforwardly to Equation 5.11. (See Section 7.2.1 for more discussion of the Euler–Poisson
equation and its application to a related variational problem.) He thus derived the following generalequation, where a and c are parameters. (See Truesdell for a more detailed explanation of Euler’sderivation [86, pp. 203–204].)79dy
a2 − c2 + x2=p.
dx(c2 − x2 )(2a2 − c2 + x2 )(5.12)

This equation is essentially the same as James Bernoulli’s general solution, Equation 5.9.
Euler then goes on to classify the solutions to this equation based on the parameters a and c.
To simplify constants, we propose the use of a single λ to replace both a and c; this formulation has
a particularly simple interpretation as the Lagrange multiplier for the straightforward variational
solution of Equation 5.10.

λ=

a2
2c2

(5.13)

Note also that this parameter has a straightforward geometrical interpretation. The tangent
angle at the inflection point (relative to the line of action) is cos−1 (1 − .5/λ). Euler also worked
with the slope of the curve at the inflection point, using the mathematically equivalent formulationof Equation 5.13.
Euler observed that there is an infinite variety of elastic curves, but that “it will be worth while
to enumerate all the different kinds included in this class of curves. For this way not only will thecharacter of these curves be more profoundly perceived, but also, in any case whatsoever offered,it will be possible to decide from the mere figure into what class the curve formed ought to be put.
We shall also list here the different kinds of curves in the same way in which the kinds of algebraic
curves included in a given order are commonly enumerated.” [20, §14]. According to the note in
Oldfather’s translation, Euler is referring to Newton’s famous classification of cubic curves [68, p.
152, Note 5]. In any case, Euler finds nine such classes, enumerated in the table below:

Euler’s
species # Euler’s Figure
1
2

6

3
4

7

5

8

6

9

7

10

8

11

9

Euler’s parameters

λ

c=0

λ=∞

0<c<a

0.5 < λstraight linec=a
λ = 0.5
rectangular elastica√.30269 < λ < .5a < c < a 1.65187√
c = a 1.65187
λ = .30269
lemnoid√√
a 1.65187 < c < a 2 .25 < λ < .30269√
c=a 2
λ = .25
syntractrix√
a 2<c0 < λ < .25
a=0

λ=0

Table 5.1. Euler’s characterization of the space of all elastica.

80

circleFigure 5.9. Euler’s elastica figures, Tabula III.

81

Figure 5.10. Euler’s elastica figures, Tabula IV.

82

Euler includes figures for six of the nine cases, reproduced here in Figures 5.9 and 5.10. Of thethree remaining, species #1 is a degenerate straight line, and species #9 is a circle. Species #3, therectangular elastica, is of special interest, so it is rather disappointing that Euler did not include
a figure for it. Possibly, it was considered already known due to Bernoulli’s publication. Note alsothat in this special case, c = a, the equation becomes equivalent to Bernoulli’s Equation 5.2.
Compare these figures with the computer-drawn Figure 5.11. It is remarkable how clearly Euler
was able to visualize these curves, even 250 years ago. It is likely that the distortions and inaccuracies
are due primarily to the draftsman engraving the figures for the book rather than to Euler himself;they display lack of symmetry that Euler clearly would have known. In fact, Euler computed anumber of values to seven or more decimal places, including the values of a and c for the lemnoid
shape (species #5).
Note that for c2 < 2a2 , or λ > .25, the inflectional cases, the equation is well-defined for−c < x < c. In these cases, the parameter c is half the width of the figure. In the non-inflectional
cases (λ < .25), the equation has no solutions at x = 0, and in Euler’s figure the curve is drawn to
the right of the y axis.
The curve described by species #7, the only non-periodic solution, was known to Poleni in
1729, and is also known as the “syntractrix of Poleni,” or, in French, “la courbe des forçats” (thecurve of convicts, or galley slaves). Euler gave its equation (in section 31) in closed form (without
citing Poleni):pcc+
y = c2 − x2 − log
2

5.8.1

√

c2 − x2.
x(5.14)

Moments

Euler was also clearly aware of the simple moment approach to the elastica, and in section
43 of the Additamentum, demonstrated its equivalence to the variational approach. To do so, hemanipulated the quadrature formulation of the elastica, Equation 5.12, grouping all the terms
involving first and second derivatives of the curve into a single term Sq(1 − p2 )−3/2 . Recall that
dp
d2 y
dy
and q = dx
= dx
p = dx2 . Thus, this term is S times the curvature, yielding an equation relating
curvature to Cartesian coordinates. Euler explains:
3:2“But −(1+pp)
is the radius of curvature R; whence, by doubling the constants β andqγ, the following equation will arise:S= α + βx − γy .
R

(5.15)

This equation agrees admirably with that which the second or direct method supplies.
For let α +βx−γy express the moment of the bending power, taking any line you pleaseas an axis, to which moment the absolute elasticity S, divided by the radius of curvature
R must be absolutely equal. Thus, therefore, not only has the character of the elastic
curve observed by the celebrated Bernoulli been most abundantly demonstrated, butalso the very great utility of my somewhat difficult formulas has been established inthis example.”83λ = 0.1

λ = 0.2

λ = 0.249

λ = 0.25

λ = 0.251

λ = 0.28

λ = 0.3027

λ = 0.35

λ = 0.4

λ = 0.5λ=1
λ=2Figure 5.11. The family of elastica solutions.

84

Figure 5.12. Bernoulli’s lemniscate.

This result certainly is confirmation of the variational technique, but in fairness it must be
pointed out that James Bernoulli was able to derive essentially the identical equation (again, notethe similarity between Equations 5.9 and 5.12) by using a combination of mechanical insight andclever integration.

5.9

Elliptic integrals

The elastica, having been present at the birth of the variational calculus, also played a majorrole in the development of another branch of mathematics: the theory of elliptic functions.
Even as the quadratures of these simple curves came to be revealed, analytic formulae for theirlengths remained elusive. The functions known by the first half of the 18th century were insufficient
to determine the length even of a curve as well-understood as an ellipse.
James Bernoulli set himself to this problem and was able to pose it succinctly (and even compute
approximate numerical values), if not fully solve it himself [81]. If the quadrature of a unit-excursion
rectangular elastica is given by the equation
Z x
y=
0

x2 dx
√
,
1 − x4(5.16)dx.
1 − x4(5.17)

then the arc length is given by:
Z x
s=

√

0

This integral is now called the “Lemniscate integral”, because of its connection with the lemniscate, another beautiful curve studied by James Bernoulli (Figure 5.12). The length of the lemniscate
is equal to that of the rectangular elastica; while Equation 5.17 gives the arc length of the elastica
as a function of the x-coordinate, this nearly identical equation relates arc length to the radial
coordinate r in the lemniscate:

85

Z r
s=

√

0

dr
.
1 − r4(5.18)

Among the lemniscate’s other representations, its implicit equation in Cartesian coordinates isa simple polynomial (no such corresponding formulation exists for the elastica),(x2 + y 2 )2 = x2 − y 2 .

(5.19)

Bernoulli approximated the integral for the arc length of the lemniscate using a series expansion
and determined upper and lower numerical bounds, but felt the calculation of them would not fall tostandard analytical techniques. He wrote, “I have heavy grounds to believe that the construction of
our curve depends neither on the quadrature nor on the rectification of any conic section.” [87]. Here,“quadrature” means the area under the curve, or the indefinite integral, and “rectification” means
the computation of the length of the curve. Bernoulli’s prediction would not turn out to be entirely
accurate; as we shall see, the curve would later be expressed in terms of Jacobi elliptic functions,which in turn are deeply related to the question of determining the arc length (rectification) of theellipse.
Fagnano took up the problem of finding the length of the lemniscate, and achieved some impressive results. Indeed, Jacobi fixes the date for the birth of elliptic functions as 23 December
1751, when Euler was asked to review Fagnano’s collected works. However, Euler had begun studyof elliptic integrals as early as 1738, when he wrote to the Bernoullis that he had “noticed a singularproperty of the rectangular elastica” having unit excursion [87]:Z 1
length · height =
0

dx
√
·
1 − x4

Z 1
0

1
x2 dx
√
= π.
4
1 − x4(5.20)

After receiving Fagnano’s work, Euler caught fire and started his remarkable research, first on
lemniscate integrals, then on the more general problem of elliptic integrals, especially the discoveryof the addition theorems for elliptic functions in the 1770s. The reader interested in more details,
including mathematical derivations, is directed to Sridharan’s delightful historical sketch [81].
These integrals would ultimately be the basis for closed-form solutions of the elastica equation, with both curvature and Cartesian coordinates given as “special functions” of the arc lengthparameter, as will be described in Section 5.12.

5.10

Laplace on the capillary – 1807

Remarkably, the elastica appears as yet another shape of the solution of a fundamental physics
problem – the capillary. Pierre Simon Laplace investigated the equation for the shape of the capillary
in his 1807 Supplément au dixième livre du Traité de mécanique céleste. Sur l’action capillaire 7 . See
I. Grattan Guinness [31, p. 442] for a thorough description of these results. In particular, Laplaceconsidered the surface of a fluid trapped between two vertical plates, and obtained the equation(722.16 in [31])7reprinted in Volume 4 of Laplace’s Complete Works [50, pp. 349–401]86z 00 = 2(αz + b−1 )(1 + z 02 )3/2 .

(5.21)

Here, z represents height, and z 0 and z 00 represent first and second derivatives with respect tothe horizontal coordinate. With suitable renaming of coordinates and substitution of Equation 5.1
for curvature in Cartesian coordinates, this equation can readily seen to be equivalent to Euler’sEquation 5.15.
Laplace also recognized that at least one instance of his equation was equivalent to the elastica.
He derives Equation 5.7 (save that Laplace writes Z for S) and noted that it is identical to theelastic curve [50, p. 379].

Figure 5.13. Maxwell’s figures for capillary action.
It is not clear when the general equivalence between the capillary surface and the elastica was
first appreciated. The special case of a single plate is published in The Elements of Hydrostatics andHydrodynamics [62, p. 32] in 1833. In any case, by 1890 James Clerk Maxwell was fully aware of it,
including it his article “Capillary Action” in the 9th edition of the Encyclopædia Britannica [57]. His
figures (reproduced here as Figure 5.13) illustrate the problem and clearly show a non-inflectional,
periodic instance of the elastica.

5.11

Kirchhoff ’s kinetic analogy – 1859

Surprisingly, the differential equation for the elastica, expressing curvature as a function of arc
length, are equivalent to those of the motion of the pendulum, as worked out by Kirchhoff in 1859
(see D’Antonio [17] and also Goss [30] for historical detail). Using simple variational techniques,taking arc length as the independent variable and the angle θ from the horizontal coordinate as thedependent, Equation 5.10 yields
θ00 + λ1 sin θ + λ2 cos θ = 0 .

87(5.22)

θr
F = mg

Fnet = mg sin θ

Figure 5.14. An oscillating pendulum.

Here, the sine and cosine arise from the specification of endpoint location constraints, which
are described as an integral of the sine and cosine of θ over the length of the curve.
This differential equation for the shape of the elastica is mathematically equivalent to that ofthe dynamics of a simple swinging pendulum. This kinetic analog is useful for developing intuition
about the solutions of the elastica equation. In particular, it’s easy to see that all solutions are
periodic, and that the family of solutions is characterized by a single parameter (modulo scaling,translation, and rotation of the system).
Consider, as shown in Figure 5.14, a weight of mass m at the end of a pendulum of length r,with angle from the vertical specified as a function of time θ(t).
Then the velocity of the mass is rθ0 (where, in this section, the 0 notation represents derivativewith respect to time), and its acceleration is rθ00 . The net force of gravity, acting with the constraintof the pendulum’s angle is mg sin θ, and thus we have
Fnet = ma = mrθ00 = mg sin θ .

(5.23)

With the substitution λ1 = −g/r, and the replacement of arc length s for time t, this equation
becomes equivalent to Equation 5.22, the equation of the elastica. Note that angle (as a functionof arc length) in the elastica corresponds to angle (as a function of time) in the pendulum, and theelastica’s curvature corresponds to the pendulum’s angular momentum.
The swinging of a pendulum is perhaps the best-known example of a periodic system. Further,
it is intuitively easy to grasp the parameter space of the system. Transformations of scaling (varying
the pendulum length) and translation (assigning different phases of the pendulum swing to t = 0)
leave the solutions essentially unchanged. Only one parameter, how high the pendulum swings,
changes the fundamental nature of the solution.
The solutions of the motion of the pendulum, as do the solutions of the elastica Equation 5.22,
form a family characterized by a single scalar parameter, once the trivial transforms of rotationand scaling are factored out. Without loss of generality, let t = 0 at the bottom of the swing(i.e. maximum velocity) and let the pendulum have unit length. The remaining parameter is thenthe ratio of the total energy of the system (which remains unchanged over time) to the potential
88

energy of the pendulum at the top of the swing (the maximum possible), in both cases counting
the potential energy at the bottom of the swing as zero. Thus, the total system energy is also equal
to the kinetic energy at the bottom of the swing.
For mathematical convenience, define the parameter λ as one fourth the top-of-swing potential energy divided by the total energy. (The constant is chosen so that λ matches the Lagrange
multiplier λ1 in Equation 5.22.)
The potential energy at the top of the pendulum is 2mgr. The kinetic energy at the bottom of
the swing is 21 mv 2 , where, as above, v = rθ0 . Thus, according to the definition above,
λ=1 2mgrg.
=
1
4 2m(rθ0 )2r(θ0 )2

In other words, assuming (without loss of generality) unit pendulum length and unit velocity
the bottom of the swing, λ simply represents the force of gravity. And this λ is the same as definedin Equation 5.13.

5.12

Closed form solution: Jacobi elliptic functions – 1880

The closed-form solutions of the elastica, worked out by Saalschütz in 1880 [76], rely heavily
on Jacobi elliptic functions. (See also [17] for more historical development, and Greenhill [33] for acontemporary presentation of these results in English, suggesting the kinetic analog was a major
motivator to deriving the closed form solutions.)
The Jacobi amplitude am(u, m) is defined as the inverse of the Jacobi elliptic integral of thefirst kind,am(u, m) = the value of φ such that
Z φ
1
p
u=
dt .
0
1 − m sin2 t(5.24)

Here, m is known as the parameter. Jacobi elliptic functions are also commonly written in terms√of the elliptic modulus k = m.
From this amplitude, the doubly periodic generalizations of the trigonometric functions follow:sn(u, m) = sin(am(u, m)) ,(5.25)

cn(u, m) = cos(am(u, m)) ,qdn(u, m) = 1 − m sin2 (am(u, m)) .

(5.26)(5.27)

Note that when m is zero, sn and cn are equivalent to sin and cos, respectively, and when m isunity, sn and cn are equivalent to tanh and sech (1/cosh) [1, §16.6, p. 571].

89Figure 5.15. Rectangular elastica as roulette of hyperbola.

5.12.1

Inflectional solutions

The closed form solutions are given separately for inflectional and non-inflectional cases. The
intrinsic form, curvature as a function of arc length, is fairly simple:√dθ= κ = 2 m cn(s, m) .
ds

(5.28)

Expressing the curve in the form of angle as a function of arc length is also straightforward:√
1
sin θ = m sn(s, m) .
2

(5.29)

Through an impressive feat of analysis, the curve can be expressed in Cartesian coordinates asa function of arc length. First, we’ll need the elliptic integral of the second kind,E(φ, k) =Z φp

1 − m sin2 θ dθ ,(5.30)

0x(s) = s − 2E(am(s, m), m) ,√y(s) = −2 m cn(s, m) .

(5.31)In Love’s presentation of these results, E(am(s, m), m) is defined in terms of an integral ofdn(u, m). The equivalence with Equation 5.30 follows through a standard identity of the Jacobi
elliptical functions ([1, p. 576], 16.26.3),
Z sE(am(s, m), m) =(dn(u, m))2 du .

(5.32)

0

The solution in terms of elliptic integrals opened the way to a few more curious results. For one,the roulette of the center of a rectangular hyperbola is a rectangular elastica [33]. More precisely, let
the hyperbola roll along a line without slipping. Then, the curve traced by its center is a rectangular90

elastica, as shown in Figure 5.15. The general term for a roulette formed from a conic section is a
Sturm’s roulette.8

5.12.2

Non-inflectional solutions

The equations in the non-inflectional case are similar.
 s

dθ
2
= κ = √ dn √ , m ,dsmm(5.33) s
1sin θ = sn √ , m .
2
m(5.34)For the Cartesian version, and for a more detailed derivation, refer to Love [52].
Truesdell considers the popularity of the elliptic function approach to be a mixed blessing formechanics and for mathematics. Where Euler’s method used direct nonlinear thinking based on
physical principles, much of the elliptic literature explored properties of special functions, which“came to be ends of research rather than means to solve a natural problem.” [87] Indeed, elliptic functions are rarely used today for computation of elastica, in favor of numerical methods.
For example, I’ve used a simple fourth-order Runge–Kutte differential equation solver (computing
Equation 5.22) to draw the figures for this work, due to its good convergence and efficiency andsimple expression in code. One particularly unappealing aspect of the elliptic approach is the sharp
split between inflectional and non-inflectional cases, while one differential equation smoothly covers
both cases.
Even so, Jacobi elliptic functions are now part of the mainstream of special functions, and fastalgorithms for computing them are well known [71]. Elliptic functions are still the method of choicefor the fastest computation of the shape of an elastica.

5.13

Max Born – 1906

In spite of the equation for the general elastica being published as early as 1695, the curves had
not been accurately plotted until Max Born’s 1906 Ph.D. thesis, “Investigations of the stability ofthe elastic line in the plane and in space” [11].
Born also constructed an experimental apparatus using weights and dials to place the elastic
band in different endpoint conditions, and used photographs to compare the equations to actual
shapes. An example setup is shown in Figure 5.16, which shows an inflectional elastica under fairly
high tension, λ ≈ 0.26.
Years later, Born wrote, “...I felt for the first time the delight of finding a theory in agreementwith measurement – one of the most enjoyable experiences I know.” [12, p. 21]. Born used the bestmodern mathematical techniques to address the problem of the elastica, and, among other things,
was able to generalize it to the three-dimensional case of a wire in space, not confined to the plane.
8

Visit http://www.mathcurve.com/courbes2d/sturm/sturm.shtml for an animated demonstration of this property.

91

Figure 5.16. Born’s experimental apparatus for measuring the elastica.

5.14

Influence on modern spline theory – 1946 to 1965

Mechanical splines made of wood or metal have long been an inspiration for the mathematicalconcept of spline (and for its name). Schoenberg’s justification for cubic splines in 1946 was a direct
appeal to the notion of an elastic strip. Schoenberg’s main contribution was to define his spline interms of a variational problem closely approximating Equation 5.10, but making the small-deflection
approximation. 9
3.1 Polynomial spline curves of order k. A spline is a simple mechanical device fordrawing smooth curves. It is a slender flexible bar made of wood or some other elastic
material. The spline is placed on the sheet of graph paper and held in place at various
points by means of certain heavy objects (called “dogs” or “rats”) such as to take theshape of the curve we wish to draw. Let us assume that the spline is so placed andsupported as to take the shape of a curve which is nearly parallel to the x-axis. If wedenote by y = y(x) the equation of this curve then we may neglect its small slope y 0 ,whereby its curvature becomes
1/R = y 00 /(1 + y 02 )3/2 ≈ y 00
The elementary theory of the beam will then show that the curve y = y(x) is a polygonalline composed of cubic arcs which join continuously, with a continuous first and second
9Note: what follows is the corrected version as appears in his selected papers [78].

92derivative10 . These junction points are precisely the points where the heavy supporting
objects are placed.
While Schoenberg’s splines are excellent for fitting the values of functions (the problem ofapproximation), it was clear that they were not ideal for representing shapes. Birkhoff and de Boor
wrote in 1965 [10] that “linearized interpolation schemes have a basic shortcoming: they are not
intrinsic geometrically because they are not invariant under rigid rotation. Physically it seems morenatural to replace linearized spline curves by non-linear splines (or “elastica”), well known amongelasticians,” citing Love [52] as an authority. They also stated the result that in a mechanical spline
constrained to pass through the control points only by pure normal forces, only the rectangularelastica is needed.
Birkhoff and de Boor also noted some shortcomings in the elastica as a replacement for themechanical spline, particularly the lack of an existence and uniqueness theory for non-linear spline
curves with given endpoints, end-slopes, and sequences of internal points. Indeed, they state that“nor does it seem particularly desirable” to have techniques for intrinsic splines approximate theelastica, and they proposed other techniques, such as Hermite interpolation by segments of Euler
spirals joined together with continuous curvature.
According to notes taken by Stanford professor George Forsythe [22], in the conference presentation of this work, they referenced Max Born’s Ph.D. thesis [11], and also described a goal of thework as providing an “automatic French curve.” The presentation must have made an impression
on Forsythe, judging from the laconic sentence, “Deep.” Forsythe, working with his colleague Erastus Lee from the Mechanical Engineering college at Stanford, would go on to analyze the nonlinear
spline considerably more deeply, deriving it from both the variational formulation of Equation 5.10
and an exploration of moments, normal forces, and longitudinal forces [51]. That publication also
included a result for energy minimization in the case where the elastica is a closed loop.

5.15

Numerical techniques – 1958 through today

The arrival of the high-speed digital computer created a strong demand for efficient algorithms
to compute the elastica, particularly to compute the shape of an idealized spline constrained topass through a sequence of control points.
Birkhoff and de Boor [10] cited several approximations to non-linear splines, including those of
Fowler and Wilson [23] and MacLaren [53] at Boeing in 1959, and described their own work along
similar lines. None of these was a particularly precise or efficient approximation.
Around the same time, and very possibly working independently, Mehlum and others designed
the Autokon system [59] using an approximation to the elastica (called the kurgla 1 algorithm,
and discussed in more detail in Section 4.9). However, this algorithm was not an accurate simulationof the mechanical spline, as it did not minimize the total bending energy of the curve, and, in fact,
could generate (approximations to) the entire family of elastica curves in service of its splines.
An objection to this technique is that the spline was not extensional; adding a new control point
coinciding with the generated spline would slightly perturb the curve. Thus, this approximate
elastica was soon replaced by the kurgla 2 algorithm implementing precisely the Euler spiral
Hermite interpolation suggested by Birkhoff and de Boor.
10

Schoenberg is indebted for this suggestion to Professor L. H. Thomas of Ohio State University.

93

A sequence of papers refining the numerical techniques for computing the elastica-based nonlinear spline followed: Glass in 1966 [26], Woodford in 1969 [90], Malcolm in 1977 [54]. Most ofthese involve discretized formulations of the problem. Interestingly, Malcolm reports that Larkindeveloped a technique based on direct evaluation of the Jacobi elliptic integrals which is nonetheless“probably quite slow” compared with the discretized approaches, “due to the large number of transcendental functions which must be evaluated.” Malcolm’s paper contains a good survey of known
numerical techniques, and is recommended to the reader interested in following this development
in more detail.
In spite of the fairly rich literature available on the elastica, at least one researcher, B. K. P.
Horn in 1981, seems to have independently derived the rectangular elastica from the principle ofminimizing the strain energy [38], going through an impressive series of derivations, and using thefull power of elliptic integral theory, to arrive at exactly the same integral as Bernoulli had derived
almost three hundred years previously.
Edwards [18] proposed a spline solution in which the elastica in the individual segments arecomputed using Jacobi elliptic functions. He also gave a global solver based on a Newton technique,each iteration of which solves a band-diagonal Jacobian matrix. A particular concern was exploring
the cases when no stable solution exists, a weakness of the elastica-based spline compared to othertechniques. Edwards claimed that his numerical techniques would always converge on a solution if
one exists.

94

Chapter 6

History of the Euler spiralThis chapter traces the history of the Euler spiral, a beautiful and useful curve known by several
other names, including “clothoid,” and “Cornu spiral.” The underlying mathematical equation is
also most commonly known as the Fresnel integral. The profusion of names reflects the fact thatthe curve has been discovered several different times, each for a completely different application:
first, as a particular problem in the theory of elastic springs; second, as a graphical computation
technique for light diffraction patterns; and third, as a railway transition spiral. Through all thenames, the curve retains its aesthetic and mathematical beauty as Euler had clearly visualized. Its
equation is related to the Gamma function, the Gauss error function (erf), and is a special case of
the confluent hypergeometric function.
The Euler spiral is defined as the curve in which the curvature increases linearly with arc length.
Changing the constant of proportionality merely scales the entire curve. Considering curvature as asigned quantity, it forms a double spiral with odd symmetry, a single inflection point at the center,as shown in Figure 4.4. According to Alfred Gray, it is “one of the most elegant of all plane curves.”
[32]

6.1

James Bernoulli poses a problem of elasticity – 1694

The first appearance of the Euler spiral is as a problem of elasticity, posed by James Bernoulli
in the same 1694 publication as his solution to a related problem, that of the elastica.
The elastica is the shape defined by an initially straight band of thin elastic material (such as
spring metal) when placed under load at its endpoints. The Euler spiral can be defined as somethingof the inverse problem; the shape of a pre-curved spring, so that when placed under load at one
endpoint, it assumes a straight line.
The problem is shown graphically in Figure 6.1. When the curve is straightened out, the moment
at any point is equal to the force F times the distance s from the force. The curvature at the pointin the original curve is proportional to the moment (according to elementary elasticity theory).

95sκ = css

M = FsF

Figure 6.1. Euler’s spiral as an elasticity problem.

Because the elastic band is assumed not to stretch, the distance from the force is equal to the arc
length. Thus, curvature is proportional to arc length, the definition of the Euler spiral.
James Bernoulli, at the end of his monumental 1694 Curvatura Laminae Elasticae (see Section
5.6), presenting the solution to the problem of the elastica, sets out a number of other problems
he feels may be addressed by the techniques set forth in that paper, for example, cases where theelastica isn’t of uniform density or thickness.
A single sentence among a list of many, poses the mechanics problem whose solution is the Euler
spiral: To find the curvature a lamina should have in order to be straightened out horizontally bya weight at one end1 .
The same year, Bernoulli wrote a note containing the integral2 entitled “To find the curve whichan attached weight bends into a straight line; that is, to construct the curve a2 = sR”.
Quia nominatis abscissa = x, applicata = y, arcu curvæ s, & posita ds constante,
radius circuli oscularis, curvedini reciproce proportionalis, est dxds : −ddy; habebitur,
ex hypothesi, hæc æquatio −aaddy = sdsdx. ...
1“Nec non qualem Lamina curvaturam habere debeat, ut ab appenso onere, vel proprio pondere, vel ab utroquesimul in rectam extedatur.” The translation here is due to Raymond Clare Archibald [4].
2
The 1694 original (Latin title “Invenire curvam, quae ab appenso pondere flectitur in rectam; h.e. construere
curvam aa = sz” is No. CCXVI of Jacob Bernoulli’s Thoughts, notes, and remarks. An expanded version was publishedin slightly expanded form as No. XX of his “Varia Posthuma,” which were collected in volume 2 of his Opera, publishedin 1744 (and available online through Google Books) [9, pp. 1084–1086]. Both works are also scheduled to be be
published in Volume 6 of Die Werke von Jacob Bernoulli.

96

Figure 6.2. Bernoulli’s construction.

97

Figure 6.3. Euler’s drawing of his spiral, from Tabula V of the Additamentum.

Translated into English (and with slightly modernized notation):
Let us call the abscissa x, the ordinate y, the arc length s, and hold ds constant. Then,the radius of the osculating circle, which is proportional to the reciprocal of the moment,
is dx ds/d2 y. Thus, we have, by hypothesis, the equation −a2 d2 y = sds dx.
The remainder of the note is a geometric construction of the curve. According to Truesdell [86,p. 109], “it is not enlightening, as it does not reveal that the curve is a spiral, nor is this indicatedby his figure.” Nonetheless, it is most definitely the equation for the Euler spiral. This construction
is illustrated in Figure 6.2. The curve of interest is ET , and the others are simply scaffolding fromthe construction. However, it is not clear from the figure that Bernoulli truly grasped the shape ofthe curve. Perhaps it is simply the fault the draftsman, but the curve ET is barely distinguishable
from a circular arc.
In summary, Bernoulli had written the equation for the curve, but did not draw its true shape,
did not compute any values numerically, and did not publish his reasoning for why the equation
was correct. His central insight was that curvature is additive; more specifically, the curvature ofan elastic band under a moment force is its curvature in an unstressed state plus the product of the
moment and a coefficient of elasticity. But he never properly published this insight. In editing his
work for publication in 1744, his nephew Nicholas I Bernoulli wrote about the equation s = −a2 /R,“I have not found this identity established” [86, p. 108].

6.2

Euler characterizes the curve – 1744

The passage introducing the Euler spiral appears in section 51 of the Additamentum3 , referring
to his Fig. 17, which is reproduced here as Figure 6.3:
3

Additamentum 1 to Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes, sive solutio problematis isoperimetrici lattissimo sensu accepti [20]. The quotes here are based on Oldfather’s 1933 translation [68],with additional monkeying by the author.

98

a
m
A
M

B

P

Figure 6.4. Reconstruction of Euler’s Fig. 17, with complete spiral superimposed.

51. Hence the figure amB, which the lamina must have in its natural state, can be
determined, so that by the force P , acting in the direction AP , it can be unfolded into
the straight line AM B. For letting AM = s, the moment of the force acting at the
point M will equal P s, and the radius of curvature at M will be infinite by hypothesis,or 1/R = 0. Now the arc am in its natural state being equal to s, and the radiusof curvature at m being taken as r, because this curve is convex to the axis AB, thequantity r must be made negative. Hence P s = Ekk/r, or rs = aa, which is the equationof the curve amB.
[20, §51]
In modern terms (and illustrated by the modern reconstruction of Euler’s drawing, Figure 6.4),
the lamina in this case is not straight in its natural (unstressed) state, as is the case in his maininvestigation of the elastica, but begins with the shape amB. At point B, the curve is held so thetangent is horizontal (i.e. point B is fixed into the wall), and a weight P is suspended from theother end of the lamina, pulling that endpoint down from point A to point a, and overall flatteningout the curve of the lamina.
The problem posed is this: what shape must the lamina amB take so that it is flattened into an
exactly straight line when the free end is pulled down by weight P ? The answer derived by Euler
appeals to the simple theory of moments: the moment at any point M along the (straightened)lamina is the force P times the distance s from A to M .
The curvature of the curve resulting from the original shape stressed by force P is equal to theoriginal curvature plus the moment P s divided by the lamina’s stiffness Ek 2 . Since this resulting
curve must be a straight line with curvature zero, the solution for the curvature of the original curveis κ = −P s(Ek 2 )−1 . Euler flips the sign for the curvature and groups all the force and elasticity
constants into one constant a for convenience, yielding 1/r = κ = s/aa.

99

From this intrinsic equation, Euler derives the curve’s quadrature (x and y as a function of thearc length parameter s), giving equations for the modern expression of the curve (formatting inthis case preserved from the original):
x=

R

ss
, &y=
ds sin. 2aaR

ss
ds cos. 2aa(6.1)

Euler then goes on to describe several properties of the curve, particularly, “Now from thefact that the radius of curvature continuously decreases the greater the arc am = s is taken, it ismanifest that the curve cannot become infinite, even if the arc s is taken infinite. Therefore thecurve will belong to the class of spirals, in such a way that after an infinite number of windings it
will roll up a certain definite point as a center, which point seems very difficult to find from thisconstruction.” [20, §52]
Euler did give a series expansion for the above integral, but in the 1744 publication is not able toanalytically determine the coordinates of this limit point, saying “Therefore analysis should gain nosmall advance, if someone were to discover a method to assign a value, even if only approximate, tothe integrals [of Equation 6.1], in the case where s is infinite; this problem does not seem unworthy
for geometers to exercise their strength.”
Euler also derived a series expansion for the integrals [20, √§53], which still remains a viable
method for computing them for reasonably small s (here b = 2a, chosen for notational convenience).

x=

6.3

s3
s7
s11
s15
−
+
−
+ &c.
1 · 3 b2 1 · 2 · 3 · 7 b6 1 · 2 · 3 · 4 · 5 · 11 b10 1 · 2 · ... · 7 · 15 b14
s9
s13
s5
+
−
+ &c.
y =s−
1 · 2 · 5 b4 1 · 2 · 3 · 4 · 9 b8 1 · 2 · ... · 6 · 13 b12(6.2)

Euler finds the limits – 1781

It took him about thirty-eight years to solve the problem of the integral’s limits. In his 1781 “Onthe values of integrals extended from the variable term x = 0 up to x = ∞”4 , he finally gave thesolution, whichphe had “recently found by a happy chance and in an exceedingly peculiar manner”,of x = y = √a2 π2 . This limit is marked with a cross in Figure 4.4, as is its mirror-symmetric twin.
The paper references a “Fig. 2”. Unfortunately, the original figure is not easy to track down. The
version appearing in Figure 6.5 is from the 1933 edition of Euler’s collected works. For reference,an accurately plotted reconstruction is shown alongside.
Euler’s technique in finding the limits is to substitute s2 /2a2 = v, resulting in these equivalences(this part of the derivation had already been done in his 1744 Additamentum [20, §54]):
4“De valoribus integralium a terminus variabilis x = 0 usque ad x = ∞ extensorum”, presented to the Academy
at Petrograd on April 30, 1781, E675 in the Enneström index. Jordan Bell has recently translated this paper into
English [21].

1002. Obtulerunt se mihi autem quondam aliae huiusmodi formulae efunctiones transcendentes
involventes, quarum valores desiderati omnes
thodos adhuc cognitas respuere videantur.
Quaesiveram scilicet eam lin
curvam, in qua radius osculi ubique recip
esset proportionalis arcui curvae, ita ut p
arcu = s et radio osculi = r essetO

rs == aa.
M

Hinc enim haud difficile est figuram cu
libero quasi manus ductu describere, quaACPquidem ea talem habere debet figuram (FiInitio nimirum curvae in A constituto
curva continuo magis incurvabitur et tan
Figure 6.5. The figure
from Euler’s “De valoribus integralium”, with modern reconstruction.
post infinitas spiras in certum punctum 0 glomerabitur, quod polum hcurvae appellare licebit. Propositum igitur mihi fuerat locum huius
AC etcoordinatarum
accuratiusquantitatem
investigare
pro Z eoque
Z ∞
∞
2
x
a
sin v
√ dv ,
sin 2 dx = √perscrutari.
2a
v
2 0
Z ∞
a
cos v
√ dv .
cos 2 dx = √
2a
v in
2 0
introducta
Hune 0 in finem
0

Z ∞

3.

(6.3)

x2

calculum

portionis

cuiusvisAM= considersr = g- function
fit 5 (which he calls ∆), defined as:(p, ut sit
To solve these
integrals, Euler
the Gamma
amplitudine
Z ∞Γ(z) =

hincque ss = 2aacp ettz−1 e−t dt

sds = aad(p

(6.4)

0

Through some manipulation, Euler solves a pair ofs fairly
general=integrals,= aV2(p2cV(p. of which the limit
point of the Euler spiral will be a special case. Assuming p = r cos α, q = r sin α, he derives

unde posita abscissa pro hoc arcu AP = x et
Γ(z) cos zαt e· cos qt dt
=
,
fore
rz
P M 0= y colligitur

Hinc iam proditZ ∞ds =z−1 −px

plicataZ ∞tz−1 e−px · sin qt dt =0Γ(z) sin zα.
rz

Euler sets q = 1, p = 0, z = 21 and derives the limit of √a2√
from the well-known value Γ( 12 ) = π.

(6.5)

pπ

2 , which follows straightforwardly

4. Hinc ergo pro polo 0 determinando requiruntur valores harum dua
fu
formularum integralium, postquam a termino <p = 0 usque ad cp oo
6.4 Relation
to theInitio
elastica
extensae.
quidem sum arbitratus hos valores aliter obtineri non
nisi approximando,
dum utraque formula successive per partes evolv
The Euler spiral can be considered something of a cousin to the elastica. Both curves weread James
deindewas aresponsible
sciliceta rp =problems.
0 usqueIn fact,<p= n forusque ad <p
cp = n,Bernoulli
primo
initially describedin terms
of elasticity

5
Euler published his discovery of the Gamma function in 1729, as De progressionibus transcendentibus seu quarum
termini generales algebraice dari nequeunt, Enneström index E19.

101

rectangular elastica
κ=xcubic parabolay’’ = xEuler spiral
κ=sFigure 6.6. Euler’s spiral, rectangular elastica, and cubic parabola.

posing both problems, and Leonhard Euler described both curves in detail about fifty years later,in his 1744 Additamentum [20, p. 276].
Both curves also appear together frequently in the literature on spline curves, starting from
Birkhoff and de Boor’s 1965 survey of nonlinear splines [10], through Mehlum’s work on the Autokon
system [59], and Horn’s independent derivation of the rectangular elastica (the “curve of least
energy”) [38].
The close relationship between the curves is also apparent in their mathematical formulations.
The simplest equation of the elastica is κ = x, while that of the Euler spiral is κ = s (here,κ represents curvature, x is a Cartesian coordinate, and s is the arc length of the curve). Thissimilarity of equation is reflected in the similarity of shape, especially in the region where the arc
is roughly parallel to the x axis, as can be seen in Figure 6.6, which shows both the rectangularelastica and Euler spiral, as well as the “cubic parabola,” the curve that results under the verysmall angle approximation κ ≈ y 00 .
Similarly, both curves can be expressed in terms of minimizing a functional. The elastica is thecurve that minimizesZE[κ(s)] =102

κ2 ds .

(6.6)

The Euler spiral is one of many solutions that minimizes the L2 -norm of the variation ofcurvature (known as the MVC, or minimum variation curve) [65]. It is, in fact, the optimal solution
when the curvatures (but not the endpoint angles) are constrained.
E[κ(s)] =

Z 

dκ 2
ds .
ds

(6.7)

It is also the curve that minimizes the L∞ norm of the curvature variation, subject to endpoint
constraints.

6.5

Fresnel on diffraction problems – 1818

Around 1818, Augustin Fresnel considered a problem of light diffracting through a slit, andindependently derived integrals equivalent to those defining the Euler spiral6 . At the time, heseemed to be unaware of the fact that Euler (and Bernoulli) had already considered these integrals,
or that they were related to a problem of elastic springs. Later, this correspondence was recognized,
as well as the fact that the curves could be used as a graphical computation method for diffraction
patterns.
φ

s
s1

x

s0
λ

Figure 6.7. Diffraction through a slit.

The following presentation of Fresnel’s results loosely follows Preston’s 1901 The Theory ofLight [72], which is among the earliest English-language accounts of the theory. Another readable
account is Houstoun’s 1915 Treatise on Light [40].
Consider a monochromatic light source diffracted through a slit. Based on fundamental principles of wave optics, the wavefront emerging from the slit is the integral of point sources at eachpoint along the slit, shown in Figure 6.7 as s0 through s1 . Assuming the wavelength is λ, the phase
φ of the light emanating from point s reaching the target on the right is
6

A. J. Fresnel, Mémoire sur la diffraction de la lumière, Annales de Chimie et de Physique, t.x., 288 1819. Henry
Crew published a translation into English in 1900 [25]

103

φ=

2π p 2x + s2 .
λ(6.8)

Assuming that s  x, apply the simplifying approximationφ≈
2π
x + 12 s2 .
λ

(6.9)

Again assuming s  x, the intensity of the wave can be considered constant for all s0 < s < s1 .
Dropping the term including x (it represents the phase of the light incident on the target, butdoesn’t affect total intensity), and choosing units arbitrarily to simplify constants, assume λ = 21 ,and then the intensity incident on the target isI=

h Z s1s0Zi2i2 h Z s1i2 h Z s1ππ 2 i2 h s1
sin s2 ds .
cos s ds +sin φ(s) ds =cos φ ds +
2
2
s0
s0
s0

(6.10)The indefinite integrals needed to compute this intensity are best known as the Fresnel integrals:
Z z

 πt2 dt ,2
Z z πt2 C(z) =dt .
cos
20
S(z) =

sin0(6.11)√
Choosing a = 1/ π, these integrals are obviously equivalent to the formula for the Cartesian
coordinates of the Euler spiral, Equation 6.1. The choice of scale factor gives a simpler limit:
S(z) = C(z) = 0.5 as z → ∞. Fresnel gives these limits, but does not justify the result [25, p. 124].
Given these integrals, the formula for intensity, Equation 6.10, can be rewritten simply as
2
2
I = S(s1 ) − S(s0 ) + C(s1 ) − C(s0 ) .

(6.12)

Fresnel included in his 1818 publication a table of fifty values (with equally spaced s) to fourdecimal places.
Alfred Marie Cornu plotted the spiral accurately in 1874 [15] and proposed its use as a graphicalcomputation technique for diffraction problems. His main insight is that the intensity I is simply
the square of the Euclidean distance between the two points on the Euler spiral at arc length s0 ands1 . Cornu observes the same principle as Bernoulli’s proposal of the formula for the integrals: Le
rayon de courbure est en raison inverse de l’arc (the radius of curvature is inversely proportional
to arc length), but he, like Fresnel, also seems unaware of Euler’s prior investigation of the integral,
or of the curve.
Today, it is common to use complex numbers to obtain a more concise formulation of the Fresnel
integrals, reflecting the intuitive understanding of the propagation of light as a complex-valued wave
Z z
C(z) + iS(z) =
0

104

π 2

ei 2 t dt .

(6.13)

Figure 6.8. Cornu’s plot of the Fresnel integrals.

Even though Euler anticipated the important mathematical results, the phrase “spiral of Cornu”
became popular. At the funeral of Alfred Cornu on April 16, 1902, Henri Poincaré had these glowing
words: “Also, when addressing the study of diffraction, he had quickly replaced an unpleasantmultitude of hairy integral formulas with a single harmonious figure, that the eye follows withpleasure and where the spirit moves without effort.” Elaborating further in his sketch of Cornu in
his 1910 Savants et écrivains, “ Today, everyone, to predict the effect of an arbitrary screen on abeam of light, makes use of the spiral of Cornu.” (Translations from the original French are mine.)
Since apparently two names were not adequate, Ernesto Cesàro around 1886 dubbed the curve“clothoide”, after Clotho (Kλωθώ), the youngest of the three Fates of Greek mythology, who spun
the threads of life, winding them around her distaff – since the curve spins or twists about its
asymptotic points. Today, judging from the number of documents retrieved by keyword from anInternet search engine, the term “clothoid” is by far the most popular7 . However, as Archibald
wrote in 1917 [4], by modern standards of attribution, it is clear that the proper name for thisbeautiful curve is the Euler spiral, and that is the name used here throughout.
7

As of 25 Aug 2008, Google search reports 17,700 results for “clothoid”, 5,660 hits for “Cornu spiral”, 935 hits for“Euler spiral”, and 17,800 for “Fresnel integral”.

105

6.6

Talbot’s railway transition spiral – 1890

The third completely independent discovery of the Euler spiral is in the context of designingrailway tracks to provide a smooth riding experience. Over the course of the 19th century, the needfor a track shape with gradually varying curvature became clear. William Rankine, in his Manualof Civil Engineering [73, p. 651], gives Mr. William Gravatt credit for the first such curve, about
1828 or 1829, based on a sine curve. The elastica makes another appearance, in a proposal about
1842 by William Froude to use circles for most of the curve, but ”a curve approximating the elastic
curve, for the purpose of making the change of curvature by degrees.”Charles Crandall [16, p. 1] gives priority for the “true transition curve” to Ellis Holbrook,in the Railroad Gazette, Dec. 3, 1880. Arthur Talbot was also among the first to approach the
problem mathematically, and derived exactly the same integrals as Bernoulli and Fresnel before
him. His introduction to “The Railway Transition Spiral” [84] describes the problem and his solution
articulately:
A transition curve, or easement curve, as it is sometimes called, is a curve of varying
radius used to connect circular curves with tangents for the purpose of avoiding theshock and disagreeable lurch of trains, due to the instant change of direction and also
to the sudden change from level to inclined track. The primary object of the transition
curve, then, is to effect smooth riding when the train is entering or leaving a curve.
The generally accepted requirement for a proper transition curve is that the degree-ofcurve shall increase gradually and uniformly from the point of tangent until the degreeof the main curve is reached, and that the super-elevation8 shall increase uniformly
from zero at the tangent to the full amount at the connection with the main curve andyet have at any point the appropriate super-elevation for the curvature. In addition tothis, an acceptable transition curve must be so simple that the field work may be easily
and rapidly done, and should be so flexible that it may be adjusted to meet the varied
requirements of problems in location and construction.
Without attempting to show the necessity or the utility of transition curves, this paper
will consider the principles and some of the applications of one of the best of thesecurves, the railway transition spiral.
The Transition Spiral is a curve whose degree-of-curve increases directly as the distancealong the curve from the point of spiral.
Thus, we have yet another concise statement of what Bernoulli and Euler wrote as rs = aa.
In his introductory figure (here reproduced as Figure 6.9), Talbot shows the spiral connecting
a straight section tangent to point A to a circular arc LH (of which DLH continues the arc, andC is the center of the circle). The remainder of the paper consists of many tables of values for thisspiral, as well as examples of its application to specific problems.
Talbot derives the basic equations in terms of the angle of direction (he uses ∆, apparently
standing for “degrees”, having already used up θ to represent the angle BAL; also, L is arc length,
not to be confused with point L in the figure): ∆ = 21 aL2 [84, p. 9]. He then expresses the x and y
8Super-elevation is the difference in elevation between the outer and inner rails, banking a moving train to reduce
the lateral acceleration felt by passengers.

106Figure 6.9. Talbot’s Railway Transition Spiral.

coordinates in terms of the simple equations dy = ds sin ∆ and dx = ds cos ∆, quickly moving toa series expansion for easy numerical evaluation, then integrating. In particular, Talbot writes [84,p.10],
y = .291 a L3 − .00000158 a3 L7 ...
x = 100 L − .000762 a2 L5 + .0000000027 a4 L9 ...

(6.14)Aside from a switch of x and y, and a change of constants representing actual physical units
used in railroad engineering, these equations are effectively identical to those derived in Euler’s
1744 Additamentum, see Equation 6.2.
Obviously, Talbot was unaware of Euler’s original work, or of the identity of his spiral tothe original problem in elasticity. Similarly, Archibald did not cite any of the railroad work in
his otherwise extremely complete 1917 survey [4]. The earliest published connection between therailway transition spiral and the clothoid I could find is in a 1922 book by Arthur Higgins [36],where it is also referred to as the “Glover’s spiral”, a reference to a 1900 derivation by James Glover[27] of results similar to Talbot’s, but with considerably more mathematical notation.

6.7

Mathematical properties

Euler’s spiral and the Fresnel integrals that generate it have a number of interesting mathematical properties. This section presents a selection.
The Fresnel integrals are closely related to the error function [1, p.301], §7.3.22:
C(z) + iS(z) =√
1+i  πerf
(1 − i)z22(6.15)The Fresnel integrals can also be considered a special case of the confluent hypergeometricfunction [1, p. 301], §7.3.25:

107

C(z) + iS(z) = z 1 F1

1 3 π 22; 2; i 2 z.

(6.16)

The confluent hypergeometric function is∞

X (a)k z kaa(a + 1) z 2+ ... =,
1 F1 (a; b; z) = 1 + z +bb(b + 1) 2!
(b)k k!

(6.17)

k=0

where (a)k and (b)k are Pochhammer symbols, also known as the “rising factorial,”(a)k =(a + k − 1)!
Γ(x + k)=.
(a − 1)!
Γ(x)(6.18)

Due to these equivalences, the Fresnel integrals are often considered part of a larger family of
related “special functions” containing the Gamma function, erf, Bessel functions, and others.

6.8

Use as an interpolating spline

Several systems and publications have proposed the use of the Euler spiral as an interpolating
spline primitive, but for some reason it has not become popular. Perhaps one reason is that many
authors seem to consider the Euler spiral nothing but an approximation to a “true” spline basedon the elastica.
An early reference to the Euler spiral for use as an interpolating spline is Birkhoff and de Boor’s
1965 survey of both linear and nonlinear splines [10]. There, they present the rectangular elastica as
the curve simulating the mechanical spline, but also point out that an exact simulation of this curve
is not “particularly desirable.” As an alternative, they suggest that “they approximate equally wellto Hermite interpolation by segments of Euler’s spirals” [10, p. 172], and cite Archibald [4] as theirsource.
Probably the earliest actual application of Euler’s spiral for splines was in the Autokon system,
developed in the beginning of the 1960’s. According to Even Mehlum [60], early versions of Autokon
used the kurgla 1 algorithm, which was a numerical approximation to the elastica. Later versions
used Euler’s spiral, having derived it as a small curvature approximation to the elastica equation.
(Note that kurgla 1 is based on the elastica, but not quite the same as a true minimal energy curve(MEC), as the curve segments can take elastica forms other than the rectangular elastica). Mehlum
writes: “There is hardly any visible difference between the curves the two algorithms produce, but
kurgla 2 is preferred because of its better computational economy. There is also a question ofstability in the kurgla 1 version if the total arc length is not kept under control. This question
disappears in the case of kurgla 2.”
Mehlum published both kurgla algorithms in 1974 [59]. It is fairly clear that he considered
the Euler spiral an approximation to true splines based on the elastica. He writes, “In Section 5
we make a ‘mathematical approximation’ in addition to the numerical, which makes the resulting
curves of Section 5 slightly different from those of Section 4 [which is based on the elastica as aprimitive]. The difference is, however, not visible in practical applications.”

108

The 1974 publication notes that the spline solution with linear variation of curvature andG2 -continuity is a piecewise Cornu spiral, and that the Fresnel integrals represent the Cartesian
coordinates for this curve, but does not cite sources for these facts.
Mehlum’s numerical techniques are based on approximating the Euler spiral using a sequenceof stepped circular arcs of linearly increasing curvature.
Stoer [82] writes that the Euler spiral spline can be considered an approximation to the “true”
problem of finding a minimal energy spline. His results are broadly similar to Mehlum’s, but hepresents his algorithms in considerably more detail, including a detailed construction of a banddiagonal Jacobian matrix. He also presents an application of Euler spirals as a smoothing (approximating), rather than interpolating spline.
Stoer also brings up the point that there may be many discrete solutions to the interpolating
spline problem, each G2 continuous and piecewise Euler spiral, based on higher winding numbers.
From these, he chooses the solution minimizing the total bending energy as the “best” one.
Later, Coope [14] again presents the “spiral spline” as a good approximation to the minimum
energy curve, gives a good Newton approximation method with band-diagonal matrices for globallysolving the splines, and notes that “for sensible end conditions and appropriately calculated chordangles convergence always occurs.” A stronger convergence result is posed as an open problem.

6.9

Efficient computation

Throughout the history of the Euler spiral, a major focus of mathematical investigation is tocompute values of the integral efficiently. Fresnel, in particular, devotes many pages to approximate
formulas for the definite integral when the limits of integration are close together [25], and used
these techniques to produce his table. Many subsequent researchers throughout the 19th century
refined these techniques and the resulting tables, including Knochenhauer, Cauchy, Gilbert, Peters,
Ignatowsky, Lommel and Peters (see [4] for more detail and citations on these results).
Today, the problem can be considered solved. At least two published algorithms provide accurate
results in time comparable to that needed to evaluate ordinary trigonometric functions.
One technique is that of the Cephes library [66], which uses a highly numerical technique,splitting the function into two ranges. For values of s < 1.6, a simple rational Chebyshev polynomial
gives the answer with high precision. For the section turning around the limit point, a similar
polynomial perturbs this approximation (valid for z  1) into an accurate value,C(z) + iS(z) ≈i i π z2
1+i
−
e 2 .
2
πz

(6.19)

Numerical Recipes [71] uses a similar technique to split the range. For the section s < 1.5
near the inflection point, the series given by Euler in 1744 (Equation 6.2) accurately computes thefunction. For the section spiraling around the limit point, the recipe calls for a continued fraction
2
based on the erf function. Thanks to equivalences involving eiz , this solution also converges toEquation 6.19 as z → ∞ [71, p. 255].
Thus, the Fresnel integrals can be considered an ordinary “special function”, and the coordinatesof the Euler spiral can be efficiently computed without fear of requiring significant resources.
109

Chapter 7

Four-parameter curves
For the basic problem of interpolating splines, two-parameter splines are ideal. However, when
greater than G2 -continuity is needed, or for designing particular shapes with straight-to-curvetransitions, a four-parameter family of primitive curves is appropriate. Where a curve from a twoparameter family is determined by its tangent angles at the endpoints, a four-parameter curve isdetermined by both tangent angles and curvatures.
Section 7.1 motivates the need for four parameters in more detail. Section 7.2 discusses a wellknown four-parameter spline, the Minimum Variation Curve, and presents a general equation forthe shape of the MVC primitive curve. Section 7.3 presents an approximation to the MVC, where
curvature is defined as a cubic polynomial of arc length. This spline has generally similar properties,
but is computationally much more tractable and easier to analyze. Finally, Section 7.4 presents atechnique for applying labels to control points representing different constraints, generalizing theidea of blended corner curves with G2 -continuity to a richer set of configurations.

7.1

Why four parameters?

Chapter 4 contained a strong argument in favor of a two-parameter family of curves for interpolating splines. Why, then, consider a four-parameter family? When are two parameters notadequate?
First, while the human visual system appears to be sensitive only to G2 -continuity, there areapplications for which higher orders of continuity are important. For example, “reflection lines”
resulting from the specular reflection of lines from a 3D surface show generally exhibit only C k−1
continuity when the underlying surface is Gk (see [39, p. 570] for a detailed discussion). Thus, toguarantee smooth reflection lines, the surface must be designed with a high order of continuity.
Another application demanding higher order continuity is the design of rail layouts for high-speed
trains. At the slower speeds characteristic of trains in the United States, G2 -continuity is adequate,
and clothoidal segments are the usual design norm, but in the rest of the world there has been
recent attention on smoother curves [45].

110

κ

splineFigure 7.1. Eccentricity 0.56199 oval with G2 and G4 continuity.

111

splineκ
Figure 7.2. The suitcase corners problem.

Figure 7.1 shows an example of a simple closed curve (an ovoid shape with eccentricity 0.56199)
drawn using both a G2 -continuous spline and a G4 -continuous spline. It is clear from the curvatureplots that the curvature variation is smoother in the G4 curve, but the points of discontinuity ofcurvature variation are not directly visible to the eye.
Second, many shapes are not best represented as a simple interpolating spline. The outline ofa letterform is not, in general, a smooth undulating curve (although such examples do exist, suchas the lowercase ‘c’ in Figure 4.6). In addition to sections of smooth curve, a letter form may also
have sharp corners, sections of straight line, and, most importantly for this discussion, smoothtransitions from straight lines to curved sections. A simple example is shown in Figure 7.2 (basedon Fig. 4.10 in Moreton’s Ph.D. thesis [64]), illustrating the problem of blending smooth corners
between sections of straight line.
In font design, similar situations arise for bracketed serifs, and for shapes such as the letter U,
where straight sections blend smoothly with curves. A typical example showing both situations isshown in Figure 7.3, a capital U from the author’s font Cecco. The entire curve (with the exceptionof the sharp corners) is G2 continuous, providing pleasing and smooth joins between the straightand curved sections.
A four-parameter curve family provides enough flexibility to meet such needs. It’s clear that asingle section of Euler spiral is not adequate for one of the corners in Figure 7.2, as the curvaturemust start out zero, increase to a value large enough for the segment to turn by a right angle, thendecrease again to zero.
Third, by annotating curves with additional constraints, it’s possible to improve locality. In particular, constraints that propagate curvature derivatives in one direction, but leave the parameters
unconstrained on one side (which we call “one-way constraints”) can isolate sections of extremechanges of curvature, so that those sections do not influence nearby regions of gentle curvaturechange.

112

Figure 7.3. Bracketed serifs and straight-to-curve joins.

For example, Figure 7.4 shows how adding one-way constraints keeps the variation of curvature
localized to just two segments of the overall spline, as opposed to the infinite extent of a pureinterpolating spline. The central point can be perturbed arbitrarily without affecting the segments
beyond the one-way constraints. In this figure, the top plots show curvature as a function of arc
length, and the bottom plots show the control points (with symbols representing the constraints)
with the interpolating curve. Note also that the curve with better locality also requires considerably
more curvature (and variation of curvature) to fit the points.
There are two examples of a four-parameter primitive curve, both straightforward generalizations of two-parameter counterparts. One is the Minimal Variation Curve, which minimizes the L2
norm of curvature variation. Another is the generalization of the Euler spiral so that curvature is anarbitrary cubic polynomial in arc length (the definition of the Euler spiral is a linear relationship).
Both of these splines are G4 -continuous and use a primitive curve which has four parameters. InκsplineFigure 7.4. One-way constraints enforce locality.

113

practice, these curves are quite similar, and one can easily be considered an approximation of theother.
The parameter space of a two-parameter spline can be understood as the tangent angles of theendpoints, relative to the chord. Given arbitrary endpoint tangents, a two-parameter curve familyfinds a smooth curve meeting those constraints. Similarly, a four-parameter curve family finds asmooth curve meeting given tangent and curvature values at the endpoints.

7.2

Minimum Variation Curve

The Minimum Variation Curve (MVC) is defined as the curve minimizing the L2 norm of thevariation in curvature. Moreton’s Ph.D. thesis [64] argued in favor of this spline to address known
limitations of the MEC. In particular, the MEC lacks roundness, but since a circular arc has zerocurvature variation, it is trivially the curve that minimizes the MVC cost functional when the input
points are co-circular.
Z l
EMVC [κ(s)] =(κ0 )2 ds .

(7.1)

0

Most presentations of the MVC, including Moreton’s, use a numerical approach, approximating the curve with some other approximation (in Moreton’s case, quintic polynomials) and using
a technique such as conjugate gradient descent to minimize the cost functional. Gumhold’s more
recent approach [34] is similar, also using conjugate gradient descent. However, a variational approach analogous to that of the MEC produces exact differential equations for the curve. The
simple Euler–Lagrange equation is inadequate for this task, as it only handles functionals written
in terms of a first derivative. Since the MVC functional is in terms of the second derivative of θ,
the more general Euler–Poisson equation is needed. Section 7.2.1 recaps the general technique fromthe calculus of variations, then Section 7.2.2 applies it to the MVC and gives the solution. One
formulation is the impressively simple κ00 = λy, which is pleasingly analogous to κ = λy as thegeneral solution to the elastica (Equation 3.13).

7.2.1

Euler–Poisson equation

Here we repeat the formulation of the general Euler–Poisson equation as given by Elsgolts [19]:Let us investigate the extreme value of the functionalv[y(x)] =Z x1
F (x, y(x), y 0 (x), . . . , y (n) (x)) dx ,

(7.2)

x0

where we consider the function F differentiable n + 2 times with respect to all arguments and weassume that boundary conditions are of the form
(n−1)

y(x0 ) = y0 , y 0 (x0 ) = y00 , . . . , y (n−1) (x0 ) = y0

;

(7.3)

(n−1)
y(x1 ) = y1 , y (x1 ) = y10 , . . . , y (n−1) (x1 ) = y1
,

(7.4)

0

114

i.e. at the boundary points are given the values not only of the function but also of its derivatives
up to the order n − 1 inclusive.
The function y = y(x) that extremizes the functional given in Equation 7.2 must be a solutionof the equation
n ∂F
d ∂F
d2 ∂F
∂F
n d
=0
−
+
+
.
.
.
+
(−1)
∂y
dx ∂y 0 dx2 ∂y 00
dxn ∂y (n)

7.2.2

(7.5)

Euler–Poisson solution of MVC

We add Lagrange multipliers to represent the constraint that the curve spans the endpoints.
Adding the Lagrange multipliers and eliminating κ in favor of θ, we seek to minimize
Z lE=(θ00 )2 + λ1 sin θ + λ2 cos θ ds .

(7.6)

0

To apply this equation to the MVC, we substitute x = s and y(x) = θ(s), so that
F (x, y, y 0 , y 00 ) = (y 00 )2 + λ1 sin y + λ2 cos y .

(7.7)

The relevant partial derivatives are∂F= λ1 cos y − λ2 sin y ,
∂y
∂F
=0,
∂y 0
∂F
= 2y 00 ,
∂y 00
d2 ∂F
= 2y 0000 .
dx2 ∂y 00

(7.8)
(7.9)
(7.10)(7.11)

Applying 7.5 and reversing the variable substitution, we derive an ordinary differential equation
for the MVC,
λ1 cos θ − λ2 sin θ + 2θ0000 = 0 .

(7.12)

Because at this point we are concerned only about the existence of the Lagrange multipliers,
not their actual values, we can ignore the constant factor of 21 and − 12 for λ1 and λ2 , respectively.
θ0000 + λ1 cos θ + λ2 sin θ = 0 .

(7.13)

Note the pleasing similarity to the elastica (MEC) equation,
θ00 + λ1 cos θ + λ2 sin θ = 0 .
115(7.14)

As a special case when λ1 and λ2 become zero, all second order polynomial spirals κ(s) =k0 + k1 s + k2 s2 are solutions to the MVC equation, as the the third derivative of κ (the fourth
derivative of θ) vanishes. These solutions are analogous to circular arcs κ(s) = k0 as solutions of
the elastica when the corresponding Lagrange multipliers become zero.
Again as in the MEC case, we can reformulate this equation solely in terms of first and higher
derivatives of θ, also eliminating the Lagrange multipliers.
First, we differentiate Equation 7.13, giving
θ(5) + θ0 (−λ1 sin θ + λ2 cos θ) .

(7.15)

To integrate Equation 7.13, we first multiply by θ0 , giving
θ0000 θ0 + θ0 (λ1 cos θ + λ2 sin θ) = 0 .

(7.16)

This equation can readily be integrated, giving
θ000 θ0 −(θ00 )2+ (λ1 sin θ − λ2 cos θ) + C = 0 .
2(7.17)

Combining with Equation 7.15,
θ(5) + θ000 (θ0 )2 − θ0 (θ00 )2 /2 + Cθ0 = 0 .

(7.18)

Rewriting in terms of κ,
κ0000 + κ00 κ2 −κ(κ0 )2+ Cκ = 0 .
2(7.19)

When the length is not constrained, the constant C becomes zero. This follows from transversality, a standard technique in the calculus of variations [19, p. 345], and is analogous to the way
the corresponding constant becomes zero in the case where the MEC is not length-constrained.
Then, the equation becomes
κ0000 + κ00 κ2 − 21 κ(κ0 )2 = 0 .

(7.20)

Equation 7.20 was also given by Brook et al [13]; their solution does not encompass additional length constraints (nonzero values of C). This equation is well suited to standard numerical
techniques such as Runge–Kutte integration for very precise computation of the shape of the curve.
The MVC spline is defined as the curve that globally minimizes the MVC metric over the spaceof all curves that interpolate the input points. As in the MEC, this curve is made of segments ofa parametrized primitive curve, joined together at the control points with continuity constraints.
The MVC primitive curve (defined by Equation 7.20) is four parameter, and it is G4 -continuous at
the control points, as opposed to two parameters and G2 for the MEC.
At the endpoints of the MVC, κ0 = 0 and κ00 = 0 (again, analogous to the endpoint conditionof the MEC, κ = 0). It’s fairly easy to see why the latter conditions hold. At an endpoint, it is116

possible to add at least points lying on a circular arc of matching tangent and curvature. The
resulting MVC spline must assign exactly that circular arc to the additional segments, as thevalue of the MVC functional itself for the arc is zero, and if there is any other curve with a lower
value for the functional, then it would have been a better solution for the original problem. And
since the MVC is G4 -continuous, then κ0 and κ00 must match the values for a circular arc, whichare zero. These endpoint conditions can also be derived straightforwardly (if abstractly) from thetransversality condition, a principle of the calculus of variation which determines scalar parametersfrom the solution space to minimize the global functional (see a textbook such as Elsgolts [19] formore details).
There is also a simple equation for the MVC relating curvature and a Cartesian coordinate,analogous to κ = λy for the MEC case. To derive it, we first assume without loss of generality that
λ1 = 0; unlike the formulation in terms of curvature alone, the new equation will not be rotationallyinvariant. We then observe that dy
ds = sin θ. Substituting that into Equation 7.13, we get
θ0000 + λ2 y 0 = 0 .

(7.21)

Integrating (the constant of integration can be set at zero, since that corresponds to a translationof the curve), substituting κ = θ0 , and adjusting the constants so λ = −λ2 , the resulting simple
formulation isκ00 = λy .

(7.22)

When λ = 0, this equation is equivalent to that of the Euler spiral, in much the same way that
λ = 0 is the value for which the elastica becomes a circle. Kimia et al [44] erroneously gave thisas the general solution to the MVC equation, correctly observing that the Euler spiral is a special
case of the MVC, but they apparently overlooked the possibility of nonzero values of λ.

7.3

Polynomial spiral spline

Another approach to four-parameter splines is to generalize the Euler spiral to an implicit
equation where curvature is a cubic polynomial function of arc length. The general equation is:κ(s) = a + bs + cs2 + ds3 .

(7.23)

This curve family is a very good approximation to the MVC. Assuming small turning angles,
y is nearly proportional to s, so substituting Equation 7.22, κ00 ≈ c1 s + c2 . Integrating twice (andbringing about two more constants of integration) yields Equation 7.23. This approximation is
analogous to that of the Euler spiral approximating the MEC, but can be expected to be evencloser, as the nonlinear terms affect only higher-order derivatives of curvature.
As in the MVC, the corresponding spline can be defined in terms of enforcing G4 -continuity
across control points. These are four constraints per control point. Since each segment requires fourconstraints to determine the four parameters of the polynomial, for an open curve there are four
additional parameters, or two for each endpoint. One approach is to set tangent and curvature
explicitly, but a more natural approach (in the interpolating spline framework) is to set κ0 = 0 and
117κ00 = 0 at the endpoints, as in the MVC. In this respect, the cubic polynomial spline is even moreclosely related to the MVC than the Euler spiral is to the MEC, because in the latter case theendpoint conditions differ; the natural end condition is constant curvature (a circular arc), while
the natural endpoint condition for the MEC is zero curvature. In the four-parameter world, theMVC endpoint conditions make sense and there is no need to adjust them.
The cubic polynomial spiral spline was first proposed by Ohlin [67], by analogy of making thepentic (polynomial) spline nonlinear in the same way that the MEC is the nonlinear counterpartof the cubic spline. Although Ohlin did express the natural variational formulation of polynomial
splines, and cited Mehlum’s work indicating that the Euler spiral spline was a consistent (extensional) approximation to the MEC, he did not explicitly mention the MVC or state that this new
spline is a good approximation to it.
The idea of a curve defined as curvature being a cubic polynomial function of arc length goesback even further. In 1936, Bloss proposed a railroad transition spiral in which the curvature isa simple cubic polynomial with respect to arc length [45]. The Euler spiral (or clothoid, as it isnearly universally known in that application domain) is popular in railroad track design. The linearchange in curvature is desirably smooth on the interior of the segment, but the literature has also
long recognized the limitations of smoothness at the joins between clothoid segments. These higher
derivatives of curvature become more significant as train speeds increase. In particular, high speed
train tracks are almost invariably banked (or “superelevated”) so that the rail on the outside of thecurve is higher than the inside, to a degree roughly proportional to the curvature. Thus, passengers
of a train traversing the tracks experience a roll acceleration roughly equal to the second derivativeof curvature. At the join of two Euler spiral segments, this roll acceleration is an impulse. Usually,
the suspension of the train provides enough cushioning to avoid an unpleasant lurch, but curves
with better continuity properties are clearly desirable for a smooth, quiet ride. The Bloss transition
spiral is a very early example of a curve with higher continuity, to avoid such impulses.

κsplineFigure 7.5. G2 spline has better locality than G4 .
If G4 -continuity is smoother than G2 , then why not use it everywhere? Unfortunately, theadditional smoothness comes at the cost of locality. As shown in Figure 7.5, the effects of perturbing
a point ripple out almost twice as far as for the Euler spiral spline. Since, for most applications,the difference in smoothness is theoretical, rather than visible to the eye, it’s better to enjoy thesuperior locality properties of the G2 -continuous spline.

118

7.4

Generalized constraints

The motivation for four-parameter splines showed a number of specific examples of blended
corners, straight-to-curve transitions, and so on. This section presents a highly general techniquethat subsumes all these examples, and provides a broad palette of curve transitions to the designer.
A single segment of cubic polynomial spline has four free parameters. In general, these can
fit any tangent and curvature constraints at the endpoints. One such example is to fit the given
tangent constraints of adjoining straight segments, as well as zero curvature, as in the blended
corners of Figure 7.2.
The approach of this section will be to start from special cases, then progressively generalize.
First, observe that the input to the solver for Figure 7.2 must distinguish which segments are straight
and which are the curved corners. We propose the notation of annotating points with half-circles.
A segment between the two straight sides of the half-circle is a straight line. A segment between
the two rounded sides is a blended corner. For such segments, the solver fits a four-parameter spline
segment to match the tangent and curvature of the adjoining straight segment.
The next generalization is to allow more than one segment in the curved section. In fact, any
G4 -continuous spline interpolating a sequence of control points has two free parameters at theendpoints. For the most natural end conditions, these are determined by constraining κ0 = 0 andκ00 = 0, but it’s also possible to determine them by constraining tangent and curvature constraints
to specific values. In the case of the right half of Figure 7.4, there are two segments bracketed
between the half-circle sides of the one-way constraints, and the tangents and curvatures are, as
above, matched to the adjoining straight segments. One important observation is that perturbing
the central point does not affect the curve beyond the one-way constraints.

κ

splineFigure 7.6. Matching nonzero curvatures.

The next generalization is to match adjoining curvatures other than zero. While straight-tocurve transitions are common, sometimes the adjoining section will have some curvature, rather
than be straight. An example is shown in Figure 7.6. The general technique for the solver is tocompute the adjoining curve segments first, treating the flat side of the one-way constraint thesame as an endpoint, then use the tangent and curvature value at that point to constrain the spline
segment on the curved side. Generalizing further, the one-way constraint is treated as a dependency,
so that the spline on the curved side depends on the tangent and curvature of the other side, but
not vice versa. In many cases, it’s possible to solve the entire spline by doing a topological sort ofthe segments, and solving in order of dependency.

119

κ

splineFigure 7.7. A one-way constraint creates a cyclic dependency.

Figure 7.8. Examples of different kinds of constraint points.

However, even this approach is not maximally general. Allowing arbitrary annotations of control
points can create a cyclic dependency, as in Figure 7.7, which even so has a meaningful solution.
More importantly, the techniques above only allow G4 -continuous joins for interior points, while,
as we have seen, Euler spiral segments (G2 -continuous joins) have better locality properties.
In the general approach, each segment between two adjacent control points is assigned a cubic
polynomial spiral segment (i.e. as defined by Equation 7.23). Note that an Euler spiral segment fits
into this framework; the parameters c and d are both zero. Equivalently, κ00 is constrained to be
zero at both endpoints. For all types of internal points (other than corners, which are considered
endpoints) are also constrained to have G2 -continuity, in other words, both tangent angle θ andcurvature κ are constrained to be equal on both sides of a control point. Different points haveadditional constraints, depending on the type of point.
The total number of degrees of freedom (and thus, number of constraints) is four times thenumber of segments. If each endpoint has two constraints and each internal point has four, then
the total number of constraints is satisfied.
The complete list of point types, and associated mathematical constraints, is given in the tablebelow.
For a given sequence of points with constraint annotations, the resulting spline is defined as thecurve that is piecewise cubic polynomial spline, such that the tangent, curvature, and higher-orderderivatives satisfy the constraints given in the table above for each point.
Figure 7.4 shows a sequence of somewhat arbitrary constraints, parallel with a curvature plot120

Type of constraint

FormulaSymbolκ0r = 0

Left endpointκ00r = 0
κ0l = 0

Right endpointκ00l = 0G2 curve pointθl = θrκ00l = 0κl = κr κ00r = 0
G4 curve pointθl = θrκ0l = κ0r

κl = κr κ00l = κ00r
One-way (curve to terminal) θl = θr

κ0r = 0κl = κr κ00r = 0
One-way (terminal to curve) θl = θr

κ0l = 0κl = κr κ00l = 0
Table 7.1. Constraints for curve points, and their corresponding formulae.
which shows the effect of the different types of constraints in the curvature domain. All segments
between two G2 curve points have linear curvature. All segments between a G2 point and an
endpoint (or between a G2 point and the straight side of a one-way constraint) have constantcurvature, i.e. the segment is a circular arc. All of the points have at least G2 continuity, although
if corner points were added, those would have only G0 .
Section 8.3 presents numerical techniques for satisfying such a system of constraints, and Chapter 9 discusses techniques for efficiently drawing the polynomial spline curves that result. With allthese components in place, the spline is both highly flexible for drawing the types of shapes likely
to arise in font design, and also entirely practical for interactive editing.

121

Chapter 8

Numerical toolbox
One of the reasons for the enduring popularity of polynomial-based splines is the ease of computation, and straightforward numerical properties such as stability. By contrast, variational techniques are considered to be extremely slow and unwieldy. For example, Moreton’s MVC solver took
about 2.5 billion CPU cycles (75 s on the hardware of the day) to compute the MVC solution to
the 7-point data set posed by Woodford [90]. Even though more refined techniques have since been
published, the prejudice continues. Kimia et al [44] state, “The major shortcomings of the elastica
are that (i) they are not scale-invariant, (ii) they are computationally expensive to derive, and (iii)they do not lead to circular completion curves when the point-tangent pairs are, in fact, co-circular(i.e., both lie on the same circle).”This chapter presents a toolbox of numerical techniques for computing sophisticated splines
quickly, and with high precision. Section 8.1 presents a polynomial approximation to a fundamental
integral that encompasses both the Euler spiral and the general polynomial spiral used as a primitivein Chapter 7. Section 8.2 presents a very efficient solver fitting an Euler spiral segment to arbitrary
endpoint tangent constraints. Section 8.3 presents a fast Newton solver that globally satisfies thecontinuity constraints for the interpolating spline. Section 8.4 presents a fast 2-D lookup tableapproach for solving arbitrary two-parameter extensional splines using a single code framework.
The combination of these techniques, along with optimized drawing code, presented in Chapter
9, yields an interpolating spline that is extremely efficient to compute and easily provides interactive
response speeds. The complete spline is implemented and released in the open source libspiro library,which has also been integrated into FontForge, a font editing tool, and Inkscape, a vector drawing
editor.
The software release accompanying this thesis also includes a JavaScript implementation, whichruns in standard Web browsers.

122

8.1

Numerical integration of polynomial spirals

Almost all calculations involving polynomial spiral curves can be reduced to an evaluation of asingle definite integral,
Z 0.5

1

2

1

3

1

4

ei(k0 s+ 2 k1 s + 6 k2 s + 24 k3 s ) ds .

spiro(k0 , k1 , k2 , k3 ) =

(8.1)

−0.5

Here, k0 ...k3 are the curvature and its first three derivatives at the center of a curve segmentof arc length 1. More generally (for an arbitrary degree polynomial),
Z 0.5

i

spiro(k) =

e

1
j+1
0≤j (j+1)! kj s

P

ds .

(8.2)

−0.5

Fixing the limits to ±0.5 simplifies the math and (as will be seen) creates symmetries which
can be exploited to improve numerical accuracy, but does not lose generality. An integral over any
limits can be expressed in terms of the spiro integral, as follows:
Z s1

1

1

2

1

3

4

ei(k0 s+ 2 k1 s + 6 k2 s + 24 k3 s ) ds =

s0
1

2

1

3

1

4

(s1 − s0 )ei(k0 s+ 2 k1 s + 6 k2 s + 24 k3 s )spiro((s1 − s0 )(k0 + k1 s + 12 k2 s2 + 61 k3 s3 ),(s1 − s0 )2 (k1 + k2 s + 12 k3 s2 ),(8.3)(s1 − s0 )3 (k2 + k3 s),(s1 − s0 )4 k3 ),where s =s0 + s1.
2Note that, as a special case of Equation 8.3, the classic Fresnel integrals also reduce readily tothis integral.
Z s
0

2

s2

eit dt = sei 8 spiro( 2s , s2 , 0, 0) .

(8.4)

The spiro integral can be thought of as analogous to the sine and cosine functions of thegeneralized polynomial spiral. This connection is clear when considering a circular arc.
sin x
= spiro(2s, 0, 0, 0) .
x

(8.5)

The absolute value of the spiro integral is the ratio of chord length to arc length of the curve
segment (note that expressing this quantity as a ratio is invariant to scaling, and it is obviously
invariant to rotation and translation as well). For a given curve segment with chord length ∆x, thecurvature and its derivatives along the curve (with arc length s normalized to the range [−0.5, 0.5])
are given by a simple formula,123κ(j) (s) = |spiro(k)| j+1 X k
i+j
si .
∆x
(i + 1)!

(8.6)

0≤i

8.1.1

Polynomial approximation of the integral

There are a number of approaches to numerically computing an integral such as Equation 8.1.
However, a straightforward numerical integration technique such as Runge–Kutte would converge
only relatively slowly, requiring significant computation time to achieve high accuracy. The approachof this section is to compute polynomial approximations of the integral, keeping all terms up to acertain order and discarding the rest.
A key insight is that the angle θ(s) is a straightforward polynomial function of s. Thus, tointegrate eiθ(s) , take the Taylor’s series expansion of ex , substitute the polynomial iθ(s) for x, andexpand symbolically, keeping terms only up to the given order.
θ(s) = k0 s +k1 2 k2 3 k3 4
s + s + s .
2
6
24

(8.7)

Below is a worked example retaining terms for eix up to quadratic, then integrating. The
approximation is1
eix = 1 + ix − x22

+ O(x3 ) .

(8.8)

Substituting x = θ(s), we geteiθ(s) = 1 + ik0 s −ik1 − k02 2s2+ O(s3 ) .

(8.9)

Integrating, we getZeiθ(s) ds = f (s) = s +ik0 2 ik1 − k02 3
s −
s
2
6

+ O(s4 ) .

(8.10)

Evaluating f (0.5) − f (−0.5), we find that all even polynomial terms (including s4 ) drop out:f (0.5) − f (−0.5) ≈ 1 +

ik1 − k02.
24

(8.11)

Thus, as promised, fixing the limits of integration to ±0.5 provides a symmetry that bothdecreases the number of terms to compute and also increases the accuracy by one degree.
Retaining all terms up to quartic yields an approximation not much more complicated. Notethat this approximation actually uses all four parameters.
∆x ≈ 1 +ik1 − k02 ik3 − 6ik02 k1 + k04 − 4k0 k2 − 3k12
+
.
24
1920

124

(8.12)

8.1.2Higher-order polynomial approximations

Higher-order polynomials offer a good way to increase the accuracy of the approximation atmodest cost. However, manually working out the coefficients becomes quite tedious. This section
describes a technique for automatically generating code for an arbitrary polynomial degree.
For functions of a single variable, the Taylor’s series expansion gives the coefficients of theapproximating polynomial simply and unambiguously (it is, of course, possible to use refinements
such as Chebysheff polynomials to redistribute the errors more evenly, but the underlying concept
remains). Here, the task is to approximate a function of four parameters. It is possible to generalize
the Taylor’s series expansion to functions of more than one parameter, but a simplistic approach
would result in a blowup of coefficients; for a function of four parameters, the number of coefficientsgrows quadratically with the order of the polynomial.
A more sophisticated approach is represent each successive xj as a vector of polynomial coefficients. Because of the symmetries of the equation (many odd terms will drop out), for j ≥ 3, it
is most efficient to compute xj by combining the vectors for xj−2 and x2 . Then, these polynomials
are substituted into the Taylor’s series expansion for ex . A fully worked out example, for order 8,
is given below.

t11 = k0
t12 = k21
t13 = k62
t14 = k243
t22 = t211
t23 = 2t11 t12
t24 = 2t11 t13 + t212
t25 = 2(t11 t14 + t12 t13 )t26 = 2t12 t14 + t213
t27 = 2t13 t14
t28 = t214
t34 = t22 t12 + t23 t11
t36 = t22 t14 + t23 t13 + t24 t12 + t25 t11
t38 = t24 t14 + t25 t13 + t26 t12 + t27 t11
t44 = t222
t45 = 2t22 t23
t46 = 2t22 t24 + t223
t47 = 2(t22 t25 + t23 t24 )t48 = 2(t22 t26 + t23 t25 ) + t224
t56 = t44 t12 + t45 t11
t58 = t44 t14 + t45 t13 + t46 t12 + t47 t11
t66 = t44 t22
t67 = t44 t23 + t45 t22
t68 = t44 t24 + t45 t23 + t46 t22
t78 = t66 t12 + t67 t11
t88 = t66 t22
t26
t66
t68
t88
22
− t24 − 896− t28 + t44 + t46 + t48 − 322560− 1658880+ 92897280
u = 1 − t24
t12
t14 160
t34
t36 4608t38 1920 t56 10752 t5855296
t78
v = 12 + 80 − 480 − 2688 − 13824 + 53760 + 276480 − 11612160(8.13)

Obviously, for higher degrees, the formulas are too complex to be worked out by hand. The
accompanying software distribution includes a utility, numintsynth.py, that synthesizes C code
for an arbitrary degree. Asymptotically, the number of lines of code grows as O(n2 ), so for any
reasonably low degree, the code is quite concise. Additionally, the structure of the dependency

125

1
n=1
n=2
n=4

0.01

n=8
n = 16n = 320.0001n = 64n = 128n = 256

Accuracy

1e-06

1e-08

1e-10
Order 2
Order 4
Order 6
Order 8Order 10
Order 12
Order 14
Order 16

1e-12

1e-14

1e-16
1e-08

1e-07
1e-06
1e-05
Time spent, in seconds, 2.33 GHz Intel Core 2 Duo

0.0001Figure 8.1. Precision vs. timing for different integration strategies.

graph permits a high level of instruction-level parallelism, making the computation very efficienton typical CPUs.

8.1.3

Hybrid approach

The higher order the polynomial, the more accurate the result. However, code complexity grows
significantly, and for small angles, accurate results are possible with a low-order polynomial.
A hybrid approach is to implement a polynomial of fixed order, and use it for parameter values
where the error is within tolerance. For parameter values outside this tolerance, subdivide. The
range [−0.5, 0.5] is split into n equal subdivisions, each is evaluated using Equation 8.3, and theresult is the sum of all such intervals.
The extreme case of an order-2 polynomial is equivalent to Euler integration. The code is simple,
but convergence is expected to be poor.
The optimum order of polynomial is determined empirically. All even orders of polynomial
between 2 and 16 were measured, and the accuracy plotted against CPU time. The results areshown in Figure 8.1.3, measured on a MacBook Pro with a 2.33 GHz Core 2 Duo processor. Eachsuccessive datapoint doubles the total number of subdivisions. This particular plot is for evaluatingspiro(1, 2, 3, 4), but other combinations of arguments all give similar results.
126

The n = 1 case is relatively faster than n > 1, because no sine and cosine functions are needed
for the subdivision. As can be seen, low order polynomials require a large number of subdivisions
to provide accurate results. However, there is diminishing return as the order of the polynomial
grows – while the number of intervals needed is decreased, the time needed to compute each one
increases. As can be seen from the graph, orders greater than 12 don’t improve overall performance
significantly.
An empirically determined estimate for the error of the order 12 spiro approximation is:
2

1|spiro12 (k0 , k1 , k2 , k3 )−spiro(k0 , k1 , k2 , k3 )| ≈ (0.006|k0 |2 +0.03|k1 |+0.03|k2 | 3 +0.025|k3 | 2 )6 (8.14)
With n subdivisions, the error scales approximately as n−12 . Thus, it should be straightforward
to compute a value for any given precision.

8.2

Determining Euler spiral parameters

One important primitive is to compute the parameters for Euler spiral segments, given a fixed
chord length and tangents relative to the chord. This primitive is needed in the global determinationof angles to satisfy the C 2 curvature continuity constraints across knots, and is also independently
of interest in vision applications, for computing a Euler spiral segment for use in shape completion.
There are quite a number of algorithms for this in the literature, the most detailed of which is
given by Kimia et al [44]. Their solution represents the curve in terms of two parameters:
κ(s) = γs + κ0, with 0 ≤ s ≤ L(8.15)

Numerically, their solution is complicated. They first construct a biarc to approximate thedesired spiral segment, then use a gradient descent optimization method for the two parameters,
computing the Euler spiral at each iteration.
Our solution is similar in many respects, but more direct; the biarc as an initial approximation
is no longer needed. There are two key insights. First, by specifying the curvature of the curve fromthe middle, as opposed to one of the endpoints, one of the parameters has a closed-form solution.
Second, it is easier to compute Euler spiral segments of fixed arc length than fixed chord length.
Once that’s done, scaling the curve to fit the original chord length constraint does not affect theangle constraint.
We represent the curve to solve as
κ(s) = k0 + k1 s, with − 0.5 ≤ s ≤ 0.5 .

(8.16)

Integrating, and fixing θ(0) = 0,
θ(s) = k0 s +

k1 2
s .
2

Again, the curve can be given in parametric form simply by integrating the angle:
127(8.17)

ys = -.5ψθ0θ1

s = .5x

Figure 8.2. Determining Euler spiral parameters.

Z sx(s) + iy(s) =eiθ(t) dt .

(8.18)

0The endpoint angles with respect to the chord can be represented in terms of the angle of thechord in this coordinate frame.

θ0 = ψ − θ(−0.5) ,(8.19)

θ1 = θ(0.5) − ψ, wherey(0.5) − y(−0.5)
ψ = arctan.
x(0.5) − x(−0.5)(8.20)(8.21)

But the formula for θ(s) is known (Equation 8.17), so

k0 k1
−
+ψ,
2
8
k0 k1
θ1 =
+
−ψ.
2
8

θ0 =

(8.22)(8.23)

From Equations 8.22 and 8.23 we can immediately solve for k0 :k0 = θ0 + θ1 .

(8.24)We have now reduced the search space to a single dimension, k1 . We will use a simple polynomial
approximation to derive an initial guess, an approximate Newton–Raphson method to derive asecond guess, and the secant method to further refine the approximation.
Using the small angle approximation eiθ ≈ 1 + iθ, Equation 8.18 becomesZ sx(s) + iy(s) ≈1 + iθ(t) dt
0

k
k1 
0 2
s + s3 .
=s+i
2
6
128

(8.25)

Combining Equations 8.21 and 8.25 (note the cancellation of even terms, just as for the integration in the previous subsection),x(.5) − x(−.5) ≈ 1 ,k1y(.5) − y(−.5) ≈,24
k1
ψ ≈ arctan.
24(8.26)(8.27)(8.28)

Applying the small angle assumption arctan θ ≈ θ and substituting Equation 8.22,
k1
k0
−
,
2
12
k0
k1
θ1 ≈
+
.
2
12

θ0 ≈

(8.29)(8.30)

Thus, we can derive a good approximation to k1 :k1 ≈ 6(θ1 − θ0 ) .

(8.31)

This approximation may be sufficient for low-precision applications. A considerably more accurate technique is the secant method, which excels at finding the roots of nearly-linear onedimensional functions. In this application, the function to be solved is angle error as a function of
k1 , with θ0 and θ1 given.
The secant method needs two initial guesses. For the first, we choose k1 = 0, simply because theangle error is trivial to calculate at that value. For the second guess, the approximation of equation
8.31 works well for θ0 + θ1 < 6 or so, but as this value approaches 2π, it significantly overestimates+θ1 3
k1 and sometimes fails to converge. The guess 6(θ1 − θ0 )(1 − θ02π
) has been found to guarantee
convergence for all θ0 + θ1 < 2π, and is also slightly more accurate, thus requiring fewer iterations
to converge to a given precision.
Complete code for this method is appealingly concise; the Python version (shown in Figure 8.3)
is only a dozen lines of code. Further, only a small handful of iterations is required to approach
floating-point accuracy across the range, and a mere two or three for small angles. Given fast codefor evaluating the spiro function, this approach is efficient indeed.

8.3

Global constraint solver

A common problem in computing interpolating splines is solving a global system of constraints.
For example, in the framework of two-parameter extensional splines, the constraints are curvature
continuity across control points. Other splines may demand a more complex set of constraints, suchas the generalized four-parameter splines of Section 7.4.
Most generally, the goal of the solver is to produce a vector of parameters, such that all thestated constraints are met. Most often, the relationship between the parameters and the corresponding constraints is nearly linear. For example, as demonstrated in Section 8.2, in the small angle
129def fit_euler(th0, th1):
k1_old = 0
e_old = th1 - th0
k0 = th0 + th1
k1 = 6 * (1 - ((.5 / pi) * k0) ** 3) * e_old
for i in range(10):x, y = spiro(k0, k1, 0, 0)
e = (th1 - th0) + 2 * atan2(y, x) - .25 * k1
if abs(e) 1e-9: break
k1_old, e_old, k1 = k1, e, k1 + (k1_old - k1) * e / (e - e_old)
return k0, k1

Figure 8.3. Python code for computing Euler spiral parameters.
approximation the relationship between endpoint tangent angles and curvature is given by linear
Equations 8.29. Thus, one viable approach is to express the relationship between the parametersand constraints as a matrix, and solve it numerically. Hobby’s spline [37] uses this approach, andachieves approximate G2 -continuity. The resulting matrix is tridiagonal, which admits a particularly
fast O(n) and simple solution.
However, this approach yields only an approximate solution, and as the bending angles increase,so do the curvature discontinuities in the resulting spline. A good general approach is to use Newton
iteration on the entire set of vector parameters, using the Jacobian matrix of partial derivatives torefine each iteration. When a solution exists, convergence is quadratic, and experience shows thatthree or four iterations tend to suffice. This general technique resurfaces a number of times in theliterature, including Malcolm [54], Stoer [82], and Edwards [18], the latter being an especially clear
and general explication of using Newton-style iteration to solve nonlinear splines.
Newton iteration is one of the most efficient techniques, but far from the only workable approach.
Since many of the splines are expressed in terms of minimizing an energy functional, a number of
researchers use a gradient descent approach, such as Glass [26], and, more recently, Moreton [64].
However, gradient descent solvers tend to be much slower than Newton approaches. For example,
Moreton’s conjugate gradient solver for MVC curves took 137 iterations, and 75 seconds, to computethe optimal solution to a 7 point spline [64, p. 99]. By contrast, the techniques in this chapter easily
achieve interactive speeds even in a JavaScript implementation that runs in standard Web browsers.

8.3.1

Two-parameter solver

For two-parameter splines, the solution is particularly easy. The vector of parameters is thetangent angle at each control point, and the constraints represent G2 -continuity across control
points. Thus, the vector to be solved consists of the difference in curvature across each control
point.
The vector of tangent angles is initialized to some reasonable values (in our implementation,
the initial angles are those assigned by Séquin, Lee, and Yen’s circle spline [80], although likely a130simpler approach would suffice). At each iteration, the curvatures at the endpoints of each segment
are computed (κli representing the curvature at the left endpoint of segment i, and κri the curvature
at the right). The differences in curvature ∆κi = κri − κli+1 form the error vector. At the same time,
the Jacobian matrix represent the partial derivatives of the error vector with respect to the tangent
angles,
Jij =

∂∆κi
.
∂θj

(8.32)

Then, a correction vector ∆θ is obtained, by solving the Jacobian matrix for the given error
vector. Since the Jacobian matrix is tridiagonal, the solution is very fast.
∆θ = J−1 ∆κ .

(8.33)

Finally, the θi values are updated for the next iteration. In all cases but the most extreme, thenew value is just θi + ∆θi . When bending angles are very large, this iteration can fail, and a second
pass is attempted, using only half ∆θi as the update; convergence is then only linear rather thanquadratic, but it is more robust.
While each segment is specified by two parameters, this solver is written in terms of one parameter per control point, the tangent angle. This approach guarantees G1 -continuity by construction,
as the same tangent is used on both sides of a control point.
One approach to determining the second parameter is to use an additional local solver, determining the two spiral parameters from the two tangent endpoints. For the Euler spiral spline, thislocal solver is described in Section 8.2. Another approach (described in more detail in Section 8.4)
is to precompute a two-dimensional lookup table with the values.
In the software distribution, the clearest and most robust implementation of this solver is inspiro.js. The top-level Newton iteration is contained in refine euler, the computation of thecoefficients of the Jacobian matrix is in get jacobian g2, and the solution of the matrix uses thegeneric band diagonal solver bandec and banbks, adapted from Numerical Recipes [71].

8.3.2

Four-parameter solver

The same general approach is used for four-parameter splines, but the layout of parameters is
a bit different. For one, the solver directly computes all parameters of each segment of the spline,
rather than using a two-stage approach. This solver implements the general constraints as describedin Section 7.4.
The inner loop for computing the error vector is essentially a spiro integral (see Section 8.1). For
the Jacobian matrix, for simplicity we use central differencing to compute the partial derivatives,
even though an analytical technique would also be feasible.
Each constraint gets one element in the error vector and one row in the Jacobian matrix. For
a constraint of the form, say, κ002 = κ003 , the entry in the error vector is the deviation from that
equality, in this case κ002 − κ003 . Each column in the Jacobian matrix represents one parameter of thepolynomial spline, i.e. there are four times as many columns as curve segments. For the system
to be well-determined, there must be an equal number of (linearly independent) constraints as131Set up structure of matrix:
column = parameter of spline
row = constraint
Initialize parameter vector to all zero
Repeat until done:
Compute error vector of deviations from constraints
Compute norm of error vector
Below total error threshold?
If so, done
Compute partial derivative matrix of constraint wrt spline parameter
Invert matrix (band-diagonal)
Compute dot product of error vector and inverse of matrix
Add to parameter vector

Figure 8.4. Pseudocode for computing four-parameter splines.
parameters, i.e. the Jacobian matrix must be square. Thus, the constraints described in 7.4 arecarefully designed to make the count come out even. In particular, for a sequence of G2 -continuous
curve points, additional constraints zero out the higher-order derivatives of the polynomial spline,
so that the resulting curve segments are segments of the Euler spiral.
The algorithm for the four-parameter solver is presented in pseudocode form in Figure 8.4.
Complete working code for the solver is included in the open-source libspiro package. This package
contains an efficient C implementation of this solver, the numerical methods being contained inspiro.c. The outer loop of the Newton iteration is contained in the spiro iter function, with thecomputation of the Jacobian matrix in compute pderivs, and the solution of the matrix in bandec11
and banbks11, again adapted from Numerical Recipes, specialized to the size of the band.
An example of the four-parameter solver at work is shown in Figure 8.5. The first iteration (i.e.
the linear approximation) is the thin red curve, and the second is in green. By the third iteration,
the constraints are satisfied to within visual tolerance.
In practice, both solvers run efficiently and robustly. For the same problem, the two-parameter
solver is a bit more robust than the four-parameter, but of course the latter is much more generaland can solve a much wider class of spline problems.

8.4

2-D Lookup table

The family of two-parameter, extensional splines admits a particularly simple and efficient implementation based on 2-D lookup tables. Since the curve segments are defined by two parameters,
for any given generator curve it is practical to precompute a lookup table indexed by the tangent
angles at the endpoints. The lookup table stores the curvatures at the endpoints (used for enforcing
curvature continuity across control points), as well as the position and tangent angle of a midpoint,used for subdivision while rendering the curve. With this precomputation done only once, drawing
132

Figure 8.5. Example of four-parameter solver.

133

4
3
2
1
0
-1

θ1

8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8

-2
-3
-4
-4

-3

-2

-1

0

1

2

3

4

θ0

Figure 8.6. Contour plot of curvature lookup table for Euler spiral.

an arbitrary spline segment is only slightly more complex than the standard de Casteljau subdivision for rendering Bézier curves; the only difference is the use of a lookup table rather than analgebraic formula to determine the midpoint.
An example curvature lookup table, for the Euler spiral, is shown in Figure 8.4. The horizontal
axis is the tangent angle, in radians, on the left side, the vertical axis is the tangent angle on theright, and the value is the resulting curvature on the left side; the other value follows by symmetry.
Moreover, the table itself is anti-symmetric, which can further reduce storage requirements. Notethat the relationship is nearly linear for small angles. The slope of that linear relationship is closely
related to locality. Because these values are so smoothly behaved, a simple representation as a twodimensional lookup table with 32 × 32 entries, and using bicubic interpolation for lookups, gives amaximum curvature error of 10−4 across the range of the function. This error scales with the fourth
power of the linear size of the table; thus, for example, if a higher precision is desired, a 64 × 64
table could reduce the error to about 6 × 10−6 .

134Chapter 9

Converting to Bézier curves
By far the most popular representation of curves is piecewise polynomial, most often with thecoefficients represented in Bernstein–Bézier form. This popularity is well earned – the representation
is simple, numerically stable, and the primitive is quite supple, capable of representing a wide rangeof curves with relatively few control points. Both cubic and quadratic Béziers are widely used. The
PostScript language [70] established cubic Béziers as a de facto standard for 2D graphics, and the
TrueType font format and Macromedia Flash are based on quadratics, presumably because therasterizer is even simpler to implement than for cubics.
There are several reasons to convert piecewise clothoid curves, or piecewise polynomial spiralcurves, into Bézier representation: interactive display using graphics API’s based on the PostScript
imaging model, creation of font files in industry-standard formats, and export of curves to CAD
and other design tools. All such tools are capable of accepting Bézier curves, but few, if any, canhandle piecewise clothoids.
The conversion to Bézier form involves a tradeoff between the conciseness of the result and theamount of computation needed. For some applications, such as interactive display, the curve will becomputed, displayed, and then discarded. A somewhat wasteful Bézier representation is perfectly
acceptable as long as the overall computation time is small. For other applications, the curve will becomputed, saved into a file such as a font, then embedded thousands of times in document files, eachof which in turn may be rendered thousands of times. In such applications, a significant investment
in computing time to reduce the size of the Bézier representation can pay off handsomely. Thissection will describe both fast and cheap methods for converting polynomial spirals into Béziers,
as well as much slower techniques, suitable for creating optimally concise Bézier representations.
There are several techniques in the literature for creating polynomial approximations of clothoid
curves. Sánchez-Reyes and Chacón [77] use a family of polynomials called s-power series, which isa fairly powerful and general technique. However, the error in their technique converges only asO(n4 ) in the number of subdivisions, while convergence of O(n6 ) is possible for cubic Béziers.

135

(x3, y3)(x2, y2)(x1, y1)(x0, y0)

Figure 9.1. A cubic Bézier.

9.1

Parameters

Assuming fixed endpoints, quadratic Béziers have two free parameters. These are exactly sufficient to determine the tangent angles at the endpoints. Thus, the problem of generating a quadraticBézier given the endpoints is fully determined; the only remaining problem is to determine the endpoints, or, in other words, how to subdivide the original curve to be fit.
Cubic Béziers have four parameters. Two of these determine the tangent angles. The other two
have no direct geometric interpretation, and can be considered free parameters useful for making
the curve fit the target more precisely.
Cubic Béziers are often described in terms of two control arms, or vectors from the endpoints
to the control points. In the typical Bézier of Figure 9.1, these vectors are (x0 , y0 ) to (x1 , y1 ) and(x3 , y3 ) to (x2 , y2 ), respectively. The tangent angles of these control arms are the tangent angles of
the resulting curve at the endpoints, so the lengths of these control arms are the free parameters.

9.2

Error metrics

Fitting of a cubic Bézier to a target curve can be seen as a simple optimization problem. Overthe parameter space, what set of parameters produces the closest fit to the target curve? However,for such a formulation to be meaningful, we must choose an error metric. What does it mean for
one fit to be “closer” than another to the target curve?
In the literature, the maximum distance deviation, or “flatness” is the most commonly used136

Figure 9.2. Two circle approximations with equal distance error.

metric. However, this metric has several shortcomings. To address these shortcomings, this section
presents other error metrics based on angle, rather than distance.
Figure 9.2 shows two approximations to a circle with equal maximum distance error. The
approximation on the left is quantized to an integer grid, while the approximation on the right issmooth but warped out of shape. Even though the actual maximum distance error from the true
circle is the same, the one on the right looks considerably smoother. Thus, it is clear this metric does
not accurately predict perceived smoothness, and optimizing to minimize it may not necessarily
produce the best results.
Thus, other error metrics seem more promising, in particular those based on angle error rather
than absolute distance error. Assume, for the sake of simplicity, that the Bézier approximation andthe target curve have identical total arc length l. Then, we propose the absolute angle metric
Z l|θapprox (s) − θ(s)| ds .

Eabs =(9.1)

0Note that this error metric is extremely closely related to the standard flatness metric. Like
the flatness metric, it is expressed in units of distance, and indeed it represents a fairly tight upperbound on the flatness. Consider the distance between points at arc length d ≤ l from the curvestart, here using complex coordinates for algebraic convenience:
Z deiθapprox (s) − eiθ(s) ds .

(9.2)

|eiθapprox (s) − eiθ(s) | ≤ |θapprox (s) − θ(s)| .

(9.3)

zapprox − z =0

But note the inequality

Note also that the integral 9.2 is bounded from above by Eabs because the absolute angle
error is everywhere positive. Thus, for all points on the original curve, there exists a point on theapproximate curve no more distant than Eabs from this point, and thus from the target curve.
137

Figure 9.3. Euler spiral approximation drawn using .03 and .001 error threshold, respectively.

However, this metric still penalizes “bumps,” or short segments of relatively large angle error
relatively lightly. Better visual results may thus be derived from minimizing the square of the angle
error, integrated over the arc length as before.

Esq =

Z l

2

(θapprox (s) − θ(s)) ds

1/2

.

(9.4)

0

√
From standard properties of norms, lEsq ≥ Eabs holds, which means that Esq also places astrong upper bound on the flatness error metric. This property implies that minimizing the Esq
metric will also yield a curve with a good flatness metric as well, although the converse is not always
true.
Figure 9.3 shows examples of this error metric in use. Both approximations of the Euler spiralsegment use a number of cubic Bézier curves, each with the approximately the same Esq error
metric. Each segment on the left has an error of approximately 0.03, and each segment on the rightapproximately 0.001 (relative to the total size of the figure being about 1). As can be seen, anerror of 0.03 results in a somewhat distorted curve, while 0.001 is essentially indistinguishable fromthe original curve. This figure also provides evidence that converting to Bézier form by solving an
optimization problem in terms of the Esq metric gives good results, even for a relatively crude error
threshold.

9.3

Simple conversion

A simple, fast conversion to cubic Bézier form is as follows. Like the s-power series, and likequadratic Béziers, the convergence is only O(n4 ), although the constant factor is considerably better
than for quadratics, and this cubic solution can easily accommodate an inflection point.
In this solution, set the length of both control arms to one-third the total arc length of the138

Figure 9.4. Three Bézier solutions: fast, equal arc length, and fully optimized. Curvatureplot is the top graph, the Bézier is on bottom.

curve. An example, fitting the section of Euler spiral from parameter t = 0.5 to t = 1.1, is shownas the leftmost instance in Figure 9.4. Both the Bézier (bottom curve) and the plot of curvature
against arc length (top curve) are shown. In the ideal Euler spiral, of course, curvature is a linear
function of arc length.

9.4

Equal arc length, equal arm length solution

A somewhat more accurate solution can be obtained by solving the length of both control arms
so that the total arc length of the Bézier matches that of the target curve. An example is shown as
the top right instance of Figure 9.4.
There is no simple formula for the total arc length of a Bézier, so for this solution we rely on
numerical integration. Since the relation between the control arm length and the total arc lengthis smooth and monotonic, a simple numerical solving technique such as the secant method gives
robust, accurate results in a few iterations.

139

9.5

Fully optimized solution

In the fully general case, the lengths of the control arms need not be equal. To find the single
best-fitting cubic Bézier is, thus, a minimization problem over a two-dimensional parameter space– the lengths of the control arms.
An example of the error metric over this parameter space is plotted in Figure 9.5. This plot is
the root mean squared angle error for all combinations of left and right control arm length, with atarget of the Euler spiral segment from t = 0.5 to t = 1.0. For clarity, only isolines are shown.
The first notable feature of this plot is that there are two local minima of the error metric. Thus,
care must be taken to choose the best global minimum. The actual curves and curvature plots ofthese two local minima are shown in Figure 9.6. In this case, the one in which the left control arm
is longer is clearly the better fit. In general, though, either local minimum may be globally optimal.
Figure 9.7 shows a fit of a different section of the Euler spiral, from t = 0.5 to t = 0.85, and thetwo local minima are nearly indistinguishable, both in overall error metric and in curvature plot.
Also plotted in Figure 9.5 is the curved line representing all values of left and right control arm
length so that the total arc length of the Bézier matches that of the target curve. Notably, error is
locally minimized across this line, suggesting a numerical approach where the outer loop searches
for the optimum ratio between left and right control arm lengths, and an inner loop determines
the total length of the control arms, given their ratio fixed, such that the arc length of the Bézier
matches that of the target curve.
As can be seen in the rightmost example in Figure 9.4, the fully optimized solution finds a Bézierthat matches the original Euler spiral segment with impressive accuracy. Even in the curvature plot,
the deviation from ideal is barely perceptible.

9.6

Global optimum

Generating an optimized Bézier outline of a shape requires not only computation of individual
segments, but also deciding the subdivision points at which the original curve is to be divided.
In general, an optimized Bézier representation of a curve has three elements:• The number of subdivisions is minimal with respect to the error metric threshold.
• The position of the subdivision points minimizes the error metric.
• Each individual segment minimizes the error metric given fixed endpoints.
Of course, the third element is the problem addressed in the previous section. This section
describes techniques to determine the subdivision points, and essentially uses the generation ofindividual Bézier segments (with computation of error metric) as an inner loop.
As before, there is a tradeoff between the amount of time spent optimizing the Bézier representation and the degree to which the representation actually meets these goals. For applications
such as quick interactive display, where the Bézier will be drawn once, it makes sense to computeit as quickly as possible, even if the number of segments is substantially above the minimum. For
140

0.4

1e-02

lines of equal error metricarc length = original curve

0.35

1e-03

right control arm length

0.3

1e-04
0.25

1e-05
1e-06

0.2

0.15

0.1

0.05

0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

left control arm length
Figure 9.5. Approximation errors over parameter space for cubic Bézier fit.

141

0.4

Figure 9.6. Local minima of cubic Bézier fitting, t=[0.5, 1.0]. Curvature plots are the top
two graphs, the Béziers are on the bottom.

142

Figure 9.7. Local minima of cubic Bézier fitting, t=[0.5, 0.85]. Curvature plots are the top
graphs, the Béziers are on the bottom.

143

thresh=0.01 optim=3 tot segs=30

Figure 9.8. Conversion into optimized cubic Béziers.

other applications, such as generating font files to be embedded in documents served over the Web,it makes sense to spend a lot of computing time up front to optimize the representation.
Figure 9.8 shows an example of an optimized conversion into cubic Bézier representation, suitable for encoding into a PostScript Type 1 font format or any of its successors, such as CFF(Compact Font Format), which can be directly embedded into a PDF file, or included in an OpenType container. This font format has the additional requirement that subdivision points appear at
all horizontal and vertical tangents. For conversion into simple vector drawings such as PostScriptor SVG (Scalable Vector Graphics), this requirement could be relaxed.

9.6.1

Hybrid error metric

The purest approach to the error metric for the curve would be to use the same metric as for
an individual segment. For a metric based on L2 -norm, this is computed easily enough by summing
the squares of the individual error metrics of the segments.
However, it is significantly easier to work with a hybrid metric in which an L∞ -norm is used
to combine the error metrics; in other words the error metric for the whole curve is the maximum
metric for each individual segment. If the metric for individual segments is also an L∞ -norm, thisis equivalent to the purest approach. However, the L2 -norm is easier to compute accurately using
numerical integration (it is a considerably smoother function), so the result is a hybrid approach.
The error metric for an individual segment is clearly monotonic with respect to the subdivision
points on either side. Otherwise, if making the segment shorter by moving a subdivision point
inwards were to increase the error metric, the new curve must clearly not be an optimum Bézierfit. Given this property and the objective of minimizing an L∞ -norm, all optimum Bézier repre-144sentations of curves have the property that each individual Bézier segment has an identical error
metric.
In theory, it would be possible to optimize for the pure L2 -norm, but in practice it would becomputationally much more expensive, and the results would be almost imperceptibly different. For
reasoning about this approach, it is perhaps simplest to accept the L2 -norm metric as a reasonably
good approximation to the L∞ -norm.

9.6.2

Simple subdivision

The simplest technique for meeting a particular error metric is to simply fit a single Bézier
segment to the curve, accept it if the error metric threshold is met, and subdivide the curve in half
if not.
Note that if the error metric is an L∞ -norm, then the same error threshold can be applied√ toboth subdivisions. If the error metric is an L2 -norm, then the threshold for each subdivision is 0.5of the original threshold.
One refinement is to subdivide at the midpoint of angle traversed, rather than at the midpointof arc length. For sections not including an inflection point, this is simply the angle halfway between
the tangent angles of the subdivision points, but for sections for which curvature changes sign, amore sophisticated approach is needed, using the integral of the absolute value of curvature as thequantity to subdivide.
Even though simple subdivision doesn’t produce fully optimum results, the number of segments
is fairly good. Using the same ‘a’ from Cecco, using full optimization for the individual Bézier segments, as described in Section 9.5, results in 34 segments. Using the less computationally expensive
equal arc length approach as described in Section 9.4 increases the count to 47. The conversions
are shown in the upper right and upper left quadrants, respectively, of Figure 9.9.

9.6.3

Smarter subdivision

Simple subdivision often results in more segments than absolutely needed. In order to produce
an actually optimum Bézier representation, we propose a smarter subdivision technique. It consistsof four nested loops. In order of outer to inner:• Choose the number of subdivisions needed (a simple approach is to count up from 1, and stop
when the error threshold is met).
• Start with an initial dumb subdivision (equal arc length is fine), then iteratively refine bychoosing the pair of adjacent segments with the greatest discrepancy of error, and finding thesubdivision point joining those segments so that the error is equal.
• A slightly modified secant technique selects the new subdivision point so that the error metricof the two segments is balanced. Simple bisection would also work but its convergence isn’t
quite as fast.
• The inner loop computes the optimized Bézier segments (along with the error metric) on
either side of the control point, using the algorithms of Section 9.5.
145

thresh=0.01 optim=0 tot segs=47

thresh=0.01 optim=1 tot segs=34

thresh=0.01 optim=2 tot segs=30

thresh=0.01 optim=3 tot segs=30

Figure 9.9. Four different levels of optimization: subdivision in half with equal arm length,
subdivision in half with fully optimized Bézier segments, greedy subdivision with fully
optimized segments, and fully optimized.

146

Figure 9.9 shows four different levels of optimization. In the top two examples, each segment is
subdivided exactly in half (by arc length) until the error metric is below the threshold. In the top
left, the Béziers are chosen to have equal arc length as the true curve, but with equal length control
arms, as in Section 9.4. In the top right, each Bézier segment is fully optimized, as in Section 9.5.
In the global optimum (the lower right corner of Figure 9.9), the locations of the subdivision
points are themselves chosen to be optimum. The lower left corner is another intermediate stage of
optimization, to be discussed in the next section.

9.6.4

Alternative optimum algorithm

The algorithm presented in the previous section gives optimum results reasonably quickly.
This section presents an alternative algorithm, which also produces optimum results, but usingslightly different techniques. It is in some ways more systematic, in particular immediately yielding
the number of segments needed, rather than relying on a search to find that number. Also, thistechnique may be faster in cases when only a minimal number of segments is needed (meeting theerror threshold), rather than the additional criterion that the error metric be minimized given that
number of segments.
Again, the central observation is that in an optimized representation, the error metric for eachsegment is the same. Given the error metric target, the algorithm walks from one endpoint of thecurve along it, until the end is reached, one segment at a time. For each segment, the left endpointis already known (again, starting at the endpoint of the curve). Bisection (or some other variantof the shooting method) determines the right endpoint such that the error metric of the segment
is equal to the target. Keep in mind that the error metric of the fit is monotonic in the position of
the right endpoint, so bisection is guaranteed to work.
Simply setting the error target to the threshold yields a minimal number of segments, but thelast segment is likely to have a much lower error metric than the target (and, thus, be shorter thanneeded). An example is shown in the lower left corner of Figure 9.9. Note that the total numberof segments (30) is equal to the fully optimized case in the lower right corner. Yet, some of thesubdivision points seem rather lopsided. For example, the subdivision point on the top curve of thebowl is too far to the right; the error on the curve to its left is exactly the target, and to its right,
much smaller. There are other similar examples.
Even so, this technique can be adapted (at the expense of more computation) to producing afully optimized solution. The key is to determine the error target such that the last segment has
the same error metric as the rest of the segments. If the target is too small, then the algorithm
will use more segments than needed. If too large, then the error metric of the last segment will besignificantly less than all the others. Thus, an outer bisection loop can determine this error target.
An additional small optimization can speed the convergence of this bisection: in trials where
the error target is too large, the error metric of the last segment yields a lower bound on the final
error target for the curve; if it is larger than the lower bound for the bisection, it can replace it.
This is most helpful especially when the guess is close to the actual value.
These results are essentially identical to those of the algorithm of the previous section, andoften slower. This algorithm requires 23.3 seconds (in a slow Python implementation) to produce
the lower left corner of Figure 9.9 and 129 seconds for the fully optimized version. The algorithm of
the previous section gives equivalent results in 48.5 seconds, and for other examples the discrepancy147

is even greater. In settings where the absolute minimum number of Bézier segments is not required,simple subdivision with fully optimized segments produces the upper right corner in 8.6 seconds,
with 34 segments, and simple subdivision with equal arm length, equal arc length produces theupper left corner, 47 segments in 0.9 seconds.
Such times are still undesirably large for interactive use. Clearly, reimplementation in a fasterlanguage such as C would improve the compute time by at least an order of magnitude, and othernumerical improvements may also be possible. But as the last phase in a font designing process thatlasts many hours, even these speeds are acceptable. The next chapter gives discusses applicationsof interpolating splines for the production of fonts in more detail.

148

Chapter 10

Applications in font design
A major motivation for this work is to develop better tools for producing fonts.
Existing font formats in the literature use many different representations of the curves. Usually
the user interface for designing the curves corresponds closely to the representation. Today, by far
the most common representation is Bézier curves, both cubic in the case of PostScript Type 1 fonts[2], and quadratic in the case of TrueType fonts [3]. (The more recent OpenType specification [69]can serve as a container for either of these two underlying formats).
However, a number of alternate formats, popular in their day, use other techniques. Metafont[37] and Ikarus [42], in particular, both used their own flavor of interpolating splines. Perhaps notcoincidentally, both these formats supported variation, or representing a parametrized family of
shapes rather than a single static outline.
Font shapes have a number of requirements beyond simple interpolating splines. The shapes ofoutlines have corners and straight line sections in addition to smooth curves. In general, the outlines
forming a letter shape are closed, although in general not entirely made from smooth points likethe inner and outer contours of a typical ‘o’.
Further, letter-shapes often contain transitions between straight and curved sections. Theseshape patterns are easily expressed using the constraints described in Section 7.4, and the resultspreserve G2 -continuity, yielding visually smooth transitions. Table 10 summarizes these constraints,
as well as the visual symbol used to represent them, used consistently in the following font outlines.
Inconsolata is an original design by the author, drawn entirely using the splines described in thisthesis. Most of the lowercase alphabet is shown in Figures 10.1 and 10.2. Most of the curves weredrawn using G4 -continuity, but using some G2 points when curvature varies more dramatically, forexample in the inner curve of the tail of the ‘g’. Of course, straight-to-curve transitions (which haveG2 -continuity) are used extensively as well, for example in the ‘u’.
The design intent is classical, and draws inspiration from such metal types as the Franklin
Gothic series by Morris Fuller Benton (see [58, p. 142] for showings of this and many other metal
fonts), as opposed to more recent designs such as TheSans Mono, which show their Bézier underlying
format rather clearly.

149

Figure 10.1. Selected glyphs from Inconsolata.

150

Figure 10.2. Selected glyphs from Inconsolata.

151

Type of constraint

SymbolLeft endpoint
Right endpoint
G2 curve point
G4 curve point
One-way (curve to terminal)
One-way (terminal to curve)
Table 10.1. Constraints used to draw font outlines.
Inconsolata is also used for the typewriter style in the formatting of this thesis.
Figures 10.3 and 10.4 show a lowercase alphabet carefully traced from ATF Baskerville Italic(again, see [58, p. 26] for a sample). These outlines generally use G4 -continuous constraints, andoften use one-way constraints to isolate sections of high curvature. Very special care was taken to
trace the outlines faithfully. Figure 10.5 shows one scan (enlarged 25× from the original) with theoutline superposed.
Figures 10.6 and 10.7 are a draft of the capital letters from Cecco, a more recent original designby the author. As opposed to Inconsolata, the curves were drawn primarily with G2 -continuous
curves. These curves are visually just as smooth as the G4 , and fewer points were needed, ingeneral, to describe them. These shapes are also classical, inspired in large part by Pacioli by
Giovanni Mardersteig [55].
Many interpolating splines in the literature are used for generating variations in fonts, rather
than a single static shape. Since fonts are desired in a wide range of weights, one of the morecommon applications is to generate these continuously. One of the most straightforward techniques,pioneered by Ikarus [42], is to start with two masters with similar structure, representing thinand thick extremes, and linearly interpolate the positions of the control points between the two.
Then, an interpolating spline is fit through these control points. Note that there are two senses ofinterpolation here – one for the control points, another for the spline.
The splines described in this thesis work especially well for such interpolation. An example
interpolation is shown in Figure 10.10, derived from the masters in Figure 10.9. Figure 10.11 moreclearly shows the nature of this interpolation, with the linear motion of the control points withrespect to the interpolation parameter represented by the equal spacing of the resulting contour
lines.
One advantage of the splines of this thesis over the Ikarus spline is that fewer control points
are needed to achieve good smoothness, as shown in Figure 10.8. As can be seen, approximately
half as many control points are needed, and the smoothness and quality of the result is favorable.
Second, the one-way constraints are more robust under interpolation than ordinary interpolating

152

Figure 10.3. Selected glyphs from Baskerville Italic.

153

Figure 10.4. Selected glyphs from Baskerville Italic.

154

Figure 10.5. Lowercase ‘w’ from Baskerville Italic, scan with outline.

splines, for such things as straight-to-curve joins, as in the upper right corner curve of this “a”. In
particular, G2 -continuity is always preserved.
By contrast, as shown in Figure 10.12, interpolation between two curves defined by multiple
cubic Bézier segments doesn’t preserve even G1 -continuity, and a “kink” can appear, even when
the source curves have a high degree of continuity (in this case, G∞ across the join, because both
source curves have been subdivided from a single Bézier curve). Even though the figure shows afairly extreme example, it is likely that unless special care is taken, small discontinuities can arise
in ordinary use (perhaps a degree or less, so barely visible). Such cases particularly arise wheninterpolating between different widths, as in this example; it requires special care to place thesubdivisions so that the interpolations are continuous.
In this particular example, the interpolation master curves are subdivided from a single Bézier
segment. If the subdivision point had been chosen at the same paramter for both, then the result
would have been at least G1 -continuous. Alternatively, if the curves had been subdivided at the
point of the same tangent angle, the result also would have been G1 -continuous. It is certainly
possible for designers to carefully choose subdivision points to achieve smooth results, as witnessedby a relatively small number of Multiple Master releases from designers such as Robert Slimbach.
However, this requirement leads to a steeper learning curve and entails ongoing difficulties in use.
Using a spline that is G2 -continuous by construction promises to be better in both respects, andmake the design of fonts with variation more accessible and productive than with previous tools.

155

Figure 10.6. Selected glyphs from Cecco.

156

Figure 10.7. Selected glyphs from Cecco.

157

(b) Changes are marked in black(c) Final IK-points
Figure 10.8. Typical letter-shape using IKARUS, with Spiro comparison.
Figure 6. Comparison of the hand-digitized result with the original

Figure 10.9. Interpolation masters.

Figure 10.10. Interpolation between extremes.

158

Figure 10.11. Interpolation between extremes, superimposed.

Figure 10.12. Interpolation between cubic Béziers fails to preserve G1 -continuity.

159

Chapter 11

Space curves
The focus of this thesis has been on curves in the plane, but much of the literature on splines
is also concerned with the generalization to space curves – a continuous mapping from arc length
to points in space, R → R3 . For example, Moreton’s investigation of the MEC and MVC [64] isequally concerned with space curves as with planar curves.

Figure 11.1. Space curves rendered into ironwork.
Space curves have many applications, including architectural ironwork. The interpolating spline
techniques of this thesis have already shown much promise for this application. Figure 11.1 shows160

ironwork space curves which have been constructed starting from planar interpolating curves. These
pieces are by Terry Ross, as is the Curve Maker software (implemented in the Ruby language asa plugin for the Google SketchUp 3D CAD tool), adapted from the algorithms described in thisthesis [74].
Many of the techniques presented here are also applicable to space curves. There are a numberof viable paths.
First, variationally defined splines are also valid in three-space. The MEC, in particular, generalizes straightforwardly. Its definition in three-space is the same as on the plane; it is the curvethat minimizes bending energy (L2 norm of curvature) that also interpolates the points.
Second, space curve splines may be defined as curve segments drawn from a parametrized space,with additional curvature constraints.
This chapter begins with a brief review of the properties of space curves (Section 11.1). Section 11.2 discusses the generalization of the MEC to space curves, including the property that
this family of curves has three parameters, in contrast to two for the planar MEC, as well as
the continuity properties of the space MEC used as an interpolating spline. Section 11.3 presents
Mehlum’s generalization of the Euler spiral to space splines, a particularly appealing approach.
Section 11.4 discusses a generalization of log-aesthetic curves to space. Finally, Section 11.5 considers the possibility of applying similar techniques to the problem of interpolating surfaces, a much
more challenging problem.

11.1

Curves in space

In Cartesian coordinates, the representation of a space curve simply adds a z coordinate to
the x and y of plane curves. The intrinsic representation also adds an additional dimension. Whilea plane curve is characterized by its curvature, a space curve is characterized by curvature andtorsion. The plane curve with constant curvature is a circle, and the corresponding space curvewith constant curvature and torsion is the helix.
An excellent introduction to the differential geometry of space curves is Chapter 2 of Kreyszig’s
Differential Geometry [48]. This section simply restates some of the elementary properties of spacecurves, for convenience.
The moving trihedron (also known as the Frenet frame) consists of three orthogonal vectors.
For a curve defined as a vector x as a function of arc length s, the tangent vector is t = dx/ds, theprincipal normal is defined asp=

1 d2 xdt/ds=.
|dt/ds|κ ds2(11.1)The tangent and principal normal are orthogonal, and both lie in the plane of the osculating
circle. The third component of the moving trihedron is the binormal, and is orthogonal to theosculating circle.
b=t×p.

161(11.2)

Torsion is the rate of change of the direction of the binormal vector.
τ = −p ·db.
ds(11.3)

An important result is that a curve is planar if and only if torsion is everywhere zero.
The concept of geometric continuity is similar in space as in the plane. In particular, a spacecurve has G2 -continuity if its osculating circle is continuous along its arc length. Torsion need not
be continuous in a G2 -continuous space curve.
The bending energy of a wire is, as in the planar case, proportional to the integral of the squareof curvature. Note that the bending energy does not directly depend on torsion. A physical wire
may also have twist, but this is not considered here; for spline applications, it may be assumed thatthe wire can untwist to find its minimum energy configuration.

11.2

Minimum energy curve in space

Likely the first actual solution to the shape of the MEC in space was by Kirchoff, who foundthat the kinetic analogy for the planar elastica (see Section 5.11) also generalizes to three-space
fairly straightforwardly. While the planar elastica corresponds to a pendulum rotating around a
fixed point, constrained to a plane, the minimum energy curve in space corresponds to a spinning
top, also gyrating around a fixed point. Goss, in his Ph.D. thesis, traces out the history of thisdevelopment [29, p. 37]. The result is also presented in Greenhill’s 1892 book on the applications ofelliptic curves [33]. Max Born’s Ph.D. thesis also contains a derivation of the differential equationsof the curve from variational principles [11], but his choice of spherical polar coordinates makes theequations fairly awkward.
A more modern derivation of the shape of the MEC in space is contained in Mehlum’s workon nonlinear splines, published in 1974 [59] and later refined [61, 60]. Most importantly, Mehlum
discovered a simple relationship between curvature and torsion for all solutions of the elastica [59,Eq. 2.21].
κ2 τ = C(11.4)

Mehlum also derived the following equation for the curvature as a function of arc length. Sincea space curve is completely characterized by its curvature and torsion, given the above relation, theshape is determined by a differential equation in one variable.
 2 0 2(κ ) + κ2 (κ2 − 2D)2 − 4ψ 2 + 4C 2 = 0,

(11.5)

where C, D, and ψ are arbitrary constants [59, Eq. 2.24]. Compare with Equation 3.10, the equation
for the elastica in the plane.

162

11.2.1

MEC as spline in 3-space

When used as a spline, the MEC has the following properties:• Each segment is an energy-minimizing space curve, constrained by the tangent vectors relative
to the chord at both endpoints.
• G2 -continuity is preserved across control points, as in the planar case. However, torsion may
be discontinuous.
• Each curve segment is controlled by three parameters. The additional parameter corresponds
to the deviation of the two tangent vectors from being planar.
Each tangent vector has two components (equivalent, for example, to the spherical polar coordinates θ and φ). This would seem to suggest that there are four parameters for each curve segment.
However, counting in this way would include one parameter for rotations around the chord, which,as a rigid body transformation, is not actually a true parameter of the underlying curve family.
Thus, the actual number of parameters is three. The relative rotation around the chord is related
to the overall torsion on the curve; when this rotation is zero, the torsion vanishes and the curve
becomes the planar MEC.
Of course, the MEC in three-space has the same drawbacks as the planar MEC, including thelack of roundness and existence. Thus, a generalization of the Euler spiral is also worth seeking.
As of now, there are two viable proposals for such a generalization: Mehlum’s spiral, discussed in
Section 11.3, and a generalization of log-aesthetic curves to space, discussed in Section 11.4.
Except in the planar case, a space MEC never has an actual inflection point, i.e. a point inwhich curvature is zero. In the nearly-planar case, the curvature drops very low, and torsion peaksin this region, so that the plane of the osculating circle swings through nearly 180◦ .
As for higher degrees of continuity, Moreton [64] described a reasonable (but computationallyintensive) techniques for computing a good approximation to the MVC, providing G4 continuity.
Likely, either an exact solution based on variational techniques, or a good approximation, would
yield similar results with much less computational cost.

11.3

Mehlum’s spiral

One beautiful proposal for generalizing the Euler spiral into space is by Even Mehlum [59]. It
can be understood both as an approximation to the space MEC (entirely analogous to the way inwhich the Euler spiral is an approximation to the planar elastica), and as a geometric constructionin which the Euler spiral is rolled up onto the surface of a sphere.
To derive the curve from the space MEC, assume that κ2  2D. Then, Equation 11.5 simplifiesto 2 0 2(κ ) − 4α2 κ2 + 4C 2 = 0 ,(11.6)

where α is a constant (it is equal to ψ 2 −D2 , combining these two constants into one). This equation,
however, has a simple analytic solution,163κ2 = α2 s2 +

C2
.
α2

(11.7)

6!
asymptote for !
"
5

!, "

4

3

2

1

0
-4

-2

0
s

2

4

Figure 11.2. Curvature as a function of arc length in the Mehlum spiral.
When C = 0, this equation simplifies to κ = αs, the equation of the Euler spiral. In general, it
is the equation of a hyperbola, as shown in Figure 11.2.
A space curve requires both curvature and torsion. For the latter, simply retain Equation 11.4
from the planar MEC. Thus, these equations fully define a family of curves suitable as the basis
for an interpolating spline. This family was used as the basis of the kurgla 2 interpolating spline
algorithm in the Autokon system [59].
Further analysis by Mehlum and Wimp [61] revealed a surprising geometric property of thiscurve – it lies entirely on the surface of a sphere. That analysis is based on a well-known intrinsic
equation relating curvature and torsion, as a necessary and sufficient condition for a space curve tolie on a sphere [61, Eq. 2.8].
τ  κ0  0
−
=0
κ
τ κ2

(11.8)For values of curvature following the function described in Equation 11.7, this relation (after a
bit of algebraic manipulation) is shown to be equivalent to Equation 11.4.
Even more recently, Mehlum discovered a purely geometric construction of this space curve(credited jointly with Professor Ronald Resch), based on rolling up an Euler spiral onto the surface
of a sphere [60]. He writes,164

Draw a Cornu spiral with wet ink in the xz plane. The curvature of the spiral is specified
by κ = αs. Place a sphere with radius R = Cα at the xz plane with the tangent pointat the origin. Roll the sphere along the Cornu spiral such that the axis of rotation
always goes through the center of the sphere and that point on the evolute of the spiralthat corresponds to the touching point of the sphere at the spiral. The sphere is always
tangent to the xz plane during the rotation.
Mehlum’s 1994 paper [60] also contains, in an impressive feat of analysis, closed-form formulas for the Cartesian coordinates of the space curve as a function of arc length, using confluent
hypergeometric functions.
This curve seems a very good choice for an interpolating space curve, both because of its
excellent practical properties (generalized from those of the Euler spiral) and for its mathematicalbeauty. I know of no other implementations than the original Autokon system.

11.4

Log-aesthetic space curves

Another fruitful avenue for exploration is the generalization of log-aesthetic curves to space [91].
Recall that the planar log-aesthetic curve is defined as a straight-line plot of curvature vs (scaled)derivative of curvature when both are plotted on a logarithmic scale. The space log-aesthetic curve
applies the analogous constraint to the torsion as well. Just as the circle is a special case of theplanar log-aesthetic curve, the helix is a special case of the corresponding space curve. There are anumber of additional parameters of this curve family, so determining these parameters remains aninteresting open problem.
The parameter space admits many different types of curves. In particular, some have both
curvature and torsion increasing as a function of arc length. These curves generally resemble acorkscrew. One of Mehlum’s main results is that, in a space curve resulting from the physicalbending of a wire, the curvature and torsion vary oppositely (see Equation 11.4). In regions of highcurvature, torsion is low, and vice versa. Perhaps this heuristic could be applied to curves from thelog-aesthetic family to reduce the dimensionality of the parameter space while preserving desirable
solutions.

11.5

Surfaces

While well beyond the scope of this thesis, some of the concepts described here also generalize
to surfaces. Moreton’s work [64] is one of the earliest presentations of an interpolating surface
defined in terms of minimizing an energy functional. These techniques are known to be extremelycomputationally intensive. Furthermore, the choice of energy functional to minimize is not obvious(just as in the planar case). Among the proposals are MES (Minimum Energy Surface), analogous
to the planar MEC, and MVS (Minimum Variation Surface), analogous to the planar MVC.
One proposal for such an energy functional to minimize is the MVScross functional of Joshi
and Séquin [41]. For the input configurations tested, the surface minimizing this functional is morebalanced and aesthetically pleasing than the pure MES and MVS functionals.
165

A very intriguing direction to explore is whether surface segments may be drawn from aparametrized family in the same way that curve segments are drawn in the planar and space-curve
cases. Such an approach holds the promise of much faster computation than iteratively evolving
a surface to reduce its energy functional. Candidates for such surfaces include a family of patches
that minimize an energy functional (analogous to the MEC), generalizations of the log-aesthetic
curve family to surfaces, and quasi-static liquids [83], which are defined as energy minimizing in
the presence of certain additional constraints.

166

Chapter 12Conclusion
Interpolating splines are a useful and important tool for designing curves. However, choosing one
from the multiplicity in the literature is not easy. While fairness is obviously a crucially important
property, there has been no universal consensus on what makes an ideal interpolating spline, or even
how fairness should be measured. In this thesis, I have argued that an “ideal” spline must share
two properties with the Minimum Energy Curve, the mathematical idealization of the mechanicalspline based on a thin strip of elastic material. First, it must be extensional, meaning that adding
an additional control point on the curve does not change the shape of the curve. Second, the spline
should be based on a two-parameter family of primitive curves, joined together at control points
with G2 -continuity. Higher degrees of continuity (also corresponding to more parameters for theprimitive curve) are possible, but come at the cost of poorer locality.
Based on evidence from existing splines in the literature, two parameter splines have better
locality than splines defined from a larger parameter space. Two-parameter splines can easily achieve
G2 -continuity, and no higher degree of continuity is required for a planar curve to appear smoothto the human visual system.
In fact, there is in general a tradeoff between fairness and locality. A degenerate spline consistingof straight lines has perfect locality but terrible fairness. More usefully, the interpolating spline basedon the log-aesthetic curve covers a range of this tradeoff, based on a scalar parameter (the exponentof the relation between curvature and arc length).
A significant result of this thesis is that all two-parameter, extensible splines can be definedin terms of a generator curve, presented in Section 4.8. As a result, for this important subclass of
G2 -continuous splines, the design space is reduced to choosing the best curve, both in terms of theaesthetic beauty of the curve itself and also the properties of the corresponding spline.
The Euler spiral is an excellent choice for this generator curve. In fact, the Euler spiral is betterin nearly every respect than the MEC, the only exception being minimization of the MEC energyfunctional itself. However, as demonstrated by human experiment (see Section 4.8.1), MEC energy
is not a good predictor of perceived fairness. The Euler spiral spline has identical locality properties,
and, unlike the MEC, is round (producing a circular arc when the control points are co-circular),and a solution nearly always exists.

167

-1 0

1
3
2
1
0 θ1
-1
-2

-3

-2

-1

0

1

2

θ0 (radians)

3

-3

13
12
11
10
9
8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8
-9
-10
-11
-12
-13

Figure 12.1. Contour plot of curvature lookup table for α = 2 log-aesthetic curve.

One promising candidate for the generating curve is the family of log-aesthetic curves, parameterized by the exponent α. At α = 1, the Euler spiral is one instance of this family. Higher values of
α exhibit better locality, but overly large values have poor perceived fairness. However, as measured
by empirical study, values of α between 1.5 and 2 have both better locality and better fairness thanthe Euler spiral. Based on these properties, it would seem like the log-aesthetic spline with a higher
exponent would be an even better interpolating spline than that based on the Euler spiral. Indeed,
for solving a spline based on a static sequence of control points, the log-aesthetic spline tends to
produce good results.
However, in interactive use, a weakness of the log-aesthetic spline emerges. For many configurations of input points, a small perturbation to the input causes a fairly significant change tothe curve, even when the angles are small. One way to analyze this phenomenon is to look at therelation between tangent angles and curvature. While it is very nearly linear for the Euler spiral(see Figure 8.4), the log-aesthetic plot has significant nonlinearities, even for very small angles, as
shown in Figure 12.1.
An excellent way to achieve better fairness (albeit at the expense of locality) is to require higher
degrees of geometric continuity. In the variational framework, the Minimum Variation Curve (MVC)
is a four-parameter spline with G4 -continuity. A strong motivation was to fix the roundness problem

168

of the MEC. While previous researchers have used numerical techniques to minimize the functional(especially Moreton [64]), one contribution of this thesis is to develop the exact general equation
for the MVC. Happily, the MVC has a very simple formula relating curvature with Cartesian
coordinates: κ00 = λy, which also has an appealing parallelism with the corresponding equation for
the MEC, κ = λy.
The MVC has a simple approximation, writing curvature as a cubic polynomial in terms ofarc length. It is entirely analogous to the Euler spiral spline, which can be seen as a simpler, morerobust approximation to the MEC. This curve admits very efficient numerical techniques (presented
in Section 8.1), for fast and robust computation. Interestingly, the cubic polynomial spiral spline
was first described by Ohlin in 1987 [67], but with little impact at the time, probably because Ohlin
did not adequately characterize its properties.
The G4 -continuous polynomial spiral spline (or “spiro” for short) is quite practical, and has
been used to draw many font outlines. Because of its poorer locality properties relative to the Euler
spiral spline, it requires more points to adequately describe letter-shapes without exhibiting theripple effect, but only about a factor of two, based on experience drawing multiple fonts (describedin detail in Chapter 10). Because G2 -continuity is sufficient for the perception of fair curves, I havetaken to using primarily the Euler spiral spline for my later font work.
For pure interpolating splines, the tradeoff between fairness and locality seems somewhat uncompromising. In practice, the biggest problem is that regions of high curvature causes ripples,
analogous to “ringing” in digital filters. Fortunately, adding explicit one-way constraints to thecontrol points isolates these regions of high curvature and prevents ripples entirely. The implementation of these constraints uses the same four-parameter curve primitive as used to generate splines
with G4 -continuity. Instead of enforcing continuity of higher derivatives of curvature, these “oneway” (asymmetrical) constraints enforce zero derivatives of curvature on one side only of the control
point, while still preserving G2 -continuity everywhere. In practice, these constraints are especially
useful for making smooth transitions from straight to curved sections, a common configuration inletter-shapes.
A major motivation of this thesis was to explore the idea of the best possible interpolating spline,
regardless of computational complexity, based on the idea that Moore’s Law would ultimately makeeven computationally extremely expensive splines practical. However, work on efficient implementation proved very fruitful, and the numerical techniques developed in this thesis are competitive
with cruder, polynomial-based splines.
Especially for font design, it is important to represent the resulting shapes in a format thatcan be widely used. Because the industry-standard font formats are all based on Bézier curves, I
developed algorithms for converting curves into this form, producing curves that are visually nearly
identical to the original spline curves (based on a carefully designed error metric), while using an
optimally small number of segments, to keep file size for the resulting fonts small.
A signifcant focus of this work was to quantify the properties of the various interpolating splines
explored in this thesis, to facilitate choosing the best one for any given application. In using these
splines for drawing a number of fonts, I also found they brought an aesthetic quality to the work.
Designing on the screen (rather than starting on paper), I found that I could readily achieve thesmooth, classical shapes I desired. Knowing that the curves have a strong mathematical foundation
added to this aesthetic sense, as did the fact that the numerical techniques were based on rapidly
computing nearly exact analytical solutions, rather than merely crude approximations. It is my169hope that in publishing this work and releasing the code under an open-source license, a wider
audience of designers will experience the pleasure of working with these superior tools.

170

Bibliography[1] Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Functions with Formulas,
Graphs, and Mathematical Tables. Dover, New York, ninth Dover printing, tenth GPO printing
edition, 1964.
[2] Adobe Systems Inc. Adobe Type 1 Font Format, 1990.
[3] Apple Computer Inc. TrueType Reference Manual, October 1996.
http://developer.apple.com/textfonts/ttrefman/index.html.
[4] Raymond Clare Archibald. Euler integrals and Euler’s spiral. American Mathematical Monthly,25(6):276–282, June 1917.
[5] Fred Attneave. Some informational aspects of visual perception.
61(3):183–193, 1954.

Psychological Review,[6] Daniel Bernoulli. The 26th letter to Euler. In Correspondence Mathématique et Physique,volume 2. P. H. Fuss, October 1742.
[7] James Bernoulli. Quadratura curvae, e cujus evolutione describitur inflexae laminae curvatura.
In Die Werke von Jakob Bernoulli, pages 223–227. Birkhäuser, 1692. Med. CLXX; Ref. UB:
L Ia 3, p 211–212.
[8] James Bernoulli. Jacobi Bernoulli, Basilieensis, Opera, volume 1. Cramer & Philibert, Geneva,
1744.
[9] James Bernoulli. Jacobi Bernoulli, Basilieensis, Opera, volume 2. Cramer & Philibert, Geneva,
1744.
[10] Garrett Birkhoff and Carl R. de Boor. Piecewise polynomial interpolation and approximation.
Proc. General Motors Symp. of 1964, pages 164–190, 1965.
[11] Max Born. Untersuchungen über die Stabilität der elastischen Linie in Ebene und Raum, underverschiedenen Grenzbedingungen. PhD thesis, University of Göttingen, 1906.
[12] Max Born. My Life and Views. Scribner, 1968.
[13] Alexander Brook, Alfred M. Bruckstein, and Ron Kimmel. On similarity-invariant fairness
measures. In Scale Space and PDE Methods in Computer Vision, number 3459 in LNCS,
pages 456–467. Springer, 2005.
[14] Ian D. Coope. Curve interpolation with nonlinear spiral splines. IMA Journal of Numerical
Analysis, 13(3):327–341, 1993.
171

[15] Marie Alfred Cornu. Méthode nouvelle pour la discussion des problémes de diffraction dans lecas d’une onde cylindrique. Journal de Physique théoretique et appliquée, pages 5–15, 1874.
[16] Charles Lee Crandall. The Transition Curve. Wiley, 1893.
[17] Lawrence D’Antonio. The fabric of the universe is most perfect: Euler’s research on elastic
curves. In Euler at 300: an appreciation, pages 239–260. Mathematical Association of America,
2007.
[18] John A. Edwards. Exact equations of the nonlinear spline. Trans. Mathematical Software,18(2):174–192, June 1992.
[19] L. Elsgolts. Differential Equations and the Calculus of Variations. Mir Publishers, Moscow,
1970.
[20] Leonhard Euler. Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes,
sive solutio problematis isoperimetrici lattissimo sensu accepti, chapter Additamentum 1. eulerarchive.org E065, 1744.
[21] Leonhard Euler. On the values of integrals extended from the variable term x = 0 up tox = ∞, 2007. Translation into English by Jordan Bell.
[22] G. Forsythe. Notes taken at the GM Research Laboratories Symposium on approximation offunctions. Special collection SC98 2–6, Stanford University, 1964.
[23] A. H. Fowler and C. W. Wilson. Cubic spline, a curve fitting routine. Technical Report Y-1400,
Oak Ridge Natl. Laboratory, September 1962.
[24] Craig G. Fraser. Mathematical technique and physical conception in Euler’s investigation of
the elastica. Centaurus, 34(3):211–246, 1991.
[25] Augustin Jean Fresnel. Memoir on the diffraction of light. In Henry Crew, editor, The Wave
Theory of Light, pages 79–144. American Book Company, 1900.
[26] J. M. Glass. Smooth curve interpolation: A generalized spline-fit procedure. BIT, 6(4):277–293,
1966.
[27] James Glover. Transition curves for railways. In Minutes of Proceedings of the Institution of
Civil Engineers, volume 140, pages 161–179, 1900.
[28] Lawrence E. Goodman and William H. Warner. Statics. Wadsworth, 1963.
[29] Victor Geoffrey Alan Goss. Snap buckling, writhing and loop formation in twisted rods. PhD
thesis, University of London, London, 2003.
[30] Victor Geoffrey Alan Goss. The history of the planar elastica: Insights into mechanics andscientific method. Science & Education, 2008.
[31] I. Grattan-Guinness. Convolutions in French Mathematics, 1800–1840. Birkhäuser, 1990.
[32] Alfred Gray. Modern Differential Geometry of Curves and Surfaces with Mathematica, chapter
3.5, Clothoids, pages 64–66. CRC Press, Boca Raton, FL, 2nd edition, 1997.
[33] Alfred George Greenhill. The Applications of Elliptic Functions. Macmillan, London, 1892.
172

[34] Stefan Gumhold. Designing optimal curves in 2D. In Proc. CEIG, pages 61–76, Sevilla/Spain,
July 2004.
[35] Peter M. Harman and Alan E. Shapiro, editors. The critical role of curvature in Newton’sdeveloping dynamics. Cambridge University Press, 2002.
[36] Arthur Lovat Higgins. The Transition Spiral and Its Introduction to Railway Curves. Van
Nostrand, 1922.
[37] John D. Hobby. Smooth, easy to compute interpolating splines. Technical Report STAN-CS85-1047, Stanford University, 1985.
[38] B. K. P. Horn. The curve of least energy. Technical Report A.I. Memo. No. 612, MIT AI Lab,
January 1981.
[39] Josef Hoschek, Dieter Lasser, and Larry L. Schumaker. Fundamentals of Computer-aided
Geometric Design. AK Peters, 1993.
[40] Robert Alexander Houstoun. A Treatise on Light. Longmans, Green and Co., 1915.
[41] Pushkar Joshi and Carlo Séquin. Energy minimizers for curvature-based surface functionals.
CAD And Applications, 4(5):607–617, 2007.
[42] Peter Karow. Digital Formats for Typefaces. URW Verlag, April 1987.
[43] Peter Karow. Digital punch cutting. Electronic Publishing, 4(3):151–170, September 1991.
[44] Benjamin B. Kimia, Ilana Frankel, and Ana-Maria Popescu. Euler spiral for shape completion.
Int. J. Comput. Vision, 54(1-3):157–180, 2003.
[45] Louis T. Klauder. A better way to design railroad transition spirals. ASCE Journal of Transportation Engineering, May 2001. Date is date submitted.
[46] Donald E. Knuth. Computer Modern Typefaces, volume E of Computers and Typesetting.
Addison-Wesley, Reading, MA, USA, 1986.
[47] Donald Ervin Knuth. Digital typography, volume 78 of CSLI lecture notes. CSLI Publications,
Stanford, Calif., 1999.
[48] Erwin Kreyszig. Differental Geometry. Dover Publications, 1991.
[49] Boris I. Kvasov. Methods of shape-preserving spline approximation. World Scientific, 2000.
[50] Pierre Simon Laplace. Œuvres complétes de Laplace, volume 4. Gauthier-Villars, 1880.
[51] Erastus H. Lee and George E. Forsythe. Variational study of nonlinear spline curves. SIAM
Rev., 15(1):120–133, 1973.
[52] A. E. H. Love. A Treatise on the Mathematical Theory of Elasticity. Dover Publications,fourth edition, 1944.
[53] D. H. MacLaren. Formulas for fitting a spline curve through a set of points. Technical Report 2,
Boeing Appl. Math. Report, 1958.

173

[54] Michael A. Malcolm. On the computation of nonlinear spline functions. SIAM J. Numer.
Anal., 14(2):254–282, April 1977.
[55] Giovanni Mardersteig. The Officina Bodoni: An Account of the Work of a Hand Press, 19231977. Edizioni Valdonega, 1980.
[56] Dan Margalit. History of curvature.
http://www3.villanova.edu/maple/misc/history of curvature/k.htm, 2003.
[57] James Clerk Maxwell. Capillary action. In Encyclopædia Britannica, pages 256–275. Henry
G. Allen, 9th edition, 1890.
[58] Mac McGrew. American Metal Typefaces of the Twentieth Century. Oak Knoll Books, 2nd,
revised edition, 1993.
[59] Even Mehlum. Nonlinear splines. Computer Aided Geometric Design, pages 173–205, 1974.
[60] Even Mehlum. Appell and the apple (nonlinear splines in space). In M. Dæhlen, T. Lyche,and L. L. Schumaker, editors, Mathematical Methods for Curves and Surfaces: Ulvik, Norway.
Vanderbilt University Press, 1994.
[61] Even Mehlum and Jet Wimp. Spherical curves and quadratic relationships for special functions.
J. Austral. Math Soc. Ser. B, 27(1):111–124, 1985.
[62] William Hallowes Miller. The Elements of Hydrostatics and Hydrodynamics. J. & J. J.
Deighton, 1831.
[63] Kenjiro T Miura, J Sone, A Yamashita, and Toru Kaneko. Derivation of a general formula ofaesthetic curves. In In: Proceedings of the Eighth International Conference on Humans andComputers, pages 166–171, 2005.
[64] Henry Moreton. Minimum Curvature Variation Curves, Networks and Surfaces for Fair FreeForm Shape Design. PhD thesis, University of California, Berkeley, Berkeley, California, 1992.
[65] Henry P. Moreton and Carlo H. Séquin. Scale-invariant minimum-cost curves: Fair and robust
design implements. Computer Graphics Forum, 12:473–484, 1993.
[66] Stephen L. Moshier. Methods and Programs for Mathematical Functions. Prentice Hall, 1989.
[67] S. C. Ohlin. Splines for engineers. In Guy Marechal, editor, Eurographics ’87, pages 555–565,
1987.
[68] W. A. Oldfather, C. A. Ellis, and Donald M. Brown. Leonhard Euler’s elastic curves. Isis,
20(1):72–160, November 1933.
[69] OpenType specification, version 1.6 draft.
http://www.microsoft.com/typography/otspec/, April 2009.
[70] Adobe Press. PostScript Language Reference Manual. Addison-Wesley Longman PublishingCo., Inc., Boston, MA, USA, 1985.
[71] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical
Recipes in C. Cambridge University Press, 2nd edition, 1992.

174

[72] Thomas Preston and Charles Jasper Joly. The Theory of Light. Macmillan, 1901.
[73] William John Macquorn Rankine. A Manual of Civil Engineering. Charles Griffin, 17th edition,
1889.
[74] Terry Ross. More about curve maker.
http://www.drawmetal.com/moreaboutcurvemaker.
[75] John Roulier and Thomas Rando. Measures of fairness for curves and surfaces. In Nickolas S.
Sapidis, editor, Designing Fair Curves and Surfaces: Shape Quality in Geometric Modeling andCumputer-Aided Design, chapter 5. SIAM, 1994.
[76] Louis Saalschütz. Der belastete Stab unter Einwirkung einer seitlichen Kraft. B. G. Teubner,
Leipzig, 1880.
[77] Javier Sánchez-Reyes and Jesús Miguel Chacón. Polynomial approximation to clothoids via
s-power series. Computer-Aided Design, 35(14):1305–1313, 2003.
[78] I. J. Schoenberg. Selected Papers, volume I. Birkhäuser, Boston, 1988. edited by Carl de Boor.
[79] Peter Selinger. Potrace: a polygon-based tracing algorithm.
http://potrace.sourceforge.net/potrace.pdf, September 2003.
[80] Carlo H. Séquin, Kiha Lee, and Jane Yen. Fair, G2- and C2-continuous circle splines for theinterpolation of sparse data points. Computer-Aided Design, 37(2):201–211, 2005.
[81] R Sridharan. Physics to mathematics: from lintearia to lemniscate – I. Resonance, pages
21–29, April 2004.
[82] Josef Stoer. Curve fitting with clothoidal splines. J. Res. Nat. Bur. Standards, 87(4):317–346,
1982.
[83] Kurt W. Swanson, Kenneth A. Brakke, and David E. Breen. Physics-based surface modelingusing quasi-static liquids. CAD And Applications, 6(6):759–779, 2007.
[84] Arthur N. Talbot. The Railway Transition Spiral. U. Illinois, 1899. Reprinted from The
Technograph No. 13.
[85] Isaac Todhunter. A History of the Theory of Elasticity and of the Strength of Materials,
volume 1. Cambridge University Press, 1886.
[86] C. Truesdell. The Rational Mechanics of Flexible or Elastic Bodies: 1638–1788. Birkhäuser,
1960. Introduction to Leonhardi Euleri Opera Omnia Vol X et XI, Series 2.
[87] C. Truesdell. The influence of elasticity on analysis: The classic heritage. Bull. AMS, 9(3):293–
310, November 1983.
[88] C. Truesdell. Der Briefwechsel von Jacob Bernoulli, chapter Mechanics, especially Elasticity,in the Correspondence of Jacob Bernoulli with Leibniz. Birkhäuser, 1987.
[89] William Whewell. Analytical Statics: A supplement to the fourth edition of an elementarytreatise on Mechanics. J. & J. J. Deighton, Cambridge, England, 1833.
[90] C. H. Woodford. Smooth curve interpolation. BIT, 9:69–77, 1969.
[91] Norimasa Yoshida, Ryo Fukuda, and Takafumi Saito. Log-aesthetic space curve segments. In
SIAM/ACM Joint Conference on Geometric and Physical Modeling, 2009.
175