<figure><img src="image_1.png" /></figure>

HAL
open science

# Rendu Temps-Réel de Nuages à l'Echelle Planétaire

Laurent Belcour

► To cite this version:

Laurent Belcour. Rendu Temps-Réel de Nuages à l'Echelle Planétaire. Graphics [cs.GR]. 2009. hal-04138703

HAL Id: hal-04138703

https://inria.hal.science/hal-04138703

Submitted on 23 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rendu Temps-Réel de Nuages à l'Echelle Planétaire

<figure><img src="image_3.png" /></figure>

Laurent Belcour sous la direction d'Eric Bruneton

Rapport de stage de Master de Recherche en Informatique Graphique, Vision et Robotique
Imprimé le 11 juin 2009
<figure><img src="image_4.png" /></figure>
Remerciements

Je tiens à remercier toute l'équipe d'Evasion pour son acceuil chaleureux et sa bonne ambiance et plus particulièrement M. Eric Bruneton.

Merci à toutes les personnes qui ont bien voulu relire ce travail permettant de le parfaire grandement et à tous ceux qui n'y ont jeté qu'un oeil mais dont l'intérêt ou les remarques m'ont permis d'avancer.

Laurent Belcour
<figure><img src="image_5.png" /></figure>
Table des matières

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Introduction</td><td>9</td></tr><tr><td>1.1</td><td>L'équipe Evasion</td><td>9</td></tr><tr><td>1.2</td><td>ProLand</td><td>9</td></tr><tr><td>1.3</td><td>Le Sujet</td><td>10</td></tr><tr><td>1.3.1</td><td>Travail d'Antoine Bouthors</td><td>10</td></tr><tr><td>1.3.2</td><td>Travail de Master</td><td>11</td></tr><tr><td>1.4</td><td>Structure du Rapport</td><td>11</td></tr><tr><td>2</td><td>Physique de la Lumière</td><td>13</td></tr><tr><td>2.1</td><td>La Lumière</td><td>13</td></tr><tr><td>2.1.1</td><td>La Radiance</td><td>13</td></tr><tr><td>2.1.2</td><td>L'Irradiance</td><td>14</td></tr><tr><td>2.1.3</td><td>L'Intensité Lumineuse</td><td>14</td></tr><tr><td>2.1.4</td><td>Contact avec un Corps</td><td>14</td></tr><tr><td>2.1.5</td><td>L'Albédo</td><td>16</td></tr><tr><td>2.2</td><td>L'Equation de Transfert de la Lumière</td><td>16</td></tr><tr><td>2.2.1</td><td>Réflection sur un Support Solide</td><td>16</td></tr><tr><td>2.2.2</td><td>Diffusion dans un Milieu Particulaire</td><td>17</td></tr><tr><td>2.2.3</td><td>Le Calcul de Radiosité</td><td>18</td></tr><tr><td>2.3</td><td>Interactions Lumineuses dans les Nuages</td><td>18</td></tr><tr><td>2.3.1</td><td>Les Nuages</td><td>18</td></tr><tr><td>2.3.2</td><td>La Dispersion de Mie</td><td>18</td></tr><tr><td>2.3.3</td><td>Transparence</td><td>20</td></tr><tr><td>2.3.4</td><td>L'influence des Nuages sur le Sol</td><td>21</td></tr><tr><td>3</td><td>État de l'Art</td><td>23</td></tr><tr><td>3.1</td><td>La Fonction de Phase</td><td>23</td></tr><tr><td>3.1.1</td><td>Utilisation d'une autre Fonction de Phase</td><td>23</td></tr><tr><td>3.1.2</td><td>Approximation de Mie</td><td>23</td></tr><tr><td>3.1.3</td><td>Utilisation de Mie</td><td>24</td></tr><tr><td>3.2</td><td>Le Scattering et l'Importance des Réflexions</td><td>24</td></tr><tr><td>3.2.1</td><td>Calcul Analytique d'un Scatter</td><td>25</td></tr><tr><td>3.2.2</td><td>Simplification du Nombre de Scatters</td><td>26</td></tr><tr><td>3.2.3</td><td>Calcul Analytique de Paquets de Scatter</td><td>26</td></tr><tr><td>3.2.4</td><td>Calculs non Analytique de Paquets de Scatter</td><td>27</td></tr><tr><td>3.3</td><td>Interaction entre Sol et Nuages</td><td>27</td></tr><tr><td>3.3.1</td><td>Shadow maps</td><td>27</td></tr></tbody></table>
<table><tbody><tr><td>3.3.2 Rayons de lumière</td><td>27</td></tr><tr><td>3.3.3 Calcul de radiosité</td><td>28</td></tr><tr><td>3.4 Transferts Radiatifs Temps-réel</td><td>28</td></tr><tr><td>3.4.1 Transferts Radiatifs Précalculés</td><td>28</td></tr><tr><td>3.4.2 Bases Sphériques</td><td>30</td></tr><tr><td>3.4.3 Lumières multiples</td><td>31</td></tr><tr><td>3.5 Bilan</td><td>32</td></tr><tr><td><strong>4 D'un Modèle Plan à une Terre Sphérique</strong></td><td><strong>33</strong></td></tr><tr><td>4.1 Quelques Définitions Utiles</td><td>33</td></tr><tr><td>4.2 Chemins Lumineux et Hauteurs Apparentes</td><td>34</td></tr><tr><td>4.2.1 Depuis la Caméra, au-Dessus des Nuages</td><td>34</td></tr><tr><td>4.2.2 Depuis la Caméra, sous la Couche Nuageuse</td><td>36</td></tr><tr><td>4.2.3 Depuis le Soleil</td><td>36</td></tr><tr><td>4.2.4 Comparaison avec le Modèle Initial</td><td>37</td></tr><tr><td>4.2.5 Conclusion</td><td>37</td></tr><tr><td>4.3 Modification du Modèle Lumineux</td><td>38</td></tr><tr><td>4.3.1 Illumination Directe</td><td>38</td></tr><tr><td>4.3.2 Single Scattering</td><td>38</td></tr><tr><td>4.3.3 Double Scattering</td><td>40</td></tr><tr><td>4.3.4 Multiple Scattering</td><td>40</td></tr><tr><td>4.3.5 Conclusion</td><td>40</td></tr><tr><td><strong>5 Les Transferts Radiatifs entre la Terre et les Nuages</strong></td><td><strong>43</strong></td></tr><tr><td>5.1 Algorithme de Bouthors et al.</td><td>43</td></tr><tr><td>5.1.1 Description de l'Algorithme</td><td>43</td></tr><tr><td>5.1.2 Restriction de l'Intégrale à un Disque</td><td>44</td></tr><tr><td>5.1.3 Comparaison avec un Sol non Plat</td><td>44</td></tr><tr><td>5.1.4 Conclusion sur les Calculs d'Erreur</td><td>47</td></tr><tr><td>5.2 Notre Algorithme</td><td>48</td></tr><tr><td>5.2.1 Simplification du Calcul Intégral</td><td>48</td></tr><tr><td>5.2.2 Choix des bases {C<sub>i</sub>}<sub>i∈N</sub> et {G<sub>i</sub>}<sub>i∈N</sub></td><td>52</td></tr><tr><td>5.2.3 Pseudo-Code</td><td>55</td></tr><tr><td>5.2.4 Validation de la Solution Proposée</td><td>59</td></tr><tr><td><strong>6 Résultats & Conclusion</strong></td><td><strong>61</strong></td></tr><tr><td>6.1 Résultats</td><td>61</td></tr><tr><td>6.1.1 Illumination des Nuages dans ProLand</td><td>61</td></tr><tr><td>6.1.2 Radiosité entre la Terre et les Nuages dans ProLand</td><td>62</td></tr><tr><td>6.2 Conclusion</td><td>62</td></tr><tr><td>6.2.1 Sur le Rendu de Nuages</td><td>62</td></tr><tr><td>6.2.2 Sur le Précalcul de Radiance</td><td>63</td></tr><tr><td><strong>A Annexes</strong></td><td><strong>65</strong></td></tr><tr><td>A.1 Erreur dans la Restriction des Calculs de Rebonds</td><td>65</td></tr><tr><td>A.2 Validation du Précalcul</td><td>68</td></tr></tbody></table>
<figure><img src="image_6.png" /></figure>
1

Introduction

Avant de commencer tout calcul, avant toute considération informatique et avant toute
formule mathématique, nous allons présenter le contexte du stage. Cette courte présenta-
tion permettra de situer le travail de ce master dans le projet auquel il appartient. Tout
d'abord, nous présenterons l'équipe dans laquelle a été effectué le stage : *Evasion*. Ensuite,
nous présenterons *ProLand*, l'environnement scientifique et logiciel dans lequel ce travail
s'insère. Pour finir cette introduction, nous présenterons le sujet du stage de Master et la
façon dont sera structuré le document (ainsi que le travail de recherche effectué).

1.1 L'équipe Evasion

Evasion est une équipe du Laboratoire d'informatique Jean Kuntzmann et de l'INRIA Rhône-Alpes. Le thème de cohésion de cette équipe est la modélisation, l'animation et la visualisation de phénomènes et objets naturels. Cette recherche comprend entre autre l'animation physique d'organes, la capture de mouvements d'êtres vivants, le rendu de données naturelles massives, etc.

Ce stage de master à été réalisé au sein d'Evasion, dans le contexte du projet ProLand dirigé par Eric Bruneton.

1.2 ProLand

ProLand est un projet de recherche visant à permettre la navigation temps-réel dans un monde constitué de données massives pouvant être générées ou modifiées procéduralement (ProLand étant la contraction de *Procedural Landscape*). Le but de ce projet est de permettre le rendu et l'édition en temps réel d'environnements très vastes (comme la Terre) à toutes les échelles de grandeur (du brin d'herbe à l'échelle spatiale).

La figure 1.1 présente deux captures issues d'un logiciel de ProLand utilisant des données satellites (élévation, couleur) de la planète Terre. Ces données proviennent de
<figure><img src="image_7.png" /><figcaption>**FIGURE 1.1** – Quelques résultats logiciels du projet Proland présentant la Terre visualisée depuis l'espace (à gauche) et du point de vue d'un avion de tourisme (à droite).</figcaption></figure>

relevés de la NASA. On peut noter que de nombreux effets atmosphériques sont rendus
car ProLand se nourrit des travaux de recherche d'Eric Bruneton et de Fabrice Neyret
notamment sur l'atmosphère [BN08].

1.3 Le Sujet

Le sujet de ce travail de master est le *rendu temps-réel de nuages à l'échelle planétaire*.
Le but de ce projet de master est d'étudier et de proposer une solution pour le rendu de
nuages sur un globe terrestre complet. A notre connaissance, une solution de rendu de
nuage à cette échelle n'a jamais été proposée. Il suit le travail de thèse d'Antoine Bouthors
sur le rendu de nuage dans l'équipe Evasion.

1.3.1 Travail d'Antoine Bouthors

Durant sa thèse dans l'équipe Evasion, Antoine Bouthors a travaillé sur le rendu
photo-réaliste de nuages. Sa thèse présente deux solutions pour le rendu de nuages :

