This page was generated from doc/euclidean/quadrangle.ipynb. Interactive online version: .

Quadrangle Interpolation

This doesn’t seem to be a very popular type of spline. We are mainly mentioning it because it is the starting point for interpolating rotations with Spherical Quadrangle Interpolation (Squad).

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

As usual, we import some helpers from utility.py and helper.py:

from utility import NamedExpression, NamedMatrix
from helper import plot_basis

Let’s start – as we have done before – by looking at the fifth segment of a spline, between \(\boldsymbol{x}_4\) and \(\boldsymbol{x}_5\). It will be referred to as \(\boldsymbol{p}_4(t)\), where \(0 \le t \le 1\).

x4, x5 = sp.symbols('xbm4:6')

Boehm [Boe82] mentions (on page 203) so-called quadrangle points:

x4bar = sp.symbols('xbarbm4^(+)')
x5bar = sp.symbols('xbarbm5^(-)')
x4bar, x5bar

$\displaystyle \left( \boldsymbol{\bar{x}}^{(+)}_{4}, \ \boldsymbol{\bar{x}}^{(-)}_{5}\right)$

t = sp.symbols('t')

def lerp(one, two, t):
    """Linear intERPolation.

    The parameter *t* is expected to be between 0 and 1.

    """
    return (1 - t) * one + t * two

Boehm [Boe82] also mentions (on page 210) a peculiar algorithm to construct the spline segment. In a first step, a linear interpolation between the start and end point is done, as well as a linear interpolation between the two quadrangle points. The two resulting points are then interpolated again in a second step. However, the last interpolation does not happen along a straight line, but along a parabola defined by the expression \(2t(1-t)\):

p4 = NamedExpression(
    'pbm4',
    lerp(lerp(x4, x5, t), lerp(x4bar, x5bar, t), 2 * t * (1 - t)))

This leads to a cubic polynomial. The following steps are very similar to what we did for cubic Bézier curves.

Basis Polynomials

b = [p4.expr.expand().coeff(x) for x in (x4, x4bar, x5bar, x5)]
b

$\displaystyle \left[ - 2 t^{3} + 4 t^{2} - 3 t + 1, \ 2 t^{3} - 4 t^{2} + 2 t, \ - 2 t^{3} + 2 t^{2}, \ 2 t^{3} - 2 t^{2} + t\right]$

plot_basis(*b, labels=(x4, x4bar, x5bar, x5))

Basis Matrix

M_Q = NamedMatrix(
    r'{M_\text{Q}}',
    sp.Matrix([[c.coeff(x) for x in (x4, x4bar, x5bar, x5)]
               for c in p4.as_poly(t).all_coeffs()]))
M_Q

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

M_Q.I

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

Tangent Vectors

pd4 = p4.diff(t)

xd4 = pd4.evaluated_at(t, 0)
xd4

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=0} = 2 \boldsymbol{\bar{x}}^{(+)}_{4} - 3 \boldsymbol{x}_{4} + \boldsymbol{x}_{5}$

xd5 = pd4.evaluated_at(t, 1)
xd5

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=1} = - 2 \boldsymbol{\bar{x}}^{(-)}_{5} - \boldsymbol{x}_{4} + 3 \boldsymbol{x}_{5}$

This can be generalized to:

\begin{align*} \boldsymbol{\dot{x}}^{(+)}_{i} &= 2 \boldsymbol{\bar{x}}^{(+)}_{i} - 3 \boldsymbol{x}_{i} + \boldsymbol{x}_{i+1}\\ \boldsymbol{\dot{x}}^{(-)}_{i} &= - \left(2 \boldsymbol{\bar{x}}^{(-)}_{i} - 3 \boldsymbol{x}_{i} + \boldsymbol{x}_{i-1}\right) \end{align*}

Quadrangle to Hermite Control Values

M_QtoH = NamedMatrix(
    r'{M_\text{Q$\to$H}}',
    sp.Matrix([[expr.coeff(cv) for cv in [x4, x4bar, x5bar, x5]]
               for expr in [
                   x4,
                   x5,
                   xd4.expr,
                   xd5.expr]]))
M_QtoH

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

M_QtoH.I.pull_out(sp.S.One / 2)

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

Quadrangle to Bézier Control Points

Since we already know the tangent vectors, it is easy to find the Bézier control points, as we have already shown in the notebook about uniform Hermite splines.

x4tilde = NamedExpression('xtildebm4^(+)', x4 + xd4.expr / 3)
x4tilde

