Algorithme de Casteljau

Considérons une courbe de Bézier définie par les points de contrôles ${\displaystyle P_{0},\dots ,P_{N}}$, où les ${\displaystyle P_{i}}$ sont des points de ${\displaystyle \mathbb {R} ^{m},m\geqslant 2}$. Ici on souhaite simplement calculer le point de paramètre ${\displaystyle t\in [0,1]}$.

Principe de l'algorithme de Casteljau

Comme on peut le voir sur l'image, en calculant les barycentres de paramètres ${\displaystyle \{t,1-t\}}$ des points de contrôle consécutifs de la courbe, puis les barycentres de même paramètres de ces barycentres et ainsi de suite itérativement, on définit de cette manière une suite de listes de points que l'on va indexer ${\displaystyle P_{i}^{j}}$, où ${\displaystyle P_{i}^{j}}$ est le barycentre de ${\displaystyle \{P_{i}^{j-1},P_{i+1}^{j-1}\}}$. La dernière liste ne contient en fait qu'un point, qui est le point de la courbe de paramètre ${\displaystyle t}$.

En pseudo-code, ceci donne:

 // Calcul des points intermédiaires
 Pour j de 1 à N faire
 |
 | Pour i de 0 à N-j faire
 | |
 | | T[i][j] = t*T[i+1][j-1] + (1-t)*T[i][j-1]
 | |
 |
 Afficher T[0][N] // Afficher (ou stocker) le point

Pour démontrer l'algorithme, il faut prouver par récurrence limitée que

${\displaystyle \forall i\in [0,n],\forall t\in [0,1],B(t)=\sum _{k=0}^{n}P_{k}^{0}B_{k}^{n}(t)=\sum _{k=0}^{n-i}P_{k}^{i}B_{k}^{n-i}(t)}$

Et d'appliquer la formule en ${\displaystyle i=n}$, ce qui nous donne le résultat directement.

La récurrence se démontre facilement en utilisant la propriété de construction des points ${\displaystyle P_{i}^{j}=(1-t)P_{i}^{j-1}+tP_{i+1}^{j-1}}$ pour séparer la somme en deux, puis en faisant un changement d'indice et en utilisant une propriété des polynômes de Bernstein ${\displaystyle B_{i}^{n}(t)=(1-t)B_{i}^{n-1}(t)+tB_{i-1}^{n-1}(t)}$ pour regrouper les deux sommes. Pour réintégrer les cas particuliers, il faut utiliser deux autres propriétés des polynômes de Bernstein: ${\displaystyle B_{n}^{n}(t)=t^{n}=tB_{n-1}^{n-1}}$ et ${\displaystyle B_{0}^{n}(t)=(1-t)^{n}=(1-t)B_{0}^{n-1}}$

L'exécution d'une étape de l'algorithme est quadratique en nombre de points de contrôle de la courbe (la boucle imbriquée donne lieu à ${\displaystyle {\mathcal {O}}(N^{2})}$ opérations).

On va ici appliquer l'algorithme précédent pour trouver le point de paramètre ${\displaystyle 1/2}$, et conserver les barycentres intermédiaires. En effet, la courbe de Bézier de points de contrôle ${\displaystyle P_{0}^{i},i\in [0,N]}$ est égale à la restriction de la courbe originale à ${\displaystyle t\leqslant 1/2}$, et la courbe de Bézier de points de contrôle ${\displaystyle P_{i}^{N-i},i\in [0,N]}$ est égale à la restriction de la courbe originale à ${\displaystyle t\geqslant 1/2}$.

On affiche ou mémorise le point de paramètre ${\displaystyle 1/2}$ (qui est ${\displaystyle P_{0}^{N}}$) et applique récursivement l'algorithme sur les deux courbes (de même nombre de points de contrôle que l'originale). On s'arrête quand un critère à choisir est vérifié (typiquement la distance entre les points inférieure à une limite).

Plutôt que le paramètre ${\displaystyle 1/2}$, on pourrait prendre un paramètre quelconque et l'algorithme fonctionnerait encore, mais le paramètre ${\displaystyle 1/2}$ est celui qui converge en moyenne le plus rapidement. Le paramètre ${\displaystyle 1/2}$ est aussi celui pour lequel les calculs sont les plus rapides quand l'on travaille en coordonnées entières : le calcul de chaque barycentre se fait par une addition et un décalage à droite pour chaque coordonnée, c'est-à-dire sans multiplication ni division.

