This page was generated from doc/rotation/catmull-rom-uniform.ipynb. Interactive online version: Binder badge.

Uniform Catmull–Rom-Like Quaternion Splines#

We have seen how to use De Casteljau’s algorithm with Slerp to create “cubic” Bézier-like quaternion curve segments. However, if we only have a sequence of rotations to be interpolated and no additional Bézier control quaternions are provided, it would be great if we could compute the missing control quaternions automatically from neighboring quaternions.

In the notebook about (uniform) Euclidean Catmull–Rom splines we have already seen how this can be done for splines in Euclidean space:

\begin{align*} \boldsymbol{\tilde{x}}_i^{(+)} &= \boldsymbol{x}_i + \frac{\boldsymbol{\dot{x}}_i}{3} \\ \boldsymbol{\tilde{x}}_i^{(-)} &= \boldsymbol{x}_i - \frac{\boldsymbol{\dot{x}}_i}{3} \end{align*}

Note that the velocity vectors \(\boldsymbol{\dot{x}}_i\) live in the same Euclidean space as the position vectors \(\boldsymbol{x}_i\). We can simply add a fraction of a velocity to a position and we get a new position in return.

Applying this to rotations is unfortunately not very straightforward. When unit quaternions are moving along the the unit hypersphere, their velocity vectors are tangential to that hypersphere, which means that the velocity vectors are generally not unit quaternions themselves. Furthermore, adding a (non-zero length) tangent vector to a unit quaternion never leads to a unit quaternion as a result.

Instead of using tangent vectors, we can introduce a (yet unknown) relative quaternion (in the global frame of reference) \(q_{i,\text{offset}}\):

\begin{align*} \tilde{q}_{i}^{(+)} &= {q_{i,\text{offset}}}^{\frac{1}{3}} \; q_i\\ \tilde{q}_{i}^{(-)} &= {q_{i,\text{offset}}}^{-\frac{1}{3}} \; q_i \end{align*}

When trying to obtain \(q_{i,\text{offset}}\), the problem is that there are many equivalent ways to write the equation for tangent vectors in Euclidean space …

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_{i-1}}{2} = \frac{(\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) + (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{2} = \frac{\boldsymbol{x}_i - \boldsymbol{x}_{i-1}}{2} + \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_i}{2} \end{equation*}

… but “translating” them to quaternions will lead to different results!

For the following experiments, let’s define three quaternions using the angles2quat() function from helper.py:

[1]:
from helper import angles2quat
[2]:
q3 = angles2quat(0, 0, 0)
q4 = angles2quat(0, 45, -10)
q5 = angles2quat(90, 0, -90)

Relative Rotations#

As a first attempt, we can try to “translate” the equation …

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_{i-1}}{2} \end{equation*}

… to unit quaternions like this:

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left(q_{i+1} {q_{i-1}}^{-1}\right)^{\frac{1}{2}} \end{equation*}

[3]:
offset_a = q3.rotation_to(q5)**(1/2)

We’ll see later whether that’s reasonable or not.

For the next few examples, we define the relative rotations associated with the the incoming and the outgoing chord:

\begin{align*} q_\text{in} &= q_i {q_{i-1}}^{-1}\\ q_\text{out} &= q_{i+1} {q_i}^{-1} \end{align*}

[4]:
q_in = q3.rotation_to(q4)
q_out = q4.rotation_to(q5)