La première "Real-time realistic illumination and shading of stratiform clouds" [BNL06] porte sur le rendu de nuages stratiformes. Elle permet un rendu convaincant des nuages (voir la figure 1.2, image de gauche) et de très bonnes performances (une centaine d'images par seconde). Les nuages sont supposés plans, tout comme le sol (pas de montagnes) ce qui constitue la majeure restriction de ce modèle.

La seconde "Interactive multiple anisotropic scattering in clouds" [BNM+08] est une méthode de rendu pour n'importe quel type de nuage. Elle permet un rendu photo-réaliste (voir la figure 1.2, image de droite), malheureusement elle est très coûteuse en calculs et n'est pas adaptée pour une scène ayant une échelle plus large que celle d'un seul nuage.
<figure><img src="image_8.png" /><figcaption>FIGURE 1.2 – Résultats des solutions de rendu de nuage de Bouthors et al.. A gauche nous présentons une image issue de Real-time realistic illumination and shading of stratiform clouds[BNL06]. A droite nous présentons une image issue de Interactive multiple anisotropic scattering in clouds[BNM+08].</figcaption></figure>

1.3.2 Travail de Master

Nous avons décidé d'étendre le modèle de rendu de nuages stratiformes de Bouthors et al. [BNL06] pour prendre en compte la forme non plane de la Terre (nous détaillerons notre choix dans la partie 3). Étant donné que cette méthode fait l'hypothèse d'un sol plat, nous avons du réétudier l'algorithme de radiosité entre le bas de la couche nuageuse et le sol. Nous avons donc séparé le travail de master en deux parties :

Modification du modèle plan du nuage en un modèle sphérique, dans laquelle nous avons cherché les limitations que la restriction planaire engendrait dans le calcul des interactions lumineuses à l'intérieur du nuage.

Étude du modèle de radiosité, qui porte sur l'étude des échanges lumineux entre le bas du nuage et le sol, quand ce dernier n'est pas plat.

1.4 Structure du Rapport

Le présent rapport de master est divisé en 5 parties.

La première partie présentera les phénomènes physiques associés au transport de la lumière dans les gaz ainsi que la réflection de la lumière sur des supports. Nous développerons les équations nécessaires au calcul de l'illumination d'une scène composée de nuages et d'éléments du sol.

Ensuite, nous présenterons les publications précédant ce travail. Nous détaillerons les méthodes qui nous ont permis de mener le travail de recherche à bien. Nous présenterons les avantages et inconvénients de chaque publication en la mettant dans le contexte de notre travail.

La troisième partie sera l'occasion pour nous de présenter le travail de validation de l'algorithme de rendu de Bouthors *et al*. dans un cadre plus général que celui proposé par les auteurs (Terre sphérique plutôt que plate).

Dans la quatrième partie, nous ferons le développement de notre contribution : un algorithme d'illumination globale temps réel avec précalcul de transferts radiatifs. Nous détaillerons en premier lieu les fondements mathématiques de notre algorithme, puis nous
proposerons un pseudo-code et pour finir nous validerons notre modèle par rapport à un simulateur.

Pour finir, nous présenterons nos résultats dans le cadre du projet ProLand, nous comparerons notre méthode à l' existant dans ProLand. Pour clore notre rapport, nous donnerons des pistes quant au travail de recherche possible sur le rendu, l'animation et la modélisation de nuage dans les milieux naturels.
"Darkness isn't the opposite
of light, it is simply its ab-
sence"

Terry Pratchett - The Light
Fantastic

2

Physique de la Lumière

Décrire l’illumination d’une scène quelconque nécessite une connaissance des réactions physiques que la lumière va avoir sur les différents objets de la scène. Pour obtenir une telle connaissance, nous devons nous intéresser de plus près à la physique de la lumière et au transport de cette dernière dans les matériaux.

2.1 La Lumière

On appelle *lumière* l'ensemble des ondes électromagnétiques perceptibles par l'oeil. Les longueurs d'ondes de la lumière vont de 380nm à 780nm. La lumière peut aussi être considérée comme un flux de particules, les *photons*. Dans la suite du document, nous utiliserons le modèle le plus simple selon le contexte, cela étant permis par la dualité onde-particule<fnref n="1" />.

De notre point de vue macroscopique, il nous est impossible de mesurer l'action d'un
seul photon sur un médium (qu'il soit solide ou gazeux). Nous devons donc définir des
grandeurs mesurant l'influence de la lumière sur des surfaces ou des volumes.

2.1.1 La Radiance

La radiance sortante ou entrante (radiance est un terme anglais désignant la luminance énergétique) décrit l'intensité lumineuse émise ou reçue par une surface dans une direction donnée divisée par l'aire apparente dans cette direction. La luminance est exprimée en watt par stéradian par mètre carré Wsr⁻¹m⁻².

La radiance est notée L.

<footnote n="1">1. La notion de photon et la dualité onde-particule a été formalisée par le physicien Albert Einstein au début du XXᵉ siècle</footnote>
## 2.1.2 L'Irradiance

L'irradiance (terme anglais désignant l'éclairement) est la puissance lumineuse reçue par une surface par unité d'aire. L'éclairement est exprimé en watt par mètre carré $Wm^{-2}$.
L'irradiance est notée $E$.

## 2.1.3 L'Intensité Lumineuse

L'intensité lumineuse (radiant intensity en anglais) décrit la quantité de lumière émise dans une direction donnée. L'intensité est exprimée en watt par stéradian $Wsr^{-1}$.
L'intensité est notée $I$.

## 2.1.4 Contact avec un Corps

Lors qu'un photon entre en contact avec un corps, le corps accumule l'énergie de la radiation. Puis, le corps va ré-émettre de l'énergie sous forme de photons dans différentes directions. On appelle la première partie l'absorption et la seconde, l'émission (on parle d'émission secondaire car elle est induite par l'absorption).

### L'Absorption

A l'échelle des photons, l'absorption représente le phénomène par lequel l'énergie du photon est prise par une autre particule. Le photon est alors détruit par l'opération.

<figure><img src="image_9.png" /><figcaption>**FIGURE 2.1 – L'Absorption**</figcaption></figure>

La figure 2.1 représente l'étape de l'absorption. Un photon (représenté sous forme ondulatoire) vient toucher une particule. L'onde s'arrête au niveau de la surface de contact et l'énergie du photon est transféré à la particule.

### L'Emission

Après avoir pris l'énergie du photon, la particule peut ré-émettre des photons (émission lumineuse comme le présente la figure 2.2) ou augmenter son agitation (donc la température). On parle généralement de scattering pour décrire l'émission secondaire dans le vocabulaire anglo-saxon. On parle de *n-scattering* pour décrire le parcours d'un photon issu de l'absorption et l'émission secondaire successive de photons par *n* particules.

Pour la plupart des substances, le taux d'absorption varie avec la longueur d'onde de la lumière incidente, ainsi un matériau de couleur noire absorbe toutes les longueurs d'onde (converties en chaleur), alors qu'un matériau de couleur blanche les réfléchit.
<figure><img src="image_10.png" /><figcaption>FIGURE 2.2 – L'Emission</figcaption></figure>

## la Fonction de Phase

La fonction donnant la probabilité de voir de la lumière dans la direction $w_1$ quand une particule est éclairée dans la direction $w_0$ est appelée fonction de phase. Cette fonction caractérise la réponse lumineuse d'une particule à un éclairement.

## Propriétés de la fonction de phase

La fonction de phase est normalisée :

$$ \int_{S^2} \int_{S^2} f(w_0, w_1, \lambda) dw_0 dw_1 = 1 $$

Dans les cas qui nous intéressent, les particules sont sphériques. Cette propriété permet de ne faire dépendre la fonction de phase que de la direction d'émission $w_1$.

## La BRDF

Dans le cas où l'absorbant est une surface, on nomme la fonction de phase par la BRDF (Bidirectional Reflectance Distribution Function en anglais).

La BRDF est définie mathématiquement grâce à la radiance sortante et l'irradiance d'une surface en un point :

$$ BRDF(x, \Omega, \Phi, \lambda) = \frac{dL(x, \Omega)}{dE(x, \Phi)} $$

Où $x$ est le point de la surface étudiée, $\Omega$ est la direction sortante depuis le point $x$, $\Phi$ est la direction entrante vers le point $x$ et $\lambda$ est la longueur d'onde du rayon lumineux entrant.

## Propriétés de la BRDF

La BRDF doit respecter le principe de conservation de l'énergie. Une surface n'étant pas une source primaire ne peut pas renvoyer plus de lumière qu'elle n'en reçoit. On a donc :

$$ \int_{\Omega \in H^2} BRDF(x, \Omega, \Phi, \lambda) d\Omega d\Phi \le 1 $$

Une BRDF bien connue est la BRDF diffuse (aussi appelée BRDF lambertienne). Cette BRDF est constante pour toute direction $\Omega$. On parle généralement d'une surface "lambertienne" (ou diffuse) quand la BRDF de cette dernière est lambertienne.
Nous noterons la BRDF de cette manière : $BRDF(x, \Omega \to \Phi, \lambda)$ pour éviter de faire une confusion entre rayon entrant et sortant.

### 2.1.5 L'Albédo

L'albédo est l'intégrale de la BRDF sur l'hémisphère. Elle varie entre 0 (toute lumière reçue est absorbée) et 1 (toue lumière reçue est réfléchie) grâce à la propriété de la conservation d'énergie. Cette grandeur sans unité permet de caractériser la réflectivité d'une surface. Plus une surface ré-émet de la lumière, plus elle a un albédo élevée.

Par exemple, les nuages ont un albédo proche de 1 car presque toute la lumière atteignant le nuage est retransmise, au contraire de la terre dont l'albédo se situe entre 0.2 et 0.4 (donc une majeure partie des rayons lumineux touchant le sol vont lui transmettre leur énergie sous forme de chaleur).

## 2.2 L'Équation de Transfert de la Lumière

Maintenant que nous avons décrit les différents phénomènes lumineux, nous avons besoin de définir des équations caractérisant le transport de la lumière dans le monde. Calculer le chemin de paquets de photons pour rendre une scène (appelé *photon mapping*) ne serait pas utilisable dans un rendu en temps réel. Nous devons donc faire nos calculs plus globalement.

### 2.2.1 Réflection sur un Support Solide

Le premier à avoir décrit une équation du transport lumineux adaptée à l'informatique graphique mais suffisamment générale pour encore être liée à la physique est James T. Kajiya [Kaj86].

<figure><img src="image_11.png" /><figcaption>FIGURE 2.3 – Réflection sur un support solide</figcaption></figure>

L'équation donnant la radiance sortante du point $x$ de la surface dans la direction $\Omega$ (ou équation de réflection) est :

$$L(x \to \Omega) = \int_{H_x} BRDF(x, \Omega \to \Phi)L(x \leftarrow \Phi)d\Phi + L_e(x \to \Omega) \quad (2.1)$$

Dans cette équation, on fait la distinction entre la *radiance entrante* $L(x \leftarrow \Phi)$ et la *radiance sortante* $L(x \to \Omega)$. Le domaine d'intégration $H_x$ est le dôme des directions partant de $x$ dont le produit scalaire avec la normale en $x$ est positif.
On fait de plus la distinction entre la radiance émise par le matériau $L_e$ et la radiance réfléchie $L_r = \int_{H_x} BRDF(x, \Phi \to \Omega)L(x \leftarrow \Phi)d\Phi$. La plupart des matériaux (bois, marbre, ...) ne sont pas des sources primaires de lumière. On a donc $L_e = 0$.

Il est possible de modifier le terme de radiance entrante $L(x \leftarrow \Phi)$ en un terme de radiance sortante $L(y \to -\Phi)$. On obtient donc une nouvelle forme de l'équation de réflection donnant la radiance sortante comme l'intégrale des radiances sortantes des autres points de la scène (équation 2.2). Il devient donc nécessaire de connaitre l'illumination de la scène pour calculer l'illumination de la scène!

$$L(x \to \Omega) = \int_{H_x} BRDF(x, \Phi \to \Omega)L(y \to -\Phi)d\Phi + L_e(x \to \Omega) \quad (2.2)$$

### 2.2.2 Diffusion dans un Milieu Particulaire

Dans le cas d'un milieu particulaire, on ne décrit pas les calculs sur une surface mais dans un volume. Kajiya et Von Herzen ont détaillé ce calcul en utilisant un volume cylindrique infinitésimal et en considérant le différentiel de l'intensité lumineuse $dI$ dans la direction du cylindre [KMH84].

<figure><img src="image_12.png" /><figcaption>FIGURE 2.4 – Volume cylindrique infinitésimal</figcaption></figure>

Nous voulons obtenir l'intensité lumineuse sortant d'un milieu particulaire. Pour cela, nous allons utiliser le calcul de conservation de l'énergie dans ce volume quand il est soumis à une intensité lumineuse arrivant en $x$ dans la direction $\Omega$ :

$$\frac{1}{\kappa\rho}s\nabla_x I(x \leftarrow \Omega) - I(x \leftarrow \Omega) + \frac{1}{4\pi}\int_S p(\Omega \to \Psi)I(x \to \Psi)d\Psi = 0 \quad (2.3)$$

Le terme $\frac{1}{\kappa\rho}s\nabla_x I(x \leftarrow \Omega)$ est l'atténuation du milieu particulaire.

$I(x \leftarrow \Omega)$ est l'intensité lumineuse entrant dans le volume de particules.

$\frac{1}{4\pi}\int_S p(\Omega \to \Psi)I(x \to \Psi)d\Psi$ est l'intensité lumineuse sortante du volume de particules.

$\kappa$ est le coefficient d'absorption du milieu, $\rho$ est la densité de particules du milieu, $p$ est la fonction de phase du milieu (obtenue à partir de celles des particules). Dans ce calcul, Kajiya utilise l'hypothèse que le milieu n'est pas un diffuseur primaire de lumière (comme une flamme).

L'équation est ici "inversée" par rapport à la réflection sur un support solide car on a l'intensité entrante en fonction de l'intensité sortante et non l'inverse.
### 2.2.3 Le Calcul de Radiosité

On a vu que l'équation de réflection de la lumière 2.1 nécessitait de connaitre la radiance de tous les points de la scène lors du calcul de la radiance d'un point de cette dernière. Cette équation n'est pas solvable telle que présentée. Une méthode de résolution est de faire un calcul de récurrence sur la radiance, en procédant par nombre de réflections de la lumière.

On peut réécrire l'équation 2.1 comme la somme d'une série de récurrence :

$$L(x \to \Omega) = \sum_{i} L^{i}(x \to \Omega) \qquad (2.4)$$

Ou l'on calcule les termes $L^i$ de la façon suivante, en prenant $L^0$ comme l'illumination directe de la scène :

$$L^{n+1}(x \to \Omega) = \int_{H_x} BRDF(x, \Phi \to \Omega)L^n(x \leftarrow \Phi)d\Phi$$

Grâce à cette méthode, on peut résoudre l'équation de réflection 2.2 de manière approchée mais convergente.

## 2.3 Interactions Lumineuses dans les Nuages

Nous avons vu précédemment une vision globale des interactions lumineuses. Nous allons maintenant l'appliquer aux nuages. Nous devons pour cela étudier ce que sont les nuages.

### 2.3.1 Les Nuages

Les nuages sont formés de particules (des amas de molécules) d'eau sphériques (appelées *droplets* en anglais) en suspension dans l'air. Ces gouttelettes d'eau proviennent de la condensation de la vapeur d'eau contenue dans l'air. Ces gouttelettes ont un diamètre de l'ordre de $10^{-5}$m, ce qui entraine une dispersion de la lumière différente de la dispersion atmosphérique qui est due aux molécules de l'atmosphère (l'azote principalement). La dispersion de la lumière peut être expliqué par la théorie de Mie.

### 2.3.2 La Dispersion de Mie

La *dispersion de Mie* apparaît lorsque la lumière entre en collision avec des particules sphériques. Elle fut décrite au début du 19ème siècle par le physicien Gustav Mie [Mie08].

La figure 2.5 permet de visualiser plusieurs éléments importants :

- Le *pic avant* ou *scattering avant* est la zone angulaire autour de 0 degrés ayant la plus forte intensité de réflexion.

- Le *lob avant* qui entraine une diffusion de la lumière dans le nuage.

- Les *pics arrière* ou *scattering arrière* sont composés de trois pics importants responsables d'effets visuels forts (changement notable d'intensité lumineuse à la surface du nuage).
<figure><img src="image_13.png" /><figcaption>**FIGURE 2.5** – Fonction de Phase de Mie pour la lumière rouge (650nm) d'après le logiciel MiePlot [Lav], tracé en coordonées polaire avec une échelle de valeur logarithmique. Le logiciel suppose un simple scattering dans un volume de droplets uniformément réparties dont le rayon varie selon une loie normale autour de 10µm avec une déviation de 1%</figcaption></figure>

## Le Scattering dans les Nuages

La compléxité de la fonction de phase de Mie associée à l'albédo élevé du nuage entraine un nombre de scatters moyens dans le nuage élevé (dépendant des propriétés du nuage). La figure 2.6 présente différents type de scattering dans un nuage.

<figure><img src="image_14.png" /><figcaption>**FIGURE 2.6** – Différents ordres de scattering dans un nuage. En jaune : le 0-scattering. En orange : le 0-scattering. En vert : 2-scattering. En bleu : le n-scattering</figcaption></figure>
## Le Phénomène de Gloire

