Illumination globale, Photon mapping
L'un des problèmes les plus compliqués en synthèse d'image est le problème de l'illumination globale. Il n'est pas difficile en soi à comprendre mais plutôt à mettre en oeuvre en un temps raisonnable. Dans la réalité l'éclairage ne va pas en ligne droite et n'est pas régie par la relation binaire "éclairé/pas éclairé". En pratique, la lumière rebondit, elle s'introduit dans les plus petits interstices, les ombres des objets ont une pénombre et ont des frontières floues. En pratique un objet ou une surface n'est pas éclairée par une source de lumière mais par une infinité de sources de lumière parce que tous les objets de la scène reçoivent de la lumière et en émettent à leur tour.
Même si l'on voulait une simulation parfaite de tous ces phénomènes on n'en aurait pas le temps de notre vivant parce que la puissance nécessaire pour calculer les chemins individuels des photons serait trop grande. On se satisfait donc de simulations partielles. Je vais proposer ici l'une des plus populaires qui est le photon mapping. Ce qui suit est un résumé/relecture du livre de Henrik Wann Jensen et j'encourage toute personne sérieuse sur le sujet du photon mapping et de l'illumination globale de se le procurer. Une autre solution partielle est la radiosité, voir le site de François Sillion pour plus d'informations. http://artis.imag.fr/%7EFrancois.Sillion/
C'est la sixième partie de notre série d'articles sur le raytracing en C++. Elle fait suite à la partie intitulée "HDR, loi de Beer et aberration chromatique".
Le lancer de photons
La première étape du photon mapping consiste à lancer des photons dans la scène. Le lancer de photons (photon tracing) est comme le lancer de rayons, sauf qu'il effectue le trajet de la lumière vers les objets à éclairer. En optique classique le chemin de la lumière est strictement symétrique ce qui explique qu'une bonne partie du code est identique avec les parties précédentes.
La différence se fait lorsque le photon intersecte une surface. À chaque intersection on décide d'enregistrer le photon avec son point d'intersection dans une structure spéciale appelée photon map (carte à photon). C'est cette photon map qui nous resservira plus tard pour notre calcul d'éclairage lors de notre phase de raytracing. Le photon tracing n'est donc qu'une phase de pré-traitement dont le seul but est d'améliorer notre calcul d'éclairage dans la partie suivante.
Notre algorithme a quelques limitations par rapport à une solution complète de photon mapping. En premier lieu, on ne stocke un photon qu'après qu'ils ont subi une réfraction ou une réflexion simple. Notre solution est donc moins globale que ce que j'ai annoncé en préambule mais ça ne change pas plus que ça le principe et les explications qui suivent. Selon le type de scène que l'on veut rendre il peut paraître plus important d'inclure les rebonds indirects après diffusion. Notre code va principalement rendre des caustiques.
* La fonction de lancer des photons est identique ou presque
* à celle de lancer de rayon. La différence est le traitement
* à effectuer en cas d'intersection avec un matériel diffus
* on ne calcule pas l'éclairage immédiatement
* mais on stocke le photon dans une structure hiérarchisé
* pour les retrouver plus facilement dans la phase d'envoi des rayons
*
*/
bool addPhotonRay(ray viewRay, scene &myScene, context myContext, photonmap& myPhotonMap)
{
int level = 0;
do {
point ptHitPoint;
vecteur vNormal;
material currentMat;
{
float distance;
if (!findIntersection(/*IN*/viewRay, /*IN*/myScene,
/*OUT*/ptHitPoint, /*OUT*/vNormal, /*OUT*/currentMat, /*OUT*/distance))
{
break;
}
if (myContext.absorption != 1.0f)
{
// la lumière transmise dans le médium transparent
// est affectée par la diffusion due aux impuretés
// la beer's law s'applique pour les faibles concentrations
float coef = powf(myContext.absorption, distance);
float fRoulette = (1.0f / RAND_MAX) * rand();
if (fRoulette > coef )
{
// Le photon s'est retrouvé absorbé avant d'arriver au point d'intersection.
// On peut juste s'arreter de stocker des photons
break;
}
}
}
bool bInside;
if (vNormal * viewRay.dir > 0.0f)
{
vNormal = -1.0f * vNormal;
bInside = true;
}
else
{
bInside = false;
}
// perturbation du vecteur normal par un bruit aléatoire
// on le renormalise derrière pour le conserver unitaire
// pour les calculs d'éclairage
if (currentMat.bump != 0.0f)
{
float noiseCoefx = noisef(0.1f * ptHitPoint.x, 0.1f * ptHitPoint.y,0.1f * ptHitPoint.z);
float noiseCoefy = noisef(0.1f * ptHitPoint.y, 0.1f * ptHitPoint.z,0.1f * ptHitPoint.x);
float noiseCoefz = noisef(0.1f * ptHitPoint.z, 0.1f * ptHitPoint.x,0.1f * ptHitPoint.y);
vNormal.x = (1.0f - currentMat.bump ) * vNormal.x + currentMat.bump * noiseCoefx;
vNormal.y = (1.0f - currentMat.bump ) * vNormal.y + currentMat.bump * noiseCoefy;
vNormal.z = (1.0f - currentMat.bump ) * vNormal.z + currentMat.bump * noiseCoefz;
float temp = vNormal * vNormal;
if (temp == 0.0f)
break;
vNormal = invsqrtf(temp) * vNormal;
}
const materialChannel & currentMatChannel = currentMat.getChannel(myContext.offset);
ray outRay = viewRay;
bool bDiffuse = LookForDiffuse(viewRay, ptHitPoint, vNormal, currentMatChannel, bInside, myContext, outRay);
if (bDiffuse)
{
// notre BRDF nous dit que le photon est soit réflechi soit réfracté
// soit absorbé par la surface diffuse..
// Pas de color bleeding dû à la diffusion.. pour l'instant
if (myScene.phMapProp.directLighting || (level > 0) )
{
// La surface n'est pas totalement transitive, on peut stocker le photon incident
// on ne le fait que pour les niveaux supérieurs à zero si l'on décide
// d'utiliser les formules classiques pour l'éclairage direct
photon currentPhoton; // {ptHitPoint, viewRay.dir, coef};
currentPhoton.pos = ptHitPoint;
currentPhoton.dir = viewRay.dir;
if (!myPhotonMap.addPhoton(currentPhoton, vNormal))
{
// Le container a photon a atteint ses limites.
// On peut juste s'arreter de stocker des photons
return false;
};
}
return true;
}
viewRay = outRay;
level++;
} while (level < 10);
return true;
}
On notera la subtile différence du traitement de l'absorption. Notre choix a été de décimer les photons en proportion de la quantité absorbée. En contrepartie notre carte de photon ne stockera que des photons d'intensité égale, cela fait un nombre de moins à stocker par photon. On lance des photons jusqu'à ce qu'il n'y ait plus de place dans le conteneur à photon.
Échantillonnage
Bon on sait comment lancer un photon mais comment évaluer le nombre, le point de départ et la direction de photons à envoyer pour notre lumière ? On va encore faire de l'échantillonnage. Choisir une direction de manière aléatoire est équivalent à prendre des points aléatoirement à l'intérieur d'une sphère à l'aide d'un crible. On utilise la distribution uniforme de la fonction rand, et l'on applique la post condition que le point se trouve à l'intérieur de la sphère (criblage). Ensuite on normalise le vecteur obtenu pour avoir une direction parfaitement aléatoire (autant que la fonction rand utilisée le permet). On aurait pu penser à d'autres solutions impliquant le choix d'angles sphériques théta et phi aléatoires mais elles ne sont pas aussi bonnes parce que ces angles n'ont pas une distribution uniforme sur la sphère.
for (;;)
{
ray viewRay;
do
{
viewRay.dir.x = (2.0f / RAND_MAX) * rand() - 1.0f;
viewRay.dir.y = (2.0f / RAND_MAX) * rand() - 1.0f;
viewRay.dir.z = (2.0f / RAND_MAX) * rand() - 1.0f;
}
while (viewRay.dir.x *viewRay.dir.x + viewRay.dir.y * viewRay.dir.y + viewRay.dir.z * viewRay.dir.z > 1.0f );
float temp = viewRay.dir * viewRay.dir;
if (temp == 0.0f)
continue;
viewRay.dir = invsqrtf(temp) * viewRay.dir;
int lightNumber = rand() % int(myScene.lgtTab.size());
viewRay.start = myScene.lgtTab[lightNumber].pos;
myContext.offset = color::OFFSET(offset);
if (!addPhotonRay(viewRay, myScene, myContext, channelPhotonMap[offset]))
{
// on a atteint la limite du stockage de photons, on peut arreter d'en ajouter
break;
};
nPhotonSent++;
}
Le nombre de photons stockés est limité par contre le nombre de photons envoyés ne l'est pas ce qui peut-etre un problème si l'on ne fait pas attention au coût global. Dans l'idéal on ne devrait envoyer des photons que dans les directions importantes ou dans celles qui vont aboutir à des photons visibles, ce qui peut-etre non trivial à determiner de par la nature même du phénomène que l'on cherche à simuler.
L'intensité de chaque photon individuel est calculé par la formule L'intensité de chaque photon individuel est calculé par la formule i = Intensité totale / nPhotonSent. Il est important de diviser par le nombre de photons envoyés et pas par le nombre de photons stockés, parce que les photons de l'éclairage direct ou ceux qui ont été décimés contribuent aussi à leur manière à l'éclairage global.
Conteneur à photons
Le stockage de nos photons se fait dans un simple std::vector, mais avec une nuance importante. Les photons sont adressables depuis une deuxième structure à base de KD-Tree. Un arbre KD est une structure d'arbre binaire qui permet de subdiviser l'espace en sous espaces de dimension équivalente. En l'occurrence ici la dimension est le nombre de photons présent dans une de ses branches. L'intérêt du KD tree est de pouvoir retrouver plus rapidement les m photons les plus proches d'un point donné sans avoir à parcourir toute la structure. Déjà il est trivial de retrouver LE plus proche, les m plus proches sont un peu moins triviaux parce qu'il faut parcourir récursivement plusieurs branches mais cela peut se faire assez rapidement avec du pruning (élagage de branches) basées sur la distance maximale.
La construction du KD tree se fait de la manière suivante. En premier lieu, on cherche le photon médian par rapport à l'axe choisi (X, Y ou Z), ce photon médian est trouvé en faisant le tri de nos photons selon la coordinée choisie et en prenant celui du milieu (en indice et non pas en distance sur l'axe). On place ce photon médian dans le premier noeud de l'arbre. Le choix de l'axe pourrait se faire en séquence (X puis Y puis Z puis on répète) mais en pratique il est plus intéressant de choisir l'axe pour lequel notre ensemble couvre la plus grande distance (si les photons sont tous sur un même plan, une subdivision selon ce plan n'apporte que peu d'information supplémentaire). Les derniers photons occupent un noeud sans fils.
C'est fait par la fonction balance(). C'est un peu fastidieux parce qu'il faut énumérer les cas mais la base est simple :
{
const photon *photonTab;
bool operator() (int p1, int p2)
{
if (n == kdnode::dimensionX)
{
return photonTab[p1].pos.x < photonTab[p2].pos.x;
}
if (n == kdnode::dimensionY)
{
return photonTab[p1].pos.y < photonTab[p2].pos.y;
}
if (n == kdnode::dimensionZ)
{
return photonTab[p1].pos.z < photonTab[p2].pos.z;
}
}
};
void photoncache::balance(vector<int>& tab, int p)
{
if (tab.size() == 0)
return;
vector<int> tabLess;
vector<int> tabMore;
if (2 * p + 1 < 100 * int(photonTab.size()))
{
kdnode node;
node.nPhotonMedian = -1;
if (tab.size() == 1)
{
node.nPhotonMedian = tab[0];
}
else
{
float minX, minY, minZ;
float maxX, maxY, maxZ;
minX = minY = minZ = numeric_limits<float>::max();
maxX = maxY = maxZ = - numeric_limits<float>::max();
tabLess.reserve(tab.size());
tabMore.reserve(tab.size());
for (int i = 0; i < int(tab.size()); ++i)
{
minX = min(minX, photonTab[tab[i]].pos.x);
minY = min(minY, photonTab[tab[i]].pos.y);
minZ = min(minZ, photonTab[tab[i]].pos.z);
maxX = max(maxX, photonTab[tab[i]].pos.x);
maxY = max(maxY, photonTab[tab[i]].pos.y);
maxZ = max(maxZ, photonTab[tab[i]].pos.z);
}
if ( maxX - minX >= maxY - minY )
{
if (maxX - minX >= maxZ - minZ)
{
node.dimensionPlane = kdnode::dimensionX;
sortingCriterion<kdnode::dimensionX> pred = {&photonTab[0]};
sort(tab.begin(), tab.end(), pred);
node.nPhotonMedian = tab[int(tab.size()) / 2];
for (int i = 0; i< int(tab.size()) / 2; ++i)
{
tabLess.push_back(tab[i]);
}
for (int i = int(tab.size()) / 2 + 1; i< int(tab.size()); ++i)
{
tabMore.push_back(tab[i]);
}
}
else
{
node.dimensionPlane = kdnode::dimensionZ;
sortingCriterion<kdnode::dimensionZ> pred = {&photonTab[0]};
sort(tab.begin(), tab.end(), pred);
node.nPhotonMedian = tab[int(tab.size()) / 2];
for (int i = 0; i< int(tab.size()) / 2; ++i)
{
tabLess.push_back(tab[i]);
}
for (int i = int(tab.size()) / 2 + 1; i< int(tab.size()); ++i)
{
tabMore.push_back(tab[i]);
}
}
}
else
{
if (maxY - minY >= maxZ - minZ)
{
node.dimensionPlane = kdnode::dimensionY;
sortingCriterion<kdnode::dimensionY> pred = {&photonTab[0]};
sort(tab.begin(), tab.end(), pred);
node.nPhotonMedian = tab[int(tab.size()) / 2];
for (int i = 0; i< int(tab.size()) / 2; ++i)
{
tabLess.push_back(tab[i]);
}
for (int i = int(tab.size()) / 2 + 1; i< int(tab.size()); ++i)
{
tabMore.push_back(tab[i]);
}
}
else
{
node.dimensionPlane = kdnode::dimensionZ;
sortingCriterion<kdnode::dimensionZ> pred = {&photonTab[0]};
sort(tab.begin(), tab.end(), pred);
node.nPhotonMedian = tab[int(tab.size()) / 2];
for (int i = 0; i< int(tab.size()) / 2; ++i)
{
tabLess.push_back(tab[i]);
}
for (int i = int(tab.size()) / 2 + 1; i< int(tab.size()); ++i)
{
tabMore.push_back(tab[i]);
}
}
}
}
if (p >= int(nodeTab.size()) )
nodeTab.resize(p + 1);
nodeTab[p] = node;
}
balance(tabLess, 2 * p);
balance(tabMore, 2 * p + 1);
}
Comme vous pouvez le voir dans ce listing, il n'y a pas de structure d'arbre à proprement parler, nos noeuds sont stockés dans un std::vector, on s'arrange juste pour que le noeud p ait deux fils 2 p et 2 p +1. Ce type de stockage implicite est intéressant en terme de taille parce que l'arbre est équilibré et donc il y aura en proportion peu de trous dans la structure.
Cela ne parait pas trop optimisé mais le coût de construction de l'arbre par rapport au reste est vraiment négligeable puisqu'on ne le fait qu'une seule fois par scène.
Digression : Pourquoi pas un octree ? Tout simplement parce que nos photons sont regroupés dans une portion congrue de l'espace et qu'un octree est basé sur une subdivision régulière de l'espace et de nombreuses branches seraient donc vides. Pourquoi pas un BSP tree ? Simplement parce qu'un BSP tree n'offre que peu d'avantage par rapport à un KD tree et parce que dans un BSP tree on stocke généralement explicitement les plans de subdivision ce qui augmente d'autant la taille de la structure (et réduit la performance). Il n'y a en plus pas de critère "naturel" pour choisir un bon plan de subdivision donc exit.
Reconstruction
La dernière phase est la reconstruction. Lorsque notre rayon vue heurte une surface, il va falloir identifier la contribution de l'éclairage direct ainsi que celui indirect de nos photons vagabonds. Pour cette dernière phase, on va trouver les m photons les plus proches puis additionner leur contribution pour obtenir l'éclairage indirect.
Comment trouver les m photons les plus proches ? Solution naïve, parcourir notre structure de photons et ajouter les photons en les triant au fur et à mesure, inconvénient cela a un coût linéaire en nombre de photons dans notre structure ce qui n'est pas si intéressant vu le nombre de calculs doit encore être multiplié par le nombre de rayons envoyés. La solution adoptée est de parcourir notre arbre KD en élagant au fur et à mesure les branches qui sont trop lointaines pour pouvoir être considérées. Un noeud est constitué d'un photon et d'un plan de subdivision de l'espace en deux parties de chaque côté du photon médian. Tout objet situé de l'autre côté de ce plan (par rapport à notre point d'intersection) est donc à plus grande distance que le photon médian. On peut donc se permettre d'élaguer si notre photon médian est plus éloigné que les m photons déjà trouvés. On peut faire mieux, comme l'influence de nos photons est une fonction décroissante de leur distance au point d'intersection (voir plus loin), il est possible de fournir dès le départ une distance maximale qui permettra d'élaguer plus vite (jump start) sans avoir à trouver m photons d'abord. Notre algorithme utilise également la distance des photons précédents pour se permettre d'élaguer plus aggressivement, ce qui permet d'économiser énormément de parcours d'arbre.
const vecteur &normal,
float fMinCos, int nPhotons,
float & maxDistSquare,
vector<pair<float, int> > & result,
int firstNode)
{
const vecteur *normalTab = &container->normalTab[0];
const photon *photonTab = &container->photonTab[0];
if (2 * firstNode + 1 < int(container->nodeTab.size()))
{
int nPhotonMedian = container->nodeTab[firstNode].nPhotonMedian ;
if (nPhotonMedian >= 0)
{
// Il y a un photon à ce noeud, evalue la distance au point d'intersection
// et stocke le s'il y a lieu
const photon * currentPhoton = &container->photonTab[nPhotonMedian];
{
vecteur vDistToPhoton = pos - currentPhoton->pos;
float fDistSquare = vDistToPhoton * vDistToPhoton;
// On ne rajoute pas un photon qui se retrouve déjà plus loin que la distance max
if ( maxDistSquare > fDistSquare)
{
// On élimine les photons qui ont été stockés sur une surface qui a peu de probabilité
// d'être la même que celle du point d'intersection. Sinon il y a risque de pollution entre
// deux surfaces proches (angle de deux murs par exemple)
if ( (normalTab[nPhotonMedian] * normal > fMinCos) && (photonTab[nPhotonMedian].dir * normal < 0.0f) )
{
if (int(result.size()) == nPhotons )
{
// le conteneur est déjà plein. Mais ça ne veut pas dire que
// le nouveau photon n'est pas plus proche que les n autres
// on compare les distance et on stocke le cas échéant
// en décalant les autres.
maxDistSquare = 0;
for (int i = 0; i < nPhotons; ++i)
{
if(result[i].first > fDistSquare)
{
float fTemp = fDistSquare;
fDistSquare = result[i].first;
result[i].first = fTemp;
int temp = nPhotonMedian;
nPhotonMedian = result[i].second;
result[i].second = temp;
}
if (result[i].first > maxDistSquare)
{
maxDistSquare = result[i].first;
}
}
}
else
{
// On ajoute à la liste qui n'est pas encore pleine
result.push_back(pair<float, int>(fDistSquare,nPhotonMedian));
}
}
}
}
float fDist;
switch (container->nodeTab[firstNode].dimensionPlane)
{
case kdnode::dimensionX:
fDist = pos.x - currentPhoton->pos.x;
break;
case kdnode::dimensionY:
fDist = pos.y - currentPhoton->pos.y;
break;
case kdnode::dimensionZ:
default:
fDist = pos.z - currentPhoton->pos.z;
break;
}
if (fDist <= 0.0f)
{
float fDistSquare = fDist * fDist;
locatePhotons(pos, normal, fMinCos, nPhotons, maxDistSquare, result, 2 * firstNode);
if ( maxDistSquare >= fDistSquare )
{
// Tout l'intérêt du KD-Tree réside dans le pruning
// si l'on a déjà trouvé n photons, pas la peine de chercher dans la liste
// des photons qui se trouvent plus loin que la distance max.
locatePhotons(pos, normal, fMinCos, nPhotons, maxDistSquare, result, 2 * firstNode + 1);
}
}
else
{
float fDistSquare = fDist * fDist;
locatePhotons(pos, normal, fMinCos, nPhotons, maxDistSquare, result, 2 * firstNode + 1);
if ( maxDistSquare >= fDistSquare )
{
// Tout l'intérêt du KD-Tree réside dans le pruning
// si l'on a déjà trouvé n photons, pas la peine de chercher dans la liste
// des photons qui se trouvent plus loin que la distance max.
locatePhotons(pos, normal, fMinCos, nPhotons, maxDistSquare, result, 2 * firstNode);
}
}
}
}
}
Une fois ces m photons trouvés il faut calculer leur contribution. On évalue la force de ces photons en utilisant une fonction d'évaluation basée sur la surface couverte par les photons, et leur distance au point d'intersection. Il y a plusieurs évaluateurs possibles. On en décrit trois. Le premier utilise une formule constante pour tous les photons et divise leur intensité par la surface couverte (on prend le photon le plus lointain et on divise par pi * distance au carré.). Cela marche bien pour les calculs d'éclairage doux comme les interreflexions on l'appelle "radiance estimator" (à défaut d'autre nom). Le deuxième utilise donne un poids plus élevé aux photons plus proche selon une fonction conique de leur distance au point d'intersection. Cela permet d'améliorer le rendu de variations rapides d'intensités. Le dernier utilise une fonction type "gaussienne" normalisée.
Le calcul d'éclairage proprement dit est quasiment identique. Dans l'éclairage direct on doit prendre en compte dans les calculs le facteur de Lambert qui fait que plus l'angle du rayon lumineux est rasant, plus un nombre de photons donné couvrira une surface importante (intensité par unité de surface est donc réduite). Dans le modèle par photon mapping cela doit en principe être simulé par la dispersion des photons et la division de la contribution totale par la surface couverte.
if ( found > 0 && nearestPhoton.maxDistSquare > 0.0f)
{
float fInvMaxDist = invsqrtf(nearestPhoton.maxDistSquare);
vector<pair<float, const photon * > >::const_iterator it = nearestPhoton.resultTab.begin();
for (;it != nearestPhoton.resultTab.end() ; ++it)
{
float fPhotonDist = sqrtf((*it).first);
float fPhotonWeight;
switch (myScene.phMapProp.eFilterType)
{
case photonMapSettings::conicEstimator:
{
static const float invKappa = 0.9f;
fPhotonWeight = 1.0f - fPhotonDist * fInvMaxDist * invKappa;
}
break;
case photonMapSettings::gaussianEstimator:
{
static const float Alpha = 1.818f;
static const float Beta = 1.953f;
float fDistRatio = fPhotonDist * fInvMaxDist;
fPhotonWeight = Alpha *
(1.0f - (1.0f - expf(-0.5f * Beta * fDistRatio * fDistRatio)) / (1.0f - expf (-Beta) ));
}
break;
default: // photonMapSettings::RadianceEstimator:
{
fPhotonWeight = 1.0f;
}
break;
}
float fLightIntensity = fPhotonWeight;
float fLightProjection = - ((*it).second->dir * vNormal) ;
if (fLightProjection > 0.0f)
{
float lambert = 1.0f;
fLightContrib += lambert * fLightIntensity
* ((1.0f - noiseCoef)* currentMatChannel.diffuse +
noiseCoef * currentMatChannel.diffuse2);
// Blinn
// La direction de Blinn est exactement à mi chemin entre le rayon
// lumineux et le rayon de vue.
// On calcule le vecteur de Blinn et on le rend unitaire
// puis on calcule le coéfficient de blinn
// qui est la contribution spéculaire de la lumière courante.
vecteur blinnDir = (*it).second->dir + viewRay.dir;
float temp = blinnDir * blinnDir;
if (temp != 0.0f )
{
float blinn = invsqrtf(temp) * max(fLightProjection - fViewProjection , 0.0f);
blinn = currentMatChannel.specular * powf(blinn, currentMatChannel.power);
fLightContrib += blinn * fLightIntensity;
}
}
}
if(fLightContrib > 0.0f)
{
switch (myScene.phMapProp.eFilterType)
{
case photonMapSettings::conicEstimator:
{
static const float invKappa = 0.9f;
static const float divider = 1.0f / ( (1.0f - 0.666666666667f * invKappa) * 3.1415927f);
fLightContrib *= myPhotonMap.getWattagePerPhoton()
*
divider * fInvMaxDist * fInvMaxDist ;
}
break;
case photonMapSettings::gaussianEstimator:
{
const float divider = 1.0f / 3.1415927f;
fLightContrib *= myPhotonMap.getWattagePerPhoton()
* divider * fInvMaxDist * fInvMaxDist ;
}
break;
default: // photonMapSettings::RadianceEstimator:
{
const float divider = 1.0f / 3.1415927f;
fLightContrib *= myPhotonMap.getWattagePerPhoton()
* divider * fInvMaxDist * fInvMaxDist ; }
break;
}
}
}
Résultats :
Dans l'image suivante on utilise le photon mapping pour l'éclairage direct ce qui est sans doute inutile, notre méthode habituelle est très rapide pour ce genre de problème et résout les équations sans artefacts. Le résultat n'est pas très satisfaisant, à cause d'un fort bruitage dans la partie éclairage direct et les frontières des ombres et des caustiques ne sont pas bien délimitées :
Ici, on utilise le modèle classique pour l'éclairage direct, et le modèle à base de photons pour les caustiques et les interreflexions. On a besoin de beaucoup moins de photons et on peut les concentrer dans les zones d'intérêt :
Enfin, pour la route, les caustiques causées par un blob introduit dans les articles précédents. Le résultat est quasiment abstrait :
Téléchargement du code source ici : raytrace_page6.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.
Conclusion
C'est donc la fin de cette série d'articles d'introduction au raytracing. J'ai essayé d'aborder le problème d'un point de vue un peu plus original que ce qu'on peut trouver sur Internet sans oublier de couvrir certaines bases. En général je considère qu'un bout de code vaut mieux qu'un long discours (les langages de programmation sont fait pour décrire des algorithmes, les langages humains ne le sont pas). Mais si vous avez des soucis pour comprendre certaines parties n'hésitez pas à m'envoyer un email.
Je ne doute pas non plus que le code est bourré d'erreurs ou que mes articles ont une orthographe qui n'est pas irréprochable donc si vous avez des suggestions, sur la qualité des explications, les fautes de grammaire et/ou d'éventuels bugs dans les programmes fournis n'hésitez pas à me contacter j'essaierai de corriger aussi vite que possible.
Signé
: LeGreg <
>
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.