M2 - Images

TP 3 - Structures Accélératrices :
Hiérarchie  d'englobants


travail à rendre :

    Pour le 16/1/2009 : en plus de transmettre vos sources à votre encadrant de TP, vous rédigerez un rapport expliquant précisement comment vous détectez et éliminez les calculs d'intersection inutiles. De même, vous préciserez les gains de performances apportés par cette méthode et vous justifierez le choix des constantes Tcube et Tobjet utilisées dans le calcul l'heuristique pour la sélection du plan de séparation. Vous préciserez également les différences entre les parcours simple et efficace.
   
   
    Vous rendrez votre rapport en pdf et une archive .tar.gz de vos sources (ainsi qu'un makefile).
    Le rapport et l'archive porterons le nom des personnes associées au tp. Vous transmettrez le rapport et l'archive en utilisant le mail de l'université en précisant vos noms et numéros d'étudiants.


    message subliminal : webmail de l'université, tar.gz, makefile, merci d'avoir lu le sujet ...

Ce TP propose d'utiliser une  hierarchie de volumes englobants pour accélérer les calculs d'intersection d'un lancer de rayons.

Partie 1 : c'est un peu lent, non ?

    La visualisation de la scène est une opération longue : a priori pour chacun des Largeur*Hauteur pixels, il faut calculer l'intersection du rayon avec les N objets de la scène, soit, au total, L*H*N calculs d'intersection (sans compter les rayons supplémentaires pour les ombres et les matières transparentes).

    Une solution classique consiste à structurer la scène (organiser les objets les uns par rapport aux autres) afin d'éliminer les calculs d'intersections inutiles. L'objectif du TP est d'utiliser une hiérarchie de volumes englobants. Un autre avantage de cette structure est de gérer correctement  les scènes dynamiques (ou animéees).

    Toutes les informations relatives à la construction, au parcours et aux optimisations sont disponibles dans l'article cité en référence :
"Ray Tracing Deformable Scenes using Dynamic Bounding Volume Hierarchies".

    Hiérarchie de volumes englobants

    Une hiérarchie  de volumes englobants est un arbre dont les noeuds représentent le volume occupé par un sous-ensemble des objets de la scène. Les noeuds internes ne stockent que le volume englobant. Les feuilles, par contre, contiennent, en plus, la liste des objets présents dans leur volume englobant. Dernière propriété, le volume englobant d'un noeud est également un englobant de ses fils. Il existe un grand nombre de hiérarchies correctes pour chaque scène, mais certaines permettront d'accélérer les calculs d'intersections en éliminant les volumes non traversés par le rayon.

    Nous nous plaçons dans un cas simple : la hiérarchie sera un arbre binaire et les volumes englobants seront des cubes alignés sur le repère de la scène.

    Construction de la hiérarchie

    L'objectif est de construire la "meilleure" hiérarchie, c'est à dire, celle qui permet de faire le minimum de calculs d'intersections pour trouver l'objet le plus proche de l'origine du rayon. L'idée principale consiste à trouver une hiérarchie qui minimise le nombre de noeuds traversés par le rayon en cours de "propagation", tout en obtenant le bon résultat.

    L'algorithme général construit la hiérarchie en partant de la racine, c'est à dire le volume occupé par l'ensemble des objets de la scène. Ensuite, il faut trouver une coupe dans ce volume qui permette de répartir les objets entre les fils gauche et droit. La construction est récursive et se termine lorsqu'il n'est plus possible de séparer les objets (il n'en reste plus que 1 ou 0).

    Le problème le plus délicat est de choisir où couper un volume englobant afin de répartir les objets entre les fils gauche et droit. La solution simple consiste à répartir "équitablement" les objets entre les deux fils et permet d'obtenir un arbre binaire équilibré. Malheureusement, cette "stratégie" est une des pires pour le calcul d'intersections. En effet, un arbre équilibré ne garanti pas la meilleure séparation des objets : c'est à dire éliminer les noeuds dont le volume n'est pas traversé par le rayon.

    Il faut donc utiliser un autre critère de séparation qui tient compte du rayon. Le plus répandu utilise la probabilité qu'un rayon traverse le fils gauche (ou droit) sachant que le père est lui aussi traversé par le rayon. Cette probabilité est proportionelle au rapport de l'aire des surfaces des volumes englobants des deux noeuds : A(fils) / A(pere), avec A(noeud), l'aire de la surface du volume englobant.

    Ce critère estime le temps d'intersection du rayon avec les objets associés au noeud sachant que l'un des fils sera également intersecté. Il suffit alors de trouver la répartition entre les fils gauche et droit qui minimise le temps d'intersection T :

