Unity Laurent Belcour

E icient Rendering of Layered Materials using an Atomic Decomposition with Statistical Operators

Keywords: CCS Concepts:, Computing methodologies →Re ectance modeling; Layered Materials, Statistical Analysis, Real-Time Rendering ACM Reference format:

Technologies

. We introduce a BSDF model to render plane-parallel layered materials using an analysis of the directional statistics of light interaction with microfacet geometry and participating media. Our model closely matches the reference and supports an arbitrary number of textured layers while being energy conserving, free from heavy per-material precomputation, and compatible with real-time constraints.

We derive a novel framework for the e cient analysis and computation of light transport within layered materials. Our derivation consists of two steps. First, we decompose light transport into a set of atomic operators that act on its directional statistics. Speci cally, our operators consist of re ection, refraction, sca ering, and absorption, whose combinations are su cient to describe the statistics of light sca ering multiple times within layered structures. We show that the rst three directional moments (energy, mean and variance) already provide an accurate summary. Second, we extend the adding-doubling method to support arbitrary combinations of such operators e ciently. During shading, we map the directional moments to BSDF lobes. We validate that the resulting BSDF closely matches the ground truth in a lightweight and e cient form. Unlike previous methods, we support an arbitrary number of textured layers, and demonstrate a practical and accurate rendering of layered materials with both an o ine and real-time implementation that are free from per-material precomputation.

INTRODUCTION

Elaborate shading models allow the realistic reproduction of the diversity and complexity of surface appearance found in nature or man-made objects. A simple way to enrich a set of materials is to coat them with one or many strata of glazing or clear coat.

is is typically what happens when objects are wet or paints are glazed. For example, applying a clear coat on a rough metal will create a distinctive hazy look. Unfortunately, coatings increase the di culty of modeling the resulting appearance and design e cient shading models for rendering. Up to now, no one has provided a comprehensive model of the Bidirectional Sca ering Distribution Function (BSDF) of coated materials. In this work, we focus on a subset of coatings called layered materials [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF]].

Layered Materials. Layered materials assume that the coating layers are plane parallel with respect to the shading normal and separated by rough interfaces. It is also assumed that geometrical optics fully describes the appearance. us, the thicknesses of layers is larger than the wavelengths of visible light (removing wave optics from our study). However, the number of layers can be arbitrary and thus light can bounce many times before exiting the material, which makes evaluation di cult. O ine productions [START_REF] Hery | Pixar s Foundation for Materials[END_REF][START_REF] Kulla | Revisiting Physically Based Shading at Imageworks[END_REF][START_REF] Langlands | Physically Based Shader Design in Arnold[END_REF] and the real-time community [START_REF] Drobot | Practical Multilayered Materials in Call of Duty: In nite Warfare[END_REF]] have expressed interest in those appearance. However, they limit themselves to only two layers evaluated with a single re ection per layer (that is paths directly re ected by the top layer (denoted R) and paths transmi ed by the top layer and re ected by the base layer (denoted T RT )). us, more general -yet e cient -models are needed by the community. Fig. 2. We express the directional statistics (energy, mean and variance) of a layered BSDF in the projected plane (a). Instead of computing the complete transport, we track statistical summary at each step. We first study the impact on directional statistics of di erent atomic operations, i.e. reflection, refraction, sca ering, and absorption. For example, rough reflection increases the variance of the incident radiance (b). We evaluate multiple sca ering between layers by combining those atomic operators with a new adding-doubling algorithm working on those statistics (c). Finally, we instantiate multiple BRDF lobes from those statistics to approximate the entire layered BSDF (d).

Technical Di culties.

A key di culty is to provide a realistic model that works with an arbitrary number of layers (possibly textured), accounts for multiple sca ering, is energy conserving, requires li le storage, has a short precomputation time, supports good importance sampling and is symmetric with respect to light transport evaluation (to be compatible with bidirectional rendering techniques). In the pursuit of some of those goals, researchers o en rely on heavy precomputation and/or tabulate much of the data [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF] to faithfully reproduce the target appearance or devise approximative light transport schemes [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] that unfortunately miss some important e ects. e main objective of our work is to reach a quality comparable to tabulated methods, while providing a lightweight solution that is e cient like the approximative methods. While we are inspired by how Jakob et al. [2014] compute multiple sca ering between layers, we use a lighter representation. Our insight is that light interaction with individual layers is simple and o en results in only stretches, compressions, and blurs of the directional distribution of light. It is the combination of interactions in the structure that creates complex appearance. us, our idea is to work at the level of individual layers where light-ma er interaction is described with simple transport operators (re ection, refraction, absorption, and sca ering) to keep derivation tractable and to use an e cient combination algorithm to build complex sca ering functions.

We make the following contributions:

• We introduce a statistical framework to estimate the energy, mean, and variance of the BSDF for any layered material con guration (Fig. 2(a)). Our framework builds upon the statistics of projected directions (Section 3) • We provide a set of atomic operators (re ection, refraction, absorption, and sca ering) within our framework and derive close approximations of light interaction with interfaces or media (Section 4 and Fig. 2(b)). • We derive a new method, similar to the one of [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF], that works on statistics and allows us to eciently evaluate the directional statistics due to multiple sca ering in the layered structure (Section 5, and Fig. 2(c)). • Using our framework, we implement a new BSDF models for o ine and real-time rendering of layered materials that is energy conserving, requires li le storage, is free from per-material precomputation, and has good importance sampling (Section 6, and Fig. 2(d)). Our o ine BSDF model is also symmetric.

Assumptions & Limitations. In this work, we assume that all interfaces are composed of specular microfacets following the GGX distribution. Due to our approximations, our model accurately reproduces the ground truth for low and moderate roughness only (see Section 7.3). Our current formulation does not support other surface models, e.g. Lambertian di use layers.

PREVIOUS WORK

Specialized models. Many models were derived for a speci c conguration of a layered structure. For example, [START_REF] Dorsey | Modeling and Rendering of Metallic Patinas[END_REF] introduced a model for aging of metals, [START_REF] Jensen | Rendering of wet materials[END_REF] derived a model accounting for a water layer on top of objects, [START_REF] Stam | An illumination model for a skin layer bounded by rough surfaces[END_REF] derived a BSDF model for human skin by layering a dielectric surface on a participating medium, [START_REF] Ershov | Rendering Pearlescent Appearance Based On Paint-Composition Modelling[END_REF] derived a model of akes in a coating, [START_REF] Dai | e Dual-microfacet Model for Capturing in Transparent Slabs[END_REF] modeled the refraction of light by two parallel rough microfacets, etc. While e cient, those models can be di cult to extend beyond their original purpose. In contrast, we propose a general formulation, that is not restricted to a particular layered con guration.

Spatial Di usion.

It is also possible to model the di usion of light in highly di using layered materials [START_REF] Donner | Light di usion in multi-layered translucent materials[END_REF]. is usually results in a Bidirectional Subsurface Sca ering Re ectance Distribution Function (BSSRDF) model as the spatial di usion of light in the medium is accounted for. We do not try to reproduce this e ect and restrict ourselves to the case in which the approximation that light enters and exits the layered structure at the same position is valid.

