Courbes et surfaces de Bézier

Avant-Propos:

Voici une introduction aux courbes et surfaces de bézier que j'ai écrit pendant mon temps libre pour un forum de programmation. Il est assez ancien et pas vraiment parfait mais je ne le modifierai plus vu que je suis passé à autre chose depuis.

Vous pouvez cependant envoyez-vos commentaires à email.

Contenu:

1- Un peu de théorie: les maths derrière les courbes et surfaces de Bézier.

2- Implantation: je discute des algos d'implantation pour le rendu.

3- Code: un exemple d'implantation d'une classe tri-patch en C++/OpenGL

4- Conclusion

Un peu de théorie:

Intro

Avis aux gens qui sont allergiques aux maths, cette partie en est pleine donc si vous ne voulez pas comprendre le pourquoi du comment des choix d'implantation passez tout de suite aux articles suivants, même si cela perd beaucoup de son intérêt. (Niveau supposé: études maths supérieures, lycéens passez votre chemin..)

Des surfaces courbes en 3D en fait c'est une imposture puisque un écran d'ordinateur et un hardware 3D n'est capable que de représenter des entités discrètes que sont des pixels et des triangles. Il reste que si l'on cherche à reproduire un monde réaliste il y a un intermédiaire entre symboliser une sphère avec 4 triangles (un tétraèdre) et une vraie sphère. Ce problème n'est pas récent puisque la solution qui est proposée ici date d'une époque où les ordinateurs balbutaient et n'étaient certainement pas encore capables d'afficher du graphisme, ne parlons pas de 3D temps réel. Donc Monsieur Bézier (un français) cherchait un moyen simple de représenter une courbe, qui soit à la fois manipulable par un non mathématicien (par exemple dans le dessin de carosseries de voitures), et que l'on puisse construire à la règle et au crayon. La solution ce sont donc les courbes appelées "de Bézier" et qui ont révolutionné bien d'autres domaines que le dessin automobile (typographie, jeux vidéo en 3D..).

Les courbes de Bézier

Les courbes de Bézier ce sont en fait de simples courbes polynomiales. Ce n'est pas un "nouveau type" de courbe mais simplement un moyen de représenter des courbes connues.

Courbe polynomiale de degré 2: X(t) = t².A + t.B + C avec t réel;

Courbe polynomiale de degré 3: X(t) = t³.A + t².B + t.C + D avec t réel;

Courbe polynomiale de degré n: X(t) = ∑ti.Ai avec t réel;

Les courbes polynomiales de degré 1 n'ont qu'un intérêt très limité puisque ce sont .. des droites!

Nota: Ces fonctions représentent un certain type de courbe. Ce ne sont pas les seules qui existent, loin de là! Par exemple, on ne peut pas obtenir de cercle "exact" à partir de courbes polynomiales, même en les mettant bout à bout. Ce qu'on peut espérer représenter ici ce sont les courbes polynomiales par morceau. Dans les faits, on ne fait pas forcément la distinction, puisqu'on peut approximer suffisamment tout type de courbe avec une courbe polynomiale par morceau. N'oublions pas non plus qu'on travaille dans un domaine (le graphisme temps réel) où l'approximation est de mise, peu importe le moyen pourvu qu'on ait le résultat qu'on cherche :).

On a des moyens de les construire: on calcule la valeur de X, pour différentes valeur de t et on relie les points. Ce n'est pas très intuitif: quel est le lien entre les coéfficients du polynome et la courbe que l'on aura à l'écran? Les coéfficients vectoriels n'ont pas de signification géométrique bien définie. De plus, impossible de tracer ce type de courbe sans avoir un calculateur pour calculer les coordonnées des points. On est loin de la facilité évoquée ci-dessus.

La solution viendra d'une réécriture des fonctions précédentes, à l'aide de fonctions intermédiaires, les fonctions de Bernstein. Cela revient à faire un changement de base dans l'espace des fonctions polynomiales de degré n.

