Algorisme de De Casteljau

En el camp matemàtic de l'anàlisi numèrica, l'algorisme de De Casteljau és un mètode recursiu per avaluar polinomis en forma de Bernstein o corbes de Bézier, que rep el nom del seu inventor Paul de Casteljau. L'algorisme de De Casteljau també es pot utilitzar per dividir una sola corba de Bézier en dues corbes de Bézier amb un valor de paràmetre arbitrari.

L'algorisme és numèricament estable^[1] en comparació amb l'avaluació directa de polinomis. La complexitat computacional d'aquest algorisme és ${\displaystyle O(dn^{2})}$, on d és el nombre de dimensions i n és el nombre de punts de control. Existeixen alternatives més ràpides.^[2][3]

Una corba de Bézier ${\displaystyle B}$ (de grau ${\displaystyle n}$, amb punts de control ${\displaystyle \beta _{0},\ldots ,\beta _{n}}$ ) es pot escriure en forma de Bernstein de la següent manera^[4]$${\displaystyle B(t)=\sum _{i=0}^{n}\beta _{i}b_{i,n}(t),}$$ on ${\displaystyle b}$ és un polinomi de base de Bernstein $${\displaystyle b_{i,n}(t)={n \choose i}(1-t)^{n-i}t^{i}.}$$ La corba en el punt ${\displaystyle t_{0}}$ es pot avaluar amb la relació de recurrència $${\displaystyle {\begin{aligned}\beta _{i}^{(0)}&:=\beta _{i},&&i=0,\ldots ,n\\\beta _{i}^{(j)}&:=\beta _{i}^{(j-1)}(1-t_{0})+\beta _{i+1}^{(j-1)}t_{0},&&i=0,\ldots ,n-j,\ \ j=1,\ldots ,n\end{aligned}}}$$Aleshores, l'avaluació de

${\displaystyle B}$ al punt ${\displaystyle t_{0}}$ es pot avaluar en ${\textstyle {\binom {n}{2}}}$ operacions. El resultat ${\displaystyle B(t_{0})}$ ve donat per $${\displaystyle B(t_{0})=\beta _{0}^{(n)}.}$$A més, la corba de Bézier ${\displaystyle B}$ es pot dividir en un punt ${\displaystyle t_{0}}$ en dues corbes amb els seus respectius punts de control: $${\displaystyle {\begin{aligned}&\beta _{0}^{(0)},\beta _{0}^{(1)},\ldots ,\beta _{0}^{(n)}\\[1ex]&\beta _{0}^{(n)},\beta _{1}^{(n-1)},\ldots ,\beta _{n}^{(0)}\end{aligned}}}$$

Interpretació geomètrica

La interpretació geomètrica de l'algorisme de De Casteljau és senzilla.

• Considereu una corba de Bézier amb punts de control ${\displaystyle P_{0},\dots ,P_{n}}$ Connectant els punts consecutius creem el polígon de control de la corba.
• Ara, subdivideix cada segment de línia d'aquest polígon amb la raó ${\displaystyle t:(1-t)}$ i connecteu els punts que obtingueu. D'aquesta manera arribeu al nou polígon amb un segment menys.
• Repetiu el procés fins que arribeu al punt únic; aquest és el punt de la corba corresponent al paràmetre ${\displaystyle t}$.

La imatge següent mostra aquest procés per a una corba de Bézier cúbica:

Cal tenir en compte que els punts intermedis que es van construir són, de fet, els punts de control per a dues noves corbes de Bézier, ambdues exactament coincidents amb l'antiga. Aquest algorisme no només avalua la corba a ${\displaystyle t}$, però divideix la corba en dues parts a ${\displaystyle t}$, i proporciona les equacions de les dues subcorbes en forma de Bézier.

La interpretació donada anteriorment és vàlida per a una corba de Bézier no racional. Per avaluar una corba de Bézier racional en ${\displaystyle \mathbf {R} ^{n}}$, podem projectar el punt a ${\displaystyle \mathbf {R} ^{n+1}}$ ; per exemple, una corba en tres dimensions pot tenir els seus punts de control ${\displaystyle \{(x_{i},y_{i},z_{i})\}}$ i pesos ${\displaystyle \{w_{i}\}}$ projectat als punts de control ponderats ${\displaystyle \{(w_{i}x_{i},w_{i}y_{i},w_{i}z_{i},w_{i})\}}$ L'algorisme continua com de costum, interpolant en ${\displaystyle \mathbf {R} ^{4}}$ Els punts quadridimensionals resultants es poden projectar de nou en tres espais amb una divisió de perspectiva.

