M1 synthèse d'images
2022

TP2 - lancer de rayons et images

Installation

vous pouvez faire le tp :

soit en installant la base de code minimaliste gKit3 sans dépendances : git clone https://forge.univ-lyon1.fr/JEAN-CLAUDE.IEHL/gkit3.git
soit en installant la base de code gKit2 complete, cf la doc en ligne : "installation"

la solution la plus simple est d'utiliser le code minimaliste et de cloner le depot gkit3.
si vous souhaitez utiliser la base de code complète, le plus simple est de travailler sous linux, les librairies sont deja installées sur les machines du nautibus.
si vous souhaitez travailler sur votre portable sous linux, il suffit d'installer les librairies.
sous windows, il faudra le faire à la main et il est plus aussi pratique d'utiliser le compilateur officiel de microsoft, cf visual studio community 2022 (pas visual studio code, ce n'est pas du tout la meme chose...), ou le code minimaliste...

si vous souhaitez un aperçu des fonctions nécessaires pour réaliser le tp, vous pouvez consulter la doc sur les points et les vecteurs. ainsi que celle sur les images + les couleurs.

Création d'un projet pour le tp

continuez à lire la doc sur l'installation, cf "étape 5 créer un projet".

avec gkit2light, la base de code complete :
pour démarrer, le plus rapide : créez un fichier tp.cpp dans le répertoire gkit2light/projets, et modifiez le fichier premake4.lua comme indiqué dans la doc.

// exemple de code test a compiler pour verifier que tout est ok
#include "color.h"
#include "image.h"
#include "image_io.h"

int main( )
{
    // cree l'image resultat
    Image image(1024, 1024);

    for(int py= 0; py < image.height(); py++)
    for(int px= 0; px < image.width(); px++)
        image(px, py)= Red();

    write_image(image, "image.png"); // par defaut en .png
    return 0;
}

générez les makefiles / projets et vérifiez que ce code d'exemple compile et s'exécute, repassez dans "étape 3 générer les projets" si ce n'est pas clair...

l'exécutable est placé dans bin/, pour tester :./bin/tp

avec gkit3, le code minimaliste :
créez un fichier dans le repertoire gkit3/projets et modifiez le fichier premake4.lua, ajoutez votre tp à la liste de projets (tp1 et tp2 sont deja presents).

générez les makefiles / projets et vérifiez que ce code d'exemple compile et s'exécute, repassez dans "étape 3 générer les projets" si ce n'est pas clair...
l'exécutable est placé dans bin/, pour tester :./bin/tp

Partie 1 : et en 2D ?

exercice 1 :
on va commencer par un cas simple, en 2D, directement dans l'image. dessinez un triangle abc.
ecrivez le test d'orientation en calculant un déterminant, par exemple.
testez tous les pixels de l'image. Vous pouvez colorier les pixels du triangle en blanc, par exemple, cf White() dans color.h

rappel : les x manières de calculer l'orientation d'un triangle / de 2 vecteurs :
cf "Optimizing Ray-Triangle Intersection via Automated Search", A.Kensler, P. Shirley,
et un rappel sur le calcul du determinant d'une matrice 2x2 et 3x3, cf wikipedia

une relation utile : si calculer le determinant d'une matrice 2x2 ne vous semble pas tres intuitif

pour un vecteur v de coordonnées (v_x, v_y), on peut construire un vecteur perpendiculaire avec les coordonnées (-v_y, v_x) et
calculer un produit scalaire en 2D : dot(u, v)= u_x*v_x+ u_y*v_y
ce qui permet de calculer l'aire (signée) d'un triangle et de faire un test d'orientation...

exercice 2 :
c'est quand meme dommage de tester l'inclusion de tous les pixels de l'image... construisez le rectangle englobant du triangle.
ne testez que ses pixels.

exercice 3 :
definissez une couleur associée à chaque sommet.
calculez les coordonnées barycentriques du pixel et utilisez l'interpolation barycentrique pour déterminer la couleur du pixel à partir des couleurs des sommets.

