M2 Image

TP2 - lancer de rayon
et Monte Carlo

L'objectif de ce tp est de mettre en place les éléments nécessaires pour une simulation "correcte" de l'éclairage direct et indirect basé sur les méthodes de Monte Carlo.

Utilisez les scenes de test suivantes, elles sont composées de très peu de geométrie (triangles), ce qui permet d'éviter de construire une structure accélératrice pour limiter les calculs de visibilité (au moins dans un premier temps ...) :

cornell box + matières + orbiter
geometry
shadow
emission
secondary
matières nécessaires pour charger les scènes de test.

Partie 1 : lancer de rayons et visibilité

exercice 1 : visualisation par lancer de rayons

Pour chaque pixel de l'image résultat, déterminez les coordonnées de l'origine du rayon dans le repère de la scène, ainsi que celles de son extrémité.
Dans quel espace est-il le plus simple de calculer l'origine et l'extrémité du rayon ?

Calculez l'intersection du rayon avec tous les triangles de la scène (cf. la fonction intersect() dans tuto_ray.cpp) et ne conservez que l'intersection valide la plus proche de l'origine du rayon.

remarques :

Ray permet de représenter et de manipuler facilement un rayon.
Hit permet de représenter simplement une intersection,
les fonctions d'intersections rayon / triangle est déjà écrite : Triangle::intersect( ).

indications : votre programme devrait ressembler à :

main
    creer une image
    charger un maillage et ses matières

    extraire les sources de lumiere (trouver les triangles dont la matière émet de la lumiere)

    pour chaque pixel de l'image :
        générer un rayon dans le repere de la scene

        // trouver le point de la scène visible pour le rayon
        pour chaque triangle :
            transformer le rayon dans le repere local de l'objet
            calculer l'intersection du rayon avec le triangle,
            conserver l'intersection si elle est valide et plus proche que la dernière trouvée

        si une intersection valide existe
            ecrire un pixel blanc dans l'image

    sauver l'image

exercice 2 : matière

Récupérez la couleur de la matière du triangle touché par chaque rayon et copiez cette couleur dans l'image.

indications : hit.object_id contient l'indice du triangle, et Mesh::triangle_material( id ) renvoie la matière du triangle d'indice id.

exercice 3 : caméra

Utilisez Orbiter::read_orbiter() pour charger la description d'un point de vue de la scène. Générez les rayons en fonction des différentes transformations définies par Orbiter.

Utilisez shader_kit pour charger votre objet, choisir un point de vue et enregistrer l'orbiter (appuyez sur C pour copier/relire un fichier, sur V pour coller / enregistrer un fichier).

./bin/shader_kit data/shaders/mesh.glsl mesh.obj

exercice 4 : ombres et éclairage direct

Choisissez un point sur chaque source de lumière et vérifiez que chaque point visible par la caméra est également visible par au moins une des sources de lumière.
Comment déterminer l'origine et la direction des rayons permettant de faire ce test ?
Les points non éclairés seront noirs.

exercice 5 : pénombres et éclairage direct

utilisez des points générés aléatoirement à la surface des sources de lumières pour évaluer l'éclairage direct.

rappel : $L_{r} (p, o) = L_{e} (p, o) + \int L_{e} (s, p) f_{r} (s, p, o) V (s, p) \frac{cos θ cos θ_{s}}{d^{2} (s, p)} d s$

$L_r(p, o) \approx L_e(p,o) + \frac{1}{N} \sum_{k=1}^{k=N} L_e(s_k, p) f_r(s_k, p, o) V(s_k, p) \frac{\cos \theta \cos \theta_{sk}}{d^2(s_k,p)} \frac{1}{pdf(s_k)}$

$L_r(p, o)= L_e(p, o) + \int L_e(s, p) f_r(s, p, o) V(s, p) \frac{\cos \theta \, \cos \theta_s}{d^2(x, y)} ds$

theta est l'angle entre la normale en p et la direction p vers s,
theta_s est l'angle entre la normale en s et la direction s vers p.

cf GI compendium eq 18 pour générer des points dans un triangle,
ou cette solution plus récente (et mieux expliquée) : "A Low-Distortion Map Between Triangle and Square", E. Heitz, 2019

exercice bonus : tonemapper le résultat

au lieu d'exporter directement la couleur résultat des calculs, appliquez une transformation à la couleur, une compression gamma 2.2, par exemple :

color= color^\frac{1}{\gamma}

pourquoi ?

la courbe rouge correspond à la conversion directe de l'énergie calculée par pixel en couleur. mais un écran standard applique une transformation aux couleurs avant de les afficher (pour raisons historiques... ça date des écrans crt...)

