\section{Determining Euler spiral parameters}
\label{solve-euler}

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 determination of 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 \cite{Kimia03}. Their
solution represents the curve in terms of two parameters:

\begin{equation}
\kappa(s) = \gamma s + \kappa 0\mbox{, with}\ 0 \leq s \leq L
\end{equation}

Numerically, their solution is complicated. They first construct a
biarc to approximate the desired 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 from the
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 the angle constraint.

\begin{figure*}
\begin{center}
\includegraphics{figs/solvecloth}
\caption{\label{solvecloth}Determining Euler spiral parameters.}
\end{center}
\end{figure*}

We represent the curve to solve as

\begin{equation}
\kappa(s) = k_0 + k_1 s\mbox{, with} -0.5 \leq s \leq 0.5\:.
\end{equation}

Integrating, and fixing $\theta(0) = 0$,

\begin{equation}
\label{cloth-theta}
\theta(s) = k_0 s + \frac{k_1}{2} s^2\:.
\end{equation}

Again, the curve can be given in parametric form simply by integrating
the angle:

\begin{equation}
\label{cloth-cesaro}
x(s) + i y(s) = \int_0^s e^{i\theta(t)}\ \mbox{d}t\:.
\end{equation}

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

\begin{align}
\label{cloth5}
\theta_0 &= \psi - \theta(-0.5)\:, \\
\theta_1 &= \theta(0.5) - \psi\mbox{, where} \\
\label{cloth5c}
\psi &= \arctan\frac{y(0.5) - y(-0.5)}{x(0.5) - x(-0.5)}\:.
\end{align}

But the formula for $\theta(s)$ is known (Equation \ref{cloth-theta}), so

\begin{eqnarray}
\label{cloth6}
\theta_0 = \frac{k_0}{2} - \frac{k_1}{8} + \psi\:, \\
\label{cloth6a}
\theta_1 = \frac{k_0}{2} + \frac{k_1}{8} - \psi\:.
\end{eqnarray}

From Equations \ref{cloth6} and \ref{cloth6a} we can immediately solve for $k_0$:

\begin{equation}
k_0 = \theta_0 + \theta_1\:.
\end{equation}

We have now reduced the search space to a single dimension, $k_1$. We
will use a simple polynomial approximation to derive an initial guess,
an approximate Newton--Raphson method to derive a second guess, and the
secant method to further refine the approximation.

Using the small angle approximation $e^{i\theta} \approx 1 + i\theta$,
Equation \ref{cloth-cesaro} becomes

\begin{equation}
\begin{split}
\label{cloth7}
x(s) + i y(s)& \approx \int_0^s 1 + i\theta(t)\ \mbox{d}t \\
& = s + i\Bigl (\frac{k_0}{2} s^2 + \frac{k_1}{6} s^3\Bigr )\:.
\end{split}
\end{equation}

Combining Equations \ref{cloth5c} and \ref{cloth7} (note the
cancellation of even terms, just as for the integration in the
previous subsection),

\begin{eqnarray}
x(.5) - x(-.5) \approx 1\:, \\
y(.5) - y(-.5) \approx \frac{k_1}{24}\:, \\
\psi \approx \arctan\frac{k_1}{24}\:.
\end{eqnarray}

Applying the small angle assumption $\arctan \theta \approx \theta$
and substituting Equation \ref{cloth6},

\begin{eqnarray}
\label{solveeuler-linear-approx}
\theta_0 \approx \frac{k_0}{2} - \frac{k_1}{12}\:, \\
\theta_1 \approx \frac{k_0}{2} + \frac{k_1}{12}\:.
\end{eqnarray}

Thus, we can derive a good approximation to $k_1$:

\begin{equation}
\label{fit-euler-approx}
k_1 \approx 6(\theta_1 - \theta_0)\:.
\end{equation}


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 one-dimensional
functions. In this application, the function to be solved is angle
error as a function of $k_1$, with $\theta_0$ and $\theta_1$ given.

The secant method needs two initial guesses. For the first, we choose
$k_1 = 0$, simply because the angle error is trivial to calculate at
that value. For the second guess, the approximation of equation
\ref{fit-euler-approx} works well for $\theta_0 + \theta_1 < 6$ or so,
but as this value approaches $2\pi$, it significantly overestimates
$k_1$ and sometimes fails to converge. The guess $6(\theta_1 -
\theta_0)(1 - \frac{\theta_0 + \theta_1}{2\pi})^3$ has been found to
guarantee convergence for all $\theta_0 + \theta_1 < 2\pi$, and is
also slightly more accurate, thus requiring fewer iterations to
converge to a given precision.

\begin{figure*}
\begin{verbatim}
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
\end{verbatim}
\caption{\label{solvecloth-py}Python code for computing Euler spiral parameters.}
\end{figure*}

Complete code for this method is appealingly concise; the Python
version (shown in Figure \ref{solvecloth-py}) 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 code for evaluating the spiro
function, this approach is efficient indeed.