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 page 3

Laurent Belcour page 7

"Every thesis begin with a single keystroke" Jorge Cham -Phd Comics 1. La notion de photon et la dualité onde-particule a été formalisée par le physicien Albert Einstein au début du XX me siècle page 13 Chapitre 2. Physique de la Lumière

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é W m -2 .

L'irradiance est notée E.

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 W sr -1 .

L'intensité est notée I.

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ésentation 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 scientifque 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é).

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.

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 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 [START_REF] Bruneton | Precomputed atmospheric scattering[END_REF].

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.

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" [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] 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.

Travail de Master

Nous avons décidé d'étendre le modèle de rendu de nuages stratiformes de Bouthors et al. [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] 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.

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 Laurent Belcour page 11

Chapitre 1. Introduction 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.

page 12

Rendu Temps-Réel de Nuages à l'Echelle Planétaire "Darkness isn't the opposite of light, it is simply its absence" 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.

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 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.

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é W sr -1 m -2 .

La radiance est notée L.

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. onde lumineuse 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.

page 14

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 2.1. La Lumière 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 :

S 2 S 2 f (w 0 , w 1 , λ)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, Ω, Φ, λ) = dL(x, Ω) dE(x, Φ)
Où x est le point de la surface étudiée, Ω est la direction sortante depuis le point x, Φ est la direction entrante vers le point x et λ 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 :

Ω∈H 2 BRDF (x, Ω, Φ, λ)dΩdΦ ≤ 1
Une BRDF bien connue est la BRDF diffuse (aussi appelée BRDF lambertienne). Cette BRDF est constante pour toute direction Ω. On parle généralement d'une surface "lambertienne" (ou diffuse) quand la BRDF de cette dernière est lambertienne.

Laurent Belcour page 15

Chapitre 2. Physique de la Lumière Nous noterons la BRDF de cette manière : BRDF (x, Ω → Φ, λ) pour éviter de faire une confusion entre rayon entrant et sortant.

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).

L'Equation 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.

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 [START_REF] Kajiya | The rendering equation[END_REF]. L'équation donnant la radiance sortante du point x de la surface dans la direction Ω (ou équation de réflection) est :

L(x → Ω) = Hx BRDF (x, Ω → Φ)L(x ← Φ)dΦ + L e (x → Ω) (2.1)
Dans cette équation, on fait la distinction entre la radiance entrante L(x ← Φ) et la radiance sortante L(x → Ω). 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.

page 16

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 2.2. L'Equation de Transfert de la Lumière

On fait de plus la distinction entre la radiance émise par le matériau L e et la radiance réfléchie L r = Hx BRDF (x, Φ → Ω)L(x ← Φ)dΦ. 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 ← Φ) en un terme de radiance sortante L(y → -P hi). 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 → Ω) = Hx BRDF (x, Φ → Ω)L(y → -Φ)dΦ + L e (x → Ω) (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 [START_REF] Kajiya | Ray tracing volume densities[END_REF]. 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 Ω :

-1 κρ s▽ x I(x ← Ω) -I(x ← Ω) + 1 4π S p(Ω → Ψ)I(x → Ψ)dΨ = 0 (2.3)
Le terme 1 κρ s▽ x I(x ← Ω) est l'atténuation du milieu particulaire. I(x ← Ω) est l'intensité lumineuse entrant dans le volume de particules. 1 4π S p(Ω → Ψ)I(x → Ψ)dΨ est l'intensité lumineuse sortante du volume de particules. κ est le coefficient d'absorption du milieu, ρ 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.

Laurent Belcour page 17

Chapitre 2. Physique de la Lumière

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écurence :

L(x → Ω) = i L i (x → Ω) (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 → Ω) = Hx BRDF (x, Φ → Ω)L n (x ← Φ)dΦ
Grâce à cette méthode, on peut résoudre l'équation de réflection 2.2 de manière approchée mais convergente.

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.

Les Nuages

Les nuages sont formés de particules (des amas de molécules) d'eau sphériques (appelées dropplets 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.

La Dispersion de Mie

La dispersion de Mie apparait 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 [START_REF] Mie | Beitrage zur Optik truber Medien, speziell kolloidaler Metallosungen[END_REF].

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).

page 18 Rendu Temps-Réel de Nuages à l'Echelle Planétaire 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. 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éments.

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.

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