la courbe verte correspond à la compression (inverse) appliquée par l'écran, gamma = 2.2. en pratique, cette transformation étale les tons sombres et compresse les tons clairs.

exercice 6 au choix : toutes les sources

Comment répartir N points à la surface de l'ensemble des sources de lumières (au lieu de N par source, solution a priori utilisee à l'exercice précédent) ?

L'idée est de choisir aléatoirement d'abord une source, puis un point sur cette source, et de recommencer N fois.
Proposez au moins 2 manières différentes de choisir une source parmi l'ensemble de sources.

Vous pouvez tester vos idées avec la scène emission.obj, prévue pour ce cas précis.
Quel comportement observez-vous dans chaque cas ?

voila un exemple, avec 8 échantillons par pixel : à gauche une solution correcte. à droite, sélection uniforme de la source.

quel est le problème ? quelle est la répartition des calculs dans les 2 cas (quelles sont les sources sélectionnées dans les 2 cas pour N échantillons ?)

(en orange : les sources de lumière, les 2 panneaux lumineux emettent le même flux)

exercice 7 au choix : tous les rebonds / éclairement global

préparation : avant de vous lancer dans la version "complète", calculez une version simplifiée de l'éclairage global : l'occultation ambiante

$L_{r} (p, o) = \int \frac{1}{π} V (p, \overset{⇀}{v}) cos θ d v$

$L_r(p, o)= \int \frac{1}{\pi} V(p, \vec{v}) \cos \theta dv$ $L_r(p, o) \approx \frac{1}{N} \sum_{k=1}^{k= N} \frac{1}{\pi} V(p, \vec{v_k}) \cos \theta_k \frac{1}{pdf(\vec{v_k})}$
utilisez des directions v uniformes, cf GI compendium eq 34 ou des directions distribuées selon $\frac{\cos \theta}{\pi}$ , cf GI compendium eq 35

remarque : les directions alétoires sont générées dans un repère local, il faut donc les transformer pour connaitre leurs coordonnées dans le repère du monde. On peut construire le changement de repère avec la normale du point p et un autre vecteur (+ 2 produits vectoriels pour construire 2 directions orthogonales). Il y a quelques subtilités dans cette construction, autant en utiliser une qui est robuste et rapide :

cf "Generating a consistently oriented tangent space", Frisvad, 2012,
et "Building an Orthonormal Basis, Revisited", Pixar, 2017

struct World
{
    World( const Vector& _n ) : n(_n)
    {
        float sign= std::copysign(1.0f, n.z);
        float a= -1.0f / (sign + n.z);
        float d= n.x * n.y * a;
        t= Vector(1.0f + sign * n.x * n.x * a, sign * d, -sign * n.x);
        b= Vector(d, sign + n.y * n.y * a, -n.y);
    }

    // transforme le vecteur du repere local vers le repere du monde
    Vector operator( ) ( const Vector& local ) const { return local.x * t + local.y * b + local.z * n; }

    // transforme le vecteur du repere du monde vers le repere local
    Vector inverse( const Vector& global ) const { return Vector(dot(global, t), dot(global, b), dot(global, n)); }

    Vector t;
    Vector b;
    Vector n;
};

voici quelques résultats pour 1, 4, 16 et 64 directions :

pour vérifier que les directions générées sont correctement distribuées, vous pouvez utiliser directions.cpp, il suffit de remplacer les fonctions direction() et pdf() par les votres et d'inspecter le résultat :

make directions config=release64 bin/directions bin/image_viewer sphere.hdr density.hdr

à gauche : la densité de directions générées,
à droite l'évaluation de la densité générée par rapport à la densité souhaitée, l'image de la sphère doit être constante, au bruit près.

exemple de problème (génération uniforme, densité souhaitée cosinus) :

Maintenant que vous savez générer des directions bien distribuées (et dans le repère de la scène), il ne "reste" plus qu'à évaluer ça :

$L_r(p, o) = L_e(p, o) + L_{direct}(p, o) + L_{indirect}(p, o)$ $L_{direct}(p, o) = \int L_e(s, p) V(p, s) f_r(s, p, o) \frac{\cos \theta \, \cos \theta_s}{d^2(p, s)} ds$ $L_{indirect}(p, o) = \int L_{direct}(q, p) f_r(q, p, o) \cos\theta dv \text{ avec } q = hit(p, \vec{v})$

Commencez par restructurer un peu votre code, pour créer les fonctions direct( ) et indirect( ). q est le point visible par p dans la direction v.

Ecrivez l'estimateur Monte Carlo de L_indirect(), quelle pdf peut-on utiliser pour l'évaluer ?

Comment choisir le nombre d'échantillons utilisé pour évaluer L_direct() et L_indirect(), afin de contrôler le temps d'exécution total de l'image ?

