TD4 - Échantillonnage

L'intégration de Monte-Carlo est au cœur de nombreux algorithmes de synthèse d'image photo-réalistes. Dans ce TD, après un échauffement en 1D, vous l'utiliserez dans le cas de l'anti-aliasing en sur-échantillonnant les pixels et dans le cas des sources étendus pour obtenir des ombres douces.

1. Intégration de Monte-Carlo 1D

L'objectif est de calculer l'intégrale de la fonction 1D : $f(x) = 5x^4$ sur l'intervalle $[0,1]$, c'est-à-dire $I = \int_0^1 5x^4 dx$. Téléchargez l'archive MonteCarlo.zip et compilez son contenu. Implémentez la fonction int_mc qui doit retourner l'estimation de l'intégrale selon l'estimateur de Monte-Carlo $\langle I \rangle$ avec une densité de probabilité uniforme $(\forall i,~pdf(x_i)=1)$, soit : $$\langle I \rangle = \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_i)} = \frac{1}{N} \sum_{i=1}^N 5x_i^4$$ Ignorez pour l'instant la variable de sortie varest.

Le main appelle int_mc avec un nombre croissant d'échantillons puis dessine avec gnuplot l'erreur d'estimation en fonction de $N$. Exécutez plusieurs fois le programme et comparez les tracés.

Calculez maintenant la variance analytique en fonction de N : $$\sigma^2_{est} = \frac{1}{N} \int_0^1 \left(\frac{f(x)}{pdf(x)} - I \right)^2 pdf(x) dx = \frac{1}{N} \int_0^1 (5x^4 - 1)^2 dx$$ Tracez sur le graphique précédent cette variance.

Estimez enfin la variance numérique (variable de sortie varest) à l'aide de l'estimateur : $$ \langle \sigma^2_{est} \rangle = \frac{1}{N-1} \left\{ \frac{1}{N} \sum_{i=1}^N \left(\frac{f(x_i)}{pdf(x_i)}\right)^2 - \left( \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_i)}\right)^2 \right\} $$ Affichez-la sur le même graphique.

Monte-Carlo 1D
Erreur d'estimation, variance analytique et numérique en fonction du nombre d'échantillons.

2. Anti-Aliasing

Les silhouettes des objets et les textures souffrent clairement de problèmes d'aliasing. Pour le résoudre, il faut intégrer la radiance sur le pixel, ce que l'on peut faire au sens de Monte-Carlo en lançant plusieurs rayons par pixel et en moyennant leur contribution.

2.1. Échantillonnage régulier

Modifiez la méthode Viewer::renderBlock pour qu'elle lance plusieurs rayons par pixel à travers une grille régulière de taille $n \times n$, où $n$ correspondra à l'attribut sampleCount de la caméra (lu depuis le fichier de scène). Vous devriez obtenir le résultat suivant avec "tw.scn" en fixant le nombre d'échantillons à 8 dans le fichier de scène :

antialised
Échantillonnage régulier – 8x8 rayons par pixel

2.2. Échantillonnage stratifié

Implémentez ensuite un échantillonnage stratifié en ajoutant du jitter à chaque rayon. Vous pouvez utiliser la fonction :

Eigen::internal::random<float>(a,b)

pour obtenir un nombre aléatoire entre a et b. Activez la stratification lorsque Camera::isStratified() renvoie true (également lu depuis le fichier de scène).

antialised
Échantillonnage stratifié – 8x8 rayons par pixel

La scène "deuxSpheres.scn" est particulièrement difficile en raison de la texture en damiers du plan infini :

2x2 8x8
32x32 128x128
Échantillonnage stratifié – 2x2, 8x8, 32x32 et 128x128 rayons par pixel

2.3. Bonus : filtre de reconstruction

