M2 - Images

TP 1 - Méthodes de monte carlo

et lancer de rayons.



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 ...) :

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_ray1.cpp) et ne conservez que l'intersection valide la plus proche de l'origine du rayon.


remarques : 

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


remarque : pour déterminer la position du point d'intersection à l'abscisse t, le plus simple est d'utiliser l'opérateur ( ) de la classe Ray. par exemple, Point p= ray(hit.t);


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.


remarque : Triangle::point(0.3, 0.3) permet de trouver le centre du triangle. pourquoi ?


exercice 5 : pénombres et éclairage direct

Utilisez des directions sur l'hémisphère, générées avec une spirale de Fibonacci, ou générées aléatoirement.



ou utilisez des points générés aléatoirement à la surface des sources de lumières.




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.

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. en pratique, cette transformation étale les tons sombres et compresse les tons clairs.



exercice 6 : 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)


Partie 2 : structures accélératrices

exercice 7 : construction d'un BVH

Construisez un arbre binaire : chaque noeud contient une boite englobante alignée sur les axes et les indices de ses 2 fils. Les feuilles contiennent l'indice d'un triangle. L'arbre est donc représenté par un vecteur de noeuds et un vecteur de triangles, l'indice de la racine et sa boite englobante.

Un algorithme de construction simple consiste à calculer (récursivement en partant de la racine) pour chaque noeud :
    la boite englobante des centres des triangles,
    l'axe le plus allongé de cette boite englobante,
    la répartition des triangles par rapport au centre de la boite sur l'axe le plus étiré.

il suffit de recommencer pour chaque sous ensemble de triangles, jusqu'a obtenir 1 seul triangle.

indications :
utilisez un functor et la fonction std::partition() pour répartir les triangles entre les fils gauche et droit. il est inutile d'utiliser un tri complet...

struct BBox { ... };
struct BVHNode { ... };
struct Triangle { ... };

unsigned int build_node( std::vector<BVHNode>& nodes, std::vector<Triangle>& triangles, const unsigned int begin, const unsigned int end )
{
    if(end - begin <= 1)
    {
        // construire une feuille qui reference le triangle d'indice begin

        // renvoyer l'indice de la feuille
    }

    // construire la boite englobante des centres des triangles d'indices [begin .. end[

    // trouver l'axe le plus etire de la boite englobante
    int axe= ...

    // couper en 2 au milieu de boite englobante sur l'axe le plus etire
    float coupe= ...

    // partitionner les triangles par rapport a la "coupe"
    Triangle *pmid= std::partition(triangles.data() + begin, triangles.data() + end, predicat(axe, coupe))
    unsigned int mid= std::distance(triangles.data(), pmid);

    // construire le fils gauche
    unsigned int left= build_node(nodes, triangles, begin, mid)

    // construire le fils droit
    unsigned int right= build_node(nodes, triangles, mid, end)

    // construire le noeud

    // renvoyer l'indice du noeud
}


predicat() est un functeur c++, une classe définissant l'operateur () :

struct predicat
{
    int axe;
    float coupe;

    predicat( const int _axe, const float _coupe ) : axe(_axe), coupe(_coupe) {}

    bool operator() ( const Triangle& t ) const
    {
        // calculer le centre du triangle
        // renvoyer vrai s'il se trouve avant la coupe...
    }
};


indication : utilisez les indices des noeuds, pas de pointeurs pour la représentation des noeuds. il sera très simple de transformer le bvh en storage buffers pour la version compute shader.

exercice 8 : parcours d'un BVH

Utilisez l'arbre construit à l'exercice précédent pour trouver l'intersection valide la plus proche de l'origine d'un rayon.

exercice 9 : et avec les heuristiques SAH ?

Commencez par vous convaincre que ça fonctionne : calculez le nombre de rayons passant par chaque noeud lors d'un parcours de N rayons et vérifiez qu'il correspond au SAH du noeud (rapport de l'aire de l'englobant du noeud par l'aire de l'englobant de la scène).

Comment générer N rayons uniformes ?


Programmez la recherche exhaustive vue en cours :

ou la version rapide :


Partie 3 : réduction de variance / multiple importance sampling

Testez les 2 stratégies de base :
    choisir un point sur l'ensemble des sources,
    choisir une direction proportionnelle au cosinus.