General models. Some models reproduce the appearance of stratied materials without the explicit evaluation of all light transport. [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] used strati ed microfacets models [START_REF] Walter | Microfacet Models for Refraction through Rough Surfaces[END_REF]] to model transmissive rough interfaces but simpli ed the transport in the layered structure by refracting query rays with respect to a single microfacet. is unfortunately fails to correctly capture the blur due to a stack of layers as shown in Fig. 13. [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF] introduced a numerical scheme to compute and store a tabulated representation of layered BRDFs that accounts for multiple sca ering between layers using the adding-doubling method [van de [START_REF] Van De Hulst | Multiple Light Sca ering[END_REF]. However, it is impractical for production scenarios or real-time rendering. One of its drawbacks is that it is not memory bound as it relies on a Fourier decomposition of the BSDFs. Consequently, a single BSDF can weight more than a gigabyte in their framework.

Real-time models. To gain e ciency, it is possible to restrict layered BSDFs to two rough interfaces and approximate the resulting transport by two BSDF lobes at the shading point (denoted bi-model here). A rst BSDF lobe accounts for the rst interface (or clearcoat) and the second lobe accounts for the remaining interfaces. For example, [START_REF] Elek | Layered Materials in Real-Time Rendering[END_REF] uses [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] microfacet interfaces with further simpli cations. Blur induced by the clear-coat is approximated by taking the max roughness. [START_REF] Guo | Rendering in Transparent Layers with Extended Normal Distribution Functions[END_REF] use two von Mises-Fisher (vMF) lobes at the surface of the object. ey approximate the convolution of vMF lobes to model the impact of the clear-coat. vMF lobes are unlikely to be applied in production as they miss the heavy tails required to model metals. Furthermore, speci c precomputations are required (e.g. precomputed environment maps and area lights) for real-time rendering. Similarly, De Rousiers et al. [2011] approximated the refraction of two successive rough interfaces using vMF lobes, sharing the same drawbacks.

Summary. In this work, we provide a new model for layered materials that permits real-time usage but without sacri cing much of the delity with respect to the reference light transport. We were inspired by the use of the adding-doubling algorithm [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF]] but we forbade the use of any parameter-dependent or per-material time-consuming tabulation. Instead, we developed a statistical analysis of light interaction with the layered structure and tailored the adding-doubling method to track the energy, mean, and variance of the BSDF. To do so, we model the interactions of light in a layered structure as operators on statistics and apply the adding-doubling method to e ciently evaluate multiple sca ering.

anks to all this, we apply our solution to an arbitrary number of textured layers and show that, for low and moderate roughnesses, our solution is close to the reference appearance and energy conserving, but at the same time lightweight and e cient. We further show how to e ciently importance sample our model and how to make it symmetric. We demonstrate our model tailored to both o ine and real-time scenarios.

DIRECTIONAL STATISTICS

Our insight is to estimate the directional statistics of a layered material for a given incidence and to inject those statistics into a BSDF for e cient shading. In this work, we restrict our statistical study to the energy e, the mean µ and the variance σ of a BSDF lobe 1 . We focus on the GGX microfacet model [START_REF] Walter | Microfacet Models for Refraction through Rough Surfaces[END_REF]], but our ideas apply to other BSDF models. We build on the property that, for isotropic microfacet BSDF models like GGX or Beckmann, the bijective mapping from the parameters of the BSDF to its directional 1 We deliberately omit the square on σ for be er readability. We study the mean and variance of the outgoing distribution in this space. We show in the right inset the variance of a GGX lobe at normal incidence with respect to roughness (in blue) as well as the variance of light bouncing twice on parallel surfaces with the same roughness (in red).

statistics enables us to instantiate a BSDF lobe from the knowledge of its energy, mean, and variance.

Statistics in the Projected Plane. To compute directional mean and variance, we need a parametrization. We use the orthographic projection of 3D coordinates of directions to the (x, ) plane (see Fig. 3). at is, we compute the mean and variance with respect to [u, ] = ω x , ω , where a x is the x-component of vector a. We experienced that in this space, an isotropic GGX peak is close to a circular shape. Löw et al. [2012, see Fig. 5] reported a similar behavior in the MERL database. We will use this property and assume that the projected BSDF is a radially symmetric function. us, we only need to track a 1D variance. Note that for grazing angles, GGX is notably skewed. However, since GGX lobes are parametrized by an energy, an incident direction and a roughness, we have no control over the higher order statistics of this distribution.

Equivalent Roughness. A BSDF lobe is not stable under convolution. us, the statistics specify an equivalent lobe that is an approximation since the shape of the spherical function is not preserved. In the following, we will display and use the equivalent roughness of a directional statistic using the mapping between roughness and variance and its inverse (see Fig. 3(b)).

ATOMIC OPERATORS ON STATISTICS

In this section, we describe statistical atomic operators approximating the e ect of speci c light transport interactions. ose operators will later be used to evaluate BSDF lobe statistics (Section 5) and instantiate an equivalent BSDF model (Section 6).

e di erent operators we use are summarized in Table 1. We will discuss them in the following order: re ection by a rough surface (Section 4.1), refraction by a rough surface (Section 4.2), transmission and sca ering in plane-parallel participating media (Section 4.3 and 4.4).

Reflection by a Rough Interface

A rough surface distributes the incoming light based on its Bidirectional Re ectance Distribution Function (BRDF). For a microfacet Table 1. We summarize the di erent atomic operators and how we approximate the outgoing energy e, mean µ, and variance σ given the incident energy e i , mean µ i , and variance σ i . For each statistic, we indicate with an exponent whether it has to be used for a reflected or transmi ed lobe (for example σ R for the variance of a reflected lobe). This notation will be reused in the adding-doubling method. In this table, we use the convention that incident and outgoing vectors leave the surface. This explains why both reflection and refraction negate the mean.

Rough Re ection

Rough Refraction Absorption Forward Scattering

energy e R = e i × FGD ∞ e T = e i × [1 -FGD ∞ ] e T = e i exp -σt h √ 1-|µ i | 2 e T =e i [ σs h √ 1-|µ i | 2 ] exp -σ t h √ 1-|µ i | 2 mean µ R = -µ i µ T = -η 12 µ i µ T = -µ i µ T = -µ i variance σ R = σ i + σ R 12 σ T = σ i η 12 + f (s × α 12 ) σ T = σ i σ T = σ i + σ
BRDF, it can be described by:

L o (ω o ) = ∫ Ω F ( ω i , h ) G(ω i , ω o ) D(h) 4 ω i , n ω o , n L i (ω i ) ω i , n dω i , (1)
where L o (resp. L i ) is the outgoing (resp. incoming) radiance, Ω is the upper hemisphere, h Energy. e amount of energy re ected by the surface is called the directional albedo. For microfacet models, it is directly de ned by the integrated Fresnel FGD:

= (ω i +ω o ) /| |ω i +ω o | | is the half-vector, ω o ,
FGD = ∫ Ω F ( ω i , h ) G(ω i , ω o ) D(h) 4 ω i , n ω o , n ω i , n dω o (2)
We will omit all parameters for FGD = FGD(ω i , α, η + iκ), where η + iκ is the complex index of refraction (IOR) at the interface. For isotropic BRDFs, the elevation is su cient to parametrize the incoming direction. We approximate the energy of the outgoing radiance by decoupling FGD from incident radiance:

∫ Ω L o (ω o )dω o FGD × ∫ Ω L i (ω i )dω i . (3)
Since no closed form exists for FGD, we precompute it in a 4D table parametrized by an elevation, a roughness, and a complex index of refraction.

Mean. Largarde and De Rousiers showed that the mean of the reected lobe is shi ed towards the normal [Lagarde and De Rousiers 2014, Section 3.1.4].

