L3 synthèse d'images
2024

TP1 - lancer de rayons et images


Installation

installez la base de code minimaliste gKit3 sans dépendances : git clone https://forge.univ-lyon1.fr/JEAN-CLAUDE.IEHL/gkit3.git

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.

gkit utiise premake pour gérer ses projets, lisez les instructions sur la page du projet.

un squelette de programme est disponible dans projets/tp2.cpp, vous pouvez le modifier pour démarrer le tp.

en résumé sous linux :
git clone https://forge.univ-lyon1.fr/JEAN-CLAUDE.IEHL/gkit3.git
cd gkit3
premake5 gmake
make tp2
bin/tp2

premake5 peut également générer un CMakeList.txt pour simplifier l'intégration avec vs code, cf premake5 cmake


Partie 1 : intersections

exercice 1 : plan
on va commencer par un cas simple, mais avant de calculer une image, on va prendre le temps de vérifier que la fonction d'intersection fournit un resultat correct / utilisable.

relisez le cours et écrivez une fonction d'intersection rayon / plan :

#include "vec.h"

float intersect_plan( /* parametres du plan */ const Point& a, const Vector& n, /* parametres du rayon */ const Point& o, const Vector& d )

{
    ...
}

rappel : le produit scalaire de 2 vecteurs u et v, déclarés comme Vector u, v; s'écrit dot(u, v), cf la doc

vous pouvez renvoyer inf si l'intersection avec le plan se trouve derriere la camera / l'origine du rayon. comment peut-on détecter ce cas ?

// infini
#include <limits>

const float inf= std::numeric_limits<float>::infinity();


comment vérifier que la fonction d'intersection calcule le bon résultat ?

on peut choisir un cas simple ou la réponse est directe : par exemple, pour le rayon o= (0 0 0), d= (0 0 -1) et un plan a= (0 0 -1) et n= (0 0 1), l'intersection doit être à la position 1.
faites le test et vérifiez que tout fonctionne :

int main( )
{
    // rayon
    Point o= Point(0, 0, 0);
    Vector d= Vector(0, 0, -1);

    // plan
    Point a= Point(0, 0, -1);
    Vector n= Vector(0, 0, 1);

    float t= intersect_plan(a, n, o, d);

    std::cout << t << std::endl;
    if(t != 1) std::cout << "oops, loupe\n";
    else       std::cout << "c'est bon...\n";

    return 0;
}

rappel : les objets devant la camera ont une coordonnée Z < 0...


exercice 2 : sphère
relisez le cours et écrivez une fonction d'intersection rayon / sphère :

#include <cmath>
#include "vec.h"

float intersect_sphere( /* parametres de la sphere */ const Point& c, const float r, /* parametres du rayon */ const Point& o, const Vector& d )
{
    float a= ...;
    float b= ...;
    float k= ...;

    float det= b*b - 4*a*k;
    ...
}

attention, il y a plusieurs cas à prendre en compte, il peut y avoir 0, 1 ou 2 intersections avec la droite infinie qui porte le rayon, mais le rayon n'est qu'une partie de cette droite... la droite peut traverser la sphère, mais pas le rayon. la fonction ne doit renvoyer que les intersections avec le rayon ou inf pour indiquer qu'il n'y a pas d'intersection.

vérifiez que votre intersection est correcte dans les 3 cas lorsque la droite infinie (qui porte le rayon) traverse la sphère :

le plus simple est de placer une sphère de rayon 1 en (0 0 0) et de modifier l'origine du rayon, pour chaque cas. la direction du rayon sera (0 0 -1), par exemple.
vérifiez également que votre fonction renvoie inf lorsque la droite et le rayon passent à cote de la sphère. il y a donc au moins 4 tests à écrire.


remarques :
vous pouvez également créer une structure sphère pour regrouper tous les paramètres d'une sphère, son centre, son rayon, et sa couleur.
vous pouvez aussi créer une structure rayon et une autre pour stocker les resultats de l'intersection : la position sur le rayon, cf t, mais aussi la couleur de l'objet, sa normale, etc. ça sera super pratique pour calculer la lumière réfléchie vers le pixel / la camera :
struct Hit
{
    float t;
    Point p;
    Vector n;
    Color color;    // ou Material material; si avez une description de matière plus complète, avec des reflets par exemple...)
};
si vous vous posez plein de questions sur comment organiser tout ce code, regardez la partie "grand nettoyage du code".


Partie 2 : camera et rayons

il ne reste plus qu'à calculer l'origine et la direction des rayons pour les différents pixels de l'image et on pourra (enfin) calculer une image complete !

