# Modélisation de la diffusion multiple de la lumière dans un réseau cubique selon une approche discrète

Lionel Simonot, Laurent Signac, Daniel Meneveaux, Mathieu Hébert

► To cite this version:

Lionel Simonot, Laurent Signac, Daniel Meneveaux, Mathieu Hébert. Modélisation de la diffusion multiple de la lumière dans un réseau cubique selon une approche discrète. [Rapport de recherche] université de Poitiers; université Jean Monnet - Saint-Etienne. 2014. hal-01077439

HAL Id: hal-01077439

https://hal.science/hal-01077439v1

Submitted on 24 Oct 2014

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

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

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

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
Rapport de recherche – octobre 2014

Modélisation de la diffusion multiple de la lumière
dans un réseau cubique selon une approche discrète

INSTITUT

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

Lionel Simonot

Institut Pprime - CNRS UPR 3346 - Département Physique et mécanique des matériaux
Université de Poitiers
lionel.simonot@univ-poitiers.fr

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

LIAS

Laurent Signac

Laboratoire d'informatique et d'automatique pour les systèmes
Université de Poitiers
laurent.signac@univ-poitiers.fr

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

Daniel Meneveaux

Laboratoire XLIM - CNRS UMR 7252
Université de Poitiers
daniel.meneveaux@univ-poitiers.fr

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

Mathieu Hébert

Laboratoire Hubert Curien - CNRS UMR 5516 - Groupe ERIS
Université de Lyon - Université Jean Monnet / Saint-Etienne
mathieu.hebert@institutoptique.fr
# Modélisation de la diffusion multiple de la lumière dans un réseau cubique selon une approche discrète