is is typically used to fetch preintegrated environment maps in video-game engines. Note that this shi is directly integrated into the FGD texture. ey showed that a good approximation of this shi is:

ω r |x = βn |x + (1 -β)ω r |x , (4)
where

ω |x = [ω.x, ω.
] is the vector in the tangential plane, and β depends on the shadowing-masking term (see [START_REF] Lagarde | Moving Frostbite to PBR[END_REF] for details). However, this shi is only important for rough con gurations and we neglect it during the estimation

f (α) α 1.1 1-α 1.1 a) 2-bounce equivalent roughness b) Linear space transformation
Fig. 4. We plot the equivalent roughness of two bounces of light on parallel layers with di erent micro-geometries using a normal incidence

|µ i | = 0.
The first layer is depicted with color, from red α = 0 to blue α = 1, and the second layer is specified by the abscissa. The equivalent roughness does not linearly depend on the input roughnesses. We found a transform (b) mapping roughness to a space of linear variance for reflection. We show our approximation in dashed.

of the BSDF. It can be applied later on, during shading, to correct misalignment of the preintegrated environment maps.

Variance. A rough surface increases the variance of the incident distribution [START_REF] Durand | A frequency analysis of light transport[END_REF]]. However, this increase is not linear with respect to roughness. It is approximately linear with respect to variance for small roughnesses (α < 0.2). For higher roughnesses, we can nd a space where rough re ections have a linear behavior (see Appendix A). at is, the variance of the light bouncing twice on parallel surfaces with roughnesses α 1 and α 2 is:

σ 12 = f (α 1 ) + f (α 2 ), (5)
where f (α) is the transform of roughness to variance. We display this transformation in Fig. 4(b). We convert roughness to linear variance using:

σ = f (a) ≈ a 1.1 (1 -a 1.1 ) , (6)
and from linear variance to roughness using the inverse. is transform does not perfectly t the linear-space mapping, but it approximates the outgoing directional statistics of two light bounces on

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 |µt| -1.0 -0.8 -0.6 -0.4 -0.2 0.0 |µt| η12 = 0.8 η12 = 1.2 η12 = 1.4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 |µi| -1.0 -0.8 -0.6 -0.4 -0.2 0.0 |µt| η12 = 0.8 η12 = 1.2 η12 = 1.4 a) α = 0.1 b) α = 0.6
Fig. 5. We display the mean of the transmi ed lobe µ t = sin(θ t ) with respect to the incident direction (specified by

|µ i | = | sin(θ i ) |) for di erent refractive index.
We also display (in dashed) the refracted direction sin(θ r ).

We can see very li le deviation of the mean transmi ed direction from the refracted direction for small roughnesses (a), but a noticeable one for higher roughnesses (b).

microfacet interfaces of moderate roughnesses (below 0.4) well, as shown in Fig. 4(a). An interactive demo of this approximation using a two bounce setup is available in our supplemental material 2 .

Energy Conservation. e classical microfacet theory only accounts for single sca ering in the microgeometry and thus omits the energy due to multiple bounces in the microgeometry. According to [START_REF] Heitz | Multiplesca ering microfacet BSDFs with the Smith model[END_REF], a scaling factor of the directional albedo approximates multiple sca ering well. us, it should have no effect on the mean and variance and only a ect the energy term. We incorporate it directly into the FGD term during precomputation. For a rough interface, we replace the previously de ned FGD term by FGD ∞ which is the result of the stochastic evaluation of Heitz et al. 's model. is new table has the same 4D parametrization.

Summary. Given an incident directional function with energy e i , mean µ i and variance σ i , the re ected energy e R , mean µ R , and variance σ R are approximated by:

e R = e i × FGD ∞ (7) µ R = -µ i (8) σ R = σ i + f (α) (9)

Refraction by a Rough Interface

A rough dielectric interface distributes light in both the upper and lower hemisphere of local directions. e later is described by the BTDF and, in the context of microfacet theory, its de nition is similar to the BRDF [START_REF] Walter | Microfacet Models for Refraction through Rough Surfaces[END_REF]]. In the following, we will note

η 12 = η 1
η 2 the ratio of refractive indices of the two media.

Energy. e energy scaling of an incident angular eld by a rough dielectric interface is given by 1 -FGD ∞ . is also ensures energy conservation between the re ected and refracted directional elds.

Mean. For rough interfaces, the transmi ed lobe peak is not the refracted incident direction. We found that, as with o -specular re ection, the mean refracted direction is roughness dependent 2 See webgl/validate reflection.html Fig. 6. The transmi ed lobe (in blue) has, for the same roughness, a di erent angular width compared to the reflected lobe (in red). We analyze them in the sca ering plane, parametrized by the sine of the outgoing angle for an incident angle of 30 degrees (c). We scaled the reflectance and transmi ance to one to compare the width of both lobes. To compute the variance we fake rough refraction using a rough reflection with a new incidence direction µ r = η 12 µ i (le , in green). The transmi ed roughness we use is α r = s × α , where α is the roughness of the microgeometry. We show that our fake refraction (le , in green) closely matches the reference curve (le , in blue), while the one of [START_REF] Kulla | Revisiting Physically Based Shading at Imageworks[END_REF] (le , in red) overestimates it.

(see Fig. 5). Its direction also depends on the refractive index. At grazing angles, and for higher refractive index ratio, the peak of the lobe di ers more from the refracted incident direction. However, we empirically found that using the purely refracted direction was su cient (see our supplemental material).

Variance. Due to Snell's law, light rays are bent when passing through a dielectric interface. Likewise, the incident variance is scaled by the ratio of the index of refraction when transmi ed. As in the re ection case, this scaled incident variance is increased by the roughness of the surface. We derive this property by considering transmission like a re ection. For the same roughness and incident direction, the refracted and transmi ed lobes have di erent widths (see Fig. 6). We use the notion of transmi ed roughness to analyze the variance of the transmi ed lobe. e idea is to fake the transmi ed lobe using a re ection from underneath the interface (see Fig. 7 (a) in green). To match the reference transmi ed lobe, we scale the roughness of our fake transmission. For stretch-invariant microfacet distributions, the scale s is the ratio of the derivatives of the transmi ed and re ected half vector's tangents:

s = 1 2 1 + η 12 ω i • n ω t • n . ( 10
)
We detail its derivation in Appendix B. [START_REF] Kulla | Revisiting Physically Based Shading at Imageworks[END_REF] derived a similar scaling factor by equating the fake re ected and transmi ed lobe peaks but it overestimates the transmi ed roughness. [START_REF] Guo | Rendering in Transparent Layers with Extended Normal Distribution Functions[END_REF] approximated a von Mises-Fisher exponent similar to our equivalent roughness. However, neither of them provide a close match to the reference BSDF. We compare our scaling factor to the one of [START_REF] Kulla | Revisiting Physically Based Shading at Imageworks[END_REF] in Fig. 7 (b). An interactive demo of the o -specular transmi ed direction and the di erent scaling factors is available in our supplemental material 3 .

Summary. Given an incident directional function with energy e i , mean µ i and variance σ i , the transmi ed energy e T , mean µ T , and variance σ T are approximated by:

e T = e i × (1 -FGD ∞ ) (11) µ T = -η 12 µ i (12) σ T = σ i η 12 + f (s × α) (13)

Volume Absorption