Définition des fonctions de Bernstein:
Les Bi,n(t) = Ci,n.ti.(1-t)(n-i), i variant entre 0 et n, constituent une base de l'espace des fonctions polynomiales de degré n (Ci,n est le coéfficient bien connu: Ci,n = n!/(i!(n-i)!) ).

Rappel: Qu'est-ce qu'une base? Cela veut dire qu'on peut obtenir toute fonction polynomiale de degré n à partir d'une combinaison linéaire de fonctions de Bernstein et que cette décomposition est unique.

(Exo: Si ca vous amuse vous pouvez prouver que c'est une base..).

Exprimée dans notre nouvelle base, nos fonctions deviennent:

Courbe polynomiale de degré 2 : P(t) = B2,2(t).P2 + B1,2(t).P1 + B0,1(t).P0 avec t réel;

Courbe polynomiale de degré n: P(t) = ∑Bi,n(t).Pi avec t réel;

Hum.. ca nous avance bien. On a remplacé une écriture compliquée par une écriture encore plus compliquée... En fait pas tant que ça. Cette écriture a bien des avantages que l'on va découvrir ici:


Courbe de Bézier avec 4 points de contrôle

Courbe de Bézier inscrite dans l'enveloppe convexe définie par ses points de contrôle

On raccorde simplement deux courbes de Bézier pour en créer une plus complexe

Approximation de la courbe par appels récursifs

Calculer tous ces milieux permet de subdiviser le problème de l'approximation de la courbe en un sous-problème d'approximation de deux sous-courbes (divide and conquer). Au fur et à mesure des récursions on a des points de contrôle qui sont de plus en plus proches de la courbe. Cela donne un algorithme de tracé, simplement à partir des points de contrôles, sans passer par la formule et qui peuvent être faits à la règle et au compas! Evidemment il est toujours possible de calculer les points en passant par la formule, et dans ce cas les points de controle sont juste une aide (non négligeable) au dessinateur.

Toutes ces propriétés ne sont intéressantes que sur l'intervalle [0,1]. En dehors de ce domaine, les informations données par les points de contrôle sont inintéressantes. C'est pour ça que quand on parle de courbes de Bézier, on parle en fait de fragments de courbes lorsque t varie entre 0 et 1.

Les surfaces de Bézier

L'idée des surfaces de Bézier est strictement la même que pour les courbes de Bézier, la différence c'est qu'on a un degré de liberté supplémentaire.

P(s,t) = ∑∑Bi,n(s).Bj,m(t).Pi,j, s et t évoluant sur R, l'ensemble des réels. (On se restreindra à s et t variant sur l'intervalle [0,1])

Les Pi,j définissent ce que l'on appelle une grille de contrôle (on les représente souvent sous forme d'une grille). L'allure locale de la surface est toujours déterminée par un point de contrôle.


Surface de Bézier définie par une grille à 16 points de contrôle

Ce qu'on remarque immédiatement c'est qu'à t=t0 constant on a simplement la définition d'une courbe de Bézier, ce qui donne des indices pour la construction de notre surface. Le point de controle Qi0(t0) de cette courbe est ∑Bj,m(t0).Pi0,j, c'est à dire le point correspondant à la valeur t sur la courbe de Bézier définie par les Pi0,j. En particulier pour t=0, les points de contrôle sont simplement les Pi0,0. Et c'est la même chose sur tous les bords de la grille: les bords de la surface sont des courbes de Bézier dont les points de contrôle sont les points de controle du bord correspondant de la grille. L'implication est la même que plus haut. C'est à dire que si l'on se donne deux grilles de contrôle avec des points en commun sur un de leur bord, alors ces deux surfaces seront parfaitement collées. Par contre ce qui est moins évident c'est l'influence des autres points de contrôle sur la qualité de ce raccord en terme de surfaces tangentes: on risque d'avoir des angles vifs sur les bords sauf dans certains cas particuliers.