Laurent Belcour page 21

État de l'Art

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 [START_REF] Craig | Absorption and scattering of light by small particles[END_REF]. 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 [START_REF] Mie | Beitrage zur Optik truber Medien, speziell kolloidaler Metallosungen[END_REF]. 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 :

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 [START_REF] Mark | Real-time cloud rendering[END_REF] l'utilisent tout comme Kajiya et Von Herzen [START_REF] Kajiya | Ray tracing volume densities[END_REF]. 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.

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. [START_REF] Nishita | Display of clouds taking into account multiple anisotropic scattering and sky light[END_REF] mais aussi par Max [START_REF] Max | Efficient Light Propagation for Multiple Anisotropic Volume Scattering[END_REF].

page 23

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 [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] 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 [START_REF] Craig | Absorption and scattering of light by small particles[END_REF] 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.

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.

page 24

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 3.2. Le Scattering et l'Importance des Réflexions

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 [START_REF] Blinn | Light reflection functions for simulation of clouds and dusty surfaces[END_REF]. 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 dîte 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 3.2 -1-scatter dans une couche uniforme de particules sphériques

Nous voulons calculer l'intensité émise dans la direction θ v à partir de l'intensité d'une source lumineuse directionnelle d'incidence θ l .

I v = H 0 s(x l (l))f (θ l + θ v )τ (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) = l cos(θ l ) et x v (l) = l cos(θv) :

I v = H 0 s( l cos(θ l ) )f (θ l + θ v )τ (x v ( l cos(θ v ) )) dl cos(θ v ) I l = f (θ l + θ v ) cos(θ l ) cos(θ l ) + cos(θ v ) (1 -τ ( H cos(θ l ) + H cos(θ v ) )I l
Cette méthode sera reprise par Bouthors [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF]. 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 [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] en réutilisant le modèle de Blinn [START_REF] Blinn | Light reflection functions for simulation of clouds and dusty surfaces[END_REF].

I v = l s(x l (l)) Ψ k f (Ψ)s(x 1,2 )f (θ l + θ v -Ψ)τ (x 2,exit )dkdΨ dl cos(θ v ) I l
Laurent Belcour page 25

Chapitre 3. État de l'Art 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 (θ l + θ v -Ψ) et f (Ψ)
sont indépendantes de k, on peut réécrire la double intégrale du second scatter comme :

τ (x 1,exit ) Ψ f (Ψ)f (θ l + θ v -Ψ)dΨ
On remarque que ce n'est rien de moins que la convolution de la fonction de phase. En notant f 2 (θ) = f * f (θ) on obtient :

I r = f 2 (θ l + θ v ) cos(θ l ) cos(θ l ) + cos(θ v ) 1 -τ ( H cos(θ l ) + H cos(θ v ) ) I l
Ce second scattering lui permet d'obtenir un arc-blanc et une gloire quasi non atténué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.

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 [START_REF] Mark | Real-time cloud rendering[END_REF]. 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].

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.

Calculs non Analytique de Paquets de Scatter

Dans [BNM + 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 rafraichissement de 2 à 10 Hz).

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 [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] [Vid08].

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 [START_REF] Williams | Casting curved shadows on curved surfaces[END_REF]. Une shadow map est une texture donnant la visibilité depuis une source lumineuse.

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 + 00]. On peut tout de même souligner que leur méthode est trop permissive, faisant apparaitre 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).

Laurent Belcour page 27

Chapitre 3. État de l'Art

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 [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] a utilisé une simplification du système d'intégration par hémicube de Sillion et Puech [START_REF] Sillion | A general two-pass method integrating specular and diffuse reflection[END_REF]. 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 ).

x

L(x → j) ≃ i L i→j sortante * form-factor i→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).

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.

Transferts Radiatifs Précalculés

Le pré-calcul des transferts radiatifs est une technique mise au point par Sloan et al. [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF]. 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 : 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).