Assuming that the interface between the medium and the outside (be it air or another medium) is already resolved (as discussed in the previous sections), we can study the changes in energy, mean and variance due to a participating media. In this Section, we cover how absorption a ects statistics. e impact of sca ering on statistics will be covered in Section 4.4. As in Monte-Carlo rendering, we treat them separately and combine them together a erwards.

Energy. e a enuation of the incident light-eld by a random medium is described by Beer-Lambert's law. It states that the incident radiance is scaled down depending on the optical depth:

L o (ω o ) = exp -σ t h ω o , n L i (ω o ), (14)
where σ t is the transmi ance cross-section and h is the depth of the layer. We approximate the average energy using the incident mean a enuation:

e o = e i exp -σ t h ω o ,n .
Mean and Variance. Since the optical depth between the bo om and the top layers depends on the incident direction, the mean and variance are also a ected (see Fig. 8(a) for the variance). We empirically found out that the impact of absorption on the mean and variance is negligible. An interactive demo of this impact and of our approximation is available in our supplemental material 4 .

Summary. Given an incident directional function with energy e i , mean µ i and variance σ i , we approximate the energy e T , mean µ T , variance σ T transmi ed by volumetric absorption by:

e T = e i exp - σ t h ω o , n (15)
µ T = -µ i (16) σ R = σ i (17)

Volume Sca ering

In this section, we restrict ourselves to the study of an optically thin homogeneous slab of height h that does not emit light. In such a case, single sca ering is predominant and we can neglect multiple sca ering. Given the light incident to the slab L i (ω), the amount of light sca ered by the medium is [START_REF] Pharr | Physically Based Rendering, Second Edition: From eory To Implementation[END_REF]:

L o (ω) = σ s ∫ h 0 exp - σ t t ω, n L s (t, ω)dt, (18)
with

L s (t, ω) = ∫ S 2 p(ω i , ω) exp - σ t (h -t) ω i , n L i (ω i )dω i , (19)
where σ s is the sca ering cross-section.

Energy and Mean. Assuming that the phase function is strongly forward sca ering (most of its energy is sca ered forward), we approximate the outgoing energy using the a enuation evaluated in the mean direction of the incoming light:

∫ Ω L o (ω)dω ≈ σ s h exp - σ t h ω o , n ∫ Ω L i (ω)dω. ( 20
)
is approximation assumes that the phase function does not "lose" energy in the backward directions. However, it is possible to account for backsca ering by modulating this energy by the amount of light the phase function sca ers forward. For ≥ 0.7, this scaling is unnecessary. As for the absorption, we empirically found that sca ering did not alter the mean signi cantly for forward phase functions.

Variance. As for rough interfaces, the incident variance increases by the width of the phase function in the forward direction [START_REF] Belcour | A local frequency analysis of light sca ering and absorption[END_REF]. We empirically found a t of the equivalent GGX roughness for the Henyey-Greenstein phase function (HG) in the forward direction (see Fig. 8). It is given by:

σ = 1 -0.8 1 1 + , ( 21
)
where is HG's anisotropy factor. However, the Henyey-Greenstein phase function can have a non-negligible backsca ering when 0.7. is creates two modes in directional statistics that need to be resolved separately in our methodology. We focused on forward sca ering only media and le the derivation of the backsca ering variance to future work. An interactive demo of the impact of HG's factor on the angular statistics and of our di erent approximations is available in our supplemental material 5 .

Summary. Given an incident directional function with energy e i , mean µ i and variance σ i , the transmi ed energy e T , mean µ T , variance σ T from volumetric sca ering are approximated by:

e T = e i σ s h ω o , n exp - σ t h ω o , n (22)
µ T = -µ i (23) σ R = σ i + σ (24)

STATISTICS WITH MULTIPLE LAYERS

e outgoing light distribution aggregates many light paths (e.g. T RT andT RRRT paths). So far, our statistical analysis models groups of paths undergoing similar transport (e.g. all T RT paths). us, we need to combine the statistics of those path groups in an e cient way. For that, we rely on the adding-doubling algorithm. We rst apply it to the case of energy (Section 5.1).

en, we apply it to compute the variance (Section 5.2). Since the mean approximately aligns with the re ected or refracted direction, it is preserved a er multiple bounces and we do not treat it.

The Adding-Doubling Method

e Adding-Doubling method (see [START_REF] Grant | Discrete space theory of radiative transfer. I. Fundamentals[END_REF] for a complete overview) allows up to express radiative transfer in plane parallel media using a discrete form.

e idea is to model transmission and re ection due to a thin homogeneous slab with linear operators and to combine them to estimate the transmission and re ection of a thick heterogeneous slab. We recall here its mathematical principle in a simple form. For that, we will assume that the interacting interfaces are purely smooth.

Adding-doubling requires the de nition of the re ection r and transmission t coe cients of radiance with respect to the light's direction of propagation: r k k+1 is the re ection coe cient of the slab for light propagating downward and r k +1 k is the re ection coe cient of the slab for light propagating upward. With no light source in the interface, the light exiting an interface upward can be expressed as (see Fig. 9):

l + o = r 12 l - i + t 21 l + i , (25)
reads: light propagating upward from the slab is the incident light propagating downward re ected by the slab, plus the incident light propagating upward transmi ed by the slab. In the same fashion, we can express the light propagating downward from the slab as:

l - o = r 21 l + i + t 12 l - i . ( 26
)
5 See webgl/validate scattering.html Fig. 9. Given the transmi ance and reflectance r 12 , t 12 , r 21 , and t 21 of a layer (le , between two media, η 1 and η 2 ). The adding-doubling method permits to calculate the sca ering properties of a stack of layers (right), accounting for the inter-reflection between them (that is all T R + T paths).

We can express the global upward re ectance summarizing light transport between two layers r 12 and r 23 as:

r 13 = r 12 + t 12 r 23 t 21 + t 12 r 2 23 r 21 t 21 + • • • = r 12 + ∞ k =0 t 12 r k +1 23 r k 21 t 21 . ( 27
)
is form accounts for all bounces between the two layers. It can be summarized using the analytic form of this arithmetic series:

r 13 = r 12 + t 12 r 23 t 21 1 -r 23 r 21 . (28)
Similarly, we can express the global upward transmi ance as:

t 31 = t 32 t 21 1 -r 23 r 21 . ( 29
)
We can as well express the downward re ectance and transmi ance as:

r 31 =
e adding algorithm iteratively expresses the re ectance and transmi ance of an increasing number of interfaces by using r 13 , t 31 , r 31 , and t 31 as the re ectance and transmi ance of a virtual interface in Equations 28 to 31. is method can also evaluate multiple sca ering in a homogeneous medium of any depth with the doubling algorithm. Here, we stack together the same layer to generate a virtual layer of twice the depth. at is, we apply Equations 28 to 31 with r 12 = r 23 , t 12 = t 23 , r 21 = r 32 , and t 21 = t 32 . By iterating this operation on an input layer of very small size (typically h = 10 -8 ), we can increase the depth of a layer to any size.

Using adding-doubling, we can now evaluate the directional albedo of the layered BSDF. is is done by replacing the di erent r and t terms with the re ectance and transmi ance terms studied in Section 4. For a rough surface, we have:

