L3 synthèse d'images
2023

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.

un squelette de programme pour le tp est disponible dans projets/tp2.cpp

générer les projets

pourquoi ? gkit compile et fonctionne sur linux, windows, mac os, ios, android et meme WebGL. Chaque système dispose de plusieurs compilateurs et environnements de travail. Il n'est pas envisageable de créer et de maintenir tous ces projets manuellement. gkit utilise donc un outil : un générateur de projet, ce qui permet de décrire les projets une seule fois et c'est l'outil (premake dans ce cas...) qui génère le projet pour votre environnement de travail.

il faut donc apprendre à générer le projet pour votre environnement de travail, en utilisant premake.

sous linux, premake4 est disponible, installez le sur votre machine, si nécessaire : sudo apt install premake4
sinon, téléchargez le dans le répertoire gkit3 : cf premake5

ouvrez un terminal, et naviguez jusqu'au répertoire contenant gKit :

rappel : commandes ls et cd pour naviguer.

windows + codeblocks

./premake5.exe codeblocks

le workspace (groupe de projets) codeblocks ainsi que les projets sont crées dans le répertoire build/, ouvrez build/gKit3.workspace.

windows + visual studio

pour générer une solution (groupe de projets) visual studio, il suffit de choisir la bonne version :

./premake5.exe vs2022

la solution visual studio ainsi que les projets sont crées dans le répertoire build/, ouvrez build/gkit3.sln.

mac os + xcode

./premake5 xcode

mac os + makefile

./premake5 gmake

le Makefile se trouve dans le répertoire de base de gKit.

linux + makefile

premake4 gmake // si premake4 est installe dans le système
./premake5 gmake // si premake5 est copie dans le repertoire de gKit

le Makefile se trouve dans le répertoire de base de gKit.

remarque : si premake5 est disponible dans les paquets de votre distribution utilisez-le ! cf premake5 gmake

linux + vscode

générez les makefiles, comme au dessus dans linux + makefile

compilez un exemple

compilez tp1, par exemple, si vous voulez vérifiez qu'une application compile et fonctionne.

utilisation des makefiles

les makefile peuvent générer les versions debug (cf utiliser un debugger comme gdb ou lldb) ou les versions release, plus rapide (2 ou 3 fois, interressant pour les projets avec beaucoup de calculs) :

les exécutables sont crées dans le répertoire gkit3/bin, pour les exécuter : 

bin/tp1

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 etre à 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 : sphere
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 sphere, 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 sphere :

le plus simple est de placer une sphere 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 sphere. il y a donc au moins 4 tests à écrire.


remarques :
vous pouvez également créer une structure sphere 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. mais ca ne sera utile que plus tard, pour calculer la lumière réfléchie vers le pixel / la camera.


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 sphere de centre (0 0 -3) et de rayon 2, calculez l'image, coloriez un pixel dont le rayon touche la sphere 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...
relisez le cours sur les astuces pour simplifier le code de comparaison de plusieurs intersections...


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 ?


exercice 3 :
vérifiez que l'intersection calculée est bien la plus proche de la camera. le plus simple et déplacer la sphère sous le plan :


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 ?

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.

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 parametre la structure qui stocke les parametres 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
}

exercice 2 : plein d'objets !
une sphere posée sur un plan est pratique pour écrire la premiere version du code, mais on est rapidement lassé de voir la meme 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 meme chose que si l'on ecrit 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'interieur de chaque pixel, ou aléatoirement) et calculer la moyenne des couleurs calcule exactement la meme 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 !