L p (s ′ ) = ρ p π s∈H ′ s Lp (s)H p (s)V p (s)ds (
On peut décomposer Lp comme la somme de la lumière directe et de la lumière indirecte Lp = L i p + L d p . L d p est la lumière de la source à l'infini et L d p 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 0 p (s) = ρ p π H p (s)V p (s)y(s) (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 i+1 p (s) = ρ p π H p (s) 1 -V p (s) M i p (3.3)
L'équation 3.3 donne la récurence. En sommant la série des M i p , 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 [START_REF] Kautz | Precomputed radiance transfer : theory and practice[END_REF].

Cette technique a été étendue à des objects articulés par Sloan et al. [START_REF] Sloan | Local, deformable precomputed radiance transfer[END_REF] en décomposant l'objet en région statiques sur lequelles on utlise des harmoniques zonales Laurent Belcour page 29

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 ∆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. [START_REF] Cabral | Bidirectional reflection functions from surface bump maps[END_REF]. 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).

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 [START_REF] Daubechies | Ten lectures on wavelets[END_REF] (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) [START_REF] Schröder | Spherical wavelets : efficiently representing functions on the sphere[END_REF]. 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) [START_REF] Bonneau | Optimal triangular haar bases for spherical data[END_REF].

Dans son rapport de Master, Lessig décrit la construction d'une base sphérique en ondelettes orthogonales et symétriques [START_REF] Lessig | SOHO : Orthogonal and Symmetric Haar Wavelets on the Sphere[END_REF].

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 [START_REF] Kinkelin | Instant radiosity, an approach for real-time global illumination[END_REF]. 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).

Laurent Belcour page 31