(à gauche : statégie 1, un point uniformément réparti sur l'ensemble des sources. à droite : stratégie 2, directions proportionnelles au cosinus)

Vérifiez que les 2 techniques convergent vers la même image.
Testez l'heuristique de mélange des 2 estimateurs vue en cours.


(à gauche : mélange des 2 stratégies. à droite (en rouge) influence de la stratégie 2)


Quel type de mélange semble le plus adapté dans ce cas ? balance ? power ? cutoff ? etc.


et avec plusieurs rebonds ? combien de stratégies pour chaque longueur de chemins ?


Partie 4 : compute shaders

c'est un peu lent, non ?
lisez les tutoriels de gKit concernant les compute shaders et les storage buffers (et les storage textures, eventuellement)

exercice 1 : et avec un shader ?

Ecrivez un compute shader calculant toutes les intersections rayons / triangles, l'équivalent de la fonction intersect() du ray_tuto.cpp

Le shader prend en entrée plusieurs storage buffers avec les rayons à tester, ainsi que les triangles.
La sortie est un vecteur avec un hit par rayon. La structure hit contient au minimum : l'indice du triangle ainsi que u et v, les coordonnées barycentriques du point d'intersection.

Comment paralléliser la boucle ? est-ce que chaque thread traite les intersections d'un triangle ? ou d'un rayon ? est-il nécessaire d'utiliser de la synchronisation ?

Comment utiliser ce shader pour calculer l'éclairage direct ? Faut-il restructurer l'algorithme de rendu ?

Utilisez la mémoire partagée pour éviter de relire chaque triangle pour chaque test d'intersection. quels gains constatez vous ?
(rappel, mesure temps cpu / gpu dans les tutos)

Quelle est la meilleure organisation mémoire pour stocker l'ensemble de triangles ? de rayons ? pour la mémoire partagée ?
(rappel : accès mémoire cohérents, quelles adresses sont manipulées lorsque tous les threads d'un groupe lisent un rayon ? ou un triangle ?)

exercice 2 : préparation, parcours sans pile

Comment transformer un bvh pour pouvoir le parcourir sans utiliser de pile ?

indication : un parcours en profondeur à une structure très (très) régulière...

Ecrivez la transformation de l'arbre.
Ecrivez également la fonction de parcours de l'arbre transformé.

exercice 3 : et avec un shader ?

Ecrivez le parcours sans pile avec un compute shader.

Le shader prend en entrée plusieurs storage buffers avec les rayons à tester, ainsi que les noeuds de l'arbre, et le vecteur de triangles.
La sortie est un vecteur de hits, contenant au minimum : l'indice du triangle ainsi que u et v, les coordonnées barycentriques du point d'intersection.


indication : l'organisation mémoire du buffer qui stocke les noeuds est très (très) importante. débrouillez vous pour que chaque noeud occupe 32 octets (ou un multiple de 32 octets, si vous n'y arrivez pas), ce qui permettra de remplir complètement les caches des processeurs graphiques).

bonus : et avec une pile ?