r 12 = FGD ∞ , t 12 = [1 -FGD ∞ ], r 21 = FGD ∞ , and t 21 = [1 -FGD ∞ ].
Total Internal Re ection. While the de ned adding-doubling algorithm permits approximating multiple sca ering between two layers, it partially misses the impact of the totally re ected light by the upper Fresnel interfaces. In the case of a smooth interface, there are angular con gurations where no light is refracted and the dielectric interface behaves like a pure mirror. is behavior is called Total Internal Re ection (TIR). Transmission t 21 does not account for the angular spread of light re ected on interface 2 → 3 and we need to bring this information back:

L t (ω o ) = ∫ Ω L i (ω i )T (ω t ) ρ(ω i , ω t , α 23 ) dω t , (32)
where

ω t = refract(ω o , η 12 ), T (ω t ) = 1 -F (ω t )
is the Fresnel transmi ance, and ρ(ω i , ω t , α 23 ) is the microfacet BRDF of layer interface 2 → 3. When TIR occurs, T (ω i ) = 0 for ω i outside of the extinction cone. A solution to compute the amount of energy lost due to TIR is to decouple its computation from the transport integral:

L t (ω o ) = ∫ Ω L i (ω i ) ρ(ω i , ω t , α 23 ) dω t × ∫ Ω (1 -F (ω t )) D(ω t , α 23 ) dω i =TIR , (33)
where D(ω i , α 23 ) is the NDF. e second integral of this product is similar to the FGD table. We thus precompute this term in a 3D table.

We show in our supplemental material 6 that this decoupling gives good results in general. During the adding-doubling, we replace r 21 , and t 21 by:

r 21 ← r 21 + (1 -TIR) × t 21 (34) t 21 ← TIR × t 21 (35)

Adding-Doubling for Variance

We reuse the idea of adding-doubling for the case of variance. As seen in Sections 4.1 and 4.4, the variance of the interaction with an interface or medium has an a ne form. Also, the variance of a weighted sum of distributions sharing the same mean is the weighted arithmetic mean of the individual variances. ose properties permit the use of the adding-doubling methodology to compute variance due to multiple sca ering between layers.

Multiple Sca ering. Given two interfaces, the unnormalized average variance σ R 13 accounting for multiple sca ering is:

σ R 13 = r 12 σ R 12 + +∞ k=0 t 12 r k+1 23 r k 21 t 21 σ T 12 + (k + 1) s R 23 + k σ R 21 21 + σ T 21 , (36)
where s R i j (resp. s T i j ) is the additional variance when re ecting (resp. transmi ing) on a rough interface between indices of refraction i and j, and i j is the transmission scaling factor. is formula can be separated into a geometric series and an arithmetico-geometric 6 See webgl/compute TIR.html series. e geometric series has the analytic form we saw for the classical adding-doubling. e arithmetico-geometric series has the following analytic form [START_REF] Riley | Mathematical methods in physics and engineering[END_REF] Similar forms can be derived for the transmi ed average variance. We detail them in Appendix C. ey will be used in the next iteration of the adding algorithm in place of σ R 12 and σ T 12 . For the case of participating media, we use σ R 12 = 0, σ R 21 = 0, σ T 12 = σ , σ T 21 = σ , 12 = 1, and 21 = 1 (see Equation 21for the de nition of σ ).

]: +∞ k =0 kr k = r (1 -r ) 2 with r ∈ [0, 1[. ( 37
Adding-Doubling. Using the multiple sca ering equations, we can build an adding-doubling algorithm. We start with an empty layer, that is e R i j = e R ji = 0, e T i j = e T ji = 1, and

σ R i j = σ R ji = σ T i j = σ T ji = 0.
For each interface in the stack (starting from the upper one to the lower one) we apply the multiple sca ering equations on energies and variances (e.g. Equation 38and 28 for the upper re ected distribution). We also need to track the refraction scaling factors i j , and ji during the process.

Limitations. Note that those derivations assume that all lobes share approximately the same mean direction. Consequently, the described adding-doubling method cannot work with multimodal BSDFs without special-case handling. For example, to handle retrore ective phase functions (e.g. HG with < 0) we need to track another set of statistics in the retro-re ective direction.

IMPLEMENTATION

Based on our ndings, we developed two BSDF models: one tailored for forward path tracing and real-time rendering, denoted the Forward model (Section 6.1); and one tailored for bidirectional light transport algorithms, denoted the Symmetric model (Section 6.2).

Forward Model

As illustrated in Fig. 10, to evaluate the BSDF, we start from a unit energy e i = 1, the incident mean

µ i = [ ω i |x , ω i | ],
a zero variance σ i = 0, and evaluate the adding-doubling method to gather intermediate variance (See Appendix C) and energy in a vector of BRDF lobes7 . For the transmi ed lobe, we record the transmi ed statistics at the end of the adding-doubling algorithm8 (Equation 49). Each entry in the vector corresponds to a BSDF lobe and contains its energy e k , mean µ k , and variance σ k . From this vector of statistics, we instantiate an approximate BSDF model consisting of a weighted sum of microfacet GGX lobes with Fresnel e k , incident direction re ect(µ k ), and roughness f -1 (σ k ). We use our fake refraction

ω i n ω i n ω r ω t ω t n ω r ω t ω i n a)
Layered structure b) Adding: rst iteration c) Adding: last iteration d) BRDF instantiation Fig. 10. For shading, we track the BRDF lobe statistics starting from the incident directions (a). We estimate iteratively a vector of outgoing lobes using the adding-doubling algorithm (b-c). Finally, we instantiate a set of BRDF lobes using the vector of statistics and the inverse mapping. This process outputs as many lobes as there are layers. It is possible to merge statistics together to reduce the number of BRDF instances.

model for the transmi ed lobe. us the re ection and transmission models are:

ρ(ω i , ω o ) = N k =0 e k ρ k (ω k , ω o , α k ) , (39)
with

           α k = f -1 (σ k ) ω k = re ect(µ k ) ρ k (ω k , ω o , α k ) = D(h)G(ω k , ω o ) 4 ω k , n ω o , n
where e k , ω k and σ k are the energy, mean, and variance for the k th lobe respectively. For be er performance, it is possible to merge some of the variances (scaled by energies and then normalized by the total energy). In our implementation, we kept all lobes as-is.

Importance Sampling & MIS. Since our models are a weighted average of multiple lobes, we randomly select one based on the energy and importance sample the visible normals [START_REF] Heitz | Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals[END_REF] based on the fake incident direction and roughness. However, this strategy is not optimal and creates re ies since the di erent lobes overlap. To avoid this we use multiple importance sampling on the di erent lobes with the balance heuristic [START_REF] Veach | Robust monte carlo methods for light transport simulation[END_REF]]. e contribution of a BRDF lobe sample becomes9 :

p = e all N i=0 e i p i N i=0 e i ρ i (ω i , ω o , α i ), (40)
where ω o is the outgoing direction selected by one of the strategies, e i is the energy for the i t h lobe, e all is the total energy, p i is the probability density function (pdf) of sampling the vNDF for roughness α i , and ρ i is the pdf of the i t h microfacet model.

Real-Time Model. For real-time scenarios, we reuse the Forward model, restricted to two or three lobes. Since, in all real-time engines, lights such as area lights and environment maps are preintegrated with respect to a GGX model, we evaluate our model in this context.

Symmetric Model

For o ine scenarios where symmetry of light transport ma ers (such as BDPT or MLT), we build an ad-hoc symmetrization from the Forward model. As [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF], we assume that the incoming and outgoing directions share the same microfacet.

