TP4 - BVH et trop de triangles

retour sur la perspective :
La construction des rayons est très simple, mais la perspective est assez brutale et déforme énormément les objets lorsqu'ils se rapprochent de la camera. Comment utiliser une ouverture plus petite que 90° ? L'image au dessus utilise une perspective plus naturelle avec une ouverture de 60°. Comment modifier la construction de l'extrémité des rayons ? Ou : comment zoomer dans l'image ?

pour les curieux : comment simuler une profondeur de champ ? ie comment rendre les objets flous s'ils ne sont pas à la distance du plan focal / du plan image ?


exercice 1 : beaucoup de triangles
Chargez une scène plus complexe que la cornell box, en utilisant read_meshio_data().

par exemple alley, vous pouvez utiliser la transformation suivante pour placer les objets devant la camera :

Translation( 180.638168, 109.031326, -1229.029297 ) * RotationX( 13 ) * RotationY( 275 );

   
Parcourez les triangles et conservez ceux dont l'aire est non nulle. oui, il y a souvent des triangles dégénérés (avec 2 sommets superposés).

rappel : scènes d'exemple, cf doc en ligne.

si vous preférez partir d'un code départ : cf scene_meshio.cpp


exercice 2 : BVH
Reprenez le principe de construction de l'arbre dans la doc en ligne.
Vous pouvez soit trier les triangles sur un axe, soit les répartir par rapport au centre de l'axe le plus étiré. Comparez votre solution avec votre voisin.

astuce : il est plus robuste de manipuler les centres des boites des triangles pour répartir.

pour les rapides : expérimentez avec une des constructions qui minimisent le coût de parcours du BVH, cf doc en ligne.

pour les curieux : peut-on déterminer que certains noeuds sont mauvais et qu'il faudrait répartir les triangles dans plus de fils / noeuds ? Que faudrait-il modifier ?
question : que pourrait-être la définition d'un "mauvais" noeud ?
indications : à partir de quelle taille / aire est-il plus intéressant de remplacer un fils par ses fils ? Ou : dans quel cas est-il intéressant de créer un noeud avec 3 fils, au lieu de 2...

si vous preférez partir d'un code départ : cf bvh_correction.cpp


exercice 3 : intersection rayon / BVH
Ecrivez le parcours du BVH permettant de visiter toutes les feuilles touchées par le rayon.


exercice 4 : intersection ordonnée
Adaptez l'ordre de parcours des fils d'un noeud en fonction de la direction du rayon. Pourquoi ? Constatez-vous une différence ?


exercice 5 : transformations BVH et rayons
Dans les tps précédents, il était probablement plus simple de pre-transformer les sommets des objets pour les placer devant la camera.
Que peut-on dire de la qualité d'un BVH construit avec des boites alignées sur les axes lorsqu'on applique une rotation à la géométrie ?

indication : sait-on estimer le nombre de droites visitant une région de la scène ? Si cette région est plus grande, comment évolue le nombre de droites ?
question : comment construire l'englobant d'une boite alignée après une rotation ? Comment évolue l'englobant ?

Si on ne souhaite pas pre-transformer les objets, il suffit de transformer l'origine et l'extrémité des rayons issus de la camera avant de calculer les intersections. Si on applique une transformation T pour passer du repère de l'objet au repère de la camera, quelle transformation appliquer à l'origine et l'extrémité du rayon pour passer du repère camera vers le repère de l'objet ?

Vérifiez que le BVH construit sur des objets non transformés est plus efficace que celui construit sur les objets pre-transformés dans le repère de la camera.

Partie 2 : estimateurs et réduction de variance

un peu de lecture pour changer !

exercice : on souhaite intégrer une somme pondérée de fonctions. on suppose que l'on connaît une distribution bien adaptée à l'intégration numérique de chaque fonction. Quelles sont les différentes stratégies pour construire un estimateur de Monte Carlo de la somme pondérée de fonctions ?

les détails sont dans :