T = Tgauche + Tdroit

avec Tfils= Tenglobant + A(fils) / A(pere) . N(fils) Tobjet

    avec :
    On peut interpréter les différentes parties de cette expression de Tfils :
    L'algorithme de construction ne teste pas tous les plans de coupe possibles, mais uniquement les plans principaux alignés sur le repère de la scène et passant par le centre du volume englobant de chaque objet. L'algorithme, présenté en détail dans l'article de référence, parcours l'ensemble d'objets dans les deux sens (du min vers le max et réciproquement, le long de chaque axe) et construit de manière incrémentale A(gauche), N(gauche), A(droit), N(droit) pour chaque plan de coupe candidat.

    Des détails supplémentaires ainsi que l'algorithme de construction se trouvent dans l'article de référence.
    Un algorithme nettement plus efficace (et plus simple) est décrit dans "On fast Construction of SAH based Bounding Volume Hierarchies".

    Algorithme de Construction simple

étape 1 : précalculer les boites englobantes de chaque primitive ainsi que leur centre.

    Le plus plus simple est d'utiliser 3 vecteurs pour représenter l'englobant :
    struct AABOX
    {
       VEC min;
       VEC max;
       VEC centre;
    };

étape 2 : construction récursive de l'arbre.

algo : construction(S : ensemble de primitives)
// construire l'ensemble de plans de coupe candidats pour chaque axe
pour axe = X, Y, Z
P = S
trier l'ensemble P par position du centre des primitives croissant, c'est à dire P[i].centre[axe] croissant