Il est bien sur également possible d'émuler un algorithme récursif en écrivant explicitement les opérations sur une pile représentée par un storage buffer.
Comment représenter les piles de tous les rayons ? comment dimensionner ces piles ?
Quelle est la meilleure organisation mémoire de l'ensemble de piles ?
(indication : accès mémoires cohérents, quelles adresses sont manipulées lorsque tous les threads d'un groupe empilent / dépilent un noeud ?)

bonus : et avec une mini pile en mémoire partagée ?

Modifiez votre shader pour qu'il n'empile un noeud que lorsque c'est strictement nécessaire, afin de limiter le nombre de noeuds à écrire en mémoire.

Afin de limiter les lectures / écritures en mémoire video, il est aussi possible de couper la pile en 2, et de conserver les niveaux proches du sommet dans la mémoire partagée de chaque processeur graphique. Modifiez votre shader pour tester l'efficacité de cette solution.





Annexe : squelette pour démarrer

cf ray_tuto.cpp


#include <cfloat>
#include <cmath>

#include "vec.h"
#include "color.h"
#include "mat.h"
#include "mesh.h"
#include "wavefront.h"
#include "image.h"
#include "image_io.h"
#include "image_hdr.h"
#include "orbiter.h"

#define EPSILON 0.00001f

//! representation d'un rayon.
struct Ray
{
    Point o;    //!< origine.
    Vector d;    //!< direction.
    float tmax;    //!< abscisse max pour les intersections valides.
   
    Ray( const Point origine, const Point extremite ) : o(origine), d(Vector(origine, extremite)), tmax(1) {}
    Ray( const Point origine, const Vector direction ) : o(origine), d(direction), tmax(FLT_MAX) {}
   
    //!    renvoie le point a l'abscisse t sur le rayon
    Point operator( ) ( const float t ) const { return o + t * d; }
};

//! representation d'un point d'intersection.
struct Hit
{
    Point p;        //!< position.
    Vector n;        //!< normale.
    float t;        //!< t, abscisse sur le rayon.
    float u, v;        //!< u, v coordonnees barycentrique dans le triangle.
    int object_id;  //! indice du triangle dans le maillage.
   
    Hit( ) : p(), n(), t(FLT_MAX), u(0), v(0), object_id(-1) {}
};

struct Triangle : public TriangleData
{
    Triangle( ) : TriangleData() {}
    Triangle( const TriangleData& data ) : TriangleData(data) {}
   
    /* calcule l'intersection ray/triangle
        cf "fast, minimum storage ray-triangle intersection"
        http://www.graphics.cornell.edu/pubs/1997/MT97.pdf

        renvoie faux s'il n'y a pas d'intersection valide, une intersection peut exister mais peut ne pas se trouver dans l'intervalle [0 htmax] du rayon. \n
        renvoie vrai + les coordonnees barycentriques (ru, rv) du point d'intersection + sa position le long du rayon (rt). \n
        convention barycentrique : t(u, v)= (1 - u - v) * a + u * b + v * c \n
    */
    bool intersect( const Ray &ray, const float htmax, float &rt, float &ru, float&rv ) const
    {
        /* begin calculating determinant - also used to calculate U parameter */
        Vector ac= Vector(Point(a), Point(c));
        Vector pvec= cross(ray.d, ac);

        /* if determinant is near zero, ray lies in plane of triangle */
        Vector ab= Vector(Point(a), Point(b));
        float det= dot(ab, pvec);
        if(det > -EPSILON && det < EPSILON)
            return false;

        float inv_det= 1.0f / det;

        /* calculate distance from vert0 to ray origin */
        Vector tvec(Point(a), ray.o);

        /* calculate U parameter and test bounds */
        float u= dot(tvec, pvec) * inv_det;
        if(u < 0.0f || u > 1.0f)
            return false;

        /* prepare to test V parameter */
        Vector qvec= cross(tvec, ab);

        /* calculate V parameter and test bounds */
        float v= dot(ray.d, qvec) * inv_det;
        if(v < 0.0f || u + v > 1.0f)
            return false;

        /* calculate t, ray intersects triangle */
        rt= dot(ac, qvec) * inv_det;
        ru= u;
        rv= v;

        // ne renvoie vrai que si l'intersection est valide (comprise entre tmin et tmax du rayon)
        return (rt <= htmax && rt > EPSILON);
    }

    //! renvoie l'aire du triangle
    float area( ) const
    {
        return length(cross(Point(b) - Point(a), Point(c) - Point(a))) / 2.f;
    }
   
    //! renvoie un point a l'interieur du triangle connaissant ses coordonnees barycentriques.
    //! convention p(u, v)= (1 - u - v) * a + u * b + v * c
    Point point( const float u, const float v ) const
    {
        float w= 1.f - u - v;
        return Point(Vector(a) * w + Vector(b) * u + Vector(c) * v);
    }

    //! renvoie une normale a l'interieur du triangle connaissant ses coordonnees barycentriques.
    //! convention p(u, v)= (1 - u - v) * a + u * b + v * c
    Vector normal( const float u, const float v ) const
    {
        float w= 1.f - u - v;
        return Vector(na) * w + Vector(nb) * u + Vector(nc) * v;
    }
};


//! representation d'une source de lumiere.
struct Source : public Triangle
{
    Color emission;     //! flux emis.
   
    Source( ) : Triangle(), emission() {}
    Source( const TriangleData& data, const Color& color ) : Triangle(data), emission(color) {}
};

// ensemble de sources de lumieres
std::vector<Source> sources;

// recuperer les sources de lumiere du mesh : triangles associee a une matiere qui emet de la lumiere, material.emission != 0
int build_sources( const Mesh& mesh )
{
    for(int i= 0; i < mesh.triangle_count(); i++)
    {
        // recupere la matiere associee a chaque triangle de l'objet
        Material material= mesh.triangle_material(i);

        if((material.emission.r + material.emission.g + material.emission.b) > 0)
            // inserer la source de lumiere dans l'ensemble.
            sources.push_back( Source(mesh.triangle(i), material.emission) );
    }

    printf("%d sources.\n", (int) sources.size());
    return (int) sources.size();
}

// verifie que le rayon touche une source de lumiere.
bool direct( const Ray& ray )
{
    for(size_t i= 0; i < sources.size(); i++)
    {
        float t, u, v;
        if(sources[i].intersect(ray, ray.tmax, t, u, v))
            return true;
    }
   
    return false;
}


// ensemble de triangles
std::vector<Triangle> triangles;

// recuperer les triangles du mesh
int build_triangles( const Mesh &mesh )
{
    for(int i= 0; i < mesh.triangle_count(); i++)
        triangles.push_back( Triangle(mesh.triangle(i)) );
   
    printf("%d triangles.\n", (int) triangles.size());
    return (int) triangles.size();
}


// calcule l'intersection d'un rayon et de tous les triangles
bool intersect( const Ray& ray, Hit& hit )
{
    hit.t= ray.tmax;
    for(size_t i= 0; i < triangles.size(); i++)
    {
        float t, u, v;
        if(triangles[i].intersect(ray, hit.t, t, u, v))
        {
            hit.t= t;
            hit.u= u;
            hit.v= v;
           
            hit.p= ray(t);      // evalue la positon du point d'intersection sur le rayon
            hit.n= triangles[i].normal(u, v);
           
            hit.object_id= i;    // permet de retrouver toutes les infos associees au triangle
        }
    }
   
    return (hit.object_id != -1);
}

// construit un repere ortho tbn, a partir d'un seul vecteur...
// cf "generating a consistently oriented tangent space"
// http://people.compute.dtu.dk/jerf/papers/abstracts/onb.html
struct World
{
    World( const Vector& _n ) : n(_n)
    {
        if(n.z < -0.9999999f)
        {
            t= Vector(0, -1, 0);
            b= Vector(-1, 0, 0);
        }
        else
        {
            float a= 1.f / (1.f + n.z);
            float d= -n.x * n.y * a;
            t= Vector(1.f - n.x * n.x * a, d, -n.x);
            b= Vector(d, 1.f - n.y * n.y * a, -n.y);
        }
    }
   
    Vector operator( ) ( const Vector& local )  const
    {
        return local.x * t + local.y * b + local.z * n;
    }
   
    Vector t;
    Vector b;
    Vector n;
};



// objet
Mesh mesh;

int main( int argc, char **argv )
{
    // init generateur aleatoire
    srand48(time(NULL));

    // lire un maillage et ses matieres   
    mesh= read_mesh("cornell.obj");
    if(mesh == Mesh::error())
        return 1;

    // extraire les sources   
    build_sources(mesh);
    // extraire les triangles du maillage
    build_triangles(mesh);
   
    // relire une camera
    Orbiter camera;
    camera.read_orbiter("orbiter.txt");

    // placer une source de lumiere   
    Point light= camera.position();
   
    // creer l'image pour stocker le resultat
    Image image(1024, 640);
  
// multi thread avec OpenMP
#pragma omp parallel for schedule(dynamic, 16)
    for(int y= 0; y < image.height(); y++)
    for(int x= 0; x < image.width(); x++)
    {
        Point o= { ... };    // origine du rayon
        Point e= { ... };    // extremite du rayon
       
        Ray ray(o, e);
        Hit hit;
        if(intersect(ray, hit))
        {
            // calculer l'eclairage direct
            Color direct= { ... };
           
            // calculer l'eclairage indirect
            Color indirect= { ... };
           
            image(x, y)= Color(Red(), 1);
        }
    }
   
    write_image(image, "render.png");
    write_image_hdr(image, "render.hdr");

    return 0;
}