relisez le cours et écrivez la boucle de parcours des pixels de l'image (déjà fait dans le squelette projets/tp2), calculez le point sur le plan image qui correspond au pixel. la droite qui passe par l'origine (0 0 0) et le point sur le plan image permettent de définir le rayon : son origine (0 0 0) et sa direction.

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

int main( )
{
    // cree l'image resultat
    Image image(512, 512);    // par exemple...

    for(int py= 0; py < image.height(); py++)
    for(int px= 0; px < image.width(); px++)
    {        
        Point o= Point(0, 0, 0);    // origine
        Point e= ...                // extremite
        Vector d= Vector(o, e);     // direction : extremite - origine

        image(px, py)= Red();

    }

    write_image(image, "image.png");
    return 0;
}

la doc sur les images et la doc sur lire/ecrire une image.


Partie 3 : image !

exercice 1 :
créez une sphère de centre (0 0 -3) et de rayon 2, calculez l'image, coloriez un pixel dont le rayon touche la sphère en rouge, cf image(px, py)= Red() ou Color(1,0,0).




exercice 2 : et avec un sol ?
ajoutez un plan a= (0 -1 0), n= (0 1 0), coloriez les intersections du plan en gris, par exemple...

   

attention : lorsqu'il y a plusieurs intersections et plusieurs objets, il faut bien renvoyer l'intersection valide la plus proche de l'origine du rayon... la bonne image est celle de gauche, pas celle de droite...


pour les curieux : et si l'image n'est pas carrée ? comment corriger la déformation ? par exemple, si on calcule une image (512 256) :



l'image est 2 fois plus large que haute, et la sphère ne ressemble plus trop à une sphère...
il faut modifier le plan image, il doit lui aussi etre 2 fois plus large que haut...
au lieu de générer le point sur la largeur du plan image entre -1 et 1, il faut utiliser -2 et 2. comment ?
relisez le cours sur la construction du plan image.



est-il possible de faire ce calcul quel que soit le rapport largeur / hauteur de l'image ? comment ?


Partie 4 : matière diffuse

exercice 1 : "soleil" et cos theta
pour commencer simplement, on va considérer que le "soleil" éclaire les objets.
relisez le cours et calculez la lumière réfléchie par les objets composés d'une matière diffuse / Lambert.

il y a 2 étapes :
quelle est la normale de l'objet au point d'intersection ? pour le plan, c'est simple, c'est la même que celle du plan. et pour la sphère ?

   

rappel : relation entre cos theta et le produit scalaire des 2 vecteurs : float cos_theta= dot(normalize(n), normalize(l))

attention : l'angle entre les 2 vecteurs doit etre compris entre 0 et 90°, donc le cosinus est compris entre 0 et 1. si les vecteurs sont de sens opposés, leur produit scalaire et le cosinus sont < 0, que faut-il faire dans ce cas ? (indication : quelle est la couleur de la sphère lorsqu'elle n'est pas éclairée ?)

indications : vous pouvez reproduire l'exemple avec Vector l= normalize(Vector(-4, 6, 1));          // direction vers le "soleil"


Intermède : grand nettoyage du code !

avant d'écrire la suite, il est conseillé de passer un peu de temps à nettoyer le code et surtout à le ré-organiser pour le rendre plus simple à manipuler. il y a probablement 2 choses à reprendre, il faut connaitre facilement la normale de l'intersection ainsi que la couleur de l'objet touché par le rayon. pour la suite du tp, on va construire plus de rayons et calculer plus d'intersections avec plus d'objets, pour calculer les ombres, les pénombres et les reflets, par exemple.

il y a bien sur plusieurs manières d'organiser tout ca. voila une solution qui fontionne bien et qui permet de faire la suite du tp sans problemes :


exercice 1 : nettoyage des fonctions d'intersection
créez une structure plan, avec les parametres d'un plan et sa couleur, créez une autre structure pour la sphère.

par exemple :

struct Plan
{
    Point a;        // point sur le plan
    Vector n;       // normale du plan
    Color color;    // couleur
};

il serait plus simple que les fonctions d'intersection renvoient toutes les infos sur une intersection : sa position sur le rayon, sa normale et la couleur de l'objet. le plus direct est d'utiliser une structure :

struct Hit
{
    float t;        // position sur le rayon, ou inf s'il n'y a pas d'intersection
    Point p;        // position du point, s'il existe
    Vector n;       // normale du point d'intersection, s'il existe

    Color color;    // couleur du point d'intersection, s'il existe

    Hit( ) : t(inf), p(), n(), color() {}     // pas d'intersection
    Hit( const float _t, const Point& _p, const Vector& _n, const Color& _color ) : t(_t), p(_p), n(_n), color(_color) {}

};

et les fonctions d'intersection prennent en paramètre la structure qui stocke les paramètres de l'objet et renvoient une structure Hit avec toutes les informations sur l'intersection :

