.
Cumulative Form#
The basic idea, as proposed by Kim et al. [KKS95] is the following:
Instead of representing a curve as a sum of basis functions weighted by its control point’s position vectors \(p_i\) – as it’s for example done with Bézier splines – they suggest to use the relative difference vectors \(\Delta p_i\) between successive control points.
These relative difference vectors can then be “translated” to local rotations (replacing additions with multiplications), leading to a form of rotation splines.
Piecewise Slerp#
As an example, they define a piecewise linear curve
\begin{equation*} p(t) = p_0 + \sum_{i=1}^n \alpha_i(t) \Delta p_i, \end{equation*}
where
\begin{align*} \Delta p_i &= p_i - p_{i - 1}\\ \alpha_i(t) &= \begin{cases} 0 & t < i - 1\\ t - i + 1 & i - 1 \leq t < i\\ 1 & t \geq i. \end{cases} \end{align*}
[1]:
def alpha(i, t):
if t < i - 1:
return 0
elif t >= i:
return 1
else:
return t - i + 1
Note
There is an off-by-one error in the paper’s definition of \(\alpha_i(t)\):
\begin{equation*} \alpha_i(t) = \begin{cases} 0 & t < i\\ t - i & i \leq t < i + 1\\ 1 & t \geq i + 1. \end{cases} \end{equation*}
This assumes that \(i\) starts with \(0\), but it actually starts with \(1\).
This “cumulative form” can be “translated” to a rotation spline by replacing addition with multiplication and the relative difference vectors by relative (i.e. local) rotations (represented by unit quaternions):
\begin{equation*} q(t) = q_0 \prod_{i = 1}^n \exp(\omega_i \alpha_i(t)), \end{equation*}
where
\begin{equation*} \omega_i = \log\left(q_{i - 1}^{-1} q_i\right). \end{equation*}
The paper uses above notation, but this could equivalently be written as
\begin{equation*} q(t) = q_0 \prod_{i = 1}^n \left(q_{i - 1}^{-1} q_i\right)^{\alpha_i(t)}. \end{equation*}
[2]:
import numpy as np
Let’s import a few helper functions from helper.py:
[3]:
from helper import angles2quat, animate_rotations, display_animation
[4]:
from splines.quaternion import UnitQuaternion
[5]:
# NB: math.prod() since Python 3.8
product = np.multiply.reduce
[6]:
[7]:
qs = [
angles2quat(0, 0, 0),
angles2quat(90, 0, 0),
angles2quat(90, 90, 0),
angles2quat(90, 90, 90),
]
[8]:
times = np.linspace(0, len(qs) - 1, 100)
[9]:
ani = animate_rotations([piecewise_slerp(qs, t) for t in times])
[10]:
display_animation(ani, default_mode='reflect')
Cumulative Bézier/Bernstein Curve#
After the piecewise Slerp, Kim, Kim and Shin (1995) show (in section 5.1) how to create a cumulative form inspired by Bézier splines, i.e. using Bernstein polynomials.
They start with the well-known equation for Bézier splines:
\begin{equation*} p(t) = \sum_{i=0}^n p_i \beta_{i,n}(t), \end{equation*}
where \(\beta_{i,n}(t)\) are Bernstein basis functions as shown in the notebook about Bézier splines.
They re-formulate this into a cumulative form:
\begin{equation*} p(t) = p_0 \tilde{\beta}_{0,n}(t) + \sum_{i=1}^n \Delta p_i \tilde{\beta}_{i,n}(t), \end{equation*}
where the cumulative Bernstein basis functions are given by
\begin{equation*} \tilde{\beta}_{i,n}(t) = \sum_{j=i}^n \beta_{j,n}(t). \end{equation*}
We can get the Bernstein basis polynomials via the function splines.Bernstein.basis() …
… and create a simple helper function to sum them up:
[12]:
from itertools import accumulate
[13]:
def cumulative_bases(degree, t):
return list(accumulate(Bernstein.basis(degree, t)[::-1]))[::-1]
Finally, they “translate” this into a rotation spline using quaternions, like before:
\begin{equation*} q(t) = q_0 \prod_{i=1}^n \exp\left(\omega_i \tilde{\beta}_{i,n}(t)\right), \end{equation*}
where
\begin{equation*} \omega_i = \log(q_{i-1}^{-1} q_i). \end{equation*}
Again, they use above notation in the paper, but this could equivalently be written as
\begin{equation*} q(t) = q_0 \prod_{i=1}^n \left(q_{i-1}^{-1} q_i\right)^{\tilde{\beta}_{i,n}(t)}. \end{equation*}
[14]:
def cumulative_bezier(qs, t):
degree = len(qs) - 1
bases = cumulative_bases(degree, t)
assert np.isclose(bases[0], 1)
return qs[0] * product([
(qs[i - 1].inverse() * qs[i])**bases[i]
for i in range(1, len(qs))
])
[15]:
times = np.linspace(0, 1, 100)
[16]:
rotations = [cumulative_bezier(qs, t) for t in times]
[17]:
ani = animate_rotations(rotations)
[18]:
display_animation(ani, default_mode='reflect')
Comparison with De Casteljau’s Algorithm#
The method described by Shoemake [Sho85] is shown in a separate notebook. An implementation is available in the class splines.quaternion.DeCasteljau:
[19]:
from splines.quaternion import DeCasteljau
[20]:
times = np.linspace(0, 1, 100)
[21]:
control_polygon = [
angles2quat(90, 0, 0),
angles2quat(0, -45, 90),
angles2quat(0, 0, 0),
angles2quat(180, 0, 180),
]
[22]:
cumulative_rotations = [
cumulative_bezier(control_polygon, t)
for t in times
]
[23]:
cumulative_rotations_reversed = [
cumulative_bezier(control_polygon[::-1], t)
for t in times
][::-1]
[24]:
casteljau_rotations = DeCasteljau([control_polygon]).evaluate(times)
[25]:
ani = animate_rotations({
'De Casteljau': casteljau_rotations,
'Cumulative': cumulative_rotations,
'Cumulative reversed': cumulative_rotations_reversed,
})
[26]:
display_animation(ani, default_mode='reflect')
Applying the same method on the reversed list of control points and then time-reversing the resulting sequence of rotations leads to an equal (except for rounding errors) sequence of rotations when using De Casteljau’s algorithm:
[27]:
casteljau_rotations_reversed = DeCasteljau([control_polygon[::-1]]).evaluate(times)[::-1]
[28]:
for one, two in zip(casteljau_rotations, casteljau_rotations_reversed):
assert np.isclose(one.scalar, two.scalar)
assert np.isclose(one.vector[0], two.vector[0])
assert np.isclose(one.vector[1], two.vector[1])
assert np.isclose(one.vector[2], two.vector[2])
However, doing the same thing with the “cumulative form” can lead to a significantly different sequence, as can be seen in the above animation.