This page was generated from doc/euclidean/hermite-uniform.ipynb. Interactive online version: .

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.

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

We load a few tools from utility.py:

from utility import NamedExpression, NamedMatrix

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, let’s look at only one specific polynomial segment, let’s say the fifth one. It goes from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\) and it is referred to as \(\boldsymbol{p}_4(t)\), 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}\), where \(0 \le t \le 1\).

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

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

$\displaystyle \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right]$

Combined with the monomial basis …

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

$\displaystyle \left[\begin{matrix}t^{3} & t^{2} & t & 1\end{matrix}\right]$

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

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

$\displaystyle \boldsymbol{p}_{4} = \boldsymbol{d}_{4} t^{3} + \boldsymbol{c}_{4} t^{2} + \boldsymbol{b}_{4} t + \boldsymbol{a}_{4}$

For more information about polynomials, see Polynomial Parametric Curves.

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

pd4 = p4.diff(t)
pd4

$\displaystyle \frac{d}{d t} \boldsymbol{p}_{4} = 3 \boldsymbol{d}_{4} t^{2} + 2 \boldsymbol{c}_{4} t + \boldsymbol{b}_{4}$

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:

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

display(x4, x5, xd4, xd5)

$\displaystyle \boldsymbol{x}_{4} = \boldsymbol{a}_{4}$

$\displaystyle \boldsymbol{x}_{5} = \boldsymbol{a}_{4} + \boldsymbol{b}_{4} + \boldsymbol{c}_{4} + \boldsymbol{d}_{4}$

$\displaystyle \boldsymbol{\dot{x}}_{4} = \boldsymbol{b}_{4}$

$\displaystyle \boldsymbol{\dot{x}}_{5} = \boldsymbol{b}_{4} + 2 \boldsymbol{c}_{4} + 3 \boldsymbol{d}_{4}$

Basis Matrix

Given an input vector of control values …

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

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

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

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

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

$\displaystyle \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right] = {M_\text{H}} \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

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

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

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right] = {M_\text{H}}^{-1} \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right]$

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

substitutions = x4, x5, xd4, xd5

control_values_H.subs_symbols(*substitutions)

$\displaystyle \left[\begin{matrix}\boldsymbol{a}_{4}\\\boldsymbol{a}_{4} + \boldsymbol{b}_{4} + \boldsymbol{c}_{4} + \boldsymbol{d}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{b}_{4} + 2 \boldsymbol{c}_{4} + 3 \boldsymbol{d}_{4}\end{matrix}\right] = {M_\text{H}}^{-1} \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right]$

… and from this equation we can directly read off the matrix coefficients of \({M_H}^{-1}\):

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

$\displaystyle {M_\text{H}}^{-1} = \left[\begin{matrix}0 & 0 & 0 & 1\\1 & 1 & 1 & 1\\0 & 0 & 1 & 0\\3 & 2 & 1 & 0\end{matrix}\right]$

The same thing for copy & paste purposes:

print(_.expr)

Matrix([[0, 0, 0, 1], [1, 1, 1, 1], [0, 0, 1, 0], [3, 2, 1, 0]])

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:

M_H

$\displaystyle {M_\text{H}} = \left[\begin{matrix}2 & -2 & 1 & 1\\-3 & 3 & -2 & -1\\0 & 0 & 1 & 0\\1 & 0 & 0 & 0\end{matrix}\right]$

Again, for copy & paste:

print(_.expr)

Matrix([[2, -2, 1, 1], [-3, 3, -2, -1], [0, 0, 1, 0], [1, 0, 0, 0]])

Now we have a new way to write the polynomial \(\boldsymbol{p}_4(t)\), given our four control values. We take those control values, left-multiply them by the Hermite basis matrix \(M_\text{H}\) (which gives us a column vector of coefficients), which we can then left-multiply by the monomial basis:

sp.MatMul(b_monomial, M_H.expr, control_values_H.name)

$\displaystyle \left[\begin{matrix}t^{3} & t^{2} & t & 1\end{matrix}\right] \left[\begin{matrix}2 & -2 & 1 & 1\\-3 & 3 & -2 & -1\\0 & 0 & 1 & 0\\1 & 0 & 0 & 0\end{matrix}\right] \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

