doc/euclidean/hermite-uniform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.\n",
    "\n",
    "[back to overview](hermite.ipynb) -\n",
    "[properties](hermite-properties.ipynb) -\n",
    "[derivation (non-uniform)](hermite-non-uniform.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uniform Cubic Hermite Splines\n",
    "\n",
    "We derive the basis matrix as well as the basis polynomials\n",
    "for cubic (= degree 3) Hermite splines.\n",
    "The derivation for other degrees is left as an exercise for the reader.\n",
    "\n",
    "In this notebook,\n",
    "we consider *uniform* spline segments,\n",
    "i.e. the parameter in each segment varies from $0$ to $1$.\n",
    "The derivation for *non-uniform* cubic Hermite splines\n",
    "can be found in [a separate notebook](hermite-non-uniform.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(order='grevlex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load a few tools from [utility.py](utility.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utility import NamedExpression, NamedMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = sp.symbols('t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are considering a single cubic polynomial segment of a Hermite spline\n",
    "(which is sometimes called a *Ferguson cubic*).\n",
    "\n",
    "To simplify the indices in the following derivation,\n",
    "let's look at only one specific polynomial segment,\n",
    "let's say the fifth one.\n",
    "It goes from $\\boldsymbol{x}_4$ to $\\boldsymbol{x}_5$\n",
    "and it is referred to as $\\boldsymbol{p}_4(t)$, where $0 \\le t \\le 1$.\n",
    "The results will be easily generalizable to an arbitrary\n",
    "polynomial segment $\\boldsymbol{p}_i(t)$\n",
    "from $\\boldsymbol{x}_i$ to $\\boldsymbol{x}_{i+1}$,\n",
    "where $0 \\le t \\le 1$.\n",
    "\n",
    "The polynomial has 4 coefficients, $\\boldsymbol{a_4}$ to $\\boldsymbol{d_4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients = sp.Matrix(sp.symbols('a:dbm4')[::-1])\n",
    "coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combined with the *monomial basis* ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_monomial = sp.Matrix([t**3, t**2, t, 1]).T\n",
    "b_monomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... the coefficients form an expression\n",
    "for our polynomial segment $\\boldsymbol{p}_4(t)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p4 = NamedExpression('pbm4', b_monomial.dot(coefficients))\n",
    "p4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information about polynomials,\n",
    "see [Polynomial Parametric Curves](polynomials.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also calculate the first derivative\n",
    "(a.k.a. velocity, a.k.a. tangent vector),\n",
    "while we are at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd4 = p4.diff(t)\n",
    "pd4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate a Hermite spline segment,\n",
    "we have to provide the value of the polynomial\n",
    "at the start and end point of the segment\n",
    "(at times $t = 0$ and $t = 1$, respectively).\n",
    "We also have to provide the first derivative at those same points.\n",
    "\n",
    "\\begin{align*}\n",
    "\\boldsymbol{x}_4 &= \\left.\\boldsymbol{p}_4\\right\\rvert_{t=0}\\\\\n",
    "\\boldsymbol{x}_5 &= \\left.\\boldsymbol{p}_4\\right\\rvert_{t=1}\\\\\n",
    "\\boldsymbol{\\dot{x}}_4 &= \\left.\\frac{d}{dt}\\boldsymbol{p}_4\\right\\rvert_{t=0}\\\\\n",
    "\\boldsymbol{\\dot{x}}_5 &= \\left.\\frac{d}{dt}\\boldsymbol{p}_4\\right\\rvert_{t=1}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call those 4 values the *control values* of the segment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluating the polynomial and its derivative\n",
    "at times $0$ and $1$ leads to 4 expressions for our 4 control values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x4 = p4.evaluated_at(t, 0).with_name('xbm4')\n",
    "x5 = p4.evaluated_at(t, 1).with_name('xbm5')\n",
    "xd4 = pd4.evaluated_at(t, 0).with_name('xdotbm4')\n",
    "xd5 = pd4.evaluated_at(t, 1).with_name('xdotbm5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(x4, x5, xd4, xd5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Matrix\n",
    "\n",
    "Given an input vector of control values ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_values_H = NamedMatrix(sp.Matrix([x4.name,\n",
    "                                          x5.name,\n",
    "                                          xd4.name,\n",
    "                                          xd5.name]))\n",
    "control_values_H.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... we want to find a way to transform those into the coefficients of our cubic polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_H = NamedMatrix(r'{M_\\text{H}}', 4, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients_H = NamedMatrix(coefficients, M_H.name * control_values_H.name)\n",
    "coefficients_H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This way, we can express our previously unknown coefficients\n",
    "in terms of the given control values.\n",
    "\n",
    "However, in order to make it easy to determine\n",
    "the coefficients of the *basis matrix* $M_H$,\n",
    "we need the equation the other way around\n",
    "(by left-multiplying by the inverse):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_values_H.expr = M_H.name.I * coefficients\n",
    "control_values_H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now insert the expressions for the control values\n",
    "that we obtained above ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "substitutions = x4, x5, xd4, xd5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_values_H.subs_symbols(*substitutions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and from this equation we can directly read off\n",
    "the matrix coefficients of ${M_H}^{-1}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_H.I = sp.Matrix(\n",
    "    [[expr.coeff(cv) for cv in coefficients]\n",
    "     for expr in control_values_H.subs_symbols(*substitutions).name])\n",
    "M_H.I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same thing for copy & paste purposes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(_.expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This transforms the coefficients of the polynomial into our control values,\n",
    "but we need it the other way round,\n",
    "which we can simply get by inverting the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, for copy & paste:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(_.expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a new way to write the polynomial\n",
    "$\\boldsymbol{p}_4(t)$,\n",
    "given our four control values.\n",
    "We take those control values,\n",
    "left-multiply them by the Hermite basis matrix $M_\\text{H}$\n",
    "(which gives us a column vector of coefficients),\n",
    "which we can then left-multiply by the monomial basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.MatMul(b_monomial, M_H.expr, control_values_H.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Polynomials\n",
    "\n",
    "However, instead of calculating from right to left,\n",
    "we can also start at the left and\n",
    "multiply the monomial basis with the Hermite basis matrix $M_\\text{H}$,\n",
    "which yields (a row vector containing) the *Hermite basis polynomials*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_H = NamedMatrix(r'{b_\\text{H}}', b_monomial * M_H.expr)\n",
    "b_H.factor().T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multiplication of this row vector\n",
    "with the column vector of control values\n",
    "again produces the polynomial $\\boldsymbol{p}_4(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the basis polynomials\n",
    "with some help from [helper.py](helper.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import plot_basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_basis(*b_H.expr, labels=sp.symbols('xbm_i xbm_i+1 xdotbm_i xdotbm_i+1'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the basis function associated with $\\boldsymbol{x}_i$\n",
    "has the value $1$ at the beginning,\n",
    "while all others are $0$ at that point.\n",
    "For this reason,\n",
    "the linear combination of all basis functions at $t=0$\n",
    "simply adds up to the value $\\boldsymbol{x}_i$\n",
    "(which is exactly what we wanted to happen!).\n",
    "\n",
    "Similarly,\n",
    "the basis function associated with $\\boldsymbol{\\dot{x}}_i$\n",
    "has a first derivative of $+1$ at the beginning,\n",
    "while all others have a first derivative of $0$.\n",
    "Therefore,\n",
    "the linear combination of all basis functions at $t=0$\n",
    "turns out to have a first derivative of $\\boldsymbol{\\dot{x}}_i$\n",
    "(what a coincidence!).\n",
    "\n",
    "While $t$ progresses towards $1$,\n",
    "both functions must relinquish their influence\n",
    "to the other two basis functions.\n",
    "\n",
    "At the end (when $t=1$),\n",
    "the basis function associated with $\\boldsymbol{x}_{i+1}$\n",
    "is the only one that has a non-zero value.\n",
    "More specifically, it has the value $1$.\n",
    "Finally,\n",
    "the basis function associated with $\\boldsymbol{\\dot{x}}_{i+1}$\n",
    "is the only one with a non-zero first derivative.\n",
    "In fact, it has a first derivative of exactly $+1$\n",
    "(the function values leading up to that have to be negative\n",
    "because the final function value has to be $0$).\n",
    "\n",
    "This can be summarized by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Matrix([[\n",
    "    b.subs(t, 0),\n",
    "    b.subs(t, 1),\n",
    "    b.diff(t).subs(t, 0),\n",
    "    b.diff(t).subs(t, 1),\n",
    "] for b in b_H.expr])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Plot\n",
    "\n",
    "To quickly check whether the matrix $M_H$ does what we expect,\n",
    "let's plot an example segment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use the same API as for the other splines,\n",
    "we can reuse the helper functions for plotting\n",
    "from [helper.py](helper.py)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import plot_spline_2d, plot_tangents_2d, plot_vertices_2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class UniformHermiteSegment:\n",
    "\n",
    "    grid = 0, 1\n",
    "\n",
    "    def __init__(self, control_values):\n",
    "        self.coeffs = sp.lambdify([], M_H.expr)() @ control_values\n",
    "\n",
    "    def evaluate(self, t):\n",
    "        t = np.expand_dims(t, -1)\n",
    "        return t**[3, 2, 1, 0] @ self.coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Note\n",
    "\n",
    "The `@` operator is used here to do\n",
    "[NumPy's matrix multiplication](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html).\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertices = [0, 0], [5, 1]\n",
    "tangents = [2, 3], [0, -2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = UniformHermiteSegment([*vertices, *tangents])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_spline_2d(s, chords=False)\n",
    "plot_tangents_2d(tangents, vertices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relation to Bézier Splines\n",
    "\n",
    "Above, we were using two positions (start and end)\n",
    "and two tangent vectors (at those same two positions) as control values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_values_H.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about using four positions (and no tangent vectors) instead?\n",
    "\n",
    "Let's use the point $\\boldsymbol{\\tilde{x}}_4$ as a \"drag point\"\n",
    "(connected to $\\boldsymbol{x}_4$) that controls the tangent vector.\n",
    "Same for $\\boldsymbol{\\tilde{x}}_5$ (connected to $\\boldsymbol{x}_5$).\n",
    "\n",
    "And since the tangents looked unwieldily long in the plot above\n",
    "(compared to the effect they have on the shape of the curve),\n",
    "let's put the drag points only at a third of the length of the tangents,\n",
    "shall we?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\tilde{\\boldsymbol{x}}_4\n",
    "&=\n",
    "\\boldsymbol{x}_4 + \\frac{\\dot{\\boldsymbol{x}}_4}{3}\n",
    "\\\\\n",
    "\\tilde{\\boldsymbol{x}}_5\n",
    "&=\n",
    "\\boldsymbol{x}_5 - \\frac{\\dot{\\boldsymbol{x}}_5}{3}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_values_B = NamedMatrix(sp.Matrix([\n",
    "    x4.name,\n",
    "    sp.Symbol('xtildebm4'),\n",
    "    sp.Symbol('xtildebm5'),\n",
    "    x5.name,\n",
    "]), sp.Matrix([\n",
    "    x4.name,\n",
    "    x4.name + xd4.name / 3,\n",
    "    x5.name - xd5.name / 3,\n",
    "    x5.name,\n",
    "]))\n",
    "control_values_B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to come up with a matrix\n",
    "that transforms our good old Hermite control values\n",
    "into our new control points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_HtoB = NamedMatrix(r'{M_\\text{H$\\to$B}}', 4, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NamedMatrix(control_values_B.name, M_HtoB.name * control_values_H.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can immediately read the matrix coefficients\n",
    "off the previous expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_HtoB.expr = sp.Matrix([\n",
    "    [expr.coeff(cv) for cv in control_values_H.name]\n",
    "    for expr in control_values_B.expr])\n",
    "M_HtoB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(_.expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse of this matrix transforms our new control points\n",
    "into Hermite control values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_BtoH = NamedMatrix(r'{M_\\text{B$\\to$H}}', M_HtoB.I.expr)\n",
    "M_BtoH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(_.expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we combine $M_H$ with this new matrix,\n",
    "we get a matrix which leads us to a new set of basis polynomials\n",
    "associated with the 4 control points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B = NamedMatrix(r'{M_\\text{B}}', M_H.name * M_BtoH.name)\n",
    "M_B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B = M_B.subs_symbols(M_H, M_BtoH).doit()\n",
    "M_B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_B = NamedMatrix(r'{b_\\text{B}}', b_monomial * M_B.expr)\n",
    "b_B.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_basis(\n",
    "    *b_B.expr,\n",
    "    labels=sp.symbols('xbm_i xtildebm_i xtildebm_i+1 xbm_i+1'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those happen to be the cubic *Bernstein* polynomials and\n",
    "it turns out that we just invented *Bézier* curves!\n",
    "See [the section about Bézier splines](bezier.ipynb)\n",
    "for more information about them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We chose the additional control points to be located\n",
    "at $\\frac{1}{3}$ of the tangent vector.\n",
    "Let's quickly visualize this\n",
    "using the example from above and $M_\\text{H$\\to$B}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = sp.lambdify([], M_HtoB.expr)() @ [*vertices, *tangents]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_spline_2d(s, chords=False)\n",
    "plot_tangents_2d(tangents, vertices)\n",
    "plot_vertices_2d(points, chords=False, markeredgecolor='purple')\n",
    "plt.annotate(r'$\\quad\\tilde{\\bf{x}}_0$', points[1])\n",
    "plt.annotate(r'$\\quad\\tilde{\\bf{x}}_1$', points[2]);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}