Partie 1 : prise en main
écrivez l'intersection rayon triangle ou rayon sphère.
Dans un premier temps, testez votre fonction d'intersection dans
un fragment shader, avec shader_kit, par exemple.
Partie 2 : compute shader
que faut-il modifier dans votre programme pour faire ce calcul
dans un compute shader ?
ou sont écrits les résultats ?
comment les afficher ?
comment passer un ensemble de triangles ou de sphères au shader ?
rappels : storage
buffers,
storage
textures
rappels : le constructeur de la classe
App
permet de préciser la version d'openGL, en plus des dimensions de
la fenetre, cf App(1024,640, 4,3)
Partie 3 : intégration numérique et shaders
L'objectif est de calculer l'éclairage ambiant en utilisant la
même solution que dans le tp précédent : générer des directions
aléatoires, tester la visibilité, etc.
Vos tests d'intersections doivent fonctionner, il ne reste plus
qu'à voir comment générer des nombres aléatoires dans un shader...
La solution la plus simple est d'utiliser un
générateur
linéaire : x
n+1= (ax
n + b) % m, et de
choisir de
"bonnes"
valeurs pour a, b, et m, (et x
0, bien sur)
utilisez celles de la glibc par exemple.
Reste un détail : comment paralléliser le générateur ? il faut
obtenir une séquence de nombres aléatoires différente par pixel ?
Une solution directe utilise une valeur x
0 différente
pour chaque pixel, mais rien ne garanti que les séquences générées
pour des pixels voisins seront indépendantes...
Il suffit de se rendre compte que lorsque l'on calcule une image
de manière séquentielle, chaque pixel utilise une séquence de
termes de la série aléatoire. Le pixel (0, 0) utilise les termes 0
à N, le pixel suivant, les termes N à 2N, etc. Il ne reste plus
qu'à trouver comment calculer la valeur de x
pN pour
obtenir les termes x
pN à x
(p+1)N de chaque
pixel p de l'image. Les générateurs de la stl ne fournissent pas
une version utilisable de cette fonctionnalité, mais bien sur, il
y a une solution.
On peut reformuler x
n= (ax
n-1 + b) % m = (a
nx
0
+ b (a
n-1) / (a-1)) % m. On peut calculer le terme x
n
directement en fonction de x
0, sans évaluer les N
termes successivement. Mais le calcul n'est pas direct, on ne peut
pas représenter la valeur a
n. Mais comme tous les
calculs sont fait modulo m, il y a quand même une solution, cf
"Random
Number Generation with Arbitrary Strides", 1994, F. Brown.
Voici une classe RNG simplifiée qui inclut cette fonctionnalité :
struct RNG
{
unsigned int x;
unsigned int x0;
RNG( const unsigned int seed ) : x(seed), x0(seed) {}
// glibc
static const unsigned int a= 1103515245;
static const unsigned int b= 12345;
static const unsigned int m= 1u << 31;
float sample( ) // renvoie un reel aleatoire dans [0 1]
{
x= (a*x + b) % m;
return float(x) / float(m);
}
unsigned int index( const size_t i ) // prepare la generation du terme i
{
unsigned int cur_mul= a;
unsigned int cur_add= b;
unsigned int acc_mul= 1u;
unsigned int acc_add= 0u;
size_t delta= i;
while(delta)
{
if(delta & 1u)
{
acc_mul= acc_mul * cur_mul;
acc_add= acc_add * cur_mul + cur_add;
}
cur_add= cur_mul * cur_add + cur_add;
cur_mul= cur_mul * cur_mul;
delta= delta >> 1u;
}
x= acc_mul * x0 + acc_add;
return x;
}
};
int main()
{
RNG rng(1);
for(int i= 0; i < 10; i++)
// implicite :
rng.index(i);
printf("%f\n",
rng.sample());
// affiche les 10 premiers termes dans
l'ordre
for(int i= 9; i >= 0; i--)
{
rng.index(i);
// affiche les memes
valeurs dans un ordre different
printf("%f\n",
rng.sample());
}
// prepare un generateur
unsigned int seed=
rng.index(9);
// recupere l'etat permettant de generer le
terme 9
RNG tmp(seed);
printf("%f\n",
tmp.sample());
// genere le terme 9
avec un autre generateur
return 0;
}
pour les curieux : on peut aussi utiliser une autre famille
de générateurs,
pcg,
une évolution des générateurs linéaires, qui ajoute une
transformation après le calcul de x
n+1. Les variantes
utilisables dans un shader (état interne sur 32bits) sont décrites
sur
pcg-random.org
Utilisez RNG::index() pour générer le premier terme de la séquence
aléatoire de chaque pixel. Pour les autres termes du pixel, il
suffit d'utiliser le générateur "normalement", cf RNG::sample().
Vous pouvez vérifier que tout fonctionne en modifiant votre tp
précédent (remplacez les générateurs de la stl par la classe RNG).
Utilisez RNG::index() pour initialiser le générateur linéaire de
chaque pixel. le shader n'évaluera que l'équivalent de
RNG::sample() pour générer les nombres aléatoires nécéssaires pour
calculer l'éclairage ambiant du pixel.
Comment transmettre au shader l'état du générateur de chaque pixel
? Ecrivez également la fonction sample() dans votre shader.
Partie 4 : c'est un peu lent ? intégration temporelle
L'intégration numérique a besoin de trop d'échantillons pour
recalculer brutalement chaque image et garder un temps d'affichage
interactif. Il est relativement simple de découper le problème :
chaque exécution du shader calcule une partie du résultat et
accumule petit à petit les N échantillons. Le temps total de
calcul ne change pas, mais l'application peut rester interactive.
Lorsque la camera (ou les objets, ou les sources...) se déplace,
il suffit de repartir de zero, sans attendre la fin du calcul des
N échantillons.
indications : creation et manipulation de textures
format "couleur" classique, 8bits par canal,
normalisé entre 0 et 1
glGenTextures(1, &m_image);
glBindTexture(GL_TEXTURE_2D, m_image);
glTexImage2D(GL_TEXTURE_2D, 0,
GL_RGBA8, window_width(), window_height(), 0,
GL_RGBA, GL_UNSIGNED_BYTE, /* data */);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAX_LEVEL, 0);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER,
GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER,
GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S,
GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T,
GL_CLAMP_TO_EDGE);
// pour
"effacer", la texture
glClearTexImage(m_image, 0, GL_RGBA, GL_UNSIGNED_BYTE,
nullptr);
// pour
selectionner la texture
glBindImageTexture(/* unit */ 0, /* texture */ m_image, /*
level */ 0, /* layered */ GL_FALSE, /* layer */ 0,
/* access */ GL_READ_WRITE, /* format */ GL_RGBA8);
// shader
layout(binding=
0, rgba8) uniform image2D image;
vec3 color=
imageLoad(image, ivec2(gl_GlobalInvocationID.xy)).rgb;
format unsigned int
32bits, 1 canal :
glGenTextures(1,
&m_seed_image);
glBindTexture(GL_TEXTURE_2D, m_seed_image);
glTexImage2D(GL_TEXTURE_2D, 0,
GL_R32UI, window_width(), window_height(), 0,
GL_RED_INTEGER, GL_UNSIGNED_INT, /* data */);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAX_LEVEL, 0);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER,
GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER,
GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S,
GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T,
GL_CLAMP_TO_EDGE);
// pour
selectionner la texture
glBindImageTexture(/* unit */ 1, /* texture */
m_seed_image, /* level */ 0, /* layered */ GL_FALSE, /*
layer */ 0,
/* access */ GL_READ_WRITE, /* format */ GL_R32UI);
layout(binding=
1, r32ui) uniform uimage2D seeds;
uint seed=
imageLoad(seeds, ivec2(gl_GlobalInvocationID.xy)).x;
indication : les shaders (compute et fragment) peuvent
lire une texture avec les fonctions imageLoad()
layout(binding= 0, rgba32f) uniform image2D
image;
vec3 color= imageLoad(image,
ivec2(gl_GlobalInvocationID.xy), 0).rgb; //
compute shader
vec3 color= imageLoad(image, ivec2(gl_FragCoord.xy),
0).rgb;
// fragment shader
pour les curieux : c'est quand meme dommage de jeter tous
les calculs lorsque la camera bouge... une idée pour continuer à
utiliser les calculs déjà fait dans les zones qui ne changent pas
ou peu, consiste à identifier les points de la scène qui sont
toujours visibles dans la nouvelle image et à re-utiliser leur
valeur. Pour les autres points, il faudra repartir de zero (ou
d'une valeur moyenne des pixels voisins qui restent valides).
Ce principe d'accumulation temporelle est utilisé quasiment
partout, cf
"TAA
/ Temporal Supersampling", M. Salvi, GDC 2016 et
"Temporal
Reprojection Anti-Aliasing in INSIDE", Playdead, GDC 2016
Partie 5 : beaucoup de triangles
Lorsque le nombre de primitives / objets augmente, il devient
rapidement nécessaire d'utiliser une structure accélératrice,
comme un BVH, pour limiter le nombre de calculs d'intersections.
Les sections suivantes vous proposent une solution pour construire
un BVH, le modifier pour pouvoir le parcourir sans utiliser de
récursion (les shaders ne permettent pas d'écrire de fonctions
récursives), écrire la fonction d'intersection iterative et enfin
utiliser l'ensemble dans un compute shader pour calculer une
image, avec un éclairage ambiant, par exemple.
Commencez par construire un BVH avec une solution simple, cf
build_node_centroids() vue en cours. Le BVH sera représenté par un
tableau de noeuds, un tableau de triangles et éventuellement,
l'indice de la racine.
rappel : il n'y a pas de pointeurs
non plus, il faut donc référencer les noeuds en utilisant leur
indice.
Vérifiez que votre arbre et votre fonction de parcours récursif
sont corrects en modifiant votre tp précédent.
Prochaine étape, préparation du parcours itératif, la fonction
d'intersection qui devra fonctionner dans un compute shader.
Il y a au moins 2 stratégies pour éliminer la récursion du
parcours : manipuler explicitement une pile de noeuds à visiter
(un tableau d'indices de noeud), ou
coudre l'arbre. L'idée
de la couture est simple : il suffit de stocker dans chaque noeud
l'indice de son successeur dans l'ordre prefixe (noeud puis fils
gauche, suivi du fils droit) ainsi que l'indice de son fils
gauche. Le parcours devient très linéaire : si le rayon touche
l'englobant du noeud, le fils du noeud sera visité, sinon le
successeur du noeud sera visité.
Ecrivez une fonction qui transforme un BVH binaire en BVH cousu.
Ecrivez la fonction d'intersection de l'arbre cousu.
Vérifiez que tout est correct en modifiant votre tp précédent.
Il ne reste plus qu'à écrire la version compute shader de la
fonction d'intersection. Utilisez 2 shader storage buffers pour
transmettre le tableau de noeuds et le tableau de triangles au
compute shader.
Modifiez votre shader de la partie 4 pour vérifier que tout
fonctionne correctement.
indications : vous aurez également besoin d'une fonction
d'intersection rayon / boite englobante alignée sur les axes, par
exemple :
struct BBox
{
Point pmin, pmax;
bool intersect( const Ray& ray, const
Vector& invd, const Float htmax ) const
{
Point rmin= pmin;
Point rmax= pmax;
if(ray.d.x < 0) std::swap(rmin.x, rmax.x);
if(ray.d.y < 0) std::swap(rmin.y, rmax.y);
if(ray.d.z < 0) std::swap(rmin.z, rmax.z);
Vector dmin= (rmin - ray.o) * invd;
Vector dmax= (rmax - ray.o) * invd;
float rtmin= std::max(dmin.z, std::max(dmin.y,
std::max(dmin.x, float(0))));
float rtmax= std::min(dmax.z, std::min(dmax.y,
std::min(dmax.x, htmax)));
return (rtmin <= rtmax);
}
// avec Vector invd(1/ray.d.x, 1/ray.d.y,
1/ray.d.z);
};
indications : il est plutot recommandé de
se "débrouiller" pour que les 2 fils d'un noeud se trouvent sur la
même ligne de cache...
Il est aussi préferable de tester directement
l'englobant des 2 fils au lieu de suivre à la lettre
l'algorithme de parcours. pourquoi ?
(combien de fils sont visités en moyenne, pour un BVH "bien
construit" ?
combien de tests rayon / englobant sont nécessaires pour visiter
le fils gauche ? pour le fils droit ? comment re-equilibrer ?)
pour les curieux : oui, le nombre d'acces mémoire est un
gros facteur de ralentissement de la fonction d'intersection.
comment le réduire ?
Une première idée consiste à compresser la
représentation d'un noeud et d'arriver à stocker 4 noeuds, ou
plus, sur une ligne de cache de 64 octets (au lieu de 2
normalement). c'est ce que propose :
"Efficient
Incoherent Ray Traversal on GPUs Through Compressed Wide BVHs",
H. Ylitie, T. Karras, S. Laine, 2017
Une autre solution consiste à limiter le nombre de fois ou un
noeud est chargé par l'ensemble des processeurs... rappel : que
se passe-t-il lorsque les données chargées par N threads
saturent le cache L1 du processeur ?
"Dual
Streaming for Hardware-Accelerated Ray Tracing", K.
Shkurko... 2017
Mais la première limitation est l'exécution cohérente des threads
: que se passe-t-il lorsque les threads ne veulent pas faire la
même chose en même temps ?
Le processeur n'a pas le choix, il est obligé
d'exécuter la totalité des instructions pour que chaque thread
obtienne le bon résultat. Par exemple, un thread veut tester
l'intersection d'un rayon et d'un englobant, et un autre thread
du meme groupe d'exécution veut tester l'intersection
d'un rayon et d'un triangle ? Comment paralléliser N parcours
rayons / BVH ? un thread par rayon et un parcours en profondeur
? ou un thread par noeud et un parcours en largeur ?
"Coherent
Ray Tracing via Stream Filtering", C. Gribble, K. Ramani,
2008
"Understanding
the Efficiency of Ray Traversal on GPUs", T. Aila, S.
Laine, 2009
"Fast Ray
Sorting and Breadth-First Packet Traversal for GPU Ray
Tracing", K. Garanzha, C. Loop, 2010
