{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bézier Splines\n",
    "\n",
    "See also https://pomax.github.io/bezierinfo/.\n",
    "\n",
    "There are several ways to get to Bézier curves, one was already shown in\n",
    "[the notebook about Hermite curves](hermite-uniform.ipynb#Relation-to-Bézier-Splines)\n",
    "(but only for cubic curves).\n",
    "\n",
    "TODO: first explain control polylines and then link to Hermite splines?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another one is the so-called De Casteljau's algorithm. (TODO: link to De Casteljau)\n",
    "\n",
    "One nice aspect of this is that the algorithm can be used for arbitrary polynomial degrees.\n",
    "\n",
    "A Bézier spline is defined by a so-called *control polyline* (or *control polygon*), which comprises a sequence of *control points*.\n",
    "Some of those control points are part of the final spline curve, others lie outside of it.\n",
    "The degree of a spline segment determines how many \"off-curve\" control points are between two \"on-curve\" control points.\n",
    "\n",
    "For example, in a cubic (degree = 3) Bézier spline there are two (= degree - 1) \"off-curve\" control points.\n",
    "\n",
    "Two equally valid viewpoints for what a Bézier spline is:\n",
    "\n",
    "* A sequence of curve segments, each defined by degree + 1 control points.\n",
    "  The first control point of a segment is the same as the last control point of the previous one.\n",
    "\n",
    "* A sequence of control points that can be used to shape the resulting curve.\n",
    "  Every degree'th control point lies on the curve and the others define the shape of the curve segments.\n",
    "\n",
    "TODO: most well-known: cubic Bézier splines (show screenshot from drawing program, e.g. Inkscape).\n",
    "The two \"off-curve\" control points are shown as \"handles\".\n",
    "\n",
    "TODO: typical set of constraints on continuity in drawing programs: C0, C1, G1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparations\n",
    "\n",
    "Before we continue, here are are few preparations for the following calculations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import stuff from the file [utility.py](utility.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utility import NamedExpression, NamedMatrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's prepare a few symbols for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t, x0, x1, x2, x3, x4 = sp.symbols('t, xbm:5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and a helper function for plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_curve(func, points, dots=30, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    times = np.linspace(0, 1, dots)\n",
    "    ax.plot(*func(points, times).T, '.')\n",
    "    ax.scatter(*np.asarray(points).T, marker='x', c='black')\n",
    "    ax.set_title(func.__name__ + ' Bézier curve')\n",
    "    ax.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to prepare for the animations we will see below.\n",
    "This is using code from the file [casteljau.py](casteljau.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from casteljau import create_animation\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "def show_casteljau_animation(points, frames=30, interval=200):\n",
    "    ani = create_animation(points, frames=frames)\n",
    "    display(HTML(ani.to_jshtml(default_mode='reflect')))\n",
    "    plt.close()  # avoid spurious figure display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degree 1, a.k.a. linear\n",
    "\n",
    "But let's start with the trivial case:\n",
    "A Bézier spline of degree 1 is just a piecewise linear curve connecting all the control points.\n",
    "There are no \"off-curve\" control points that could bend the curve segments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that we have two control points, $\\boldsymbol{x}_0$ and $\\boldsymbol{x}_1$ ...\n",
    "\n",
    "... linear equation ...:\n",
    "\n",
    "\\begin{equation}\n",
    "\\boldsymbol{p}_{0,1}(t) = \\boldsymbol{x}_0 + t (\\boldsymbol{x}_1 - \\boldsymbol{x}_0)\n",
    "\\end{equation}\n",
    "\n",
    "... in other words ... this is called *affine combination*, but we don't really have to worry about it ...\n",
    "\n",
    "\\begin{equation}\n",
    "\\boldsymbol{p}_{0,1}(t) = (1 - t) \\boldsymbol{x}_0 + t \\boldsymbol{x}_1\n",
    "\\end{equation}\n",
    "\n",
    "... with $t \\in [0, 1]$ (which is called *uniform*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: show change of variables for *non-uniform* curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we will be needing quite a bunch of those affine combinations, let's create a helper function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affine_combination(one, two):\n",
    "    return (1 - t) * one + t * two"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define the equation in SymPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p01 = NamedExpression('pbm_0,1', affine_combination(x0, x1))\n",
    "p01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1 = [p01.expr.expand().coeff(x.name).factor() for x in (x0, x1)]\n",
    "b1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doesn't look like much, but those are the Bernstein bases for degree 1 (<https://en.wikipedia.org/wiki/Bernstein_polynomial>).\n",
    "\n",
    "It doesn't get much more interesting if we plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.plot(*b1, (t, 0, 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to convert this to coefficients for the monomial basis $[t, 1]$ instead of the Bernstein basis functions, you can use this matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B1 = NamedMatrix(\n",
    "    r'{M_\\text{B}^{(1)}}',\n",
    "    sp.Matrix([[c.coeff(x) for x in (x0, x1)]\n",
    "               for c in p01.expr.as_poly(t).all_coeffs()]))\n",
    "M_B1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying this matrix leads to the coefficients of the linear equation mentioned in the beginning of this section\n",
    "($\\boldsymbol{p}_{0,1}(t) = t (\\boldsymbol{x}_1 - \\boldsymbol{x}_0) + \\boldsymbol{x}_0$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.MatMul(M_B1.expr, sp.Matrix([x0, x1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you ever need that, here's the inverse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B1.I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anywho, let's calculate points on the curve by using the Bernstein basis functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear(points, times):\n",
    "    \"\"\"Evaluate linear Bézier curve (given by two points) at given times.\"\"\"\n",
    "    return np.column_stack(sp.lambdify(t, b1)(times)) @ points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [\n",
    "    (0, 0),\n",
    "    (1, 0.5),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_curve(linear, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_casteljau_animation(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I know, not very exciting. But it gets better!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degree 2, a.k.a. quadratic\n",
    "\n",
    "Consider three control points, $\\boldsymbol{x}_0$, $\\boldsymbol{x}_1$ and $\\boldsymbol{x}_2$ ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the affine combinations of the first two points from above ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and we do the same thing for the second and third point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p12 = NamedExpression('pbm_1,2', affine_combination(x1, x2))\n",
    "p12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we make another affine combination of those two results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p02 = NamedExpression('pbm_0,2', affine_combination(p01.expr, p12.expr))\n",
    "p02"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bernstein basis functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b2 = [p02.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2)]\n",
    "b2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.plot(*b2, (t, 0, 1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B2 = NamedMatrix(\n",
    "    r'{M_\\text{B}^{(2)}}',\n",
    "    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2)]\n",
    "               for c in p02.expr.as_poly(t).all_coeffs()]))\n",
    "M_B2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B2.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic(points, times):\n",
    "    \"\"\"Evaluate quadratic Bézier curve (given by three points) at given times.\"\"\"\n",
    "    return np.column_stack(sp.lambdify(t, b2)(times)) @ points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [\n",
    "    (0, 0),\n",
    "    (0.2, 0.5),\n",
    "    (1, -0.3),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_curve(quadratic, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_casteljau_animation(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some more insight, let's look at the first derivative of the curve (i.e. the tangent vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v02 = p02.diff(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... at the beginning and the end of the curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v02.evaluated_at(t, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v02.evaluated_at(t, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that the tangent vector at the beginning and end of the curve is parallel to the line\n",
    "from $\\boldsymbol{x}_0$ to $\\boldsymbol{x}_1$ and\n",
    "from $\\boldsymbol{x}_1$ to $\\boldsymbol{x}_2$, respectively.\n",
    "The length of the tangent vectors is twice the length of those lines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might have already seen that coming, but it turns out that the last line in de Casteljau's algorithm ($\\boldsymbol{p}_{1,2}(t) - \\boldsymbol{p}_{0,1}(t)$ in our case) is exactly half of the tangent vector (at any given $t \\in [0, 1]$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(v02.expr - 2 * (p12.expr - p01.expr)).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you are wondering, the factor 2 comes from the degree 2 of our quadratic curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degree 3, a.k.a. cubic\n",
    "\n",
    "Consider four control points, $\\boldsymbol{x}_0$, $\\boldsymbol{x}_1$, $\\boldsymbol{x}_2$ and $\\boldsymbol{x}_3$ ...\n",
    "\n",
    "By now, the pattern should be clear: We take the result from the first three points from above and affine-combine it with the result for the three points $\\boldsymbol{x}_1$, $\\boldsymbol{x}_2$ and $\\boldsymbol{x}_3$.\n",
    "\n",
    "Combination of $\\boldsymbol{x}_2$ and $\\boldsymbol{x}_3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p23 = NamedExpression('pbm_2,3', affine_combination(x2, x3))\n",
    "p23"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combination of $\\boldsymbol{x}_1$, $\\boldsymbol{x}_2$ and $\\boldsymbol{x}_3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p13 = NamedExpression('pbm_1,3', affine_combination(p12.expr, p23.expr))\n",
    "p13"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combination of $\\boldsymbol{x}_0$, $\\boldsymbol{x}_1$, $\\boldsymbol{x}_2$ and $\\boldsymbol{x}_3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p03 = NamedExpression('pbm_0,3', affine_combination(p02.expr, p13.expr))\n",
    "p03"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bernstein bases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b3 = [p03.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2, x3)]\n",
    "b3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: show that those are the same Bernstein bases as in the notebook about Hermite splines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.plot(*b3, (t, 0, 1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B3 = NamedMatrix(\n",
    "    r'{M_\\text{B}^{(3)}}',\n",
    "    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2, x3)]\n",
    "               for c in p03.expr.as_poly(t).all_coeffs()]))\n",
    "M_B3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B3.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cubic(points, times):\n",
    "    \"\"\"Evaluate cubic Bézier curve (given by four points) at given times.\"\"\"\n",
    "    return np.column_stack(sp.lambdify(t, b3)(times)) @ points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [\n",
    "    (0, 0.3),\n",
    "    (0.2, 0.5),\n",
    "    (0.1, 0),\n",
    "    (1, 0.2),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_curve(cubic, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_casteljau_animation(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, let's look at the derivative (i.e. the tangent vector) of the curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v03 = p03.diff(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... at the beginning and the end of the curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v03.evaluated_at(t, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v03.evaluated_at(t, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that the tangent vector at the beginning and end of the curve is parallel to the line\n",
    "from $\\boldsymbol{x}_0$ to $\\boldsymbol{x}_1$ and\n",
    "from $\\boldsymbol{x}_2$ to $\\boldsymbol{x}_3$, respectively.\n",
    "The length of the tangent vectors is three times the length of those lines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now see that the last line in de Casteljau's algorithm ($\\boldsymbol{p}_{1,3}(t) - \\boldsymbol{p}_{0,2}(t)$ in this case) is exactly a third of the tangent vector (at any given $t \\in [0, 1]$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(v03.expr - 3 * (p13.expr - p02.expr)).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the factor 3 comes from the degree 3 of our curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now know the tangent vectors at the beginning and the end of the curve, and obviously we know the values of the curve at the beginning and the end:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p03.evaluated_at(t, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p03.evaluated_at(t, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these four pieces of information, we can find a transformation from the four Bézier control points to the two control points and two tangent vectors of Hermite splines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_BtoH = NamedMatrix(\n",
    "    r'{M_\\text{B$\\to$H}}',\n",
    "    sp.Matrix([[expr.coeff(cv) for cv in [x0, x1, x2, x3]]\n",
    "               for expr in [\n",
    "                   x0,\n",
    "                   x3,\n",
    "                   v03.evaluated_at(t, 0).expr,\n",
    "                   v03.evaluated_at(t, 1).expr]]))\n",
    "M_BtoH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can simply invert this if we want to go in the other direction, from Hermite to Bézier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_BtoH.I.pull_out(sp.S.One / 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, those are the same matrices as shown in the [notebook about uniform cubic Hermite splines](hermite-uniform.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: show tangent vectors for non-uniform case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degree 4, a.k.a. quartic\n",
    "\n",
    "Consider five control points, $\\boldsymbol{x}_0$, $\\boldsymbol{x}_1$, $\\boldsymbol{x}_2$, $\\boldsymbol{x}_3$ and $\\boldsymbol{x}_4$ ...\n",
    "\n",
    "More combinations!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p34 = NamedExpression('pbm_3,4', affine_combination(x3, x4))\n",
    "p24 = NamedExpression('pbm_2,4', affine_combination(p23.expr, p34.expr))\n",
    "p14 = NamedExpression('pbm_1,4', affine_combination(p13.expr, p24.expr))\n",
    "p04 = NamedExpression('pbm_0,4', affine_combination(p03.expr, p14.expr))\n",
    "p04"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kinda long, but anyway, let's try to extract the Bernstein bases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b4 = [p04.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2, x3, x4)]\n",
    "b4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.plot(*b4, (t, 0, 1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B4 = NamedMatrix(\n",
    "    '{M_B^{(4)}}',\n",
    "    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2, x3, x4)]\n",
    "               for c in p04.expr.as_poly(t).all_coeffs()]))\n",
    "M_B4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_B4.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quartic(points, times):\n",
    "    \"\"\"Evaluate quartic Bézier curve (given by five points) at given times.\"\"\"\n",
    "    return np.column_stack(sp.lambdify(t, b4)(times)) @ points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [\n",
    "    (0, 0),\n",
    "    (0.5, 0),\n",
    "    (0.7, 1),\n",
    "    (1, 1.5),\n",
    "    (-1, 1),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_curve(quartic, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_casteljau_animation(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For completeness' sake, let's look at the derivative (i.e. the tangent vector) of the curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v04 = p04.diff(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... at the beginning and the end of the curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v04.evaluated_at(t, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v04.evaluated_at(t, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By now it shouldn't be surprising that the tangent vector at the beginning and end of the curve is parallel to the line\n",
    "from $\\boldsymbol{x}_0$ to $\\boldsymbol{x}_1$ and\n",
    "from $\\boldsymbol{x}_3$ to $\\boldsymbol{x}_4$, respectively.\n",
    "The length of the tangent vectors is four times the length of those lines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line in de Casteljau's algorithm ($\\boldsymbol{p}_{1,4}(t) - \\boldsymbol{p}_{0,3}(t)$ in this case) is exactly a fourth of the tangent vector (at any given $t \\in [0, 1]$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(v04.expr - 4 * (p14.expr - p03.expr)).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the factor 4 comes from the degree 4 of our curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arbitrary Degree\n",
    "\n",
    "We could go on doing this for higher and higher degrees, but this would get more and more annoying.\n",
    "\n",
    "Luckily, there is a closed formula available to calculate Bernstein polynomials for an arbitrary degree $n$!\n",
    "\n",
    "\\begin{equation}\n",
    "b_{i,n}(x) = {n \\choose i} x^i \\left( 1 - x \\right)^{n - i}, \\quad i = 0, \\ldots, n.\n",
    "\\end{equation}\n",
    "\n",
    "with the *binomial coefficient* ${n \\choose i} = \\frac{n!}{i!(n - i)!}$.\n",
    "\n",
    "TODO: link to proof?\n",
    "\n",
    "TODO: show Bernstein polynomials for \"quintic\" etc.?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_casteljau_animation([\n",
    "    (0, 0),\n",
    "    (-1, 1),\n",
    "    (-0.5, 2),\n",
    "    (1, 2.5),\n",
    "    (2, 2),\n",
    "    (2, 1.5),\n",
    "    (0.5, 0.5),\n",
    "    (1, -0.5),\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}