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

# Uniform Cubic Hermite Splines

TODO: Image for the 1D case (points and slopes)?

TODO: Image for the 2D case (points and tangent vectors)? Probably combine two 1D examples?

TODO: Hermite's two-point interpolation formula?


```python
%matplotlib inline
from IPython.display import display
```


```python
import sympy as sp
sp.init_printing(order='rev-lex')
```


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


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


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


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


```python
p = b_monomial.dot(coefficients)
p
```

We want to provide the value of the polynomial at time $t_0 = 0$ and $t_1 = 1$ (start and end point).
We also want to provide the first derivative (a.k.a. velocity, a.k.a. tangent vector) at the start and at the end.

\begin{align}
\boldsymbol{x}_0 &= \boldsymbol{p}(0)\\
\boldsymbol{x}_1 &= \boldsymbol{p}(1)\\
\boldsymbol{\dot{x}}_0 &= \boldsymbol{p}'(0)\\
\boldsymbol{\dot{x}}_1 &= \boldsymbol{p}'(1)
\end{align}

The curve segment defined by these 4 values is sometimes called *Ferguson cubic*.

Velocity = Tangent Vector = First Derivative:


```python
velocity = p.diff(t)
velocity
```


```python
x0 = NamedExpression('xbm0', p.subs(t, 0))
x1 = NamedExpression('xbm1', p.subs(t, 1))
xd0 = NamedExpression('xdotbm0', velocity.subs(t, 0))
xd1 = NamedExpression('xdotbm1', velocity.subs(t, 1))
```


```python
display(x0, x1, xd0, xd1)
```

Given an input vector of control values ...


```python
control_values_H = NamedMatrix(sp.Matrix([x0.name, x1.name, xd0.name, xd1.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
```

And the other way round:


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


```python
substitutions = x0, x1, xd0, xd1
```


```python
control_values_H.subs(substitutions)
```


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


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

This transforms the coefficients into our control values, but we need it the other way round:


```python
M_H
```


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

Those are called Hermite basis functions:


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


```python
sp.plot(*b_H.expr, (t, 0, 1));
```

TODO: describe properties of Hermite basis functions

## Plotting


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


```python
def hermite_curve(control_values):
    coeffs = sp.lambdify([], M_H.expr)() @ control_values
    
    def evaluate(t):
        t = np.expand_dims(t, -1)
        return t**[3, 2, 1, 0] @ coeffs
        
    return evaluate
```


```python
t_values = np.linspace(0, 1, 15)
```


```python
def plot_hermite(control_values, t_values, ax=None):
    if ax is None:
        ax = plt.gca()
    control_values = np.asarray(control_values)
    c = hermite_curve(control_values)
    xyuv = np.column_stack(np.vsplit(control_values, 2))
    # Move beginning of second arrow:
    xyuv[1, :2] -= xyuv[1, 2:]
    ax.quiver(*xyuv.T, angles='xy', scale_units='xy', scale=1, color='lightgrey')
    ax.plot(*c(t_values).T, '.')
    ax.axis('equal')
```


```python
plot_hermite([[0, 0], [5, 2], [1, 2], [1, -2]], t_values)
plot_hermite([[5, 2], [6, 1], [1, -2], [5, 1]], t_values)
```

TODO: same length, different angles?


```python
plot_hermite([[0, 0], [1, 0], [1, 1], [1, -1]], t_values)
plot_hermite([[0, 0], [1, 0], [0, 1], [0, -1]], t_values)
```

## Relation to Bézier Splines

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


```python
coefficients_H
```

What about using four positions instead?

Let's use the point $\boldsymbol{\hat{x}}_0$ as "handle" connected to $\boldsymbol{x}_0$.
Same for $\boldsymbol{\hat{x}}_1$ and $\boldsymbol{x}_1$.

And since the tangents looked unwieldily long in the examples before, let's make the handles only a third of the length of the tangents:


```python
x_hat0 = NamedExpression('xhatbm0', x0.name + xd0.name / 3)
x_hat1 = NamedExpression('xhatbm1', x1.name - xd1.name / 3)
```


```python
control_values_B = NamedMatrix(
    sp.Matrix([x0.name, x_hat0.name, x_hat1.name, x1.name]),
    sp.Matrix([x0.name, x_hat0.expr, x_hat1.expr, x1.name]))
control_values_B
```


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


```python
sp.Eq(control_values_B.name, M_HtoB.name * control_values_H.name)
```


```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)
```


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


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


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


```python
M_B = M_B.subs([M_H, M_BtoH])
M_B
```


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


```python
sp.plot(*b_B.expr, (t, 0, 1));
```

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

3bef90e2a5d6d5356aa6632546d839d26dc93d17
ipynb
Jupyter Notebook
doc/hermite-uniform.ipynb
mgeier/splines
f54b09479d98bf13f00a183fd9d664b5783e3864
MIT
doc/hermite-uniform.ipynb
mgeier/splines
f54b09479d98bf13f00a183fd9d664b5783e3864
MIT
doc/hermite-uniform.ipynb
mgeier/splines
f54b09479d98bf13f00a183fd9d664b5783e3864
MIT
Qwen/Qwen-72B
1. YES
2. YES