bonus, exercice 4 :
et avec 2 triangles, de couleur différentes, en partie superposés ?
comment choisir la couleur d'un pixel sur lequel se dessinent les 2 triangles ?

remarque : on peut aussi melanger les 2 couleurs...

Partie 2 : version 3D

exercice 1 : camera
générez le vecteur op / le rayon passant par l'origine de la camera et le pixel (px, py) dans le plan image.

votre programme devrait ressembler à quelque chose comme ça :

#include "color.h"
#include "image.h"
#include "image_io.h"

int main( )
{

    // cree
        l'image resultat

Image image(1024, 1024);

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

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

{
// rayon

   
            float x= ...;

float y= ...;
float z= -1;

        Point p= Point(x, y, z);

        Point o=
          Point(0, 0, 0);

                Vector d= Vector(o,
        p);    // ou Vector d= p - o; // si vous
        preferrez...

        // triangle
        Point a= ...;
        Point b= ...;
        Point c= ...;

        Vector n= normalize(
        cross( Vector(a, b), Vector(a, c) ) );

        image(px, py)=
        White();

}

   
        write_image(image, "image.png"); // par defaut en .png

return 0;
}

exercice 2 : test d'inclusion
ecrivez le test d'inclusion direct, vérifiez simplement que le vecteur op est du bon cote des 3 aretes...

attention à l'orientation des faces du tetraedre !
pour valider/debugguer le test, utilisez une configuration connue et verifiez que les signes des volumes des 4 tetraedres oabc, oabp, obcp et ocap sont identiques.

relisez la section 3.2 de "3D Rasterization", pour les notations précises.
les auteurs proposent de numéroter les sommets abc d'un triangle, ce qui simplifie les notations, par exemple, pour calculer les normales des faces du tetraedre pour le test d'orientation :

p₀= a, p₁= b, p₂= c

n_i= cross( p_(i+2)%3 , p_(i+1)%3 ), ou n_i= cross( p_(i+2)%3 - o, p_(i+1)%3 - p_(i+2)%3 ) ou n'importe quelle autre combinaison qui preserve l'orientation...
V_i(u)= dot( n_i, u )

lambda_i= V_i(d) / ( V₀(d) + V₁(d) + V₂(d) )

et le test d'inclusion de la direction op pour le pixel p s'écrit directement, il suffit d'evaluer les 3 V_i( op ). s'ils sont tous les 3 positifs, la droite op passe par le triangle et le pixel p est dans la projection le triangle dans l'image...

pour le triangle de l'exemple ci-dessus, on peut également afficher le signe de V_i pour chaque face du tetradre :

le signe de chaque volume, positif en blanc, negatif en rouge...

exercice 3 : et les coordonnées barycentriques...
évaluez lambda_i pour chaque sommet.
même question que précédemment, associez une couleur à chaque sommet et interpolez-la pour dessiner le triangle.

exercice 4 : et le Z ?
utilisez l'interpolation pour calculer la coordonnée z pour un pixel à l'intérieur du triangle.

exercice 5 : et maintenant avec plusieurs triangles !
meme question que precedemment, mais utilisez un ZBuffer pour conserver la couleur du triangle le plus proche.

attention ! le Z interpolé à la question précédente est la coordonnée dans l'espace camera, pas une distance / profondeur et comme la camera regarde dans la direction -Z, c'est le point avec la plus grande coordonnée Z qui est le proche de la camera...

data/robot.obj, 3200 triangles, image 1024x1024, test sur tous les pixels de l'image, 2s5
ne dessiner que les triangles bien orientes est plus rapide : 1s2...

