{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Courbes de Bézier\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "19 mars 2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Les polynômes de Bernstein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout entier naturel $n$ et tout entier $k\\in[\\![0,n]\\!]$, notons\n",
    "\n",
    "$$B_k^n=\\binom n kX^k(1-X)^{n-k}$$\n",
    "\n",
    "Les polynômes $B_k^n$ sont les *polynômes de Bernstein*. Ces polynômes sont inéressants à plus d'un titre, ils interviendront dans tout le notebook. Commençons par une propriété facile.\n",
    "\n",
    "**Proposition.** Pour tout $n\\in\\mathbb N$,\n",
    "\n",
    "$$\\sum_{k=0}^n B_k^n=1$$\n",
    "\n",
    "**Démonstration.**\n",
    "\n",
    "$$\\sum_{k=0}^n B_k^n=\\sum_{k=0}^n \\binom n k X^k(1-X)^{n-k}=(X+(1-X))^n=1\\ \\diamond$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici un résultat moins évident.\n",
    "\n",
    "**Proposition.** Soit $n\\in\\mathbb N$. La famille $(B_k^n)_{0\\le k\\le n}$ est une base de l'espace vectoriel $\\mathbb R_n[X]$ des polynômes de degré inférieur ou égal à $n$.\n",
    "\n",
    "**Démonstration.** Montrons-le par récurrence sur $n$. Pour $n=0$, c'est évident car $B_0^0=1$. Soit $n\\ge 1$. Supposons la propriété vraie pour $n-1$. Remarquons que si $1\\le k\\le n-1$, on a\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "\\left(B_k^n\\right)'&=&k\\binom n k X^{k-1}(1-X)^{n-k}-(n-k)\\binom n k X^k(1-X)^{n-k-1}\\\\\n",
    "&=&n\\binom {n-1} {k-1} X^{k-1}(1-X)^{n-k}-n\\binom {n-1} k X^k(1-X)^{n-k-1}\\\\\n",
    "&=&nB_{k-1}^{n-1}-nB_k^{n-1}\n",
    "\\end{array}$$\n",
    "\n",
    "De plus,\n",
    "\n",
    "$$\\left(B_0^n\\right)'=-nB_0^{n-1}$$\n",
    "\n",
    "et\n",
    "\n",
    "$$\\left(B_n^n\\right)'=nB_{n-1}^{n-1}$$\n",
    "\n",
    "Soient $\\lambda_0,\\ldots,\\lambda_n$ $n+1$ réels. Supposons que\n",
    "\n",
    "$$(1)\\ \\sum_{k=0}^n\\lambda_kB_k^n=0$$\n",
    "\n",
    "En dérivant et en divisant par $n$, on obtient\n",
    "\n",
    "$$-\\lambda_0B_0^{n-1}+\\lambda_nB_{n-1}^{n-1}+\\sum_{k=1}^{n-1}\\lambda_k\\left(B_{k-1}^{n-1}-B_k^{n-1}\\right)=0$$\n",
    "\n",
    "Ceci peut encore s'écrire, avec un petit changement d'indice,\n",
    "\n",
    "$$-\\lambda_0B_0^{n-1}+\\lambda_nB_{n-1}^{n-1}-\\sum_{k=1}^{n-1}\\lambda_k B_k^{n-1}+\\sum_{k=0}^{n-2}\\lambda_{k+1}B_{k}^{n-1}=0$$\n",
    "\n",
    "Par l'hypothèse de récurrence, la famille $(B_k^{n-1})_{0\\le k\\le n-1}$ est libre. On en déduit que pour $0\\le k\\le n-1$, le coefficient de $B_k^{n-1}$ dans l'égalité ci-dessus est nul. Il en résulte facilement que\n",
    "\n",
    "$$\\lambda_0=\\ldots=\\lambda_n$$\n",
    "\n",
    "En reportant dans $(1)$, on en déduit que\n",
    "\n",
    "$$\\lambda_0\\sum_{k=0}^n B_k^n=\\lambda_0=0$$\n",
    "\n",
    "La famille $(B_k^n)_{0\\le k\\le n}$ est donc libre. Comme elle est de cardinal $n+1=\\dim\\mathbb R_n[X]$, c'est une base de $\\mathbb R_n[X]$. $\\diamond$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `binom` prend en paramètres deux entiers $n$ et $k$. Elle renvoie le coefficient binomial $\\binom n k$. Plus exactement, elle en renvoie une approximation flottante."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binom(n, k):\n",
    "    p = 1\n",
    "    for i in range(k):\n",
    "        p *= (n - i) / (i + 1)\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([binom(10, k) for k in range(11)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `bernstein` renvoie $B_k^n(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bernstein(n, k, t):\n",
    "    return binom(n, k) * t ** k * (1 - t) ** (n -k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons les courbes des polynômes de Bernstein sur $[0,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subdi(a, b, n):\n",
    "    d = (b - a) / n\n",
    "    return [a + k * d for k in range(n + 1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bernstein(n, rng=None):\n",
    "    ts = subdi(0, 1, 500)\n",
    "    if rng == None: rng = range(n + 1)\n",
    "    for k in rng:\n",
    "        bs = [bernstein(n, k, t) for t in ts]\n",
    "        plt.plot(ts, bs, 'k')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bernstein(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici les mêmes courbes, en supprimant celles de $B_0^n$ et $B_n^n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bernstein(20, range(1, 20))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Courbes de Bézier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $n\\in\\mathbb N$. Soient $P_0,\\ldots,P_n$ $n+1$ points du plan. Soit $t\\in[0,1]$. Considérons le point $B_{P_0,\\ldots,P_n}(t)$ défini par\n",
    "\n",
    "$$\\overrightarrow{OB_{P_0,\\ldots,P_n}(t)}=\\sum_{k=0}^nB_k^n(t)\\overrightarrow{O P_k}$$\n",
    "\n",
    "On définit ainsi une fonction\n",
    "\n",
    "$$B_{P_0,\\ldots,P_n}:[0,1]\\longrightarrow \\mathbb R^2$$\n",
    "\n",
    "Lorsqu'il n'y a pas d'ambiguïté sur les points $P_i$, nous noterons plus simplement $B(t)$ au lieu de $B_{P_0,\\ldots,P_n}(t)$. Pour alléger encore les notations, nous écrirons\n",
    "\n",
    "$$B(t)=\\sum_{k=0}^nB_k^n(t)P_k$$\n",
    "\n",
    "Les points $P_0,\\ldots,P_n$ sont appelés *points de contrôle*. Lorsque $t$ décrit le segment $[0,1]$, le point $B(t)$ décrit une courbe qui s'appelle la *courbe de Bézier* associée aux points de contrôle $P_0,\\ldots,P_n$.\n",
    "\n",
    "Les courbes de Bézier possèdent des propriétés tout à fait remarquables que nous allons explorer. Avant tout, remarquons les évidences. L'abscisse et l'ordonnée de $B(t)$ sont des fonctions polynomiales en $t$, combinaisons linéaires des $B_k^n$ qui forment une base de $\\mathbb R_n[X]$. Il est donc possible d'obtenir de cette manière, par un choix judicieux des points de contrôle, *toutes* les fonctions de $[0,1]$ vers $\\mathbb R^2$ dont l'abscisse et l'ordonnée sont des fonctions polynomiales.\n",
    "\n",
    "Remarquons aussi que $B(0)=P_0$ et $B(1)=P_n$. En revanche, sauf coïncidences, les points $P_1,\\ldots,P_{n-1}$ n'appartiennent pas à la courbe de Bézier.\n",
    "\n",
    "Regardons ce qui se passe pour $n=0$ et $n=1$.\n",
    "\n",
    "- Si $n=0$, on a pour tout $t\\in[0,1]$, $B(t)=P_0$. La fonction $B$ est constante.\n",
    "- Si $n=1$, on a pour tout $t\\in[0,1]$, $B(t)=(1-t)P_0+tP_1$. La courbe de Bézier est le segment $[P_0,P_1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `bezier0` prend en paramètres une liste `Ps` de points du plan (représentés chacun par un couple de flottants) et un réel $t\\in[0,1]$. Elle renvoie $B(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bezier0(Ps, t):\n",
    "    x = 0\n",
    "    y = 0\n",
    "    n = len(Ps) - 1\n",
    "    for k in range(n + 1):\n",
    "        xk, yk = Ps[k]\n",
    "        m =  bernstein(n, k, t)\n",
    "        x += m * xk\n",
    "        y += m * yk\n",
    "    return (x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 2), (4, 1), (4, 0)]\n",
    "print(bezier0(Ps, 0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Tracé dezs courbes de Bézier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `plot_bezier1` prend en paramètre une liste de points du plan. Si le paramètre `pts` vaut `True` elle affiche ces points. Elle affiche également la courbe de Bézier dont ces points sont les points de contrôle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bezier1(Ps, pts=True):\n",
    "    ts = subdi(0, 1, 500)\n",
    "    Qs = [bezier0(Ps, t) for t in ts]\n",
    "    xs = [x for (x, y) in Qs]\n",
    "    ys = [y for (x, y) in Qs]\n",
    "    plt.plot(xs, ys, 'k')\n",
    "    if pts:\n",
    "        for (x, y) in Ps:\n",
    "            plt.plot([x], [y], 'or')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 2), (4, 1), (4, 0)]\n",
    "plot_bezier1(Ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici une fonction `cycle` qui renvoie une liste de $n+1$ points situés sur une « spirale »."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cycle(n):\n",
    "    Ps = []\n",
    "    for k in range(n + 1):\n",
    "        t = - math.pi / 2 + 2 * k * math.pi / n\n",
    "        Ps.append((math.cos(t) / (k + 1), math.sin(t) / (k + 1)))\n",
    "    return Ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bezier1(cycle(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Quelques remarques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deux remarques concernant notre fonction `bezier0` :\n",
    "\n",
    "- Lorsque $n$ est grand les coefficients binomiaux deviennent **très** grands. En fait, **trop** grands. Les erreurs d'arrondis sur les sommes deviennent prépondérandes et `bezier0` renvoie des résultats faux. Essayez donc de prendre $n=200$ ci-dessus. Vous verrez que la courbe a « disparu ».\n",
    "- Le calcul de $B(t)$ nécessite $n$ additions et, pour chaque terme, le calcul d'un coefficient binomial. La complexité de ce calcul est $O(n^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Premières propriétés"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Barycentres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelons la formule\n",
    "\n",
    "$$B(t)=\\sum_{k=0}^nB_k^n(t)P_k$$\n",
    "\n",
    "Pour tout $t\\in[0,1]$, $B_k^n(t)\\ge 0$ et nous avons vu que\n",
    "\n",
    "$$\\sum_{k=0}^n B_k^n(t)=1$$\n",
    "\n",
    "Ainsi, le point $B(t)$ est un barycentre à coefficients positifs des points de contrôle $P_k$. Il en résulte que la courbe de Bézier est incluse dans *l'enveloppe convexe* de l'ensemble $\\{P_0,\\ldots,P_n\\}$, c'est à dire la plus petite partie convexe du plan qui contient les points de contrôle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Tangentes aux extrémités"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout $k\\in[\\![0,n]\\!]$, $0$ est racine d'ordre $k$ du polynôme $B_k^n$. On a donc pour tout $k\\ge 2$, $\\left(B_k^n\\right)'(0)=0$. De plus,\n",
    "\n",
    "$$\\left(B_0^n\\right)'(t)=-n(1-t)^{n-1}$$\n",
    "\n",
    "et \n",
    "\n",
    "$$\\left(B_1^n\\right)'(t)=n(1-t)^{n-1}-n(n-1)t(1-t)^{n-2}=n(1-t)^{n-2}(1-nt)$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$B'(0)=-nP_0+nP_1=n\\overrightarrow{P_0P_1}$$\n",
    "\n",
    "Ainsi, la tangente à la courbe de Bézier au point $P_0=B(0)$ est dirigée par le vecteur $\\overrightarrow{P_0P_1}$.\n",
    "\n",
    "Des calculs analogues montrent que la tangente à la courbe de Bézier au point $P_n=B(1)$ est dirigée par le vecteur $\\overrightarrow{P_{n-1}P_n}$.\n",
    "\n",
    "Nous commençons à comprendre pourquoi les points $P_k$ sont appelés des *points de contrôle*, en tout cas pour 4 de ces points.\n",
    "\n",
    "- Le point $P_0$ est l'origine de la courbe de Bézier.\n",
    "- Le point $P_1$ contrôle la tangente à l'origine.\n",
    "- Le point $P_n$ est l'extrémité de la courbe de Bézier.\n",
    "- Le point $P_{n-1}$ contrôle la tangente à l'extrémité.\n",
    "\n",
    "Remarquons que si $n=3$ nous avons tout compris. Évidemment, si $n\\ge 4$ il reste encore beaucoup à dire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `plot_bezier2` trace la courbe de Bézier, les points de contrôle, et les tangentes aux extrémités de la courbe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bezier2(Ps):\n",
    "    n = len(Ps) - 1\n",
    "    ts = subdi(0, 1, 500)\n",
    "    Qs = [bezier0(Ps, t) for t in ts]\n",
    "    xs = [x for (x, y) in Qs]\n",
    "    ys = [y for (x, y) in Qs]\n",
    "    plt.plot(xs, ys, 'k')\n",
    "    for (x, y) in Ps:\n",
    "        plt.plot([x], [y], 'or')\n",
    "    a, b = Ps[0]\n",
    "    c, d = Ps[1]\n",
    "    plt.plot([a, c], [b, d], 'b')\n",
    "    a, b = Ps[n - 1]\n",
    "    c, d = Ps[n]\n",
    "    plt.plot([a, c], [b, d], 'b')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 2), (4, 1), (4, 0)]\n",
    "plot_bezier2(Ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bezier2(cycle(8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. L'algorithme de De Casteljau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons dans cette section discuter d'un algorithme qui permet le calcul des courbes de Bézier sans passer par des coefficients du binôme. Il s'agit de *l'algorithme de De Casteljau*. Cet algorithme a été développé à la fin des années 1950 par Paul de Faget de Casteljau, alors qu'il travaillait chez Citroën. Il a été « redécouvert » par Pierre Bézier au début des années 1960. Ce dernier travaillait chez Renault (il faut savoir que les recherches de De Casteljau étant protégées par le secret industriel, Bézier ne pouvait pas connaître les résultats de De Casteljau)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $t\\in[0,1]$. Pour $k=0,\\ldots,n-1$, notons $Q_k$ le barycentre de $P_k$ et $P_{k+1}$ affectés des coefficients $1-t$ et $t$. On a donc\n",
    "\n",
    "$$Q_k = (1-t)P_k+tP_{k+1}$$\n",
    "\n",
    "Considérons maintenant la courbe de Bézier associée aux points de contrôle $Q_0,\\ldots,Q_{n-1}$. On a le résultat suivant.\n",
    "\n",
    "**Proposition** $B_{Q_0,\\ldots,Q_{n-1}}(t)=B_{P_0,\\ldots,P_{n}}(t)$.\n",
    "\n",
    "**Démonstration.** On a\n",
    "\n",
    "$$B_{Q_0,\\ldots,Q_{n-1}}(t)=\\sum_{k=0}^{n-1}\\binom {n-1}k t^k(1-t)^{n-1-k}((1-t)P_k+tP_{k+1})=\\sum_{k=0}^{n-1}\\binom {n-1}k t^{k}(1-t)^{n-k}P_k+\\sum_{k=0}^{n-1}\\binom {n-1}k t^{k+1}(1-t)^{n-1-k}P_{k+1}$$\n",
    "\n",
    "Dans la deuxième somme (notons-la $S_2$), faisons le changement d'indice $k'=k+1$. Il vient\n",
    "\n",
    "$$S_2=\\sum_{k=1}^{n}\\binom {n-1}{k-1} t^{k}(1-t)^{n-k}P_k$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$B_{Q_0,\\ldots,Q_{n-1}}(t)=(1-t)^nP_0+t^nP_n+\\sum_{k=1}^{n-1}\\left(\\binom {n-1}k +\\binom{n-1}{k-1}\\right)t^{k}(1-t)^{n-k}P_k$$\n",
    "\n",
    "Remarquons que\n",
    "\n",
    "$$\\binom {n-1}k +\\binom{n-1}{k-1}=\\binom n k$$\n",
    "\n",
    "et donc\n",
    "\n",
    "$$B_{Q_0,\\ldots,Q_{n-1}}(t)=(1-t)^nP_0+t^nP_n+\\sum_{k=1}^{n-1}\\binom {n}kt^{k}(1-t)^{n-k}P_k$$\n",
    "\n",
    "On reconnaît dans les deux termes isolés du membre de droite deux quantités qui peuvent être intégrées à la somme. Finalement :\n",
    "\n",
    "$$B_{Q_0,\\ldots,Q_{n-1}}(t)=\\sum_{k=0}^{n}\\binom {n}kt^{k}(1-t)^{n-k}P_k=B_{P_0,\\ldots,P_{n}}(t)\\ \\diamond$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Itération"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quel est l'intérêt du calcul que nous venons de faire ? La courbe de Bézier $B_{Q_0,\\ldots,Q_{n-1}}$ est associée à $n$ points de contrôle, alors que la courbe $B_{P_0,\\ldots,P_{n}}$ est associée à $n+1$ points. Nous pouvons recommencer en prenant les barycentres de deux points $Q_k$ successifs pour fabriquer une courbe de Bézier $B_{R_0,\\ldots,R_{n-2}}$ associée à $n-1$ points de contrôle et telle que \n",
    "\n",
    "$$B_{R_0,\\ldots,R_{n-2}}(t)=B_{Q_0,\\ldots,Q_{n-1}}(t)=B_{P_0,\\ldots,P_{n}}(t)$$\n",
    "\n",
    "Et ainsi de suite. Après $n$ itérations nous obtenons une courbe de Bézier $B_S$ associée à un unique point de contrôle $S$. On a alors\n",
    "\n",
    "$$B_{P_0,\\ldots,P_{n}}(t)=B_S(t)=S$$\n",
    "\n",
    "Ainsi, par des calculs de barycentres successifs de deux points, on est capables de déterminer $B_{P_0,\\ldots,P_{n}}(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `casteljau` prend en paramètres une liste `Ps` de points du plan et un réel $t\\in[0,1]$. Elle renvoie $B(t)$. Elle appelle une fonction `bary` dont nous parlerons juste après."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def casteljau(Ps, t):\n",
    "    n = len(Ps) - 1\n",
    "    for i in range(n):\n",
    "        Ps = bary(Ps, t)\n",
    "    return Ps[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `bary` prend en paramètres une liste `Ps` de points du plan et un réel $t\\in[0,1]$. Si `Ps`$=[P_0,\\ldots,P_n]$, la fonction renvoie $[Q_0,\\ldots,Q_{n-1}]$ où $Q_k$ est le barycentre de $P_k$ et $P_{k+1}$ affectés des coefficients $1-t$ et $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bary(Ps, t):\n",
    "    n = len(Ps)\n",
    "    Qs = []\n",
    "    for k in range(n - 1):\n",
    "        a, b = Ps[k]\n",
    "        c, d = Ps[k + 1]\n",
    "        Qs.append([(1 - t) * a + t * c, (1 - t) * b + t * d])\n",
    "    return Qs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [[0, 0], [1, 2], [2, 2], [4, 1], [4, 0]]\n",
    "print(casteljau(Ps, 0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un petit dessin vaut mieux qu'un long discours ... La fonction `demo_casteljau` prend en paramètres une liste `Ps` de points de contrôle et un réel $t\\in[0,1]$. Elle affiche la construction de $B(t)$ par l'algorithme de Casteljau."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def demo_casteljau(Ps, t):\n",
    "    n = len(Ps) - 1\n",
    "    ts = subdi(0, 1, 500)\n",
    "    Qs = [bezier0(Ps, t) for t in ts]\n",
    "    xs = [x for (x, y) in Qs]\n",
    "    ys = [y for (x, y) in Qs]\n",
    "    plt.plot(xs, ys, 'k')\n",
    "    for (x, y) in Ps:\n",
    "        plt.plot([x], [y], 'or')\n",
    "    for i in range(n):\n",
    "        xs = [x for (x, y) in Ps]\n",
    "        ys = [y for (x, y) in Ps]\n",
    "        plt.plot(xs, ys, '--k')\n",
    "        plt.plot(xs, ys, 'ok')\n",
    "        Ps = bary(Ps, t)\n",
    "    x, y = casteljau(Ps, t)\n",
    "    plt.plot([x], [y], 'or')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 2), (4, 1), (4, 0)]\n",
    "demo_casteljau(Ps, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [[0, 0], [1, 2], [2, 2], [4, 1], [4, 0]]\n",
    "demo_casteljau(cycle(5), 0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remarque.** Sur les dessins ci-dessus, il semble que la tangente à la courbe de Bézier au point rouge soit l'un des segments tracés. C'est effectivement le cas, nous y reviendrons plus loin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Comment utiliser l'algorithme de De Casteljau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Utilisé naïvement, l'algorithme de De Casteljau calcule $B(t)$ en temps $O(n^2)$. Il n'y a donc aucun gain de temps par rapport à l'algorithme initial utilisant la formule. Notons tout de même que nous avons réussi à éliminer les problèmes liés aux erreurs cumulées. Il reste maintenant à voir comment on utilise précisément cet algorithme. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mettons en place quelques notations. Notons tout d'abord $P_0^0,\\ldots,P_n^0$ les points de contrôle.\n",
    "\n",
    "Par un premier calcul de barycentres, on obtient les points $P_0^1,\\ldots,P_{n-1}^1$ où $P_i^1$ est le barycentre de $P_i^0$ et $P_{i+1}^0$ affectés des coefficients $1-t$ et $t$.\n",
    "\n",
    "Un second calcul de barycentres fournit les points $P_0^2,\\ldots,P_{n-2}^2$ où $P_i^2$ est le barycentre de $P_i^1$ et $P_{i+1}^1$ affectés des coefficients $1-t$ et $t$.\n",
    "\n",
    "Etc. On obtient ainsi une famille $(P_i^k)$ de points du plan, où $0\\le k\\le n$ et $0\\le i\\le n-k$. Le point $B(t)$ est le point $P_0^n$. Nous avons le résultat suivant.\n",
    "\n",
    "**Proposition.** La courbe de Bézier de points de contrôle $P_0,\\ldots,P_n$ est la réunion de deux courbes de Bézier :\n",
    "\n",
    "- La courbe de points de contrôle $P_0^0,P_0^1,\\ldots,P_0^n$.\n",
    "- La courbe de points de contrôle $P_0^n,P_1^{n-1},\\ldots,P_n^0$.\n",
    "\n",
    "Pour prouver cette proposition, montrons d'abord le lemme suivant.\n",
    "\n",
    "**Lemme.** Pour tout $k\\in[\\![0,n]\\!]$, pour tout $i\\in[\\![0,n-k]\\!]$,\n",
    "\n",
    "$$P_i^k=\\sum_{j=0}^kB_j^k(t)P_{i+j}^0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Démonstration du lemme.** Montrons la propriété par récurrence sur $k$.\n",
    "\n",
    "Pour $k=0$, c'est immédiat. Soit $k\\in[\\![0,n-1]\\!]$. Supposons que pour tout $i\\in[\\![0,n-k]\\!]$,\n",
    "\n",
    "$$P_i^k=\\sum_{j=0}^kB_j^k(t)P_{i+j}^0$$\n",
    "\n",
    "Soit $i\\in[\\![0,n-k-1]\\!]$. On a\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "P_{i}^{k+1}&=&(1-t)P_i^k+tP_{i+1}^k\\\\\n",
    "&=&(1-t)\\sum_{j=0}^kB_j^k(t)P_{i+j}^0+t\\sum_{j=0}^kB_{j}^k(t)P_{i+j+1}^0\\\\\n",
    "&=&(1-t)\\sum_{j=0}^kB_j^k(t)P_{i+j}^0+t\\sum_{j=1}^{k+1}B_{j-1}^k(t)P_{i+j}^0\\\\\n",
    "&=&(1-t)B_0^k(t)P_i^0+tB_k^k(t)P_{i+k+1}^0+\\sum_{j=1}^k((1-t)B_j^k(t)+tB_{j-1}^k(t))P_{i+j}\n",
    "\\end{array}$$\n",
    "\n",
    "Remarquons que\n",
    "\n",
    "$$(1-t)B_j^k(t)+tB_{j-1}^k(t)=\\left(\\binom k j +\\binom k {j-1}\\right)t^{j}(1-t)^{k+1-j}=\\binom {k+1} jt^{j}(1-t)^{k+1-j}$$\n",
    "\n",
    "De plus,\n",
    "\n",
    "$$(1-t)B_0^k(t)=(1-t)^{k+1}=B_0^{k+1}(t)$$\n",
    "\n",
    "et\n",
    "\n",
    "$$tB_k^k(t)=t^{k+1}=B_{k+1}^{k+1}(t)$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$P_{i}^{k+1}=\\sum_{j=0}^{k+1}B_j^{k+1}(t)P_{i+j}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Démonstration de la proposition.** On a pour tout $k\\in[\\![0,n]\\!]$,\n",
    "\n",
    "$$P_0^k=\\sum_{j=0}^kB_j^k(t)P_{j}^0$$\n",
    "\n",
    "Soit $u\\in[0,1]$. On a\n",
    "\n",
    "$$\\sum_{k=0}^nB_k^n(u)P_0^k=\\sum_{k=0}^nB_k^n(u)\\sum_{j=0}^kB_j^k(t)P_{j}^0=\\sum_{j=0}^nS_jP_{j}^0$$\n",
    "\n",
    "où\n",
    "\n",
    "$$S_j=\\sum_{k=j}^nB_j^k(t)B_k^n(u)$$\n",
    "\n",
    "Que vaut $S_j$ ? On a\n",
    "\n",
    "$$S_j=\\sum_{k=j}^nB_j^k(t)B_k^n(u)=\\sum_{k=j}^n\\binom k jt^j(1-t)^{k-j}\\binom n ku^k(1-u)^{n-k}$$\n",
    "\n",
    "Remarquons que\n",
    "\n",
    "$$\\binom j k \\binom k n=\\binom n j\\binom{n-j}{k-j}$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "S_j&=&\\binom n j\\sum_{k=j}^n\\binom{n-j}{k-j}t^j(1-t)^{k-j}u^k(1-u)^{n-k}\\\\\n",
    "&=&\\binom n j\\sum_{k=0}^{n-j}\\binom{n-j}{k-j}t^j(1-t)^ku^{k+1}(1-u)^{n-k-j}\\\\\n",
    "&=&\\binom n j(tu)^j\\sum_{k=0}^{n-j}\\binom{n-j}k(u(1-t))^k(1-u)^{n-j-k}\\\\\n",
    "&=&\\binom n j(tu)^j(1-u+u(1-t))^{n-j}\\\\\n",
    "&=&\\binom n j(tu)^j(1-tu)^{n-j}\\\\\n",
    "&=&B_j^n(tu)\n",
    "\\end{array}$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$\\sum_{k=0}^nB_k^n(u)P_0^k=\\sum_{j=0}^nB_j^n(tu)P_{j}^0=B(tu)$$\n",
    "\n",
    "Pour $u=0$, $tu=0$. Pour $u=1$, $tu=t$. On obtient donc la partie de la courbe de $B$ pour les valeurs du paramètre entre 0 et $t$.\n",
    "\n",
    "où $B$ est la fonction de Bézier associée aux points de contrôle $P_0,\\ldots,P_n$.\n",
    "\n",
    "On montre de même que\n",
    "\n",
    "$$\\sum_{k=0}^nB_k^n(u)P_{n-k}^k=B(1-(1-t)u)$$\n",
    "\n",
    "Pour $u=0$, $1-(1-t)u=1$. Pour $u=1$, $1-(1-t)u=t$. On obtient donc la partie de la courbe de $B$ pour les valeurs du paramètre entre $t$ et $1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remarque.** Nous pouvons maintenant confirmer ce que nous avions remarqué un peu plus haut au sujet des tangentes à la courbe de Bézier. Le point $B(t)$ est l'extrémité de la courbe de Bézier associée aux points de contrôle $P_0^0$, $P_0^1$, ... $P_0^n$. La tangente à gauche à la courbe en ce point est donc dirigée par le vecteur $\\overrightarrow{P_0^{n-1}P_0^n}$. De même, la tangente à droite est dirigée par le vecteur $\\overrightarrow{P_{n-1}^1P_n^0}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour obtenir la courbe de Bézier complète, il suffit d'appeler récursivement notre algorithme sur deux « moitiés » de courbe. L'algorithme s'arrête lorsqu'un certain critère est vérifié. Dans les applications pratiques, on travaille avec des points à coordonnées entières représentant des pixels et on prend $t=\\frac 1 2$. Ainsi, tous les calculs peuvent être réalisés avec des entiers. Dans ce cas, un bon critère d'arrêt est que les points de contrôle soient à des distances inférieures à 1. En effet, comme la courbe de Bézier est contenue dans l'enveloppe convexe de l'ensemble des points de contrôle, elle est alors incluse dans un pixel.\n",
    "\n",
    "Nous travaillons quant à nous avec des points dont les coordonnées sont des flottants. Pour simplifier, dans la fonction `bezier` ci-dessous, le critère choisi est le nombre $r$ d'appels récursifs, passé en paramètre à la fonction. Lorsque $r=0$, la fonction renvoie les deux points de contrôle extrêmes : ces deux points, nous le savons, font partie de la courbe de Bézier. Cette fonction prend également un paramètre optionnel $t$, fixé par défaut à $\\frac 1 2$. Il s'agit du paramètre utilisé pour les calculs de barycentres.\n",
    "\n",
    "Un appel à `bezier(Ps, r)` renvoie ainsi une liste de $2^r+1$ points. Dans la pratique, $r=9$ suffit, puisqu'on obtient environ 500 points (c'est ce nombre de points que nous prenons depuis le début du notebook pour tracer les courbes). Si la courbe de Bézier est vraiment très compliquée, on peut être amenés à prendre une valeur de $r$ plus grande."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bezier(Ps, r, t=0.5):\n",
    "    if r == 0: return [Ps[0], Ps[-1]]\n",
    "    n = len(Ps) - 1\n",
    "    Qs1 = [Ps[0]]\n",
    "    Qs2 = [Ps[-1]]\n",
    "    for k in range(n):\n",
    "        Ps = bary(Ps, t)\n",
    "        Qs1.append(Ps[0])\n",
    "        Qs2.append(Ps[-1])\n",
    "    Qs2.reverse()\n",
    "    Ps1 = bezier(Qs1, r - 1)\n",
    "    Ps2 = bezier(Qs2, r - 1)\n",
    "    return Ps1 + Ps2[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 3), (4, 1), (5, 3)]\n",
    "len(bezier(Ps, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `plot_bezier` affiche la courbe de Bézier avec notre nouvel algorithme. Elle prend en paramètres la liste des points de contrôle et le nombre $r$ d'appels récursifs à effectuer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bezier(Ps, r):\n",
    "    Qs = bezier(Ps, r)\n",
    "    xs = [x for (x, y) in Qs]\n",
    "    ys = [y for (x, y) in Qs]\n",
    "    plt.plot(xs, ys, 'k')\n",
    "    for (x, y) in Ps:\n",
    "        plt.plot([x], [y], 'or')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ps = [(0, 0), (1, 2), (2, 3), (4, 1), (5, 3)]\n",
    "plot_bezier(Ps, 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bezier(cycle(20), 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Les courbes de Bézier dans le monde « réel »"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une recherche rapide sur Internet permet de trouver facilement un grand nombre d'applications des courbes de Bézier. Voici, en quelques lignes, une porte d'entrée pour les recherches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Quelques applications concrètes des courbes de Bézier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les courbes de Bézier sont utilisées par un grand nombre d'applications logicielles. Pour n'en citer que quelques unes :\n",
    "\n",
    "- Le logiciel [Blender](https://www.blender.org)\n",
    "- Le format graphique [SVG](https://www.w3.org/Graphics/SVG/)\n",
    "- Le programme de manipulation d'images [Gimp](https://www.gimp.org/)\n",
    "- Les polices de caractères TrueType\n",
    "- Les polices de caractères utilisées par le système de composition de documents $\\LaTeX$ : le langage **MetaFont**, qui est un langage utilisé pour créer des polices de caractères vectorielles, utilise des courbes de Bézier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 B-splines, NURBS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Des généralisations des courbes de Bézier existent. Citons en particulier les B-splines, et une généralisation des B-splines, les NURBS (non uniform rational B-Splines)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Surfaces de Bézier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les *surfaces de Bézier* sont une généralisation à deux dimensions des courbes de Bézier. Celles-ci sont en particulier utilisées pour le rendu d'objets 3D. Plutôt que de décomposer l'objet à représenter en triangles, on le décompose en surfaces de Bézier. Cette technique a ses avantages, et aussi quelques inconvénients que nous ne détaillerons pas ici."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