Le phénomène de Gloire (Glory en anglais) est la conséquence directe du pic arrière central de la fonction de Mie (voir la figure 2.5). Lorsque la caméra (ou l'oeil) est dans la direction anti-solaire, on peut voir ce phénomène. L'ouverture de la Gloire est directement dépendante de la taille des particules d'eau et peut varier de 5 à 20 degrés.

## Le Phénomène d'Arc Blanc

Le phénomène d'Arc Blanc (Fogbow en anglais) est la conséquence directe des deux pics de part et d'autre du pic arrière de la fonction de Mie (voir la figure 2.5). De même qu'une Gloire, un Arc Blanc est variable selon la taille des particules d'eau dans le nuage.

<figure><img src="image_15.png" /><figcaption>**FIGURE 2.7** – Une gloire et un arc blanc vus depuis le Golden Gate (photographie de Mila Zinkova)</figcaption></figure>

La figure 2.7 présente la photographie des deux effets optiques cités précédement. Ces deux effets sont liés par la fonction de Mie et donc apparaissent simultanément.

### Le Phénomène de Silver Lining

Le phénomène de Silver Lining est la conséquence du pic avant de la fonction de phase. Cet effet ce caractérise par une intensité lumineuse plus importante sur le contour du nuage. Le nombre de scattering moyen pour un photon traversant le nuage en bordure étant faible, on peut ressentir l'effet non atténué du pic avant de la fonction de Mie.

### 2.3.3 Transparence

La transparence est la traversée de rayon lumineux à travers un médium sans être absorbé. Dans le cas des nuages, on définira la transparence comme la probabilité q'un rayon lumineux ne soit pas dévié lors de son passage dans le nuage. En considérant
la densité de probabilté des particules du médium comme uniforme, on obtient que la
transparence peut être exprimée comme suit :

$$
\tau(H) = e^{-\kappa H}
$$

Où *H* est la longueur du chemin à travers le nuage et *κ* est le coefficient d'extinction (l'inverse du chemin moyen que parcours un photon avant d'être scatterisé).

On définit aussi la probabilité qu’un rayon lumineux soit dévié a une distance d’exac-
tement *Hm* dans le nuage :

$$
s(H) = \kappa e^{-\kappa}
$$

2.3.4 L'influence des Nuages sur le Sol

Il est évident que les nuages occultent la lumière venant du Soleil et ombrent le sol.
Mais les nuages, à cause de leur albédo très élevée, reflètent aussi la lumière venant du sol.
Cela étant, les nuages stoppent la lumière arrivant sur la Terre, mais permettent aussi aux
régions sous elles de ne pas être trop sombres (en réfléchissant la lumière provenant des
zones éclairées). La figure 2.8 est une représentation de ces deux phénomènes (occlusion
et illumination indirecte du sol).

<figure><img src="image_16.png" /><figcaption>FIGURE 2.8 – Influence des Nuages sur l’Illumination de la Terre. Certains rayons lumineux sont réfléchis par le nuage et ne parviennent jamais au sol. D’autres atteignent la Terre et son réfléchis vers le bas de la couche nuageuse puis réfléchis de nouveaux.</figcaption></figure>
<figure><img src="image_17.png" /></figure>
## 3.1 La Fonction de Phase

Les données physiques et l'analyse microscopique ont permis de décrire la réaction de corps nuageux à l'excitation de rayons lumineux. Borhen et Huffman ont fait l'étude des *scatters* se produisant lorsqu'un rayon lumineux frappe une particule [BH04]. C'est ainsi que nous savons que les particules d'un nuage réagissent avec une fonction de phase particulière : *la fonction de Mie* [Mie08]. A cause des effets variés de la *fonction de Mie*, il n'est pas aisé de l'utiliser. Dans le cadre du rendu de nuages, plusieurs types de fonction de phase sont utilisés :

### 3.1.1 Utilisation d'une autre Fonction de Phase

La première méthode consiste à utiliser une fonction de phase plus simple correspondant à un phénomène physique proche comme la *fonction de Rayleigh*.

La *fonction de Rayleigh* est une fonction de phase pour les aérosols très fins (présent dans l'atmosphère par exemple). Cette dernière ne permet pas de capturer tous les effets de la *la fonction de Mie* (comme l'*arc blanc* par exemple). Harris et Lastra [HL01] l'utilisent tout comme Kajiya et Von Herzen [KVH84]. Grâce à cette approximation, on obtient un rendu de nuages avec de très bonnes performances car la fonction peut être évaluée analytiquement très rapidement. Néanmoins, on peut constater que bon nombres d'effets optiques disparaissent. En effet, dans le cas présent, il n'apparait plus ni d'>*arc blanc*, ni de *gloire* car la *fonction de Rayleigh* n'a pas de pics arrières.

### 3.1.2 Approximation de Mie

Une seconde méthode consiste à utiliser une approximation de *Mie* comme la *fonction de Henyey-Greenstein* qui n'est pas symétrique et rend bien compte du *scattering avant* fort. Cette fonction a été utilisée par Nishita *et al*. [NDN96] mais aussi par Max [Max94].
<figure><img src="image_18.png" /><figcaption>**FIGURE 3.1** – La fonction de phase de Rayleigh tracée par le logiciel MiePlot [Lav] pour la longueur d'onde du rouge (650nm) en coordonnées polaires le rayon arrivant dans la direction $180^\circ \rightarrow 0^\circ$.</figcaption></figure>

Elle gomme totalement les effets de scattering arrière responsables de l'arc blanc et de la gloire. Une autre méthode d'approximation de la fonction de Mie est l'utilisation de gaussiennes comme le font Premoze et al. où ils cherchent à associer des gaussiennes avec des résultats expérimentaux [PAT$^+$04]. On peut voir dans cette méthode une volonté de s'approcher d'un modèle physique complexe avec des éléments mathématiques bien contrôlés. Il est à noter que cette méthode n'a pas été conçue uniquement pour le rendu de nuages et donc que le réalisme dans le cas particulier des nuages n'était pas la priorité de ces chercheurs.

### 3.1.3 Utilisation de Mie

La dernière méthode consiste à utiliser un pré-calcul de la fonction de Mie. Cette fonction dépendant de la taille des particules, de leur densité, etc, il est impossible (car trop coûteux en calcul) d'évaluer la fonction de phase en temps-réel. Cependant, en se restreignant à certains types de nuages (pour que les valeurs de densité et de taille ne varient pas trop), on peut utiliser une fonction de Mie ne prenant en argument que l'angle entre le rayon incident et le rayon réfléchi. Bouthors [BNL06] utilise une telle technique tout comme Riley et al. [REK$^+$04]. L'intérêt est que l'on peut rendre ainsi tous les effets optiques des nuages, mais elle demande des efforts de calculs supplémentaires ainsi que de stocker les valeurs de l'évaluation de Mie car cette fonction ne peut pas être calculée en temps réel. A titre d'exemple, l'algorithme de Borhen et Huffman [BH04] implémenté par Philip Laven [Lav] prend plusieurs minutes de calcul sur un ordinateur actuel pour obtenir la fonction de phase de Mie telle qu'à la figure 2.5.

## 3.2 Le Scattering et l'Importance des Réflexions

Si la fonction de phase permet facilement de savoir la distribution des rayons lumineux après un scatter dans le milieu, il devient très compliqué d'intégrer l'ensemble des scatters dans une simulation temps-réel.
3.2.1 Calcul Analytique d'un Scatter

Blinn a décrit une méthode simple pour calculer un scatter dans le cadre d'une couche
de particules [Bli82]. Dans ce papier, il décrit comment obtenir une solution analytique à
l'équation de scattering dans le cas d'un volume de particules sphériques uniformément
réparties. Il fait l'approximation que la géométrie du volume est dite slab, c'est-à-dire que
pour tout rayon lumineux entrant dans le volume, la surface d'entrée ainsi que la surface
de sortie du volume est plane. La figure 3.2 présente un rayon lumineux entrant dans le
slab et sortant par la face d'entrée après un scatter.

<figure><img src="image_19.png" /><figcaption>FIGURE 3.2 – 1-scatter dans une couche uniforme de particules sphériques</figcaption></figure>

Nous voulons calculer l'intensité émise dans la direction $\theta_v$ à partir de l'intensité d'une source lumineuse directionnelle d'incidence $\theta_l$.

$$
I_v = \int_{0}^{H} s(x_l(l)) f(\theta_l + \theta_v) \tau(x_v(l)) dx_l(l) I_l
$$

$x_l$ est la distance parcourue par le rayon lumineux dans le gaz avant d'être scatterisé,
$x_v$ est la distance parcourue par le rayon lumineux dans le gaz après le scatter.

Dans le cadre d'une surface dite slab (les surfaces sont planes infinies et parallèles),
on peut obtenir un calcul analytique de cette intégrale grâce aux valeurs de $x_l(l) = \frac{l}{\cos(\theta_l)}$
et $x_v(l) = \frac{l}{\cos(\theta_v)}$ :

$$
\begin{align*}
I_v &= \int_0^H s\left(\frac{l}{\cos(\theta_l)}\right) f(\theta_l + \theta_v) \tau\left(x_v\left(\frac{l}{\cos(\theta_v)}\right)\right) \frac{dl}{\cos(\theta_v)} I_l \\
&= f(\theta_l + \theta_v) \frac{\cos(\theta_l)}{\cos(\theta_l) + \cos(\theta_v)} \left(1 - \tau\left(\frac{H}{\cos(\theta_l)} + \frac{H}{\cos(\theta_v)}\right)\right) I_l
\end{align*}
$$

Cette méthode sera reprise par Bouthors [BNL06]. Cependant elle ne permet de calculer qu'un scattering simple, qui si il permet d'obtenir un arc-blanc et une gloire incomplète ne permet pas d'obtenir une illumination réaliste du nuage (qui serait beaucoup trop sombre avec seulement ce terme).

Il est possible de faire un calcul analytique des ordres de scatter plus importants.
Bouthors fait le calcul analytique d'un scattering d'ordre 2 [BNL06] en réutilisant le
modèle de Blinn [Bli82].

$$
I_v = \int_l s(x_l(l)) \int_{\Psi} \int_k f(\Psi) s(x_{1,2}) f(\theta_l + \theta_v - \Psi) \tau(x_{2,exit}) dk d\Psi \frac{dl}{\cos(\theta_v)} I_l
$$
<figure><img src="image_20.png" /><figcaption>FIGURE 3.3 – 2-scatter dans une couche uniforme de particules sphériques</figcaption></figure>

$x_{1,2}$ est la distance parcourue par le rayon lumineux dans le nuage après le premier scatter et avant le second et $x_{2,exit}$ est la distance parcourue par le rayon lumineux après le second scatter et avant de sortir du nuage.

Comme $f(\theta_l + \theta_v - \Psi)$ et $f(\Psi)$ sont indépendantes de $k$, on peut réécrire la double intégrale du second scatter comme :

$$ \tau(x_{1,exit}) \int_{\Psi} f(\Psi) f(\theta_l + \theta_v - \Psi) d\Psi $$

On remarque que ce n'est rien de moins que la convolution de la fonction de phase. En notant $f^2(\theta) = f * f(\theta)$ on obtient :

$$ I_r = f^2(\theta_l + \theta_v) \frac{\cos(\theta_l)}{\cos(\theta_l) + \cos(\theta_v)} \left( 1 - \tau \left( \frac{H}{\cos(\theta_l)} + \frac{H}{\cos(\theta_v)} \right) \right) I_l $$

Ce second scattering lui permet d'obtenir un arc-blanc et une gloire quasi non atteintués. Mais ne permet pas d'obtenir l'illumination complète du nuage car elle est dépend fortement des ordres de scatters plus important.

Dans le cadre des nuages, la surface d'entrée n'est pas plane, mais Bouthors et al. supposent qu'elle l'est localement et utilisent la fonction de phase de Mie pour montrer que les faibles ordres de scatter ne s'écartent pas assez pour invalider cette approximation.

### 3.2.2 Simplification du Nombre de Scatters

Une méthode pour simplifier les calculs de scatters consiste à ne prendre en compte qu'un petit ordre maximum de scatters (entre 2 et 5). Cette méthode est utilisée par Harris et Lastra [HL01]. Mais elle reste trop simpliste car les nuages ont un albédo très élevé (proche de 1) qui rend la diffusion par scatter très importante. Bouthors a montré les ordres élevés de scattering ne sont pas négligeables dans l'illumination du nuage [BNM$^{+08}$].

### 3.2.3 Calcul Analytique de Paquets de Scatter

Une autre méthode consiste à donner une expression analytique de l'intégrale de scatter pour des ensembles d'ordres de scattering. Bouthors utilise cette technique dans Real-time realistic illumination and shading of stratiform clouds [BNL06], ainsi que Riley et al. [REK$^{+04}$].

Bouthors et al. simplifient les calculs des ordres supérieurs (supérieurs à 2) de scatters grâce à un modèle 1D se trouvant sur l'axe du rayon lumineux.
### 3.2.4 Calculs non Analytique de Paquets de Scatter

Dans [BNM<sup>+</sup>08], Bouthors *et al*. calculent les chemins de scatter à partir de surfaces de collection de lumière permettant une estimation du chemin le plus probable à partir d'une étude préliminaire du scattering multiple. Mais ces calculs sont très coûteux et ne sont mis un place que pour un nuage (taux de rafraîchissement de 2 à 10 Hz).

## 3.3 Interaction entre Sol et Nuages

Maintenant que nous avons vu un tour d'horizon des méthodes de rendu de nuages, nous allons nous intéresser aux interactions entre lesdits nuages et le sol. Même si ce phénomène n'est généralement que peu traité, il n'est pas négligeable lors de fortes variations d'albédo du sol [BNL06] [Vid08].

### 3.3.1 Shadow maps

La plupart du temps, les interactions entre les nuages et le sol se bornent à l'action des nuages sur le sol : l'ombre portée. Cette ombre peut être calculée en utilisant une *shadow map* décrite par Williams [Wil78]. Une *shadow map* est une texture donnant la visibilité depuis une source lumineuse.

<figure><img src="image_21.png" /><figcaption>FIGURE 3.4 – La *shadow map* permet de calculer l'illumination directe d'une scène</figcaption></figure>

Cette technique permet de calculer très rapidement l'illumination directe d'une scène. Le coût de calculs étant très faible, cette méthode est très utilisée dans les jeux vidéos comme Far Cry par exemple [Wen06].

### 3.3.2 Rayons de lumière

Les rayons lumineux émergents des trouées nuageuses (l'effet Tyndall) ne sont pas des interactions entre le sol et les nuages, mais ils sont un effet tellement connu qu'il est difficile de ne pas en parler. A notre connaissance, très peu de papiers proposent une telle méthode. Néanmoins, Dobashi *et al.* proposent un modèle qui simule ces rayons [DKY<sup>+</sup>00]. On peut tout de même souligner que leur méthode est trop permissive, faisant apparaître des raies lumineuses dès qu'il existe une séparation entre deux nuages. On sait pourtant que ces rayons apparaissent dans des conditions atmosphériques bien précises (présence de suspensions dans l'atmosphère).
### 3.3.3 Calcul de radiosité

Une méthode pour obtenir une vraie interaction entre le sol et les nuages est de faire un calcul de radiosité entre ces deux éléments cependant celle-ci est très coûteuse. Nous avons vu que les termes de l'équation 2.4 contiennent une intégrale sur tous les éléments de la scène. Pour palier à ce gouffre de performance, Bouthors [BNL06] a utilisé une simplification du système d'intégration par hémicube de Sillion et Puech [SP89]. En utilisant une texture en *mipmap* (texture contenant des version basses résolutions de l'originale), l'intégrale est approchée par la somme des moyennes des radiances sortantes de disques multipliées par la visibilité associée (appelée *form factor*).

$$ \int_x L(x \to j) \approx \sum_i (L_{\text{sortante}}^{i \to j} * \text{form-factor}_{i \to j}) $$

On notera quand même que cette méthode repose sur le fait que les éléments étudiés sont homogènes et n'ont pas de variations trop importantes dans les zones où ils interagissent. Dans le cadre de l'interaction en le bas de la couche nuageuse et le sol, nous voulons pouvoir avoir une variation de l'albédo du sol ainsi qu'une couche nuageuse réaliste (aux contours cassants).

## 3.4 Transferts Radiatifs Temps-réel

La solution que nous proposons dans ce document pour calculer les interactions entre sol et nuage fait parti de la famille des techniques de précalul de transferts radiatifs (*Precomputed Radiance Transfert* ou PRT en anglais). L'idée de ces techniques est de remplacer le calcul d'intégrale de la radiosité par un produit matrice-vecteur :

$$ r = l * T $$

Dans ce calcul, *l* est le vecteur décrivant la lumière de la scène, *T* est la matrice de transfert et *r* est le vecteur décrivant la radiance finale de la scène.

### 3.4.1 Transferts Radiatifs Précalculés

Le pré-calcul des transferts radiatifs est une technique mise au point par Sloan et al. [SKS02]. Elle a pour but de pré-calculer les inter-réflexions lumineuses d'un objet soumis à une lumière distante. Pour cela, on décompose la lumière distante (généralement décrite par une texture comme le montre la figure 3.5) dans une base sphérique bien connue (généralement les harmoniques sphériques).

On obtient un vecteur de coefficients *r* représentant l'illumination de la scène.

Pour chaque point caractéristique de l'objet *p*, on va précalculer le vecteur de transfert $m_p$ qui est la décomposition dans la base de la lumière de la fonction de transfert. La fonction de transfert $M_p(s)$ ($s$ est une direction) est la solution de radiosité de l'équation :

$$ L_p(s') = \frac{\rho_p}{\pi} \int_{s \in H'_s} \bar{L}_p(s) H_p(s) V_p(s) ds \quad (3.1) $$
<figure><img src="image_22.png" /><figcaption>FIGURE 3.5 – Map d'environnement. Cette texture décrit la lumière arrivant au centre de la scène. En assumant que la lumière vient de l'infini, on peut supposer que toutes les zones de l'objet sont éclairées par cette map.</figcaption></figure>

Où $H_p(s) = \max(N_p.s, 0)$ ($N_p$ est la normale en $p$) est le noyau hémisphérique en $p$ (qui correspond à l'influence d'un objet dans la direction $s$ sur $p$) et $V_p(s)$ est la visibilité de l'infini (0 si on rencontre un objet dans cette direction, 1 sinon).

On peut décomposer $\bar{L}_p$ comme la somme de la lumière directe et de la lumière indirecte $\bar{L}_p = L_p^i + L_p^d$. $L_p^i$ est la lumière de la source à l'infini et $L_p^d$ est l'illumination de $p$ par son voisinage sur l'object (qui a été illuminé par la source à l'infini).

Nous allons résoudre l'équation de radiosité en sommant les termes d'illumination

$$M_p^0(s) = \frac{\rho_p}{\pi} H_p(s) V_p(s) y(s) \quad (3.2)$$

Où $y(s)$ est une fonction de la base sur la sphère.

L'équation 3.2 donne le terme d'illumination directe du point $p$ par la la source lumineuse $y$. Pour obtenir l'interaction lumineuse de l'object sur lui même, nous devons effectuer une récurrence sur le calcul d'illumination :

$$M_p^{i+1}(s) = \frac{\rho_p}{\pi} H_p(s) (1 - V_p(s)) M_p^i \quad (3.3)$$

L'équation 3.3 donne la récurrence. En sommant la série des $M_p^i$, on obtient l'illumination globale du point $p$ quand la source lumineuse est $y$. Nous devons effectuer ce calcul pour toutes les fonctions de la base sur la sphère.

Lors du rendu, il suffit de projeter la texture de lumière distance sur notre base sphérique et effectuer l'intégrale dans la base sphérique. Dans le cadre des *harmoniques sphériques* (voir la section 3.4.2), l'intégrale se résume faire le produit scalaire des vecteurs coefficients. Cette propriété à justifiée l'utilisation des harmoniques sphériques pour les *PRT*.

Les *PRT* sont maintenant très utilisées dans les applications temps-réel comme les jeux vidéos. *Halo 2* du studio *Bungie* a été le premier jeu à en utiliser [KSL05].

Cette technique a été étendue à des objects articulés par *Sloan et al.* [SLS05] en décomposant l'objet en région statiques sur lequelles on utilise des harmoniques zonales
<figure><img src="image_23.png" /><figcaption>**FIGURE 3.6** – Comparaison entre le modèle d'illumination globale par environnement map de Ramamoorthi et al. [RH01] (à gauche) et la PRT de Sloan et al. [SKS02] (image de Sloan et al.)</figcaption></figure>

(insensibles au rotations). De plus, d'autres bases que les harmoniques sphériques initialement utilisées ont été appliquées pour permettre de faire des calculs sur toutes les fréquences [XJF+08][TS06]. Par exemple, en utilisant des bases en Ondelettes, des bases de fonctions constantes par morceau, on peut récupérer les interactions lumineuses de hautes fréquences (les ombres cassantes).

L'intérêt de cette méthode est de prendre en compte des effets très complexes comme les *softs-shadows*, les *inter-réflexions* (la figure 3.6 présente certains des ces effets), les *caustiques*, ... L'inconvénient majeur est la nécessité de faire des pré-calculs très longs (plusieurs heures pour une modèle simple). De plus, le PRT ne fonctionne correctement que pour des objets rigides et ne prend pas en compte les interactions lumineuses entre deux objets.

## 3.4.2 Bases Sphériques

Pour la PRT, le choix de la base est déterminant. En effet, on cherche à tout prix à réduire le nombre de coefficients à calculer lors de la projection. Il nous faut donc une base adaptée au type de signal à représenter (plus généralement une carte décrivant l'environnement lumineux appelée *environnement map* en anglais). On utilise des bases sphériques pour décrire ces environnements.

### Harmoniques Sphériques

Les harmoniques sphériques (*spherical harmonics* ou SH en anglais) sont des polynômes sur la sphère qui ont un laplacien nul $\Delta y = 0$. On décrit généralement les harmoniques sphériques comme l'équivalent pour les fonctions sur la sphère des séries de Fourier. Les SH formant une base orthonormée et capturant bien les basses fréquences d'un signal, on les utilise pour simplifier les calculs d'intégral des rayons lumineux.

Une propriété des HS est que l'intégrale du produit de deux fonctions décomposées sur cette base est le produit scalaire des vecteurs coefficients, cela est due à l'orthonormalité de cette base.

L'utilisation des HS pour le graphisme a été proposé par Cabral *et al.* [CMS87]. L'avantage est de transformer l'intégrale en un simple produit scalaire. L'inconvénient des SH est leur fonction de filtre passe-bas. En effet, ces fonctions ne capturent les hautes fréquences que pour des indices élevés (tout comme les séries de Fourier).
<figure><img src="image_24.png" /><figcaption>**FIGURE 3.7** – Représentation de certaines fonctions en harmoniques sphériques. Le bleu représente les valeurs positives, le rouge les négatives. La distance par rapport au centre de la sphère donne la valeur de l'HS</figcaption></figure>

## Fonctions Constantes par Morceaux

Pour pallier au problème de filtrage que crée les harmoniques sphériques, Xu *et al.* [XJF+08] ont proposé d'utiliser des fonctions constantes par morceau sur la sphère comme base. Cette technique permet d'obtenir des résultats très convaincants et diminue fortement l'erreur de troncature des fréquences. Néanmoins il est nécessaire de connaitre l'allure de la fonction que l'on veut décomposer sur la sphère (par exemple la BRDF).

## Ondelettes Sphériques

Les bases en ondelettes sont devenus le nouveau standard de compression d'un signal. Elles permettent la reconstruction de détails locaux grâce au traitement multirésolution, le tout pour un faible nombre de coefficients. Dans la suite, nous ne parlerons que d'ondelettes de seconde génération comme les décrit Daubechies [Dau92] (c'est à dire un couple composé d'une fonction de référence et d'une fonction de subdivision de l'espace). Cependant, les ondelettes sont difficilement applicables à la sphère.

Schröder et Sweldens ont décrit une méthode pour construire une base d'ondelettes sur une sphère sans utiliser la projection d'une base en ondelettes sur une autre surface (cette méthode souffre d'artefacts) [SS95]. D'autres méthodes de subdivision, ou d'autres bases ont étés explorées. On peut citer le travail de Bonneau qui a proposé une base semi-orthogonale (orthogonale pour une subdivision infinie) [Bon99].

Dans son rapport de Master, Lessig décrit la construction d'une base sphérique en ondelettes orthogonales et symétriques [LF08].

### 3.4.3 Lumières multiples

Une façon de se passer du calcul de radiosité est de rajouter des sources lumineuses de faible intensité pour simuler es transferts lumineux "seconds". Chaque source lumineuse va émettre des sources de lumières virtuelles (*Virtual Point Lights* en anglais) aléatoirement. Ces *VPL* sont ensuite atténuées pour que le résultat de l'illumination soit consistant avec la conservation d'énergie. Kinkelin et Liensberger utilisent cette idée [KL08]. La figure 3.8 présente un rendu de scène en utilisant ces sources de lumière virtuels (*virtual point lights*) dans le document d'origine).
<figure><img src="image_25.png" /><figcaption>**FIGURE 3.8** – Rendu d'une scène en utilisant l'algorithme de radiosité instantanée de Kinkelin et Leinsberger [KL08]. Les sources de luminière secondaires sont représentées dans la scène par des sphère jaunes.</figcaption></figure>

Mais si cette méthode permet de donner l'impression d'une radiosité, elle nécessite de créer un nombre important de luminières additionnelles pour simuler des interactions complexes ce qui alourdit énormément les calculs. On nuancera toutefois ce propos en notant que les algorithmes de *lightcuts* [HPB07] et de *row-column sampling* [WFA<sup>+</sup>05] permettent d'obtenir des temps de calcul très raisonnable (de l'ordre de la seconde).

## 3.5 Bilan

Nous avons vu différentes techniques de rendu de nuages. De cette étude ressort que pour le rendu de nuages stratiformes dans ProLand, la première méthode de Bouthors et al. [BNL06] permettra de très bonne performances ainsi qu'un rendu convaincant de nuages de type stratiformes. Cependant, cette méthode n'est pas conçue pour une couche nuageuse ayant la forme d'une surface sphérique. Il nous a donc fallut faire une adaptation et la valider. De plus, le rendu de radiosité n'est pas correct pour les sols non plats. Il nous faut donc trouver d'autres méthode de calcul.

Les techniques de PRT actuelles ne permettent pas de faire le calcul de l'interaction lumineuse entre une couche de nuage dynamique et un sol statique car la lumière est considérée comme distante et l'objet est statique (ou statique par morceau). Nous allons voir comment nous avons créer un algorithme de PRT sur deux slabs permettant de prendre en compte les nuages dynamiques.
"I have always read that the
world, both land and water,
was spherical"

Christophe Colomb - Letters
to the Queen

4

D'un Modèle Plan à une Terre Sphérique

Dans ce document, nous étendons le modèle de nuages de Bouthors [BNL06]. Ce modèle proposait un rendu temps-réel de nuage en s'appuyant sur la simplification suivante : pour une échelle locale, la Terre peut être approché par un plan. Cette approximation permet bon nombre de simplifications de calculs. Notre but étant de faire un rendu temps-réel de nuages à l'échelle de la planète, nous allons modifier ce modèle pour passer d'une Terre plane à une Terre sphérique telle qu'utilisée dans ProLand. Dans ce chapitre, nous allons tout d'abord présenter quelques notions géométriques sur la terre, en faisant nous aussi quelques approximations, puis nous ferons une étude comparative du modèle de Bouthors avec un modèle terrestre sphérique.

4.1 Quelques Définitions Utiles

La Terre peut être approximé par une sphère. Cette simplification quoique erronée (er-
reur de plus de 22 km à l'équateur) permet une simplification des calculs plutôt salutaire,
c'est pourquoi ce modèle sphérique est utilisé dans ProLand.

Dans la suite du texte, on considérera le modèle terrestre suivant : la Terre est une sphère de rayon R. Elle est entourée d'une atmosphère. Dans cette atmosphère on distingue une couche nuageuse dont la base aura toujours la même distance au centre la Terre, notée H₁. La couche nuageuse se finit à la distance H₂ ≥ H₁. La figure 4.1 permet de visualiser cette disposition. On parlera aussi du terme h = H₂ – H₁ hauteur du nuage et du terme H = H₁ – R la distance entre la surface de la Terre et la base de la couche nuageuse.

Les valeurs numériques utilisées dans la suite du texte (cas d'une couche nuageuse stratiforme) sont les suivantes : R = 6,35.10⁶m, H = 2500m. La valeur de H₂ est variable (selon la position sur la couche nuageuse), mais on peut souligner que h est de l'ordre de grandeur de la centaine de mètres.
<figure><img src="image_26.png" /><figcaption>**FIGURE 4.1 – Modèle géométrique de la Terre et de la couche nuageuse**</figcaption></figure>

## 4.2 Chemins Lumineux et Hauteurs Apparentes

Dans cette partie nous allons tout d'abord calculer la hauteur apparente du nuage depuis la caméra, puis celle depuis le Soleil. La hauteur apparente (notée $H_v$) est la longueur de nuage traversée par un rayon lumineux virtuel partant depuis soit la caméra, soit le Soleil.

### 4.2.1 Depuis la Caméra, au-Dessus des Nuages

<figure><img src="image_27.png" /><figcaption>**FIGURE 4.2 – Les différentes configuration lorsque que le rayon provient d'au-dessus des nuages**</figcaption></figure>

On peut distinguer deux scénarios pour un rayon partant de la caméra (en dehors du fait qu'il ne rencontre aucun nuage), comme le montre la figure 4.2 un rayon peut :

- traverser la couche nuageuse et toucher la terre (4.2 figure de gauche).

- traverser la couche nuageuse une ou deux fois et sortir de l'atmosphère (4.2 figures de droite).

<figure><img src="image_28.png" /><figcaption>**FIGURE 4.3 – 1<sup>ère</sup> configuration**</figcaption></figure>
Le premier cas est atteint pour $0 \le \theta \le \sin^{-1}(\frac{R}{H_2})$ où $\theta$ est l'angle entre le rayon entrant dans la couche nuageuse et la normale de la couche nuageuse au point de contact.

On peut cependant trouver une équation plus précise de la longueur du rayon dans le nuage en utilisant le théorème d'Al-Kashi sur le triangle ABC de la figure 4.4, le point A étant le point d'incidence, le point B étant le centre de la Terre et le point C étant le point de sortie du nuage du rayon partant depuis la caméra.

<figure><img src="image_29.png" /><figcaption>FIGURE 4.4 – Utilisation du théorème d'Al-Kashi sur le triangle ABC</figcaption></figure>

D'après ce théorème :

$$H_1^2 = H_v^2 + H_2^2 - 2 \times H_2 \times H_v \times \cos(\theta)$$

On trouve donc une équation polynômiale du second degré en $H_v$ à résoudre :

$$H_v^2 + H_v(-2H_2\cos(\theta)) + (H_2^2 - H_1^2) = 0$$

On trouve que le déterminant $\Delta = H_1^2 - H_2^2\sin^2(\theta)$ est bien positif pour la région angulaire testée. De plus, la solution du problème doit être positive et bien inférieure à $H_2$ (sinon le rayon irait à l'intérieur de la Terre!), cela permet d'éliminer une solution. On a donc :

$$H_v = H_2\cos(\theta) - \sqrt{H_1^2 - H_2^2\sin^2(\theta)}$$

Le second cas est atteint pour $\sin^{-1}(\frac{R}{H_2}) \le \theta \le \frac{\pi}{2}$, dans ce cas nous devons distinguer deux cas :

- le cas où la couche nuageuse est traversée en un seul morceau $\theta \ge \sin^{-1}(\frac{H_1}{H_2})$ (figure 4.5 dessin de gauche) et

- le cas où elle est traversée en deux fois $\theta < \sin^{-1}(\frac{H_1}{H_2})$ avec une interface atmosphérique (figure 4.5 dessin de droite).

Quand on traverse en une seule fois la couche nuageuse, la longueur du rayon est facilement calculable dans le triangle rectangle ABD où D est le milieu de AC.

$$H_v = 2H_2\cos(\theta)$$

Dans l'autre cas, on peut se ramener au même triangle que dans la première configuration. On a donc $H_v = 2H_2\cos(\theta) - \sqrt{H_1^2 - H_2^2\sin^2(\theta)}$.
<figure><img src="image_30.png" /><figcaption>FIGURE 4.5 – Cas d'un rayon lumineux ne touchant pas la Terre</figcaption></figure>

### 4.2.2 Depuis la Caméra, sous la Couche Nuageuse

Cette configuration a le mérite de ne pas avoir de cas "dégénéré" où un rayon issu de la caméra peut traverser deux fois une couche nuageuse. Nous pouvons réutiliser les calculs que nous avons faits en 4.2.1 car la configuration est analogue (voir figure 4.6)

La longueur apparente traversée dans cette configuration est donc :

$$H_v = H_2\cos(\theta) - \sqrt{H_1^2 - H_2^2\sin^2(\theta)}$$

<figure><img src="image_31.png" /><figcaption>FIGURE 4.6 – Caméra sous les nuages</figcaption></figure>

On vérifie que cette formule est valide pour tout $\theta$ provenant d'une position de la caméra sous les nuages, regardant en direction des nuages car cet angle ne dépasse pas $\sin^{-1}(\frac{R}{H_2})$.

### 4.2.3 Depuis le Soleil

Le Soleil n'est pas considéré comme une source ponctuelle, mais comme une source directionnelle (placée à l'infini) confère la figure 4.7.

On voit apparaître une configuration un peu spéciale du soleil dans lequel ce dernier grâce au pic de phase avant (*forward scattering*) peut éclairer le dessous d'un nuage après que ses rayons lumineux aient traversés l'atmosphère (et donc une première couche nuageuse) comme le montre la figure 4.7.

Les rayons lumineux étant calculés de la même façon que les rayons provenant de la caméra, les remarques géométriques données en 4.2.1 restent vraies pour le Soleil.
<figure><img src="image_32.png" /><figcaption>**FIGURE 4.7 –** *Le Soleil comme source directionnelle*</figcaption></figure>

4.2.4 Comparaison avec le Modèle Initial

Nous avons tracé les courbes des hauteurs apparentes pour $\theta \in [0; \pi/2]$ en utilisant l'un et l'autre des modèles, la figure 4.8 présente ce tracé. Si pour des angles faibles $< 0.4$ radians ces deux modèles coïncident, la simplification du sol à un plan induit une "explosion" de la fonction de profondeur (qui est logique dans le cas d'un *slab* infini, mais totalement irréaliste pour la Terre).

<figure><img src="image_33.png" /><figcaption>**FIGURE 4.8 –** *Comparaison des hauteurs apparentes, la caméra au-dessus des nuages*</figcaption></figure>

En ce qui concerne le cas où la caméra est sous le nuage, il y a un ensemble d'angles qui sont impossibles à atteindre pour le modèle sphérique de la Terre ($\theta$ ne peut pas être plus grand que $sin^{-1}(R/H_2)$).

Cependant, dans la zone d'angles atteignables (voir figure 4.9), les deux fonctions sont conformes.

4.2.5 Conclusion

Nous avons vu que le modèle géométrique de Bouthors *et al.* entraine une erreur pour la hauteur apparente du nuage pour des rayons provenant du haut de la couche nuageuse. Dans la suite, nous utiliserons ces longueurs modifiées.
<figure><img src="image_34.png" /><figcaption>FIGURE 4.9 – Comparaison des longueurs apparentes, la caméra en-dessous des nuages</figcaption></figure>

## 4.3 Modification du Modèle Lumineux

Dans cette partie, nous comparons le modèle lumineux du slab avec un modèle lumineux pour une sphère couverte d'une couche nuageuse.

Comme Bouthors, nous allons étudier les différents "ordres" de scattering de façon indépendantes. Cela nous permettra de faire une comparaison valable avec son modèle.

### 4.3.1 Illumination Directe

La transparence décrot exponentiellement en fonction de la hauteur apparente depuis la caméra :

$$ \tau(H_v) = e^{-\kappa H_v} $$

$\kappa$ étant le coefficient d'extinction du nuage.

L'illumination directe du Soleil à travers les nuages aussi appelé forward scattering est l'atténuation exponentielle du rayon lumineux scatterisé $\kappa H_v$ fois pour une déviation inférieure à 5°. $\kappa H_v$ est le nombre moyen de scattering ayant lieu sur une distance $H_v$.

Ces deux termes, transparence et illumination directe du soleil forment l'illumination directe.

Cette illumination ne dépendant pas de la forme du nuage, seulement de la hauteur apparente $H_v$ que nous avons déjà modifiée. Nous n'avons aucun changement à faire sur ce modèle.

### 4.3.2 Single Scattering

Bouthors calcule l'intensité totale sortant du slab dans une direction $\bar{v}$ en intégrant tous les rayons lumineux arrivant dans la direction $\bar{l}$. Le pourcentage de l'intensité transmise après un scatter dans la direction $\bar{v}$ est :
$$I_1 = \int_0^{H_v} s(l) P(\theta_{vl}) \tau(x) dl$$

Dans cette équation $s(l) = \kappa e^{-\kappa H_v}$ est la probabilité de rencontrer une particule après un chemin de distance $l$ dans le nuage. $\tau(x) = e^{-\kappa H_v}$ est la probabilité de parcourir la distance $x$ sans rencontrer de particule. $P(\theta_{vl})$ est la fonction de phase donnant la probabilité de rebondir d'un angle $\theta_{vl}$ après le choc avec la particule. $l$ est la distance d'intégration, elle représente la distance d'un rayon lumineux ayant effectué un rebond (ou scatter) dans la couche nuageuse et qui en ressort. $x$ est la distance parcourue par un rayon lumineux avant d'effectuer un rebond pour une longueur de sortie $l$ ($x$ est en fait une fonction de $l$, $\theta_v$ et $\theta_l$). La figure 4.10 montre ces distances.

<figure><img src="image_35.png" /><figcaption>FIGURE 4.10 – Intégration du single scattering selon $l$</figcaption></figure>

Dans le cas d'une couche nuageuse autour d'une Terre sphérique, l'équation permettant de déterminer la distance $x(l)$ se complique au point de ne plus permettre l'intégration de façon analytique. En effet, $x(l)$ peut se calculer comme suit :

$$ 
\begin{aligned}
C &= A + \vec{A}B + \vec{B}C \\
  &= (0, H_2) + l \times (\sin(\theta_v), -\cos(\theta_v)) + x \times (\sin(\theta_l), \cos(\theta_l))
\end{aligned}
$$

Or $C$ est sur le cercle de rayon $H_2$ :

$$C_x^2 + C_y^2 = H_2^2$$

On peut donc en déduire une équation polynômiale en $x$ :

$$x^2 + 2x(H\sin(\theta_l) - l\cos(\theta_v + \theta_l)) + (l^2 - 2H_2l\cos(\theta_v))$$

La solution d'une telle équation serait un polynôme $x(l)$ de la forme $al^2 + bl + c + \sqrt{dl^2 + el + f}$. La fonction $\tau(x)$ n'a pas de primitive analytique (c'est l'exponentielle de la somme d'un polynôme et de la racine d'un polynôme). Pour utiliser une telle solution, il faudrait ruser en utilisant soit un développement limité, soit un modèle de Terre apportant une solution plus simple (par exemple une parabole).
<figure><img src="image_36.png" /><figcaption>FIGURE 4.11 – Single Scattering pour un slab (graphe de gauche) et pour une Terre sphérique (graphe du milieu) pour une hauteur de nuage de 250m. L'erreur est présentée dans le graphe de droite</figcaption></figure>

Nous avons comparé le pourcentage de lumière obtenue en utilisant notre nouvelle distance *x* avec le pourcentage utilisé par Bouthors. Les deux courbes sont présentées dans les figures 4.11.

Ces calculs montrent que l'approximation par un slab plan est valable. Ceci est expliqué par le fait que l'épaisseur du nuage est très faible par rapport au rayon de la Terre.

### 4.3.3 Double Scattering

L'équation utilisée par Bouthors pour le *scattering* double est très proche de celle du *scattering* simple dans le sens où il se restreint au cas d'un scattering simple avec une fonction de phase égale à la convolution de la fonction de phase avec elle même.

Nos résultats de 4.3.2 étant indépendants de la fonction de phase, nous en tirerons les mêmes conclusions.

### 4.3.4 Multiple Scattering

Le *multiple scattering* est indépendant de la configuration géométrique étudiée étant donné qu'elle se base sur un modèle 1D de *slab*. Ces transmissions/réflexions étant calculées uniquement sur un rayon, seule la longueur du rayon est prise en compte, au contraire du *scattering* de premier et second ordre où la forme de la couche nuageuse entrait dans l'intégrale.

Cette vision du *multiple scattering* est très simpliste. On pourrait même penser qu'elle n'est pas réaliste du tout car elle ne prend en compte qu'une infime partie des rayons. Pourtant elle donne des résultats très convaincants visuellement car il ne faut pas oublier que le scattering lumineux se fait essentiellement dans l'axe de la lumière (voir la fonction de phase de Mie).

### 4.3.5 Conclusion

Nous avons décidé de conserver le modèle d'illumination de Bouthors car après analyse, il reste toujours cohérent dans le cas de la Terre en utilisant les nouvelles hauteur apparentes $H_v$.

On a tout de même pu constater des erreurs visuelles au niveau du coucher de soleil, lorsque la caméra est au-dessus des nuages. Nous n'avons pas pu déterminer si ces erreurs
provenaient du modèle trop simpliste d'illumination du nuage ou du modèle d'illumination d'atmosphère utilisé pour l'illumination diffuse par le ciel [BN08]. Ces deux éléments étant très liés et vu le temps que nous avons passé à valider le modèle, nous avons préféré laisser ceci de côté.
<figure><img src="image_37.png" /></figure>
# 5

Les Transferts Radiatifs entre la Terre et les
Nuages

Dans ce chapitre, nous allons nous intéresser aux interactions lumineuses entre le sol et la couche basse des nuages. La lumière du soleil frappant sur le sol est en partie ré-émise par ce dernier pour aller éclairer les nuages qui à leur tour vont éclairer le sol et ainsi de suite. Ce genre d'interaction lumineuse est appelée *Illumination Globale* car l'illumination d'un point de la scène ne dépend pas uniquement des sources lumineuses présentes, mais aussi des objets présents dans la scène.

Nous allons tout d'abord nous intéresser à la solution proposée par Bouthors [BNL06] *et al.*. Puis nous expliciterons l'algorithme que nous utilisons pour la simulation des rebonds de lumière.

## 5.1 Algorithme de Bouthors *et al.*

Dans cette partie, nous allons nous concentrer sur l'algorithme de radiosité proposé par Bouthors *et al.* [BNL06]. Nous détaillerons en premier lieu cet algorithme. Puis nous étudierons les impacts des hypothèses prises sur l'erreur commise par rapport à un calcul de radiance classique.

### 5.1.1 Description de l'Algorithme

Le calcul de l'illumination se décompose en deux parties :

- Le calcul de l'illumination des nuages par les rayons du Soleil.

- Le calcul de l'illumination du sol et du bas des nuages par radiosité.
## Illumination des Nuages par le Soleil

Pour calculer l'illumination des nuages par le soleil, Bouthors et al. décomposent les ordres de scattering en 4 groupes comme nous l'avons vu au chapitre 4 :

- Le non scattering qui calcule l'illumination par transparence.

- Le scattering simple qui reprend le calcul de Blinn [Bli82].

- Le scattering double qui étend le principe de Blinn pour calculer un scattering double analytiquement.

- Le scattering multiple regroupant tous les termes non calculés.

## Calcul de Radiosité

Pour faire son calcul de radiosité, Bouthors utilise l'algorithme suivant. Pour chaque surface (nuage et sol) il utilise deux textures : une texture T contenant l'intensité lumineuse accumulée et une autre dT contenant l'intensité lumineuse sortante.

Tant que les textures d'accumulation ne convergent pas (généralement 2 itérations pour une albédo au sol faible, jusqu'à 8 quand l'albédo au sol est élevée), on fait le calcul suivant :

Pour chaque point de la texture, on utilise les 3 niveaux de mimap ayant une résolution de H, 2H et 8H par pixel de dT (H est la hauteur du sol au nuage). Chaque pixel du centre des niveaux de mimap est multiplié par une combinaison des form factors de ces niveaux ce qui donne une approximation de l'intégrale du calcul de l'irradiance 2.1. On ajoute ce résultat à la texture d'accumulation.

### 5.1.2 Restriction de l'Intégrale à un Disque

Dans son algorithme, Bouthors et al. intègrent sur un pseudo disque de rayon 8H. On peut comprendre que cette restriction permet de se limiter au 3ème niveau de la mimap, mais ce calcul est-il bien correct ?

Calculons l'erreur commise entre l'intégrale de la visibilité sur une hémisphère et l'intégrale sur la restriction à atan(8) :

$$ \frac{\left| \int_{0}^{2\pi} \int_{0}^{atan(8)} \cos(\theta)\sin(\theta)d\theta d\varphi - \int_{0}^{2\pi} \int_{0}^{\pi/2} \cos(\theta)\sin(\theta)d\theta d\varphi \right|}{\int_{0}^{2\pi} \int_{0}^{\pi/2} \cos(\theta)\sin(\theta)d\theta d\varphi} = 1.5\% $$

Cette approximation est légitime, l'erreur commise est somme toute minime. On peut donc considérer que la zone d'influence d'un point sur un nuage est le disque au sol de rayon 8H. Dans le cas des nuages stratiformes H ≈ 2,5km. On a donc un disque d'influence de rayon 20km. Comme nous l'avons vu dans le chapitre 4, l'approximation de la Terre par un plan est valide pour cette distance, faible par rapport à l'ordre de grandeur de la Terre. On pourra donc aussi se limiter à l'intégration sur un disque de rayon supérieur à 8H pour le calcul de l'irradiance d'un point.

### 5.1.3 Comparaison avec un Sol non Plat

La plus importante restriction de l'algorithme de Bouthors et al. est le fait qu'il considère le sol comme plat. Cette considération permet de ne s'occuper que de faire une
integration sur les textures. Dans la réalité, les déformations du sol (montagnes, val-
lées, gorges, ...) vont créer des discontinuités dans l'émission lumineuse du sol. Les zones
d'ombres ne reflèteront pas autant de lumière que les zones éclairées, les zones pentues ne
diffuseront pas exclusivement vers les nuages, etc.

Les interactions nuage-nuage

Dans le cas d'une couche nuageuse sphérique, il peut y avoir des interactions de la
couche nuageuse avec elle même comme le montre la figure 5.1. Est-il important de les
prendres en compte ? On sait qu'à cause de l'albédo élevé du nuage on ne pouvait pas
négliger les termes de réflexions multiples [BNM+08] [Vid08]. En est-il de même pour des
interactions extérieures ? Pour vérifier cela, nous avons calculé le rapport entre l'irradiance
de la partie visible du nuage et de l'hémisphère totale à partir d'un point du bas du nuage.

<figure><img src="image_38.png" /><figcaption>FIGURE 5.1 – Schéma d’intégration de la radiance du nuage vers le nuage. La zone en jaune correspond à la zone d’influence de la couche nuageuse vers elle même</figcaption></figure>

En prenant un point quelconque sur le bas de la couche nuageuse, le vecteur normal
pointant vers le centre de la Terre, la zone angulaire pointant sur des nuage est $\theta =$
$[asin(\frac{R}{H_1}), \frac{\pi}{2}]$ :

$$
\frac{\int_{\varphi \in [0,2\pi]} \int_{\theta \in [asin(\frac{R}{H_1}), \frac{\pi}{2}]} cos(\theta)sin(\theta)d\theta d\varphi}{\int_{\varphi \in [0,2\pi]} \int_{\theta \in [0, \frac{\pi}{2}]} cos(\theta)sin(\theta)d\theta d\varphi} = 2 \int_{\theta \in [asin(\frac{R}{H_1}), \frac{\pi}{2}]} cos(\theta)sin(\theta)d\theta = 2a
$$

En prenant les valeurs numériques de la partie 4 et pour différentes valeurs de $H_1$ (donc différents $H$), que vaut le terme $a$ ?

<table>
  <tr>
    <td><i>H</i></td>
    <td>500m</td>
    <td>1500m</td>
    <td>2500m</td>
  </tr>
  <tr>
    <td><i>a</i></td>
    <td>0, 8.10<sup>-3</sup></td>
    <td>2, 3.10<sup>-3</sup></td>
    <td>3, 9.10<sup>-3</sup></td>
  </tr>
</table>

On voit bien qu'on peut négliger le terme des interactions nuage-nuage dans nos calculs de radiance, car la lumière provenant du reste des nuages vers un point sur le bas de la couche nuageuse représente moins d'un pour cent du total de la lumière reçue par le nuage.

Les interactions sol-sol

Dans le cas d'un sol non plat, nous devons estimer si l'influence du sol sur lui même
lors d'un calcul de transfert radiatif est importante. Pour cela, nous allons calculer l'erreur
commise par un modèle de sol plat lorsque l'on a une pente d'angle $\alpha$ dans le calcul de la
radiance reçue par un point de la pente.
<figure><img src="image_39.png" /><figcaption>**FIGURE 5.2** – Calcul simple pour un sol non plat</figcaption></figure>

Nous avons modélisé l'expérience suivante : nous comparons le calcul fait sur un sol plat avec celui fait pour une montagne comme le montre la figure 5.2 (la zone d'intégration est le demi-cercle rouge).

Le calcul de l'irradiance reçue par la base de la pente donne dans le cas non plat :

$$
\begin{align*}
E &= \int_{\omega \in H^{2+}} L(\omega) N.\omega d\omega \\
&= \int_{0}^{\pi} \int_{\frac{\pi}{2}-\alpha}^{\frac{\pi}{2}} L_{sol} \cos(\theta) \sin(\theta) d\theta d\varphi + \int_{0}^{\pi} \int_{0}^{\frac{\pi}{2}-\alpha} L_{nuage} \cos(\theta) \sin(\theta) d\theta d\varphi \\
&\quad + \int_{\pi}^{2\pi} \int_{0}^{\frac{\pi}{2}} L_{nuage} \cos(\theta) \sin(\theta) d\theta d\varphi \\
&= \frac{\pi}{2} \left( I_{sol} \sin^2\alpha + L_{nuage} (\cos^2\alpha + 1) \right)
\end{align*}
$$

On a supposé que les intensités émises par la couche nuageuse et le sol sont constantes.
De plus, en supposant que ces deux support reçoivent exactement la même quantité de
lumière, l'équation de l'erreur fait ressortir les albédos respectifs :

$$
\epsilon = \frac{\left| \frac{1}{2} (a_{sol} \sin^2\alpha + a_{nuage} (\cos^2\alpha + 1)) - a_{nuage} \right|}{a_{nuage}}
$$

On a tracé cette fonction d'erreur pour des angles allant de 0 à π/2 en utilisant les valeurs a<sub>sol</sub> = 0.3 et a<sub>nuage</sub> = 0.9. Le résultat est montré à la figure 5.3. On remarque que l'erreur peut atteindre 30% dans des cas extêmes.

On peut faire remarquer que notre expérience ne permet de savoir que l'erreur commise
pour un point du sol, elle ne reflète pas l'erreur répercutée sur un nuage (pour cela voir
plus loin). Néanmoins, on peut remarquer qu'il devient nécessaire d'ajouter les termes
d'interaction sol-sol dans nos calculs quand nous avons de fortes déclivités. Il faut aussi
noter que ce test n'est qu'une version très simplifiée et peut être vu comme une erreur
minimum. En effet, d'autres sources d'erreur vont s'ajouter comme la présence d'autres
montagnes qui bloqueront la visibilité du ciel.

Interactions sol-nuage

Une autre forme d'erreur que génère l'utilisation d'un modèle de Terre plane est la
non prise en compte de l'influence de la pente du sol dans le calcul de l'illumination sur
<figure><img src="image_40.png" /><figcaption>**FIGURE 5.3** – Fonction d'erreur entre un sol non plat et un sol plat, interaction sol-sol</figcaption></figure>

un nuage. La lumière du soleil étant directionnelle, il existe des zones d'ombre dans les reliefs. Ces dernières vont moins contribuer à l'éclairement des nuages.

<figure><img src="image_41.png" /><figcaption>**FIGURE 5.4** – Schéma d'intégration pour un sol non plat, l est la demie-largeur de la montagne, $\beta$ est l'angle d'incidence des rayons du Soleil</figcaption></figure>

Nous avons évalué l'erreur commise par la non prise en compte de ce phénomène dans la simulation décrite par la figure 5.4. Dans cette dernière, nous devons prendre en compte l'incidence du soleil $\beta$ et la pente $\alpha$.

On voit que lorsque le soleil est peu rasant, on obtient une erreur très faible, qui pourrait être négligée. Or dans le cas où l'on a des montagnes importantes et à pente forte (on peut penser à l'Himalaya), les ombres seront vite importantes, même pour des angles d'incidence des rayons du Soleil peu éloignés de l'angle de la montagne.

### 5.1.4 Conclusion sur les Calculs d'Erreur

Nous avons analysé les sources d'erreur que pouvait avoir l'algorithme de radiosité de Bouthors et al.. Il apparaît clairement que nous devons gérer les interactions lumineuses du sol avec lui même, car ces dernières auront un impact sur l'irradiance de la couche nuageuse.
<figure><img src="image_42.png" /><figcaption>**FIGURE 5.5** – Fonctions d'erreur entre un sol non plat et un sol plat dans l'interaction sol-nuage. Le graphe de gauche présente l'erreur obtenue pour une distance entre le sol et le nuage et une largeur de montagne constante, pour plusieurs valeurs de pente. Le graphe de sortie présente l'erreur obtenue pour une distance entre le sol et la couche nuageuse et une pente constante, pour plusieurs valeurs de largeur de montagne.</figcaption></figure>

## 5.2 Notre Algorithme

### 5.2.1 Simplification du Calcul Intégral

Nous allons chercher à simplifier le calcul des transferts radiatifs entre une couche de nuages et le sol. Nous avons choisi de nous inspirer de la méthode de pré-calcul de transferts de radiance de Sloan et al. qui permet d'effectuer un éclairage complexe en temps réel [SKS02].

#### Définition d'un Rebond Lumineux

On appellera par la suite un rebond lumineux entre les nuages et le sol le transfert de la lumière émise par une partie de la couche nuageuse, absorbée puis ré-émise par le sol allant illuminer les nuages. La figure 5.6 illustre cette définition. Nous devons bien surprendre en compte que la lumière est diffusée par les deux couches, ce n'est pas un seul rayon qui est ré-émis, mais une infinité de rayons.

<figure><img src="image_43.png" /><figcaption>**FIGURE 5.6** – Un rebond lumineux, $L_C^n$ est la radiance réémise par la couche nuageuse après n rebonds, $L_G^{n+1}$ est la radiance réémise par le sol après n rebonds et $L_C^{n+\frac{1}{2}}$ est la radiance entrante dans le nuage après n rebonds</figcaption></figure>

On définit le demi-rebond lumineux nuage-sol (respectivement sol-nuage) comme le transfert de la lumière émise par le nuage (resp. le sol) reçue par le sol (resp. le nuage).
Dans la figure 5.6, le demi-rebond nuage-sol est le rayon vert et le demi-rebond sol-nuage est le rayon rouge.

## Transfert entre éléments de bases

Soit {$C_i$}$_{i \in \mathbb{N}}$ une base de fonctions permettant de décrire les fonctions de la couche nuageuse. Soit {$G_i$}$_{i \in \mathbb{N}}$ une base de fonctions permettant de décrire les fonctions du sol.

On définit deux fonctions : $L_{G,i}$ et $L_{C,i}$ qui sont respectivement la radiance reçue par le sol lorsque la radiance émise par les nuages est égale à l'élément $i$ de la base {$C_i$} et la radiance reçue par le nuage lorsque la radiance émise par le sol est égale à l'élément $i$ de la base {$G_i$}. On peut "projeter" ces deux fonctions respectivement sur {$G_i$} et {$C_i$}. Cela donne :

$$L_{G,i}(x) = \sum_{j \in \mathbb{N}} g_{i,j} G_i(x) \qquad (5.1)$$

$$L_{C,i}(x) = \sum_{j \in \mathbb{N}} c_{i,j} C_i(x) \qquad (5.2)$$

Dans ces définitions, $x$ représente une position respectivement sur le sol et sur le bas de la couche nuageuse.

Il nous est donc possible de rendre compte de l'illumination du sol par une fonction des nuages décomposable dans la base des nuages (et réciproquement). Nous allons maintenant voir le calcul de récurrence sur le nombre de rebonds lumineux.

## Principe de superposition

Un propriété fondamentale dans la suite des calculs est le fait que l'intensité lumineuse est cumulative. C'est à dire que lorsque l'on a deux sources lumineuses, l'intensité lumineuse en un point éclairé par l'une et l'autre des sources est la somme de l'intensité de ce point éclairé par une source puis l'autre (dans le cadre d'une lumière non polarisée).

Cela est du au fait que l'énergie absorbée aux photons s'accumule.

## Calcul des rebonds lumineux sans prise en compte des inter-réflexions du sol

Nous voulons calculer notre équation de radiosité comme montré en 3.3.3. Pour cela, nous devons tout d'abord définir un schéma de récurrence sur la radiance sortante du nuage après $n$ rebonds : $L_C^n$.

Supposons que nous avons connaissance de la fonction de radiance du nuage vers le sol après $n$ rebonds. Cette fonction est décomposable sur la base du nuage. On peut donc écrire :

$$L_C^n(x) = \sum_{i \in \mathbb{N}} r_i^n C_i(x)$$

Par application du principe de superposition, la radiance reçue par le sol (égale à la radiance émise par le nuage) est décomposable en fonction des $L_{G,i}$ de la manière suivante :
$$L_G^{n+1}(x) = \sum_{j \in \mathbb{N}} r_j^n L_{G,j}(x)$$

Or les $L_{G,i}$ sont décomposables sur la base des $G_i$. On a donc :

$$L_G^{n+1}(x) = \sum_{j \in \mathbb{N}} r_j^n \sum_{i \in \mathbb{N}} g_{j,i} G_i(x)$$

En faisant une interversion de signe de sommation, on obtient l'égalité suivante :

$$L_G^{n+1}(x) = \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) G_i(x)$$

On voit que l'intensité lumineuse transmise après n rebonds au sol peut être donnée par une multiplication vecteur-matrice entre la matrice de transmission entre les nuages et le sol et le vecteur de radiance après n rebonds : $\vec{b} = r^n \mathbf{G}$

Continuons les calculs pour obtenir l'intensité dans le nuage. Grâce à la décomposition de $L_G$ sur la base $\{G_i\}$, on peut calculer les coefficients de $L_C^{n+\frac{1}{2}}$, la radiance entrante sur le nuage, dans la pseudo-base $\{L_{C,i}\}$ :

$$L_C^{n+\frac{1}{2}}(x) = \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) L_{C,i}(x)$$

De même que pour le sol, nous allons décomposer cette somme de fonction dans la base du nuage :

$$L_C^{n+\frac{1}{2}}(x) = \sum_{k \in \mathbb{N}} \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) c_{i,k} C_k(x)$$

Tout comme pour les coefficients de l'intensité lumineuse transmise au sol, nous pouvons matriciellement calculer les coefficients de transmission du sol au nuage dans la base du nuage $a_i$ par l'opération suivante : $\vec{a} = r^n \mathbf{G} \mathbf{C}$

Nous devons maintenant obtenir les coefficients de la radiance transmise par le nuage au sol. Pour cela, nous devons multiplier la fonction de radiance entrante $L_C^{n+\frac{1}{2}}$ par la fonction de réflectance $R$, cette fonction est décomposable sur la base du nuage :

$$R(x) = \sum_{i \in \mathbb{N}} r_i C_i(x)$$

On obtient donc l'équation :

$$\begin{align*} L_C^{n+1}(x) R(x) &= \sum_{k \in \mathbb{N}} \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) c_{i,k} C_k(x) \times \sum_{l \in \mathbb{N}} r_l C_l(x) \\ &= \sum_{l \in \mathbb{N}} \sum_{k \in \mathbb{N}} \left( \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) c_{i,k} \right) C_k(x) C_l(x) \end{align*}$$

