M2 - Images


TP5 - lancer de rayons
et
intégration numérique




Utilisez les scenes de test suivantes, elles sont composées de très peu de triangles, ce qui permet, dans un premier temps, de se concentrer sur les principes du lancer de rayons et la réflexion de la lumière. Dans un second temps, il faudra construire un arbre pour aller plus vite et traiter des objets plus complexes...

Partie 1 : visualisation

exercice 1 : visualisation par lancer de rayon

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 le code de départ) et ne conservez que l'intersection valide la plus proche de l'origine du rayon.

Il faudra déplacer les objets pour qu'ils soient visibles par la camera, comment ?


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 / le plus proche de la camera
        pour chaque objet :
            transformer le rayon dans le repere local de l'objet

            pour chaque triangle :
                calculer l'intersection du rayon et du 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

un code de départ est fourni en annexe. relisez la section premake et les projets dans la doc pour créer un nouveau projet avec gKit.

exercice 2 : matière et couleur

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. utilisez diffuse_color.




indication :
vous pouvez utiliser shader_kit pour choisir un point de vue, l'enregistrer dans un fichier texte (appuyez sur C pour copier/relire un fichier, sur V pour coller / enregistrer un fichier) et le recharger ensuite dans votre programme, cf Orbiter::read_orbiter( ) :
./bin/shader_kit data/shaders/mesh.glsl mesh.obj


Partie 2 : occultation ambiante

remarque : utilisez firefox pour visualiser correctement les équations... chrome/chromium ne gèrent pas MathML.


L'idée est de calculer la luminance réfléchie par un point associé à une matière diffuse éclairée par un ciel uniforme, pour toutes les directions qui peuvent éclairer ce point :
I=1πL(x,v)V(x,v)cosθdvI= \frac{1}{\pi} \int L(x, v) \, V(x, v) \, \cos \theta \, dvV(x, v) vaut 1 si le ciel est visible dans la direction v, et 0 sinon.
dans le cas d'un ciel uniforme, L(x, v) est une constante arbitraire, 1 par exemple.

Comment estimer la valeur de la somme ? Il suffit de choisir N directions uniformes et de calculer la moyenne :
I1N1πV(x,vi)cosθiI \approx \frac{1}{N} \frac{1}{\pi}\sum V(x, v_i) \, \cos \theta_iavec vi, la ième direction et cos thetai le cosinus de l'angle entre la normale du point p et la ième direction.
pour déterminer la valeur de V(x, vi), il faut construire un rayon (origine x, direction vi) et calculer ses intersections avec les objets de la scène.


exercice 1 : spirale de Fibonacci

Comment choisir N directions uniformes ?
la solution la plus simple est d'utiliser la spirale de Fibonacci :
cosθi=1-2i+12N,ϕi=2π[iΦ]2π(iΦ-iΦ),Φ=(5+1)2cos \theta_i= 1 - \frac{2i +1}{2N}, \phi_i= 2 \pi \left[\frac{i}{\Phi}\right] \equiv 2 \pi \left( \frac{i}{\Phi} - \left\lfloor \frac{i}{\Phi} \right\rfloor \right), \Phi= \frac{(\sqrt 5 + 1)}{2}
et connaissant les angles, il suffit de retrouver les coordonnées de vi :

