This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.

[back to overview](hermite.ipynb)

# Uniform Cubic Hermite Splines

We derive the basis matrix as well as the basis polynomials
for cubic (= degree 3) Hermite splines.
The derivation for other degrees is left as an exercise for the reader.

In this notebook,
we consider *uniform* spline segments,
i.e. the parameter in each segment varies from $0$ to $1$.
The derivation for *non-uniform* cubic Hermite splines
can be found in [a separate notebook](hermite-non-uniform.ipynb).


```python
import sympy as sp
sp.init_printing(order='grevlex')
```

We load a few tools from [utility.py](utility.py):


```python
from utility import NamedExpression, NamedMatrix
```


```python
t = sp.symbols('t')
```

We are considering a single cubic polynomial segment of a Hermite spline
(which is sometimes called a *Ferguson cubic*).

To simplify the indices in the following derivation,
we are looking at the fifth polynomial segment $\boldsymbol{p}_4(t)$
from $\boldsymbol{x}_4$ to $\boldsymbol{x}_5$, where $0 \le t \le 1$.
The results will be easily generalizable to an arbitrary
polynomial segment $\boldsymbol{p}_i(t)$
from $\boldsymbol{x}_i$ to $\boldsymbol{x}_{i+1}$.

The polynomial has 4 coefficients, $\boldsymbol{a_4}$ to $\boldsymbol{d_4}$.


```python
coefficients = sp.Matrix(sp.symbols('a:dbm4')[::-1])
coefficients
```

Combined with the *monomial basis* ...


```python
b_monomial = sp.Matrix([t**3, t**2, t, 1]).T
b_monomial
```

... the coefficients form an expression
for our polynomial segment $\boldsymbol{p}_4(t)$:


```python
p4 = NamedExpression('pbm4', b_monomial.dot(coefficients))
p4
```

For more information about polynomials,
see [Polynomial Parametric Curves](polynomials.ipynb).

Let's also calculate the first derivative
(a.k.a. velocity, a.k.a. tangent vector),
while we are at it:


```python
pd4 = p4.diff(t)
pd4
```

To generate a Hermite spline segment,
we have to provide the value of the polynomial
at the start and end point of the segment
(at times $t = 0$ and $t = 1$, respectively).
We also have to provide the first derivative at those same points.

\begin{align*}
\boldsymbol{x}_4 &= \left.\boldsymbol{p}_4\right\rvert_{t=0}\\
\boldsymbol{x}_5 &= \left.\boldsymbol{p}_4\right\rvert_{t=1}\\
\boldsymbol{\dot{x}}_4 &= \left.\frac{d}{dt}\boldsymbol{p}_4\right\rvert_{t=0}\\
\boldsymbol{\dot{x}}_5 &= \left.\frac{d}{dt}\boldsymbol{p}_4\right\rvert_{t=1}
\end{align*}

We call those 4 values the *control values* of the segment.

Evaluating the polynomial and its derivative
at times $0$ and $1$ leads to 4 expressions for our 4 control values:


```python
x4 = p4.evaluated_at(t, 0).with_name('xbm4')
x5 = p4.evaluated_at(t, 1).with_name('xbm5')
xd4 = pd4.evaluated_at(t, 0).with_name('xdotbm4')
xd5 = pd4.evaluated_at(t, 1).with_name('xdotbm5')
```


```python
display(x4, x5, xd4, xd5)
```

## Basis Matrix

Given an input vector of control values ...


```python
control_values_H = NamedMatrix(sp.Matrix([x4.name,
                                          x5.name,
                                          xd4.name,
                                          xd5.name]))
control_values_H.name
```

... we want to find a way to transform those into the coefficients of our cubic polynomial.


```python
M_H = NamedMatrix(r'{M_\text{H}}', 4, 4)
```


```python
coefficients_H = NamedMatrix(coefficients, M_H.name * control_values_H.name)
coefficients_H
```

This way, we can express our previously unknown coefficients
in terms of the given control values.

