L3 synthèse d'images
2022

TP1 - lancer de rayons et images


Installation

vous pouvez faire le tp :

sur les machines du nautibus sous linux :
les librairies de la base de code complete gkit2light sont deja installées, il n'y rien a faire. il n'y a pas de différence entre les 2 versions de la base de code.

sur votre machine sous linux :
il suffit d'installer les librairies si vous souhaitez utiliser la base de code complete gkit2light. cf la doc / installation.

sur votre machine sous windows : le plus simple est d'utiliser la base de code minimaliste gkit3 et le compilateur officiel de microsoft
cf visual studio community 2022 (pas visual studio code, ce n'est pas du tout la meme chose...)
vous ne devez installer les librairies que si vous voulez vraiment utiliser la base de code complete, cf la doc / installation (mais encore une fois ce n'est pas necessaire pour ce cours)


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 3 : générer les projets et "étape 5 : créer un projet".
pour démarrer, le plus rapide : créez un fichier tp.cpp dans le répertoire `projects`, 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

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 verifier 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.


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, 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(1024, 1024);    // 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 -5) 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 -5 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 le 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 : "frontale" et cos theta
pour commencer simplement, on va considérer que la camera éclaire aussi 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 3 é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 ?


exercice 2 : et avec une vraie lumière ?
complétez le calcul en divisant par le carré de la longueur de l, le vecteur entre le point d'intersection et la source de lumière.
ajustez l'émission de la source pour obtenir une image pas complètement noire ! ou enregistrez l'image au format .hdr et utilisez image_viewer...

rappel : le carré de la distance entre 2 points se calcule avec : float distance2(a, b); et le carré de la longueur d'un vecteur : float length2(v), cf la doc



les images précédentes sont obtenues avec un plan a= (0 -1 0), n= (0 1 0) et une sphère c= (0 0 -3), r= 2
si la sphère est plus loin de la camera, elle apparaitra toute noire, pareil pour le plan... et si la sphère est plus proche, elle sera toute rouge / sur-exposée...
et oui un point qui éclaire des objets est comparable à une bougie... et il est assez difficile de placer la lumière pour obtenir une image "propre", ni toute noire ni entièrement sur-exposée...

indication : si vous galérez trop à placer la source de lumière ou les objets, vous pouvez remplacer la division par le carré de la distance entre le point et la source, par la distance entre le point et la source, ce sera plus simple à ajuster...


exercice 3 : qu'est ce qui change pour une source de lumière placée ailleurs ?
placer la source de lumière à une autre position.


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, pareil 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
    Vector n;       // normale du point d'intersection, s'il existe
    Color color;    // couleur du point d'intersection, s'il existe
};

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 )
{
    ...

    if(t < 0)
        return { inf, plan.n, plan.color };  // l'intersection n'est pas valide / derriere l'origine du rayon
    else
        return { t, 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...


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 )
{
    for(unsigned i= 0 ; i < scene.spheres.size(); i++)
    {
        // tester la ieme sphere
        Hit h= intersect(scene.spheres[i], o, d);
        if(h ... )
    }

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

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

vous pouvez également ajouter la source de lumière (sa position 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 :
struct Hit
{
    float t;        // position sur le rayon ou inf
    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
    }


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 position de la source, et de calculer les intersections. s'il y a une intersection le point est à l'ombre !

origine du rayon : p, le point d'intersection
direction du rayon : l, le vecteur entre p et la position de la source.

indications : 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, t > 0 et t < 1


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 etre 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= Vector(o, source);

ce qui produit une image bizarre, pleine de défauts :


exemple de défauts sur les ombres...

on va utiliser :
Point o= p + epsilon * n;    // avec epsilon= 0.001 par exemple...
Vector d= Vector(o, source);


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 : reflets

exercice 1 :
relisez le cours et construisez le vecteur réfléchi dans la direction miroir connaissant la direction du rayon incident, le point d'intersection, et sa normale.





exercice 2 :
calculez les intersections avec les objets de la scène dans la direction miroir. quelle est l'origine du rayon ? sa direction ?
attention : pensez à décoller l'origine du rayon de la surface, sinon il y aura les mêmes problèmes que pour les ombres...

on se contente d'utiliser la couleur de l'objet touché par ce rayon pour vérifier que la direction et les intersections sont correctes.

exercice 3 :
quelle couleur faut-il utiliser pour le reflet ? le miroir réfléchit son environnement. si un objet réfléchi dans le miroir est à l'ombre, le miroir réfléchit l'ombre, sinon il réfléchit la couleur de l'objet éclairé... donc pour obtenir la "bonne" couleur, il faut calculer la lumière réfléchie par l'objet visible dans la direction miroir...

il est sans doute nécessaire de restructurer le code du tp et de créer une fonction Color eclairage_direct( const Point& p, const Vector& n, const Color& color, const Scene& scene ) qui regroupe les calculs de lumière et d'ombre. il ne reste plus qu'à l'utiliser.
 
exercice 4 : et les coefficients de Fresnel ?
relisez le cours et multipliez la couleur de l'objet visible dans la direction miroir par le coefficient de reflexion de Fresnel. vous pouvez utiliser la formule complète ou l'approximation pour calculer le coefficient à partir de l'indice de réfraction de la matière de l'objet, ou la couleur F0.

rappel : l'eau à un indice de réfraction de 1.3, le verre de 1.5, et un F0 de 0.02 (eau) ou 0.04 (verre).

pour les curieux : et les métaux ? les métaux sont des conducteurs, et la lumière (une onde electro-magnétique) se comporte différemment avec les métaux. les formules complètes sont encore plus pénibles à écrire et à valider. on peut utiliser directement l'appoximation, manipuler F0 est plus simple que de faire les calculs avec les indices de réfraction complexes (avec une partie imaginaire) retrouvez le tableau avec les bonnes couleurs dans le cours (cf fer, cuivre, etc, colonne RGB lineaire).


exercice 5 : description de matière.
modifiez la description des objets pour ajouter l'indice de réfraction (ou F0, selon le cas) en plus de la couleur diffuse. certains objets seront des miroirs, les autres seront diffus. il faut également adapter le calcul de la couleur de l'objet...

pour les curieux : et si l'objet visible dans la direction miroir est aussi un miroir ? pour connaitre la lumière qu'il réfléchit, il faut aussi calculer la couleur de l'objet qu'il réfléchit... le calcul devient récursif.


Partie 7 : 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.

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, les générateurs de nombres aléatoires par exemple. une solution simple qui fonctionne sans défaut consiste à utiliser un générateur par 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> u;
        // chaque thread a une copie privee de ces 3 variables

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

        {   
            // generer les nombres aleatoires avec u(rng)
            float ux= u(rng);
            float uy= u(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 !


Partie 8 : pénombres ! avec plein de sources...

il y a plusieurs manières d'avoir un éclairage un peu plus naturel dans les images, mais il va falloir utiliser plus de sources et faire plus de calculs...

exercice 1 : points dans une grille
pour démarrer, le plus simple est de modifier le code pour utiliser un tableau de sources. relisez le cours pour disposer des points dans une grille / panneau.


(cliquez pour une version plus grande)

écrivez une fonction qui génère n points dans une grille passée en paramètre :
std::vector<Source>& genere_sources( const Point& a, const Vector& u, const Vector& v, const Color& emission, const int n );
il est également plus simple de fixer l'emission totale du panneau lumineux et de la répartir sur les n points. comme ça vous aurez à peu près la même image quelque soit le nombre de points...

ensuite, dans la boucle de calcul des pixels de l'image, il suffit de parcourir ce tableau de sources de lumière pour calculer la lumière réfléchie par un point d'intersection, connaissant sa position, sa normale et les parametres de sa matière... on peut aussi écrire une fonction :
Color eclairage_points( /* point d'intersection */ const Point& p, const Vector& n, /* matiere */ const Color& color, /* sources */ const std::vector<Source>& points )

exercice 2 : et avec des points aléatoires, qu'est ce qui change ?
au lieu de générer les points au début du programme comme dans l'exercice précédent, il faut les générer à chaque fois que l'on calcule la lumière réfléchie par un objet. la fonction qui calcule cette lumière réfléchie doit donc connaitre les parametres du panneau lumineux et le nombre de points à générer. il faudra aussi passer en paramètre un générateur de nombre aléatoire.
Color eclairage_panneau( /* parametres du point d'intersection */ const Point& p, const Vector& n, /* matiere */ const Color& color,
    /* parametres du panneau */ const Point& a, const Vector& u, const Vector& v, const Color& emission, const int n, std::default_random_engine& rng )
indication : c'est sans doute une bonne idee de créer une struture pour stocker tous les parametres qui decrivent un panneau lumineux.


(cliquez pour une version plus grande)


Partie xx : éclairage ambient

relisez le cours.
comment constuire un ensemble de sources de lumière ou de directions disposées sur un dome ?

   



Partie xx : flou de profondeur

Partie xx : bouger la camera ?

Partie xx : réfraction et objets transparents