We further use the property that the refracted rays will share the same half-vector (see Appendix E for proof). It is thus possible to use our method with respect to the di erence vector θ d . is virtually injects the adding-doubling roughness on the half-vector normal.

e resulting statistics might not match the Forward model, since the FGD term is evaluated using cos(θ d ) and not cos(θ i ) and the transmi ed roughnesses will slightly di er. is doesn't however introduce an over-blur of the resulting lobe, and we found that the results of the Symmetric model were o en closer to the reference (see our supplemental material).

e Symmetric model uses the same importance sampling strategy as the Forward model.

RESULTS

In this section, we demonstrate the use of our di erent models in both o ine (Section 7.1) and real-time (Section 7.2) contexts. All results are computed on a 16 core i7 processor with 32GB of RAM and an Nvidia GTX 980 graphics card. All computations are using parallel computing as much as possible. For our o ine method, we precomputed the FGD using the complex index of refraction, resulting in a 4D table, precomputed at a resolution of 64 4 (making it 64MB in size). We used linear interpolation to evaluate it. For the real-time method, we used the available FGD 2D texture that performs a split sum of Schlick's Fresnel [START_REF] Karis | Real Shading in Unreal Engine 4[END_REF]]. We also precomputed the 3D TIR table using a 64 3 resolution, making it 1MB in size. Unless noted, we use the Symmetric model for our o ine results.

O line Rendering

Validation. We validated our method against a stochastic ground truth. See Appendix D and our supplemental material for more details. We show that our model captures the color saturation due to multiple sca ering (Fig. 11) and is energy conserving (Fig. 12).

anks to the adding method, we do not lose energy due to multiple sca ering between layers. Even though the re ectance and transmi ance terms are approximate, our model is close to the saturated color of the reference. We also show in the F M scene (Fig. 13) that we correctly predict the rough appearance of frosted paints. Here, we replicate a gure of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF]: a rough dielectric interface on top of a smooth conducting one. e method of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] Fig. 12. Our model is energy conserving (as long as the FGD texture is). In this scene, we show that for a complex layered structure with five dielectric layers (alternating η = 1.3 and η = 1) on top of a pure mirror, our model does not lose energy as opposed to the one of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF].

roughness between layers: the base layer still behaving as a mirror on the global re ectance. Note that [START_REF] Elek | Layered Materials in Real-Time Rendering[END_REF] mitigates this issue by using the maximum of roughness of the top and bo om layers. Please refer to our supplemental document for more comparisons with ground truth and the model of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] In the D scene (Fig. 14), we replicated a gure from Jakob et al. [2014] to show how our method handles participating media.

e dragon gure consists of a golden base on top of which we apply a di using medium with increasing optical depth. We used a forward phase function = 0.7, and used σ a = 0 and σ s = 1.

anks to our statistical analysis, we easily predict the di usion and a enuation of light rays inside this structure. However, to simplify our implementation, we only account for single sca ering within media layers (see Section 7.3).

In the P (Fig. 15) and R (Fig. 16) scenes, we display textured assets. In P , we use a single texture for the roughness of some of the layers and in R we combined three textures together. Compared to Jakob et al. [2014], our method is not limited by precomputation of discrete values of the parameters and resolves a BSDF on the y. Consequently, our method handles e ciently any textured combination of inputs for both surface and media layers. For example, in the R scene, three parameters are textured. To illustrate, the method of Jakob et al. [2014] requires 1.7GB of storage for the M P example alone and its computation takes 21min and requires up to 15.5GB of RAM. Texturing this asset would require an enormous amount of disk storage for the di erent BSDF les and they might not even t in memory for rendering. Fig. 13. We reproduce the F M example from [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF]. While [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] misses the correlation between layers, we stay close to the reference image computed using a stochastic evaluation of the layered structure.

No media d = 0.1 d = 0.2 Fig. 14. The D scene illustrates the ability of our method to add absorption and single sca ering. Here, we start from a specular gold dragon and add a layer of sca ering particles on top (with = 0.7). As we increase the optical depth of this layer, the dragon appears rougher and darker due to sca ering and absorption.

Symmetric model Forward model Symmetric model

Two layers with conducting base ree layers with absorption Fig. 15. The P scene illustrates our ability to work with textured assets. In the first example, the plate is composed of a two layers material with η 0 = 2, α 0 = 0 for the first layer and η 1 = 1 + 0.5i, α 1 ∈ [0, 1] for the second layer. We display side-by-side our symmetric and forward model. In the second example, the plate is composed of three layers with η 0 = 1.8, α 0 ∈ [0, 1], d 1 = 1, σ a = [1, 0.7, 0.2] , and η 2 = 1, α 0 =∈ [0, 1]. We used a specialized plugin for dielectrics since Mitsuba separates the two cases.

Peformance. We evaluated the performance of our code base. We report the timings of our o ine method compared to the method of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] in Table 2. Our method's costs are primarily due to the adding-doubling calculation and to the linear interpolation of tables.

e la er cost is alleviated on graphics hardware.

Real-Time Rendering

We implemented our real-time BRDF model inside a commercial 3D engine using forward rendering to validate the practicability of our implementation. Note that our model applies to deferred rendering as well. In our implementation, we restricted the number

η 0 η 0 κ 1 κ 1 α 0 α 0 η = η 0 , α 0 η = η 0 , α 0 η = 1 + iκ 1 , α 1 η = 1 + iκ 1 , α 1
Fig. 16. The R scene illustrates the ability of our method to e iciently work with textured assets. In this scene, we used a 2-layer structure where η 0 , α 0 , and κ 1 are texture mapped. We set η 1 = 1 and α 1 = 0. of layers to three, and outgoing lobes to two. We also xed the second interface to be a participating medium. We used the same adding code for the Mitsuba implementation and for the real-time rendering one. We show a comparison between our real-time implementation, our Forward model, and the reference in Fig. 17. Our real-time implementation matches the reference appearance albeit with the technical limitations of using a simple Fresnel form, and environment map pre-integration. We measure the running time of our shader using a full-screen quad at a resolution of 1920 × 1080 pixels. Here, the full frame takes between 1.9ms to 2.1ms to render. We compare this running time to the standard shader of the engine. In this case, the full frame takes between 1.7ms to 2.0ms to render. e B scene (Fig. 18) demonstrates our real-time implementation. It depicts a car fully covered with our layered BRDF model. We rendered it at a resolution of 1920×1080 during a material design session that we show in our supplemental video.

Limitations & Failure Cases

Accuracy with respect to roughness. Our model approximates multiple sca ering lobes using equivalent GGX models. While this faithfully predicts the behavior of layered materials for relatively smooth surfaces, it becomes increasingly inaccurate when most of the layers are very rough. We illustrate the lost of accuracy of our model in Fig. 19. Here we can see that for multiple layers with Table 2. We report the cost of our method compared to the one of [START_REF] Weidlich | Arbitrarily Layered Micro-facet Surfaces[END_REF] for three scene: F M (Fig. 13, at 1024 spp), W F (Fig. 12, at 1024 spp), and D (Fig. 14, at 512 spp). F M and W F use the direct integrator while D uses the path integrator. We report the timings of our method without linear interpolation of the FGD and TIR texture in parenthesis. roughnesses close to one, our model signi cantly diverges from the ground truth.

Color shi . Because of total internal re ection, the color of the multiple sca ering between layers might not be the same as the transmi ed color. is happens for metals with a strong color shi at grazing angles. In such a case, our method will not be able to reproduce the exact tint of the material. is is noticeable in the F M example (Fig. 22) where at grazing angles the base layer takes on a green tint that our model misses.