Hit intersect( const Plan& plan, const Point& o, const Vector& d )
{
    float t= ... ;
    Point p= ... ;


    if(t < 0)
        return Hit();                          // renvoie une intersection non valide / par defaut. l'intersection n'est pas valide / derriere l'origine du rayon

    else
        return Hit(t, p, plan.n, plan.color);  // renvoie la position de l'intersection, + la normale et la couleur du plan

}

remarque : oui bien sur, on peut définir une structure pour le rayon...

on peut encore se simplifier la vie en ajoutant un paramètre aux fonctions d'intersection : la valeur max de t, ie quelles sont les valeurs possibles de t : entre 0 et tmax. les intersections en dehors de cet intervalle ne sont pas valides. ca permet de simplifier encore un peu la logique du reste du code pour garder l'intersection la plus proche de la camera.

pourquoi tmax ?
les rayons sont soit des demi droites infinies et tmax= inf, soit des segments, par exemple pour vérifier qu'un point est à l'ombre, ou pour vérifier qu'il existe une intersection plus proche que celle que l'on vient de trouver.

Hit intersect( const Plan& plan, const Point& o, const Vector& d, const float tmax )
{
    float t= ... ;
    Point p= ... ;
   
    if(t < 0 || t > tmax)

        return Hit();                           // renvoie une intersection non valide / par defaut. l'intersection n'est pas valide / derriere l'origine du rayon, ou trop loin

    else

        return Hit(t, p, plan.n, plan.color);   // renvoie la position de l'intersection, + la normale et la couleur du plan
}


remarque : si vous êtes à l'aise en c++, n'hésitez pas à utiliser des fonctions membres... par exemple :

struct Plan
{
    Point a;        // point sur le plan
    Vector n;       // normale du plan
    Color color;    // couleur

    Hit intersect( const Point& o, const Vector& d );
};

vous pouvez également créer une classe de base pour les objets et rendre la fonction intersect() virtuelle, ce qui simplifiera les calculs d'intersections avec des objets différents, par exemple des sphères et des plans.
struct Objet
{
    Objet() {}
    virtual ~Objet() {}

    virtual Hit intersect( const Point& o, const Vector& d, const float tmax );
};

struct Plan : public Objet
{
   
...
};

struct Sphere : public Objet
{
   
...
};




exercice 2 : plein d'objets !
une sphère posée sur un plan est pratique pour écrire la première version du code, mais on est rapidement lassé de voir la même image...
il n'est pas très compliqué de stocker plusieurs sphères dans un std::vector :

struct Scene
{
    std::vector<Sphere> spheres;
    Plan plan;
};

et d'écrire une fonction d'intersection qui renvoie l'intersection valide la plus proche trouvée en testant toutes les sphères + le plan :

Hit intersect( const Scene& scene, const Point& o, const Vector& d )
{
    Hit plus_proche;
    plus_proche.t= inf;
    for(unsigned i= 0 ; i < scene.spheres.size(); i++)

    {
        // tester la ieme sphere
        Hit h= intersect(scene.spheres[i], o, d, plus_proche.t);
        if(h ... )
            plus_proche= ... ;

    }

    // et bien sur, on n'oublie pas le plan...
    Hit h= intersect(scene.plan, o, d, plus_proche.t);
    if(h ...)
        plus_proche= ...;

    return plus_proche;

}

et voila !
ce travail préparatoire va permettre de passer à la suite très facilement !

vous pouvez également ajouter une source de lumière (ie la direction vers le soleil et son emission) dans la structure scene et même prévoir une structure Source, ainsi qu'un tableau de sources de lumières...



pour les curieux :
on peut se simplifier encore la vie en utilisant une fonctionnalité du c++, en surchargeant l'operateur conversion vers bool de la structure Hit :
struct Hit
{
    float t;        // position sur le rayon ou inf
    Point p;
    Vector n;       // normale du point d'intersection, s'il existe
    Color color;    // couleur du point d'intersection, s'il existe

    operator bool() { return t >= 0 && t < inf; } // renvoie vrai si l'intersection existe et quelle est valide...
};


à quoi ça sert ? et pourquoi ? on va pouvoir écrire les tests comme ca :

    Hit h= intersect( ... )
    if(h)
    {
        // faire un truc avec l'intersection h
    }

`if(h)` appelle la fonction que l'on vient d'ajouter dans la structure hit, c'est la même chose que si l'on écrit le test à chaque fois `if(hit.t >= 0 && hit.t < inf)` mais c'est plus simple et ça peut éviter certaines erreurs.

ou encore plus simple :

    if(Hit h= intersect( ... ))
    {
        // faire un truc avec l'intersection h
    }


remarque :
on peut aussi définir l'opérateur de comparaison entre 2 structures Hit pour vérifier laquelle est la plus petite...