pour les curieux : on peut charger les positions des sommets des triangles d'un objet 3d au format .obj / wavefront avec :

    #include "mesh_io.h"

    const char *filename= "data/robot.obj";         std::vector<Point> positions;
    if(!read_positions(filename, positions))
        return "erreur";
    printf("%d triangles\n", int(positions.size() / 3));

    // deplace tous les sommets devant la camera
    for(unsigned i= 0; i < positions.size(); i++)
        positions[i]= positions[i] + Vector(0, -2, -4);     // a ajuster en fonction de l'objet...

    // englobant des points, verifier qu'ils sont bien devant la camera...
    Point pmin= positions[0];
    Point pmax= positions[0];
    for(unsigned i= 1; i < positions.size(); i++)
    {
        pmin= min(pmin, positions[i]);
        pmax= max(pmax, positions[i]);
    }
    printf("bounds [%f %f %f]x[%f %f %f]\n", pmin.x, pmin.y, pmin.z, pmax.x, pmax.y, pmax.z);

    // parcours tous les triangles     for(unsigned i= 0; i +2 < positions.size(); i+= 3)
    {
        Point p[3]= {
            positions[ i    ],
            positions[ i +1 ],
            positions[ i +2 ]
        };

        ...
    }

exercice 6 : ne pas dessiner les triangles mal orientés...
c'est quand meme dommage de dessiner les triangles qui sont à l'arrière de l'objet, vu que la surface de l'objet est fermée (en général), ils ne sont jamais visibles dans l'image ! comment détecter l'orientation du triangle dans l'image ? et ne dessiner que les triangles orientés vers la camera ?

pour les curieux :
c'est quand meme dommage de tester tous les pixels de l'image, non ? peut on calculer le rectangle englobant de la projection du triangle dans tous les cas ?
indication : comment se projette un triangle dont au moins un sommet se trouve derriere la camera ?

au lieu de tester directement les pixels, on commence par tester l'englobant d'un bloc de 32x32 pixels (par exemple). comment déterminer rapidement (en ne testant que les 4 coins...) qu'un bloc est à l'exterieur du triangle ? quelle est la taille de bloc la plus efficace ?

l'exemple precedent avec les 3200 triangles de data/robot.obj est dessine un peu plus vite : < 0.1s... mais le test sur les blocs fait au moins 4 fois trop de tests... on peut aller nettement plus vite en reflechissant un peu plus pour eviter de refaire plusieurs fois le meme calcul. on peut gagner encore du temps en se rendant compte que l'on peut aussi interpoler les volumes entre 2 pixels voisins ! ie si on calcule les 3 volumes sur un pixel, on peut evaluer les volumes du pixel voisin sans tout recalculer...

Partie 3 : j'aime pas les tetraedres...

l'objectif de cette partie est de construire une fonction d'intersection rayon / triangle en re-utilisant les calculs réalisés aux questions précédentes.

on peut calculer l'intersection en 2d (dans un plan en tout cas) :

commencer par calculer t, la position du point d'intersection entre le plan qui porte le triangle et le rayon. calculez les coordonnées de ce point d'intersection p(t).
rappel : le triangle abc se trouve dans le plan qui passe par a et dont la normale est celle du triangle, n= cross( ab, ac ), par exemple.

il ne reste plus qu'à évaluer les coordonnées barycentriques, c'est à dire les aires normalisées des triangles abp, bcp, cap
rappel : on est quand meme en 3d... donc orientation(abp)= dot(n, cross(ab,ap) )

on peut également utiliser les calculs de la version 3d. mais comment déterminer t, la position de l'intersection sur le rayon ?

indication : quelle est la relation entre V₀(d) et V₀(td) ? ou pourquoi les lambda sont-il corrects dans les calculs précédents ?
indication : quelle est la relation entre V₀(d) et V₀(op₀) ?
indication : une propriete du produit scalaire permet de repondre : V(td) == dot( td, n ) == t dot( d, n ) == t V(d)...