• Initialisation: affecter le tableau des points de contrôle dans les ${\displaystyle P_{i}^{0}}$
• Voir que le critère d'arrêt n'est pas vérifié: il y a plusieurs possibilités dont :
  • Si on fait tous les calculs avec des nombres entiers, un choix peut être de s'arrêter lorsque tous les points sont confondus. (c'est-à-dire que la courbe n'est représentée que par un seul pixel).
  • Si le calcul n'est pas fait sur des nombres entiers, on peut s'arrêter quand les points sont distants d'une distance inférieure à une valeur choisie.
  • On effectue M itérations puis on relie les points obtenus (c'est-à-dire que l'on trace le polygone de Bézier), M étant déterminé empiriquement.
• Calcul des points intermédiaires de l'algorithme : il y a deux méthodes possibles qui donnent finalement les mêmes points :
  • Méthode constructive (en construisant une suite de milieux): On définit itérativement les points ${\displaystyle P_{0}^{1},\dots ,P_{N-1}^{1},P_{0}^{2},\dots ,P_{N-2}^{2},\dots ,P_{0}^{N-1},P_{1}^{N-1},P_{0}^{N}}$ tels que ${\displaystyle P_{i}^{j}}$ soit le milieu de ${\displaystyle [P_{i}^{j-1}P_{i+1}^{j-1}]}$
  • Méthode matricielle (en construisant directement les points comme barycentre, les coefficients étant donnés par les matrices de Casteljau) :

${\displaystyle {\begin{pmatrix}P_{0}^{0}\\\vdots \\P_{0}^{N}\end{pmatrix}}=D_{0}^{N}*{\begin{pmatrix}P_{0}^{0}\\\vdots \\P_{N}^{0}\end{pmatrix}}}$ et ${\displaystyle {\begin{pmatrix}P_{0}^{N}\\\vdots \\P_{N}^{0}\end{pmatrix}}=D_{1}^{N}*{\begin{pmatrix}P_{0}^{0}\\\vdots \\P_{N}^{0}\end{pmatrix}}}$

• Mémorisation: on mémorise le point ${\displaystyle P_{0}^{N}}$ qui est sur la courbe.
• Appel récursif: on appelle l'algorithme sur les deux courbes de Bézier intermédiaires définies par les points ${\displaystyle P_{0}^{i},i\in [0,n]}$ d'une part, ${\displaystyle P_{i}^{n-i},i\in [0,n]}$ d'autre part.

Voici un exemple d'implémentation de l'algorithme en pseudo-code avec pour critère d'arrêt l'égalité des points (on travaille donc sur des entiers) et la construction constructives pour calculer les points intermédiaires :

 Entrée : tableau T[0][0…N] des coordonnées des points de contrôle.

 Si T[0][0] = T[0][1] = … = T[0][N]	//si le critère d'arrêt est vérifié on s'arrête
 alors
 |
 | fin
 |
 Sinon  //pour dessiner
 |
 | // Calcul des points intermédiaires
 | Pour i de 1 à N faire
 | |
 | | Pour j de 0 à N-i faire
 | | |
 | | | T[i][j] = milieu de T[i-1][j] T[i-1][j+1]
 |
 | Afficher T[N][0] // Afficher (ou stocker) le point milieu
 |
 | // Construction des courbes restreintes
 | Pour i de 0 à N faire
 | |
 | | T'[i] = T[i][0]
 | | T"[i] = T[N-i][i]
 |
 | // Appel récursif
 | de_Casteljau(T')
 | de_Casteljau(T")

Si on prend comme critère d'arrêt un nombre d'appels récursifs constant, comme le nombre de points de contrôle est constant pendant les appels récursifs reste constant et qu'à chaque étape de récursion on double le nombre de courbes étudiées, la complexité de l'algorithme est en ${\displaystyle {\mathcal {O}}(N^{2}*2^{M})}$ avec ${\displaystyle M}$ le nombre de récursions (mais il est linéaire en nombre de points calculés car pour ${\displaystyle M}$ itérations il y a ${\displaystyle 2^{M}}$ points calculés).

• Claude Brezinski, Méthodes numériques de base, Techniques de l'ingénieur (lire en ligne), p. 7
• Jean-Louis Merrien, Analyse numérique avec MATLAB, Dunod, coll. « Sciences Sup / Exercices et problèmes », 2007 (lire en ligne), p. 89-90
• Alain Yger et Jacques-Arthur Weil, Mathématiques appliquées L3 : Cours complet avec 500 tests et exercices corrigés, Pearson Education France, 2009 (lire en ligne), p. 115