Revenons sur ce problème de surfaces tangentes. Celles-ci sont définies par un simple vecteur normal à la surface, qui est égal au produit vectoriel de la dérivée selon s et selon t de la définition de la surface. (dans le cas non dégénéré où ces deux dérivées sont non nulles). Il n'y a pas de formules magiques pour interpoler ces dérivées sur toute la surface, il parait difficile d'échapper au calcul de dérivées. Cela a une importance si l'on a besoin de calculer les normales en chaque point de la surface (calcul d'éclairage, environment map etc..).

La surface tangente n'a de sens que sur un des points de la surface, donc elle n'a pas de sens en tout point de contrôle. Si l'on a besoin de spécifier tout de même un vecteur normal pour chaque point de la grille de contrôle, même pour les points de contrôle qui ne sont pas sur la surface, on peut imaginer faire une interpolation de ce vecteur normal à partir des 4 coins de la grille. Si l'on a suffisamment subdivisé notre surface et que l'on est très proche de la surface exacte, alors cette inexactitude est probablement négligeable d'un point de vue visuel.

Implantation:

Je n'ai pas trop envie de traiter le cas général, j'espère que quand vous voudrez implanter des surfaces ou des courbes plus complexes, les éléments que j'ai détaillés dans cet article vous suffiront. Ce qui suit concerne donc les surfaces et les courbes de degré 2, aussi appelées quadratiques (bicoze le terme de plus haut degré est un carré).

Courbes de Bezier:

Courbe polynomiale de degré 2 : P(t) = t².P2 + 2.t.(1-t).P1 + (1-t)².P0 avec t évoluant sur [0,1];

Pour tracer la courbe de Bézier en fonction de ses points de controles, on peut utiliser la méthode Divide and Conquer (diviser pour régner) ou la méthode plus directe de calcul de la courbe à des points prédéfinis. Dans un cas, on fournit au programme les points de contrôle et le degré de récursion, dans l'autre cas, on fournit les points de contrôle et le nombre de points intermédiaires que l'on veut sur la courbe.

Première solution: on a une grille de coefficients qui donnent les coordonnées de nouveaux points de contrôle en fonction des anciens. 

     P0  P1  P2
    ------------
P'0|  1   0   0;
P'1| 1/2 1/2  0;
P'2| 1/4 2/4 1/4;
P'3|  0  1/2 1/2;
P'4|  0   0   1;

A partir de 3 points de contrôle on en calcule 5. Soit on a atteint le niveau de complexité souhaité et on trace ces points à l'écran, soit on continue le processus avec les 3 premiers points de contrôle et les 3 derniers points de contrôle. On peut y arriver avec une simple procédure récursive. On remarque que les coéfficients restent les mêmes et donc si on traite toujours la même famille de courbe, il est possible de coder ces coéfficients en dur.