Chapitre 3. État de l'Art 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 + 05] permettent d'obtenir des temps de calcul très raisonnable (de l'ordre de la seconde).

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. [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF] 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.

page 32

Rendu Temps-Réel de Nuages à l'Echelle Planétaire "I have always read that the world, both land and water, was spherical"' Christophe Colomb -Letters to the Queen

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

Dans ce document, nous étendons le modèle de nuages de Bouthors [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF]. 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.

Quelques Définitions Utiles

La Terre peut être approximé par une sphère. Cette simplification quoique erronée (erreur 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 1 . La couche nuageuse se finit à la distance H 2 ≥ H 1 . La figure 4.1 permet de visualiser cette disposition. On parlera aussi du terme h = H 2 -H 1 hauteur du nuage et du terme H = H 1 -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 6 m, H = 2500m. La valeur de H 2 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.

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. 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 :

Depuis la Caméra, au-Dessus des Nuages

-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).

H 2 1 = H 2 v + H 2 2 -2 × H 2 × H v × cos(θ)
On trouve donc une équation polynômiale du second degré en H v à résoudre :

H 2 v + H v (-2H 2 cos(θ)) + (H 2 2 -H 2 1 ) = 0
On trouve que le déterminant ∆ = H 2 1 -H 2 2 sin 2 (θ) 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(θ) -H 2 1 -H 2 2 sin 2 (θ)
Le second cas est atteint pour sin -1 ( R H 2 ) ≤ θ ≤ π 2 , dans ce cas nous devons distinguer deux cas :

-le cas où la couche nuageuse est traversée en un seul morceau θ ≥ sin -1 ( H 1 H 2 ) (figure 4.5 dessin de gauche) et -le cas où elle est traversée en deux fois θ < sin -1 ( 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(θ)
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(θ) -H 2 1 -H 2 2 sin 2 (θ).

Laurent

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 : On vérifie que cette formule est valide pour tout θ 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 ( R H 2 ).

H v = H 2 cos(θ) -H 2 1 -H 2 2 sin 2 (θ) θ h

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 apparaitre 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.

Comparaison avec le Modèle Initial

Nous avons tracé les courbes des hauteurs apparentes pour θ ∈ [0; π/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). 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 (θ 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.

Conclusion

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.

Illumination Directe

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

τ (H v ) = e -κHv
κ é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é κH v fois pour une déviation inférieure à 5°. κ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.

Single Scattering

Bouthors calcule l'intensité totale sortant du slab dans une direction v en intégrant tous les rayons lumineux arrivant dans la direction l. Le pourcentage de l'intensité transmise après un scatter dans la direction v est : Dans cette équation s(l) = κe -κHv est la probabilité de rencontrer une particule après un chemin de distance l dans le nuage. τ (x) = e -κHv est la probabilité de parcourir la distance x sans rencontrer de particule. P (θ vl ) est la fonction de phase donnant la probabilité de rebondir d'un angle θ 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, θ v et θ l ). La figure 4.10 montre ces distances. 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 :

C = A + AB + BC = 0, H 2 + l × sin(θ v , -cos θ v ) + x × sin(θ l ), cos(θ l )
Or C est sur le cercle de rayon H 2 :

C 2 x + C 2 y = H 2 2
On peut donc en déduire une équation polynômiale en x :

x 2 + 2x Hsin(θ l ) -lcos(θ v + θ l ) + l 2 -2H 2 lcos(θ v )
La solution d'une telle équation serait un polynôme x(l) de la forme al 2 + bl + c + dl 2 + el + f . La fonction τ (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). 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.

Laurent

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.

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).

Conclusion

Nous avons décidé de conserver le modèle d'illumination de

Algorithme de Bouthors et al.

Dans cette partie, nous allons nous concentrer sur l'algorithme de radiosité proposé par Bouthors et al. [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF]. 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.

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é. -Le non scattering qui calcule l'illumination par transparence.

-Le scattering simple qui reprend le calcul de Blinn [START_REF] Blinn | Light reflection functions for simulation of clouds and dusty surfaces[END_REF].

-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 mipmap 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 mipmap 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.

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 eme niveau de la mipmap, 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 à 2atan(8) :

| 2π 0 atan(8) 0 cos(θ)sin(θ)dθdϕ - 2π π/2 0 cos(θ)sin(θ)dθdϕ| 2π π/2 0 cos(θ)sin(θ)dθdϕ = 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.

Comparaison avec un Sol non Plat

-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

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 θ = [asin( R H 1 ), π 2 ] :

ϕ∈[0,2π] θ∈[asin( R H 1 ), π 2 ] cos(θ)sin(θ)dθdϕ ϕ∈[0,2π] θ∈[0, π 2 ] cos(θ)sin(θ)dθdϕ = 2 θ∈[asin( R H 1 ), π 2 ] cos(θ)sin(θ)dθ = 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 ? H 500m 1500m 2500m a 0, 8.10 -3 2, 3.10 -3 3, 9.10 -3

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 α dans le calcul de la radiance reçue par un point de la pente. Le calcul de l'irradiance reçue par la base de la pente donne dans le cas non plat :

Laurent

E = ω∈H 2+ L(ω)N.ωdω = π 0 π 2 π 2 -α L sol cos(θ)sin(θ)dθdϕ + π 0 π 2 -α 0 L nuage cos(θ)sin(θ)dθdϕ + 2π π π 2 0 L nuage cos(θ)sin(θ)dθdϕ = π 2 I sol sin 2 α + L nuage cos 2 α + 1
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 :

ǫ = 1 2 a sol sin 2 α + a nuage (cos 2 α + 1) -a nuage a nuage
On a tracé cette fonction d'erreur pour des angles allant de 0 à π/2 en utilisant les valeurs a sol = 0.3 et a nuage = 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 page 46

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 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 β et la pente α.

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.

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 apparait 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.

Laurent

Notre Algorithme

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 [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF].

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 sur prendre 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.

L n C L n+ 1 2 C L n+1 C L n+1 G

Transfert entre éléments de bases

Soit {C i } i∈N une base de fonctions permettant de décrire les fonctions de la couche nuageuse. Soit {G i } i∈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) = j∈N g i,j G i (x)
(5.1)

L C,i (x) = j∈N c i,j C i (x) (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 n C . 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 n C (x) = i∈N r n i 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 :

Laurent Belcour page 49

Chapitre 5. Les Transferts Radiatifs entre la Terre et les Nuages

L n+1 G (x) = j∈N r n j L G,j (x)
Or les L G,i sont décomposables sur la base des G i . On a donc :

L n+1 G (x) = j∈N r n j i∈N g j,i G i (x)
En faisant une interversion de signe de sommation, on obtient l'égalité suivante :

L n+1 G (x) = i∈N j∈N r n j g j,i 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 : b = r n 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

n+ 1 2 C
, la radiance entrante sur le nuage, dans la pseudo-base {L C,i } :

L n+ 1 2 C (x) = i∈N j∈N r n j g j,i 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 n+ 1 2 C (x) = k∈N i∈N j∈N r n j g j,i 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 : a = r n GC 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 n+ 1 2 C par la fonction de réflectance R, cette fonction est décomposable sur la base du nuage :

R(x) = i∈N r i C i (x)
On obtient donc l'équation :

L n+1 C (x)R(x) = k∈N i∈N j∈N r n j g j,i c i,k C k (x) × l∈N r l C l (x) = l∈N k∈N i∈N j∈N r n j g j,i c i,k C k (x)C l (x)
Pour pouvoir appliquer une récurrence sur ce calcul nous devons décomposer cette fonction dans la base C i : page 50

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 5.2. Notre Algorithme

r n+1 m = l∈N k∈N i∈N j∈N r n j g j,i c i,k < C k C l , C m >
On peut aussi représenter ce calcul de façon tensorielle, en posant

Γ i,j k =< C i Cj, C k >. r n+1 = r n GC × 1 r × 2 Γ (5.3)
Γ étant un tenseur d'ordre 3, nous devons spécifier dans notre opérateur de multiplication 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 = r n × 1 r × 2 Γ GC (5.4)
Ou Γ GC est le produit de la matrice GC avec le tenseur Γ selon la première dimension du tenseur uniquement. Γ GC m i,k = j GC i,j Γ m j,k 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 -ref lectance. On a donc :

L 0 C (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 0 G (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 b n = r n G + b n-1 S.

Pré-calcul des matrices

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

L G,i (u) = C i x(θ, ϕ) cos(θ)V u (θ, ϕ)dθdϕ
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∈N est orthonormée, on a : Laurent Belcour page 51

Chapitre 5. Les Transferts Radiatifs entre la Terre et les Nuages

g i,j = < L G,i , G j > = L G,i (u)G j (u)du = C i x(θ, ϕ) cos(θ)V u (θ, ϕ)dθdϕ G j (u)du
De même pour les coefficients c i,j . Il est à noter que les coefficients du tenseur Γ sont indépendant de la forme du sol.

Résultat

Pour calculer les grandeurs qui nous intéressent : [START_REF] Vidal | Illumination réaliste de terrains en temps-réel[END_REF] 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 et L C , nous devons sommer tous les L n G et L n C . Vidal a montré [
L G = i∈N L i G (5.5) L C = i∈N L i C (5.6)

Choix des bases {C

i } i∈N et {G i } i∈N

Utilisation de bases sphériques

Nous allons montrer dans cette partie que les calculs précédents sur une sphère complè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 :

πR l + 1 = 1km ⇒ l = πR 1km - 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 2
max ⇒ N = 401 441 296 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.

page 52

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 5.2. Notre Algorithme

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. (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).

Erreur 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 ≥ R r . On sait de plus, grâce à la section 5.1. 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.5km.

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. 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 < 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.

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 Laurent Belcour page 55

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 Terre = 510 067 420km 2 . La surface de notre disque de calcul est grâce à l'approximation planaire : s disque = π(32H) 2 . On obtient donc qu'il faut s Terre s 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'éxé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 N uage 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(G) = n × m × N c taille(C) = m × n × N c taille(S) = m × m × N c taille(Sol) = m × N c taille(N uage) = n × 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 : page 56

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 5.2. Notre Algorithme

taille(N ) = n × N c taille(S) = m × 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. La matrice ProjN est égale à r × 2 Γ 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ée de gèrer le cas ou le nuage est dynamique.

Algorithme de rendu

Le coût de cet algorithme pour chaque disque de calcul est O(K × max(m, n) 2 ). 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).

Laurent

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∈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∈N . De même pour chaque élément de la base {G i } i∈N on connait son influence sur les éléments de la base {C i } i∈N . On connait donc l'influence de {G i } i∈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∈N . On peut donc accélérer les calculs grâce à cela. Cette méthode est assez proche de ce que fait Snyder [START_REF] Snyder | Code generation and factoring for fast evaluation of low-order spherical harmonic products and squares[END_REF]. 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 [START_REF] Bouthors | Real-time realistic illumination and shading of stratiform clouds[END_REF]. 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%).

page 58

Rendu Temps-Réel de Nuages à l'Echelle Planétaire 5.2. Notre Algorithme

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 [START_REF] Soler | Fast calculation of soft shadow textures using convolution[END_REF]. On remarque bien ce resultat par l'effet de flou qui apparait 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.

Laurent Belcour page 59

"Every cloud has a silver lining." an English idiom 6

Résultats & Conclusion

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.

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. 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.

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.

Conclusion

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. 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.

Laurent

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 16H l'erreur commise sur le demi-disque de calcul est inférieure à 2%.

Laurent Belcour page 67

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ène 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. 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 [START_REF] Soler | Fast calculation of soft shadow textures using convolution[END_REF]. On remarque bien ce resultat par le floutage qui apparait 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.

page 68 Rendu Temps-Réel de Nuages à l'Echelle Planétaire A.2. Validation du Précalcul Laurent Belcour page 69

Annexe A. Annexes

  A.1 Erreur dans la Restriction des Calculs de Rebonds . . . . . . . . . . . . . . A.2 Validation du Précalcul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . page 6 Rendu Temps-Réel de Nuages à l'Echelle Planétaire

Figure 1 . 1 -

 11 Figure 1.1 -Quelques résultat 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).

Figure 1 . 2 -

 12 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].

Figure 2

 2 Figure 2.1 -L'Absorption

Figure 2 .

 2 Figure 2.2 -L'Emission

Figure 2 . 3 -

 23 Figure 2.3 -Réflection sur un support solide

Figure 2 . 4 -

 24 Figure 2.4 -Volume cylindrique infinitésimal

Figure 2 . 5 -

 25 Figure 2.5 -Fonction de Phase de Mie pour la lumière rouge (650nm) d'après le logicielMiePlot[Lav], tracé en coordonées polaire avec une échelle de valeur logarithmique.Le logiciel suppose un simple scattering dans un volume de dropplets uniforméments réparties dont le rayon varie selon une loie normale autour de 10µm avec une déviation de 1%

Figure 2 . 6 -

 26 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

Figure 2 . 7 -

 27 Figure 2.7 -Une gloire et un arc blanc vus depuis le Golden Gate (photographie de Mila Zinkova)

Figure 3 . 1 -

 31 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 polaire le rayon arrivant dans la direction 180 → 0.

Figure 3 .

 3 Figure 3.3 -2-scatter dans une couche uniforme de particules sphériques

Figure 3 . 4 -

 34 Figure 3.4 -La shadow map permet de calculer l'illumination directe d'une scène

  Figure 3.5 -Map d'environement. 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.

Figure 3 . 6 -

 36 Figure 3.6 -Comparaison entre le modèle d'illumination globale par environement map de Ramamoorthi et al. [RH01] (à gauche) et la PRT de Sloan et al. [SKS02] (image de Sloan et al.)

  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

Figure 3 . 8 -

 38 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.

Figure 4 . 1 -

 41 Figure 4.1 -Modèle géométrique de la Terre et de la couche nuageuse

Figure 4 . 2 -

 42 Figure 4.2 -Les différentes configuration lorsque que le rayon provient d'au-dessus des nuages

Figure 4 Figure 4 . 4 -

 444 Figure 4.3 -1 ere configuration

Figure 4 . 5 -

 45 Figure 4.5 -Cas d'un rayon lumineux ne touchant pas la Terre

Figure 4 . 6 -

 46 Figure 4.6 -Caméra sous les nuages

  Figure 4.7 -Le Soleil comme source directionnelle

  Hauteurs relatives pour un rayon au-dessus de la surfaceHauteur relative pour un slab Hauteur relative pour la terre

Figure 4 . 8 -

 48 Figure 4.8 -Comparaison des hauteurs apparentes, la caméra au-dessus des nuages

  Figure 4.9 -Comparaison des longueurs apparentes, la caméra en-dessous des nuages

  )P (θ vl )τ (x)dl

Figure 4 . 10 -

 410 Figure 4.10 -Intégration du single scattering selon l

  Scattering pour le cercle et le slab

Figure 4 .

 4 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

  Les Transferts Radiatifs entre la Terre et les Nuages 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 :

  La plus importante restriction de l'algorithmede 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 page 44 Rendu Temps-Réel de Nuages à l'Echelle Planétaire 5.1. Algorithme de Bouthors et al. intégration sur les textures. Dans la réalité, les déformations du sol (montagnes, vallé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-nuageDans 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][START_REF] Vidal | Illumination réaliste de terrains en temps-réel[END_REF]. 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 5. 1

 1 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

Belcour page 45 Chapitre 5 .Figure 5 . 2 -

 45552 Figure 5.2 -Calcul simple pour un sol non plat

5. 1 .

 1 Figure 5.3 -Fonction d'erreur entre un sol non plat et un sol plat, interaction sol-sol

Figure 5 . 4 -

 54 Figure 5.4 -Schéma d'intégration pour un sol non plat, l est la demie-largeur de la montagne, β est l'angle d'incidence des rayons du Soleil

Figure 5 . 5 -

 55 Figure 5.5 -Fonctions d'erreur entre un sol non plat et un sol plat dans l'interaction solnuage. 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.

Figure 5 . 6 -

 56 Figure 5.6 -Un rebond lumineux, L n C est la radiance réémise par la couche nuageuse après n rebonds, L n+1 G est la radiance réémise par le sol après n rebonds et L n+ 1 2

Figure 5 . 7 -

 57 Figure 5.7 -Les disques de rendu et de calcul pour le sol et la couche nuageuse

Figure 5 . 9 -

 59 Figure 5.9 -Zoom sur l'erreur pour 4 rebonds et R c = 16H

Figure 5 . 10 -

 510 Figure 5.10 -Décomposition du disque unité par angle solide constant

Figure 5 . 11 -

 511 Figure 5.11 -Exemple de map de nuages et des compressions obtenues pour 15 coefficients

  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 d i s q u e de rendu v i s i b l e dN <-I l l u m i n a t i o n N u a g e ( ) N <-[ 0 ] dS <-I l l u m i n a t i o n S o l ( ) S <-[ 0 ] pour k dans [ 1 ; K] tempS <-G * dN + S * dS S

Belcour page 57 Chapitre 5 .

 575 Les Transferts Radiatifs entre la Terre et les Nuages 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 p i x e l à a f f i c h e r d i s q u e <-d é t e r m i n e r l e d i s q u e de rendu a s s o c i é à p i x e l ( u , v ) <-d é t e r m i n e r p o s i t i o n p i x e l s i p i x e l dans nuage Nuage <-r é c u p è r e c o e f (N, d i s q u e ) retourne I l l u m i n a t i o n B ( u , v ) + é v a l u a t i o n de Nuage ( u , v , Nuage ) sinon S o l <-r é c u p è r e c o e f ( S , d i s q u e ) retourne ShadowMap ( u , v ) + é v a l u a t i o n de S o l ( u , v S o l ) f i n s i f i n 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).