Basis Polynomials

However, instead of calculating from right to left, we can also start at the left and multiply the monomial basis with the Hermite basis matrix \(M_\text{H}\), which yields (a row vector containing) the Hermite basis polynomials:

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

$\displaystyle {b_\text{H}}^{T} = \left[\begin{matrix}\left(t - 1\right)^{2} \cdot \left(2 t + 1\right)\\- t^{2} \cdot \left(2 t - 3\right)\\t \left(t - 1\right)^{2}\\t^{2} \left(t - 1\right)\end{matrix}\right]$

The multiplication of this row vector with the column vector of control values again produces the polynomial \(\boldsymbol{p}_4(t)\).

Let’s plot the basis polynomials with some help from helper.py:

from helper import plot_basis

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 specifically, 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:

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

$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$

Example Plot

To quickly check whether the matrix \(M_H\) does what we expect, let’s plot an example segment.

import matplotlib.pyplot as plt
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.

from helper import plot_spline_2d, plot_tangents_2d, plot_vertices_2d

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

Note

The @ operator is used here to do NumPy’s matrix multiplication.

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

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

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:

control_values_H.name

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

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

Let’s use the point \(\boldsymbol{\tilde{x}}_4\) as a “drag point” (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 drag points 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*}

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

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{\tilde{x}}_{4}\\\boldsymbol{\tilde{x}}_{5}\\\boldsymbol{x}_{5}\end{matrix}\right] = \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{4} + \frac{\boldsymbol{\dot{x}}_{4}}{3}\\\boldsymbol{x}_{5} - \frac{\boldsymbol{\dot{x}}_{5}}{3}\\\boldsymbol{x}_{5}\end{matrix}\right]$

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

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

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

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{\tilde{x}}_{4}\\\boldsymbol{\tilde{x}}_{5}\\\boldsymbol{x}_{5}\end{matrix}\right] = {M_\text{H$\to$B}} \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

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

M_HtoB.expr = sp.Matrix([
    [expr.coeff(cv) for cv in control_values_H.name]
    for expr in control_values_B.expr])
M_HtoB

$\displaystyle {M_\text{H$\to$B}} = \left[\begin{matrix}1 & 0 & 0 & 0\\1 & 0 & \frac{1}{3} & 0\\0 & 1 & 0 & - \frac{1}{3}\\0 & 1 & 0 & 0\end{matrix}\right]$

print(_.expr)

Matrix([[1, 0, 0, 0], [1, 0, 1/3, 0], [0, 1, 0, -1/3], [0, 1, 0, 0]])

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

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

$\displaystyle {M_\text{B$\to$H}} = \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 1\\-3 & 3 & 0 & 0\\0 & 0 & -3 & 3\end{matrix}\right]$

print(_.expr)

Matrix([[1, 0, 0, 0], [0, 0, 0, 1], [-3, 3, 0, 0], [0, 0, -3, 3]])

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.

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

$\displaystyle {M_\text{B}} = {M_\text{H}} {M_\text{B$\to$H}}$

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

$\displaystyle {M_\text{B}} = \left[\begin{matrix}-1 & 3 & -3 & 1\\3 & -6 & 3 & 0\\-3 & 3 & 0 & 0\\1 & 0 & 0 & 0\end{matrix}\right]$

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

$\displaystyle {b_\text{B}}^{T} = \left[\begin{matrix}- t^{3} + 3 t^{2} - 3 t + 1\\3 t^{3} - 6 t^{2} + 3 t\\- 3 t^{3} + 3 t^{2}\\t^{3}\end{matrix}\right]$

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 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}\):

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

plot_spline_2d(s, chords=False)
plot_tangents_2d(tangents, vertices)
plot_vertices_2d(points, chords=False, markeredgecolor='purple')
plt.annotate(r'$\quad\tilde{\bf{x}}_0$', points[1])
plt.annotate(r'$\quad\tilde{\bf{x}}_1$', points[2]);