// énumérer les plans de coupe, tous les P[i].centre[axe]
pour i = [ 0 .. N(P) [
répartir les primitives placées avant P[i].centre[axe] dans S1
répartir les primitives placées après P[i].centre[axe] dans S2

// la position d'une primitive P[j] est déterminée par la position du centre de sa boite englobante, P[j].centre
// si P[j].centre[axe] < P[i].centre[axe] alors ajouter la primitive P[j] à S1
// si P[j].centre[axe] > P[i].centre[axe] alors ajouter la primitive P[j] à S2

calculer la boite englobante de S1 et de S2
calculer l'aire des englobants des noeuds : A(S1), A(S2) et A(S)
calculer T, le coût du plan candidat

// conserver le candidat de coût minimum
si T < Tmin
Tmin = T
imin = i
axemin = axe
fin si
fin pour
fin pour

construire les fils gauche et droit avec la répartition déterminée pour Tmin,
    dans le fils gauche : ensemble S1, les primitives P[j] telles que P[j].centre[axemin] < P[imin].centre[axemin]
    dans le fils droit : ensemble S2, les autres primitives P[j] telles que P[j].centre[axemin] > P[imin].centre[axemin]

// continuer récursivement la construction des fils gauche et droit
construction(S1)
construction(S2)
fin


    Q1. Quel est le critère d'arrêt de la construction récursive ?
    Q2. Quelle est la complexité de cet algorithme ?
    Q3. Proposez une version incrémentale de l'algorithme afin de répartir efficacement les primitives avant et après P[i].centre[axe], le plan candidat.
        Quelle est la complexité de cette variante de l'algorithme ?

    Parcours "simple"

    Une fois la hiérarchie construite, il ne reste plus qu'à la parcourir afin de calculer l'intersection du rayon avec les objets qu'elle organise.
    L'algorithme effectue un parcours en profondeur de l'arbre et ne traite que les noeuds pouvant intersecter le rayon :

    si le noeud est une feuille
       renvoyer l'intersection du rayon avec les objets associés à la feuille
    fin si
    si le rayon intersecte le volume englobant du fils gauche
       descendre à gauche
    fin si
    si le rayon intersecte le volume englobant du fils droit
       descendre à droite
    fin si
   
    Si les volumes englobants des deux fils se chevauchent, ils seront testés tous les deux, dans un ordre indéterminé, ce qui est particulièrement inefficace.



    résumé : Structures accélératrices et lancer de rayons

    rappels : intersection rayon - boite englobante alignée sur les axes
    explications plus complètes, sans cas particulier : "An Efficient and Robust Ray-Box Intersection Algorithm" + code.tar.gz
   

Partie 2 : parcours "efficace"

    Nous ne sommes interressés que par l'intersection la plus proche de l'origine du rayon. Dans le cas ou les deux fils sont visités par le parcours, il est préferable de tester en premier les intersections avec le fils le plus proche de l'origine du rayon. L'autre fils ne sera testé que si le rayon traverse le premier fils sans intersection.

    Il existe plusieurs solutions, la plus simple est de commencer par le fils dont le volume englobant se projette le plus près de l'origine du rayon. La solution proposée dans l'article de référence précalcule ce résultat et l'exploite directement sans refaire de calculs, ce qui permet de gagner énormément de temps. C'est, de plus, très simple à faire :

    par construction, le fils gauche d'un noeud correspond aux primitives dont le centre est placé à gauche du plan de coupe : P[j].centre[axemin] < P[imin].centre[axemin]
    donc si la direction du rayon est positive sur l'axe de coupe, direction[axemin] > 0, il suffit de parcourir le fils gauche, puis le fils droit.
    sinon, (le rayon a une direction < 0), il suffit de parcourir le fils droit puis le fils gauche, relisez le résumé sur le parcours de rayon, si ce n'est pas clair.

    Q1. Ecrivez les fonctions nécessaires pour construire une hiérarchie de volumes englobants, ainsi que les parcours simple et efficace.
    Q2. Quels gains obtient-on par rapport à la version naïve qui intersecte la totalité des objets de la scène ?
    Q3. Il est encore possible d'éliminer certains noeuds du parcours, lesquels ?

Partie Subsidiaire

    Il est possible de construire un arbre de meilleure qualité tres simplement, la méthode est présentée dans "Early Split Clipping for Bounding Volume Hierarchies".

    Q1. modifiez la construction de votre hierarchie afin d'utiliser la modification présentée dans l'article.
    Q2. quels sont les gains apportés ?

    Q3. modifiez la construction de votre hierarchie en utilisant l'algorithme rapide "On fast Construction of SAH based Bounding Volume Hierarchies"

Annexes : Références


"Ray Tracing Deformable Scenes using Dynamic Bounding Volume Hierarchies"
Ingo Wald, Solomon Boulos, and Peter Shirley
ACM Transactions on Graphics 26(1), 2007
pdf movie(mov)


"On fast Construction of SAH based Bounding Volume Hierarchies"
Ingo Wald
Eurographics/IEEE Symposium on Interactive Ray Tracing, 2007
pdf


"Early Split Clipping for Bounding Volume Hierarchies"
M. Ernst, G. Greiner
pdf



Annexes : outils


    Finissez d'installer sdlkit, dans le repertoire objtoy se trouve une librairie supplémentaire permettant de charger des objets 3D au format maya .obj
    cd sdlkit/objtoy
    make
    make install

    scenes de test pour objtoy


    La documentation de objtoy explique comment parcourir un objet (représenté par la structure MODEL), comment parcourir l'ensemble de faces décrivant la surface de l'objet, retrouver la position des sommets d'une face, etc.    [Documentation et exemple]

    Il est plus simple de trianguler la scene après l'avoir chargée : model_set_triangles(MODEL **),

    exemple : charger un modele et le trianguler :

    #include "model.h"

    MODEL *model;
    if(model_is_obj(nom_de_fichier) && model_load_obj(&model, nom_de_fichier) < 0)
       erreur de chargement

    // triangule les faces du modele
    model_set_triangles(&model);


    exemple : parcourir toutes les faces d'un modele :

    int model_intersect( MODEL *model,
        const VERTEX origine, const VEC direction,
        float *t )
    {
       const int faces_n= model_get_faces_n(model);
       int i;
       for(i= 0; i < faces_n; i++)
       {
            // recuperer la ieme face
            const FACE *face= model_get_face_ptr(model, i);

            // recuperer les sommets de la face / triangle
            const float *v0= model_face_get_vertex_ptr(model, face, 0);
            const float *v1= model_face_get_vertex_ptr(model, face, 1);
            const float *v2= model_face_get_vertex_ptr(model, face, 2);

            // calculer l'intersection avec le triangle v0, v1, v2 et renvoyer la plus petite intersection trouvee
       }

         return -1;
    }


    Annexe : intersection d'un triangle et d'un rayon :

    "Fast, minimum storage ray-triangle intersection"
    T. Moller, B. Trumbore,
    Journal of Graphics Tools, 1997
    pdf

    version sdlkit :

    /* intersection d'un rayon et du triangle vert0, vert1, vert2
       renvoie 1 en cas d'intersection et t, l'abcisse de l'intersection,
       renvoie 0, sinon.
     */
    static
    int intersect( const VERTEX vert0, const VERTEX vert1, const VERTEX vert2,
        const VERTEX origine, const VEC direction,
        float *t )
    {
    VEC edge1, edge2, tvec, pvec, qvec;
    float det, inv_det;
    float u, v;

    /* find vectors for two edges sharing vert0 */
    vec3_sub(edge1, vert1, vert0);
    vec3_sub(edge2, vert2, vert0);

    /* begin calculating determinant - also used to calculate U parameter */
    vec3_cross(pvec, direction, edge2);

    /* if determinant is near zero, ray lies in plane of triangle */
    det= vec3_dot(edge1, pvec);
    if (det > -EPSILON && det < EPSILON)
        return 0;

    inv_det= 1.0f / det;

    /* calculate distance from vert0 to ray origin */
    vec3_sub(tvec, origine, vert0);

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

    /* prepare to test V parameter */
    vec3_cross(qvec, tvec, edge1);

    /* calculate V parameter and test bounds */
    v= vec3_dot(direction, qvec) * inv_det;
    if(v < 0.0f || u + v > 1.0f)
        return 0;

    /* calculate t, ray intersects triangle */
    *t= vec3_dot(edge2, qvec) * inv_det;
    return 1;
    }


    /* intersecte un rayon avec une face d'un modele (charge et triangule avec objtoy)
     */
    int triangle_intersect( const MODEL *model, const int face_id,
        const VERTEX origine, const VEC direction,
        float *t )
    {
    const FACE *face= model_get_face_ptr(model, face_id);
    const float *v0= model_face_get_vertex_ptr(model, face, 0);
    const float *v1= model_face_get_vertex_ptr(model, face, 1);
    const float *v2= model_face_get_vertex_ptr(model, face, 2);

    return intersect(v0, v1, v2, origine, direction, t);
    }

    Annexe : enregistrer une image HDR

    installez imgtoy dans un répertoire
    cd XXX
    make
    make install

    vous pouvez utiliser ~/local/bin/viewhdr pour visualiser les images hdr et les convertir en tga.

    exemple de création, remplissage et écriture d'une image HDR :

    #include <assert.h>
    #include "img.h"
    #include "hdr.h"

    void img_set_pix( IMG *img, const int x, const int y,
        const float r, const float v, const float b )
    {
        assert(img != NULL);
        assert(x >= 0 && x < img->largeur);
        assert(y >= 0 && y < img->hauteur);
   
        const unsigned int pix= (img->hauteur -1 - y) * img->largeur * img->channels + x * img->channels;
        img->dataf[pix]= r;
        img->dataf[pix +1]= v;
        img->dataf[pix +2]= b;
    }

    ... quelque part dans le programme :

    /* creer l'image resultat */
    IMG *img= new_img_datafloat(largeur, hauteur);

    /* ecrire des pixels dans l'image */
    int y, x;
    for(y= 0; y < hauteur; y++)
       for(x= 0; x < largeur; x++)
          img_set_pixel(img, x, y, 3.0f, 2.0f, 1.0f);

    /* sauver l'image sous image.hdr */
    hdr_write(img, "image.hdr");
    free_img(img);