à cogiter pour la suite : à partir du moment ou l'on commence à jouer avec les distributions des échantillons, il faut être préparé aux défauts introduits dans certaines situations... on a peut être grandement réduit la variance dans la majorité des cas, mais il reste sans doute des situations mal gérées. Ces situations mal échantillonnées, avec une distribution peu adaptée à la fonction intégrée peuvent introduire des évaluations avec des pondérations aberrantes qui cassent complètement la convergence de l'estimateur. Il y a deux familles de solutions, soit toujours utiliser une stratégie de défense, ie uniforme, dans les estimateurs "multiple importance sampling" au prix de la vitesse de convergence, soit réduire a posteriori l'impact de ces pondérations aberrantes, soit ignorer le problème...

L'article ci-dessous illustre très bien le problème de pondération aberrante et propose de construire plusieurs estimateurs au lieu d'un seul et d'utiliser leur médiane comme estimateur final. C'est plutôt élégant, mais peu adapté lorsque l'on se limite à un faible nombre d'échantillons.
"Firefly removal in Monte Carlo rendering with adaptive Median of meaNs", J. Buisine, S. Delepoulle, C. Renaud, 2021


pour les curieux : l'utilisation des textures comporte également une intégration, une solution alternative est présentée dans :
"Filtering After Shading With Stochastic Texture Filtering", M. Pharr, B. Wronski, M. Salvi, M. Fajardo, I3D 2024.

---

Annexe : génération de directions pour le modèle Blinn Phong

Le modèle de matière se présente comme une mixture de 2 termes, une réflexion diffuse ainsi qu'une réflexion spéculaire :

\[
    f_r( \vec{d}, \vec{l} )= \frac{diffuse}{\pi} + specular \times \frac{8+\alpha}{8 \pi} ( \cos \theta_h )^\alpha \\
\]

rappel : la construction du modèle se trouve sur la doc en ligne.

En pratique, dans l'équation de rendu, le modèle de matière est toujours pondéré par le cosinus de l'angle entre la normale en p et la direction \( \vec{l} \) :
\[
    L_r(p, \vec{o})= L_e(p, \vec{o}) + \int_{\Omega} L_i(p, \vec{l}) \, f_r( \vec{o}, \vec{l} ) \cos \theta dl
\]
On peut ré-ecrire l'équation de rendu explicitement en décomposant les 2 termes du modèle de matière :
\[
    L_r(p, \vec{o})= L_e(p, \vec{o}) + \int_{\Omega} L_i(p, \vec{l}) \, \frac{diffuse}{\pi} \, \cos \theta dl + \int_{\Omega} L_i(p, \vec{l}) \, \frac{8+\alpha}{8 \pi} ( \cos \theta_h )^\alpha \, \cos \theta dl
\]

On peut écrire directement un estimateur utilisant une distribution de directions proportionnelle à \( \cos \theta \) pour l'intégration de la partie diffuse, par contre, il va falloir cogiter un peu pour construire une distribution proportionnelle à \( (\cos \theta_h)^\alpha \cos \theta \) pour la partie spéculaire.

En consultant le recueil de distributions GI Compendium, P. Dutre, 2003, il est assez direct d'identifier une distribution de directions proportionnelle à \( \cos \theta \) présentée en eq 35 :
\[
\begin{eqnarray*}
    p( \vec{v} ) & = & \frac{\cos \theta}{\pi} \\
    \cos \theta & = & \sqrt{u_1} \\
    \sin \theta & = & \sqrt{1 - (\cos \theta)^2} \\
    \phi & = & 2\pi u_2
\end{eqnarray*}
\]


On peut également trouver, cf eq 36, une distribution adaptée à \( (\cos \theta)^n \) :
\[
\begin{eqnarray*}
    p( \vec{v} ) & = & \frac{n+1}{2\pi} (\cos \theta)^n\\
    \cos \theta & = & u_1^{\frac{1}{n+1}} \\
    \sin \theta & = & \sqrt{1 - (\cos \theta)^2} \\
    \phi & = & 2\pi u_2
\end{eqnarray*}
\]

Une solution assez simple consiste à générer la direction \( \vec{h} \) avec cette distribution et ensuite à construire la direction \( \vec{l} \) qui doit être symétrique de \( \vec{o} \) par rapport à \( \vec{h} \).
rappel : \( \vec{h} \) doit être la direction bissectrice de  \( \vec{o} \) et \( \vec{l} \) dans le modèle Blinn-Phong.
\[
    \vec{l} = 2(\vec{o} \cdot \vec{h}) \vec{h} - \vec{o}
\]

