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

Non-Uniform Cubic Hermite Splines#

We have already derived uniform cubic Hermite splines, where the parameter \(t\) ranges from \(0\) to \(1\).

When we want to use non-uniform cubic Hermite splines, and therefore arbitrary ranges from \(t_i\) to \(t_{i+1}\), we have (at least) two possibilities:

  • Do the same derivations as in the uniform case, except when we previously evaluated an expression at the parameter value \(t=0\), we now evaluate it at the value \(t=t_i\). Of course we do the same with \(t = 1 \to t = t_{i+1}\).

  • Re-scale the non-uniform parameter using \(t \to \frac{t - t_i}{t_{i+1} - t_i}\) (which makes the new parameter go from \(0\) to \(1\)) and then simply use the results from the uniform case.

The first approach leads to more complicated expressions in the basis matrix and the basis polynomials, but it has the advantage that the parameter value doesn’t have to be re-scaled each time when evaluating the spline for a given parameter (which might be slightly more efficient).

The second approach has the problem that it doesn’t actually work correctly, but we will see that we can make a slight adjustment to fix that problem (spoiler alert: we will have to multiply the tangent vectors by \(\Delta_i\)).

The class splines.CubicHermite is implemented using the second approach (because its parent class splines.Monomial also uses the re-scaling approach).

We show the second approach here, but the first approach can be carried out very similarly, with only very few changed steps. The appropriate changes are mentioned below.

[1]:
from pprint import pprint
import sympy as sp
sp.init_printing(order='grevlex')
[2]:
from utility import NamedExpression, NamedMatrix

To simplify the indices in the following derivation, we are again looking at the fifth polynomial segment \(\boldsymbol{p}_4(t)\) from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\), where \(t_4 \le t \le t_5\). 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 \(t_i \le t \le t_{i+1}\).

[3]:
t, t4, t5 = sp.symbols('t t4:6')
[4]:
coefficients = sp.Matrix(sp.symbols('a:dbm4')[::-1])
b_monomial = sp.Matrix([t**3, t**2, t, 1]).T
b_monomial.dot(coefficients)
[4]:
$\displaystyle \boldsymbol{d}_{4} t^{3} + \boldsymbol{c}_{4} t^{2} + \boldsymbol{b}_{4} t + \boldsymbol{a}_{4}$

We use the humble cubic polynomial (with monomial basis) to represent our curve segment \(\boldsymbol{p}_4(t)\), but we re-scale the parameter to map \(t_4 \to 0\) and \(t_5 \to 1\):

[5]:
p4 = NamedExpression('pbm4', _.subs(t, (t - t4) / (t5 - t4)))

If you don’t want to do the re-scaling, simply un-comment the next line!

[6]:
#p4 = NamedExpression('pbm4', b_monomial.dot(coefficients))

Either way, this is our polynomial segment …

[7]:
p4
[7]:
$\displaystyle \boldsymbol{p}_{4} = \frac{\boldsymbol{d}_{4} \left(t - t_{4}\right)^{3}}{\left(- t_{4} + t_{5}\right)^{3}} + \frac{\boldsymbol{c}_{4} \left(t - t_{4}\right)^{2}}{\left(- t_{4} + t_{5}\right)^{2}} + \frac{\boldsymbol{b}_{4} \left(t - t_{4}\right)}{- t_{4} + t_{5}} + \boldsymbol{a}_{4}$

… and it’s derivative/velocity/tangent vectors:

[8]:
pd4 = p4.diff(t)
pd4
[8]:
$\displaystyle \frac{d}{d t} \boldsymbol{p}_{4} = \frac{3 \boldsymbol{d}_{4} \left(t - t_{4}\right)^{2}}{\left(- t_{4} + t_{5}\right)^{3}} + \frac{\boldsymbol{c}_{4} \cdot \left(2 t - 2 t_{4}\right)}{\left(- t_{4} + t_{5}\right)^{2}} + \frac{\boldsymbol{b}_{4}}{- t_{4} + t_{5}}$

The next steps are very similar to what we did in the uniform case, except that we use \(t_4\) and \(t_5\) instead of \(0\) and \(1\), respectively.

[9]:
x4 = p4.evaluated_at(t, t4).with_name('xbm4')
x5 = p4.evaluated_at(t, t5).with_name('xbm5')
xd4 = pd4.evaluated_at(t, t4).with_name('xdotbm4')
xd5 = pd4.evaluated_at(t, t5).factor().with_name('xdotbm5')