Procédure Trace (P0,P1,P2) {
Soit:
  Calcule(P'0,P'1,P'2,P'3,P'4);
  Trace (P'0,P'1,P'2); //appel récursif à Trace
  Trace (P'2,P'3,P'4);
Soit: //si on a atteint le niveau désire 
  Relie (P0,P1,P2); // on relie les points passés en paramètre
}

Voilà c'est on ne peut plus simple, non?

Deuxième solution: on calcule une grille de coéfficients qui donnent les coordonnées de points de la courbe en fonction des points de contrôle:

     P0    P1    P2  
    -----------------
P'0|   1     0     0;
P'1|  9/16  3/8  1/16;
P'2|  1/4   1/2   1/4;
P'3|  1/16  3/8  9/16;
P'4|   0     0     1;

Là encore on a 5 nouveaux points en fonction des 3 points de contrôle mais leur signification est différente. Tout d'abord ce sont des points de la courbe et non pas de nouveaux points de contrôle. D'autre part avec la procédure récursive, on calcule un certain nombre de points intermédiaires qui n'apparaitront pas à l'écran alors qu'avec cette méthode on est certain de ne calculer que ce qui est nécessaire. Par contre si l'on désire faire varier le niveau de complexité et de décomposition, il faut calculer et stocker un peu plus de coéfficients. Un jeu de coéfficient par nouveau point.

Procedure Trace(P0,P1,P2) {
  Calcule(P'0,...,P'n);
  Relie(P'0,..P'n);
}

Voilà, ça reste encore très simple, n'est-ce pas?

Surfaces de Bezier:

L'optique est la même que pour les courbes. On exprime les nouveaux points en fonction des points de contrôle. La différence évidemment c'est qu'il sont beaucoup plus nombreux. 9 points de contrôle pour une surface de degré 2 sur 2. On peut se ramener au problème des courbes en effectuant d'abord un calcul selon s de nouveaux points de contrôle, puis un calcul selon t en introduisant les points calculés précédemment. On peut également faire un calcul direct mais cela impose de hardcoder ou de calculer à la volée toute une floppée de coéfficients.

Dans le programme qui suit j'ai décidé, de recourir à la méthode récursive: à chaque étape on a neuf points de contrôle et on doit en calculer 25. J'ai parlé juste avant de l'indépendance des deux problèmes verticaux et horizontaux; j'adopte donc une méthode qui n'est pas optimale en temps de calcul (ceux qui le désirent peuvent hardcoder les coéfficients pour gagner des multiplications/additions). A partir de 3 points de contrôle horizontaux, j'en calcule 5 nouveaux qui définissent évidemment la même courbe de Bézier horizontale. Au total j'ai donc 3*5 = 15 points de contrôle horizontaux qui peuvent également être vus comme des points de contrôle verticaux. Pour trois points de contrôle verticaux calculés précédemment, je peux calculer 5 nouveaux points de contrôle qui définissent une même courbe de Bézier verticale. J'ai donc au final 5*5 = 25 nouveaux points de contrôle qui représentent une subdivision de ma surface en 4 sous-surfaces de Bézier. Et je réapplique le processus à chaque sous surface jusqu'au niveau de subdivision désiré.

Il ne suffit généralement pas de calculer les coordonnées des points, mais aussi des choses utiles au tracé, comme la couleur, les coordonnées de texture ou le vecteur normal à la surface. Ces calculs ne sont pas totalement innocents. Un méthode ou une autre peuvent avoir des répercussions sur l'aspect de la surface.

La couleur, la solution la plus simple c'est de se choisir 4 coordonnées de couleur aux quatre coins et de décider que la couleur du point (s,t) sera une interpolation bilinéaire en s et en t des coordonnées de couleur des coins.

Les coordonnées de texture, on peut tout simplement considérer que les coordonnées s et t sont des coordonnées de texture. On peut aussi faire le choix de coordonnées de texture implicites: mappage de texture plan, sphérique, cylindrique.. Les coordonnées de texture seraient donc définies de manière "externe" à la surface. Se pose alors le problème de la cohérence des coordonnées de texture par rapport aux mouvements éventuels de la surface, mais ça dépasse le cadre strict des surfaces de Bézier.

Les normales, difficile d'échapper au calcul des normales en chaque point de la surface. Dans le cas du degré 2, on a une idée des tangentes grâce aux points de contrôle: la normale est dans la direction du produit vectoriel des dérivées selon s et selon t en un point de la surface. Or ces dérivées valent respectivement (P1-P0), (P3-P0) en (0,0), (P1-P2),(P5-P2) en (1,0) etc.. Pour ce qui est de la normale en un point de contrôle, cela n'a pas de sens malheureusement, on peut fournir une "pseudo-normale" en interpolant linéairement par exemple les 4 normales aux 4 coins de la grille.

Optimisation/Utilité:

Si l'on est dans le cas d'une application temps réel, il n'est pas raisonnable de refaire tous ces calculs à chaque affichage et donc on peut envisager de stocker les coordonnées des nouveaux points en mémoire et de réutiliser ces calculs au prochain affichage (si le degré de complexité n'a pas changé).

Quel est l'intérêt de passer par ces fonctions, si l'on doit tout de même les transformer et les stocker sous forme de polygones au final? Et bien, il faut peser le pour et le contre. Pour l'édition, le plus est indéniable, cela augmente la créativité des formes si l'artiste peut utiliser tout une classe de formes rondes dans ses créations. Il suffit de regarder les maps de Quake 3 pour s'en rendre compte. Par ailleurs, si ces formes sont utiles pour la création, c'est également une forme avancée de compression de données, plutôt que de stocker chaque point d'une surface courbe dans un fichier autant ne stocker que les points de contrôle qui ont servi à sa création. La décompression (ou tesselation par référence aux tesselles d'une mosaïque) est un processus relativement rapide. D'autant que si le programme doit s'adapter à des configurations plus ou moins puissantes, il poussera cette décompression plus ou moins loin, ce qui évite d'avoir à stocker différents niveaux de complexité pour chaque modèle 3D. Dans un futur proche, on peut même imaginer que le hardware soit capable de traiter cette décomposition en temps réel, on n'aurait qu'à lui fournir les points de contrôle et il ferait la tesselation comme, il y a peu encore, il faisait la rasterisation. Imaginons une machine avec une bonne puissance de calcul mais une mémoire faible (Playstation 2 pour ne pas la citer), dans ce cas on utiliserait la puissance disponible au mieux en faisant effectuer par son unité de calcul vectoriel programmable toutes les opérations de tesselation, tout en ne surchargeant pas la mémoire de données inutiles.

Des surfaces de Bézier en lancer de rayon:

Supposons que l'on cherche à afficher des surfaces de Bézier avec une méthode par lancer de rayons.

Première solution: On traite la surface de Bézier comme une simple surface polygonale. On applique donc l'algorithme de tesselation défini ci dessus. On a en prime une enveloppe convexe de la surface, fonction des points de contrôle, qui peut accélérer le calcul de collision surface/rayon. Pour ce qui est de l'affichage, on est ramené au cas des surfaces polygonales que je ne traiterai pas ici.

Deuxième solution: On veut conserver l'écriture implicite de la courbe, et pouvoir calculer de manière exacte le point d'intersection et la normale en ce point. Calcul de l'intersection surface/droite: on cherche u, s, t tels que X(s,t)=M(u) avec M(u) = O+u.V. On est ramené à un système de trois équations à trois inconnues, de degré au plus deux. On peut le résoudre exactement et on obtient s, t et u d'où les coordonnées du point d'intersection et la valeur de la normale en ce point.

Code:

Que serait un article de programmation sans un peu de code. Le programme que j'ai réalisé dans le cadre de cet article essaie d'implanter une classe tripatch de surfaces de Bézier définies par une grille 3 sur 3 de points de contrôles. Il vous donnera les bases nécessaires j'espère pour passer à des surfaces plus compliquées ou pour implanter des optimisations indispensables dans le monde temps réel. C'est un peu rapide mais c'est une simple illustration de l'article pas un programme complet.

Pour utiliser la classe vous devez récuperer les deux fichiers:

tripatch.h

tripatch.cpp

On commence par définir une classe vecteur3D, ainsi que quelques opérations de base sur ces vecteurs:

class vecteur3D {
public:
    GLfloat x,y,z;
    vecteur3D() {}
    vecteur3D(GLfloat _x, GLfloat _y, GLfloat _z ) :
        x(_x), y(_y), z(_z) {}
    const vecteur3D operator +  (const vecteur3D & _vect) const {
        return (vecteur3D(x+_vect.x, y+_vect.y, z+_vect.z));
    }
    const vecteur3D operator -  (const vecteur3D & _vect) const {
        return (vecteur3D(x-_vect.x, y-_vect.y, z-_vect.z));
    }
    const vecteur3D operator /  (GLfloat scal) const {
        return (vecteur3D(x/scal, y/scal, z/scal));
    } 
    GLfloat normecarree() const {
        return (x*x+y*y+z*z);
    }
    const vecteur3D operator ^ (const vecteur3D & _vect) const {
        return (vecteur3D(y*_vect.z-z*_vect.y, 
                          z*_vect.x - x*_vect.z, 
                          x*_vect.y - y*_vect.z));
    }
};

inline const vecteur3D operator * (GLfloat scal, const vecteur3D& _vect) {
    return (vecteur3D(scal * _vect.x, scal * _vect.y, scal * _vect.z));
}

inline const vecteur3D normalize (const vecteur3D& _vect) {
    GLfloat _norm = _vect.normecarree();
    if (_norm==0)
        return _vect;
    else
        return (_vect/(sqrt(_norm)));
}

Cette implémentation est vraiment basique, il faudra penser à l'étoffer un peu ou à adapter mon code à votre propre version de classe vecteur3D.

On continue par implanter une classe vertex, avec ses attributs position, vecteur normal, couleur, coordonnées de texture..

class vertex {
public:
    vecteur3D vect;
    vecteur3D norm;
    vecteur3D rgb;
    vecteur3D uvw;
    vertex() {}
    vertex(vecteur3D _vect, vecteur3D _norm, vecteur3D _rgb, vecteur3D _uvw):
        vect(_vect), norm(_norm), rgb(_rgb), uvw(_uvw) {}
};

Enfin la classe tripatch qui représente la surface de Bézier à 9 points de contrôle:

class TriPatch {
private:
    //La grille de contrôle
    vertex grid[9]; 
    GLint type;
    GLint nsize; 
    vector<vertex> vertexArray; 
    // ce conteneur doit contenir les données temporaires
    // issues de la tesselation
    GLint currentLevel;
    void _compute(GLint _level);
    void _compute2(vertex **_array);
public:
    TriPatch():currentLevel(-1), type(0) {};
    TriPatch(GLint _type):currentLevel(-1), type(_type) {}
    // renvoie une référence pour manipuler le vertex de la grille
    vertex & getVertex(GLint i){ return grid[i]; }
    
    // accesseurs du type
    inline void setType(GLint _type){type=_type; currentLevel=-1;}
    inline GLint getType() { return type; }
    // On trace au niveau désiré et on recalcule si besoin
    void draw(GLint _level=0); 
};

Résultat:

Voici ce que donne un objet tripatch instancié dans un simple framework GLUT:


Tripatch après une subdivision 2.

Résumé, conclusion:

On a donc vu qu'on dispose avec les surfaces de Bézier de surfaces courbes faciles à appréhender et que l'on peut approximer de manière rapide (linéaire en le nombre de points..). Les intérêts de telles surfaces sont nombreux: stockage des courbes de manière compacte, subdivision à la demande en fonction des besoins et du niveau de détail désiré. Je crois même que les hardware de prochaine génération comme le GeForce3 sauront "tesseler" ces courbes sans aide du processeur principal, ce qui augmente d'autant leur intérêt.

Quelles sont les limites de cette méthode?

Nous avons abordé un peu plus haut le problème lié aux raccords entre surfaces, problème qui est encore plus fort en cas d'animation des meshes puisque un raccord dans une position ne sera plus bon à la position suivante. Il y a aussi un problème lié à la topologie: les patchs qu'ils soient tri- quadri- ou n présentent tous la même topologie, celui d'un "drap" ou "nappe"; vous voyez l'idée. L'idée de rapprocher tout objet qu'il soit sphère, donut, ou monstre géant de Quake à un simple mouchoir que l'on tord dans tous les sens est un peu saugrenue; même si pour y arriver on peut mettre les pièces de mouchoir bout à bout, le fait est que les coutures sont apparentes lorsqu'on commence à bouger!

Certains de ces problèmes sont résolus par une surclasse de ce type de surfaces, aussi appelées subdivision surface, mais c'est un autre sujet. (peut-etre un prochain article?)

Références:

A titre d'info, les articles suivants parmi d'autres m'ont aidé:

Curved surface using Bézier Patch:
http://www.gamasutra.com/features/19990611/bezier_01.htm

Unraveling Bézier spline:
http://www.gamedev.net/reference/articles/article888.asp

Implementing Curved Surface Geometry:
http://gamasutra.com/features/20000530/sharp_01.htm

Sinon la littérature sur ce genre de sujet est légion sur le net.

 

Partner websites : LEGREG | GRAPHICS | GRAPHISME | PHOTOGRAPHY | OUT OF MY MIND | ANIMATION MENTOR | GREEN LIVING | VOXEL | RAY TRACING