Pour pouvoir appliquer une récurrence sur ce calcul nous devons décomposer cette fonction dans la base $C_i$ :
$$r_m^{n+1} = \sum_{l \in \mathbb{N}} \sum_{k \in \mathbb{N}} \left( \sum_{i \in \mathbb{N}} \left( \sum_{j \in \mathbb{N}} r_j^n g_{j,i} \right) c_{i,k} \right) < C_k C_l, C_m >$$

On peut aussi représenter ce calcul de façon tensorielle, en posant $\Gamma_{i,j}^k = \langle C_i C_j, C_k \rangle$.

$$r^{n+1} = (\vec{r}^n \mathbf{G} \mathbf{C}) \times_1 (\vec{r} \times_2 \mathbf{\Gamma}) \quad (5.3)$$

$\Gamma$ étant un tenseur d'ordre 3, nous devons spécifier dans notre opérateur de multipli- cation sur quelle dimension se fait la multiplication (c'est en fait une contraction au sens tensoriel). Il est possible de contracter cette expression en utilisant les interversions de signe de sommation. Une forme plus compacte de 5.3 est :

$$r^{n+1} = \vec{r}^n \times_1 (\vec{r} \times_2 \Gamma_{GC}) \quad (5.4)$$

Ou $\Gamma_{GC}$ est le produit de la matrice $\mathbf{G} \mathbf{C}$ avec le tenseur $\Gamma$ selon la première dimension du tenseur uniquement. $\Gamma_{GC}^{m, i, k} = \sum_j \mathbf{G} \mathbf{C}_{i,j} \Gamma_{j,k}^m$ est toujours un tenseur d'ordre 3.

