il est de temps d'écrire un peu de code pour tester tout ça. la structure générale va ressembler au code de départ de principes du lancer de rayons :
- charger une scène,
- extraire les triangles et leurs matières, éventuellement construire un BVH pour accélérer les intersections rayons / scène, cf construction de BVH optimal, SAH et parcours efficace,
- extraire les triangles dont la matière émet de la lumière, ie les sources de lumière,
- charger une camera.
on peut charger des objets stockés dans un fichier .obj (cf read_mesh() ou read_positions() ou read_meshio_data()) ou .gltf (cf read_gltf_scene()). les données de ces fichiers sont organisées différemment, et c'est une assez bonne idée de cacher ce détail dans une structure Scene pour éviter de rendre le reste du code dépendant du format de fichier...
au final, on va écrire une fonction qui parcourt tous les pixels de l'image, génère un rayon et utilise un estimateur Monte Carlo pour calculer la couleur du pixel. par exemple :
la ligne #pragma omp ... permet d'utiliser tous les coeurs de calcul du processeur pour calculer plus rapidement, cf le lancer de rayons, ça rame ? ou pas ?.
il ne manque plus que la partie Monte Carlo... on va avoir besoin de générer des nombres aléatoires entre 0 et 1, au minimum :
mais il faut manipuler à la fois le générateur et la distribution, ce qui n'est pas toujours très pratique. il suffit de les regrouper dans une structure :
on peut profiter des fonctions sample() et sample_range() pour vérifier que les nombres aléatoires ne sont pas égaux à 1 ou n, ce qui pourrait provoquer quelques problèmes numériques sinon... par exemple, choisir un indice de source dans un tableau de n éléments avec int s= rng.sample_range(sources.size()), on ne veut surtout pas que s soit plus grand que sources.size()...
remarque : on aurait aussi pu utiliser une std::uniform_int_distribution<int>(0, n-1) dans sample_range()...
dernière subtilité, le code de rendu est multi-threadé et les threads ne doivent pas modifier les mêmes variables simultanément, par exemple les générateurs aléatoires. le plus simple, dans ce cas, est tout simplement de créer un générateur local, ou privé pour chaque thread :
et le code Monte Carlo, ie l'évaluation de la couleur du pixel aura toujours la même forme :
et voila : on a tous les éléments pour évaluer un estimateur Monte Carlo...
éclairage ambiant
pour évaluer l'éclairage ambiant, la scène éclairée par un ciel uniforme, on utilise la formulation sur les directions de l'équation de rendu, et son estimateur, cf Monte Carlo et équation de rendu :
\[ L_r(p, \vec{o}) = \int_{\Omega} \frac{k}{\pi} V(p, \vec{l}) L_i(p, \vec{l}) \cos \theta \, dl \]
pour calculer l'éclairage ambiant, on considère que \( L_i(p, \vec{l})= 1\), ce qui simplifie un peu plus la formulation :
\[ \begin{eqnarray*} L_r(p, \vec{o}) &= & \int_{\Omega} \frac{k}{\pi} V(p, \vec{l}) \cos \theta \, dl\\ &\approx & \frac{1}{N} \sum_j^N \frac{k}{\pi} V(p, \vec{l_j}) \cos \theta \frac{1}{p(\vec{l_j})} \end{eqnarray*} \]
on va utiliser des directions générées selon la densité :
\[ \begin{eqnarray*} p(\vec{l})&= & \frac{1}{2\pi}\\ \cos \theta &= & u_1\\ \phi &= & 2\pi \, u_2\\ \vec{l} &= & ( \cos \phi \sin \theta, \, \sin \phi \sin \theta, \, \cos \theta)\\ \mbox{avec }\sin \theta &= & \sqrt{1 - \cos \theta \cos \theta} \end{eqnarray*} \]
une fois que l'on connait une direction \( \vec{l} \), que reste-t-il à évaluer ?
- la visibilite, \( V(p, \vec{l})\),
- le cosinus, \( \cos \theta \),
- la matière, \( \frac{k}{\pi}\), noté brdf habituellement,
- et bien sur la densité de proba \( p(\vec{l}) \).
si le ciel est visible, ie il n'y a pas d'intersection avec de la géométrie dans la direction \( \vec{l} \), il reste à calculer :
- le cosinus, \( \cos \theta \), on utilise la relation \( \cos \angle(\vec{n}, \vec{l})= \mbox{dot}( normalize(n), normalize(l) )\).
il faut aussi se rappeller que dans l'équation de rendu, on ne s'interresse qu'aux directions de l'hemisphere au dessus de la surface du point p.
le plus simple est de remplacer le cosinus, s'il est négatif, par 0 : \( \cos \angle(\vec{n}, \vec{l})= \mbox{std::max}(\mbox{float}(0), \mbox{dot}( normalize(pn), normalize(l) )) \) - la matière, \( \frac{k}{\pi}\), noté brdf habituellement,
- et la densité \(p(\vec{l}) = \frac{1}{2\pi} \).
la matière associée au point p et sa normale sont déjà connues, on recupère ces valeurs après les calculs d'intersection entre le rayon du pixel et la géométrie de la scène, cf les variables p, pn, et pmaterial dans le fragment de code au dessus.
pour constuire le rayon pour évaluer la visibilité \( V(p, \vec{l})\), il faut aussi penser à décoller l'origine du rayon de la surface de l'intersection pour éviter des gros défauts dans l'image, cf précision numérique et lancer de rayons
au final, c'est assez direct en ré-utilisant les fragments de code précédents :
mais il reste un problème : la direction aléatoire \( \vec{l} \) est construite dans un repère local, il faut la transformer dans le repère de la scène... la construction de la variable aléatoire utilise une convention : l'axe Z est aligné sur la normale à la surface, pour construire un changement de repère il faut connaitre 3 axes, on n'en connait qu'un seul...
indication : si on connaissait la normale \( \vec{n} \) et une autre direction \( \vec{d} \), non alignée avec la normale, il suffirait de 2 produits vectoriels pour construire 2 autres vecteurs orthogonaux à la normale : \( \vec{x} = \vec{d} \times \vec{n} \), puis \( \vec{y}= \vec{n} \times \vec{x} \) avec \( \vec{z}= \vec{n} \), ce qui définit une base xyz et permet de transformer la direction \( \vec{l} \) vers la scène.
et connaissant une base xyz, comment transformer la direction ? il faut se rappeler la définition des coordonnées d'un vecteur dans une base : ie \( x= \vec{v} \cdot \vec{x} \ldots \) et qu'un vecteur s'exprime comme la somme de ses projections sur les axes de la base : \( \vec{v}= x \, \vec{x} + y \, \vec{y} + z \, \vec{z}\). du coup \( l_x \, \vec{x} + l_y \, \vec{y} + l_z \, \vec{z}\) correspond bien à la direction \( \vec{l} \) dans le repère de la scène (si les axes ont des coordonnées dans le repère de la scène bien sur).
on peut essayer de construire une direction non alignée avec \( \vec{n} \) en modifiant les coordonnées de \( \vec{n} \) , mais ca n'est pas toujours très robuste. par contre, il existe une construction complètement différente, à base de quaternions, qui est plus rapide : cf
"Building an orthonormal basis from a 3d unit vector without normalization", J.Frisvad, 2012.
mais cette solution souffre aussi d'un problème numérique qui a été corrigé depuis :
"Building an Orthonormal Basis, Revisited" Pixar, 2017
et qui se code sans problèmes :
et voila ! on a tous les éléments pour calculer l'image, il suffit de modifier le fragment qui calcule la couleur du pixel :
alors c'est faisable ?
bien sur, on peut créer des fonctions utilitaires pour générer une direction et évaluer une densité de proba. mais cet exemple est sans doute plus lisible rédigé comme ça.
éclairage direct
même démarche, on part de la formulation du problème, on écrit l'estimateur Monte Carlo et on choisit une densité de proba pour les échantillons.
on utilise la formulation de l'équation de rendu sur la surface des sources de lumière, cf Monte Carlo et éclairage direct :
\[ \begin{eqnarray*} L_r(p, o) &= & \int_{S} \frac{k}{\pi} V(p, q) L_i(p, q) \cos \theta_p \frac{\cos \theta_q}{||\vec{pq}||^2} \, dq\\ &\approx & \frac{1}{N} \sum_j^N \frac{k}{\pi} V(p, q_j) L_i(p, q_j) \cos \theta_p \frac{\cos \theta_{q_j}}{||\vec{pq_j}||^2} \, \frac{1}{p(q_j)} \end{eqnarray*} \]
pour générer des points sur les sources de lumière, on peut utiliser cette densité de proba :
\[ p(q) = \frac{1}{n} \frac{1}{aire} \]
c'est la 1ère solution proposée dans Monte Carlo et éclairage direct : on choisit uniformément une source parmi n et ensuite un point sur cette source. connaissant le point \( q \) et sa densité, comment évaluer la fonction intégrée ? on peut évaluer directement les cosinus et le carré de la distance entre p et q. il reste \( L_i(p, q) \), la lumière émise par q qui éclaire p, c'est tout simplement l'emission de la source de lumière que l'on a choisit pour générer q. le dernier terme, la visibilité entre p et q s'évalue directement en vérifiant l'existence d'une intersection sur le rayon entre p et q.
exercice pour les curieux : comment utiliser l'autre solution présentée dans Monte Carlo et éclairage direct, ie choisir une source en fonction de son aire ? que faut-il modifier ?
on va aussi simplifier la représentation des sources : ce sont les triangles de la scène dont la matière émet de la lumière. on peut les décrire avec une structure qui regroupe quelques propriétés utiles :
avec ces éléments, le code de génération d'un point sur une source (un triangle dans ce cas) est sans surprise :
il ne reste plus qu'à finir d'écrire l'estimateur complet :
mêmes remarques qu'au dessus, il faut penser à décoller l'origine et l'extrémité du rayon des surfaces pour éviter des défauts d'intersections. autre subtilité, il faut faire bien attention à l'orientation des vecteurs pour évaluer les cosinus.
on pourrait aussi créer une fonction sample( ) dans Source qui génère un point sur la source pour rendre tout ça un peu plus lisible... par exemple :
l'estimateur complet est plus compact, du coup :