However, in order to make it easy to determine
the coefficients of the *basis matrix* $M_H$,
we need the equation the other way around
(by left-multiplying by the inverse):


```python
control_values_H.expr = M_H.name.I * coefficients
control_values_H
```

We can now insert the expressions for the control values
that we obtained above ...


```python
substitutions = x4, x5, xd4, xd5
```


```python
control_values_H.subs_symbols(*substitutions)
```

... and from this equation we can directly read off
the matrix coefficients of ${M_H}^{-1}$:


```python
M_H.I = sp.Matrix(
    [[expr.coeff(cv) for cv in coefficients]
     for expr in control_values_H.subs_symbols(*substitutions).name])
M_H.I
```

The same thing for copy & paste purposes:


```python
print(_.expr)
```

This transforms the coefficients of the polynomial into our control values,
but we need it the other way round,
which we can simply get by inverting the matrix:


```python
M_H
```

Again, for copy & paste:


```python
print(_.expr)
```

## Basis Polynomials

Multiplying the monomial basis with this matrix
yields the *Hermite basis polynomials*:


```python
b_H = NamedMatrix(r'{b_\text{H}}', b_monomial * M_H.expr)
b_H.factor().simplify().T
```

Let's plot the basis polynomials
with some help from [helper.py](helper.py):


```python
from helper import plot_basis
```


```python
plot_basis(*b_H.expr, labels=sp.symbols('xbm_i xbm_i+1 xdotbm_i xdotbm_i+1'))
```

Note that the basis function associated with $\boldsymbol{x}_i$
has the value $1$ at the beginning,
while all others are $0$ at that point.
For this reason,
the linear combination of all basis functions at $t=0$
simply adds up to the value $\boldsymbol{x}_i$
(which is exactly what we wanted to happen!).

Similarly,
the basis function associated with $\boldsymbol{\dot{x}}_i$
has a first derivative of $+1$ at the beginning,
while all others have a first derivative of $0$.
Therefore,
the linear combination of all basis functions at $t=0$
turns out to have a first derivative of $\boldsymbol{\dot{x}}_i$
(what a coincidence!).

While $t$ progresses towards $1$,
both functions must relinquish their influence
to the other two basis functions.

At the end (when $t=1$),
the basis function associated with $\boldsymbol{x}_{i+1}$
is the only one that has a non-zero value.
More concretely, it has the value $1$.
Finally,
the basis function associated with $\boldsymbol{\dot{x}}_{i+1}$
is the only one with a non-zero first derivative.
In fact, it has a first derivative of exactly $+1$
(the function values leading up to that have to be negative
because the final function value has to be $0$).

This can be summarized by:


```python
sp.Matrix([[
    b.subs(t, 0),
    b.subs(t, 1),
    b.diff(t).subs(t, 0),
    b.diff(t).subs(t, 1),
] for b in b_H.expr])
```

## Example Plot

To quickly check whether the matrix $M_H$ does what we expect,
let's plot an example segment:


```python
import numpy as np
```

If we use the same API as for the other splines,
we can reuse the helper functions for plotting
from [helper.py](helper.py).


```python
from helper import plot_spline_2d, plot_tangents_2d
```


```python
class UniformHermiteSegment:

    grid = 0, 1

    def __init__(self, control_values):
        self.coeffs = sp.lambdify([], M_H.expr)() @ control_values

    def evaluate(self, t):
        t = np.expand_dims(t, -1)
        return t**[3, 2, 1, 0] @ self.coeffs
```

<div class="alert alert-info">

Note

