précision numérique et lancer de rayons

Le lancer de rayon fonctionne plutot bien :

mais dès que l'on calcule des ombres, on tombe sur un probleme assez pénible :

les images sont pleines de défauts... pourquoi ?

les float utilisés pour les calculs ne sont que des approximations, plus ou moins précises des nombres réels, et le point d'intersection exact peut se trouver n'importe ou dans une région autour de la valeur des float qui le représente :

et lorsque l'on teste la visibilité de la source de lumière, on utilise ce point d'intersection comme origine du rayon de test. Selon la position de ce point, le résultat sera correct (s'il est au-dessus de l'objet intersecté par le rayon), ou pas (lorsqu'il est sous la surface). Une solution classique décale l'origine du rayon de test pour l'éloigner de la surface de l'objet :

mais il faut choisir de quelle distance éloigner l'origine du rayon de test... une solution simple utilise une constante, 1/1000 par exemple :

Point p= { ... }; // intersection du rayon 
Vector n= { ... }; // normale au point d'intersection
 
// decale l'origine du rayon vers la source de lumiere
const float epsilon= float(1.0 / 1000.0);
Point origine= p + epsilon * n;

mais cette valeur, ne fontionnera pas avec toutes les scenes... pourquoi ?

parce que la précision des float est relative, pas absolue... plus un float est grand, moins il est précis, et la taille de la boite / la distance nécessaire pour décaler l'origine du rayon de test doit augmenter. La cornell box est inclut dans un cube de taille 10, alors que le manege est dans un cube de taille > 10000... et il faut un décalage très différent dans ces 2 cas... C'est meme pire, pour le manege ou le grand huit en arriere plan, il faut aussi des valeurs différentes...

pour vous convaincre, executez ce petit code qui affiche le premier float qui ne peut pas représenter exactement valeur + 1/1000 :

float epsilon= 1.0 / 1000.0;
for(int i= 0; i < 1000000; i++)
{
    if(float(i) + epsilon == float(i))
    {
        printf("%.15f + %.15f = %.15f\n", float(i), epsilon, float(i) + epsilon);
        break;
    }
}

et les valeurs plus grandes que 32000 ne permettent pas de représenter exactement 32000 + 1/1000...

les float 32bits représentent les réels en utilisant un entier 24bits et un exposant sur 8bits, ils ne peuvent donc représenter exactement que \( 2^{32} \) réels... chaque float représente donc un intervalle de réels, qui grandit avec la valeur du float. on peut connaitre la taille de cet intervalle en utilisant std::numeric_limits<float>::epsilon(), qui représente l'intervalle pour la valeur 1. pour les autres valeurs, il suffit de calculer \( valeur * (1 + epsilon) \), pour obtenir le prochain float représentable :

voila la taille des intervalles pour les valeurs < 100000 :

for(int i= 0; i < 100000; i+= 100)
{
    float f= float(i);
    float fpmax= f * (1 + std::numeric_limits<float>::epsilon());
    printf("%f %f\n", f, fpmax - f);
}

pour vous convaincre que ces valeurs ont vraiment un sens, vérifiez que si l'on ajoute une valeur inférieure à epsilon le résultat de l'addition est faux... par exemple :

{
    float x= 1;
    x= x + std::numeric_limits<float>::epsilon();
    assert(x > 1); // le resultat doit etre representable exactement
}
{
    float x= 1;
    x= x + std::numeric_limits<float>::epsilon() / 2;
    assert(x == 1); // le resultat n'est pas representable. on ajoute une valeur trop petite...
}

ou pour des valeurs > 1 :

{
    float x= 100000;
    x= x * (1 + std::numeric_limits<float>::epsilon());
    assert(x > 100000);
}
{
    float x= 100000;
    x= x * (1 + std::numeric_limits<float>::epsilon() / 2);
    assert(x == 100000);
}

et alors ?

comment utiliser la taille de l'intervalle / la precision d'un float pour décaler l'origine du rayon de test ?

première étape : estimer la precision du point d'intersection, le plus simple est de garder l'intervalle de la plus grande valeur des coordonnées x, y, z du point :

float epsilon_point( const Point& p )
{
    // plus grande erreur
    float pmax= std::max(std::abs(p.x), std::max(std::abs(p.y), std::abs(p.z)));
    
    // evalue l'epsilon relatif du point d'intersection
    float pe= pmax * std::numeric_limits<float>::epsilon();
    return pe;
}

et il ne reste plus qu'à décaler l'origine du rayon, le long de la normale de la surface :

Point p= { ... }; // intersection du rayon 
Vector pn= { ... }; // normale au point d'intersection
 
Point origine= p + epsilon_point(p) * pn;

et voila le resultat :

à gauche, sans déclage, à droite, avec... ce n'est pas suffisant pour régler le problème, pourquoi ?

tout simplement, parce que le point d'intersection est le résultat d'un calcul et que toutes les valeurs utilisées sont aussi arrondies par les float, et l'erreur réelle du point d'intersection est plus importante, elle dépend aussi des erreurs associées à toutes les valeurs utilisées et aux opérations utilisées. au minimum, le point d'intersection est calculé par \( p = o + t \vec{d} \) mais \( o, t\) et \( \vec{d} \) sont aussi représentés par des float...

il est possible de calculer assez précisement comment se propagent toutes ces erreurs, mais une solution très simple consiste à augmenter l'erreur estimée par une constante, voila quelques exemples :

p + 1 * epsilon_point(p) * pn	p + 2 * epsilon_point(p) * pn
p + 4 * epsilon_point(p) * pn	p + 32 * epsilon_point(p) * pn

et cette solution p + 32 * epsilon_point(p) * pn marche correctement sur des scènes très différentes :

et si l'extrémité du rayon est aussi sur un objet ?

il faut decaler l'origine (le point p) et l'extrémité (le point q) du rayon de test :

Point p= { ... };
Vector pn= { ... };
 
Point q= { ... }:
Vector qn= { ... };
 
Point origine= p + 10 * epsilon_point(p) * pn;
Point extremite= q + 10 * epsilon_point(q) * qn;
Ray test(origine, extremite);
if(Hit hit= intersect(ray ... )) 
    ... ;

c'est quand meme moche ? non ?

un peu oui, on a simplement remplacé une constante absolue 1/1000 par une constante relative \( 10*\epsilon(x) \) ... qui a tout de meme l'avantage d'etre simple à estimer et raisonnablement robuste.

PBRT explique en détail comment estimer les erreurs lors des calculs et comment modifier le lancer de rayon en fonction : Managing rounding errors