## État initial

La lumière traversant le nuage et allant illuminer le sol est la lumière transmise par le ciel multipliée par la transmission du nuage. Or la transmission du nuage est égale à 1 – reflectance. On a donc :

$$L_C^0(x) = L_{sky}(1 - R(x))$$

La lumière initialement au sol peut être calculée grâce à une shadow map résultant de l'occlusion des nuages.

$$L_G^0(x) = L_{sky} SD(x)$$

Où $SD(x)$ est la fonction de shadow map déterminant si le sol est éclairé par le soleil ou non (occlusion des nuages, occlusion du sol sur lui même).

## Ajout de l'influence du sol sur lui même

Pour ajouter dans nos calculs, l'influence du sol sur lui même, nous devons ajouter au terme du sol au rebond $n$ le terme $n-1$ multiplié par une matrice de passage. Cela nous donne $\vec{b}^n = \vec{r}^n \mathbf{G} + \vec{b}^{n-1} \mathbf{S}$.

## Pré-calcul des matrices

Le pré-calcul des matrices $\mathbf{G}$ et $\mathbf{C}$ et du tenseur $\Gamma$ permet l'accélération des calculs. Pour cela, nous avons besoin de réécrire l'intégrale de transfert :

$$L_{G,i}(u) = \iint C_i(x(\theta, \varphi)) \cos(\theta) V_u(\theta, \varphi) d\theta d\varphi$$

