doc/rotation/catmull-rom-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 rotation splines](index.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uniform Catmull--Rom-Like Quaternion Splines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen how to use\n",
    "[De Casteljau's algorithm with Slerp](de-casteljau.ipynb)\n",
    "to create \"cubic\" Bézier-like quaternion curve segments.\n",
    "However, if we only have a sequence of rotations to be interpolated\n",
    "and no additional Bézier control quaternions are provided,\n",
    "it would be great if we could compute\n",
    "the missing control quaternions automatically\n",
    "from neighboring quaternions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the\n",
    "[notebook about (uniform) Euclidean Catmull--Rom splines](../euclidean/catmull-rom-uniform.ipynb#Using-Bézier-Segments)\n",
    "we have already seen how this can be done for splines in Euclidean space:\n",
    "\n",
    "\\begin{align*}\n",
    "\\boldsymbol{\\tilde{x}}_i^{(+)}\n",
    "&= \\boldsymbol{x}_i + \\frac{\\boldsymbol{\\dot{x}}_i}{3} \\\\\n",
    "\\boldsymbol{\\tilde{x}}_i^{(-)}\n",
    "&= \\boldsymbol{x}_i - \\frac{\\boldsymbol{\\dot{x}}_i}{3}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the velocity vectors $\\boldsymbol{\\dot{x}}_i$\n",
    "live in the same Euclidean space as the position vectors $\\boldsymbol{x}_i$.\n",
    "We can simply add a fraction of a velocity to a position\n",
    "and we get a new position in return."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying this to rotations is unfortunately not very straightforward.\n",
    "When unit quaternions are moving along the the unit hypersphere,\n",
    "their velocity vectors are tangential to that hypersphere,\n",
    "which means that the velocity vectors are generally\n",
    "not unit quaternions themselves.\n",
    "Furthermore, adding a (non-zero length) tangent vector to a unit quaternion\n",
    "never leads to a unit quaternion as a result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of using tangent vectors,\n",
    "we can introduce a (yet unknown)\n",
    "[relative quaternion (in the global frame of reference)](quaternions.ipynb#Relative-Rotation-(Global-Frame-of-Reference))\n",
    "$q_{i,\\text{offset}}$:\n",
    "\n",
    "\\begin{align*}\n",
    "\\tilde{q}_{i}^{(+)}\n",
    "&=\n",
    "{q_{i,\\text{offset}}}^{\\frac{1}{3}} \\; q_i\\\\\n",
    "\\tilde{q}_{i}^{(-)}\n",
    "&=\n",
    "{q_{i,\\text{offset}}}^{-\\frac{1}{3}} \\; q_i\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When trying to obtain $q_{i,\\text{offset}}$,\n",
    "the problem is that there are many equivalent ways\n",
    "to write the equation for tangent vectors in Euclidean space ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_i =\n",
    "\\frac{\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_{i-1}}{2} =\n",
    "\\frac{(\\boldsymbol{x}_i - \\boldsymbol{x}_{i-1}) +\n",
    "(\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_i)\n",
    "}{2} =\n",
    "\\frac{\\boldsymbol{x}_i - \\boldsymbol{x}_{i-1}}{2} +\n",
    "\\frac{\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_i}{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... but \"translating\" them to quaternions will lead to different results!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the following experiments,\n",
    "let's define three quaternions using the `angles2quat()` function\n",
    "from [helper.py](helper.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import angles2quat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q3 = angles2quat(0, 0, 0)\n",
    "q4 = angles2quat(0, 45, -10)\n",
    "q5 = angles2quat(90, 0, -90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relative Rotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a first attempt,\n",
    "we can try to \"translate\" the equation ...\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_i =\n",
    "\\frac{\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_{i-1}}{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... to unit quaternions like this:\n",
    "\n",
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left(q_{i+1} {q_{i-1}}^{-1}\\right)^{\\frac{1}{2}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_a = q3.rotation_to(q5)**(1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll see later whether that's reasonable or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the next few examples,\n",
    "we define the\n",
    "[relative rotations](quaternions.ipynb#Relative-Rotation-(Global-Frame-of-Reference))\n",
    "associated with the the incoming and the outgoing chord:\n",
    "\n",
    "\\begin{align*}\n",
    "q_\\text{in} &= q_i {q_{i-1}}^{-1}\\\\\n",
    "q_\\text{out} &= q_{i+1} {q_i}^{-1}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_in = q3.rotation_to(q4)\n",
    "q_out = q4.rotation_to(q5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next equation ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_i =\n",
    "\\frac{(\\boldsymbol{x}_i - \\boldsymbol{x}_{i-1}) +\n",
    "(\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_i)\n",
    "}{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... can be \"translated\" to unit quaternions like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left(q_\\text{out} q_\\text{in} \\right)^{\\frac{1}{2}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_b = (q_out * q_in)**(1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that this is actually equivalent to the previous one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_b - offset_a).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the Euclidean case,\n",
    "the order doesn't matter,\n",
    "but in the quaternion case ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left(q_\\text{in} q_\\text{out} \\right)^{\\frac{1}{2}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_c = (q_in * q_out)**(1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... there is a (quite large!) difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_b - offset_c).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the equation ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_i =\n",
    "\\frac{\\boldsymbol{x}_i - \\boldsymbol{x}_{i-1}}{2} +\n",
    "\\frac{\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_i}{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... we can try another pair of equations ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left({q_\\text{out}}^{\\frac{1}{2}} {q_\\text{in}}^{\\frac{1}{2}} \\right)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_d = (q_out**(1/2) * q_in**(1/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left({q_\\text{in}}^{\\frac{1}{2}} {q_\\text{out}}^{\\frac{1}{2}} \\right)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_e = (q_in**(1/6) * q_out**(1/6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... but they are also non-symmetric:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_e - offset_d).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try a slightly more involved variant,\n",
    "where the order of $q_\\text{in}$ and $q_\\text{out}$\n",
    "can actually be reversed:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left(\n",
    "{q_\\text{out}}\n",
    "{q_\\text{in}}^{-1}\n",
    "\\right)^{\\frac{1}{2}}\n",
    "{q_\\text{in}}\n",
    "= \\left(\n",
    "{q_\\text{in}}\n",
    "{q_\\text{out}}^{-1}\n",
    "\\right)^{\\frac{1}{2}}\n",
    "{q_\\text{out}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_f = (q_out * q_in**-1)**(1/2) * q_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_g = (q_in * q_out**-1)**(1/2) * q_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_g - offset_f).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is nice to have symmetric behavior,\n",
    "but the curvature of the unit hypersphere still causes an error.\n",
    "We can check that by scaling down the components before the calculation\n",
    "(leading to a smaller curvature) and scaling up the result:\n",
    "\n",
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\left(\\left(\n",
    "{q_\\text{out}}^{\\frac{1}{10}}\n",
    "{q_\\text{in}}^{-\\frac{1}{10}}\n",
    "\\right)^{\\frac{1}{2}}\n",
    "{q_\\text{in}}^{\\frac{1}{10}}\n",
    "\\right)^{10}\n",
    "= \\left(\\left(\n",
    "{q_\\text{in}}^{\\frac{1}{10}}\n",
    "{q_\\text{out}}^{-\\frac{1}{10}}\n",
    "\\right)^{\\frac{1}{2}}\n",
    "{q_\\text{out}}^{\\frac{1}{10}}\n",
    "\\right)^{10}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_h = ((q_out**(1/10) * q_in**(-1/10))**(1/2) * q_in**(1/10))**10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_i = ((q_in**(1/10) * q_out**(-1/10))**(1/2) * q_out**(1/10))**10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_h - offset_i).xyzw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_j = ((q_out**(1/100) * q_in**(-1/100))**(1/2) * q_in**(1/100))**100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_k = ((q_in**(1/100) * q_out**(-1/100))**(1/2) * q_out**(1/100))**100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_j - offset_k).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we choose a larger scaling factor,\n",
    "the the error caused by curvature becomes smaller\n",
    "(as we will see in the next section).\n",
    "However, the numerical error gets bigger.\n",
    "We cannot scale down the components arbitrarily,\n",
    "but there is a different mathematical tool that we can use,\n",
    "which boils down to the same thing,\n",
    "as we'll see in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tangent Space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *logarithmic map* operation\n",
    "transforms a unit quaternion\n",
    "into a vector that's a member of the\n",
    "tangent space at the identity quaternion (a.k.a. $\\boldsymbol{1}$).\n",
    "In this tangent space\n",
    "-- which is a flat, three-dimensional Euclidean space --\n",
    "we can add and scale components without worrying about curvature.\n",
    "Using the *exponential map* operation,\n",
    "the result can be projected back onto the unit hypersphere.\n",
    "This way,\n",
    "we can take the equation for the tangent vector in Euclidean space ...\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\dot{x}}_i =\n",
    "\\frac{(\\boldsymbol{x}_i - \\boldsymbol{x}_{i-1}) +\n",
    "(\\boldsymbol{x}_{i+1} - \\boldsymbol{x}_i)\n",
    "}{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and \"translate\" it into unit quaternions ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "q_{i,\\text{offset}}\n",
    "\\overset{?}{=}\n",
    "\\exp\\left(\n",
    "\\frac{\\ln(q_\\text{in}) + \\ln(q_\\text{out})}{2}\n",
    "\\right)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from splines.quaternion import UnitQuaternion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_l = UnitQuaternion.exp_map((q_in.log_map() + q_out.log_map()) / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This approach is implemented in the\n",
    "[splines.quaternion.CatmullRom](../python-module/splines.quaternion.rst#splines.quaternion.CatmullRom)\n",
    "class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare this to the variants from the previous section:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_l - offset_f).xyzw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_l - offset_h).xyzw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (offset_l - offset_j).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Increasing the scaling factor from the previous section\n",
    "will get us closer and closer,\n",
    "but only until the numerical errors eventually take over."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After all those more or less successful experiments,\n",
    "let's show an example with actual rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def offset(q_1, q0, q1):\n",
    "    q_in = q0 * q_1.inverse()\n",
    "    q_out = q1 * q0.inverse()\n",
    "    return UnitQuaternion.exp_map((q_in.log_map() + q_out.log_map()) / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the \n",
    "[DeCasteljau](../python-module/splines.quaternion.rst#splines.quaternion.DeCasteljau)\n",
    "class to create a Bézier-like curve from the given control points,\n",
    "using\n",
    "[canonicalized()](../python-module/splines.quaternion.rst#splines.quaternion.canonicalized)\n",
    "to avoid angles greater than 180 degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from splines.quaternion import DeCasteljau, canonicalized"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, some helper functions from [helper.py](helper.py)\n",
    "will come in handy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import animate_rotations, display_animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't want to worry about end conditions here,\n",
    "so let's create a closed curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_closed_curve(rotations):\n",
    "    rotations = list(canonicalized(rotations + rotations[:2]))\n",
    "    control_points = []\n",
    "    for q_1, q0, q1 in zip(rotations, rotations[1:], rotations[2:]):\n",
    "        q_offset = offset(q_1, q0, q1)\n",
    "        control_points.extend([\n",
    "            q_offset**(-1/3) * q0,\n",
    "            q0,\n",
    "            q0,\n",
    "            q_offset**(1/3) * q0])\n",
    "    control_points = control_points[-2:] + control_points[:-2]\n",
    "    segments = list(zip(*[iter(control_points)] * 4))\n",
    "    return DeCasteljau(segments)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotations = [\n",
    "    angles2quat(0, 0, 180),\n",
    "    angles2quat(0, 45, 90),\n",
    "    angles2quat(90, 45, 0),\n",
    "    angles2quat(90, 90, -90),\n",
    "    angles2quat(180, 0, -180),\n",
    "    angles2quat(-90, -45, 180),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = create_closed_curve(rotations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.linspace(0, len(rotations), 200, endpoint=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ani = animate_rotations(s.evaluate(times))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_animation(ani, default_mode='loop')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shoemake's Approach\n",
    "\n",
    "In section 4.2,\n",
    "<cite data-cite-t=\"shoemake1985animating\">Shoemake (1985)</cite>\n",
    "provides two function definitions:\n",
    "\n",
    "\\begin{align*}\n",
    "\\operatorname{Double}(p, q) &= 2 (p \\cdot q) q - p\\\\\n",
    "\\operatorname{Bisect}(p, q) &= \\frac{p + q}{\\|p + q\\|}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def double(p, q):\n",
    "    return 2 * p.dot(q) * q - p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisect(p, q):\n",
    "    return (p + q).normalized()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given three successive key quaternions\n",
    "$q_{n-1}$, $q_n$ and $q_{n+1}$,\n",
    "these functions are used to compute control quaternions\n",
    "$b_n$ (controlling the incoming tangent of $q_n$) and\n",
    "$a_n$ (controlling the outgoing tangent of $q_n$):\n",
    "\n",
    "\\begin{align*}\n",
    "a_n &=\n",
    "\\operatorname{Bisect}(\\operatorname{Double}(q_{n-1}, q_n), q_{n+1})\\\\\n",
    "b_n &=\n",
    "\\operatorname{Double}(a_n, q_n)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is unclear where these equations come from,\n",
    "we only get a little hint:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> For the numerically knowledgeable, this construction\n",
    "approximates the derivative at points of a sampled function by\n",
    "averaging the central differences of the sample sequence.\n",
    ">\n",
    "> ---<cite data-cite-t=\"shoemake1985animating\">Shoemake (1985)</cite>, footnote on page 249"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shoemake_control_quaternions(q_1, q0, q1):\n",
    "    \"\"\"Shoemake's control quaternions.\n",
    "    \n",
    "    Given three key quaternions, return the control quaternions\n",
    "    preceding and following the middle one.\n",
    "    \n",
    "    Actually, the great arc distance of the returned quaternions to q0\n",
    "    still has to be reduced to 1/3 of the distance\n",
    "    to get the proper control quaternions (see the note below).\n",
    "    \n",
    "    \"\"\"\n",
    "    a = bisect(double(q_1, q0), q1)\n",
    "    b = double(a, q0).normalized()\n",
    "    return b, a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalization of $b_n$ is not explicitly mentioned in the paper,\n",
    "but even though the results have a length very close to `1.0`,\n",
    "we still have to call `normalized()` to turn the\n",
    "[Quaternion](../python-module/splines.quaternion.rst#splines.quaternion.Quaternion)\n",
    "result into a\n",
    "[UnitQuaternion](../python-module/splines.quaternion.rst#splines.quaternion.UnitQuaternion)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b, a = shoemake_control_quaternions(q3, q4, q5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are close (but by far not identical) to\n",
    "the tangent space approach from above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (a - offset_l * q4).xyzw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(map(abs, (b - offset_l.inverse() * q4).xyzw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Note\n",
    "\n",
    "Shoemake's result has to be scaled by $\\frac{1}{3}$,\n",
    "just as we did with $q_{i,\\text{offset}}$ above:\n",
    "\n",
    "> A simple check proves the curve touches $q_n$ and $q_{n+1}$ at\n",
    "its ends. A rather challenging differentiation shows it is\n",
    "tangent there to the segments determined by $a_n$ and\n",
    "$b_{n+1}$. However, as with Bézier's original curve, the\n",
    "magnitude of the tangent is three times that of the\n",
    "segment itself. That is, we are spinning three times\n",
    "faster than spherical interpolation along the arc.\n",
    "Fortunately we can correct the speed by merely\n",
    "truncating the end segments to one third their original\n",
    "length, so that $a_n$ is closer to $q_n$ and $b_{n+1}$ closer to\n",
    "$q_{n+1}$.\n",
    ">\n",
    "> ---<cite data-cite-t=\"shoemake1985animating\">Shoemake (1985)</cite>, section 4.4: \"Tangents revisited\"\n",
    "\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}