exemple de résultats pour 1, 4, 16 et 64 directions :

bonus : le bruit, c'est moche...

colonne de gauche : calculs classiques avec 16 directions,
colonne de droite : mêmes calculs, mais le bruit est différent, l'erreur est la même... elle est répartie de manière moins visible / génante.

cf "Distributing Monte Carlo Errors as a Blue Noise in Screen Space by Permuting Pixel Seeds Between Frames", Heitz, Belcour, 2019

il est aussi possible de filtrer le bruit avec des gaussiennes / "filtres classiques" pour éliminer une partie du bruit :

cf "Spatiotemporal Variance-Guided Filtering: Real-Time Reconstruction for Path-Traced Global Illumination", C. Schied, 2017,
et sa suite "Gradient Estimation for Real-Time Adaptive Temporal Filtering", C.Schied 2018
(cette méthode est utilisée dans q2vkpt, quake 2 version illumination globale)
cf "Blockwise Multi-Order Feature Regression for Real-Time Path Tracing Reconstruction", M. Kokseka, 2019

Annexe : code de départ

cf tuto_ray.cpp

int main( const int argc, const char **argv )

{

    const char *mesh_filename=
      "cornell.obj";

    const char *orbiter_filename=
      "cornell_orbiter.txt";

    if(argc > 1) mesh_filename=
      argv[1];

    if(argc > 2) orbiter_filename=
      argv[2];

    printf("%s: '%s' '%s'\n", argv[0],
      mesh_filename, orbiter_filename);

// creer l'image resultat

Image image(1024, 640);

// charger un objet

    Mesh mesh= read_mesh(mesh_filename);

if(mesh.triangle_count() == 0)

        // erreur de
      chargement, pas de triangles

return 1;

    // creer l'ensemble de triangles /
      structure acceleratrice

BVH bvh(mesh);

Sources sources(mesh);

// charger la camera

Orbiter camera;

   
      if(camera.read_orbiter(orbiter_filename))

        // erreur,
      pas de camera

return 1;

    // recupere les transformations
      view, projection et viewport pour generer les rayons

Transform model= Identity();

Transform view= camera.view();

    Transform projection=
      camera.projection(image.width(), image.height(), 45);

    Transform viewport=
      Viewport(image.width(), image.height());

    auto cpu_start=
      std::chrono::high_resolution_clock::now();

    // parcourir tous les pixels de
      l'image

    // en parallele avec openMP, un
      thread par bloc de 16 lignes

#pragma omp parallel for schedule(dynamic, 1)

    for(int py= 0; py <
      image.height(); py++)

{

        // nombres
      aleatoires, version c++11

       
      std::random_device seed;

        // un
      generateur par thread... pas de synchronisation

       
      std::default_random_engine rng(seed());

        // nombres
      aleatoires entre 0 et 1

       
      std::uniform_real_distribution<float> u01(0.f, 1.f);

        for(int px=
      0; px < image.width(); px++)

{

           


















      Color color= Black();

           


















      // generer le rayon pour le pixel (x, y)

           


















      float x= px + u01(rng);

           


















      float y= py + u01(rng);

           


















      Point o= { ... } // origine dans l'image

           


















      Point e= { ... } // extremite dans l'image

           


















      Ray ray(o, e);

           


















      // calculer les intersections

           


















      if(Hit hit= bvh.intersect(ray))

               


















      const TriangleData& triangle=
      mesh.triangle(hit.triangle_id);          


















      // recuperer le triangle

               


















      const Material& material=
      mesh.triangle_material(hit.triangle_id);     
      // et sa matiere

               


















      Point p= point(hit,
      ray);              


















      // point d'intersection

               


















      Vector pn= normal(hit,
      triangle);       // normale
      interpolee du triangle au point d'intersection

               


















      // retourne la normale pour faire face a la camera / origine du
      rayon...

               


















      if(dot(pn, ray.d) > 0)

               


















      // accumuler la couleur de l'echantillon

               


















      float cos_theta= std::max(0.f, dot(pn, normalize(-ray.d)));

               


















      color= color + 1.f / float(M_PI) * material.diffuse * cos_theta;

           


















      image(px, py)= Color(color, 1);

}

    auto cpu_stop=
      std::chrono::high_resolution_clock::now();

    int cpu_time=
      std::chrono::duration_cast<std::chrono::milliseconds>(cpu_stop
      - cpu_start).count();

    printf("cpu  %ds %03dms\n",
      int(cpu_time / 1000), int(cpu_time % 1000));

// enregistrer l'image resultat

write_image(image, "render.png");

    write_image_hdr(image,
      "render.hdr");

return 0;

}

TP2 - lancer de rayon et Monte Carlo