$\displaystyle \boldsymbol{\tilde{x}}^{(+)}_{4} = \frac{2 \boldsymbol{\bar{x}}^{(+)}_{4}}{3} + \frac{\boldsymbol{x}_{5}}{3}$

x5tilde = NamedExpression('xtildebm5^(-)', x5 - xd5.expr / 3)
x5tilde

$\displaystyle \boldsymbol{\tilde{x}}^{(-)}_{5} = \frac{2 \boldsymbol{\bar{x}}^{(-)}_{5}}{3} + \frac{\boldsymbol{x}_{4}}{3}$

M_QtoB = NamedMatrix(
    r'{M_\text{Q$\to$B}}',
    sp.Matrix([[expr.coeff(cv) for cv in (x4, x4bar, x5bar, x5)]
               for expr in [
                   x4,
                   x4tilde.expr,
                   x5tilde.expr,
                   x5]]))
M_QtoB.pull_out(sp.S.One / 3)

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

M_QtoB.I.pull_out(sp.S.One / 2)

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

The inverse matrix can be used for converting from Bézier control points to quadrangle points:

NamedMatrix(
    sp.Matrix([x4, x4bar, x5bar, x5]),
    M_QtoB.I.expr * sp.Matrix([x4, x4tilde.name, x5tilde.name, x5]))

$\displaystyle \left[\begin{matrix}\boldsymbol{x}_{4}\\\boldsymbol{\bar{x}}^{(+)}_{4}\\\boldsymbol{\bar{x}}^{(-)}_{5}\\\boldsymbol{x}_{5}\end{matrix}\right] = \left[\begin{matrix}\boldsymbol{x}_{4}\\- \frac{\boldsymbol{x}_{5}}{2} + \frac{3 \boldsymbol{\tilde{x}}^{(+)}_{4}}{2}\\- \frac{\boldsymbol{x}_{4}}{2} + \frac{3 \boldsymbol{\tilde{x}}^{(-)}_{5}}{2}\\\boldsymbol{x}_{5}\end{matrix}\right]$

We can generalize the equations for the outgoing and incoming quadrangle points:

\begin{align*} \boldsymbol{\bar{x}}_i^{(+)} &= \frac{3}{2} \boldsymbol{\tilde{x}}_i^{(+)} - \frac{1}{2} \boldsymbol{x}_{i+1}\\ \boldsymbol{\bar{x}}_i^{(-)} &= \frac{3}{2} \boldsymbol{\tilde{x}}_i^{(-)} - \frac{1}{2} \boldsymbol{x}_{i-1} \end{align*}

The two equations are also shown by Boehm [Boe82] on page 203.

Non-Uniform Parameterization

Just like cubic Bézier splines, the shape of a segment (i.e. the image) is fully defined by its four control points. Re-scaling the parameter does not change the shape, but it changes the speed and therefore the tangent vectors.

t4, t5 = sp.symbols('t4:6')

p4nu = p4.subs(t, (t - t4) / (t5 - t4)).with_name(
    r'\boldsymbol{p}_\text{4,non-uniform}')

pd4nu = p4nu.diff(t)

pd4nu.evaluated_at(t, t4)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_\text{4,non-uniform}}\right\rvert_{t=t_{4}} = \frac{2 \boldsymbol{\bar{x}}^{(+)}_{4}}{- t_{4} + t_{5}} - \frac{3 \boldsymbol{x}_{4}}{- t_{4} + t_{5}} + \frac{\boldsymbol{x}_{5}}{- t_{4} + t_{5}}$

pd4nu.evaluated_at(t, t5)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_\text{4,non-uniform}}\right\rvert_{t=t_{5}} = - \frac{2 \boldsymbol{\bar{x}}^{(-)}_{5}}{- t_{4} + t_{5}} - \frac{\boldsymbol{x}_{4}}{- t_{4} + t_{5}} + \frac{3 \boldsymbol{x}_{5}}{- t_{4} + t_{5}}$

This can be generalized to:

\begin{align*} \boldsymbol{\dot{x}}^{(+)}_{i,\text{non-uniform}} &= \frac{2 \boldsymbol{\bar{x}}^{(+)}_{i} - 3 \boldsymbol{x}_{i} + \boldsymbol{x}_{i+1} }{ \Delta_i }\\ \boldsymbol{\dot{x}}^{(-)}_{i,\text{non-uniform}} &= -\frac{2 \boldsymbol{\bar{x}}^{(-)}_{i} - 3 \boldsymbol{x}_{i} + \boldsymbol{x}_{i-1} }{ \Delta_{i-1} } \end{align*}