The `@` operator is used here to do
[NumPy's matrix multiplication](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html).

</div>


```python
vertices = [0, 0], [5, 1]
tangents = [2, 3], [0, -2]
```


```python
s = UniformHermiteSegment([*vertices, *tangents])
```


```python
plot_spline_2d(s, chords=False)
plot_tangents_2d(tangents, vertices)
```

## Relation to Bézier Splines

Above, we were using two positions (start and end)
and two tangent vectors (at those same two positions) as control values:


```python
control_values_H.name
```

What about using four positions (and no tangent vectors) instead?

Let's use the point $\boldsymbol{\tilde{x}}_4$ as a "handle"
(connected to $\boldsymbol{x}_4$) that controls the tangent vector.
Same for $\boldsymbol{\tilde{x}}_5$ (connected to $\boldsymbol{x}_5$).

And since the tangents looked unwieldily long in the plot above
(compared to the effect they have on the shape of the curve),
let's put the handles only at a third of the length of the tangents,
shall we?

\begin{align*}
\tilde{\boldsymbol{x}}_4
&=
\boldsymbol{x}_4 + \frac{\dot{\boldsymbol{x}}_4}{3}
\\
\tilde{\boldsymbol{x}}_5
&=
\boldsymbol{x}_5 - \frac{\dot{\boldsymbol{x}}_5}{3}
\end{align*}


```python
control_values_B = NamedMatrix(sp.Matrix([
    x4.name,
    sp.Symbol('xtildebm4'),
    sp.Symbol('xtildebm5'),
    x5.name,
]), sp.Matrix([
    x4.name,
    x4.name + xd4.name / 3,
    x5.name - xd5.name / 3,
    x5.name,
]))
control_values_B
```

Now let's try to come up with a matrix
that transforms our good old Hermite control values
into our new control points.


```python
M_HtoB = NamedMatrix(r'{M_\text{H$\to$B}}', 4, 4)
```


```python
NamedMatrix(control_values_B.name, M_HtoB.name * control_values_H.name)
```

We can immediately read the matrix coefficients
off the previous expression.


```python
M_HtoB.expr = sp.Matrix([[expr.coeff(cv) for cv in control_values_H.name]
                         for expr in control_values_B.expr])
M_HtoB.pull_out(sp.S.One / 3)
```


```python
print(_.expr)
```

The inverse of this matrix transforms our new control points
into Hermite control values:


```python
M_BtoH = NamedMatrix(r'{M_\text{B$\to$H}}', M_HtoB.I.expr)
M_BtoH
```


```python
print(_.expr)
```

When we combine $M_H$ with this new matrix,
we get a matrix which leads us to a new set of basis polynomials
associated with the 4 control points.


```python
M_B = NamedMatrix(r'{M_\text{B}}', M_H.name * M_BtoH.name)
M_B
```


```python
M_B = M_B.subs_symbols(M_H, M_BtoH).doit()
M_B
```


```python
b_B = NamedMatrix(r'{b_\text{B}}', b_monomial * M_B.expr)
b_B.T
```


```python
plot_basis(
    *b_B.expr,
    labels=sp.symbols('xbm_i xtildebm_i xtildebm_i+1 xbm_i+1'))
```

Those happen to be the cubic *Bernstein* polynomials and
it turns out that we just invented *Bézier* curves!
See [the section about Bézier splines](bezier.ipynb)
for more information about them.

We chose the additional control points to be located
at $\frac{1}{3}$ of the tangent vector.
Let's quickly visualize this
using the example from above and $M_\text{H$\to$B}$:


```python
points = sp.lambdify([], M_HtoB.expr)() @ [*vertices, *tangents]
```


```python
import matplotlib.pyplot as plt
```


```python
plot_spline_2d(s, chords=False)
plot_tangents_2d(tangents, vertices)
plt.scatter(*points.T, marker='X', color='black')
plt.annotate(r'$\quad\tilde{\bf{x}}_0$', points[1])
plt.annotate(r'$\tilde{\bf{x}}_1\quad$', points[2], ha='right');
```

aee899ea3f1bc7ab2f9318649dca3f4ba044e100
ipynb
Jupyter Notebook
doc/euclidean/hermite-uniform.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2021-06-21T11:41:02.000Z
2022-02-15T19:18:47.000Z
doc/euclidean/hermite-uniform.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2018-07-09T17:47:29.000Z
2022-01-28T19:55:55.000Z
doc/euclidean/hermite-uniform.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2020-05-15T18:34:42.000Z
2021-12-28T06:29:54.000Z
Qwen/Qwen-72B
1. YES
2. YES