Mais bien sur, comme on transforme \( \vec{h} \), la direction générée, pour obtenir \( \vec{l} \), il faut aussi évaluer le changement de mesure pour déterminer la densité de probabilité de \( \vec{l} \) :
\[
    p( \vec{l} ) = p( \vec{h} ) \frac{1}{4 \cos \angle( \vec{o}, \vec{h} )} = p( \vec{h} ) \frac{1}{4 \cos \angle( \vec{l}, \vec{h} )}
\]

c'est sans doute plus lisible en utilisant directement le produit scalaire des directions :
\[
    p( \vec{l} ) = p( \vec{h} ) \frac{1}{4(\vec{o} \cdot \vec{h})} = p( \vec{h} ) \frac{1}{4(\vec{l} \cdot \vec{h})}
\]

pour les curieux : les détails du calcul du changement de variable entre \( \vec{h} \) et \( \vec{l} \) sont dans Pbrt, chapitre 9, par exemple.

Dernier "détail", à cause de la transformation de \( \vec{h} \) vers \( \vec{l} \), \( \vec{l} \) peut assez souvent se trouver sous la surface, ie \( \vec{n} \cdot \vec{l} < 0 \), dans ce cas, la densité de proba doit être nulle et le modèle de matière noir.


Bilan on vient de construire une distribution de directions proportionnelle au terme\(  ( \cos \theta_h )^\alpha \) du modèle Blinn-Phong. La génération n'est pas parfaite, ie certaines directions générées sont aberrantes et on ne tient pas compte du \( \cos \theta \). Mais malgré ces défauts cette solution est largement préférable à une distribution uniforme ou simplement proportionnelle à \( \cos \theta \) et devrait permettre de très fortement réduire la variance dans l'exemple précédent, ie l'image avec la boule rouge. A vérifier !


exercice : calculez l'albedo / la reflectance hémisphérique de vos matières afin de vérifier le code calcul que vous venez d'écrire et de vérifier que votre modèle de matière est correctement normalisé. La définition se trouve dans GI Compendium eq 64, ou Pbrt, chaptire 4  :
\[
    \rho= \frac{1}{\pi} \int_\Omega \int_\Omega f_r( \vec{o}, \vec{l} ) \cos \theta_o \cos \theta_l \, dl do \mbox{ avec } \theta_o : \angle(\vec{n}, \vec{o}) \mbox{ et } \theta_l : \angle(\vec{n}, \vec{l})
\]

ou pour une seule direction \( \vec{o} \) :
\[
    \rho( \vec{o} ) = \int_\Omega f_r( \vec{o}, \vec{l} ) \cos \theta_l \, dl \mbox{ avec } \theta_l : \angle(\vec{n}, \vec{l})
\]

vérifiez que vos estimations convergent bien vers la même valeur quelque soit la distribution utilisée pour générer les directions, ie uniforme, \( \cos \) ou \( \cos^n \). Ces estimations devraient être inférieures ou égales à 1 dans tous les cas.

indications : écrivez systématiquement un couple de fonctions pour chaque distribution, une pour générer une direction, l'autre pour évaluer sa densité de probabilité. Ce sera plus simple pour utiliser les combinaisons d'estimateurs de "multiple importance sampling" ou "importance resampling" ou pour tester plusieurs distributions, tout simplement. par exemple :

Vector sample_cos( const float u1, const float u2 ) { ... }
float pdf_cos( const Vector& v ) { ... }

Vector sample_uniform( const float u1, const float u2 ) { ... }
float pdf_uniform( const Vector& v ) { ... }

exercice : pour intégrer tout le modèle, ie les 2 termes, diffus et spéculaire, comment construire une seule distribution correcte / adaptée ?
indication : la distribution sera très probablement une mixture, comment écrire l'estimateur Monte Carlo dans ce cas ?


exercice : comment utiliser "importance resampling" pour construire une distribution proportionnelle à \( (\cos \theta_h)^\alpha \cos \theta \) ? Même question pour une distribution proportionelle à \( f_r( \vec{o}, \vec{l} ) \cos \theta \) ?


pour les curieux : pourrait-on utiliser un estimateur différent dans ce cas ? Quel autre estimateur Monte Carlo permet de normaliser la densité de probabilité a posteriori ? (ce n'est pas celui utilisé par "importance resampling"...)