Pour obtenir les coefficients $g_{i,j}$, nous devons projeter sur l'élément $G_j$ de la base du sol. Quand la base $\{G_i\}_{i \in \mathbb{N}}$ est orthonormée, on a :
$$
\begin{align*}
g_{i,j} &= \langle L_{G,i}, G_j \rangle \\
&= \int L_{G,i}(u) G_j(u) du \\
&= \int \left( \iint C_i(x(\theta, \varphi)) \cos(\theta) V_u(\theta, \varphi) d\theta d\varphi \right) G_j(u) du
\end{align*}
$$

De même pour les coefficients $c_{i,j}$. Il est à noter que les coefficients du tenseur $\Gamma$ sont indépendant de la forme du sol.

Résultat

Pour calculer les grandeurs qui nous intéressent : $L_G$ et $L_C$, nous devons sommer tous les $L_G^n$ et $L_C^n$. Vidal a montré [Vid08] qu'il était possible de ne calculer que 4 rebonds pour obtenir une convergence de l'illumination de la couche nuageuse et du sol.

$$
L_G = \sum_{i \in \mathbb{N}} L_G^i \qquad (5.5)
$$

$$
L_C = \sum_{i \in \mathbb{N}} L_C^i \qquad (5.6)
$$

5.2.2 Choix des bases {$C_i$}$_{i\in\mathbb{N}}$ et {$G_i$}$_{i\in\mathbb{N}}$

Utilisation de bases sphériques

Nous allons montrer dans cette partie que les calculs précédents sur une sphère com-
plète sont impossibles pour une application en temps réel. Nous allons pour cela utiliser
l'exemple d'une base en harmoniques sphériques.

Le calcul sur une base en harmoniques sphérique a l'effet d'un filtre *passe-bas* sur les variations de la fonction approximée. Ici, nous faisons un calcul simple pour savoir quelle résolution est nécessaire pour avoir une fonction fiable au kilomètre près pour une fonction définie sur la Terre.

Nous allons nous intéresser à la résolution latitudinale, car la résolution longitudinale dépend de cette dernière (et est du même ordre de grandeur). Avec une résolution d'au moins 1 km, nous faisons le calcul suivant :

$$
\frac{\pi R}{l+1} = 1\mathrm{km} \Rightarrow l = \frac{\pi R}{1\mathrm{km}} - 1
$$

On trouve que pour le rayon maximum de la terre, $l = 20036$. Calculons maintenant
le nombre de coefficients nécessaire. On sait que pour chaque $l$ on doit calculer $2l + 1$
coefficients ($m$ varie de $-l$ à $l$). On a donc :

$$
N = l_{max}^2 \Rightarrow N = 401\,441\,296 \text{ coefficients}
$$

Il serait impossible de faire des calculs en temps réel avec une telle base. D'une part la taille requise pour stocker les pré-calculs (qui doivent être exprimés dans la base de la sphère) serait immense, d'autre part l'algorithme ne serait jamais temps réel si nous avons besoin de faire des calculs sur des vecteurs à prêt de $10^8$ éléments.
Restriction des calculs à des disques

Pour rendre notre algorithme adapté au temps réel, nous allons faire nos calcul sur
une subdivision de la Terre.

Nous découpons la sphère en *disques*. Ces disques vont nous permettre d'accélérer les calculs. En effet, on a vu que l'on pouvait restreindre le calcul de l'intégrale de l'irradiance à un disque (section 5.1.2). Nous allons étendre cette approche en restreignant le calcul complet de tous les rebonds lumineux à un disque. On nommera par la suite ce disque le *disque de calcul*.

<figure><img src="image_44.png" /><figcaption>FIGURE 5.7 – Les disques de rendu et de calcul pour le sol et la couche nuageuse</figcaption></figure>