En general, les operacions sobre una corba (o superfície) racional són equivalents a les operacions sobre una corba no racional en un espai projectiu. Aquesta representació com a "punts de control ponderats" i pesos sovint és convenient a l'hora d'avaluar corbes racionals.^[5]

En avaluar una corba de Bézier de grau n en un espai tridimensional amb n + 1 punts de control P _i $${\displaystyle \mathbf {B} (t)=\sum _{i=0}^{n}\mathbf {P} _{i}b_{i,n}(t),\ t\in [0,1]}$$ amb $${\displaystyle \mathbf {P} _{i}:={\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}},}$$ Dividim la corba de Bézier en tres equacions separades

$${\displaystyle {\begin{aligned}B_{1}(t)&=\sum _{i=0}^{n}x_{i}b_{i,n}(t),&t\in [0,1]\\[1ex]B_{2}(t)&=\sum _{i=0}^{n}y_{i}b_{i,n}(t),&t\in [0,1]\\[1ex]B_{3}(t)&=\sum _{i=0}^{n}z_{i}b_{i,n}(t),&t\in [0,1]\end{aligned}}}$$

que avaluem individualment mitjançant l'algorisme de De Casteljau.

Volem avaluar el polinomi de Bernstein de grau 2 amb els coeficients de Bernstein $${\displaystyle {\begin{aligned}\beta _{0}^{(0)}&=\beta _{0}\\[1ex]\beta _{1}^{(0)}&=\beta _{1}\\[1ex]\beta _{2}^{(0)}&=\beta _{2}\end{aligned}}}$$ en el punt t₀.

Comencem la recursió amb $${\displaystyle {\begin{aligned}\beta _{0}^{(1)}&&=&&\beta _{0}^{(0)}(1-t_{0})+\beta _{1}^{(0)}t_{0}&&=&&\beta _{0}(1-t_{0})+\beta _{1}t_{0}\\[1ex]\beta _{1}^{(1)}&&=&&\beta _{1}^{(0)}(1-t_{0})+\beta _{2}^{(0)}t_{0}&&=&&\beta _{1}(1-t_{0})+\beta _{2}t_{0}\end{aligned}}}$$ i amb la segona iteració la recursió s'atura amb $${\displaystyle {\begin{aligned}\beta _{0}^{(2)}&=\beta _{0}^{(1)}(1-t_{0})+\beta _{1}^{(1)}t_{0}\\\ &=\beta _{0}(1-t_{0})(1-t_{0})+\beta _{1}t_{0}(1-t_{0})+\beta _{1}(1-t_{0})t_{0}+\beta _{2}t_{0}t_{0}\\\ &=\beta _{0}(1-t_{0})^{2}+\beta _{1}2t_{0}(1-t_{0})+\beta _{2}t_{0}^{2}\end{aligned}}}$$ que és el polinomi de Bernstein esperat de grau 2.^[6]

Aquí teniu exemples d'implementacions de l'algorisme de De Casteljau en diversos llenguatges de programació.

deCasteljau :: Double -> [(Double, Double)] -> (Double, Double)
deCasteljau t [b] = b
deCasteljau t coefs = deCasteljau t reduced
 where
 reduced = zipWith (lerpP t) coefs (tail coefs)
 lerpP t (x0, y0) (x1, y1) = (lerp t x0 x1, lerp t y0 y1)
 lerp t a b = t * b + (1 - t) * a

def de_casteljau(t: float, coefs: list[float]) -> float:
 """De Casteljau's algorithm."""
 beta = coefs.copy() # values in this list are overridden
 n = len(beta)
 for j in range(1, n):
 for k in range(n - j):
 beta[k] = beta[k] * (1 - t) + beta[k + 1] * t
 return beta[0]

public double deCasteljau(double t, double[] coefficients) {
 double[] beta = coefficients;
 int n = beta.length;
 for (int i = 1; i < n; i++) {
 for (int j = 0; j < (n - i); j++) {
 beta[j] = beta[j] * (1 - t) + beta[j + 1] * t;
 }
 }
 return beta[0];
}