Rapport Algorithmique et complexité |
L’objectif du projet est de réaliser un algorithme de lancer de rayon permettant la visualisation de surfaces implicites de type f(x,y,z)=0.
Ce projet comporte donc deux notions ; d’une part, l’algorithme proprement dit, d’autre part, la résolution d’équation de surfaces implicites.
Dans un premier temps nous étudierons le concept de lancer de rayon, puis nous étudierons les surfaces algébriques et leur méthode de résolution.
Nous conclurons avec l’implémentation de notre moteur de rendu programmé en ocaml ainsi que l’analyseur syntaxique (ocamlyacc et lex) utilisé pour définir nos surfaces.
Historiquement, le lancer de rayon est la première méthode de rendu qui a permis d’approcher le réalisme photographique. Le terme de ’lancer de rayons’ à été introduit en 1968 par Arthur Appel.
Le lancer de rayon (ray tracing en anglais) est une technique de rendu en synthèse d’image simulant le parcours inverse de la lumière de la scène vers l’œil.
Cette technique simple reproduit les phénomènes physiques que sont la réflexion et la réfraction. Une mise en œuvre naïve du lancer de rayon ne peut rendre compte d’autres phénomènes optiques tels que les caustiques (taches lumineuses créées à l’aide d’une lentille convergente par exemple) et la dispersion lumineuse (la radiosité s’attaque à ce problème).
En revanche, contrairement à d’autres algorithmes de synthèse d’image, elle permet de définir mathématiquement les objets à représenter et non pas seulement par une multitude de facettes.
Le ray-tracing est entièrement fondé sur une vision particulaire.
Le principe de la vision est le suivant : Les sources de lumière (lampes, soleil...) émettent des rayons lumineux. Ces rayons viennent frapper les objets qui en réfléchissent une partie, en fonction de leurs propriétés physiques ("leur matériau"). Parmi tous les rayons lumineux, une partie va (directement ou après réflexion sur les objets) parvenir à l’œil et former une image sur la rétine.
Il serait illusoire de vouloir appliquer ce principe informatiquement. Une infime partie de la lumière émise parvient à l’oeil de l’utilisateur. La majorité de la lumière "se perd" dans la nature. La technique du lancer de rayon tente de reconstituer le parcours inverse de la lumière, en partant de l’oeil et en allant vers les sources.
Pour cela, on place une image virtuelle dans la scène devant l’observateur (c’est cette image qui s’affichera sur votre écran). Pour chacun des pixels qui constituent l’image, on lance un rayon partant du point de vue (l’observateur) passant par le pixel (son centre ou un autre point au hasard). La couleur de pixel traversé va être déterminée en suivant le cheminement du rayon lancé jusqu’aux sources lumineuses de la scène 3D.
Pour chaque rayon lancé de l’œil, on calcule le premier objet qu’il intersecte, ce qui implique un parcours de l’ensemble des faces de la scène.
Une fois le premier objet trouvé, il faut déterminer la couleur du pixel à afficher. Cette couleur dépend du matériau de la face, et bien sûr de l’éclairage. Aussi faut-il déterminer quelles sont les sources lumineuses qui éclairent le point. Pour cela, on relance des rayons issus du point considéré vers chacune des sources lumineuses de la scène. Si ce rayon parvient à la source sans intersecter d’autre objet, alors la source contribue à l’éclairage du point. Sinon, un objet masque la source, et un effet d’ombre est engendré naturellement.
La figure 1-2 est une vue de dessus du cheminement du rayon. Le rayon R1 provient de l’oeil, disons pour le pixel (x,y). Il heurte la face F1 au point p. Pour déterminer l’éclairage du pixel (x,y). Les rayons lancés de p montrent que la source S1 va contribuer à l’éclairage de P, mais que S2 est masquée par F2. La face F2 va porter son ombre sur F1.
Dans le cas où le rayon lancé depuis le point de vue n’intersecte aucun objet, le pixel (x,y) considéré prend une valeur "de fond", qui peut être soit une couleur fixe, soit la couleur du pixel (x,y) d’une image.
La méthode du lancer de rayons permet d’obtenir des résultats de très grande qualité et est capable de recréer des phénomènes physiques que les autres méthodes de rendu ne peuvent reproduire. Mais le prix à payer est un temps de calcul bien plus long qu’avec le tampon de profondeur (z-buffer), qui fait office de référence en terme de vitesse.
Développé depuis de nombreuses années, les algorithmes de lancer de rayons ont fait l’objet de nombreuses recherches qui ont abouti à des améliorations notables de la méthode.
La grande difficulté de cet algorithme est de mettre en place une méthode permettant de tester les intersections avec des objets de la scène.
Notre projet se référant aux surfaces implicite de type f(x,y,z)=0, il faut par conséquent résoudre les équations et trouver les solutions pour déterminer les points d’impact. Ceci sera détaillé dans la partie surface algébrique.
Une surface est dite algébrique lorsqu’elle possède une équation cartésienne polynomiale à coefficients réels.
Le degré (ou parfois l’ordre) d’une surface algébrique d’équation P(x,y,z)=0 est le degré du polynôme P, supposé sans facteur carré .
C’est aussi le nombre de points d’intersection, comptés avec les multiplicités, de la surface avec une droite quelconque, dans l’espace projectif complexe (et dans la pratique, le nombre maximum de points d’intersection de la surface avec une droite dans l’espace affine réel est un minorant du degré et donne très souvent le degré).
La notion de degré est projective, à savoir que toute homographie transforme une surface algébrique en une surface algébrique de même degré.
La surface est dite décomposable si elle est réunion de surfaces algébriques de degré plus petit et indécomposable dans le cas contraire, c’est-à-dire quand le polynôme P est irréductible dans R[X,Y,Z].
Les surfaces de degré 1 sont les plans, de degré 2 les quadriques, de degré 3 les surfaces cubiques, de degré 4 les surfaces quartiques.
Exemple de surface algébrique de degré 8 : l’oloïde.
Exemple de surface algébrique de degré 9 : la surface d’Enneper.
Toute surface rationnelle est algébrique, mais la réciproque est fausse.
Le problème majeur pour l’utilisation des surfaces algébriques dans le moteur de rendu en lancer de rayon est de définir les points d’intersections entre elles et les rayons projetés.
Le test d’intersection entre une surface et un rayon correspond à trouver les solution pour f(x,y,z)=0.où x,y,z sont les coordonnées du point d’intersection.
Les équations de ces surfaces étant en majeure partie de degré élevé, la recherche des racines pour f( )=0 devient particulièrement délicate sans une méthode appropriée.
Nous avons donc choisi d’utilisé la méthode de Casteljau qui utilise les polynômes dans la base de Bernstein. Nous définirons l’utilisation de cette méthode dans la partie résolution de ce rapport.
Un autre problème pour appliquer le lancer de rayon sur ce type de surface est de définir la normale au point d’intersection avec l’objet de manière à réémettre un rayon à partir du point d’intersection.
Pour ce type de surface, on calcule le gradient de l’équation de la surface et puis on génère le vecteur en appliquant les coordonnées du point d’intersection à ce gradient.
Pour définir le gradient, il suffit de calculer les dérivées partielles de l’équation.
Comme expliqué précédemment, le calcul des intersections nécessite l’utilisation des bases de Bernstein et de la méthode de Casteljau appliqué aux courbes de Bezier. Voici l’explication de cette méthode issue du cours de D.Michelucci.
1) La base de Bernstein permet de résoudre l’équation algébrique f(x)=0 avec x dans [0..1].
Les racines de f(x) dans [1, +∞] sont les inverses des racines du polynôme g(X)=0 avec X=1/x. Si f(x)=f0 + f1 x + … fn xn, alors g(X)=f0 xn + f1 Xn−1 + … fn (trivial).
Les racines négatives de f(x) sont les opposées des racines de h(X) avec X=−x. Les coefficients de h sont: h2i=f2i et h2i+1 = − f2i+1.
Il suffit donc de résoudre f(x)=0 avec x dans [0..1]; on convertit d’abord la forme canonique dans la forme de Bernstein, en utilisant le fait que (n st le degré du polynôme):
| xk = |
| C(k,i) / C(k,n) Bin(x) |
o C(i,n) est le nombre de façons de choisir i objets parmi n:
C(0,n)=1, C(i>n,n)=0, C(1,n)=n, C(i,n)=C(i−1,n)×(n−i+1)/i.
Soit f(x)=∑pi Bi(n)(x) la forme de Bernstein de f; alors la courbe de Bézier a pour points de contrôle: (i/n, pi) pour i de 0 à n.
La méthode isoler( pi, x0, x1), pour isoler les racines de f, avec initialement x0=0 et x1=1 est alors:
Si l’intervalle des pi ne contient pas 0, f n’a pas de racine dans l’intervalle courant x0, x1.
Si l’intervalle courant x0, x1 est d’amplitude inférieure à є=10−6, alors x0+x1)/2 est (une approximation d’une) racine de f;
sinon on coupe la courbe de Bézier en 2 par la méthode de Casteljau:
elle fournit les points de contrôle gi de la courbe entre 0 et 1/2 (x entre x0 et xm=(x0+x1)/2, et les points de contrôle di de la courbe entre 1/2 et 1 (x entre m et x1).
isoler( gi, x0, xm) et isoler( di, xm, x1) sont appelées récursivement. Au fur et à mesure, les racines trouvées sont stockées.
2) Les racines de f(x)=0 sont dans l’intersection entre l’enveloppe convexe des points de contrôle Pi=(xi=i/n,yi=pi) et l’axe y=0. On ne coupe plus bêtement en 2 la courbe de Bézier.
Donner une méthode de calcul de l’enveloppe convexe en 2D.
3) Cette méthode peut calculer l’intersection entre une droite
| x=x0+at, y=y0+bt, z=z0+ct |
et une surface implicite algébrique d’équation:
| F(x,y,z)= |
|
|
| Fαβγ xαyβzγ=0 |
Après substitution de x, y, z:
|
en utilisant le fait que: ∑i Pi ∑j Qj= ∑i ∑j Pi Qj.
Initialiser les fd à 0, pour d de 0 à n (le degré de la surface), et pour chaque monôme Fα,β,γxαyβzγ de F, faire: pour i=0… α, pour j=0… β, pour k=0… γ, incrémenter le coefficient fi+j+k du monome ti+j+k du polynome f(t) de la quantité
| Fα,β,γC(i, α) C(j, β) C(k,γ) x0α−iy0β−jz0γ−k aibj ck |
Ensuite, résoudre l’équation en t par la méthode présente.
Dans cette partie, nous allons vous présenter comme nous avons implémenté ces concepts pour créer notre moteur de rendu en lancer de rayon.
Nous avons choisi de développer le moteur en langage Ocalm, langage créer par l’INRIA couramment utilisé dans des applications nécessitant une grande quantité de calcul.
Notre programme permet d’afficher dans un fichier au format ppm le rendu du lancer de rayon de surfaces algébriques passées en paramètre.
Le moteur permet également d’appliquer les transformations spatiales usuelles aux surfaces comme la rotation, translation,mise à l’échelle ainsi que le choix de la couleur de l’objet.
La qualité du rendu n’étant pas une priorité, nous n’avons pas implémenter de modèle d’illumination complexe ni les principe de réflexion/réfraction.
De manière à passer ces équations en paramètre nous avons décidé d’utilisé l’analyseur lexical et syntaxique d’ocaml : Ocamllex Ocamlyacc , ce qui permet à un utilisateur de saisir, selon une syntaxe précise, l’équation de son choix ainsi que les valeurs de translation/rotation/mise à l’échelle et couleur à lui appliquer .
Nous allons tout d’abord présenter les différentes fonctions utilisé pour le programme et nous finirons sur la définition de l’analyseur. Vous trouverez le détail des fonctions dans le source du projet.
Pour débuter, il nous à fallu implémenter plusieurs outils mathématiques pour les calculs matriciels ou vectoriels.
Nous avons par conséquent créé les outils suivants :
(*| ) : Multiplication entre un scalaire et un vecteur.
(+| ) : Addition entre un scalaire et un vecteur.
(-| ) : Soustraction entre un scalaire et un vecteur
(*| | ) : Multiplication entre un scalaire et une matrice
(+| |) : Addition des composantes de deux matrices
(-| |) : Soustraction des composantes de deux matrices
dot : Produit scalaire de deux vecteurs.
unitise : Normalisation d’un vecteur
Pi : Valeur de pi
Choices : Permet de déterminer les combinaisons d’un ensemble
Puiss : Détermine la puissance d’un nombre
Cosinus : Détermine le cosinus
Sin : Détermine le sinus
Milieu : Calcul le milieu d’une distance
Gradient :Calcule le gradient d’une équation.
mat_mul : Calcul le produit de deux matrices.
Tous ces outils mathématiques nous permettrons d’effectuer tous les calculs nécessaire sur les vecteur et matrices intervenant dans le projet.
Nous allons passer à présent aux fonctions de mise en œuvre du lancer de rayon. Le projet se défini en plusieurs partie : la partie du moteur et la partie du calcul d’intersections. Pour présenter le programme, nous allons dérouler le fonctionnement du moteur en précisant chaque étape. Notre fonction principal let( ) consiste à écrire dans un fichier ppm le résultat du lancer de rayon en ayant en paramètre une ou plusieurs équations définies par notre analyseur syntaxique(nous reviendrons sur son fonctionnement plus tard). Cette fonction génère une série de rayon pour chaque pixel de l’image. Ces rayons sont crées avec une origine et une direction, puis on lance le calcul d’intersection avec les équations de la scène. Donc pour chaque rayon, on lance la fonction ray_trace avec en paramètre la lumière de la scène(light), le rayon courant ,ainsi que la scène( groupe d’équations).
Ray_trace lance la fonction intersection et en fonction du résultat, affiche la couleur du fond si aucun point n’est trouvé ou la couleur de l’objet si intersection renvoie un point. La couleur de l’objet est déterminée par la normal au point d’intersection et la direction de la lumière avec la formule intensité lumineuse au point = cos(normal au point, direction lumière).
La fonction intersection permet de renvoyer les coordonnés du point d’intersection entre un objet de la scène(équation) et un rayon s’il existe en le construisant à partir des racines (solutions) de f(x,y,z)=0 ainsi que la normal au point à l’aide de la fonction gradient.
Ces racines sont obtenues à partir de la fonction root qui est lancé dans le calcul d’intersection.
Root permet de récupérer la première solution trouvée dans la liste des racines ou de renvoyer l’infinie si la liste est vide.
La liste des racines est crée en lançant la fonction hit_ray_canonic qui créer un polynôme à partir du rayon et de l’équation et lance la fonction solve pour le résoudre. C’est la fonction solve qui permet de tester les racines.
La fonction solve lance deux fonctions qui sont root_neg_canonic et root_pos_canonic en concaténant les résultats.
Les deux fonctions root_neg_canonic et root_pos_canonic lance la fonction root_0_1_canonic sur le polynôme qui elle lance root_0_1_bernstein sur le polynôme en base de bernstein (à l’aide de la fonction to_bernstein)
Root_0_1_bernstein lance Decasteljau sur le polynôme dans la base de bernstein.
C’est la fonction Decasteljau qui calcule les racines par approximation. Comme préciser dans la sous partie résolution de la partie surface algébrique.
Ainsi en remontant jusqu’à la scène principale, on obtient les racines qui sont transformées en coordonnées à afficher.
Pour les transformations sur les équations (rotation,translation,mise à l’échelle), nous multiplions une matrice identité par les matrices de rotation/translation/mise à l’échelle en coordonnées homogènes nous récupérons donc une matrice de transformation globale que nous multiplions à l’origine et à la direction du rayon.
Comme dit précédemment, nous avons utilisé les modules lex et yacc d’ocaml. De manière à ce que l’utilisateur puisse saisir les équations qu’il souhaite et qu’il puisse définir les transformation et la couleur des objets de sa scène. Nous avons donc défini une série de mot et une syntaxe pour écrire les paramètres. Nous donnerons la syntaxe d’utilisation en fin de rapport. Pour utiliser le fichier de paramètre, il suffit d’écrire le nom du fichier en paramètre en lançant le programme. Dans ce fichier, il faudra écrire les équations et leur paramètre.
Syntaxe d’utilisation :
Exemple commenté :
Equation x*x+y*y+z*z(−1) [r 45,45,45;t 10,10,10;s0.5,0.5,0.5 ;c0.5,0.5,0.5]
=>équation de sphère
Pour définir une équation, on la note entre accolade, en précisant les valeur de x,y,z par leur degré et leur coefficient entre parenthèse .
exemple 22x =>(2)*x*x
−3xy2 =>(−3)*x*y*y
Pour définir les transformations et la couleur, on écrit entre crochet avec soit r pour rotation
T pour translation, s pour mise à l’échelle et c pour couleur.
Rotation est en degré et la couleur est entre 0et1 pour les composantes rvb .
On sépare chaque paramètres par un point-virgule et chaque valeur par une virgule.
Pour conclure, nous vous présentons les difficultés rencontrées ainsi que le future améliorations du projet .
Pour commencer, nous avons eu des difficultés principalement avec les méthode de résolution de polynôme de degré élevé. Il nous a fallut comprendre comment utiliser les polynômes dans la base de Bernstein pour ensuite se servir de la méthode de Casteljau pour approximé les solutions.
L’algorithme général du lancer de rayon n’a pas été la partie la plus délicate, la toile étant fournis d’exemple précis.
Un autres problème, était l’apprentissage du langage ocaml pour lequel nous n’étions guère familiarisé. Mais les exemples disposé en cours ainsi que sur la toile nous ont permis de mener à bien le projet.
Dernier problème, l’utilisation de l’analyseur syntaxique et lexical pour transcrire les équations dans notre code. N’étant encore que novice dans le domaine, il nous a fallu un certains temps pour comprendre le fonctionnement et la syntaxe .
Pour finir nous dirions dans l’ensemble que les fonctionnalités que nous avions décidées d’implémenter on été respectés.
Nous allons amélioré cette version du lancer de rayons en ajoutant la gestion de la réflexion et la transparence et l’utilisation d’arbre CSG sur des primitive standard pour le projet en Imagerie.
Flying Frog Consultancy
Université de Franche Compté
Université de Franche Compté
Projet Lancer de Rayon en c++
Exemple du lancer de rayon en c++ avec Bernstein
Calcul des racines réelles d’un polynôme
Ce document a été traduit de LATEX par HEVEA