On peut encore améliorer les résultats en pondérant la contribution des rayons correspondant à un pixel par un filtre de reconstruction, par exemple par un noyau gaussien (centré sur le milieu du pixel). Pour plus de flexibilité, vous pouvez coder une classe Filter qui donnera les poids ainsi que les décalages (par rapport au centre du pixel) associés à chaque rayon.

3. Sources étendues

Les sources ponctuelles et directionnelles sont rarement physiquement réalistes. La majorité des sources lumineuses est surfacique. Calculez les ombres portées pour de telles sources étendues nécessite cependant d'intégrer leur visibilité depuis le point d'évaluation. Pour ce faire, vous utiliserez à nouveau un échantillonnage de Monte-Carlo.

Parcourez la classe AreaLight représentant une source lumineuse étendue rectangulaire, potentiellement texturée. Analysez en particulier la méthode :

Color3f AreaLight::intensity(const Point3f &x, const Point3f &y) const
calculant l'intensité émise depuis le point $\vec{y}$ sur la source vers le point $\vec{x}$ de la scène.

3.1. Échantillonnage de Monte-Carlo

Modifiez l'intégrateur de Whitted pour tirer un rayon d'ombre dans une direction aléatoire vers la source étendue. Pour déterminer cette position, utilisez les vecteurs unitaires AreaLight::uVec et AreaLight::vVec ainsi que la position et la taille de la source. Notez que l'on peut détecter qu'une Light est une AreaLight avec un dynamic_cast :

1
2
3
4
5
6
7
8
9
if(dynamic_cast<const AreaLight*>(lights[i]))
{
   const AreaLight* light = dynamic_cast<const AreaLight*>(lights[i]);
   // source étendue
}
else
{
   // lampe ponctuelle ou directionnelle
}

Vous devriez obtenir ce type de résultat pour la scène "tw_area.scn" :

arealight
Échantillonnage stratifié – 8x8 rayons par pixel

Testez également avec la scène "killeroo_area.scn" :

killeroo   killeroo  
Échantillonnage stratifié source de taille 0.5x0.5 (gauche) et 0.1x0.1 (droite) – 8x8 rayons par pixel

3.2. Source texturée

On veut moduler l'intensité de la source étendue par la texture "light_source.exr".

Pour cela, modifiez AreaLight::intensity pour accéder à l'attribut Bitmap m_texture. Ajoutez également la ligne suivante au fichier de scène "tw_area.scn" :

<light type="areaLight">
  <string name="texture" value="light_source.exr" />
  <!-- le reste est identique -->
</light>

Avant de tester, essayez de prédire quel sera le résultat ...afficher le résultat :

coloredLight   coloredLight2
Échantillonnage stratifié – 16x16 rayons par pixel– source étendue texturée

3.3. Bonus : échantillonnage stratifié

Les ombres douces générées sont très bruitées car, si l'échantillonnage des pixel est stratifié, celui de la source étendue est purement aléatoire. À titre de comparaison, implémentez l'échantillonnage stratifié de la lampe en tirant $n \times n$ rayons d'ombres sur une grille avec du jittering. Pour que la comparaison des ombres soit juste, n'utilisez qu'un rayon par pixel.

Vous devriez obtenir ce résultat pour "tw_area.scn" :

arealight
Échantillonnage stratifié – 16x16 rayons d'ombre (1 rayon par pixel)

...et pour "killeroo_area.scn" :

killeroo
Échantillonnage stratifié – 32x32 rayons d'ombre (1 rayon par pixel)

Il est évidemment possible de combiner les deux approches, puis d'ajuster le nombre d'échantillons par pixel et par source étendue pour obtenir la qualité souhaitée.

killeroo
Échantillonnage stratifié – 8x8 rayons par pixel et d'ombre

4. Ambiante Occlusion

