Flou, Fresnel et Blobs

Je rappelle que le code qui suit n'est pas particulièrement optimisé, il sert principalement d'illustration aux différentes notions qui sont abordées dans cette suite d'articles. Je suis désolé aussi pour les commentaires succincts et l'absence de "projets" (pour ceux qui ont l'habitude de travailler sur des projets avec Visual Studio) en contrepartie je m'arrange pour que tout soit compilable avec une simple ligne de commande et ne nécessite pas trop de librairies externes.

Aussi je peux paraitre sauter du coq à l'âne mais la synthèse d'image est un très vaste sujet qu'on ne peut pas couvrir avec dix pages d'articles et de code. Il s'agit donc de "browser" certains sujets suivant mon humeur et mon temps libre disponible. Les livres complets sur le sujet ne manquent pas et j'encourage les gens qui sont sérieux sur la volonté d'apprendre de compléter cette lecture par l'achat de quelques bouquins sur l'image de synthèse (et d'apprendre l'anglais accessoirement).

C'est la quatrième partie de notre série d'articles sur le raytracing en C++. Elle fait suite à la partie intitulée "Textures".

Modèle de caméra et profondeur de champ

Ce qui suit est une simplification d'un système optique tel qu'il existe dans l'oeil ou dans un appareil photographique.

Une caméra simplifiée est constituée d'un centre optique qui est le point où les rayons lumineux convergent, et d'un plan de projection. Ce plan de projection peut être placé n'importe où il peut même être situé dans la scène puisqu'il s'agit d'un calcul sur ordinateur et non pas d'un vrai appareil photo avec pellicule. Il dispose également d'un système optique qui est là pour faire converger les rayons qui passent à travers le diaphragme. Il faut en effet comprendre qu'une des grosses caractéristiques des caméras et qu'ils ont besoin de faire passer un faisceau de lumière pour que la pellicule ou la rétine soit suffisamment éclairées pour être marquées. Le système optique doit donc être configuré pour que le point d'intérêt de la scène soit net sur le plan de projection. C'est souvent mutuellement exclusif c'est à dire qu'on ne peut pas avoir deux points nets qui soient proche et lointain à la fois. En photographie, on appelle cela la profondeur de champ (depth of field). Notre modèle simplifié de caméra permet d'obtenir le même effet pour plus de réalisme (mais notre avantage par rapport aux photographes c'est que l'on peut s'en passer).

Le flou entrainé par la profondeur de champ est un peu lourd à mettre en oeuvre en terme de calcul. On va utiliser ce qu'on appelle une méthode stochastique. C'est à dire qu'on va faire intervenir du pseudo-hasard pour prendre en compte la complexité du problème. Cela consiste à imaginer que l'on suive un photon dans notre système optique et que ce photon ait une déviation aléatoire du centre optique causé par l'ouverture du diaphragme. Le système optique est une boite noire. Moins on en sait dessus mieux on se porte. Tout ce que l'on veut savoir c'est à quel distance de notre centre les points de la scène seront perçus comme nets et de quelle direction notre photon venait avant de heurter notre plan de projection afin de suivre son chemin inverse. On va utiliser les propriétés géométriques du système optique. Premièrement tout rayon passant directement par le centre optique n'est pas dévié. Deuxièmement tout rayon passant par un point du plan "net" de la scène, va converger sur un unique point du plan de projection, c'est la définition du plan "net".