Contrairement à Bouthors *et al.* (qui travaillent pixel par pixel), nous cherchons à calculer le transfert de radiance pour tout un domaine. Nous allons donc chercher à définir un *disque de rendu* inclus dans le *disque de calcul* pour lequel les valeurs des illuminations calculées sont proches des valeurs calculées sur un domaine non fini (disque de calcul de rayon infini). Comme montré dans la figure 5.7, à chaque couple disque de calcul du sol et disque de rendu du sol, on associe le couple de disques pour la couche nuageuse juste au-dessus du premier couple.

Pour pouvoir récupérer l'illumination des couches nuageuse et terrestre, nous devons faire attention à ce que la partition de la Terre en *disques de rendu* (et non *disques de calcul*) décrive bien toute la sphère.

Nous avons décrit les surfaces de calcul et de rendu comme des disques pour des raisons de clarté et pour simplifier les calculs qui vont suivre. Cependant, on laisse le choix à l'implémentation de la forme de la surface (tant que la surface de rendu est incluse dans le disque de rendu et que le disque de calcul est inclus dans la surface de calcul).

Erfour due au disque de calcul

Dans cette partie, nous allons essayer de déterminer des conditions sur la taille du disque de rendu et du disque de calcul. Dans la suite nous noterons $R_r$ le rayon du disque de rendu et $R_c$ le rayon du disque de calcul.

On sait déjà que le disque de rendu est inclus dans le disque de calcul : $R_c \geq R_r$.

On sait de plus, grâce à la section 5.1.2, que le rayon du disque de calcul doit être au minimum de 8H (H étant la distance du sol au nuage) pour le calcul d'un demi-rebond.
Lorsque l'on effectue le calcul d'un demi-rebond suivant, la zone de calcul est étendue par le calcul du premier demi-rebond.

<figure><img src="image_45.png" /><figcaption>**FIGURE 5.8 –** *erreur engendrée par le calcul d'un rebond, R*<sub>c</sub> *est le rayon du disque de calcul et R*<sub>r</sub> *est le rayon du disque de rendu*</figcaption></figure>

En effet, prenons un point sur la bordure du disque de calcul. Lorsque l'on calcule le premier demi-rebond, ce point reçoit de la lumière de la zone hors du disque de calcul représentée par le rayon rouge de la figure 5.8. Quand on effectuera le second demi-rebond, une partie de l'énergie des points de la bordure n'aura pas été calculée et donc on aura une erreur d'illumination.

Nous avons donc étudié quel impact les rayons $R_c$ et $R_r$ avaient sur l'erreur d'illumination commise après $n$ rebonds. Pour cela, nous avons calculé le résultat de l'illumination après plusieurs rebonds entre une couche nuageuse uniforme d'albédo 0.9 et un sol plat d'albédo 0.3 tous deux infinis. Les deux couches sont espacées d'une distance $H = 2.5$km.

Il est facile de calculer le résultat du calcul de radiosité quand on ne restreint pas la zone de calcul car dans ce cas les illuminations du sol et du nuage sont constantes sur les deux plans.

<figure><img src="image_46.png" /><figcaption>**FIGURE 5.9 –** Zoom sur l'erreur pour 4 rebonds et $R_c = 16H$</figcaption></figure>

La figure 5.9 présente un agrandissement d'une coupe de la fonction d'erreur d'illumination à une zone où cette dernière est inférieure à 4%. On voit que pour $R_r < 8H = \frac{1}{2}R_c$ on a une erreur inférieure à 1%.
Nous avons fait des calculs sur plusieurs valeurs de $R_c$ et $H$. Ces résultats sont présentés dans l'annexe A.1. Ces résultats confirment qu'une erreur acceptable est atteinte pour $R_r < \frac{1}{2}R_c$ et que $R_c$ doit être bien supérieur à $8H$ (rayon pour un demi-rebond). Nous avons choisi d'utiliser un rayon de calcul deux fois plus grand que le rayon de rendu ($R_c = 2R_r$) par la suite.

## Utilisation d'une base en ondelettes

Il est clair que les bases en harmoniques hémisphériques entraineront un trop grand nombre de coefficients car nous devons capturer la forme du nuage qui est un signal haute fréquence. Nous avons choisi de nous tourner vers les bases en ondelettes car elles permettent de capturer les hautes fréquences. Les bases en ondelettes ont de plus l'avantage de présenter un meilleur rapport compression/résolution que les bases issues de l'analyse harmonique (Fourier, harmoniques sphériques, ...).

Nous allons ici présenter deux bases très simples qui permettent de mettre en évidence deux caractéristiques sur lesquelles nous avons travaillé. Ces deux bases sont :

**La base objet** suppose que l'on travail sur domaine constant et que l'on veut juste modifier les valeurs à l'intérieur du domaine (par exemple si on veut recolorer un objet dans une image, on n'a pas besoin de connaitre le reste de l'image, on peut utiliser une base décrivant uniquement et entièrement l'objet). On obtient un ensemble qui n'est plus une base de l'espace, mais une base de l'objet.

**La base à angles constants** consiste à décomposer la surface par une projection de la décomposition d'une sphère par angle solide constant. On peut déplacer le centre de la sphère pour obtenir une décomposition différente.

Voici un exemple : supposons que nous étudions la décomposition d'un disque unité par cette méthode. On placera le centre de la sphère à une distance $z$ de la surface de telle façon qu'il se projette sur le centre du disque. La figure 5.10 présente l'idée de la décomposition.

<figure><img src="image_47.png" /><figcaption>FIGURE 5.10 – Décomposition du disque unité par angle solide constant</figcaption></figure>

Cette décomposition permet d'avoir une base ayant une résolution fine au centre du disque et une résolution plus grossière au fur et à mesure qu'on approche de la bordure car l'erreur commise sur les bords du disque est moins importante que celle au milieu.

### 5.2.3 Pseudo-Code

Le programme de calcul de radiosité peut être décomposé en deux parties. La première partie sera lancée avant le processeur de fragments (*pixel shader*) et servira à calculer les
<figure><img src="image_48.png" /><figcaption>FIGURE 5.11 – Exemple de map de nuages et des compressions obtenues pour 15 coefficients</figcaption></figure>

coefficients, sur la base de fonctions, de chaque disque de rendu visible. Ensuite, dans
le pixel shader, on cherchera pour chaque point quel est le jeu de coefficients associés
(ainsi que la position dans la base associée). Une simple évaluation de la fonction donnera
l’illumination du point. Le calcul des réflexions est donc amorti sur les pixels.

Structure de données

Nous devons stocker les données pré-calculées pour chaque disque de calcul. La surface de la Terre est d'environ $s_{\text{Terre}} = 510\ 067\ 420\text{km}^2$. La surface de notre disque de calcul est grâce à l'approximation planaire : $s_{\text{disque}} = \pi(32H)^2$. On obtient donc qu'il faut $\frac{s_{\text{Terre}}}{s_{\text{disque}}} = 80\ 000$ disques à stocker en mémoire (par disque nous entendons jeu de coefficients à stocker). Notons $N_c$ le nombre de disques.

Pour chaque disque, nous devons stocker 4 textures de pré-calcul et une texture qui sera calculée à l'exécution du programme. Les textures de pré-calcul sont : le passage de la base du nuage à celle du sol **G**, le passage de la base du sol à celle du nuage **C**, le passage de la base du sol à celle du sol **S**, les coefficients décrivant le sol dans la base du sol **Sol**. La texture "dynamique" est la texture des coefficients décrivant le nuage dans la base du nuage **Nuage**.

En prenant *n* comme le nombre de coefficients décrivant le nuage dans la base du nuage et *m* comme le nombre de coefficients décrivant le sol dans la base du sol, on obtient que les matrices précédemment décrites ont des tailles :

$$
taille(\mathbf{G}) = n \times m \times N_c
$$

$$
taille(\mathbf{C}) = m \times n \times N_c
$$

$$
taille(\mathbf{S}) = m \times m \times N_c
$$

$$
taille(\mathbf{Sol}) = m \times N_c
$$

$$
taille(\text{Nuage}) = n \times N_c
$$

A cette structure de donnée servant de base à nos calculs, nous devons ajouter des
textures hautement dynamiques (qui seront modifier à chaque rafraichissement écran). Ces
textures servent au calcul temporaire et à stocker les résultats. On notera ces textures *N*
pour la texture du nuage et *S* pour la texture du sol. Ce sont des textures de coefficients
dans les bases respectives. On a donc :
$$ 
\mathit{taille}(N) = n \times N_c
$$

$$ 
\mathit{taille}(S) = m \times N_c
$$

## Algorithme de pré-calcul

Dans le pré-calcul, nous calculons les matrices de passage, ainsi que la matrice du sol. Cette partie n'est pas contrainte à se dérouler en temps réel, on pourra donc ici chercher à faire le moins d'erreur possible en utilisant les formules exactes et non pas les approximations de certaines fonctions.

## Algorithme de pré-rendu

Dans le pré-rendu, nous devons calculer la projection du nuage sur la base de chaque disques de calcul visible. On n'a besoin de le faire qu'une seule fois par texture de nuages, cela peut donc être fait sur le processeur de l'ordinateur et non sur la carte graphique.

## Algorithme de rendu

Lors du rendu (dans le *pixel shader*), nous devons déterminer pour chaque pixel du nuage et du sol quel est son intensité lumineuse. Avant d'entrer sur le *pixel shader*, nous devons calculer les matrices *S* et *N*. On restreint ce calcul aux seuls disques de rendu visibles :

pour chaque disque de rendu visible

dN <- IlluminationNuage()

N <- [0]

dS <- IlluminationSol()

S <- [0]

pour k dans [1; K]

tempS <- G * dN + S * dS

S <- S + tempS

tempN <= C * dS

N <- N + tempN

dS <- tempS

dN <- ProjN * tempN

fin pour

fin pour

La matrice *ProjN* est égale à $\vec{r} \times_2 \Gamma_{GC}$ qui peut être calculée avant la boucle de calcul sur les rebonds lumineux. Cela évite de rajouter cette multiplication matrice-tenseur dans la boucle mais n'empêche pas de devoir stocker les 3 dimensions du tenseur, cela étant interdit par la nécessité de gérer le cas ou le nuage est dynamique.

Le coût de cet algorithme pour chaque disque de calcul est $O(K \times max(m, n)^2)$.
Pour chaque pixel, nous devons déterminer dans quel disque il se trouve, puis nous devons faire l'évaluation des coefficients en ce pixel :

**pour chaque pixel à afficher**
disque $\langle-$ déterminer le disque de rendu associé à pixel
(u, v) $\langle-$ déterminer position pixel

**si pixel dans nuage**
Nuage $\langle-$ récupère coef (N, disque)
retourne IlluminationB(u, v)
+ évaluation de Nuage (u, v, Nuage)

**sinon**
Sol $\langle-$ récupère coef (S, disque)
retourne ShadowMap(u, v)
+ évaluation de Sol (u, v, Sol)

**fin si**
**fin pour**

Le coût des fonctions d'évaluation est $O(n)$ pour évaluation de Nuage et $O(M)$ pour évaluation de Sol dans le pire des cas (évaluation de la somme des coefficients fois la fonction de la base associée).

La fonction ShadowMap donne l'intensité lumineuse du sol calculée à partir de l'intensité du soleil et de l'épaisseur de la couche de nuage (illumination directe avec atténuation).

La fonction IlluminationB donne l'intensité lumineuse calculée par l'algorithme modifié d'Antoine Bouthors sur le bas du nuage (sans le calcul de radiosité).

## Accélération des calculs

Pour tout élément de la base {$C_i$}$_{i \in N}$ on regarde si le coefficient associé est nul ou non. S'il est nul, on n'a pas besoin d'en tenir compte dans les calculs. Soit $C'$ l'ensemble des fonctions $C_i$ n'ayant que des coefficients non nul. Pour chaque élément de $C'$ on connait son influence sur les éléments de la base {$G_i$}$_{i \in N}$. De même pour chaque élément de la base {$G_i$}$_{i \in N}$ on connait son influence sur les éléments de la base {$C_i$}$_{i \in N}$. On connait donc l'influence de {$G_i$}$_{i \in N}$ sur $C'$. Il est possible de ne faire les calculs qu'uniquement sur les éléments à coefficients non nul de la base {$C_i$}$_{i \in N}$. On peut donc accélérer les calculs grâce à cela. Cette méthode est assez proche de ce que fait Snyder [Sny06]. Elle provient de la simple constatation que l'ont peut faire les calculs sur une restriction des bases soit à cause de considération d'orthogonalité, soit à cause de non influence d'éléments nuls.

On peut penser utiliser l'algorithme de Bouthors lorsque le sol est à peu près plat [BNL06]. Cela permet de réduire le nombre de zones à calculer par pré-calcul de radiance.

Il est de plus possible d'accélérer les calculs en ajoutant de la multi-résolution au niveau des disques de calcul. En effet, si on fait varier le rayon des disques de rendu et de calcul en fonction de la distance de la caméra au nuage tout en gardant le même nombre de coefficients par disque, on réduit le nombre de disques à calculer ce qui augmente les performances. Cela nécessite de stocker plusieurs niveaux de pavage de disques, comme une mipmap (qui augmente la taille nécessaire en mémoire de 33%).
5.2.4 Validation de la Solution Proposée

Pour valider notre solution, nous avons comparé les résultats de l’illumination globale avec ceux d’un raytraceur codé pour l’occasion. Nous avons fait ces tests sur des scènes spécifiquement choisies :

Nuage uniforme et sol plat : cette scène est la plus simple de nos tests et permet de d'évaluer si les coefficients de précalcul sont correctement calculés. On peut de plus vérifier l'erreur due à la restriction des calculs à un disque car on connait les valeurs analytiques de radiance sortante pour le nuage et pour le sol.

Nuage uniforme de forme carré et sol plat : cette scène permet de tester que l'influence du nuage sur le sol est correctement gérée. Il a été montré que l'illumination du sol et du nuage peut être approché grâce à des convolutions des textures de radiance initiales [SS98]. On remarque bien ce résultat par l'effet de flou qui apparaît dans les images au niveau des contours du nuage.

Nuage uniforme et montagne au centre : cette scène permet de mesurer l'influence de changement d'albédo et de relief du sol sur le nuage.

Nuage uniforme et sol entouré de montagnes avec variation d'albédo : cette scène permet de mesurer l'influence de changement d'albédo et de relief du sol sur le sol car les montagnes vont illuminer le centre de la scène.

Nous avons vérifié que les résultats obtenus par notre algorithme tendait vers le résultat du simulateur, voir l'annexe A.2.
<figure><img src="image_49.png" /></figure>
"Every cloud has a silver linning."
an English idiom

6

Résultats & Conclusion

**6.1 Résultats**

Dans cet ultime chapitre, nous détaillons les résultats de l’implémentation de notre travail recherche effectué lors de ce master. Nous conclurons en présentant les recherches possibles pour permettre d’aller au delà des limitations que nous avons.

**6.1.1 Illumination des Nuages dans ProLand**

Nous avons implémenté la méthode d’illumination des nuages d’Antoine Bouthors en appliquant les modifications que nous avons développées au cours du chapitre 4. Les effets complexes de lumière des nuages comme des *arcs blancs*, des *gloires*, des *silver linning*, etc.

<figure><img src="image_50.png" /><figcaption>FIGURE 6.1 – Comparaison entre les deux shaders (nouveau en haut)</figcaption></figure>
La figure 6.1 présente une vue du sol lors d'un coucher de soleil. Dans cette figure, nous n'avons pas fait le calcul de radiosité pour présenter uniquement les effets visuels de scattering. On voit bien que l'ancien shader ne rendait pas le *silver linning* car il combinait un calcul de radiance du ciel avec la profondeur du nuage. On peut noter dans les deux images des artefacts qui viennent de la faible résolution de la texture de nuage losqu'elle est vue depuis le sol.

Le *QTAttenuation* permettant les effets visuels des nuages ne ralenti pas le logiciel *ProLand*. Lors de nos tests le taux de rafraîchissement des deux méthodes était le même, cela étant du au fait que le modèle de Bouthors *et al.* contient une majorité de précalculs rendant les calculs très rapides.