L'ambiant occlusion est une technique de rendu faisant l'hypothèse qu'une surface (diffuse) reçoit une illumination uniforme depuis toutes les directions et que la visibilité est le seul effet important. Certaines portions de la surface vont recevoir moins de lumière que d'autres en raison des occlusions, et donc apparaître plus sombres. Formellement, l'ambiant occlusion $L$ au point $\vec{x}$ est l'intégrale de la fonction de visibilité $V$ sur l’hémisphère supérieur centré en ce point :

$$L(\vec{x}) = \int_{\Omega}V(\vec{x}, \vec{x}+\alpha\omega)\,\frac{\cos\theta}{\pi}\,\mathrm{d}\omega$$

où $\theta$ est l'angle entre la direction $\omega$ et la normale au point $\vec{x}$. La variable $\alpha$ permet d'ajuster la distance maximale jusqu'à laquelle l'effet d'occlusion est considérée.

4.1. Échantillonnage d'un hémisphère

Pour calculer cette intégrale au sens de Monte-Carlo, il vous faut échantillonner l'hémisphère au point $\vec{x}$ orienté selon la normale $\vec{n}$. Pour cela, vous allez compléter la classe Warp qui permet de transformer un échantillonnage uniforme du carré unité $[0,1[ \times [0,1[$ dans un autre domaine.

Pour vous échauffer, implémenter les fonctions Warp::squareToUniformDisk et Warp::squareToUniformDiskPdf en utilisant le cours. Pour déboguer votre code, utilisez l'utilitaire warptest qui vous permet de visualiser l'effet des fonctions de warping sur un nuage de points et de les valider par un test statistique du $\chi^2$. Vous devriez obtenir :

squareToUniformDisk   squareToUniformDiskPdf
Échantillonnage uniforme de l'hémisphère autour du pôle $(0,0,1)$

Effectuez la dérivation et implémentez ensuite la fonction Warp::squareToUniformHemisphere qui effectue la transformation vers la distribution uniforme sur l'hémisphère autour du pôle $(0,0,1)$. Complétez également Warp::squareToUniformHemispherePdf qui retourne la $pdf$ associée. Vous devriez obtenir :

squareToUniformHemisphere   squareToUniformHemispherePdf
Échantillonnage uniforme de l'hémisphère autour du pôle $(0,0,1)$

Implémentez ensuite les fonctions Warp::squareToCosineHemisphere et Warp::squareToCosineHemispherePdf effectuant la transformation vers la distribution sur l'hémisphère pondérée par le cosinus (cosine-weighted hemisphere sampling).

squareToCosineHemisphere   squareToCosineHemispherePdf
Échantillonnage pondéré par le cosinus de l'hémisphère autour du pôle $(0,0,1)$

4.2. Calcul de l'AO

Complétez l'intégrateur AO évaluant la fonction de visibilité en lançant m_sampleCount rayons d'ombre dans des directions aléatoires $\vec{p}$ sur l'hémisphère. Pour cela, utilisez la fonction Warp::squareToUniformHemisphere précédente et transformez la direction aléatoire $\vec{d}$ obtenue dans l'espace local vers l'espace monde : $$ \vec{p} = \vec{d}.x() \, \vec{u} + \vec{d}.y() \, \vec{v} + \vec{d}.z() \, \vec{n}$$ où $\vec{u}$ et $\vec{v}$ sont deux vecteurs orthogonaux dans le plan tangent à $\vec{n}$.

Moyennez les réponses de la fonction de visibilité et n'oubliez pas de pondérez par la $pdf$ adéquate. Vous pouvez faire l'hypothèse que l'occlusion est un effet global ($\alpha = \infty$). Testez également avec Warp::squareToCosineHemisphere lorsque m_cosineWeighted est à True. Vous devriez obtenir ces deux résultats pour "killeroo_ao.scn" :

AO   AO  
Échantillonnage uniforme (à gauche) et pondéré par le cosinus (à droite) avec 256 échantillons

Selon vous, quelle est la stratégie d'échantillonnage la plus efficace en termes de convergence et de performances ? Pourquoi ?