La modélisation de la diffusion multiple de la lumière doit prendre en compte la localisation spatiale des éléments diffusants, et la répartition angulaire des flux donnée par la diffusion simple. Le problème général étant complexe, de nombreuses méthodes de résolution ont été proposées selon des hypothèses simplificatrices. Nous présentons une méthode originale basée sur la discrétisation des éléments diffusants supposés formés un réseau cubique, et sur le découpage en six flux de la lumière diffusée par une seule particule. Cette méthode à six flux avec discrétisation sur les trois dimensions spatiales peut être considérée comme une extension tridimensionnelle de la méthode de Kubelka [Kubelka48] (méthode à deux flux avec discrétisation sur une seule dimension spatiale). Avec un nombre de paramètres réduits pour décrire la diffusion simple (au minimum trois paramètres indépendants : probabilités d'absorption, de diffusion vers l'avant et vers l'arrière), ce nouveau modèle permet d'accéder à la répartition de la lumière diffusée latéralement dans le matériau.

## I. Introduction

Les matériaux translucides naturels ou synthétiques (tissus biologiques, peau humaine, émail des dents, végétaux, neige, cire, marbre, papier...) sont très fréquents et la simulation réaliste de la diffusion de la lumière par ces matériaux est une problématique importante dans de nombreux domaines (synthèse d'image, télédétection, prédiction couleur en impression, mesure et modélisation de l'apparence visuelle pour la cosmétique ou dans un but marketing...).

Ces matériaux sont constitués de particules ou centres diffusants. Après interaction avec une de ces particules, la lumière peut être en partie absorbée et en partie diffusée (diffusion simple), puis à nouveau diffusée par une particule voisine (diffusion multiple). L'aspect translucide des matériaux s'explique par la diffusion multiple dans l'épaisseur du matériau mais surtout latéralement. D'un point de vue métrologique, pour décrire la tache de diffusion, il faut déterminer la PSF (Point Spread Function) du matériau, c'est-à-dire la répartition latérale du flux diffusé pour un rayon lumineux incident impactant la surface en un point. Il est parfois plus pratique de manipuler la fonction de transfert de modulation, transformée de Fourier de la PSF. Cet outil est en particulier utilisé en impression [Rogers97] avec des PSF pour les papiers déterminés la plupart du temps de manière empirique sous la forme d'une Gaussienne [Yule67] ou d'une fonction exponentielle [Wakeshima68]. Une description plus générale est donnée par la BSSRDF (bidirectional surface scattering reflectance distribution function [Nicodemus77]), rapport entre la luminance en réflexion du matériau et l'éclairement reçu par celui. La BSSRDF dépend de la longueur d'onde, des directions d'incidence et d'observation (chacune définie par les deux angles, zénithal et azimutal) et des positions du point d'impact et du point de sortie de la lumière. En l'absence de diffusion latérale, point d'impact et point de sortie sont confondus, la PSF est impulsionnelle et la BSSRDF se résume à la BRDF (bidirectional reflectance distribution function).

Pour une modélisation réaliste des matériaux, les interactions aux interfaces doivent bien sûr être prises en compte comme conditions aux limites. Mais cet aspect n'est pas abordé ici : nous considérons exclusivement les phénomènes d'absorption et de diffusion dans le volume
du matériau. En outre, nous supposons qu'il n'y a pas de relation de phase entre deux rayons
diffusés consécutifs (régime de diffusion multiple incohérente). Nous supposons enfin que la
lumière incidente est non polarisée et négligeons tout effet de polarisation par diffusion.

La diffusion simple décrit à la fois la proportion de lumière absorbée et diffusée par une particule unique mais aussi, par le biais de la fonction de phase, la répartition angulaire du flux diffusé. Dans l'approche électromagnétique de Mie [Bohren83, Born99], la diffusion simple est traduite comme un phénomène de diffraction d'une onde plane par une particule sphérique. Mais les matériaux réels présentent souvent une forte dispersion en taille et en forme des particules diffusantes. C'est pourquoi, plutôt que la théorie de Mie, des fonctions de phase empiriques et paramétriques (e.g. Henyey-Greenstein [Henyey41]) sont davantage utilisées afin de décrire le comportement moyen en diffusion simple de ces particules.

Modéliser la diffusion multiple revient à faire un bilan des flux absorbés et diffusés en y introduisant des caractéristiques de la diffusion simple. Il y a deux grands moyens pour y parvenir : soit simuler ces événements de manière probabiliste par des techniques de type Monte-Carlo, soit résoudre l'équation de transfert radiatif (ETR) [Chandrasekhar60]. La première est souvent extrêmement couteuse en temps de calcul mais parfois la seule solution envisageable notamment dans le cas de géométries complexes. Pour des géométries simples, en particulier pour des matériaux homogènes arrangés en une superposition de couches, la résolution de l'ETR est privilégiée. Cette dernière a une expression intégro-différentielle qui prend en compte la position spatiale du flux mais aussi les directions des flux incidents et diffusés. L'ETR étant complexe à résoudre dans le cas général, il existe un très grand nombre de techniques de résolution et d'approximations.

Une première approche est de ne pas tenir compte de la diffusion latérale et donc de ramener
le problème à une seule dimension spatiale selon la profondeur (en z).

<figure><img src="image_6.png" /><figcaption>**Figure 1 :** Quelques grandes familles de résolution de l'ETR sans prise en compte de la diffusion latérale (schéma non exhaustif)</figcaption></figure>
Sous cette forme, une première technique de résolution consiste à discrétiser les directions des flux en N angles solides et à intégrer en z. Ce modèle à N-flux a été initialement proposé par Mudgett [Mudgett71] dans le cas d'une diffusion isotrope en azimut, et généralisé par Stamnes [Stamnes88] sous le nom de méthode aux ordonnées discrètes. Pour N=2, on considère seulement deux flux hémisphériques, l'un montant, l'autre descendant. La méthode à deux flux correspond donc au modèle de Kubelka-Munk établi dès 1931 [Kubelka31] sous la forme d'un système de deux équations différentielles couplées, et dont la résolution offre une relation analytique très simple. Une autre technique de résolution, « the adding-doubling method » [van de Hulst81, Prahl95], consiste au contraire en une discrétisation en épaisseur de la couche et en une intégration angulaire des flux. Enfin, il est possible de discrétiser à la fois en z et en direction. Proposée par Kubelka [Kubelka48], cette méthode consiste à modéliser le matériau en une superposition de couches diffusantes. La diffusion multiple y est alors décrite en termes de réflexions et transmissions multiples par ces couches. Ce modèle à deux flux, discret en z, peut être formalisé en utilisant des approches de type chaines de Markov [Hébert07] ou marche aléatoire [Simon03]. La correspondance entre l'approche discrète de Kubelka et l'approche continue de Kubelka-Munk a été établie et discutée [Vöge07, Hébert08]. Enfin, des extensions du modèle à deux flux discret ont été proposées afin de modéliser des milieux particulaires en associant la reflectance et transmittance d'une couche aux fractions de diffusion vers le haut et vers le bas d'une particule, déduites de la diffusion simple [Melamed63, Mandelis90, Shkuratov99, Simonot08]. Ces modèles sont basés sur une forte approximation qui consiste à simplifier le problème tridimensionnel de la diffusion par une collection de particules en un problème unidimensionnel de couches diffusantes comme illustré sur la figure 2.

<figure><img src="image_7.png" /><figcaption>**Figure 2** : Schémas de diffusion multiple de la lumière (a) dans un milieu particulaire (avec *d*, le diamètre moyen des particules diffusantes) (b) à travers une superposition de couches (avec *d*, l'épaisseur moyenne des couches diffusantes) d'après [Shkuratov99].</figcaption></figure>

L'autre approche de résolution de l'ETR consiste à tenir compte de la diffusion latérale et doit donc prendre en compte l'aspect tridimensionnel du milieu. Sous cette forme, l'ETR présente cinq variables indépendantes (les variables spatiales x, y et z et les deux variables angulaires zénithale et azimutale).
<figure><img src="image_8.png" /><figcaption>**Figure 3** : Quelques grandes familles de résolution de l'ETR avec prise en compte de la diffusion latérale (schéma non exhaustif)</figcaption></figure>

Afin de réduire le nombre de variables indépendantes, l'approximation de diffusion [Ishimaru78] est l'approche la plus courante. Celle-ci s'applique plutôt à des matériaux fortement diffusants pour lesquels la distribution de la lumière a tendance à devenir isotrope quand le nombre d'événements de diffusion à prendre en compte augmente. Ceci reste valable même si la source de lumière incidente et la diffusion simple sont très anisotropes. Des hypothèses simplificatrices sur l'expression de la luminance à l'intérieur du matériau permettent alors de simplifier l'ETR. La méthode du dipôle [Pharell92] consiste à positionner deux points sources de part et d'autre de l'interface supérieure : cette technique permet la résolution finale de l'ETR en tenant compte des conditions aux limites. Cette méthode (approximation de diffusion et méthode du dipôle) a été introduite par Jensen dans le domaine de la synthèse d'image en 2001 [Jensen01] et a révolutionné la façon de représenter des matériaux translucides. Il est intéressant de noter que ces modèles en synthèse d'image ont été étendus à la représentation de systèmes multicouches. Des relations identiques à celles de l'approche discrète de Kubelka [Donner05] ou de l'approche continue de Kubelka-Munk [Hasan10] ont été obtenues mais dans l'espace fréquentiel (les multiplications de réflectances ou de transmittances devenant des produits de convolution).

Une autre méthode est de discrétiser les flux dans les six directions cartésiennes. L'ETR se simplifie alors en six équations différentielles couplées en x, y et z dont la résolution est souvent numérique, et est plus ou moins simple selon les hypothèses retenues. Cette approche proposée par Chu et Churchill en 1955 [Chu55] a été reprise, avec quelques variantes dans les hypothèses de diffusion, pour des domaines d'application variés : milieu particulaire dans l'infrarouge [Emslie73], tissus biologiques éclairés par un laser [Yoon87], impression sur papier [Mourad01], estimation de l'exposition dans un photoréacteur [Brucato06]. Ce modèle à six flux est une extension à trois dimensions du modèle unidimensionnel de Kubelka-Munk. Il faut noter qu'une simplification plus radicale consiste à utiliser directement le modèle de Kubelka-Munk non en z mais sur la seule dimension radiale afin de déduire la PSF de papiers pour l'impression [Oittinen82]. Cette approche a été largement discutée et adaptée de manière empirique par rapport à des mesures obtenues [Engeldrum95, Arney03, Coppel11].

Nous proposons ici une technique de résolution qui consiste en une généralisation à trois dimensions du modèle discret unidimensionnel de Kubelka. Il s'agit donc de discrétiser les flux à la fois spatialement et en direction. Pour être exploitable facilement, ce modèle nécessite de choisir une distribution régulière des éléments diffusants. L'arrangement choisi
impose le nombre de directions de flux à considérer. En choisissant un réseau cubique, nous proposons ainsi un modèle à six flux avec une discrétisation en x, y et z.

Nous allons dans la partie suivante présenter le principe général du modèle. Dans la partie III, le découpage en six flux de la diffusion simple est explicité avec l'exemple de la fonction de phase d'Henyey-Greenstein. La détermination analytique de la réflectance totale par la méthode à deux flux de Kubelka, en tenant compte de la diffusion latérale est exposée dans la partie IV. Enfin, la répartition de la réflectance sur l'interface supérieure est calculée de manière littérale puis numérique pour un nombre fini d'événements de diffusion (en pratique 20). La technique de résolution est présentée dans la partie V ainsi que les résultats obtenus.
## II. Description générale du modèle

Le tableau 1 résume les principales notations utilisées par la suite.

### Diffusion simple / Probabilités pour une particule unique

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

### Diffusion multiple pour une couche

$r_L$ Probabilité de réflectance totale pour une couche

$t_L$ Probabilité de transmittance totale pour une couche

### Diffusion multiple pour un matériau semi-infini / Probabilités de réflectance sur l'interface supérieure pour N événements de diffusion pris en compte

$$N = 1 \rightarrow \text{diffusion simple}$$

$$N = \infty \rightarrow \text{diffusion multiple totale}$$

$r_N(i, j)$ Probabilité de réflectance au point de coordonnées $(i, j)$ avec $(i, j) \in \mathbb{Z}^2$

$\langle r_N(i, j) \rangle_{|i|+|j|=d}$ Probabilité moyenne des réflectances $r_N(i, j)$ telles que $|i|+|j|=d$

$r_N(d)$ Probabilité de réflectance à la distance $d$ du point d'impact $(0,0)$ avec $d \in \mathbb{N}$

$r_N$ Probabilité de réflectance totale sur l'ensemble de la surface

<caption>**Tableau 1 :** Notations pour les différentes probabilités</caption>
Nous décrivons un matériau diffusant opaque semi-infini en divisant le demi-espace en un réseau cubique discret aux nœuds duquel se trouvent des éléments diffusants supposés identiques (Figure 4-a). Connaissant la concentration volumique $\rho$ de centres diffusants, le pas du réseau est alors égal à $1/\sqrt[3]{\rho}$.

<figure><img src="image_10.png" /><figcaption>**Figure 4** : (a) Réseau cubique au sommet duquel se trouvent les éléments diffusants. Exemple d'un chemin de diffusion avec les probabilités associées.</figcaption></figure>

(b) Probabilités associées à la diffusion simple : absorption $k$, diffusion vers l'avant $f$, diffusion vers l'arrière $b$, diffusion latérale $4l$ avec $k+f+b+4l=1$.

On note $(i, j)$, les coordonnées des éléments diffusants sur l'interface supérieure et on suppose que la lumière est incidente au point $(0,0)$.

Lorsque la lumière rencontre un élément diffusant, elle peut être absorbée selon la probabilité $k$ ou diffusée selon la probabilité $s$ avec $k+s=1$. Pour décrire la diffusion simple, on divise l'espace en six directions (Figure 4-b). La lumière peut être diffusée vers l'avant selon la probabilité $f$, vers l'arrière selon la probabilité $b$ ou vers l'une des quatre directions latérales selon la probabilité $l$. La symétrie généralement supposée de la diffusion simple par rotation autour de la direction d'incidence permet de choisir une probabilité $l$ identique pour les quatre directions latérales. On a alors :

$$f+b+4l=s \qquad (1)$$

A partir du découpage de la diffusion simple, le principe consiste à comptabiliser l'ensemble des chemins possibles et leurs probabilités associées pour que la lumière entrant au point $(0,0)$ en ressorte en un point quelconque $(i, j)$ de l'interface supérieure. Un exemple est donné sur la figure 4-a. On peut de proche en proche déterminer les expressions des reflectances $r_N(i, j)$ après $N$ événements de diffusion. L'obtention de la réflectance $r_\infty(i, j)$ sans approximation ($N = \infty$) ne semble pas possible de manière analytique. Le calcul récursif des expressions littérales de $r_N(i, j)$ est explicité et appliqué pour $N = 20$ dans la partie V tandis que les expressions littérales de $r_S(i, j)$ sont présentées dans l'annexe.
# III. Découpage en six flux de la diffusion simple

Les modèles à deux flux ou à six flux nécessitent préalablement de découper le flux diffusé par une particule unique. Ceci revient à intégrer angulairement la fonction de phase caractéristique de la diffusion simple.

## III.1. Fonction de phase : définitions et propriétés

Le tableau 2 et la figure 5 précisent les notations utilisées pour définir les directions d'incidence **i** et d'observation **o**.

### Vecteurs direction

**i** direction d'incidence

**o** direction d'observation

**n** direction de référence

### Pour chaque direction **v**

$\theta_v$ angle zénithal

$\phi_v$ angle azimutal

$d\omega_v$ angle solide élémentaire autour de la direction **v**

### Produits scalaires

$$ \mu = i \cdot o $$

$$ \mu_i = i \cdot n = \cos \theta_i $$

$$ \mu_o = o \cdot n = \cos \theta_o $$

$\mu_{lim} = \cos \theta_{lim} = \frac{N-2}{N}$ limitation du cône pour un
découpage à N flux (cf. Figure 5-b)

<caption>Tableau 2 : Notations des différentes directions</caption>
<figure><img src="image_11.png" /><figcaption>**Figure 5** : (a) Schéma explicitant les directions d'incidence **i** et d'observation **o** par rapport à une direction de référence **n**. (b) Schéma explicitant les angles solides considérés pour l'incidence, pour la diffusion vers l'avant et pour la diffusion vers l'arrière par rapport à une direction **n** de référence.</figcaption></figure>

L'indicatrice de diffusion $f(\mathbf{i},\mathbf{o})$ d'une particule est définie comme le rapport entre l'intensité $I_d(\mathbf{i},\mathbf{o})$ diffusée dans une direction $\mathbf{o}$ pour une direction incidente $\mathbf{i}$ et le flux total diffusé $F_d$ :

$$f(\mathbf{i}, \mathbf{o}) = \frac{I_d(\mathbf{i}, \mathbf{o})}{F_d} \qquad (2)$$

Par définition, l'indicatrice de diffusion est sans dimension, s'exprime en $\text{sr}^{-1}$ et vérifie la relation suivante :

$$\int_{4\pi} f(\mathbf{i}, \mathbf{o}) d\omega_o = 1 \quad \forall \mathbf{i} \qquad (3)$$

La fonction $f(\mathbf{i},\mathbf{o})$ peut donc être interprétée comme une probabilité de diffusion. La fonction de phase $P(\mathbf{i},\mathbf{o})$ correspond à une normalisation différente, rapport entre l'intensité $I_d(\mathbf{i},\mathbf{o})$ diffusée et l'intensité $F_d / 4\pi$ diffusée pour un comportement isotrope, toutes choses égales par ailleurs. On a alors :

$$P(\mathbf{i}, \mathbf{o}) = 4\pi f(\mathbf{i}, \mathbf{o}) \qquad (4)$$

La fonction de phase est sans dimension, sans unité et :

$$\int_{4\pi} P(\mathbf{i}, \mathbf{o}) d\omega_o = 4\pi \quad \forall \mathbf{i} \qquad (5)$$

La plupart du temps, les particules considérées (typiquement des sphères) présentent une symétrie telle que la fonction de phase ne dépend que du produit scalaire $\mu = \mathbf{i}\cdot\mathbf{o}$. L'équation (5) pour la fonction $P(\mu)$ s'écrit alors :

$$\frac{1}{2} \int_{\mu=-1}^{1} P(\mu) d\mu = 1 \qquad (6)$$
### III.2. Différentes fonctions de phase

La théorie de Mie est l'application des lois de l'électromagnétisme (équations de Maxwell) pour déterminer la diffusion de la lumière par une sphère d'indice de réfraction donné, éclairée par une onde incidente plane [Bohren83]. Pour des particules non absorbantes (d'indice réel) et de diamètres petits par rapport à la longueur d'onde (théorie de Rayleigh), la fonction de phase est relativement isotrope. Plus le diamètre augmente, plus la théorie de Mie prévoit une fonction de phase constituée de nombreux lobes et une diffusion vers l'avant plus importante. Lorsque le diamètre devient plus grand que la longueur d'onde, la lumière diffusée peut être déterminée dans le cadre de l'optique géométrique en considérant les réflexions/réfractions à la surface de la particule transparente [Simonot08]. La diffusion vers l'avant y est majoritaire. Dans une approche radiométrique, on peut aussi déterminer la fonction de phase d'une particule opaque Lambertienne. Dès les origines de la photométrie [Lambert1772], Lambert décrit comment une planète réfléchit la lumière du soleil et donne l'expression de cette fonction de phase permettant typiquement de calculer l'éclairement indirect apporté par la lune sur la terre. Dans ce cas, la rétro-diffusion est majoritaire.

La plupart des matériaux présentant une forte dispersion en forme et en taille des particules diffusantes, ce sont davantage des fonctions de phase décrivant le comportement « moyen » en diffusion simple qui sont utilisées. Il s'agit de fonctions de phase empiriques, décrites par un nombre réduit de paramètres. Une des plus usitées est celle de Henyey-Greenstein [Henyey41] dépendant du seul paramètre g avec $-1 \le g \le 1$ :

$$P(\mu) = \frac{1 - g^2}{(1 + g^2 - 2g\mu)^{3/2}} \quad (7)$$

g est appelé paramètre d'anisotropie et peut varier de $g = -1$ pour une diffusion strictement vers l'arrière à $g = 1$ pour une diffusion strictement vers l'avant. Le cas $g = 0$ correspond à une diffusion isotrope. Pour des valeurs de g opposées, les fonctions de phase sont symétriques par rapport au plan perpendiculaire à la direction d'incidence (Figure 6).

<figure><img src="image_12.png" /><figcaption>**Figure 6** : Diagramme polaire de la fonction de phase d'Henyey-Greenstein (Equation (7)) pour $g = -0,5$, $g = 0$ et $g = 0,5$.</figcaption></figure>
### III.3. Découpage en deux ou six flux

On découpe l'espace en $N$ parties égales d'angle solide $4\pi/N$. On en déduit la valeur $\mu_{lim}$ délimitant les angles solides (Figure 5-b) :

$$ \mu_{lim} = \cos \theta_{lim} = \frac{N-2}{N} \qquad (8) $$

Pour deux flux (N=2), l'espace est découpé en deux hémisphères tels que $\mu_{lim} = 0$.

Pour six flux (N=6), afin de conserver les symétries du réseau cubique, il faudrait en toute rigueur définir des angles solides coniques à base carrée. Pour simplifier les expressions, sans perdre la généralité du propos, nous choisissons pour les angles solides avant et arrière des cônes tels que $\mu_{lim} = 2/3$ (soit $\theta_{lim} \approx 48,2^{\circ}$), tandis que la diffusion latérale se déduit de la relation (1).

Pour une direction incidente collimatée, les fractions de diffusion vers l'avant $(f/s)_c$ et vers l'arrière $(b/s)_c$ se déterminent par simple intégration sur $\mu$ :

$$ (f/s)_c = \frac{1}{2} \int_{\mu=\mu_{lim}}^{1} P(\mu) d\mu \qquad (9) $$

$$ (b/s)_c = \frac{1}{2} \int_{\mu=-1}^{-\mu_{lim}} P(\mu) d\mu \qquad (10) $$

Pour la fonction de phase d'Henyey-Greenstein, les intégrales (9) et (10) conduisent aux expressions analytiques suivantes :

$$ (f/s)_c = \frac{1-g^2}{2g} \left( \frac{1}{1-g} - \frac{1}{(1+g^2-2g\mu_{lim})^{1/2}} \right) \qquad (11) $$

$$ (b/s)_c = \frac{1-g^2}{2g} \left( -\frac{1}{1+g} + \frac{1}{(1+g^2+2g\mu_{lim})^{1/2}} \right) \qquad (12) $$

La figure 7 présente le découpage de la diffusion simple dans le cas de la fonction de phase d'Henyey-Greenstein en deux flux (Figure 7-a) et six flux (Figure 7-b). La fraction de diffusion vers l'avant augmente continument avec le paramètre $g$ tandis que la diffusion latérale est maximale dans le cas isotrope ($g=0$) et diminue continument quand $|g|$ augmente.
<figure><img src="image_13.png" /><figcaption>**Figure 7** : Découpage de la diffusion simple pour la fonction de phase d'Henyey-Greenstein en fonction du paramètre *g* pour une lumière incidente collimatée. (a) Découpage en deux flux (vers l'avant (f/s)$_c$, vers l'arrière (b/s)$_c$) (b) Découpage en six flux (vers l'avant (f/s)$_c$, vers l'arrière (b/s)$_c$ et latéralement 4(l/s)$_c$).</figcaption></figure>

Mais pour les modèles développés, le flux incident n’est pas collimaté mais diffus. La lumière est supposée arriver de manière équiprobable dans le cône d’incidence d’angle solide $4\pi/N$. Le découpage du flux diffusé est plus complexe et nécessite d’exprimer les directions **i** et **o** en relation de l’axe photométrique nommé **n**, selon la relation de géométrie sphérique suivante :

$$ \mu = \mu_i \mu_o + \sqrt{1 - \mu_i^2} \sqrt{1 - \mu_o^2} \cos \varphi \quad (13) $$

avec $\mu_i = i \cdot n$, $\mu_o = o \cdot n$ et $\varphi = \varphi_o - \varphi_i$, l'angle azimutal entre les deux directions **o** et **i**.

La direction **n** est prise comme la direction incidente moyenne. Pour la première diffusion, il s'agit de la normale de la surface supérieure (dirigée vers le bas). Pour les événements de diffusion suivants, il s'agit d'une arrête du réseau cubique.

On introduit les fonctions de diffusion $h_f(\mu_i)$ et $h_b(\mu_i)$, probabilité de diffusion dans l'angle solide vers l'avant, respectivement vers l'arrière, pour une incidence $\mu_i$ donnée.

$$ h_f(\mu_i) = \frac{1}{4\pi} \int_{\mu_o=\mu_{lim}}^{1} \int_{\varphi=0}^{2\pi} P(\mu_i \mu_o + \sqrt{1-\mu_i^2} \sqrt{1-\mu_o^2} \cos \varphi) d\mu_o d\varphi \quad (14) $$

$$ h_b(\mu_i) = \frac{1}{4\pi} \int_{\mu_o=-1}^{-\mu_{lim}} \int_{\varphi=0}^{2\pi} P(\mu_i \mu_o + \sqrt{1-\mu_i^2} \sqrt{1-\mu_o^2} \cos \varphi) d\mu_o d\varphi \quad (15) $$

Enfin, nous obtenons les fractions *f/s* de diffusion vers l'avant et *b/s* de diffusion vers l'arrière en calculant, par intégration angulaire, la valeur moyenne des fonctions de diffusion $h_f(\mu_i)$ et $h_b(\mu_i)$ pour toutes les directions incluses dans l'angle solide d'incidence.

$$ f/s = \frac{\int_{\mu_i=-1}^{-\mu_{lim}} h_f(\mu_i) d\mu_i}{1 - \mu_{lim}} \quad (16) $$
$$b/s = \frac{\int_{\mu_i = -1}^{-\mu_{lim}} h_b(\mu_i) d\mu_i}{1 - \mu_{lim}} \quad (17)$$

On déduit la diffusion latérale en utilisant une nouvelle fois la relation (1).

La figure 8 est l'équivalent pour un flux incident diffus de la figure 7. Les variations avec le paramètre g sont similaires mais les intégrations angulaires supplémentaires pour un flux diffus (Eqs. (14)-(17)) a tendance à atténuer l'anisotropie de la fonction de phase par rapport à l'équivalent pour une incidence collimatée. A titre d'exemple, pour g = 0,5, le découpage en deux flux est de 83% pour la diffusion vers l'avant et 17% pour la diffusion vers l'arrière pour une lumière incidente collimatée. Ce découpage devient 69% et 31% pour un éclairage incident diffus. Dans le cas du six flux, ce rééquilibrage se fait essentiellement en faveur de la diffusion latérale. Pour g = 0,5, le découpage est de 52%, 4% et 44% pour la diffusion respectivement vers l'avant, vers l'arrière et latéralement dans le cas d'une incidence collimatée. Ce découpage devient 41%, 5% et 54% pour un éclairage diffus.

<figure><img src="image_14.png" /><figcaption>**Figure 8** : Découpage de la diffusion simple pour la fonction de phase d'Henyey-Greenstein en fonction du paramètre g pour un flux diffus. (a) Découpage en deux flux (vers l'avant f/s, vers l'arrière b/s) (b) Découpage en six flux (vers l'avant f/s, vers l'arrière b/s et latéralement 4l/s).</figcaption></figure>

Par la suite, nous présentons les résultats pour trois cas typiques :

- Forte diffusion vers l'arrière (g = -0,5) : f = 0,05s, b = 0,41s, 4l = 0,54s

- Diffusion isotrope (g = 0) : f = b = l = s / 6

- Forte diffusion vers l'avant (g = 0,5) : f = 0,41s, b = 0,05s, 4l = 0,54s
# IV. Détermination de la réflectance totale par la méthode à deux flux

La réflectance totale $r_∞$, rapport entre le flux total réfléchi et le flux incident, pourrait se calculer ainsi :

$$r_{\infty} = \sum_{i=-\infty}^{+\infty} \sum_{j=-\infty}^{+\infty} r_{\infty}(i, j) \qquad (18)$$

Mais n'ayant pas accès à la répartition de la réflectance $r_∞(i, j)$, il est nécessaire de procéder autrement pour déterminer la réflectance totale $r_∞$. Pour cela, le réseau cubique (figure 4-a) est décrit comme une superposition de couches identiques (Figure 9).

<figure><img src="image_15.png" /><figcaption>**Figure 9 :** Superposition infinie de couches identiques</figcaption></figure>

La réflectance $r_L$ et la transmittance $t_L$ d'une couche est déterminée selon le schéma de la figure 10 représentant les événements de diffusion successifs ayant lieu dans la couche mais sans préoccupation de la distribution latérale des flux diffusés.

<figure><img src="image_16.png" /><figcaption>**Figure 10 :** Détermination de la réflectance $r_L$ et de la transmittance $t_L$ d'une couche seule (un carré représente le ième événement de diffusion)</figcaption></figure>
Pour un flux incident provenant de la couche, lors du ième événement de diffusion (i ≥ 2), la lumière émerge en réflectance dans une direction latérale selon une probabilité l. Il en est de même pour la transmittance. La lumière diffusée dans les quatre autres directions selon la probabilité f+b+2l participe aux événements de diffusion ultérieurs à l'intérieur de la couche. En sommant l'ensemble des probabilités de rétrodiffusion, respectivement de transmission, on obtient la réflectance, respectivement la transmittance, de la couche :

$$
r_L = b + 4l^2 \sum_{i=2}^{+\infty} (b+f+2l)^{i-2} = b + \frac{4l^2}{1-(b+f+2l)} \quad (19)
$$

$$
t_L = f + 4l^2 \sum_{i=2}^{+\infty} (b + f + 2l)^{i-2} = f + \frac{4l^2}{1 - (b + f + 2l)} \quad (20)
$$

Le matériau est finalement modélisé par une succession infinie de couches identiques dont on distingue la couche supérieure de réflectance $r_L$ et transmittance $t_L$ du reste des couches de réflectance notée $r_{\infty-}$. La figure 9 montre les réflexions multiples entre la couche supérieure et le fond diffusant, conduisant à une série géométrique similaire au modèle à deux flux de Kubelka [Kubelka48].

$$
r_{\infty} = r_{L} + t_{L}^{2} r_{\infty-} \sum_{i=2}^{+\infty} (r_{L} r_{\infty-})^{i-2} = r_{L} + \frac{t_{L}^{2} r_{\infty-}}{1 - r_{L} r_{\infty-}} \quad (21)
$$

Comme le nombre de couches est infini, la réflectance $r_\infty$ n'est pas influencée par l'addition
ou la soustraction d'une couche, donc $r_\infty = r_{\infty-}$. On aboutit à l'équation suivante :

$$
r_{\infty}^{2} - 2Ar_{\infty} + 1 = 0 \tag{22}
$$

avec

$$
A = \frac{1 + r_L^2 - t_L^2}{r_L} \tag{23}
$$

et dont la seule solution positive et inférieure à 1 est :

$$
r_{\infty} = A - \sqrt{A^2 - 1} \quad (24)
$$

La figure 11 présente l'évolution de la réflectance totale $r_\infty$ et de la diffusion simple $r_1=b$
avec la probabilité d'absorption $k$ pour les trois situations choisies. De manière attendue, la
réflectance $r_\infty$ augmente quand, pour une répartition de diffusion donnée, l'absorption
diminue et quand, pour une probabilité d'absorption donnée, la rétrodiffusion ($r_1=b$) augmente. La figure 11 permet en outre de quantifier la part croissante de la diffusion
multiple quand l'absorption diminue.
Forte diffusion vers l'arrière

$f = 0,05s, b = 0,41s, 4l = 0,54s$

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

Diffusion isotrope

$f = b = l = s/6$

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

Forte diffusion vers l'avant

$f = 0,41s, b = 0,05s, 4l = 0,54s$

<figure><img src="image_19.png" /><figcaption>**Figure 11** : Réflectance totale $r_\infty$ (et réflectance due à la diffusion simple $r_1 = b$) en fonction de la probabilité d'absorption $k$ pour trois situations différentes</figcaption></figure>
# V. Résolution et résultats

## V.1. Méthode de résolution

Une fois déterminées ou choisies les probabilités $f, b, l$ de diffusion simple, le principe du calcul consiste à propager de proche en proche les probabilités de diffusion. On peut déterminer ainsi les expressions littérales des réflectances $r_N(i, j)$ de manière récursive. A titre d'exemple, l'annexe donne les expressions obtenues pour $r_5(i, j)$ et nous présentons par la suite les résultats obtenus pour $N=20$.

Pour passer de l'étape $n-1$ à l'étape $n$, cela implique un nombre de nœuds proportionnel à $n^3$. Le temps de calcul total à l'étape $N$ varie donc en $\sum_{n=1}^{N} n^3 = (N(N+1)/4)^2$ soit un ordre de grandeur variant en $N^4$. Cette variation se vérifie pour un traitement entièrement numérique du problème, c'est un peu moins vrai pour le traitement littéral employé ici : la gestion de formules littérales de plus en plus complexes rend difficile l'estimation du temps de calcul. Bien sûr la détermination des expressions littérales est sensiblement plus longue qu'un calcul entièrement numérique à partir des valeurs choisies pour les paramètres $k, b, f$ et $l$. Mais une fois ces expressions littérales obtenues, il suffit de les appliquer afin d'obtenir les valeurs numériques $r_N(i, j)$ pour n'importe quel jeu de paramètres choisi.

Pour estimer la précision obtenue en fonction de $N$, la réflectance totale $r_\infty$ (Eqs. (19)-(24)) est comparée à la réflectance totale approchée $r_N$ pour $N$ événements de diffusion pris en compte :

$$r_N = \sum_{i=-N}^{N} \sum_{j=-N}^{N} r_N(i, j) \qquad (26)$$

La figure 12 présente cette comparaison en fonction de la probabilité d'absorption dans le cas d'une diffusion isotrope ($f=b=l=s/6$). Un écart de moins de 1% sur la réflectance est obtenu pour $k > 0,75$ pour la diffusion simple $r_1$, $k > 0,3$ pour $r_5$ et $k > 0,08$ pour $r_{20}$. Le cas $N=20$ est donc suffisant pour obtenir la réflectance totale avec une très bonne précision même pour des faibles valeurs d'absorption mais reste insuffisant pour traiter le cas limite d'un matériau parfaitement non absorbant ($k=0$).

<figure><img src="image_20.png" /><figcaption>**Figure 12** : Réflectance totale $r_\infty$ et valeur approchée $r_N$ (avec $N=1, 5$ ou $20$ événements de diffusion) en fonction de la probabilité d'absorption $k$ pour une diffusion isotrope $f=b=l=s/6$</figcaption></figure>
Les exemples présentés par la suite considèrent N=20 événements de diffusion et une probabilité d'absorption k=0,1.

## V.2. Vers une répartition continue à symétrie circulaire

Par construction du modèle, la répartition de la réflectance $r_N(i, j)$ est discrète. Autrement dit, la réflectance est en théorie nulle en dehors des points $(i, j)$ du réseau. Une interpolation de la réflectance est nécessaire afin d'obtenir une répartition continue. En pratique, nous choisissons de traiter à part la diffusion simple (probabilité $b$) en conservant son aspect impulsionnel tandis que la tâche de diffusion liée à la diffusion multiple ($N \ge 2$) est interpolée.

En outre, le choix d'un réseau cubique induit une tâche de diffusion sur la surface supérieure qui n'est pas à symétrie circulaire comme attendu mais sous forme d'une figure dont les dimensions dans les directions diagonales sont plus petites que celles le long des axes (Figure 13). En effet, pour atteindre le point de coordonnées $(i, j)$, la distance minimale parcourue par la lumière est $|i| + |j|$ et non $\sqrt{i^2 + j^2}$.

<figure><img src="image_21.png" /><figcaption>**Figure 13** : Exemple de répartition de la réflectance $r_{20}(i,j)$ obtenue pour $k=0,1$ et une diffusion simple isotrope ($f=b=l=s/6$). Cette répartition en principe discrète est ici interpolée. Image en niveau de gris (échelle logarithmique), contour en rouge pour une réflectance de 0,01. Le réseau cubique choisi implique une tâche de diffusion qui n'est pas à symétrie circulaire.</figcaption></figure>

Pour résoudre ce problème, une solution envisageable serait de reprendre l'ensemble du calcul en choisissant un autre réseau, par exemple hexagonal compact, qui respecterait davantage la symétrie circulaire que le réseau cubique. L'inconvénient d'une telle solution est qu'il faudrait découper la diffusion simple en un nombre beaucoup plus important de directions (douze directions à définir pour un réseau hexagonal compact contre seulement six pour le réseau cubique).
Nous choisissons de conserver la simplicité du réseau cubique et de contraindre la symétrie circulaire en définissant la réflectance $r_N(d)$ où l'entier positif $d$ correspond à la distance au point d'incidence (0,0). On associe à $r_N(d)$, la moyenne des probabilités $r_N(i, j)$ telles que $|i|+|j|=d$, et notée $\langle r_N(i, j) \rangle_{|i|+|j|=d}$ :

$$r_N(d) = C_N \langle r_N(i, j) \rangle_{|i|+|j|=d} \qquad (27)$$

où $C_N$ est une constante de normalisation afin que la réflectance totale $r_\infty$ soit conservée dans les deux écritures. Or avec le découpage cartésien, chaque probabilité $r_N(i, j)$ est associée à une surface carrée centrée en $(i, j)$ et d'aire unité (Figure 14.a). Avec le découpage à symétrie circulaire, les probabilités $r_N(d)$ sont associées à des couronnes circulaires de largeur unité telles que représentées sur la figure 14.b, sauf pour la probabilité $r_N(0)$ à laquelle est associée le disque de rayon $1/2$.

<figure><img src="image_22.png" /><figcaption>**Figure 14 :** (a) Découpage cartésien associé aux réflectances $r_N(i,j)$ (b) Découpage à symétrie circulaire associé aux réflectances $r_N(d)$</figcaption></figure>

La constante de normalisation est alors :

$$C_N = \frac{r_\infty - b}{\pi \left( \frac{1}{4} (r_N(0,0) - b) + \sum_{d=1}^{N} 2d \langle r_N(i, j) \rangle_{|i|+|j|=d} \right)} \qquad (28)$$

La diffusion simple étant traitée à part, l'expression (27) est différente pour $d=0$ :

$$r_N(0) = C_N (r_N(0,0) - b) + b\delta(0) \qquad (29)$$

où $b\delta(0)$ est l'impulsion de probabilité $b$ au point d'impact ($d=0$) due à la diffusion simple.

Il faut noter qu'il nous a semblé intéressant de normaliser en imposant la réflectance totale $r_\infty$ car celle-ci est connue (Eqs. (19)-(24)), plutôt que la réflectance totale approchée $r_N$ (Eq.(26)).

A titre d'exemple, la figure 15 présente une comparaison de l'évolution de $r_N(2,0)$ et $r_N(2)$ en fonction de $N$. Il faut au moins trois événements de diffusion ($N \ge 3$) pour que les
réflectances $r_N(2,0)$ et $r_N(2)$ soient non nulles. L'évolution des réflectances avant et après normalisation diffère : $r_N(2,0)$ croît à chaque fois que $N$ est impair (condition nécessaire pour que la lumière diffusée émerge en un point $(i,j)$ avec $i+j$ pair). La normalisation (Eqs. (27)-(28)) induit d'une part que la valeur asymptotique est légèrement différente, d'autre part que la réflectance $r_N(2)$ n'est plus strictement croissante mais oscille autour de la valeur asymptotique. La figure 14 montre ainsi un avantage supplémentaire apporté par la normalisation : la convergence plus rapide avec $N$.

<figure><img src="image_23.png" /><figcaption>**Figure 14** : Evolution de $r_N(2,0)$ et $r_N(2)$ en fonction de $N$ pour une diffusion simple isotrope ($f=b=l=s/6$) et $k=0,1$.</figcaption></figure>

### V.3. Distributions de la réflectance $r_{20}(d)$

Pour déterminer les réflectances $r_N(d)$, le protocole de calcul est finalement le suivant :

- détermination des relations littérales de $r_N(i, j)$ pour $N$ donné (ici $N=20$);

- valeurs numériques de $r_N(i, j)$ pour les paramètres $k, b, f$ et $l$ choisis ou déterminés selon le calcul décrit dans la partie III;

- valeurs numériques de $r_N(d)$ selon les équations (27)-(29).

La figure 15 donne la répartition $r_{20}(d)$ pour trois situations typiques mais présentant la même probabilité d'absorption $k=0,1$ par particule. D'un point de vue quantitatif, la différence est surtout liée à la diffusion simple de probabilité $b$. Mais il est notable que les allures, amplitude et étendue, de la tâche de diffusion sont sensiblement différentes pour les trois situations traitées. Ainsi une réflectance à longue distance du point d'impact est favorisée dans le cas d'une forte diffusion vers l'avant (cas de la plupart des matériaux biologiques par exemple), tandis que c'est le contraire pour la réflectance à courte distance du point d'impact.
<figure><img src="image_24.png" /><figcaption>**Figure 15** : Répartition de la réflectance $r_{20}(d)$ pour $k=0,1$ et trois situations. Valeurs calculées pour $d$ entiers et interpolées pour $d$ non entiers. Les probabilités de diffusion simple $b$ et de réflectance totale $r_\infty$ sont rappelées pour comparaison. Les taches de diffusion sont centrées sur le centre des carrés dont les côtés ont une distance normalisée de 10. Représentation avec la diffusion simple (1ère ligne) et une échelle logarithmique en niveau de gris, sans la diffusion simple (2ème ligne) et une échelle linéaire en niveau de gris et courbes correspondantes (3ème ligne).</figcaption></figure>
Rappelons que la distance *d* est ici sans dimension et qu'il faut diviser par la racine cubique de la densité volumique des éléments diffusants pour donner une échelle métrique. Ainsi, toutes choses égales par ailleurs, la densité volumique n'influence pas l'allure de la tâche de diffusion mais est reliée à son étendue : moins le matériau est dense, plus la tâche de diffusion est large.

## Conclusion

Nous avons proposé un modèle de diffusion multiple complet et alternatif aux modèles existants. La méthodologie est décomposée en trois parties quasiment indépendantes : le découpage de la diffusion simple en six flux, le calcul de la réflectance totale sans approximation et enfin la détermination de la distribution de la réflectance sur l'interface supérieure du matériau. Ce dernier calcul donne une valeur approchée mais satisfaisante à condition de considérer un nombre suffisamment important d'événements de diffusion (et par conséquent un temps de calcul important) et/ou une probabilité d'absorption par particule qui ne soit pas trop faible.

Remarquons que des extensions du modèle peuvent être proposées pour une distribution des éléments diffusants autre que cubique (mais cela implique de redéfinir également le découpage de la diffusion simple), pour des éléments diffusants non identiques, avec des probabilités de diffusion différentes selon les directions latérales... Si nous avons choisi les hypothèses présentées ici (réseau cubique d'éléments diffusants identiques présentant la même probabilité de diffusion dans les quatre directions latérales), ce n'est pas seulement parce que la mise en œuvre du modèle y est la plus simple mais davantage parce que ces hypothèses minimisent le nombre de paramètres décrivant le modèle rendant ce dernier plus directement utilisable.

Le modèle a été présenté ici pour une couche infinie. Il peut être généralisé, sans difficulté apparente, à une épaisseur finie et permettre ainsi d'obtenir la répartition de la lumière diffusée sur les deux faces du matériau, autrement dit en réflectance et en transmittance. L'effet d'une interface plane (i.e. d'un changement d'indice de réfraction) peut également être pris en compte en appliquant la correction de Saunderson [Saunderson42]. Par rapport aux méthodes à deux flux, ce modèle à six flux présente l'indéniable avantage de permettre l'accès à la répartition latérale du flux réfléchi, sans hypothèse restrictive sur la diffusion simple. En revanche, tout comme les méthodes à deux flux, le modèle ne décrit pas la répartition angulaire du flux réfléchi : la lumière est supposée être diffusée de manière équiprobable dans l'angle solide de sortie considéré. Pour s'affranchir de cette difficulté, une solution est de découpler la modélisation de la diffusion de volume et de la diffusion de surface. Le comportement Lambertien du modèle présenté est en effet davantage adapté à la diffusion de volume tandis qu'un modèle de BRDF par exemple de type microfacette [Simonot09] peut être utilisé pour la diffusion de surface, qui est souvent le principal responsable de la répartition angulaire du flux réfléchi.
# Annexe : Expressions littérales de $r_5(i, j)$

Les expressions littérales de $r_N(i, j)$ peuvent s'écrire sous la forme d'une somme de termes $Kf^{i_f} b^{i_b} l^{i_l}$ signifiant qu'il y a K chemins ressortant au point $(i, j)$ et ayant subi $i_f$ diffusions vers l'avant, $i_b$ diffusions vers l'arrière et $i_l$ diffusions latérales.

Le tableau suivant donne les expressions littérales de $r_5(i, j)$ pour $(i, j) \in [0, 4]^2$. Les expressions dans les autres quadrants se déduisent par symétrie et $r_N(i, j) = 0$ quand $|i| + |j| \ge N$.

<table><tbody><tr><td rowspan="5">j</td><td>4</td><td>$f^3 l^2$</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3</td><td>$f^2 l^2$</td><td>$l^3 (2f^2 + fl)$</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2</td><td>$r_5 (0,2)$</td><td>$l^3 (2f + l)$</td><td>$2l^3 (f^2 + fl + l^2)$</td><td>0</td><td>0</td></tr><tr><td>1</td><td>$r_5 (0,1)$</td><td>$r_5 (1,1)$</td><td>$l^3 (2f + l)$</td><td>$l^3 (2f^2 + fl)$</td><td>0</td></tr><tr><td>0</td><td>$r_5 (0,0)$</td><td>$r_5 (1,0)$</td><td>$r_5 (2,0)$</td><td>$f^2 l^2$</td><td>$f^3 l^2$</td></tr><tr><td colspan="2"></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td colspan="7" style="text-align: center;">i</td></tr></tbody></table>

avec

$$r_5 (0,0) = b + 4bl^2 + bf^2 + bf^4 + 8bf^2l^2 + 12bl^4 + 8b^2l^3 + 12b^2fl^2 + b^3f^2 + 4b^3l^2 + 8fl^4 + 8l^5$$

$$r_5 (1,0) = r_{\infty}^{(5)} (0,1) = l^2 (1 + 4bf + 4bl + b^2 + f^2 + 2l^2)$$

$$r_5 (2,0) = r_{\infty}^{(5)} (0,2) = l^2 (f + 4bf^2 + 4bfl + 3bl^2 + 2b^2f + f^3 + 4fl^2 + 4l^3)$$

$$r_5 (1,1) = 2l^3 (1 + 6bf + 5bl + 4fl + 2b^2 + f^2)$$
Références

[Arney03] J.S. Arney, J. Chauvin, J. Nauman, P.G. Anderson, *Kubelka-Munk Theory and MTF of Paper*, J. Imaging Sci. Technol. **47** (2003) 339-345.

[Bohren83] C.F. Bohren, D.R. Huffman, *Absorption and scattering of light by small particles*, Wiley-Interscience (1983).

[Born99] M. Born, E. Wolf, *Principle of Optics*, Pergamon, 7th Edition (1999).

[Brucato06] A. Brucato, A.E. Cassano, F. Grisafi, G. Montante, L. Rizzuti, G. Vella, *Estimating Radiant Fields in Flat Heterogeneous Photoreactors by the Six-Flux Model*, American Institute of Chemical Engineers **52** (2006) 3882-3890.

[Chandrasekhar60] S. Chandrasekhar, *Radiative transfert*, Dover, New-York (1960).

[Chu55] C. Chu, S. Churchill, *Numerical solution of problems in multiple scattering of electromagnetic radiation*, J. Phys. Chem. **59** (1955) 855-863.

[Coppel11] L.G. Coppel, M. Neuman, P. Edström, *Lateral light scattering in paper - MTF simulation and measurement*, Optics Express **19** (2011) 25181-25187.

[Donner05] C. Donner, H.W. Jensen, *Light diffusion in multi-layered translucent materials*, In SIGGRAPH '05: ACM SIGGRAPH 2005, New York (2005) 1032-1039.

[Emslie73] A.G. Emslie, J.R. Aronson, *Spectral reflectance and emittance of particulate materials 1:Theory*, Appl. Opt. **12** (1973) 2563-2572.

[Engeldrim95] P.G. Engeldrum, B. Pridham, *Application of Turbid Medium Theory of Paper Spread Function Measurements*, TAGA Proceedings (1995) 339-352.

[Farell92] T. J. Farell, M.S. Patterson, B. Wilson, *A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo*, Med. Phys. **19** (1992) 879-888.

[Hašan10] M. Hašan, M. Fuchs, W. Matusik, H. Pfister, S. Rusinkiewicz, *Physical Reproduction of Materials with Specified Subsurface Scattering*, ACM Transactions on Graphics (Proc. SIGGRAPH) **29** (2010) 4 (10pp).

[Hébert07] M. Hébert, R. Hersch, J.-M. Becker, *Compositional reflectance and transmittance model for multilayer specimens*, J. Opt. Soc. Am. A **24** (2007) 2628-2644.

[Hébert08] M. Hébert, J.-M. Becker, *Correspondence between continuous and discrete two-flux models for reflectance and transmittance of diffusing layers*, J. Opt. A: Pure Appl. Opt. **10** (2008) 035006 (6pp).

[Henyey41] L.G. Henyey, J.L. Greenstein, J.L., *Diffuse Radiation in the Galaxy*, Astrophys. J. **93** (1941) 70-83.
[van de Hulst81] H.C. van de Hulst, *Light scattering by small particles*, Dover Publications (1981) 200–227.

[Ishimaru78] A. Ishimaru, *Wave Propagation and Scattering in Random Media*, Academic Press, New York (1978).

[Jensen01] H.W. Jensen, S.R. Marschner, M. Levoy, P. Hanrahan, *A practical model for subsurface light transport*, In Proceedings of ACM SIGGRAPH 2001, ACM Press/Addison-Wesley Publishing Co., New York, Computer Graphics Proceedings (2001) 511–518.

[Kubelka31] P. Kubelka, F. Munk, *Ein Beitrag zur Optik der Farbanstriche*, Zeitschrift für technische Physik **12** (1931) 593-601.

[Kubelka48] P. Kubelka, *New contributions to the optics of intensely light-scattering material, part I*, J. Opt. Soc. Am. **38** (1948) 448-457.

[Lambert1772] J.H. Lambert, *Beschreibung einer mit dem Calaunischen Wachse ausgemalten Farbenpyramide*, (Haude and Spener ed., Berlin 1772)

[Mandelis90] A. Mandelis, F. Boroumand, H. van den Bergh, *Quantitative diffuse reflectance spectroscopy of large powders: the Melamed model revisited*, Appl. Opt. **29** (1990) 2853-2860.

[Melamed63] N.T. Melamed, *Optical properties of powders: Part I. Optical absorption coefficients and the absolute value of the diffuse reflectance*, Journal of Applied Physics **34** (1963) 560-570.

[Mourad01] S. Mourad, P. Emmel, K. Simon, R. Hersch, *Extending Kubelka-Munk's Theory with Lateral Light Scattering*, In D. Weiss and M. Yuasa, editors, NIP 17: International Conference on Digital Printing Technologies, Fort Lauderdale, Florida, September October 5, (2001) 469-473.

[Mudgett71] P.S. Mudgett, L.W. Richards, *Multiple scattering calculations for technology*, App. Opt. **10** (1971) 1485-1502.

[Nicodemus77] F.E. Nicodemus, J.C. Richmond, J. J. Hsia, I. W. Ginsber, T. Limperis, *Geometrical consideration and nomenclature for reflectance*, J. Res. Natl. Bur. Stand. **160** (1977) 1-52.

[Oittinen82] P. Oittinen, *Limits of microscopic print quality*, in Advances in Printing Science and Technology **16**, W. H. Banks, ed. Pentech, London (1982) 121–128.

[Prahl95] S.A. Prahl, *The adding-doubling method* in “Optical-thermal response to laser irradiated tissue”, A. J. Welch and M. J. C. van Gemert Ed. (New York, 1995) 101-129.

[Rogers97] G. Rogers, *Optical dot gain in a halftone print*, J. Imaging Sci. Technol. **41** (1997) 643-656.

[Saunderson42] J.L. Saunderson, *Calculation of the color pigmented plastics*, J. Opt. Soc. Am. **32** (1942) 727-736.

[Shkuratov99] Y.G. Shkuratov, L. Starukhina, H. Hoffmann, G. Arnold, *A model of spectral albedo of particulate surfaces: implication to optical properties of the Moon*, Icarus **137** (1999) 235-246.
[Simon03] K. Simon, B. Trachsler, *A Random Walk Approach for Light Scattering in Material*, Discrete Mathematics and Theoretical Computer Science AC (2003) 289–300.

[Simonot08] L. Simonot, M. Hébert, R. Hersch, H. Garay, *Ray scattering model for spherical transparent particles*, JOSA A **25** (2008) 1521-1534.

[Simonot09] L. Simonot, *A photometric model of diffuse surfaces described as a distribution of interfaced Lambertian facets*, Applied Optics **48** (2009) 5793-5801.

[Stamnes88] K. Stamnes, S. Chee Tsay, W. Wiscombe, K. Jayaweera, *Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layer media*, Appl. Opt. **27** (1988) 2502-2510.

[Vöge07] M. Vöge, K. Simon, *The Kubelka-Munk and Dyck paths*, Journal of Statistical Mechanics: Theory and Experiment, (2007) P02018.

[Wakeshima68] H. Wakeshima, T. Kunishi and S. Kaneko, *Light scattering in paper and its effect on halftone reproduction*, J. Opt. Soc. Amer. **58** (1968) 272-273.

[Yoon87] G. Yoon, A.J. Welch, M. Motamedi, M.J.C. van Gemert, *Development and application of three-dimensional light distribution model for laser irradiated tissue*, IEEE Journal of Quantum Electronics **23** (1987) 1721-1732.

[Yule67] J.A.C. Yule, D.J. Howe, J.H. Altman, *The Effect of the Spread-Function of Paper on Halftone Reproduction*, Tappi Journal **50** (1967) 337-344.