doc/euclidean/quadrangle.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 Euclidean splines](index.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quadrangle Interpolation\n",
    "\n",
    "This doesn't seem to be a very popular type of spline.\n",
    "We are mainly mentioning it because it is\n",
    "the starting point for interpolating rotations with\n",
    "[Spherical Quadrangle Interpolation (Squad)](../rotation/squad.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": [
    "As usual, we import some helpers from\n",
    "[utility.py](utility.py) and [helper.py](helper.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utility import NamedExpression, NamedMatrix\n",
    "from helper import plot_basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start\n",
    "-- as we have done before --\n",
    "by looking at the fifth segment of a spline,\n",
    "between $\\boldsymbol{x}_4$ and $\\boldsymbol{x}_5$.\n",
    "It will be referred to as $\\boldsymbol{p}_4(t)$,\n",
    "where $0 \\le t \\le 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x4, x5 = sp.symbols('xbm4:6')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<cite data-cite-t=\"boehm1982cubics\">Boehm (1982)</cite>\n",
    "mentions (on page 203) so-called\n",
    "*quadrangle points*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x4bar = sp.symbols('xbarbm4^(+)')\n",
    "x5bar = sp.symbols('xbarbm5^(-)')\n",
    "x4bar, x5bar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = sp.symbols('t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lerp(one, two, t):\n",
    "    \"\"\"Linear intERPolation.\n",
    "\n",
    "    The parameter *t* is expected to be between 0 and 1.\n",
    "\n",
    "    \"\"\"\n",
    "    return (1 - t) * one + t * two"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<cite data-cite-t=\"boehm1982cubics\">Boehm (1982)</cite>\n",
    "also mentions (on page 210)\n",
    "a peculiar algorithm\n",
    "to construct the spline segment.\n",
    "In a first step,\n",
    "a linear interpolation between the start and end point is done,\n",
    "as well as a linear interpolation between the two quadrangle points.\n",
    "The two resulting points are then interpolated again in a second step.\n",
    "However, the last interpolation does not happen along a straight line,\n",
    "but along a parabola defined by the expression $2t(1-t)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p4 = NamedExpression(\n",
    "    'pbm4',\n",
    "    lerp(lerp(x4, x5, t), lerp(x4bar, x5bar, t), 2 * t * (1 - t)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads to a cubic polynomial.\n",
    "The following steps are very similar\n",
    "to what we did for\n",
    "[cubic Bézier curves](bezier-de-casteljau.ipynb#Degree-3-(Cubic))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = [p4.expr.expand().coeff(x) for x in (x4, x4bar, x5bar, x5)]\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_basis(*b, labels=(x4, x4bar, x5bar, x5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_Q = NamedMatrix(\n",
    "    r'{M_\\text{Q}}',\n",
    "    sp.Matrix([[c.coeff(x) for x in (x4, x4bar, x5bar, x5)]\n",
    "               for c in p4.as_poly(t).all_coeffs()]))\n",
    "M_Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_Q.I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tangent Vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd4 = p4.diff(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd4 = pd4.evaluated_at(t, 0)\n",
    "xd4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xd5 = pd4.evaluated_at(t, 1)\n",
    "xd5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be generalized to:\n",
    "\n",
    "\\begin{align*}\n",
    "\\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)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadrangle to Hermite Control Values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_QtoH = NamedMatrix(\n",
    "    r'{M_\\text{Q$\\to$H}}',\n",
    "    sp.Matrix([[expr.coeff(cv) for cv in [x4, x4bar, x5bar, x5]]\n",
    "               for expr in [\n",
    "                   x4,\n",
    "                   x5,\n",
    "                   xd4.expr,\n",
    "                   xd5.expr]]))\n",
    "M_QtoH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_QtoH.I.pull_out(sp.S.One / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadrangle to Bézier Control Points\n",
    "\n",
    "Since we already know the tangent vectors,\n",
    "it is easy to find the Bézier control points,\n",
    "as we have already shown in the notebook about\n",
    "[uniform Hermite splines](hermite-uniform.ipynb#Relation-to-Bézier-Splines)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x4tilde = NamedExpression('xtildebm4^(+)', x4 + xd4.expr / 3)\n",
    "x4tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x5tilde = NamedExpression('xtildebm5^(-)', x5 - xd5.expr / 3)\n",
    "x5tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_QtoB = NamedMatrix(\n",
    "    r'{M_\\text{Q$\\to$B}}',\n",
    "    sp.Matrix([[expr.coeff(cv) for cv in (x4, x4bar, x5bar, x5)]\n",
    "               for expr in [\n",
    "                   x4,\n",
    "                   x4tilde.expr,\n",
    "                   x5tilde.expr,\n",
    "                   x5]]))\n",
    "M_QtoB.pull_out(sp.S.One / 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_QtoB.I.pull_out(sp.S.One / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse matrix can be used for converting\n",
    "from Bézier control points to quadrangle points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NamedMatrix(\n",
    "    sp.Matrix([x4, x4bar, x5bar, x5]),\n",
    "    M_QtoB.I.expr * sp.Matrix([x4, x4tilde.name, x5tilde.name, x5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generalize the equations\n",
    "for the outgoing and incoming quadrangle points:\n",
    "\n",
    "\\begin{align*}\n",
    "\\boldsymbol{\\bar{x}}_i^{(+)} &=\n",
    "\\frac{3}{2} \\boldsymbol{\\tilde{x}}_i^{(+)} -\n",
    "\\frac{1}{2} \\boldsymbol{x}_{i+1}\\\\\n",
    "\\boldsymbol{\\bar{x}}_i^{(-)} &=\n",
    "\\frac{3}{2} \\boldsymbol{\\tilde{x}}_i^{(-)} -\n",
    "\\frac{1}{2} \\boldsymbol{x}_{i-1}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two equations are also shown by\n",
    "<cite data-cite-t=\"boehm1982cubics\">Boehm (1982)</cite>\n",
    "on page 203."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-Uniform Parameterization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like\n",
    "[cubic Bézier splines](bezier-non-uniform.ipynb),\n",
    "the shape of a segment (i.e. the [image](https://en.wikipedia.org/wiki/Image_(mathematics))) is fully defined by its four control points.\n",
    "Re-scaling the parameter does not change the shape,\n",
    "but it changes the speed and therefore the tangent vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t4, t5 = sp.symbols('t4:6')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p4nu = p4.subs(t, (t - t4) / (t5 - t4)).with_name(\n",
    "    r'\\boldsymbol{p}_\\text{4,non-uniform}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd4nu = p4nu.diff(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd4nu.evaluated_at(t, t4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd4nu.evaluated_at(t, t5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be generalized to:\n",
    "\n",
    "\\begin{align*}\n",
    "\\boldsymbol{\\dot{x}}^{(+)}_{i,\\text{non-uniform}} &=\n",
    "\\frac{2 \\boldsymbol{\\bar{x}}^{(+)}_{i} - 3 \\boldsymbol{x}_{i} + \\boldsymbol{x}_{i+1}\n",
    "}{\n",
    "\\Delta_i\n",
    "}\\\\\n",
    "\\boldsymbol{\\dot{x}}^{(-)}_{i,\\text{non-uniform}} &= \n",
    "-\\frac{2 \\boldsymbol{\\bar{x}}^{(-)}_{i} - 3 \\boldsymbol{x}_{i} + \\boldsymbol{x}_{i-1}\n",
    "}{\n",
    "\\Delta_{i-1}\n",
    "}\n",
    "\\end{align*}"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}