<figure><img src="image_51.png" /><figcaption>**FIGURE 6.2 –** *Gloire et Arc Blanc du nouveau shader*</figcaption></figure>

La figure 6.2 présente une vue de l'espace, prise lorsque le soleil est de dos à la caméra. On voit sur cette capture d'écran le phénomène de scattering arrière avec un arc blanc et une gloire.

## 6.1.2 Radiosité entre la Terre et les Nuages dans ProLand

Il nous n'avons pas eu le temps d'implémenter l'algorithme de radiosité sur la base logicielle de ProLand. Néanmoins, nous avons pu valider ce modèle par rapport à un simulateur.

## 6.2 Conclusion

### 6.2.1 Sur le Rendu de Nuages

Nous avons utilisé une méthode permettant de rendre de façon réaliste et temps-réelle les nuages stratiformes. Nous ne nous sommes pas intéressé à d'autres types de nuages qu'il serait bon de pouvoir rendre dans ProLand pour augmenter le réalisme. Dans le cas des cumulus, des stratus et autres nuages, il faudrait surement changer le modèle géométrique et donc la méthode de rendu. Nos nuages sont actuellement stockés sous forme de textures
haute-définition ce qui peut être gênant dans le cadre d'animation de nuages. On pourrait
imaginer utiliser une autre méthode de stockage comme les 'particules' d'Harris et al.
[HL01]

### 6.2.2 Sur le Précalcul de Radiance

Nous avons proposé une méthode de calcul de radiosité temps-réel entre deux surface sphériques grâce à un précalcul de l'interaction entre les décompositions des radiances et à des considérations géométriques sur les zones d'influences entre ces deux objects. Notre méthode nécessite un stockage de nos précalculs, qui force à faire l'équilibre entre stockage, qualité et fréquence de rafraichissement.

Nous n'avons pas eu le temps d'approfondir réellement le choix des bases utilisés. On pourrait penser qu'utiliser des harmoniques hémisphériques comme base pour le disque du sol permettrait d'obtenir des résultats plus lisses, néanmoins nous avons vu qu'il sera dur de les utiliser sur le disque des nuages sans avoir à stocker énormément de coefficients.

L'idée d'utiliser des connaissance géométriques pour simplifier les calculs et précalculs pourrait être étendue à une scène quelconque.
<figure><img src="image_52.png" /></figure>
A.1 Erreur dans la Restriction des Calculs de Rebonds

Les figures A.1, A.2 et A.3 présentent une coupe de la fonction d'erreur lorsque le disque de calcul est respectivement de rayon 8H, 10H et 16H pour une distance des nuages au sol de 500m.

<figure><img src="image_53.png" /><figcaption>FIGURE A.1 – Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 8H (avec H = 500m)</figcaption></figure>

Les figures A.4, A.5 et A.6 présentent une coupe de la fonction d'erreur lorsque le disque de calcul est respectivement de rayon 8H, 10H et 16H pour une distance des nuages au sol de 2500m.

En comparant les graphes ayant les disques de calcul de même rayon, on remarque qu'ils sont similaire. Il n'y a donc pas d'influence de la hauteur séparant les deux plans.

En revanche, on remarque que pour un rayon de 8H, l'erreur minimum est déjà de 1.5%. On approche une erreur de près de 10% pour un rayon de 4H. Pour un rayon de
<figure><img src="image_54.png" /><figcaption>**FIGURE A.2 –** Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 10H (avec $H = 500m$)</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**FIGURE A.3 –** Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 16H (avec $H = 500m$)</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>**FIGURE A.4 –** Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 8H (avec $H = 2500m$)</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>**FIGURE A.5** – Fonction d’erreur lors de l’intégration sur un disque de calcul de rayon 10H (avec $H = 2500m$)</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>**FIGURE A.6** – Fonction d’erreur lors de l’intégration sur un disque de calcul de rayon 16H (avec $H = 2500m$)</figcaption></figure>

16H l’erreur commise sur le demi-disque de calcul est inférieure à 2%.
## A.2 Validation du Précalcul

Dans cette partie, nous présentons les résultats de validation de notre modèle de précalul de radiance. Nous détaillons les 4 scénarios que nous avons choisi pour la validation et pour chaque scénario nous présentons les graphes issus d'un simulateur et de notre modèle côte-à-côte.

Nos calculs, que ce soit pour le simulateur ou pour notre algorithme, se bornent au calcul de radiance sur un domaine non infini. Toutes les scènes sont bornées au disque de calcul, que nous avons simplifié en carré de calcul pour pouvoir utiliser des textures comme entrées de nos programmes.

Les figures que nous présentons par la suite sont des coupes des scènes étudiées.

**Nuage uniforme et sol plat** : cette scène est la plus simple de nos tests et permet de d'évaluer si les coefficients de précalcul sont correctement calculés. On peut de plus vérifier l'erreur due à la restriction des calculs à un disque car on connait les valeurs analytiques de radiance sortante pour le nuage et pour le sol.

Soleil

Nuage

Sol

**FIGURE A.7 –** Premier test de validation. Le nuage est considéré comme constant et recouvrant toute la scène, il est donc la seule source secondaire de lumière de la scène.

La figure A.7 présente la configuration de la scène vue de côté et la figure A.8 présente les graphes de la radiance sortante totale issue du simulateur et de notre programme.

**Nuage uniforme de forme carré et sol plat** : cette scène permet de tester que l'influence du nuage sur le sol est correctement gérée. Il a été montré que l'illumination du sol et du nuage peut être approché grâce à des convolutions des textures de radiance initiales [SS98]. On remarque bien ce résultat par le floutage qui apparaît dans les images au niveau des contours du nuage.

La figure A.9 présente la configuration de la scène vue de côté et la figure A.10 présente les graphes de la radiance sortante totale issue du simulateur et de notre programme.

**Nuage uniforme et montagne au centre** : cette scène permet de mesurer l'influence du relief du sol sur le nuage (il n'y a pas de modification de l'albédo de la montagne). La figure A.11 présente la configuration de la scène vue de côté et la figure A.12 présente les graphes de la radiance sortante totale issue du simulateur et de notre programme.
<figure><img src="image_59.png" /></figure>

**FIGURE A.8** – Résultat du premier test de validation. A gauche : la radiance sortante du sol à la sortie du simulateur. A droite : la radiance sortante du sol lors d'un calcul sur une base objet de 64 composants.

**Nuage uniforme et sol entouré de montagnes** : cette scène permet de mesurer l'influence de changement d'albédo et de relief du sol sur le sol car les montagnes vont illuminer le centre de la scène.

La figure A.13 présente la configuration de la scène vue de côté et la figure A.14 présente les graphes de la radiance sortante totale issue du simulateur et de notre programme.
Soleil

Nuage

Sol

<figure><img src="image_60.png" /><figcaption>FIGURE A.9 – Second test de validation</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>FIGURE A.10 – Résultat du second test de validation. A gauche : la radiance sortante du sol à la sortie du simulateur. A droite : la radiance sortante du sol lors d'un calcul sur une base objet de 64 composants.</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>FIGURE A.11 – Troisième test de validation</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>**FIGURE A.12** – Résultat du troisième test de validation. A gauche : la radiance sortante du sol à la sortie du simulateur. Au milieu : la radiance sortante du sol lors d'un calcul sur une base objet de 64 composants. On peut voir que l'influence de la montagne est plus large pour notre modèle car la fréquence des éléments du sol est inférieur à celle capturée par notre base.</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>FIGURE A.13 – Quatrième test de validation</figcaption></figure>
<figure><img src="image_65.png" /><figcaption>FIGURE A.14 – Résultat du troisième test de validation. A gauche : la radiance sortante du nuage à la sortie du simulateur. A droite : la radiance sortante du nuage lors d'un calcul sur une base objet de 256 composants. On peut voir que l'influence de la montagne est plus large pour notre modèle car la fréquence des éléments du sol est inférieur à celle capturée par notre base. (Le schéma de gauche et de droite ne renvoient pas la même radiance aux extrêmités car le simulateur tronque une partie de la texture d'élévation et d'albédo)</figcaption></figure>
Bibliographie

[BH04] Craig F Bohren and Donald R Huffman. *Absorption and scattering of light by small particles*. Wiley-VCH, cop., 2004.

[Bli82] James F. Blinn. Light reflection functions for simulation of clouds and dusty surfaces. In *SIGGRAPH '82 : Proceedings of the 9th annual conference on Computer graphics and interactive techniques*, pages 21–29, New York, NY, USA, 1982. ACM.

[BN08] Éric Bruneton and Fabrice Neyret. Precomputed atmospheric scattering. *Comput. Graph. Forum*, 27(4):1079–1086, June 2008. Special Issue : Proceedings of the 19th Eurographics Symposium on Rendering 2008.

[BNL06] Antoine Bouthors, Fabrice Neyret, and Sylvain Lefebvre. Real-time realistic illumination and shading of stratiform clouds. 2006.

[BNM+08] Antoine Bouthors, Fabrice Neyret, Nelson Max, Eric Bruneton, and Cyril Crassin. Interactive multiple anisotropic scattering in clouds. In *ACM Symposium on Interactive 3D Graphics and Games, i3D, February, 2008*, pages 173–182, New York, NY, USA, 2008. ACM.

[Bon99] Georges-Pierre Bonneau. Optimal triangular haar bases for spherical data. In *VIS '99 : Proceedings of the conference on Visualization '99*, pages 279–284, Los Alamitos, CA, USA, 1999. IEEE Computer Society Press.

[CMS87] Brian Cabral, Nelson Max, and Rebecca Springmeyer. Bidirectional reflection functions from surface bump maps. *SIGGRAPH Comput. Graph.*, 21(4):273–281, 1987.

[Dau92] Ingrid Daubechies. *Ten lectures on wavelets*. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1992.

[DKY+00] Yoshinori Dobashi, Kazufumi Kaneda, Hideo Yamashita, Tsuyoshi Okita, and Tomoyuki Nishita. A simple, efficient method for realistic animation of clouds. In *SIGGRAPH '00 : Proceedings of the 27th annual conference on Computer graphics and interactive techniques*, pages 19–28, New York, NY, USA, 2000. ACM Press/Addison-Wesley Publishing Co.

[HL01] Mark J Harris and Anselmo Lastra. Real-time cloud rendering. Technical report, Chapel Hill, NC, USA, 2001.

[HPB07] Miloš Hašan, Fabio Pellacini, and Kavita Bala. Matrix row-column sampling for the many-light problem. In *SIGGRAPH '07 : ACM SIGGRAPH 2007 papers*, page 26, New York, NY, USA, 2007. ACM.

[Kaj86] James T. Kajiya. The rendering equation. In *SIGGRAPH '86 : Proceedings of the 13th annual conference on Computer graphics and interactive techniques*, pages 143–150, New York, NY, USA, 1986. ACM.
[KL08] Martin Kinkelin and Christian Liensberger. Instant radiosity, an approach for real-time global illumination. *Forschungsseminar aus Computergraphik und Digitaler Bildverarbeitung*, 2008.

[KSL05] Jan Kautz, Peter-Pike Sloan, and Jaakko Lehtinen. Precomputed radiance transfer: theory and practice. In *SIGGRAPH '05: ACM SIGGRAPH 2005 Courses*, page 1, New York, NY, USA, 2005. ACM.

[KVH84] James T. Kajiya and Brian P Von Herzen. Ray tracing volume densities. *SIGGRAPH Comput. Graph.*, 18(3):165-174, 1984.

[Lav] Philip Laven. Mieplot. http://www.philiplaven.com/mieplot.htm.

[LF08] Christian Lessig and Eugene Fiume. SOHO: Orthogonal and Symmetric Haar Wavelets on the Sphere. *ACM Trans. Graph.*, 27(1), 2008.

[Max94] Nelson L. Max. Efficient Light Propagation for Multiple Anisotropic Volume Scattering. In *Fifth Eurographics Workshop on Rendering*, pages 87-104, Darmstadt, Germany, 1994.

[Mie08] Gustav Mie. Beitrage zur Optik truber Medien, speziell kolloidaler Metallosungen. *Annalen der Physik*, 330(3), 1908.

[NDN96] Tomoyuki Nishita, Yoshinori Dobashi, and Eihachiro Nakamae. Display of clouds taking into account multiple anisotropic scattering and sky light. In *SIGGRAPH '96: Proceedings of the 23rd annual conference on Computer graphics and interactive techniques*, pages 379-386, New York, NY, USA, 1996. ACM.

[PAT+04] Simon Premoze, Michael Ashikhmin, Jerry Tessendorf, Ravi Ramamoorthi, and Nayar Shree. Practical Rendering of Multiple Scattering Effects in Participating Media. In *Eurographics Symposium on Rendering 2004*, pages 363-374, 2004.

[REK+04] Kirk Riley, David S. Ebert, Martin Kraus, Jerry Tessendorf, and Charles D. Hansen. Efficient rendering of atmospheric phenomena. In Alexander Kessler and Henrik Wann Jensen, editors, *Rendering Techniques*, pages 374-386. Eurographics Association, 2004.

[RH01] Ravi Ramamoorthi and Pat Hanrahan. An efficient representation for irradiance environment maps. In *SIGGRAPH '01: Proceedings of the 28th annual conference on Computer graphics and interactive techniques*, pages 497-500, New York, NY, USA, 2001. ACM.

[SKS02] Peter-Pike Sloan, Jan Kautz, and John Snyder. Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments. In *SIGGRAPH '02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques*, pages 527-536, New York, NY, USA, 2002. ACM.

[SLS05] Peter-Pike Sloan, Ben Luna, and John Snyder. Local, deformable precomputed radiance transfer. In *SIGGRAPH '05: ACM SIGGRAPH 2005 Papers*, pages 1216-1224, New York, NY, USA, 2005. ACM.

[Sny06] John Snyder. Code generation and factoring for fast evaluation of low-order spherical harmonic products and squares. Technical report, Microsoft Research, One Microsoft Way Redmond, WA 98052-6399, USA, 2006.
[SP89] F. Sillion and C. Puech. A general two-pass method integrating specular and diffuse reflection. In *SIGGRAPH '89 : Proceedings of the 16th annual conference on Computer graphics and interactive techniques*, pages 335–344, New York, NY, USA, 1989. ACM.

[SS95] Peter Schröder and Wim Sweldens. Spherical wavelets: efficiently representing functions on the sphere. In *SIGGRAPH '95 : Proceedings of the 22nd annual conference on Computer graphics and interactive techniques*, pages 161–172, New York, NY, USA, 1995. ACM.

[SS98] Cyril Soler and François Sillion. Fast calculation of soft shadow textures using convolution. In *Computer Graphics Proceedings, pages 321–332*, Jul 1998. Annual Con.

[TS06] Yu-Ting Tsai and Zen-Chung Shih. All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation. In *SIGGRAPH '06 : ACM SIGGRAPH 2006 Papers*, pages 967–976, New York, NY, USA, 2006. ACM.

[Vid08] Vincent Vidal. Illumination réaliste de terrains en temps-réel. Master's thesis, Université Joseph Fourier, June 2008.

[Wen06] Carsten Wenzel. Real-time atmospheric effects in games. In *SIGGRAPH '06 : ACM SIGGRAPH 2006 Courses*, pages 113–128, New York, NY, USA, 2006. ACM.

[WFA<sup>+</sup>05] Bruce Walter, Sebastian Fernandez, Adam Arbree, Kavita Bala, Michael Doniak, and Donald P. Greenberg. Lightcuts: a scalable approach to illumination. *ACM Trans. Graph.*, 24(3):1098–1107, 2005.

[Wil78] Lance Williams. Casting curved shadows on curved surfaces. In *SIGGRAPH '78 : Proceedings of the 5th annual conference on Computer graphics and interactive techniques*, pages 270–274, New York, NY, USA, 1978. ACM.

[XJF<sup>+</sup>08] Kun Xu, Yun-Tao Jia, Hongbo Fu, Shimin Hu, and Chiew-Lan Tai. Spherical piecewise constant basis functions for all-frequency precomputed radiance transfer. *IEEE Transactions on Visualization and Computer Graphics*, 14(2):454–467, 2008.