To simplify things, we define a new symbol \(\Delta_4 = t_5 - t_4\), representing the duration of the current segment. However, we only use this for simplifying the display, further calculations are still carried out with \(t_i\).

[10]:
delta = {
    t5 - t4: sp.Symbol('Delta4'),
}
[11]:
display(x4, x5, xd4.subs(delta), xd5.subs(delta))
$\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} = \frac{\boldsymbol{b}_{4}}{\Delta_{4}}$
$\displaystyle \boldsymbol{\dot{x}}_{5} = \frac{\boldsymbol{b}_{4} + 2 \boldsymbol{c}_{4} + 3 \boldsymbol{d}_{4}}{\Delta_{4}}$

Basis Matrix#

In contrast to the uniform case, where the same basis matrix could be used for all segments, here we need a different matrix for each segment.

[12]:
M_H = NamedMatrix(r'{M_{\text{H},4}}', 4, 4)
[13]:
control_values_H = NamedMatrix(
    sp.Matrix([x4.name, x5.name, xd4.name, xd5.name]),
    M_H.name.I * coefficients)
control_values_H
[13]:
$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{\dot{x}}_{4}\\\boldsymbol{\dot{x}}_{5}\end{matrix}\right] = {M_{\text{H},4}}^{-1} \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right]$
[14]:
substitutions = x4, x5, xd4, xd5
[15]:
control_values_H.subs_symbols(*substitutions).subs(delta)
[15]:
$\displaystyle \left[\begin{matrix}\boldsymbol{a}_{4}\\\boldsymbol{a}_{4} + \boldsymbol{b}_{4} + \boldsymbol{c}_{4} + \boldsymbol{d}_{4}\\\frac{\boldsymbol{b}_{4}}{\Delta_{4}}\\\frac{\boldsymbol{b}_{4} + 2 \boldsymbol{c}_{4} + 3 \boldsymbol{d}_{4}}{\Delta_{4}}\end{matrix}\right] = {M_{\text{H},4}}^{-1} \left[\begin{matrix}\boldsymbol{d}_{4}\\\boldsymbol{c}_{4}\\\boldsymbol{b}_{4}\\\boldsymbol{a}_{4}\end{matrix}\right]$
[16]:
M_H.I = sp.Matrix([
    [expr.expand().coeff(c) for c in coefficients]
    for expr in control_values_H.subs_symbols(*substitutions).name])
M_H.I.subs(delta)
[16]:
$\displaystyle {M_{\text{H},4}}^{-1} = \left[\begin{matrix}0 & 0 & 0 & 1\\1 & 1 & 1 & 1\\0 & 0 & \frac{1}{\Delta_{4}} & 0\\\frac{3}{\Delta_{4}} & \frac{2}{\Delta_{4}} & \frac{1}{\Delta_{4}} & 0\end{matrix}\right]$
[17]:
pprint(_.expr)
Matrix([
[       0,        0,        0, 1],
[       1,        1,        1, 1],
[       0,        0, 1/Delta4, 0],
[3/Delta4, 2/Delta4, 1/Delta4, 0]])
[18]:
M_H.factor().subs(delta)
[18]:
$\displaystyle {M_{\text{H},4}} = \left[\begin{matrix}2 & -2 & \Delta_{4} & \Delta_{4}\\-3 & 3 & - 2 \Delta_{4} & - \Delta_{4}\\0 & 0 & \Delta_{4} & 0\\1 & 0 & 0 & 0\end{matrix}\right]$
[19]:
pprint(_.expr)
Matrix([
[ 2, -2,    Delta4,  Delta4],
[-3,  3, -2*Delta4, -Delta4],
[ 0,  0,    Delta4,       0],
[ 1,  0,         0,       0]])

Basis Polynomials#

[20]:
b_H = NamedMatrix(r'{b_{\text{H},4}}', b_monomial * M_H.expr)
b_H.factor().subs(delta).simplify().T
[20]:
$\displaystyle {b_{\text{H},4}}^{T} = \left[\begin{matrix}\left(t - 1\right)^{2} \cdot \left(2 t + 1\right)\\t^{2} \cdot \left(- 2 t + 3\right)\\\Delta_{4} t \left(t - 1\right)^{2}\\\Delta_{4} t^{2} \left(t - 1\right)\end{matrix}\right]$

Those are the non-uniform (cubic) Hermite basis functions. Not surprisingly, they are different for each segment, because generally the values \(\Delta_i\) are different in the non-uniform case.

Example Plot#

To quickly check whether the matrix \(M_{H,4}\) does what we expect, let’s plot an example segment.

[21]:
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:

