\chapter{Converting to B\'ezier curves}
\label{chap-tobez}
By far the most popular representation of curves is piecewise
polynomial, most often with the coefficients represented in
Bernstein--B\'ezier form. This popularity is well earned -- the
representation is simple, numerically stable, and the primitive is
quite supple, capable of representing a wide range of curves with
relatively few control points. Both cubic and quadratic B\'eziers are
widely used. The PostScript language \cite{PLRM1} established cubic
B\'eziers as a de facto standard for 2D graphics, and the TrueType
font format and Macromedia Flash are based on quadratics, presumably
because the rasterizer is even simpler to implement than for cubics.
There are several reasons to convert piecewise clothoid curves, or
piecewise polynomial spiral curves, into B\'ezier 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 B\'ezier curves, but few, if any, can handle
piecewise clothoids.
The conversion to B\'ezier form involves a tradeoff between the
conciseness of the result and the amount of computation needed. For
some applications, such as interactive display, the curve will be
computed, displayed, and then discarded. A somewhat wasteful B\'ezier
representation is perfectly acceptable as long as the overall
computation time is small. For other applications, the curve will be
computed, saved into a file such as a font, then embedded thousands of
times in document files, each of which in turn may be rendered
thousands of times. In such applications, a significant investment in
computing time to reduce the size of the B\'ezier representation can
pay off handsomely. This section will describe both fast and cheap
methods for converting polynomial spirals into B\'eziers, as well as
much slower techniques, suitable for creating optimally concise
B\'ezier representations.
There are several techniques in the literature for creating polynomial
approximations of clothoid curves. S\'anchez-Reyes and
Chac\'on \cite{Sanchez-ReyesC03} use a family of polynomials called
\emph{s-power series,} which is a fairly powerful and general
technique. However, the error in their technique converges only as
$O(n^4)$ in the number of subdivisions, while convergence of $O(n^6)$
is possible for cubic B\'eziers.
\section{Parameters}
Assuming fixed endpoints, quadratic B\'eziers have two free
parameters. These are exactly sufficient to determine the tangent
angles at the endpoints. Thus, the problem of generating a quadratic
B\'ezier 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.
\begin{figure*}
\begin{center}
\includegraphics[scale=.6]{figs/simplebez}
\caption{\label{simplebez}A cubic B\'ezier.}
\end{center}
\end{figure*}
Cubic B\'eziers 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 B\'eziers are often described in terms of two \emph{control
arms}, or vectors from the endpoints to the control points. In the
typical B\'ezier of Figure \ref{simplebez}, these vectors are $(x_0,
y_0)$ to $(x_1, y_1)$ and $(x_3, y_3)$ to $(x_2, y_2)$, 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.
\section{Error metrics}
Fitting of a cubic B\'ezier to a target curve can be seen as a simple
optimization problem. Over the 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 used metric. However, this metric has several
shortcomings. To address these shortcomings, this section presents
other error metrics based on angle, rather than distance.
\begin{figure*}
\begin{center}
\includegraphics[scale=.8]{figs/quantcircle}
\hspace{10mm}
\includegraphics[scale=.8]{figs/warpedcircle}
\caption{\label{errorcompare}Two circle approximations with equal
distance error.}
\end{center}
\end{figure*}
Figure \ref{errorcompare} 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 is
smooth 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 B\'ezier approximation and the target
curve have identical total arc length $l$. Then, we propose the
\emph{absolute angle metric}
\begin{equation}
E_{\mbox{\scriptsize abs}} = \int_0^l |\theta_{\mbox{\scriptsize approx}}(s) - \theta(s)|\; ds\:.
\end{equation}
Note 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 upper bound
on the flatness. Consider the distance between points at arc length $d
\leq l$
from the curve start, here using complex coordinates for algebraic
convenience:
\begin{equation}
\label{int_abs_error}
z_{\mbox{\scriptsize approx}} - z = \int_0^d e^{i\theta_{\mbox{\scriptsize approx}}(s)} - e^{i\theta(s)}\; ds\:.
\end{equation}
But note the inequality
\begin{equation}
|e^{i\theta_{\mbox{\scriptsize approx}}(s)} - e^{i\theta(s)}| \leq |\theta_{\mbox{\scriptsize approx}}(s) - \theta(s)|\:.
\end{equation}
Note also that the integral \ref{int_abs_error} is bounded from above
by $E_{\mbox{\scriptsize abs}}$ because the absolute angle error is everywhere positive.
Thus, for all points on the original curve, there exists a point on
the approximate curve no more distant than $E_{\mbox{\scriptsize abs}}$ from this point,
and thus from the target curve.
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 \emph{square} of the angle
error, integrated over the arc length as before.
\begin{equation}
E_{\mbox{\scriptsize sq}} = \Big (\int_0^l (\theta_{\mbox{\scriptsize approx}}(s) - \theta(s))^2\; ds \Big )^{1/2}\:.
\end{equation}
From standard properties of norms, $\sqrt{l} E_{\mbox{\scriptsize sq}} \geq
E_{\mbox{\scriptsize abs}}$
holds, which means that $E_{\mbox{\scriptsize sq}}$ also places a strong upper bound on
the flatness error metric. This property implies that minimizing
the $E_{\mbox{\scriptsize sq}}$ metric will also yield a curve with a good flatness
metric as well, although the converse is not always true.
\begin{figure*}
\begin{center}
\includegraphics[scale=.35]{figs/draw_cornu_9e4}
\hspace{10mm}
\includegraphics[scale=.35]{figs/draw_cornu_1e6}
\caption{\label{l2-angle-fig}Euler spiral approximation drawn using
.03 and .001 error threshold, respectively.}
\end{center}
\end{figure*}
Figure \ref{l2-angle-fig} shows examples of this error metric in
use. Both approximations of the Euler spiral segment use a number of
cubic B\'ezier curves, each with the approximately the same
$E_{\mbox{\scriptsize sq}}$ error metric. Each segment on the left has
an error of approximately 0.03, and each segment on the right
approximately 0.001 (relative to the total size of the figure being
about 1). As can be seen, an error of 0.03 results in a somewhat
distorted curve, while 0.001 is essentially indistinguishable from the
original curve. This figure also provides evidence that converting to
B\'ezier form by solving an optimization problem in terms of the
$E_{\mbox{\scriptsize sq}}$ metric gives good results, even for a
relatively crude error threshold.
\section{Simple conversion}
A simple, fast conversion to cubic B\'ezier form is as follows. Like
the s-power series, and like quadratic B\'eziers, the convergence is
only $O(n^4)$, 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 the curve. An example, fitting the section of Euler
spiral from parameter $t = 0.5$ to $t = 1.1$, is shown as the leftmost
instance in Figure \ref{threebez}. Both the B\'ezier (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.
\begin{figure*}
\begin{center}
\includegraphics[scale=.35]{figs/fastbez}
\includegraphics[scale=.35]{figs/eqarcbez}
\includegraphics[scale=.35]{figs/goodbez}
\caption{\label{threebez}Three B\'ezier solutions: fast, equal arc length, and fully optimized. Curvature plot is the top graph, the B\'ezier is on bottom.}
\end{center}
\end{figure*}
\section{Equal arc length, equal arm length solution}
\label{sec-equal-arclength}
A somewhat more accurate solution can be obtained by solving the
length of both control arms so that the total arc length of the
B\'ezier matches that of the target curve. An example is shown as the
top right instance of Figure \ref{threebez}.
There is no simple formula for the total arc length of a
B\'ezier, so for this solution we rely on numerical integration. Since
the relation between the control arm length and the total arc length is
smooth and monotonic, a simple numerical solving technique such as the
secant method gives robust, accurate results in a few iterations.
\section{Fully optimized solution}
\label{sec-full-opt}
In the fully general case, the lengths of the control arms need not be
equal. To find the single best-fitting cubic B\'ezier 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 \ref{isolines}. This plot is the root mean squared angle error for all
combinations of left and right control arm length, with a target 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 of these
two local minima are shown in Figure \ref{localminima}. 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 \ref{localminima2} shows a fit of a different section
of the Euler spiral, from $t = 0.5$ to $t = 0.85$, and the two local
minima are nearly indistinguishable, both in overall error metric and
in curvature plot.
\begin{figure*}
\begin{center}
\includegraphics[scale=.8]{figs/isolines}
\caption{\label{isolines}Approximation errors over parameter space for cubic B\'ezier fit.}
\end{center}
\end{figure*}
\begin{figure*}
\begin{center}
\includegraphics[scale=.8]{figs/localminima}
\caption{\label{localminima}Local minima of cubic B\'ezier fitting,
t=[0.5, 1.0]. Curvature plots are the top two graphs, the B\'eziers
are on the bottom.}
\end{center}
\end{figure*}
\begin{figure*}
\begin{center}
\includegraphics[scale=.8]{figs/localminima2}
\caption{\label{localminima2}Local minima of cubic B\'ezier fitting,
t=[0.5, 0.85]. Curvature plots are the top graphs, the B\'eziers are
on the bottom.}
\end{center}
\end{figure*}
Also plotted in Figure \ref{isolines} is the curved line representing all
values of left and right control arm length so that the total
arc length of the B\'ezier 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 B\'ezier matches that of the target curve.
As can be seen in the rightmost example in Figure \ref{threebez}, the
fully optimized solution finds a B\'ezier that matches the original
Euler spiral segment with impressive accuracy. Even in the curvature
plot, the deviation from ideal is barely perceptible.
\section{Global optimum}
Generating an optimized B\'ezier 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 B\'ezier representation of a curve has three
elements:
\begin{itemize}
\item The number of subdivisions is minimal with respect to the error
metric threshold.
\item The position of the subdivision points minimizes the error
metric.
\item Each individual segment minimizes the error metric given fixed
endpoints.
\end{itemize}
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 of individual
B\'ezier segments (with computation of error metric) as an inner loop.
As before, there is a tradeoff between the amount of time spent
optimizing the B\'ezier representation and the degree to which the
representation actually meets these goals. For applications such as
quick interactive display, where the B\'ezier will be drawn once, it
makes sense to compute it as quickly as possible, even if the number
of segments is substantially above the minimum. For 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.
\begin{figure*}
\begin{center}
\includegraphics[scale=.5]{figs/cecco_a}
\includegraphics[scale=.5]{figs/cecco_a_3_2}
\caption{\label{optbez1}Conversion into optimized cubic B\'eziers.}
\end{center}
\end{figure*}
Figure \ref{optbez1} shows an example of an optimized conversion into
cubic B\'ezier 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
PostScript or SVG (Scalable Vector Graphics), this requirement could
be relaxed.
\subsection{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
$L^2$-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^\infty$-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^\infty$-norm, this is equivalent to the purest
approach. However, the $L^2$-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 B\'ezier fit. Given this property and the objective of
minimizing an $L^\infty$-norm, all optimum B\'ezier representations of
curves have the property that each individual B\'ezier segment has an
identical error metric.
In theory, it would be possible to optimize for the pure $L^2$-norm,
but in practice it would be computationally much more expensive, and
the results would be almost imperceptibly different. For reasoning
about this approach, it is perhaps simplest to accept the $L^2$-norm
metric as a reasonably good approximation to the $L^\infty$-norm.
\subsection{Simple subdivision}
The simplest technique for meeting a particular error metric is to
simply fit a single B\'ezier 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^\infty$-norm, then the same
error threshold can be applied to both subdivisions. If the error
metric is an $L^2$-norm, then the threshold for each subdivision is
$\sqrt{0.5}$ of the original threshold.
One refinement is to subdivide at the midpoint of angle traversed,
rather than at the midpoint of 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, a more sophisticated approach is needed,
using the integral of the absolute value of curvature as the
quantity 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 B\'ezier segments, as
described in Section \ref{sec-full-opt}, results in 34 segments. Using
the less computationally expensive equal arc length approach as
described in Section \ref{sec-equal-arclength} increases the count to
47. The conversions are shown in the upper right and upper left
quadrants, respectively, of Figure \ref{optbez2}.
\subsection{Smarter subdivision}
Simple subdivision often results in more segments than absolutely
needed. In order to produce an actually optimum B\'ezier
representation, we propose a smarter subdivision technique. It
consists of four nested loops. In order of outer to inner:
\begin{itemize}
\item Choose the number of subdivisions needed (a simple approach is
to count up from 1, and stop when the error threshold is met).
\item Start with an initial dumb subdivision (equal arc length is
fine), then iteratively refine by choosing the pair of adjacent
segments with the greatest discrepancy of error, and finding the
subdivision point joining those segments so that the error is equal.
\item A slightly modified secant technique selects the new subdivision
point so that the error metric of the two segments is
balanced. Simple bisection would also work but its convergence isn't
quite as fast.
\item The inner loop computes the optimized B\'ezier segments (along
with the error metric) on either side of the control point, using
the algorithms of Section \ref{sec-full-opt}.
\end{itemize}
\begin{figure*}
\begin{center}
\includegraphics[scale=.5]{figs/cecco_a_0_2} \hspace{4mm}
\includegraphics[scale=.5]{figs/cecco_a_1_2}
\vspace{4mm}
\includegraphics[scale=.5]{figs/cecco_a_2_2} \hspace{4mm}
\includegraphics[scale=.5]{figs/cecco_a_3_2}
\caption{\label{optbez2}Four different levels of optimization:
subdivision in half with equal arm length, subdivision in half with
fully optimized B\'ezier segments, greedy subdivision with fully optimized
segments, and fully optimized.}
\end{center}
\end{figure*}
Figure \ref{optbez2} 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 B\'eziers are chosen to have equal arc length as the true
curve, but with equal length control arms, as in Section
\ref{sec-equal-arclength}. In the top right, each B\'ezier segment is
fully optimized, as in Section \ref{sec-full-opt}.
In the global optimum (the lower right corner of Figure \ref{optbez2}),
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.
\subsection{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 using slightly 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, this technique may be
faster in cases when only a minimal number of segments is needed
(meeting the error 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 each segment is the same. Given the error metric
target, the algorithm walks from one endpoint of the curve along it,
until the end is reached, one segment at a time. For each segment, the
left endpoint is already known (again, starting at the endpoint of the
curve). Bisection (or some other variant of 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 the last segment is likely to have a much
lower error metric than the target (and, thus, be shorter than
needed). An example is shown in the lower left corner of Figure
\ref{optbez2}. Note that the total number of segments (30) is equal to
the fully optimized case in the lower right corner. Yet, some of the
subdivision points seem rather lopsided. For example, the subdivision
point on the top curve of the bowl 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 a fully 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 be significantly 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, and often slower. This algorithm requires 23.3
seconds (in a slow Python implementation) to produce the lower left
corner of Figure \ref{optbez2} 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 discrepancy is
even greater. In settings where the absolute minimum number of
B\'ezier 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 the upper left corner, 47 segments in 0.9 seconds.
Such times are still undesirably large for interactive use. Clearly,
reimplementation in a faster language such as C would improve the
compute time by at least an order of magnitude, and other numerical
improvements may also be possible. But as the last phase in a font
designing process that lasts many hours, even these speeds are
acceptable. The next chapter gives discusses applications of
interpolating splines for the production of fonts in more detail.