on vient de voir dans Monte Carlo et équation de rendu que l'on peut reformuler l'équation de rendu sur \( S \), l'aire des sources de lumières, et écrire un estimateur Monte Carlo :
\[ \begin{eqnarray*} L_r(p, o) &= & L_e(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 & L_e(p, o) + \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*} \]
mais, bien sur, il reste quelques détails à régler pour échantilloner / choisir des points sur plusieurs sources de lumière...
comme d'habitude pour estimer cette intégrale, cf intégration numérique et Monte Carlo, il faut définir une densité de proba \( p(q) \) et générer des points selon cette densité...
que peut-on faire ? on va commencer par le plus direct : choisir une source parmi les \( n \) sources existantes puis un point sur la source sélectionnée. la densité de proba est le produit des 2 :
\[ p(q) = \frac{1}{n} \frac{1}{aire} \]
ie on choisit d'abord une source uniformément parmi les \( n \) avec une variable aléatoire discrète, puis un point sur la source sélectionnée avec une variable aléatoire continue.
le code est direct, on suppose que l'on a décrit chaque source, des carrés, par un point et 2 aretes :
et il ne reste plus qu'à sélectionner une source et un point en utilisant 3 nombres aléatoires uniformes entre 0 et 1, comme les autres fois :
et voila, sur une scène avec plusieurs sources de lumières :
pour \( N = 16, 64 \) :
 |  |
argh ! mais pourquoi la partie droite de l'image est toute moche ?
il suffit de regarder ou se trouvent les sources de lumière et comment sont distribués les points \( q \) que l'on vient de construire pour comprendre le problème. il y a 2 panneaux lumineux qui éclairent la scène, mais celui de gauche est découpé en 100 petits panneaux. il y a 101 sources de lumière au total...
et en observant comment sont distribués les points sur les sources (les points en rouge dans les vignettes), on se rend vite compte que presque tous les points sont choisis sur les 100 petits panneaux. ce qui n'est pas très surprenant finalement...
 |  |
c'est exactement ce qu'il faut faire dans la partie gauche de l'image, quand les points \( p \) sont proches des petits panneaux. mais pas du tout adapté pour l'autre partie de l'image lorsque les points \( p \) sont proches du grand panneau. on vient de placer tous les échantillons utilisés pour estimer l'intégrale dans une petite partie du domaine qui ne contribue pas beaucoup au résultat... et on peut constater que l'estimateur ne converge pas très bien.
il va falloir construire une densité de proba différente pour régler ce problème...
éclairage direct uniforme
dans l'exemple précédent, on souhaiterait plutôt répartir autant de points sur les 2 panneaux, ou sélectionner un panneau / une source de lumière en fonction de son aire. comment faire ?
il va falloir construire une variable aléatoire discrète pour sélectionner une source. mais bien sur, au lieu d'utiliser une sélection uniforme, on va en profiter pour tenir compte de l'aire de chaque source, pour sélectionner 2 fois plus souvent une source 2 fois plus grande... on a déjà vu un exemple simple dans intégration numérique et Monte Carlo, pour sélectionner un intervalle 3 fois plus souvent qu'un autre.
pour construire la variable aléatoire discrète, on va imposer \( p(i)= \frac{aire_i}{\sum_j aire_j} \), ie l'aire normalisée de la ième source (ne pas oublier qu'une densité de proba est normalisée, donc il faut diviser l'aire de la source par la somme des aires des sources), puis construire la fonction de répartition et enfin la parcourir pour sélectionner la source qui correspond à un nombre aléatoire uniforme entre 0 et 1. pour générer un point sur cette source, on utilise la même solution qu'au dessus.
au final, la densité de proba s'écrit comme un produit de 2 densités, choisir une source en fonction de son aire et choisir un point sur la source sélectionnée :
\[ p(q) = \frac{aire}{\sum_j aire_j} \frac{1}{aire} = \frac{1}{\sum_j aire_j} \]
le code est sans surprise, en plusieurs étapes : d'abord sommer les aires des sources pour normaliser la densité de proba, puis construire la fonction de répartition (ou plutot l'histogramme cumulé vu que c'est une variable aléatoire discrète), c'est l'application directe de la définition :
et voila, il ne reste plus qu'à chercher l'indice de la source associée à une valeur aléatoire uniforme entre 0 et 1 :
et on peut faire exactement la même chose sans stocker la fonction de répartition dans un vecteur dynamique, pour gagner un peu de temps pendant le calcul de l'image (ie les allocations dynamiques dans un code multi-threadé sont encore plus lentes que d'habitude, donc si on peut les éliminer du code, c'est mieux...)
relisez intégration numérique et Monte Carlo si ce n'est pas clair, ie la fonction de répartition, le parcours... l'exemple n'utilise que 2 valeurs, c'est plus simple pour se faire une idée.
et alors ? pour \( N = 16, 64 \) :
 |  |
un peu mieux non ?
voila ou se trouvent les points \( q \) pour un point \( p \) dans la partie gauche de l'image, ie devant les petits panneaux et pour un autre devant le grand panneau :
 |  |
on a bien réussi à répartir les points sur les 2 panneaux et l'estimateur converge mieux, le bruit est mieux réparti dans l'image... mais peut être que si on arrivait à placer un nombre de points en fonction de l'aire des sources vues depuis le point \( p \) (ie en fonction du nombre de pixels blancs dans les vignettes) se serait encore mieux...
et si les sources sont des triangles ?
il faut trouver comment choisir uniformément un point sur l'aire d'un triangle. la construction n'est pas directe, mais il y a une solution en utilisant les coordonnées barycentriques, cf principes du lancer de rayons. en construisant 3 densités de proba continues, une pour chaque coordonnée barycentrique, leurs fonctions de répartition et en les inversant pour enfin générer un point dans un triangle abc avec 2 nombres aléatoires uniformes \( u_1, u_2\) entre 0 et 1 :
\[ \begin{eqnarray*} p(q) &= & \frac{1}{aire}\\ \alpha &= & 1 - \sqrt{u_1}\\ \beta &= & (1 - u_2) \sqrt{u_1}\\ \gamma &= & u_2 \sqrt{u_1}\\ q &= & \alpha a + \beta b + \gamma c \end{eqnarray*} \]
ce qui se code directement :
solution alternative
il existe également une construction géométrique, sans doute plus simple à comprendre, cf
"A Low-Distortion Map Between Triangle and Square", E.Heitz, 2019.
l'idee est différente : on génére un point dans le quadrilatère englobant le triangle, et on écrase ce quadrilatère, ie on pousse le sommet externe sur la diagonale du triangle, pour le transformer en triangle. le code est assez simple aussi :