[22]:
from helper import plot_spline_2d, plot_tangents_2d

The following code re-scales the parameter with t = (t - begin) / (end - begin). If you did not re-scale \(t\) in the derivation above, you’ll have to remove this line.

[23]:
class HermiteSegment:

    def __init__(self, control_values, begin, end):
        array = sp.lambdify([t4, t5], M_H.expr)(begin, end)
        self.coeffs = array @ control_values
        self.grid = begin, end

    def evaluate(self, t):
        t = np.expand_dims(t, -1)
        begin, end = self.grid
        # If you derived M_H without re-scaling t, remove the following line:
        t = (t - begin) / (end - begin)
        return t**[3, 2, 1, 0] @ self.coeffs
[24]:
vertices = [0, 0], [5, 1]
tangents = [2, 3], [0, -2]

We can simulate the uniform case by specifying a parameter range from \(0\) to \(1\):

[25]:
s1 = HermiteSegment([*vertices, *tangents], 0, 1)
[26]:
plot_spline_2d(s1, chords=False)
plot_tangents_2d(tangents, vertices)
../_images/euclidean_hermite-non-uniform_43_0.svg

But other ranges should work as well:

[27]:
s2 = HermiteSegment([*vertices, *tangents], 2.1, 5.5)
[28]:
plot_spline_2d(s2, chords=False)
plot_tangents_2d(tangents, vertices)
../_images/euclidean_hermite-non-uniform_46_0.svg

Utilizing the Uniform Basis Matrix#

If you did not re-scale \(t\) in the beginning of the derivation, you can use the matrix \(M_{H,i}\) to calculate the monomial coefficients of each segment (as shown in the example code above) and be done with it. The following simplification only applies if you did re-scale \(t\).

If you did re-scale \(t\), the basis matrix and the basis polynomials will look very similar to the uniform case, but they are not quite the same. This means that simply re-scaling the parameter is not enough to correctly use the uniform results for implementing non-uniform Hermite splines.

However, we can see that the only difference is that the components associated with \(\dot{\boldsymbol{x}}_4\) and \(\dot{\boldsymbol{x}}_5\) are simply multiplied by \(\Delta_4\). That means if we re-scale the parameter and multiply the given tangent vectors by \(\Delta_i\), we can indeed use the uniform workflow.

Just to make sure we are actually telling the truth, let’s check that the control values with scaled tangent vectors …

[29]:
control_values_H_scaled = sp.Matrix([
    x4.name,
    x5.name,
    (t5 - t4) * xd4.name,
    (t5 - t4) * xd5.name,
])
control_values_H_scaled.subs(delta)
[29]:
$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\Delta_{4} \boldsymbol{\dot{x}}_{4}\\\Delta_{4} \boldsymbol{\dot{x}}_{5}\end{matrix}\right]$

… really lead to the same result as when using the uniform basis matrix:

[30]:
sp.Eq(
    sp.simplify(M_H.expr * control_values_H.name),
    sp.simplify(sp.Matrix([
        [ 2, -2,  1,  1],
        [-3,  3, -2, -1],
        [ 0,  0,  1,  0],
        [ 1,  0,  0,  0],
    ]) * control_values_H_scaled))
[30]:
$\displaystyle \text{True}$

The following line will fail if you did not rescale \(t\):

[31]:
assert _ == True

To make a long story short, to implement a non-uniform cubic Hermite spline segment, we can simply re-scale the parameter to a range from \(0\) to \(1\) (by substituting \(t \to \frac{t - t_i}{t_{i+1} - t_i}\)), multiply both given tangent vectors by \(\Delta_i = t_{i+1} - t_i\) and then use the implementation of the uniform cubic Hermite spline segment.

Another way of looking at this is to consider the uniform polynomial segment \(\boldsymbol{u}_i(t)\) and its tangent vector (i.e. first derivative) \(\boldsymbol{u}'_i(t)\). If we want to know the tangent vector after substituting \(t \to \frac{t - t_i}{\Delta_i}\), we have to use the chain rule (with the inner derivative being \(\frac{1}{\Delta_i}\)):

\begin{equation*} \frac{d}{dt} \boldsymbol{u}_i\!\left(\frac{t-t_i}{\Delta_i}\right) = \frac{1}{\Delta_i} \boldsymbol{u}'_i\!\left(\frac{t-t_i}{\Delta_i}\right). \end{equation*}

This means the tangent vectors have been shrunk by \(\Delta_i\)! If we want to maintain the original lengths of our tangent vectors, we can simply scale them by \(\Delta_i\) beforehand.