\chapter{Two-parameter splines}
\label{two-chapter}

The broad goal of this work is to find the \emph{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:

\begin{quote}
\label{two-param-def}
\textbf{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 the endpoints of the segment.
\end{quote}

Examples of two-parameter splines include the Minimum Energy Curve,
the Euler spiral spline, Mehlum's {\sc kurgla 1}, Hobby's spline,
Ikarus, and S\'equin's circle spline. Note that the count of
parameters excludes the scaling, rotation, and translation required to
get the curve segment to coincide with the two endpoints -- those
additional parameters are uniquely determined by the constraints 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 $G^2$-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 have satisfied the old constraints as well. In addition,
both splines are based on cutting segments from a fixed
\emph{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.

\section{Extensional two-parameter splines have a generating curve}
\label{extensional-equivalence}

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 from that
generator curve, subject to scaling, rotation, translation, and mirror
image transformations to make the endpoints of the cut segments align
with the endpoints in the spline. Thus, any spline in 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 the relevant
equations, one realizes that the schemes of \cite{MacLaren58} and
\cite{Fowler62} do not approximate mechanical splines more closely
than other curves -- nor does it seem particularly desirable to have
them do so. For example, they approximate equally well to Hermite
interpolation by segments of \emph{Euler's spirals,} joined together
with continuous curvature.'' \cite[p. 171]{Birkhoff65}. The Euler
spiral is but one of a family of desirable curves for constructing
extensional interpolating splines.

\newtheorem{theorem}{Theorem}
\begin{theorem}
In an extensional, $G^2$-continuous, two-parameter spline, there
exists a curve such that for any spline segment (parametrized by
angles $\theta_0$ and $\theta_1$) there exist two points on the curve
($s_0$ and $s_1$) 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.
\end{theorem}

\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=0.9]{figs/two_continue.pdf}
\caption{\label{two-continue}Incremental construction of the generator
  curve of a two-parameter spline.}
\end{center}
\end{figure*}

A typical section of curve is shown in Figure \ref{two-continue}. For
some given endpoint angles, the generated curve segment is the section from
$s_0$ to $s_2$. Because of the extensionality, for any $s_1$ along
this curve, the subdivided curve from $s_0$ to $s_1$ (which is the
result of applying $\theta_0$ and $\theta_1$ endpoint constraints) must
be consistent with the original segment. Modulo scaling, the curvature
$\kappa$ and its derivatives $\kappa'$ and $\kappa''$ must of course
also match between these two segments.

The proof of this theorem revolves around the quantity
$\kappa'/\kappa^2$, which is also the limit of
$6(\theta_1-\theta_0)/(\theta_0+\theta_1)^2$ as $\theta_0$ and
$\theta_1$ both approach zero, i.e. as the length of the curve segment
becomes infinitesimal (this limit follows from Equation
\ref{fit-euler-approx}, which is derived in detail in Section
\ref{solve-euler}). This quantity 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 \emph{second} derivative of curvature, suitably normalized for
scale invariance, $\kappa''/\kappa^3$, is uniquely determined given a
particular value of $\kappa'/\kappa^2$. If there are two curve
segments with different values for this second derivative, then it
is possible to cut subsegments from each with identical $\theta_0$ and
$\theta_1$ values. However, the value of $\kappa''/\kappa^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 $\theta_0,\theta_1$ pair -- a
contradiction!

Thus, there is a process for constructing the generator curve. Start
with arbitrary $\kappa$ and $\kappa'$ values sampled from some point
on the curve, then continually derive $\kappa''$ based on the lookup
process suggested by the above relation: find a curve segment that has
the correct value for $\kappa'/\kappa^2$, then sample
$\kappa''/\kappa^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
$\theta_0,\theta_1$ pair, there exists a segment cut from this curve
such that it is equal to the spline segment. Determine
$\kappa'/\kappa^2$ for one side ($s_0$) 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
$\kappa''/\kappa^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 $\theta_0$ and $\theta_1$. We conjecture
that this assumption follows from the assumption of $G^2$-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 \ref{sec-mvc}) is extensional but
does not fit into a two-parameter space, and also does not have a
single generator curve.


% --- old text follows ---

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 $G^1$-continuity, determining a single tangent angle
at each control point suffices. Thus, combining a two-parameter curve
family with a technique 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. Parametrizing by chord length approximates the
MEC more closely but does not reduce the dimensionality of the
parameter space.

Similarly, the scale invariant MEC (dubbed SI-MEC in \cite{Moreton93})
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 of individual segments from the two of the pure MEC to
the three admitted by the elastica.


% --- begin old stuff

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
$(x_i, y_i)$ in order, the one that minimizes the MEC energy
functional, restated here:

\begin{equation}
E_{\mbox{\scriptsize MEC}}[\kappa(s)] = \int_0^l \kappa^2 ds
\end{equation}

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 them with 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 \cite{ForsytheLee}. A brief summary of their results follows:

\begin{itemize}
\item Each segment between knots is a MEC with $\lambda = 0.5$.
\item Curvature is continuous across knots ($G^2$-continuity).
\item Curvature is zero at endpoints (in the open case).
\item The spline is extensional.
\end{itemize}

For each segment between knots, the elastica with $\lambda$ fixed has
two remaining parameters, corresponding to the tangent angles at the
endpoints. Recall that Figure~\ref{mecrange-fig} shows instances over
the entire two-dimensional space in which such solutions exist, as
well as the boundary of that solution space. In particular, in the
case of symmetrical endpoint tangents, the curve is \emph{not} a
circular arc, so the MEC spline does not have the roundness
property. Forsythe and Lee found this unappealing, and propose adding
an arc length penalty term to the functional,

\begin{equation}
E[\kappa(s)] = \int_0^l \kappa(s)^2 + c\ ds\:.
\end{equation}

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

\section{Properties of the MEC spline}
\label{mec-props}

\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=0.35]{figs/quadrant.pdf}
\caption{\label{mec-roundness}MEC roundness failure.}
\end{center}
\end{figure*}

The most important property the MEC spline is lacking is
roundness. For three control points placed in an equilateral triangle
(shown in Figure~\ref{mec-roundness}, with the spline drawn in a thick
line, and the circle 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 deviation decreases to about $1.009$, and
asymptotically converges as $O(n^4)$ with the number of points.

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 $\lambda=0.5$
is $\pi/2$, and it is not much better for any of the non-symmetric
solutions.

\section{SI-MEC}
\label{simec}

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 S\'equin \cite{Moreton93}. This spline is
defined as the curve minimizing the SI-MEC functional

\begin{equation}
E[\kappa(s)] = l\int_0^l \kappa(s)^2 \ ds\: .
\end{equation}

The functional is termed ``scale invariant'' because its
value is invariant as the entire curve is scaled 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, in fact, 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, an appealing physical
interpretation is that of a mechanical spline with a bit of tension
pulling at the endpoints, making the overall curve slightly shorter.
Such a curve is generally \emph{not} two-parameter 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 than
the corresponding pure MEC \cite{Moreton93}.

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 $G^2$-continuity and
similar locality properties with the MEC.

% All of this text is really subsumed in Section \ref{kurgla1-sec}.
%\section{Piecewise SI-MEC}
%
%Another practical approach to improving the MEC is to use segments
%that each minimize the SI-MEC functional between any adjacent pair of
%control points (given the tangent constraints at the endpoints), and
%enforce $G^2$ continuity across the control points. A careful reading
%of Mehlum's description \cite{Mehlum74} shows that this is the spline
%approximated by his {\sc kurgla} 1 algorithm.
%
%The MEC, piecewise SI-MEC, and Euler spiral splines are all very
%closely related, and are essentially identical in the small-angle
%case. For larger bending angles at the control points, the difference
%is in the way that curvature varies:
%
%\begin{itemize}
%
%\item Curvature in the MEC increases linearly with the distance from
%  the line tangent to the inflection point.
%
%\item Curvature in the piecewise SI-MEC increases linearly along the
%  chord.
%
%\item Curvature in the Euler spiral spline increases linearly along
%  the arclength.
%
%*** really needs a figure ***
%
%\end{itemize}
%
%Since any dependence on the chord breaks extensibility (adding a new
%point changes the chords), the arclength is the most robust measure.
%For any given set of control points, the Euler spiral spline will have
%a slightly larger MEC functional, but almost all of that can be traced
%to the fact that the scaling of the MEC functional rewards a longer
%total curve length.

% TODO: should probably reinstate cubic lagrange section

\section{Cubic Lagrange (with length parametrization)}
\label{cubic-lagrange-sec}

The cubic Lagrange spline is defined as follows (see
\cite[p. 23]{Kvasov00} for more detail), given a sequence of input
data $(t_i, z_i)$, where $z_i$ can be interpreted as the vector $(x_i,
y_i)$:

\begin{equation}
L_{i,3}(t) = \sum_{j = i - 1}^{i + 2} z_j \frac{\omega_{i-1,3}(t)}{(t
  - t_j)\omega'_{i- 1, 3}(t_j)},\ \ \ \omega_{i-1,3}(t) = \prod_{j=i-1}^{i
  + 2}(t - t_j)
\end{equation}

To obtain the Lagrange polynomial normalized by chord lengths, set
successive $t_{i+1}$ so that $t_{i+1} - t_i = |z_{i + 1} - z_i|$. In
other words, the parametrization is the same as the chord length. (Set
$t_0 = 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 are small. 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 length parametrization.

Without chord length parametrization, the cubic Lagrange spline is
extensional, but prone to wild 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 the basis of the Ikarus spline, discussed in the
next section.

\section{Ikarus (biarc)}
\label{ikarus-sec}

The spline used by Peter Karow's Ikarus system \cite{Karow87} 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 elements of 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
\ref{cubic-lagrange-sec}) is fit through the
control points. Then the polynomial curve is discarded (but the
tangents preserved), and, for each segment, a curve defined by biarcs
is substituted. An arbitrary biarc has one additional degree of
freedom, corresponding 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 points
must be explicitly marked in the input.

% TODO: probably move the comparison to the fonts chapter.
\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=1.25]{figs/karow-a.pdf}
\caption{\label{ikarus-fig}Typical letter-shape using Ikarus.}
\end{center}
\end{figure*}

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 tangents are
symmetrical. Continuity is only $G^1$, 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 \cite{Karow91} recommends one control
point for every 30 degrees of arc.  A typical example
is shown in Figure \ref{ikarus-fig} (reproduced from
\cite{Karow91}). See also Figure \ref{ikarus-fig-compare} for a
comparison with the splines of this thesis.


\section{Euler's spiral}

One particularly well-studied two-parameter spline is Euler's spiral,
also known as the spiral of Cornu or the clothoid. More precisely, it is
the spline formed by joining piecewise Euler's spiral segments with
curvature ($G^2$) continuity at the knots.

\begin{figure*}[tb]
\begin{center}
\includegraphics[scale=0.8]{figs/cornu.pdf}
\caption{\label{euler-spiral-fig}Euler's spiral}
\end{center}
\end{figure*}


Euler's spiral, plotted in Figure~\ref{euler-spiral-fig} over the
range $[-5, 5]$, is most easily characterized by its intrinsic
definition: curvature is linear in arc length.

\begin{equation}
\kappa(s) = \frac{s}{2}\: .
\end{equation}

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

\begin{equation}
\label{fresnel-int-eq}
x(s) + iy(s) = \int_0^s e^{it^2} dt\: .
\end{equation}

This integral is also commonly expressed with a scale factor applied
to the parameter. The resulting curve is similar; its size is
$\sqrt{\frac{2}{\pi}}$ 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.

\begin{equation}
\label{fresnel-rescaled}
C(s) + iS(s) = \int_0^s e^{i\pi t^2 / 2} dt\: .
\end{equation}

\subsection{Euler's spiral as a spline primitive}

\begin{figure*}
\begin{center}
\includegraphics[scale=.7]{figs/clothmap}
\caption{\label{clothmap}Euler's spiral spline primitive over its parameter space.}
\end{center}
\end{figure*}

Euler's spiral is readily adapted for use as a two-parameter primitive
for a spline. Given tangent angles at the endpoints, find parameters
$t_0$ and $t_1$ from which to ``cut'' a segment from the spiral, so
that the resulting segment meets the endpoint constraints (for now,
assume that such parameters can be found; an efficient algorithm will
be given in Section~\ref{solve-euler}).

The behavior of this primitive is shown in
Figure~\ref{clothmap}. Compare to Figure~\ref{mecrange-fig}, showing the
corresponding 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.

\subsection{Variational formulation of Euler's spiral}

The spline literature has generally been split into two camps:
variational techniques, and the use of explicitly stated curves such
as the Euler spiral. Ironically, Euler's spiral was born of the study
of 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, showing that
the Euler spiral does indeed satisfy some fundamental variational
properties.

First, as noted by Kimia et al \cite{Kimia03}, 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:

\begin{equation}
E_{MVC}[\kappa(s)] = \int_0^l \Big(\frac{d\kappa}{ds}\Big)^2ds
\end{equation}

This result can be readily derived by treating the problem as an Euler
variational problem expressed in terms of $\kappa$, i.e. $F(x, y, y')
= y'^2$, with $s$ for $x$ and $\kappa$ for $y$ in the terms of
Equation~\ref{var-statement}, as stated in
Chapter~\ref{elastica-chapter}. The corresponding Euler--Lagrange equation is
particularly simple:

\begin{equation}
\kappa'' = 0
\end{equation}

The solution family is thus $\kappa = k_0 + k_1s$, which are simple
similarity transformations of the basic 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 $\kappa$ rather than $\theta$, the boundary
conditions are written in terms of curvature rather than tangent
angle. Similarly, the additional terms corresponding to
Equation~\ref{elastica-pos} representing the {\em positional}
endpoint constraints are missing. Recall that without these
constraints, the only solution to the MEC equation is a straight line.

% todo - needs fixing
So this result only tells us that the Euler spiral is the MVC-optimum
curve satisfying {\em curvature} constraints on the endpoints, but
with unconstrained endpoint positions and tangent angles. As we shall
explore more deeply in Section~\ref{sec-eulerpoisson-mvc}, 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 of piecewise Euler spiral segments has $G^2$-continuity,
the same as a MEC spline, while the solution to the MVC functional
with positional constraints for internal knots has $G^4$-continuity.

So the Euler spiral is a reasonably good first-order approximation to
the MEC functional, and is a special case in the solution space of the
MVC functional, but might there be a functional for which it is
actually \text{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 which happens
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 over the entire arc length,

\begin{equation}
E_{\mbox{\scriptsize MAXVC}}[\kappa(s)] = \max\Big
|\frac{d\kappa}{ds}\Big |\: .
\end{equation}

This functional is related to the MVC; while the MVC is the 2-norm of
curvature variation, $E_{\mbox{\scriptsize MAXVC}}$ corresponds to the
$\infty$-norm.

\subsection{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~\ref{euler-c-example}. Estienne
is a distinguished book face designed by the eminent English printer
George W. Jones, and released by Linotype in
1930 \cite[p. 133]{McGrew93}. 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, the sharp 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 to zero crossings on the curvature plot and marked on
the outline with short lines.

\begin{figure*}[tb]
\begin{center}
\includegraphics[scale=0.8]{figs/c_g2.pdf}
\caption{\label{euler-c-example}Closed Euler example.}
\end{center}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[scale=0.8]{figs/g_g2.pdf}
\caption{\label{euler-g-example}Complex Euler example.}
\end{center}
\end{figure*}

A more complex example is shown in Figure~\ref{euler-g-example}, 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
$G^2$-continuous smooth knots and ``corner points'', with only
$G^0$-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. Note that sections with no smooth
knots are straight lines (there is one at the top of the ``ear''), and
sections with a single smooth knot are circular arcs (two of these
make up the rest of the ear). More generally, each spline segment
bordered by a corner point is a circular arc, with zero variation of
curvature.

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 an
interactive editor, the correspondence between the input (the
placement of knots and control points) seems considerably more
intuitive than the standard B\'ezier curve editing techniques.

\begin{figure*}
\begin{center}
\includegraphics[scale=0.8]{figs/pure_u.pdf}
\caption{\label{pure-u}Difficulty with straight-to-curved joins.}
\end{center}
\end{figure*}

The trickiest curves to draw using pure Euler's spiral splines are
transitions between straight lines and curves, as typical in the shape
of a ``U.'' The difficulty is illustrated in
Figure~\ref{pure-u}. Essentially, 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 the curved 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 \ref{sec-constraints}.

\subsection{Existence and uniqueness of Euler spiral spline}
\label{euler-existence-sec}

One of the serious shortcomings of the MEC spline is lack of existence
for all configurations of control points (see Section
\ref{props-existence-sec} 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.

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
\cite{Stoer82}.\footnote{Although, the reader is warned, this paper
  contains numerous errors, including in the statement of the crucial
  Theorem 2.25.} Stoer considered a slightly weaker formulation of the
Euler spiral spline: each segment is piecewise an Euler spiral
segment, and there is $G^2$-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 sequence of control points, there exist
infinitely many assignments of Euler spiral parameters such that the
resulting segments have $G^2$-continuity.

The proof is dependent on the following theorem:

\begin{quote}
\textbf{Stoer's theorem 2.25} (restated): For all points $P_0$ and
$P_1$, tangent angles $\theta_0$, and curvatures $\kappa_0$, there
exist countably many different Euler spiral segments connecting $P_0$
to $P_1$ such that the tangent at $P_0$ is $\theta_0$ and the
curvature at $P_0$ is $\kappa_0$.
\end{quote}

\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=1.0]{figs/stoer_fig6.pdf}
\caption{\label{stoer-fig6}Multiple windings for Euler spiral segment
  connecting two points (from \cite{Stoer82}).}
\end{center}
\end{figure*}

Stoer's proof is based on the fact that an Euler spiral segment with
gently varying curvature may wind around arbitrarily many times
before reaching $P_1$. In the limit, the total winding number of the
spiral segment increases by approximately $2\pi$ for each successive
solution. Such multiple windings are shown in Figure \ref{stoer-fig6},
reproduced from \cite[p. 329]{Stoer82}.

The result that there are infinitely many piecewise-Euler spiral,
$G^2$-continuous interpolating splines follows readily from Stoer's
theorem 2.25. Start at the first point, choosing arbitrary tangent and
curvature, then choose one of the solutions guaranteed by the theorem,
and repeat until the final point is reached.

With arbitrary choices, such a procedure will tend to produce splines
with excessive winding. Stoer proposed using the MEC energy to select
a single ``best'' spline from the myriad alternatives. Stoer also gave
an algorithm based on the standard band-diagonal Jacobian matrix, and
claimed that it gave good results in general, but did not include a
proof that the algorithm achieves a global optimum in all cases.

In addition, using the global minimum of MEC energy will also have an
effect on the curvature at the endpoints. Since minimizing MEC over
all curves assigns zero curvature to endpoints (see Section
\ref{mec-interp-sec}), 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 methods of
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
\ref{sec-global-solver} converge on a solution (and a solution with no
excess windings, as well) for reasonable input sequences of control
points, but can fail to converge when turning angles become
excessive.

\section{Hobby's spline}
\label{hobby-sec}

A significant interpolating spline, especially in font design, is
Hobby's spline as implemented in Metafont \cite{Hobby85}. It is perhaps
easiest to understand this spline as a very computationally efficient
approximation to the Euler spiral spline of the previous section.

The definition of the primitive segment, with endpoint tangent angles
$\theta_0$ and $\theta_1$, is

\begin{equation}
\begin{array}{l}
a = \sqrt{2}, \\
b = \frac{1}{16}, \\
c = (3 - \sqrt{5}) / 2, \\
\alpha = a(\sin\theta_0 - b\sin\theta_1)(\sin\theta_1 -
b\sin\theta_0)(\cos\theta_0-\cos\theta_1), \\
\rho = \frac{2 + \alpha}{1 + (1-c)\cos\theta_0 + c\cos\theta_1}, \\
\sigma = \frac{2 - \alpha}{1 + (1-c)\cos\theta_1 + c\cos\theta_0}, \\
x_0, y_0 = 0, 0, \\
x_3, y_3 = 1, 0, \\
x_1, y_1 = \frac{\rho}{3\tau}\cos\theta_0,
\frac{\rho}{3\tau}\sin\theta_0, \\
x_2, y_2 = 1 - \frac{\sigma}{3\tau}\cos\theta_1,
\frac{\sigma}{3\tau}\sin\theta_1.
\end{array}
\end{equation}

The $(x_i, y_i)$ are the coordinates for a standard cubic B\'ezier
curve. The adjustable parameter $\tau$ 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' \cite[p. 311]{Knuth:1986:CMT}.

\begin{figure*}
\begin{center}
\includegraphics[scale=.7]{figs/hobby}
\caption{\label{hobby-fig}Hobby spline primitive over its parameter space.}
\end{center}
\end{figure*}

The behavior of the Hobby spline primitive over its parameter space
can be seen in Figure~\ref{hobby-fig} (see Figure~\ref{mecrange-fig}
for comparison to the MEC). It is defined for all
values of $\theta_0$ and $\theta_1$, but the figure only shows values
up to $2.156$ radians (123.5 degrees).

\begin{figure*}
\begin{center}
\includegraphics[scale=.7]{figs/hobbyk}
\caption{\label{hobbyk-fig}Curvature plot of Hobby spline over its
  parameter space.}
\end{center}
\end{figure*}

Similarly, the curvature plot is shown in Figure~\ref{hobbyk-fig}. 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 $G^2$-continuous at the control points, the
Hobby spline uses a slightly different criterion. Hobby introduces the
concept of \emph{mock curvature}, which is the linear approximation
(taking only the first-order term of the exact equation) of endpoint
curvature as a function of the endpoint angles $\theta_0$ and
$\theta_1$. 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
\emph{linear} problem. In fact, since each tangent angle affects two
neighboring segments, and each segment in turn contributes to the mock
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 $\tau$, including $\tau=1$. Thus, Hobby's is one
of the few ``advanced'' splines for which there is a solution for
\emph{all} configurations of 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(n^6)$
as opposed to the MEC, which is only $O(n^4)$. Similarly, the
extensionality and continuity of curvature are only approximate, and
only good at small angles. However, as noted above, it is an
exceedingly 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 used in Potrace \cite{potrace03}. While Potrace is
not an interpolating spline (rather, it is a tool for automatically
converting 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 B\'ezier curves
(inherently four parameters), with mathematical equations to determine
the additional two parameters.

\section{Splines from fixed generating curves}
\label{generator-curve-sec}

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 match the 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 \ref{extensional-equivalence}, 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, $\lim_{s \rightarrow \infty} \kappa'(s) / \kappa(s)^2 =
0$. The Euler spiral meets this condition because $\kappa$ tends to
infinity, while $\kappa'$ is fixed. Another possibility would be a
curve tending to fixed curvature in the asymptote, in which case
$\kappa'$ would tend to $0$. The MEC, by contrast, clearly does not
have this property.

The envelope of endpoint chord angles (as plotted in
Figure~\ref{mecrange-fig} for the MEC) also has an effect. A small
envelope yields failure to exist. A generating curve without an
inflection (such as a parabola) excludes all antisymmetric endpoint
angle configurations from this envelope, even for tiny angles.

\subsection{Log-aesthetic curve splines}
\label{sec-aesthetic}

A very promising family of candidates for generator curves is the
\emph{aesthetic curve,} also known as the \emph{log-aesthetic
  curve} \cite{Miura05}. These curves are defined as having a straight
line when $\kappa'/\kappa$ is plotted against $\kappa$ on a log-log
scale. In addition to a number of special cases (including the
equiangular spiral and Nielsen's spiral), the general formula for
aesthetic curves is $\kappa = s^\alpha$, where the exponent $\alpha$
is the parameter that distinguishes the members of this family. An
intuitive motivation for the aesthetic 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.

\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=0.6]{figs/aesthetic.pdf}
\caption{\label{aesthetic-examples-two}Test instrument for fairness user study.}
\end{center}
\end{figure*}

The literature on aesthetic curves doesn't tend to include curves with
inflection points -- on a log-log plot, these are off the
scale. However, since it is well known that the aesthetic curve with
$\alpha=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 $\kappa=\mbox{sign}(s)|s|^\alpha$ 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
\ref{aesthetic-examples-two}.

As described in Section
\ref{empirical-fairness-sec}, an empirical user study establishes that
MEC energy is not an accurate predictor of perceived fairness. MEC
energy is minimized at $\alpha=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 $\alpha=1.78$, that hypothesis is refuted.

% note: hand-cropped to remove whitespace
\begin{figure*}[tbh]
\begin{center}
\includegraphics[scale=0.6]{figs/locality.pdf}
\caption{\label{aesthetic-locality}Perturbation fall-off rate as a
  function of aesthetic curve exponent.}
\end{center}
\end{figure*}

The exponent of the aesthetic curve used as the generator curve in a
spline also has a profound effect 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 for
a single perturbed point, a quantitative measurement of locality, is
plotted against the exponent of the aesthetic curve in Figure
\ref{aesthetic-locality}. For $\alpha=1$, which represents the Euler
spiral, the amplitude 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 approximation and solving for $\kappa_0\theta_1 =
\kappa_1\theta_0$ for the endpoints of a segment, which yields a value
of $1/(2+\sqrt{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.

\section{Piecewise SI-MEC splines, or Kurgla 1}
\label{kurgla1-sec}

The Autokon system of Mehlum \cite{Mehlum74} was inspired by the goal
of simulating a real, mechanical spline, but its {\sc kurgla} 1 algorithm
actually implemented something a bit different, with both advantages
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 from the line of action. His solver approximated
the integration of this intrinsic equation as a sequence of 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 angle of the line of action with respect to the
control point chord, and one for the amplitude of the curvature
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 {\sc kurgla} 1
did not attempt it.

In particular, it nailed down the \emph{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 tangents at the endpoints, and the remaining global problem of
choosing tangents for $G^2$-continuity yields to simple Newton-style
solution.

As shown in Section \ref{simec-section}, an elastica segment with the
line of action perpendicular to the chord is the one minimizing the
SI-MEC energy for the angle constraints at those endpoints. Thus, the
{\sc kurgla} 1 algorithm effectively computes a piecewise SI-MEC.

One consequence is that {\sc 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 consequence is considerably better
robustness than the pure MEC, because the envelope of parameter space
is approximately 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, {\sc 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 {\sc kurgla} 2 algorithm is the earliest known
implementation of the Euler spiral spline.\footnote{Mehlum wrote to
  George Forsythe on July 9, 1971, mentioning ``Fresnel integral
  splines'' as recent work. The work was published in 1974
  \cite{Mehlum74}.}

\begin{figure*}[thb]
\begin{center}
\includegraphics[scale=0.8]{figs/mec_euler_compare.pdf}
\caption{\label{mec-euler-compare}MEC, SI-MEC, and Euler spiral compared.}
\end{center}
\end{figure*}

The MEC, SI-MEC, and Euler spiral are closely related. One way to
visualize this relation is 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 \ref{mec-euler-compare} shows these gradients of
curvature visually -- along each of the curves, tick marks meter the
rate of change of curvature, and the tick mark is angled so that 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. 

\section{Circle splines}
\label{cspline-sec}

\begin{figure*}[thb]
\begin{center}
\includegraphics[scale=0.6]{figs/cspline.pdf}
\caption{\label{cspline-fig}Circle spline extensionality failure.}
\end{center}
\end{figure*}

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
S\'equin, Lee, and Yin \cite{DBLP:journals/cad/SequinLY05}.

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 $\kappa'$ at the endpoints,
and that the curvature is identical to that of the circle used to
determine tangents, so the resulting spline has $G^3$-continuity.

Unfortunately, even though the spline has good locality and
$G^3$-continuity, it is not particularly smooth for all inputs, and
adding additional on-curve control points can cause serious deviations
from extensionality, in many cases adding inflections. Figure
\ref{cspline-fig} clearly shows the lack of extensionality; the two
paths differ only by the addition of the central point in the
right-hand figure.

\section{Summary}

This 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 of a 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 most common such metric, bending energy, does not coincide
accurately with user preference.

\newcommand{\rateexc}{$+++$}
\newcommand{\ratevg}{$++$}
\newcommand{\rategood}{$+$}
\newcommand{\ratefair}{$-$}
\newcommand{\ratepoor}{$--$}
\newcommand{\ratebad}{$---$}
\newcommand{\sml}{\footnotesize}

\vspace{5mm}
\begin{table}[ht]
\label{twoparam-table}
%\[
\begin{tabular}{ll|lcccc}
% \sml in every cell is inelegant, but gets the job done.
\sml \textrm{Curve type} & \sml \textrm{Tangents} & \sml \textrm{Continuity}
& \sml \textrm{Extensionality} & \sml \textrm{Roundness} & \sml \textrm{Robustness} &
\sml \textrm{Locality} \\
\hline
\sml MEC & \sml $G^2$ & \sml $G^2$ & \sml \rateexc & \sml \ratepoor & \sml \ratepoor & \sml \rategood \\
\sml Euler spiral & \sml $G^2$ & \sml $G^2$ & \sml \rateexc & \sml
\rateexc & \sml \ratevg & \sml \rategood \\
\sml log-aesthetic & \sml $G^2$ & \sml $G^2$ & \sml \rateexc & \sml \rateexc & \sml \ratefair & \sml \ratevg
\\
\sml Hobby & \sml $\approx G^2$ & \sml $\approx G^2$ & \sml \ratefair & \sml \ratefair & \sml \rateexc &
\sml \rategood \\
\sml Biarc (Ikarus) & \sml Lagrange & \sml $G^1$ & \sml \ratepoor & \sml \rateexc & \sml \ratevg &
\sml \rategood \\
\sml {\sc kurgla} 1 & \sml $G^2$ & \sml $G^2$ & \sml \ratefair & \sml \rateexc & \sml \rategood &
\sml \rategood \\
\sml Circle & \sml local & \sml $G^3$ & \sml \ratebad & \sml \rateexc & \sml \rateexc & \sml \rateexc
\end{tabular}
%\]
\caption{Comparison of properties of two-parameter splines.}
\end{table}

Table \ref{twoparam-table} summarizes the properties of several splines in this
two-parameter family. All of these are known and appear in the
literature.

\begin{itemize}
\item 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 $G^2$
  continuity is preserved across each control point. (Section \ref{mec-interp-sec})

\item Euler spiral is the interpolating spline for which each segment
  between control points is a section of an Euler spiral,
  i.e. curvature is linear in arc length. Like the MEC, the tangents
  are chosen for $G^2$-continuity. (Section \ref{euler-interp-sec})

\item log-aesthetic is the two-parameter extensible spline based on
  log-aesthetic curves with an exponent near $\alpha=2$. (Section
  \ref{sec-aesthetic})

\item Ikarus (Section \ref{ikarus-sec}).

\item Hobby is a piecewise cubic spline developed for the Metafont
  system. It is best understood as a polynomial approximation to the
  Euler spiral spline, with better robustness but lower performance on
  other metrics such as smoothness and extensionality \cite{Hobby85}
  (Section \ref{hobby-sec}).

\item {\sc kurgla} 1 (Section \ref{kurgla1-sec}).

\item Circle spline (Section \ref{cspline-sec}).

\end{itemize}

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.