dans les 2 cas, peut-on éviter de faire certains calculs pour gagner un peu de temps ?
(on va calculer des millions d'intersections rayon / triangle pour construire une image, voir beaucoup plus en ajoutant les ombres, les reflets, etc)
les aires ou les volumes utilises pour faire les tests d'orientation sont-ils tous necessaires ?

indication : par construction, les coordonnées barycentriques d'un point dans le triangle sont entre 0 et 1 et leur somme est égale à 1...

vérfiez que la fonction d'intersection rayon / triangle dessine le meme triangle que dans la partie 1...

pour les rapides : chargez un objet et dessinez- le !!

pour les curieux : une solution alternative

lisez la section 2.1 "Ray Triangle Intersection" de l'article "Realtime ray tracing of dynamic scenes on an FPGA chip".

exercice 1 :
construisez la transformation du triangle vers le triangle unitaire, cf section 2.1.1
vérifiez que le triangle transformé est bien unitaire... et que le rayon est bien transforme lui aussi.

est-il nécessaire de stocker les sommets du triangle ? ou : que faut -ilstocker comme informations pour representer un triangle ? un rayon ?

indication : cf les differents constructeurs de Transform dans mat.h

exercice 2 :
ecrivez la fonction d'intersection, cf section 2.1.2

quel est le resultat de cette fonction d'intersection ?

comment obtenir la position du point d'intersection ?
comment calculer les coordonnées baycentriques ? et interpoler des attributs de sommets ?

exercice 3 :
meme question que précédemment, associez 3 couleurs aux sommets et interpolez-les pour dessiner le triangle.

exercice 4 :
et le ZBuffer ? est-il necessaire d'utiliser un ZBuffer pour obtenir le triangle le plus proche de la camera pour chaque pixel, ou peut-on structurer le code differement ?

Partie 4 : une image complète (à peu près...)

maintenant que vous avez une fonction d'intersection rayon / triangle qui fonctionne, il ne reste plus qu'à construire une premiere image.
le principe est toujours le meme, pour chaque pixel de l'image, on construit la direction op, dit le rayon, on determine quel est l'objet visible et il ne reste plus qu'à calculer une couleur pour le pixel.

pour charger les triangles, utilisez l'extrait de code de la partie précédente pour charger geometry.obj et materials.mtl
il y a peu de triangles dans cette scene, et les calculs seront plus rapides que sur le robot et ses 3000 triangles. il faudra construire un BVH pour accélérer les calculs d'intersections, mais ce sera pour plus tard.

les fonctions read_positions(), read_materails(), etc et leurs docs sont dans src/mesh_io.h

ecrivez une fonction utilitaire qui calcule les intersections avec un ensemble de triangles et renvoie la plus proche, ce sera plus pratique pour la suite :

Hit intersect( const Point& o, const Vector& d, const std::vector<Point>& positions )

cette fonction peut aussi renvoyer la position du point d'intersection et eventuellement la normale du triangle, si l'intersection existe.
par exemple :

struct Hit
{
    float t;
    float u,v;
    int triangle_id;    // indice du triangle, pour retrouver la matiere associee au triangle, par exemple

    Point p;            // p= o + td
    Vector n;           // normale du triangle
};

on peut mettre -1 dans triangle_id ou t pour indiquer qu'il n'y a pas d'intersection. vous pouvez aussi ajouter un bool, tout simplement.

exercice 1 :
comment calculer une couleur pour chaque pixel ? au minimum, il faut reproduire le comportement d'une matière diffuse, ie la matiere réfléchit la même quantité de lumière dans toutes les directions. reste à calculer quelle quantité de lumière éclaire le point d'intersection à la surface de l'objet, et cette quantité, la lumière incidente, dépend de l'orientation de surface.

la lumière incidente varie en fonction de cos theta, le cosinus de l'angle entre la normale et la direction vers la source.
vous pouvez recalculer la normale du triangle sur lequel se trouve le point d'intersection.

rappel : pour évaluer le cosinus on utilise une relation avec le produit scalaire :

cos theta= dot( normalize(n),
      normalize(l) )

exercice 2 : et les ombres ?
le point d'intersection sera à l'ombre, si un objet se trouve entre la source de lumière et le point d'intersection, sinon il sera éclairé.

il suffit de construire un nouveau rayon entre p et la position de la source, et de calculer les intersections. s'il y a une intersection le point est à l'ombre ! et le pixel sera noir.

attention : on peut très bien trouver une intersection sur le rayon, mais plus loin que la source, ie avec t > 1, elle n'est pas valide dans ce cas, on ne s'interresse qu'aux intersections entre p et la source, ie t > 0 et t < 1

pensez également à décoller l'origine du rayon du triangle, sinon vous aurez une image pleine de défauts. cf doc