vi(x,y,z)=(cosϕisinθi,sinϕisinθi,cosθi),sinθi=1-cos2θiv_i(x, y, z)= (\cos \phi_i \sin \theta_i, \sin \phi_i \sin \theta_i, \cos \theta_i), \sin \theta_i= \sqrt{1 - \cos^2 \theta_i}
reste une dernière étape : les directions sont construites dans un repère arbitraire, il faut les connaître dans le repère du monde (pour générer les rayons et faire les tests d'intersections). l'axe Z de la spirale correspond à la normale du point, il ne reste plus qu'à "trouver" ou construire 2 autres vecteurs orthogonaux... une solution consiste à créer aléatoirement un vecteur (qui ne devrait pas être colinaire à la normale). ensuite, 2 produits vectoriels permettent de construire 2 autres vecteurs orthogonaux.

une construction très élégante et robuste est détaillée dans

#include <cmath>

// b1, b2, n sont 3 axes orthonormes.
void branchlessONB(const Vector &n,
Vector &b1, Vector &b2)
{
    float sign = std::copysign(1.0f, n.z);
    const float a = -1.0f / (sign + n.z);
    const float b = n.x * n.y * a;
    b1 = Vector(1.0f + sign * n.x * n.x * a, sign * b, -sign * n.x);
    b2 = Vector(b, sign + n.y * n.y * a, -n.y);
}
connaissant les 3 axes dans le repère du monde, comment calculer les coordonnées de la direction vi (dans le repère du monde) ?


exercice 2 : directions aléatoires

comparez les résultats obtenus avec la spirale de Fibonacci, la version perturbée de la spirale et les directions uniformes, cf GI Compendium eq 34 / cours.
Faut-il modifier l'estimation du résultat ?

mêmes questions avec des directions dont la densite (pdf) est cosθπ\propto \frac{cos \theta}{\pi} ? cf GI Compendium


exercice 3 : comparaison avec le tp4

comment formuler la stratégie du tp4 / partie 3 / exercice 3 pour calculer l'occultation ambiante ?
comparez avec vos résultats précédents.




(scènes de test avec 256 directions sur la spirale de Fibonacci perturbée)

Partie 3 : structures accélératrices

La fonction intersect() calcule l'intersection d'un rayon avec chaque triangle de la scène, ce qui est très long. Il est possible d'organiser les triangles (en fonction de leur position dans l'espace) de manière à limiter le nombre de tests d'intersections. L'idée est de grouper les triangles proches et de ne tester que les groupes qui peuvent intercepter le rayon. En recommençant cette opération sur les groupes de triangles, jusqu'à n'obtenir qu'un seul groupe, on construit un arbre.

Une solution simple pour construire cet arbre fonctionne dans l'autre sens : connaissant l'ensemble des triangles de la scène et leur englobant (une boite alignée sur les axes fonctionne très bien) comment répartir les triangles en 2 sous ensembles, afin de construire un arbre binaire ? La encore, une solution simple consiste à déterminer l'axe le plus étiré de la boite englobante et son centre, puis à répartir les triangles par rapport au plan passant par le centre de l'englobant. Tous les triangles qui se trouvent avant le plan forment le premier groupe (et le fils gauche du noeud), les autres triangles forment le 2ieme groupe (et le fils droit du noeud). Cet algorithme est appliqué récursivement à chaque fils jusqu'à n'obtenir qu'un seul triangle. Cet algorithme est une adaptation du quick sort couramment utilisé pour trier des réels (choisir un pivot, répartir les éléments se trouvant avant et après le pivot, recommencer).

Ce type d'arbre d'englobants s'appelle un BVH (Bounding Volume Hierarchy) et les BVHs sont la solution la plus utilisée (et la plus efficace) actuellement. Un noeud interne de l'arbre stocke un englobant (les points extremes de la boite englobante) ainsi que les indices de ses 2 fils. Une feuille stocke l'indice d'un triangle. L'arbre est donc représenté par un tableau (cf std::vector) de noeuds, un tableau de triangles et l'indice de la racine. Pour simplifier le code, on utilise la même structure pour les noeuds et les feuilles.

détails pratiques :
comment déterminer la position d'un triangle par rapport au plan de séparation ? un triangle peut être entièrement à gauche, entièrement à droite, ou à cheval... le plus simple est de ne considérer qu'un seul point (en général, le centre de la boite englobante du triangle) pour faire ce test. si le centre est avant la séparation sur l'axe choisi, le triangle sera affecté au fils gauche. au final, l'algorithme de construction ne tient compte que des boites englobantes des triangles, on peut l'utiliser avec n'importe quel objet / primitive.

L'algorithme de construction de partition comporte plusieurs cas particuliers, pour gagner du temps, il est recommandé d'utiliser la version fournie par la STL, cf std::partition. mais il faut définir le prédicat permettant de comparer la position du centre d'une primitive et le pivot. c'est une structure définissant l'opérateur ( ).

voici un résumé de l'algorithme de construction :

#include <algorithm>

struct BBox { ... };
struct Node { ... };

struct Primitive
{
    BBox bounds;
    Point center;
    int triangle_id;
};