Multiple sca ering in media. Our prototypes are currently limited to single sca ering from participating media. is is a design choice to keep the implementation simple to read. Accounting for multiple sca ering would require to add another direction in the addingdoubling algorithm, i.e. the retro-re ective direction. is would allow us to track sca ering between media interfaces and to add the missing retro-re ective look of dust.

Interfaces. Due to the derivation of the linear space for GGX, our model only accounts for a single type of NDF for interfaces. We cannot accurately predict the case of multiple interfaces with mixed NDF geometries (such as Beckmann [START_REF] Walter | Microfacet Models for Refraction through Rough Surfaces[END_REF], GTR [START_REF] Burley | Physically-based shading at disney[END_REF]], Student-t [Ribardière et al. 2017], etc.). Note, however, that a linear space can be derived for any NDF as long as all interfaces are described with such an NDF. Similarly, we cannot accurately reproduce the e ect of a Lambertian layer. is is due to α = 0.9 Fig. 19. For higher roughnesses, our model lose its accuracy compared to a stochastic ground truth. In this example, we use a rough dielectric layer, η = 2, on top of a rough conducting base, η = 0.01 + i, with equal roughnesses α 1 = α 2 = α . Here, we did not account for multiple sca ering in the microfacet geometry and isolate the error due to approximating multiple sca ering between layers.

the fact that the resulting distribution will always be a GGX lobe. Another strong limitation that is shared by all the works on layered materials is that we cannot work with non-parallel interfaces. While using a di erent normal to evaluate Fresnel for the base and coat layers is typically done in production and works with our method, it is incorrect with respect to the resolution of multiple sca ering.

Skewness at grazing angles. Due to the extinction cuto of the Fresnel transmission, transmi ed lobes at grazing angles are anisotropic and skewed. Since we used radially symmetric NDFs for our analysis, we cannot model this behavior correctly.

DISCUSSION & FUTURE WORK

Complex Appearance. While we restricted ourselves to geometrical optics for our study by enforcing the thickness of layers to be greater than the visible light wavelengths, it is possible to incorporate iridescence e ects into our model. Since we rely on the evaluation of the integrated Fresnel term during the adding-doubling, we can use any model that replaces the Fresnel term. For example, it is possible to replace one of the evaluations of the FGD term by the model of [START_REF] Belcour | A Practical Extension to Microfacet eory for the Modeling of Varying Iridescence[END_REF]. Similarly, we could incorporate gli ery appearance into our framework by using one of the available glint models [START_REF] Wenzel | A comprehensive framework for rendering layered materials[END_REF]. Note that in such a case, the correlation between the discrete evaluation of glints and the upper layers' Fresnel will be missed.

Anisotropy & Spatial Di usion. While we restricted our analysis to isotropic GGX BSDFs, our analysis could be applied to track anisotropic lobes. Here, the tracking of variance in the addingdoubling algorithm would require tracking eight scalars instead of four. It would be interesting to see if the statistical analysis can be performed on the spatial components as well, in order to track a subsurface sca ering pro le.

Other Applications. Since our model is fast to evaluate and doesn't require a per-material pre-computation, it could be used to perform nonlinear ing of data. Also, it would be interesting to see if our layered materials model can be used for inverse design. ere, artists could design the desired look using many lobes and extract a set of layers with speci ed index of refraction and roughnesses. is could be interesting for the design and manufacture of real life objects.

CONCLUSION

We introduced a novel statistical analysis of layered materials that builds on an atomic decomposition of light transport to track the energy, mean and variance of the layered BSDF. Our analysis is versatile and can account for both surface and volumetric light-ma er interaction. Furthermore, we leveraged the power of the addingdoubling method on this statistical representation to compute the layered BSDF accounting for multiple sca ering in the layered structure. We demonstrated the applicability and accuracy of our representation with o ine and real-time implementations showing an e cient, energy conserving and symmetric BSDF model for layered materials. Our model still requires precomputing a 4D and a 3D table. But, given the regularity of this data, it seems possible to nd good and e cient approximations that would make our model much faster. Although expanding our model to account for anisotropy or spatial di usion is still open, we believe that our method is a rst step towards e cient and lightweight BSDF/BSSRDF models for layered materials.

A FINDING GGX'S LINEAR SPACE

To nd the function f (α) that transform the roughness of the GGX lobe into a space where multiple bounces behave linearly, we start from the 2-bounces case where both surface have the same roughness α:

α 12 =f -1 [f (α) + f (α)] (41) =f -1 [2f (α)] , (42)
It follows that for this con guration:

f (α 12 ) = 2f (α). (43)
Using our numerical evaluation of the variance w/r roughness for the 1-bounce and 2-bounces case, we can extract α 12 as a function of α. For that, we invert numerically the function variance w/r roughness for the 1-bounce case and evaluate the roughness corresponding to the 2-bounces variance. Fig. 20, illustrate this step.

Using this knowledge, we can extract the linear space transformation (see Fig. 4 le inset) up to a scaling factor. We iterate Equation 43 from a small value for both α and f (α) (we used 10 -3 ) and iteratively evaluate α 12 from the graph and f (α 12 ) using the current function value.

B ROUGHNESS SCALING FOR TRANSMISSION

Here we derive the following Jacobians:

t = d tan(h t ) dθ t = η 12 ω t • n η 12 ω i • n + ω t • n (44) r = d tan(h r ) dθ o = 1 2 (45)
Remember that:

tan(h t ) = η 1 sin(θ i ) + η 2 sin(θ t ) η 1 cos(θ i ) + η 2 cos(θ t ) .
It follows that:

t = d tan(h t ) dθ t = η 2 cos(θ t ) η 1 cos(θ i ) + η 2 cos(θ t ) + η 2 sin(θ t ) [η 1 sin(θ i ) + η 2 sin(θ t )] [η 1 cos(θ i ) + η 2 cos(θ t )] 2 = η 12 [cos(θ i ) cos(θ t ) + sin(θ i ) sin(θ t )] + 1 [η 12 cos(θ i ) + cos(θ t )] 2
and using the equality sin(θ i ) sin(θ t ) = -1 η 12 sin 2 (θ t ), we obtain: Assuming that we are in the fake re ection con guration, we equate θ o = θ t and obtain: Fig. 21. We validated our adding method (in blue) with respect to a Monte-Carlo simulation of the layered structure (in red). We can predict quite closely the mean variance curve of a 2-layers structure when varying one of the roughness. Due the limitation of our fi ing of the linear space transformation, our method is unable to accurately predict the correct mean roughness for higher input roughnesses. In this figure, we used η 0 = 1.4, α 0 ∈ [0, 1], η 1 = 1 + i, α 1 = 0.3 for the two layers configuration and η 0 = 1.4, α 0 ∈ [0, 1], η 1 = 1.2, α 1 = 0.1, η 2 = 1.4, α 2 = 0.1, et a 3 = 1 + i, α 3 = 0.01 for the four layers configuration. Both configurations use ω i = [0, 0, 1].

C ADDING-DOUBLING

We express here the transmi ed average variance in its unnormalized form σT Note that when computing those intermediate variances, we will have to keep track of the global scaling factor 1i and i1 due to all interfaces from the top layer to the current bo om layer.