Figure 6 . 1 -

 61 Figure 6.1 -Comparaison entre les deux shaders (nouveau en haut)

Figure 6 . 2 -

 62 Figure 6.2 -Gloire et Arc Blanc du nouveau shader

6. 2 . 1

 21 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 page 62 Rendu Temps-Réel de Nuages à l'Echelle Planétaire 6.2. Conclusion 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]

Figure A. 1 -

 1 Figure A.1 -Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 8H (avec H = 500m)

Figure A. 2 -Figure A. 3 -FigureFigure A. 5 -Figure A. 6 -

 2356 Figure A.2 -Fonction d'erreur lors de l'intégration sur un disque de calcul de rayon 10H(avec H = 500m)

Figure A. 7 -

 7 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.

Figure A. 8 -

 8 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.

Figure A. 9 -

 9 Figure A.9 -Second test de validation

Figure A. 12 -

 12 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.

Figure A. 14 -

 14 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)

Table des matières

 des Bibliographie 75

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.

  son passage dans le nuage. En considérant

		2.3. Interactions Lumineuses dans les Nuages
	la densité de probabilté des particules du médium comme uniforme, on obtient que la
	transparence peut être exprimée comme suit :
	τ (H) = e -κ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) = κe -κ
	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).
		Lumière provenant du Soleil
	Nuage
		Terre
	Figure 2.8 -
	page 20	Rendu Temps-Réel de Nuages à l'Echelle Planétaire

  Bouthors car après analyse, il reste toujours cohérent dans le cas de la Terre en utilisant les nouvelles hauteurs apparentes H Modification du Modèle Lumineux 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é. 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.

	Les Transferts Radiatifs entre la Terre et les
		Nuages
	Laurent Belcour	page 41

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 page 40

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

  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. Les Transferts Radiatifs entre la Terre et les Nuages Lorsque l'on effectue le calcul d'un demi-rebond suivant, la zone de calcul est étendue par le calcul du premier demi-rebond. Erreur engendrée par le calcul d'un rebond, R c est le rayon du disque de calcul et R r est le rayon du disque de rendu

	Chapitre 5. R c R r
	Figure 5.8 -
	Laurent Belcour	page 53

Remerciements