int build_node( std::vector<Node>& nodes, std::vector<Primitive>& primitives, const int begin, const int end )
{
    if(end - begin <= 1)
    {
        // construire une feuille qui reference la primitive d'indice begin, et la boite englobante du triangle associee a la primitive...
        // renvoyer l'indice de la feuille
   
}

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

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

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

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

    // verifier que la partition n'est pas degeneree (toutes les primitives du meme cote de la separation)
    assert(mid != begin);
    assert(mid != end);
    // remarque : il est possible que 2 (ou plus) primitives aient le meme centre,
    // dans ce cas, la boite englobante des centres est reduite à un point, et la partition est forcement degeneree
    // une solution est de construire une feuille,
    // ou, autre solution, forcer une repartition arbitraire des primitives entre 2 les fils, avec mid= (begin + end) / 2

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

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

    // construire un noeud interne
    // quelle est sa boite englobante ?
    // renvoyer l'indice du noeud

}


et le prédicat de la partition :

struct predicat
{
    int axe;
    float coupe;

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

    bool operator() ( const Primitive& p ) const 
    {
        // renvoyer vrai si le centre de l'englobant de la primitive se trouve avant la coupe...
    }
};

exercice 1 : BVH

construisez un BVH avec les triangles d'un Mesh.

exercice 2 : parcours

écrivez l'algorithme d'intersection rayon / BVH.
vous aurez également besoin de calculer l'intersection d'un rayon et de l'englobant de chaque noeud, utilisez une solution robuste comme celle vue en cours :
// utilitaire : renvoie les coordonnees min / max de 2 points
inline Point min( const Point& a, const Point& b ) { return Point( std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z) ); }
inline Point max( const Point& a, const Point& b ) { return Point( std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z) ); }


struct
BBox
{
    Point pmin;
    Point pmax;
   
    bool intersect( const Ray& ray, const float htmax, float& rtmin, float& rtmax ) const
    {
        Vector invd= Vector(1.f / ray.d.x, 1.f / ray.d.y, 1.f / ray.d.z);
        // remarque : il est un peu plus rapide de stocker invd dans la structure Ray, ou dans l'appellant / algo de parcours, au lieu de la recalculer à chaque fois
       
        Point rmin= pmin;
        Point rmax= pmax;
        if(ray.d.x < 0) std::swap(rmin.x, rmax.x);    // le rayon entre dans la bbox par pmax et ressort par pmin, echanger...
        if(ray.d.y < 0) std::swap(rmin.y, rmax.y);
        if(ray.d.z < 0) std::swap(rmin.z, rmax.z);

        Vector dmin= (rmin - ray.o) * invd;          
// intersection avec les plans -U -V -W attachés à rmin
        Vector dmax= (rmax - ray.o) * invd;           // intersection avec les plans +U +V +W attachés à rmax
        rtmin= std::max(dmin.x, std::max(dmin.y, std::max(dmin.z, 0.f)));        // borne min de l'intervalle d'intersection
        rtmax= std::min(dmax.x, std::min(dmax.y, std::min(dmax.z, htmax)));      // borne max
       
        // ne renvoie vrai que si l'intersection est valide (l'intervalle n'est pas degenere)
        return (rtmin <= rtmax);
    }
};

exercice 3 : accélération

utilisez le BVH et l'algorithme de parcours pour calculer l'occultation ambiante d'une scène "ouverte", cf partie 4. Vous pourrez trouver plusieurs scènes au formatwavefront .obj sur la page de M. McGuire : Meshes ou dans les forums de 3drender.com. pensez à vérifier les matières dans le fichier .mtl associé, elles sont souvent "bizarres".

rappel :
vous pouvez utiliser shader_kit pour visualiser la scène, choisir un point de vue et l'enregistrer (en appuyant sur C).


(cf 3drender.com, challenge 22, TheCarnival, 256 directions sur la spirale de Fibonacci perturbée, à gauche avec des normales interpolées, à droite normales géométriques des triangles)


(cf 3drender.com, Quixel challenge ~8M triangles, 256 directions sur la spirale de Fibonacci perturbée, normales interpolées)

exercice 4 : parcours ordonné

modifiez votre algorithme de parcours pour choisir dans quel ordre visiter les fils de chaque noeud. pourquoi ?

exercice 5 : encore plus vite ?

au lieu de ne tester qu'un plan pour répartir les triangles en 2 groupes, il est possible de trouver une meilleure répartition en testant plusieurs positions et de garder la meilleure. Mais il faut d'abord définir un critère de comparaison :



Annexe : code de départ

cf ray_tuto.cpp