We validated our computation of the mean variance for various layer con guration (see Fig. 21. We found that our approximation of the mean variance is quite close to the ground truth for interfaces with low and moderate roughnesses (from 0 to 0.5). For higher roughnesses, we experience higher deviation from the ground truth.

is is primarily due to our approximation of GGX's linear space transformation, we only ed well the low to moderate part of the curve.

D COMPARING OFFLINE RESULTS

In Fig. 22, we provide the comparison of our method to a ground truth. To generate this later model, we stochastically evaluated the layer material by explicitly tracing the rays in the structure.

E HALF VECTOR OF TRANSMITTED RAYS

We derive the half vector of two rays wi and ω o refracted by a planar interface of normal the half vector h = ω i +ω o |ω i +ω o | . e expression of both refracted rays is [START_REF] Walter | Microfacet Models for Refraction through Rough Surfaces[END_REF]]:

i t = η cos(θ d ) -1 + η (cos(θ d ) -1) h -ηω i ( 52
)
o t = η cos(θ d ) -1 + η (cos(θ d ) -1) h -ηω o , ( 53
)
where η is the ratio of index of refraction, and θ d is the di erence angle in the half-vector parametrization. If we express the unnormalized half vector of i t , and o t , we obtain:

i t + o t = 2 η cos(θ d ) -1 + η (cos(θ d ) -1) h -η [ω i + ω o ] (54)
= 2 η cos(θ d ) -1 + η (cos(θ d ) -1) -η h.

e normalized half vector is thus h. Fig. 22. We produce various appearance using two layers (one dielectric and one metallic). We show that our method faithfully reproduces the apparent roughness of those materials and, to a certain degree, the color saturation due to multiple sca ering. For each example, we report the RMSE ∆ between our symmetric model and the reference.

Fig. 3 .

 3 Fig.3. Our statistical analysis is done on the projective plane of directions [u, ] = ω x , ω (le inset). We study the mean and variance of the outgoing distribution in this space. We show in the right inset the variance of a GGX lobe at normal incidence with respect to roughness (in blue) as well as the variance of light bouncing twice on parallel surfaces with the same roughness (in red).

  n denotes the dot product, and F , G, and D describe the BRDF: F is the Fresnel term, D the Normal Distribution Function (NDF) and G the associated shadowing/masking.

  ance or re ectance )

  Fig. 7. To compute the variance we fake rough refraction using a rough reflection with a new incidence direction µ r = η 12 µ i (le , in green). The transmi ed roughness we use is α r = s × α , where α is the roughness of the microgeometry. We show that our fake refraction (le , in green) closely matches the reference curve (le , in blue), while the one of[START_REF] Kulla | Revisiting Physically Based Shading at Imageworks[END_REF] (le , in red) overestimates it.

Fig. 8 .

 8 Fig.8. We display the equivalent roughness of both absorption (a) and sca ering (b) of a participating layer with the Henyey-Greenstein (HG) phase function. For the absorption, we used varying optical depths (from h = 0 in red to h = 1 in blue). For both plots, we vary the anisotropy parameter ∈ [0, 1]. Our approximation of the HG roughness (in red) closely matches the reference curve (in blue).

  Fig. 11. Thanks to the adding algorithm, we estimate the color saturation due to multiple sca ering between interfaces. Here we show our method (le ) compared to two references: one computing all bounces between interfaces (middle) and one computing only one reflection per layer (right).

Fig. 17 .

 17 Fig. 13 Fig. 12 Fig. 14 Ours 46s (24s) 1.8m (35s) 1.84m (1.64m) Weidlich [2007] 17s 20s 1.60m

Fig. 20 .

 20 Fig. 20. Using the mapping between roughness and variance for the 1bounce and 2-bounces cases (where both roughnesses are equal (a)), we can extract the apparent roughness of the 2-bounces configuration (b).

  49, α = 0.01 η = 1 + i[1 0.1 0.1], α = 0.01 η = 1.49, α = 0.1 η = 1 + i[1 0.1 0.1], α = 0.01 η = 1.49, α = 0.01 η = 1 + i[1 0.1 0.1], α = 0.1 η = 1.49, α = 0.1 η = 1 + i[1 0.1 0.1], α = 0.1 η = 1.49, α = 0.1 η = [0.143 0.373 1.444] +i[3.983 2.387 1.602] α = 0.1

  =cos(θ t ) η 12 cos(θ i ) + cos(θ t ) which is Equation 44. For the second Jacobian, we use the de nition of the re ected tangent:tan(h r ) = sin(θ t ) + sin(θ o ) cos(θ t ) + cos(θ o ) ) [cos(θ t ) + cos(θ o )] + sin(θ o ) [sin(θ t ) + sin(θ o )] [cos(θ t ) + cos(θ o )] 2.

		,
	we derive:
	r =	d tan(h r ) dθ o
	=	cos(θ o

t = η 12 [cos(θ i ) cos(θ t )] + 1sin 2 (θ t ) [η 12 cos(θ i ) + cos(θ t )] 2

  13 : σT 13 = t 12 t 23 23 σ T 12 + s T 23 + t 12 r 23 r 21 t 23 23 σ T normalized average variance, that will later be used in the remaining adding-doubling algorithm, is the unnormalized variance divided by the energy of light transmi ed through both interfaces. It can be expressed as: contains two elements. First the variance of the purely transmi ed lobe 23 σ T 12 +s T 23 (e.g. the TT paths) to which adds an additional variance due to multiple sca ering 23 s R 23 + σ R 21 (e.g. the T (RR) + T paths). Similarly, we nd the following expression for σ R 31 and σ T 31 :

					σ R 31 = r 32 s R 32 +	t 32 t 23 r 21 1 -r 21 r 23
					× s T 23 + 23 s T 32 + σ R 21 + s R 23 + σ R 21	r 21 r 23 1 -r 21 r 23	. (51)
				12 + s R 23 + s R 21 + s T 23 + • • •	(46)
	which simpli es to:
	σT 13 =	t 12 t 23 1 -r 23 r 21	23 σ T 12 + s T 23
			+ 23 s R 23 + σ R 21	r 23 r 21 1 -r 23 r 21	.	(47)
		σ T 13 = σT 13	1 -r 23 r 21 t 12 t 23	.	(48)
	We can further simplify this later equation to get:
	σ T 13 = 23 σ T 12 + s T 23 + 23 s R 23 + σ R 21	r 23 r 21 1 -r 23 r 21	.	(49)
	is later equation σ T 31 = 21 s T 32 + σ T 21 + 21 s R 23 + σ R 21	r 23 r 21 1 -r 23 r 21	(50)

e

ACM Transactions on Graphics, Vol. 37, No. 4, Article 73. Publication date: August 2018.

see our code mitsuba/plugins/layered forward.cpp, line 105

see our code mitsuba/plugins/layered dielectric.cpp, line 258 ACM Transactions on Graphics, Vol. 37, No. 4, Article 73. Publication date: August 2018.

see our code mitsuba/plugins/layered forward.cpp, line 453

ACKNOWLEDGEMENTS

e author thank Pascal Barla, Emmanuel Turquin, Simon Premože, and Romain Pacanowski for their comments during the early stage of writing, Kenneth Vanhoey for the video voice-over, and Stephen Hill and Derek Nowrouzezahrai for their English proof-reading. e R , T and B models were designed by MasterX-eon1001, Optic, and LSFMT from www.blendswap.com respectively. e F model is designed by Slavik from www.thingiverse.com.