The next equation …

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{(\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) + (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{2} \end{equation*}

… can be “translated” to unit quaternions like this:

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left(q_\text{out} q_\text{in} \right)^{\frac{1}{2}} \end{equation*}

[5]:
offset_b = (q_out * q_in)**(1/2)

We can see that this is actually equivalent to the previous one:

[6]:
max(map(abs, (offset_b - offset_a).xyzw))
[6]:
1.1102230246251565e-16

In the Euclidean case, the order doesn’t matter, but in the quaternion case …

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left(q_\text{in} q_\text{out} \right)^{\frac{1}{2}} \end{equation*}

[7]:
offset_c = (q_in * q_out)**(1/2)

… there is a (quite large!) difference:

[8]:
max(map(abs, (offset_b - offset_c).xyzw))
[8]:
0.2563304531880035

Based on the equation …

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{\boldsymbol{x}_i - \boldsymbol{x}_{i-1}}{2} + \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_i}{2} \end{equation*}

… we can try another pair of equations …

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left({q_\text{out}}^{\frac{1}{2}} {q_\text{in}}^{\frac{1}{2}} \right) \end{equation*}

[9]:
offset_d = (q_out**(1/2) * q_in**(1/2))

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left({q_\text{in}}^{\frac{1}{2}} {q_\text{out}}^{\frac{1}{2}} \right) \end{equation*}

[10]:
offset_e = (q_in**(1/6) * q_out**(1/6))

… but they are also non-symmetric:

[11]:
max(map(abs, (offset_e - offset_d).xyzw))
[11]:
0.20225984693486293

Let’s try a slightly more involved variant, where the order of \(q_\text{in}\) and \(q_\text{out}\) can actually be reversed:

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left( {q_\text{out}} {q_\text{in}}^{-1} \right)^{\frac{1}{2}} {q_\text{in}} = \left( {q_\text{in}} {q_\text{out}}^{-1} \right)^{\frac{1}{2}} {q_\text{out}} \end{equation*}

[12]:
offset_f = (q_out * q_in**-1)**(1/2) * q_in
[13]:
offset_g = (q_in * q_out**-1)**(1/2) * q_out
[14]:
max(map(abs, (offset_g - offset_f).xyzw))
[14]:
1.1102230246251565e-16

It is nice to have symmetric behavior, but the curvature of the unit hypersphere still causes an error. We can check that by scaling down the components before the calculation (leading to a smaller curvature) and scaling up the result:

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \left(\left( {q_\text{out}}^{\frac{1}{10}} {q_\text{in}}^{-\frac{1}{10}} \right)^{\frac{1}{2}} {q_\text{in}}^{\frac{1}{10}} \right)^{10} = \left(\left( {q_\text{in}}^{\frac{1}{10}} {q_\text{out}}^{-\frac{1}{10}} \right)^{\frac{1}{2}} {q_\text{out}}^{\frac{1}{10}} \right)^{10} \end{equation*}

[15]:
offset_h = ((q_out**(1/10) * q_in**(-1/10))**(1/2) * q_in**(1/10))**10
[16]:
offset_i = ((q_in**(1/10) * q_out**(-1/10))**(1/2) * q_out**(1/10))**10
[17]:
max(map(abs, (offset_h - offset_i).xyzw))
[17]:
2.1094237467877974e-15
[18]:
offset_j = ((q_out**(1/100) * q_in**(-1/100))**(1/2) * q_in**(1/100))**100
[19]:
offset_k = ((q_in**(1/100) * q_out**(-1/100))**(1/2) * q_out**(1/100))**100
[20]:
max(map(abs, (offset_j - offset_k).xyzw))
[20]:
1.4277468096679513e-13

If we choose a larger scaling factor, the the error caused by curvature becomes smaller (as we will see in the next section). However, the numerical error gets bigger. We cannot scale down the components arbitrarily, but there is a different mathematical tool that we can use, which boils down to the same thing, as we’ll see in the next section.

Tangent Space#

The logarithmic map operation transforms a unit quaternion into a vector that’s a member of the tangent space at the identity quaternion (a.k.a. \(\boldsymbol{1}\)). In this tangent space – which is a flat, three-dimensional Euclidean space – we can add and scale components without worrying about curvature. Using the exponential map operation, the result can be projected back onto the unit hypersphere. This way, we can take the equation for the tangent vector in Euclidean space …

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{(\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) + (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{2} \end{equation*}

… and “translate” it into unit quaternions …

\begin{equation*} q_{i,\text{offset}} \overset{?}{=} \exp\left( \frac{\ln(q_\text{in}) + \ln(q_\text{out})}{2} \right) \end{equation*}

[22]:
offset_l = UnitQuaternion.exp_map((q_in.log_map() + q_out.log_map()) / 2)

This approach is implemented in the splines.quaternion.CatmullRom class.

Let’s compare this to the variants from the previous section:

[23]:
max(map(abs, (offset_l - offset_f).xyzw))
[23]:
0.01742323752655639
[24]:
max(map(abs, (offset_l - offset_h).xyzw))
[24]:
0.000167758442754129
[25]:
max(map(abs, (offset_l - offset_j).xyzw))
[25]:
1.6769343111344703e-06

Increasing the scaling factor from the previous section will get us closer and closer, but only until the numerical errors eventually take over.

Example#

After all those more or less successful experiments, let’s show an example with actual rotations.

[26]:
def offset(q_1, q0, q1):
    q_in = q0 * q_1.inverse()
    q_out = q1 * q0.inverse()
    return UnitQuaternion.exp_map((q_in.log_map() + q_out.log_map()) / 2)

We’ll use the DeCasteljau class to create a Bézier-like curve from the given control points, using canonicalized() to avoid angles greater than 180 degrees.

Also, some helper functions from helper.py will come in handy.

[28]:
from helper import animate_rotations, display_animation

We don’t want to worry about end conditions here, so let’s create a closed curve.

[29]:
def create_closed_curve(rotations):
    rotations = list(canonicalized(rotations + rotations[:2]))
    control_points = []
    for q_1, q0, q1 in zip(rotations, rotations[1:], rotations[2:]):
        q_offset = offset(q_1, q0, q1)
        control_points.extend([
            q_offset**(-1/3) * q0,
            q0,
            q0,
            q_offset**(1/3) * q0])
    control_points = control_points[-2:] + control_points[:-2]
    segments = list(zip(*[iter(control_points)] * 4))
    return DeCasteljau(segments)
[30]:
rotations = [
    angles2quat(0, 0, 180),
    angles2quat(0, 45, 90),
    angles2quat(90, 45, 0),
    angles2quat(90, 90, -90),
    angles2quat(180, 0, -180),
    angles2quat(-90, -45, 180),
]
[31]:
s = create_closed_curve(rotations)
[32]:
import numpy as np
[33]:
times = np.linspace(0, len(rotations), 200, endpoint=False)
[34]:
ani = animate_rotations(s.evaluate(times))
[35]:
display_animation(ani, default_mode='loop')

Shoemake’s Approach#

In section 4.2, Shoemake [Sho85] provides two function definitions:

\begin{align*} \operatorname{Double}(p, q) &= 2 (p \cdot q) q - p\\ \operatorname{Bisect}(p, q) &= \frac{p + q}{\|p + q\|} \end{align*}

[36]:
def double(p, q):
    return 2 * p.dot(q) * q - p
[37]:
def bisect(p, q):
    return (p + q).normalized()

Given three successive key quaternions \(q_{n-1}\), \(q_n\) and \(q_{n+1}\), these functions are used to compute control quaternions \(b_n\) (controlling the incoming tangent of \(q_n\)) and \(a_n\) (controlling the outgoing tangent of \(q_n\)):

\begin{align*} a_n &= \operatorname{Bisect}(\operatorname{Double}(q_{n-1}, q_n), q_{n+1})\\ b_n &= \operatorname{Double}(a_n, q_n) \end{align*}

It is unclear where these equations come from, we only get a little hint:

For the numerically knowledgeable, this construction approximates the derivative at points of a sampled function by averaging the central differences of the sample sequence.

Shoemake [Sho85], footnote on page 249

[38]:
def shoemake_control_quaternions(q_1, q0, q1):
    """Shoemake's control quaternions.

    Given three key quaternions, return the control quaternions
    preceding and following the middle one.

    Actually, the great arc distance of the returned quaternions to q0
    still has to be reduced to 1/3 of the distance
    to get the proper control quaternions (see the note below).

    """
    a = bisect(double(q_1, q0), q1)
    b = double(a, q0).normalized()
    return b, a

Normalization of \(b_n\) is not explicitly mentioned in the paper, but even though the results have a length very close to 1.0, we still have to call normalized() to turn the Quaternion result into a UnitQuaternion.

[39]:
b, a = shoemake_control_quaternions(q3, q4, q5)

The results are close (but by far not identical) to the tangent space approach from above:

[40]:
max(map(abs, (a - offset_l * q4).xyzw))
[40]:
0.013831724198409168
[41]:
max(map(abs, (b - offset_l.inverse() * q4).xyzw))
[41]:
0.018852903209093046

Note

Shoemake’s result has to be scaled by \(\frac{1}{3}\), just as we did with \(q_{i,\text{offset}}\) above:

A simple check proves the curve touches \(q_n\) and \(q_{n+1}\) at its ends. A rather challenging differentiation shows it is tangent there to the segments determined by \(a_n\) and \(b_{n+1}\). However, as with Bézier’s original curve, the magnitude of the tangent is three times that of the segment itself. That is, we are spinning three times faster than spherical interpolation along the arc. Fortunately we can correct the speed by merely truncating the end segments to one third their original length, so that \(a_n\) is closer to \(q_n\) and \(b_{n+1}\) closer to \(q_{n+1}\).

Shoemake [Sho85], section 4.4: “Tangents revisited”