En partant d'un point de projection on remarque donc que suivant le chemin suivi dans le système de lentille on va voir plusieurs points différents de la scène. On n'a pas besoin de beaucoup plus pour avoir un flou convaincant et facilement configurable (d'abord sa force et ensuite la distance où les objets apparaissent nets). Meme si après cela on pourrait avoir envie de simuler des appareillages réels avec un impact plus important sur la forme du flou créé.

On a rajouté au fichier scène des informations de perspective. Celle qui nous intéresse c'est la dispersion dans le code ci dessus. Si la dispersion est non nulle, on prend le rayon sortant du centre optique et on cherche la position du point net qui est le point visé par le rayon sur le plan net (situé à la distance "clearPoint"). Tous les rayons passant par ce point visé arriveront par construction sur notre point de projection. Par contre leur direction diffèrent, par exemple ici on perturbe légèrement le point de départ du rayon, en le faisant pointer sur le point net on obtient la nouvelle direction. On fait ça un certain nombre de fois par pixel, nombre défini par la variable "complexity". Si l'on veut un effet de flou convaincant et sans trop de bruit dans le résultat il faut un grand nombre de rayons par pixels. Il est sans doute possible de réduire ce nombre mais on ne s'intéresse pas trop à l'optimisation pour l'instant. Le code est déjà bien assez compliqué comme ça.

La "distance" virtuelle du plan de projection est déterminée par le FOV et par la taille du viewport.

Voici le résultat quand on place le plan net à distance de la sphère du milieu:

Et voilà le résultat quand le plan net est situé à grande distance (par rapport à la taille de l'ouverture), seul le paysage est net dans ce cas :

Pour obtenir ce niveau de qualité il a fallu envoyer un très grand nombre de rayons. Il a fallu plus d'une demi heure pour calculer chaque image. C'est beaucoup et il est évidemment possible d'optimiser le nombre de rayons à envoyer. Mais cela n'est pas notre priorité pour l'instant.

Isosurfaces (blobs)

Une isosurface est une surface qui est défini dans un champ de potentiel variable (et continu!) par une valeur de potentiel donné. On suppose ce potentiel scalaire et qui a les propriétés suivantes :

- les champs sont additifs : le champ créé par trois sources est égal au champ créé par une source multiplié par trois. En pratique dans la vie réelle, les champs magnétiques, électriques et gravitationnels sont aussi additifs c'est donc un domaine assez connu. Si on voulait représenter une surface isogravitationnelle approchée on pourrait utiliser ma méthode.

- la valeur du potentiel pour un point est défini comme l'inverse du carré de la distance à ce point. Grace à la propriété ci dessus, on a donc facilement la valeur du potentiel pour n points. Il suffit d'ajouter leur valeur.

- On peut approximer comme un cochon, vous pourriez etre en mal de différencier à l'oeil nu un blob de son approximation.

Voici une représentation d'isopotentielles en 2D. Il s'agit de trois attracteurs dans un champ gravitationel. En rouge on a tracé plusieurs isopotentielles pour des valeurs du champ différent et en bleue on a tracé la direction du gradient. Ce sont en fait les lignes de support de la force de gravitation.

Les caractéristiques de ce problème sont posées par définition. Il n'est pas la peine de chercher la justification de ces propriétés, elles sont données.

Passons maintenant à notre code : On va chercher à éviter la résolution d'équations polynomiales du nième degré qui dépasse le cadre de notre série d'article. À la place on va adopter une méthode approchée mais très suffisante pour nos scènes limitées. Chaque centre découpe l'espace en zones avec des sphères concentriques. En dehors de la sphère la plus extérieure c'est comme si le centre n'existait pas. C'est une optimisation rapide, en dehors de la zone maximale on fait comme si on ne subit l'influence que des N-1 sphères qui restent.

Voilà comment un rayon traverse notre champ de potentiel : la fonction de potentiel est nulle au départ du rayon les coefficients sont zéro, zéro, zéro. Puis il heurte la première sphère, la fonction de potentiel devient une fonction du premier degré en r^2, r étant la distance qui le sépare du centre de la première sphère et qui est une fonction linéaire de t. En remplaçant donc t dans la fonction de potentiel on obtient la variation de notre fonction de potentiel en fonction de t. On peut savoir si le rayon heurte notre isosurface en résolvant l'équation fonction de potentiel = notre seuil. C'est une équation du second degré en t donc facile à résoudre. On fait ça pour chaque nouvelle zone traversée. Une nouvelle zone est définie sur notre rayon comme l'intersection avec une nouvelle sphère qui change l'équation. On a de bonnes caractéristiques qui font que notre surface isopotentielle approchée sera continue.

Notre travail consiste donc à étiqueter chaque entrée/sortie d'une sphère d'influence de trier les points d'intersection en fonction de t, et de résoudre une à une les équations posées sur le trajet du rayon, jusqu'à trouver notre rayon d'intersection le plus proche.

Voici ci dessus le code qui fait le calcul d'intersection du rayon avec toutes les sphères d'influence. On peut les faire toutes une à une mais en pratique on fait pour chacune des sources (qui a dix sphères d'influences dans le code) le test d'intersection de la plus éloignée à la plus rapprochée, ainsi si une sphère éloignée n'est pas intersectée ce n'est pas la peine de chercher plus proche. S'il y a intersection, alors on a probablement deux valeurs de t qui correspondent aux deux intersections. On va placer ces valeurs de t dans une liste triée (ici une map) qui comprend en association avec des valeurs, la variation de l'expression du champ de potentiel à ces frontières. À cause de la symétrie du problème, la variation en sortie est toujours l'opposée de celle en entrée. Les variations sont exprimées en terme de coéfficients polynomiaux a, b et c.

Une fois que l'on a notre liste triée avec les variations de l'expression du champ, on va parcourir cette liste depuis l'extérieur jusqu'à ce qu'on trouve un point d'intersection intéressant (dans notre intervalle d'intersection possible). J'espère que le code parle de lui-même.

Une fois que l'on a trouvé le point d'intersection il reste encore la normale à la surface nécessaire pour les calculs d'éclairage et de réflexion/refraction. Vous verrez parfois suggéré d'approcher la normale comme une moyenne pondérée de normales des sphères d'influence. Mais on peut faire beaucoup mieux, comme il s'agit d'une surface d'isopotentielle, une des belles propriétés de ces surfaces est que la normale en un point est dans la direction du gradient du champ de potentiel.

Je vous épargne le calcul du gradient d'un potentiel en 1/R^2, voici le code résultant :

Voici le résultat avec trois "sources ponctuelles" :

Fresnel et les volumes transparents

Certains matériaux ont la propriété de transmettre la lumière (le verre, l'eau) mais affectent sa progression ce qui conduit a deux phénomènes. Le premier est la réfraction, le rayon lumineux est propagé mais est dévié de sa trajectoire initiale selon une formule formulée par Descartes (ou Snell pour les anglophones).

En gros si ThetaI est l'angle d'incidence par rapport à la normale à la surface alors on a ThetaT, l'angle de refraction (angle du rayon transmis) qui est tel que n2 * sin(Thetat) = n1 * sin(ThetaI). C'est la fameuse loi de Snell-Descartes.

n est appelée indice de refraction du milieu. C'est une variable sans unité calculée en fonction de celle du vide qui par définition est égale à 1. n1 est du coté de l'incidence, n2 est du coté de la transmission. Bien évidemment la loi est entièrement symétrique (en optique classique il n'y a pas de différence entre un rayon arrivant et un rayon partant). Cela veut dire que ThetaT et ThetaI sont interchangeables.

Symétriques à une exception près. si n1 > n2, alors il y a des angles ThetaI pour lesquels il n'y a pas d'angle ThetaT solution. Au delà d'un certain angle d'incidence, il n'y a donc tout simplement pas de rayonnement transmis.

Le deuxième effet est la reflexion d'une partie de la lumière. Cette reflexion est toujours dans la direction refletée par la normale par contre son intensité varie en fonction de la quantité de lumière transmise par réfraction.

Quelle est la quantité de lumière transmise et la quantité de lumière réfléchie ? C'est là que Fresnel intervient. L'effet de Fresnel est quelque chose d'assez complexe qui fait intervenir la nature ondulatoire de la lumière. Jusqu'à présent on a traité les surfaces réfléchissantes comme de simples miroirs. Le reflet d'un miroir est du à l'effet metalique des électrons excités qui renvoient intégralement le rayon lumineux dans la direction réfléchie.

Cependant la plupart des surfaces ne sont pas totalement réflechissantes comme un miroir. La plupart du temps elles réfletent un peu la lumière à certains angles et deviennent très réfléchissantes à des angles bas. C'est le cas de matériaux comme l'eau ou le verre. Mais également de nombreux autres matériaux qui ne sont pas considérés comme "transparents". Une feuille de papier réflète une partie de la lumière incidente, celle qui arrive aux angles très bas et elle diffuse le reste (selon la formule de Lambert), de meme que le bois, la céramique, etc. La quantité de lumière réfléchie dans une direction en fonction de la direction d'incidence est appelée réflectance.

Pour une surface qui réagit à un champ électromagnétique incident il y a deux réflectances en fonction de la polarité de la lumière. La polarité c'est l'orientation dans lesquelles varient les champs électriques et magnétiques de chaque photon. On sépare les photons en orientation horizontale par rapport à la surface et les photons en orientation perpendiculaire à l'horizontale (pas nécessairement verticale). Pour une lumière non polarisée de base, il y a moitié moitié. Notre calcul se fiche de la polarité et fait une simple moyenne de la réflectance calculée pour les photons horizontaux et les photons perpendiculaires.

On doit modifier légèrement le code qui calcule l'éclairage et le reste. D'abord si on a un coéfficient de réflexion et de transmission tous les deux égaux à un, on se fiche du calcul d'éclairage classique, on considère que toute lumière qui arrive est soit transmise soit réflechie. Ensuite on doit décider de la suite à donner à notre algorithme de raytracing. Avant on n'avait un choix limité, soit on terminait la course, soit on réfléchissait le rayon et on continuait le raytracing à partir du point courant. Ici un troisième choix s'offre à nous qui est de continuer la course du rayon, mais cette fois ci de l'autre coté de la surface, le rayon transmis. On peut traiter cela récursivement, en utilisant la pile du C++, ou bien décider d'une méthode alternative qui est de ne suivre qu'un seul choix par rayon. Comment prendre en compte les deux effets dans ce cas ? J'ai pris le parti d'implémenter ici une méthode aléatoire dite de "roulette russe". Chaque rayon suivra une piste unique et a chaque branchement lancera un dé pour savoir quel chemin il va prendre. Le choix de la roulette russe fait que à partir d'un point de départ proche les rayons peuvent diverger fortement en intensité et donnera un effet fortement bruité à notre image. Pour éviter cet effet bruité on est censé faire plus d'un rayon par pixel pour prendre en compte les nombreuses variations de parcours. Cela tombe bien on faisait déjà cela pour rendre l'effet de flou de la profondeur de champ. On peut combiner aisément les deux.

Pour implanter cela correctement, il faut utiliser ce qu'on appelle importance sampling. C'est un principe similaire à celui qui pourrait arriver dans la réalité si l'on suivait le chemin individuel de chaque photon (et le raytracing est vraiment le calcul du chemin inverse du photon). Les coefficients calculés ci dessus, sont donc des probabilités qu'un photon décide de partir dans une direction ou dans une autre.

Transcrit en C++, cela donne:

Comme pour la profondeur de champ, il faut envoyer un certain nombre de rayons pour éviter l'effet bruité, voilà ce que ça donne avec un blob et une densité de réfraction de 2.

Voici le résultat du programme avec les notions abordées dans cette page :

Téléchargement du code source ici : raytrace_page4.zip Vous n'avez pas besoin de librairie particulière pour le compiler et le faire tourner, seulement un compilateur C++ livré avec la standard C++ library. Testé avec le compilateur de GCC 3.2 et visual C++ .net 2003.

Direction page 5 : "HDR, loi de Beer et aberration chromatique".

Quick Navigation :

page 1 : "Premiers pas".
page 2 : "Matériaux spéculaires et post processing".
page 3 : "Textures".
page 4 : "Flou Fresnel et Blobs".
page 5 : "HDR loi de beer et aberration chromatique".
page 6 : "Photon mapping".

Plus d'articles et commentaires : Retour au journal.