Partie 5 : ombres

exercice 1 :
maintenant que l'on connait l'objet visible pour chaque pixel, il ne reste plus qu'à déterminer s'il est éclairé ou à l'ombre. comment faire ?
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 est éclairé.

il suffit de construire un nouveau rayon entre p et la source, et de calculer les intersections. s'il y a une intersection valide avec t > 0 le point est à l'ombre !

origine du rayon : p, le point d'intersection
direction du rayon : l, la direction de la source de lumière

exemple de défauts sur les ombres...



exercice 2 : et des ombres propres ? c'est possible ?
bien sur ! il suffit de décaler l'origine du rayon, il faut le séparer de la surface, les calculs avec des floats ne sont pas exacts et un point d'intersection ne se trouve jamais exactement sur la surface, il peut être au dessus ou en dessous de la surface de l'objet, et lorsque l'on teste les intersections pour savoir si p est éclairé, on trouve presque toujours une intersection, mais avec p lui meme...
pour régler ce problème, il suffit de décaler l'origine du rayon en s'éloignant de la surface dans la direction de la normale.

en gros, au lieu de construire le rayon comme ça :
Point o= p;
Vector d= l;

ce qui produit une image bizarre, pleine de défauts, comme au dessus, on va utiliser :
Point o= p + epsilon * n;    // avec epsilon= 0.001 par exemple...
Vector d= l;


mais ça se corrige facilement !

pour les curieux : plus d'explications sur ce problème et une meilleure solution dans la doc en ligne : "precision numérique et lancer de rayons"


Partie 6 : anti-aliasing

relisez le cours et générez plusieurs rayons par pixel et calculez une image "lisse".
ou, calculez une image plus grande et utilisez opencv pour la filtrer et la réduire à la bonne résolution. les explications sur opencv et les filtres sont dans le cours d'analyse / traitement d'images.

rappel : oui, les 2 manières de faire sont équivalentes. le filtre le plus simple est une moyenne. générer plusieurs rayons dans un pixel (sur une grille à l'intérieur de chaque pixel, ou aléatoirement) et calculer la moyenne des couleurs calcule exactement la même chose...


(cliquez pour une version plus grande)
   

(cliquez pour une version plus grande)

Intermède : plus vite !

tous les ordinateurs disposent de plusieurs coeurs de calculs, ce qui permet de réaliser les calculs plus vite. comment les utiliser ? il faut écrire un programme parallèle... cf le cours de système.
le plus simple est d'utiliser openMP, il suffit d'ajouter une ligne dans le programme principal autour de la boucle qui calcule les pixels de l'image et de compiler avec les bonnes options :
#pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    for(int px= 0; px < image.width(); px++)
    {
        ...
        image(px, py)= { ... };
    }

pour compiler, c'est déjà prévu (sous linux) : make config=release tp, au lieu de make tp (qui génère la version debug du programme).

pour visual studio, il faut cocher l'option openmp dans les propriétes du projet.
pour mac, c'est plus compliqué, il faut installer g++ avec brew, le gestionnaire de paquets, le compilateur apple ne supporte pas openmp...

temps d'exécution avant (cf make tp) : time bin/tp, 1m34s
temps d'exécution openMP (cf make config=release tp) : time bin/tp, 8s

pas mal pour une seule ligne changée dans le programme ! (12 coeurs sur la machine de test et un facteur d'accélération de 11.75, c'est plutot bien !)

qu'est ce qui a changé ? chaque ligne de l'image est calculée par un thread en parallèle sur les coeurs de calcul du processeur. regardez le moniteur systeme pendant le calcul d'une image, tous les coeurs devraient être actifs...

par contre, il va falloir faire attention aux variables globales. chaque ligne de l'image est maintenant calculée par un thread, et les threads n'aiment pas trop modifier les mêmes variables en meme temps, les générateurs de nombres aléatoires, par exemple. une solution simple qui fonctionne sans défaut consiste à utiliser un générateur différent pour chaque ligne de l'image :
#pragma omp parallel for schedule(dynamic, 1)
    for(int py= 0; py < image.height(); py++)
    {

        // initialiser un generateur de nombre aléatoire pour chaque ligne / chaque thread
        std::random_device hwseed;
        std::default_random_engine rng( hwseed() );
        std::uniform_real_distribution<float> uniform;
        // chaque thread a une copie privee de ces 3 variables

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

        {   
            // generer les nombres aleatoires avec uniform( rng )
            float ux= uniform( rng );
            float uy= uniform( rng );
       
            ...

            image(px, py)= { ... };
        }
    }

remarque : et l'image ? oui, c'est bien une variable globale, mais chaque thread modifie une ligne différente, du coup... pas de problèmes !