Keywords: microfacettes, rendu, photoréalisme, anticrénelage, filtrage, MIP map microfacet theory, rendering, photorealism, antialiasing, filtering, MIP map

La synthèse d'images dites photoréalistes nécessite d'évaluer numériquement la manière dont la lumière et la matière interagissent physiquement, ce qui, malgré la puissance de calcul impressionnante dont nous bénéficions aujourd'hui et qui ne cesse d'augmenter, est encore bien loin de devenir une tâche triviale pour nos ordinateurs. Ceci est dû en majeure partie à la manière dont nous représentons les objets : afin de reproduire les interactions subtiles qui mènent à la perception du détail, il est nécessaire de modéliser des quantités phénoménales de géométries. Au moment du rendu, cette complexité conduit inexorablement à de lourdes requêtes d'entrées-sorties, qui, couplées à des évaluations d'opérateurs de filtrage complexes, rendent les temps de calcul nécessaires à produire des images sans défaut totalement déraisonnables. Afin de pallier ces limitations sous les contraintes actuelles, il est nécessaire de dériver une représentation multiéchelle de la matière.

Dans cette thèse, nous construisons une telle représentation pour la matière dont l'interface correspond à une surface perturbée, une configuration qui se construit généralement via des cartes d'élévations en infographie.

Nous dérivons notre représentation dans le contexte de la théorie des microfacettes (conçue à l'origine pour modéliser la réflectance de surfaces rugueuses), que nous présentons d'abord, puis augmentons en deux temps.

Dans un premier temps, nous rendons la théorie applicable à travers plusieurs échelles d'observation en la généralisant aux statistiques de microfacettes décentrées. Dans l'autre, nous dérivons une procédure d'inversion capable de reconstruire les statistiques de microfacettes à partir de réponses de réflexion d'un matériau arbitraire dans les configurations de rétroréflexion. Nous montrons comment cette théorie augmentée peut être exploitée afin de dériver un opérateur général et efficace de rééchantillonnage approximatif de cartes d'élévations qui (a) préserve l'anisotropie du transport de la lumière pour n'importe quelle résolution, (b) peut être appliqué en amont du rendu et stocké dans des MIP maps afin de diminuer drastiquement le nombre de requêtes d'entréessorties, et (c) simplifie de manière considérable les opérations de filtrage par pixel, le tout conduisant à des temps de rendu plus courts. Afin de valider et démontrer l'efficacité de notre opérateur, nous synthétisons des images photoréalistes anticrenelées et les comparons à des images de référence. De plus, nous fournissons une implantation C++ complète tout au long de la dissertation afin de faciliter la reproduction des résultats obtenus.

Nous concluons avec une discussion portant sur les limitations de notre approche, ainsi que sur les verrous restant à lever afin de dériver une représentation multiéchelle de la matière encore plus générale.

opportunity to experience this intense adventure in the first place. Thank you Jean-Claude for introducing me to the world of computer graphics, to OpenGL and the challenging problem of level-of-detail during my bachelor degree. Thank you Pierre for spending so much time on the deadlines with me, and for always being available for discussions during my stay at LIGUM; this year and a half spent in Montreal was unforgettable. Thank you Victor for sharing with me your fascinating stories on tilings and quasicrystals, and for always introducing me to the people you would invite over at LIRIS. My next thoughts go to my colleague and friend Eric Heitz, who has been by far my most important collaborator since the beginning on my PhD. To you, Eric, I must say that it has been both a pleasure and privilege to work with you; it really feels to me that no problem is too hard when we are working together. Next, I wish to thank several colleagues I have had the chance to meet during my PhD and who also had a great influence on my work. Stephen Hill: for giving me the opportunity to be a Most objects we see in everyday life can be described in two steps: first by shape, then by appearance; Figure 1.1 illustrates their differences. Computer graphics applications tend to follow this approach, relying on polygon meshes for the former, and physically based shading for the latter; Figure 1.2 shows a few such combinations.

When high-frequency detail is desired, these representations are further augmented with texture maps that either affect shape (e.g., displacement maps) or appearance (e.g., normal maps); Figure 1.3 shows a few common texture mapping techniques.

Despite being a remarkably effective solution for detail management, this combination of descriptors has led to a challenging problem for physically based rendering: we now have to somehow process thousands (or even more) of polygons and/or texels per pixel just to properly capture visibility and direct lighting effects. Since there usually are no analytic solutions to this problem, it must be solved numerically with, inexorably, a large number of samples (recall that the sampling rate should be at least twice the maximum frequency component of the scene according to the Nyquist sampling theorem). In most cases unfortunately, the computational cost induced by the Nyquist sampling theorem is too high to be practical. As a consequence, the integration ends up being undersampled for the sake of reasonable rendering times, ultimately leading in final images to unacceptable artifacts such as aliasing or noise.

In order to deal with this specific issue, computer graphics applications commonly rely on more or less sophisticated screen-space filters, ranging from fully automatic (e.g., SMAA [START_REF] Jimenez | SMAA: Enhanced Morphological Antialiasing[END_REF]) to fully manual (e.g., an artist working overtime in Photoshop). As can be evidenced by the quality of current photorealistic renderings, this approach works reasonably well in practice. However, it does tend to defeat some the main benefits of physically based rendering. In particular, screen-space filters do not guarantee shading consistency across pixels and time, thus requiring extensive tweaking to achieve satisfying results. Such limitations can quickly bottleneck a rendering pipeline, so it makes sense to look for alternative solutions.

Given the current evolution of computational power and the innovations in higher resolution / higher framerate displays, it seems that brute-force supersampling will remain out of reach for the foreseeable future (note: it will also be interesting to see how screen-space filters evolve, as their performance is intimately linked to display characteristics). In the meantime, rather than trying to adapt sampling rates to geometric complexity, there are potentially great benefits to be had by doing the opposite, i.e., adapting geometry to fixed amounts of samples;

this can be achieved through the use of precomputed level-of-detail schemes.

Finding generic and lightweight object representations that can maintain physically based principles at all scales is quite difficult, however. The difficulty mainly resides in the mathematics involved, which are mostly nonlinear, highly dimensional, and sensitive to discretization. The main goal of this thesis is to motivate the conjecture that, despite difficulty, viable solutions for precomputed level-of-detail schemes can be derived. We demonstrate this by introducing an effective prefiltering scheme for displacement mapped surfaces of arbitrary appearance, which works similarly to MIP mapping [START_REF] Williams | Pyramidal Parametrics[END_REF]. Our prefiltering scheme is derived within an augmented microfacet theory framework. Originally, microfacet theory is a mathematical framework designed for describing appearance through a bidirectional reflectance distribution function (BRDF) in a biscale approach;

microfacet BRDFs are the object of Chapter 2. In order to make this framework applicable to displacement mapped surface filtering, we extend it in two ways. First, we generalize the theory to multiple scales; this is the purpose of Chapter 3. Second, we introduce an inverted theory, i.e., a framework to recover a microfacet configuration given an input material; this is the purpose of Chapter 4. Our extensions of microfacet theory can then be used to derive a multiscale representation for displaced surfaces that is general, yet simple, robust, and straightforward to implement in a physically based rendering pipeline. To emphasize this last point, a complete state-of-the-art C++ implementation of a multiscale microfacet BRDF library is presented progressively throughout the manuscript, along with numeric validation tests to help the interested reader in better understanding and reproducing the results of this thesis.

Thesis Overview

More formally, this thesis is structured as follows (note: each chapter starts with a visual summary of its contents).

• In the remainder of this chapter, we make a more rigorous mathematical treatment of the problem of filtering displacement mapped surfaces of arbitrary reflectance. We take advantage of this effort to recall some of the fundamental rules that apply to physically based rendering, to introduce mathematical notations, and to describe a general approach that we employ to solve numerical problems.

• In Chapter 2, we review previous and concurrent work on BRDF modelling with microfacet theory. While microfacet BRDFs can be found in the microfacet literature since the 1960s, some of its most powerful properties, e.g., importance sampling with the distribution of visible normals and roughness invariance, were introduced concurrently to the development of this thesis.

• In Chapter 3, we introduce noncentral microfacet theory, a generalization of microfacet theory to derive the BRDF of a rough specular microsurface at any scale. This contribution allows us to transpose the difficult problem of displacement map filtering into that of evaluating a microfacet BRDF. By approximating the microsurface with a Gaussian random process, we derive an analytic displacement filtering scheme that can be precomputed into MIP maps. We refer to this scheme as LEADR mapping.

• In Chapter 4, we introduce inverted microfacet theory, a framework to extract state-of-the-art microfacet BRDF components from arbitrary materials. This contribution allows us to convert any material into a microfacet configuration, which may then benefit from useful microfacet BRDF properties, e.g., roughness authoring and fast importance sampling. We also show how to approximate the extracted microfacets into Gaussian or GGX random processes, which are very popular in the industry.

• In Chapter 5, we show how to combine the contributions of the thesis to derive a general multiscale representation for displaced surfaces that is lightweight, efficient, and authorable. We discuss benefits and limitations of this representation with respect to previous work, and identify future extensions required to derive a more general multiscale representation for matter.

(1) 3D coordinate system.

(2) Camera model (located at point p).

(3) The Rendering Equation [START_REF] Kajiya | The Rendering Equation[END_REF].

Theoretical Background

Camera Model A photorealistic renderer is concerned with the problem of generating a purely synthetic image that is indistinguishable from a real photograph. In the simplest cases, this problem is approached by discretizing a continuous distribution of incident radiance L > 0 into a set of pixels that constitute the final image, thereby mimicking the behavior of a pinhole camera sensor with instantaneous exposure. We restrict our mathematical treatment of light transport to geometric optics, i.e., light propagates as straight rays that can be reflected or absorbed by geometry throughout the 3D scene; this is a common assumption in physically based rendering. We can determine the radiance incident to any 3D point p ∈ R 3 by measuring the amount of light rays that pass through this point. We can visualize such light rays by mapping the radiance incident to p on a unit sphere centered on p; this is illustrated in Figure 1. 4 (2). Any portion of this sphere can be used to define the image plane of a virtual pinhole camera.

Photorealistic Pixels Since we wish to generate a discrete image, the frequencies exhibited by the incident radiance distribution function L that exceed the Nyquist limit should be removed in order to avoid artifacts. This is why the discretization process of L is typically done after low-pass filtering. Mathematically, letting W ∈ R denote the low-pass filter (e.g., box or Gaussian) weight associated to a given screen pixel, the intensity

I > 0 of the pixel is an inner product [Hec89, BN12] I = ΔW L o •W (o) • dω o , (1.1)
where, letting Ω denote the unit 3D sphere, ΔW ⊆ Ω is the filter's angular support on the sphere, o is the unit direction associated to the infinitesimal solid angle dω o , and L o ≥ 0 denotes the value of L at direction o, i.e., L o = L(o). Note that throughout this thesis, we manipulate unit directions either from their Cartesian components, or their spherical coordinates, which we link according to the following convention

o = ⎡ ⎢ ⎢ ⎣ sin θ o • cos φ o = x o sin θ o • sin φ o = y o cos θ o = z o ⎤ ⎥ ⎥ ⎦ ∈ Ω. (1.2)
The geometry of the coordinate system used throughout this thesis is illustrated in Figure 1.4 (1). We reserve the θ and φ notations to denote elevation and azimuthal angles, respectively.

The Rendering Equation In order to evaluate the incident radiance distribution function L during the lowpass filtering process, i.e., Equation (1.1), a physically based renderer computes the rendering equation [Kaj86]

L o = Ω + L i • f r • cos θ i • dω i . (1.3)
This equation gives the amount of light reflected towards o as a measure of the interactions occurring between an incident lighting distribution function L i ≥ 0 and a BRDF f r ≥ 0 over the set of unit directions defined on the hemisphere Ω + aligned with the tangent space of an infinitesimal surface element that projects into the pixel; the geometric configuration of the rendering equation is illustrated in Figure 1.4 (3). It is interesting to note that Equation (1.1) and Equation (1.3) look very similar; as such, the projected BRDF term of Equation (1.3), i.e., f r • cos θ i , can be interpreted as a filter applied to the incident radiance distribution function L i .

Displacement Mapped Surfaces The filter's angular support ΔW defines the footprint of the pixel in the 3D scene. Assuming that this footprint projects onto a displacement mapped surface (if not, split the footprint in two parts), then we know from the Nyquist sampling theorem that the amount of samples required to compute Equation (1.1) without artifacts should be at least twice the number of texels; the same logic applies to Equation (1.3). Hence, the further away we are from displacement mapped surfaces, the greater the rendering cost.

The practical contribution of this thesis is to show that this costly filtering operation can be accelerated with precomputations. In particular, we show that in the case where displacement mapped surfaces fill the pixels' footprints, Equation (1.1) can be approximated by the rendering equation directly by using a multiscale microfacet BRDF, effectively reducing a 4D integration problem into a 2D integration problem, which, furthermore, can be measured with efficient importance sampling; these ideas are illustrated in Figure 1.4 (4, 5), and will be formulated mathematically throughout this thesis.

Practical Considerations

Numerical Integrations Throughout this thesis, we often rely on numerical integrations to solve and/or validate certain equations. Whenever we face an integral of the form

I = D f (x) • dx, (1.4)
where D denotes the domain of integration, we always proceed in two steps. First, we reformulate the integral in the unit segment [0, 1). If D = [0, +∞) for instance, then we can first perform the change of variable x = tan θ (dx = sec 2 θ • dθ) so that the integration domain becomes [0, π/2)

I = π/2 0 f (tan θ) • sec 2 θ • dθ. (1.5)
We apply next a simple scaling operation, i.e., x = uπ/2 (dx = π 2 • du) to arrive at the unit segment

I = π 2 1 0 f (tan(uπ/2)) • sec 2 (uπ/2) • du. (1.6)
Note that this technique is also known as transform sampling. Whenever the transformation involves quantile functions, we say that the integral is importance sampled according to the distributions from which the quantile functions originate. Once the transformation complete, we rely on Euler's method to compute the integration.

Recalling that computing the integral of a function boils down to finding the area under its curve, Euler's technique consists in discretizing the integral into a sum of rectangular areas. The width of the rectangles is constant and set to one over the number of rectangles, while their heights are given by the function evaluated at their lower-left corner, which defines the location of the samples. This gives

I = 1 N π 2 N-1 ∑ i=0 f (tan(u i π/2)) • sec 2 (u i π/2), (1.7)
where u i = i/N. This sum is particularly simple to implement on a computer, as illustrated in Listing 1.1, which provides a pseudocode to evaluate the above equation with N = 100 samples. Whenever it is of interest to compute an integral numerically, the series of substitutions that leads to such configurations will be provided. Listing 1.1 -Simple numeric integration in C++.

Implementing Bidirectional Reflectance Distribution Functions Due to its role in the rendering equation, i.e., Equation (1.3), the BRDF is a central component of a physically renderer. Thoroughly described by Nicodemus et al. [NRH + 77], the BRDF implements the concept that a material can be described from the moment we know the fraction of light it reflects for any pair of unit directions involving an emitter and a receiver direction. We reserve the notation i to refer to the former and o for the latter. The emitter and receiver direction couple is used as the standard parameterization for the BRDF, and an implementation should be capable of returning its value for any i and o located on the hemisphere Ω + . Alternatively, it is also convenient to evaluate the BRDF with the Rusinkiewicz parameterization [START_REF] Szymon | A New Change of Variables for Efficient BRDF Representation[END_REF], which involves the halfway vector h ∈ Ω + and a difference vector d ∈ Ω + . The Rusinkiewicz parameters are linked to the standard parameters through the

relations ⎧ ⎪ ⎨ ⎪ ⎩ h = i + o i + o d = Rot y (-θ h ) • Rot z (-φ h ) • i ⇔ ⎧ ⎨ ⎩ i = Rot z (+φ h ) • Rot y (+θ h ) • d o = 2(i • h)h -i (1.8)
Finally, a BRDF implementation should expose a transform sampling procedure. This ability is crucial for solving the rendering equation quickly, as the incident radiance function L i is generally not known in advance.

In general, the transform sampling procedure accounts for the cosine term too, so that we end up importance sampling the projected BRDF term f r • cos θ i . All these features are present in the C++ interface of Listing 1.2, which we rely on to implement a multiscale microfacet BRDF library in the following chapters. Listing 1.2 -C++ interface of a BRDF.

To make the C++ interface more intuitive to the reader, Listing 1.3 provides an implementation for a Lambertian material of (wavelength-dependent) reflectance parameter

k d ∈ [0, 1], i.e., f r = k d π .
(1.9)

Note that such BRDFs are very limited in terms of artistic control and plausibility compared to microfacet BRDFs, which are the main focus of this thesis. Chapter 2

State-of-the-art Microfacet Theory

(1) BRDF at macroscopic scale.

(2) BRDF at microscopic scale.

(3) Geometry of mirror reflection.

(4) Visible microfacets.

(5) Illuminated microfacets.

(6) Microfacets involved in reflection. Beckmann

φ α = 0 • φ α = 45 • φ α = 90 • φ α = 135 • Figure 2
.6 -Effect of anisotropy orientation on a microfacet BRDF.

Torrance-Sparrow Equation

Microfacet theory [START_REF] Torrance | Theory for Off-Specular Reflection From Roughened Surfaces[END_REF] is a framework that models the reflectance of rough materials. It defines the BRDF as a global response of a microscopic surface, referred hereafter as microsurface, composed of microfacets that act as Fresnel mirrors when radiated by light rays. What makes this approach particularly elegant is that, under the assumption that single-bounce mirror reflection dominates on the microsurface, it transposes the problem of evaluating the response of the BRDF into that of determinining the amount of visible and illuminated microfacets whose orientations form a mirror-like reflection between given viewing and illumination directions, which is mathematically tractable; Figure 2.1 (1, 2, 3) illustrates this idea. While determining the unoccluded areas of an arbitrary surface is a difficult problem in general, some solutions can be derived for certain types of surfaces. For such surfaces, the BRDF can be derived analytically from three main terms, according to the

microfacet BRDF equation [WMLT07] f r = F • D • G 4 • cos θ i • cos θ o . (2.1)
We follow standard computer graphics notation for the various terms: F ∈ [0, 1] is the microfacet Fresnel term, D ≥ 0 the microfacet normal distribution function (NDF) of the microsurface, and G ∈ [0, 1] the geometric attenuation factor (GAF) of the microsurface. We now discuss each of these terms in detail. Microfacet Fresnel Term The microfacet Fresnel term F describes how microfacets act when a light ray hits them; it can be seen as the fundamental interaction between light and microsurfaces. In contrast to the other terms of Equation (2.1), F may vary across wavelengths and is thus responsible for explaining the color of the material. In general, F is parameterized by Rusinkiewicz's difference angle

F = F(θ d ).
(2.2)

Various analytic expressions for F have been proposed in the computer graphics literature. Among the most popular models, we find the Fresnel equation for unpolarized light, which was introduced in computer graphics by Cook and Torrance [START_REF] Cook | A Reflectance Model for Computer Graphics[END_REF], as well as its approximation derived by Schlick [START_REF] Schlick | An Inexpensive BRDF Model for Physically-based Rendering[END_REF]; their expressions are given in Table 2.I. When computational performance is critical, such as for video games, Schlick's approximation is often preferred over the unpolarized formulation, due to the simpler arithmetic operations involved.

Finally, note that if we set F = 1 at any light wavelength, then we can simulate a microsurface composed of ideal mirrors. In such cases, the microfacet BRDF takes on the analytic form

f r,id = D • G 4 • cos θ i • cos θ o . ( 2
F(θ d ) = 1 F(θ d ) = 1 2 (g -c) 2 (g + c) 2 1 + (c(g + c) -1) 2 (c(g -c) + 1) 2 F(θ d ) = F 0 + (1 -F 0 )(1 -c) 5
Notation η: index of refraction Microfacet Normal Distribution Function The microfacet NDF term D gives the proportion of microfacets oriented towards the halfway vector, i.e., the normals that connect the emitter and receiver directions through

θ d : difference angle c = cos θ d g 2 = η 2 + c 2 -1 F 0 = (η -1) 2 /(η + 1) 2
mirror reflection i = 2(o • h)h -o (2.4) h = i + o i + o (2.5)
Note that these relations are illustrated in Figure 2.1 (3). The microfacet NDF is a directional distribution, defined on the hemisphere Ω + oriented towards the up direction of the tangent frame

D = D(h).
(2.6)

The measure of D over Ω + gives the area of the microsurface [START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF], and should be finite to be geometrically valid, i.e.,

Ω + D(h) • dω h ∈ [1, +∞).
(2.7)

In addition, D should be normalized such that the projection of the microfacets onto the tangent plane has unit area [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF][START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF], i.e.,

Ω + D(h) • cos θ h • dω h = 1. (2.8)
Various analytic expressions of D can be found in the computer graphics literature [WMLT07, BSH12, LKYU12].

In the general case, analytic directional distributions are not straightforward to manipulate or design. In contrast to unit directions, the slopes of the microsurface, i.e., its derivatives, are much easier to study, as they live in the R 2 plane. In the Ω + set, normals and slopes are linked through the bijection

h = ⎡ ⎣ -tan θ h • cos φ h = xh -tan θ h • sin φ h = ỹh ⎤ ⎦ , h ∈ R 2 , (2.9)
whose inverse is

h = 1 x2 h + ỹ2 h + 1 ⎡ ⎢ ⎢ ⎣ -xh -ỹh 1 ⎤ ⎥ ⎥ ⎦ , h ∈ Ω + .
(2.10)

The bijection between normals and slopes is illustrated geometrically in Figure 2.7. From this bijection, we may define a microfacet NDF from a bivariate slope probability distribution function (PDF) P ≥ 0 such that The secant term is the Jacobian that converts the measure of microfacet slope probability into a microfacet normal distribution [START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF]. As long as P is a normalized PDF, i.e.,

D(h) = P( h) • sec 4 θ h . ( 2
R 2 P( h) • d h = 1,
(2.12) Equation (2.8) holds. Among the most popular microfacet slope distributions, we find the Beckmann distribution [START_REF] Beckmann | Shadowing of Random Rough Surfaces[END_REF], and the GGX distribution [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF]. The Beckmann distribution is simply a Gaussian PDF, while the GGX distribution emerges from the curvature of ellipsoids [START_REF] Trowbridge | Average Irregularity Representation of a Rough Surface for Ray Reflection[END_REF]; their respective expressions are provided in and the bivariate t distribution with two degrees of freedom in the mathematics literature [START_REF] Jones | Student's Simplest Distribution[END_REF].

(1) Set of normals.

(2) Normals seen as points.

(3) Normals on the lower hemisphere.

(4) The slopes are obtained by stereographic projection.

Figure 2.7 -Bijection between slopes and normals. (1) A set of normals defined within Ω + can also be seen as (2) a set of points lying on the boundary of Ω + . In order to compute the slope associated to each normal, we first negate the z component of these points, which creates (3) a set of points that lie on the boundary of Ω -.

(4) The stereographic projection of each of these points onto the plane (shown in green) yields the associated slope; note that this transformation is nonlinear.

Microfacet Geometric Attenuation Factor The microfacet GAF term G is responsible for modelling occlusion effects occurring on the microsurface. Occlusions can be due to shadowing, i.e., the emitter direction, and to masking, i.e., the receiver direction; this idea is illustrated in Figure 2.1 (4,5). Note that, although their appellation is different, computing the amount of either shadowing or masking boils down to determining the hidden portions of the microsurface from a specific direction. This problem is extremely difficult to solve for arbitrary surfaces, and solutions only exist for very specific surface profiles. In current state-of-the-art models, the GAF is given by the correlated bistatic Smith shadowing function [HMD + 14]

G(i, o) = G 1 (i) • G 1 (o) G 1 (i) + G 1 (o) -G 1 (i) • G 1 (o) , (2.13) Beckmann GGX P( h) = 1 π exp -x2 h -ỹ2 h P( h) = 1 π 1 (1 + x2 h + ỹ2 h ) 2
Notation h ∈ R 2 : slope of the halfway vector. where

G 1 ∈ [0, 1] denotes the Smith monostatic shadowing function [Smi67, Hei14] G 1 (k) = cos θ k Ω + kh • D(h) • dω h . (2.14)
Note that we use the notation kh to express vector dot products, and the underlining to denote a quantity clamped to zero if it is negative, e.g., kh = max(0, k • h). The Smith model assumes that the microsurface structure is independent from the GAF term, which essentially means that the microsurface can be approximated by a collection of flakes [START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF], as illustrated in Figure 2.8. While this approximation may appear to be quite crude, it matches experiments very well in practice [START_REF] Bourlier | Effect of Correlation Between Shadowing and Shadowed Points on the Wagner and Smith Monostatic One-dimensional Shadowing Functions[END_REF][START_REF] Bourlier | One-and Two-dimensional Shadowing Functions for Any Height and Slope Stationary Uncorrelated Surface in the Monostatic and Bistatic Configurations[END_REF][START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF]. We also mention here two less realistic GAF models, namely the ∨-groove model of Torrance and Sparrow [START_REF] Torrance | Theory for Off-Specular Reflection From Roughened Surfaces[END_REF], and the normal map model.

The former assumes that the microsurface is composed exclusively of symmetrical ∨-groove cavities, while the latter models the microsurface as a normal mapped surface, which is equivalent to neglecting occlusion effects;

both models are illustrated in Figure 2.8, and their equations are given in Table 2.IV. Figure 2.3 shows materials generated with different G terms, while F and D are fixed.

(1) Visible Smith microfacets.

(2) Visible ∨-groove microfacets.

(3) Visible normal map microfacets.

(4) Illuminated Smith microfacets.

(5) Illuminated ∨-groove microfacets. (6) Illuminated normal map microfacets.

Importance Sampling

Distribution of Visible Slopes Physically based renderers employ BRDFs for shading. On such platforms, shading is typically computed by solving the rendering equation [START_REF] Kajiya | The Rendering Equation[END_REF], which we already mentioned in the previous chapter (see Equation (1.3))

L o = Ω + L i • f r • cos θ i • dω i .
In the general case, the integral is computed numerically because the incident radiance L i is not known in advance. When relying on microfacet BRDFs, we can compute this integral efficiently by transforming uniformly distributed samples into the space of visible microfacet normals with a quantile transform [START_REF] Heitz | Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals[END_REF]. The PDF

Beckmann

GGX

G 1 (k) = 2 1 + erf(ν) + exp(-ν 2 ) ν √ π G 1 (k) = 2 1 + √ 1 + ν -2 Notation k ∈ Ω + : incident direction, of elevation angle θ k ∈ [0, π/2] ν = cot θ k
Table 2.III -The Smith masking term of analytic microfacet slope distributions.

Smith

∨-groove

Normal map

G = G 1 (i) • G 1 (o) G 1 (i) + G 1 (o) -G 1 (i) • G 1 (o) G = min 1, 2 z i • z h (i • h) , 2 z o • z h (o • h) G = z o • z h (o • h) Notation i ∈ Ω + : emitter direction, with elevation z i ≥ 0 o ∈ Ω + : receiver direction, with elevation z o ≥ 0 h ∈ Ω + : halfway direction, with elevation z h ≥ 0 Table 2.IV -Various microfacet GAF terms.
The Smith model is considered as the most plausible.

P vis > 0 of the distribution of visible microfacet slopes is defined as [HD14]

P vis ( h; k) = kh • P( h) • G 1 (k) cos θ h • cos θ k . (2.15)
Note that in order to simplify notations, we will make the direction parameter k implicit, thus writing, e.g.,

P vis ( h) instead of P vis ( h; k).
Quantile Functions Let P 2 ≥ 0 denote the marginal PDF

P 2 ( xh ) = ∞ -∞ P vis ( xh , ỹh ) • d ỹh .
(2.16)

The PDFs P vis and P 2 are linked through the relation

P vis ( xh , ỹh ) = P 2 ( xh ) • P 3 ( ỹh | xh ) ⇒ P 3 ( ỹh | xh ) = P vis ( xh , ỹh ) P 2 ( xh ) ,
where P 3 ≥ 0 is the PDF of ỹh conditioned on xh . Let F 2 ∈ [0, 1] and F 3 ∈ [0, 1] respectively denote the cumulative distribution function of P 2 and P 3 , i.e.,

F 2 ( xh ) = xh -∞ P 2 ( x) • d x, F 3 ( ỹh | xh ) = ỹh -∞ P 3 ( ỹ| xh ) • d ỹ,
and

Q 2 = F -1 2 and Q 3 = F -1 3 their respective quantile functions. Given realizations u 1 , u 2 ∈ [0, 1] of two
independent uniform variates, the variates obtained by the quantile transformation

h = ⎡ ⎣ Q 2 (u 1 ) Q 3 (u 2 |Q 2 (u 1 )) ⎤ ⎦
are distributed according to P vis [START_REF] Heitz | Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals[END_REF]. Note that if we transform these slopes into normals according to Equation (2.10), then the resulting normals are distributed according to the distribution of visible microfacet normals [START_REF] Heitz | Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals[END_REF], which has the PDF

D vis (h; k) = kh • D(h) • G 1 (k) cos θ k .
(2.17 Solving the Rendering Equation In order to compute the rendering equation numerically, we first reexpress the integral as a measure in terms of infinitesimal solid angles oriented towards halfway vectors. We accomplish this with the change of variable

ω i = 2(ω h • ω o )ω h -ω o (dω i = 4 • oh • dω h [WMLT07]), which yields L o = Ω + L i • f r • cos θ i • 4 • oh • dω h .
(2.18)

We then express this measure in slope space with the change of variable

ω h = ω h ( h) (dω h = cos 3 θ h • d h) L o = R 2 L i • f r • cos θ i • 4 • oh • cos 3 θ h • d h (2.19) = ∞ -∞ ∞ -∞ L i • f r • cos θ i • 4 • oh • cos 3 θ h • d xh • d ỹh .
(2.20)

From this double integral, we can then apply quantile transformations. We start with the x slope component.

Letting xh = Q 2 (u 1 ) (d xh = Q 2 (u 1 ) • du 1 = 1 P 2 (Q 2 (u 1 )) • du 1 ), we get L o = ∞ -∞ 1 0 L i • f r • cos θ i • 4 • oh • cos 3 θ h P 2 (Q 2 (u 1 )) • du 1 • d ỹh . (2.21)
We then apply the second quantile substitution. Letting ỹh

= Q 3 (u 2 | xh ) (d ỹh = Q 3 (u 2 | xh ) • du 2 , and so d ỹh = 1 P 3 (Q 3 (u 2 )|Q 2 (u 1 )) • du 2 ), we get L o = 1 0 1 0 L i • f r • cos θ i • 4 • oh • cos 3 θ h P 2 (Q 2 (u 1 )) • P 3 (Q 3 (u 2 )|Q 2 (u 1 )) • du 1 • du 2 , (2.22)
which, using the fact that P vis = P 2 • P 3 , as well as the definition of the microfacet BRDF, i.e., Equation (2.1), leads after simplifications to

L o = 1 0 1 0 L i • F(θ d ) • G(i, o) G 1 (o) • du 1 • du 2 .
(2.23)

We can then apply Euler's rule to compute the integral numerically.

Controlling Roughness

We can control the roughness of a microfacet BRDF intuitively by scaling the slopes of the microsurface, i.e., its derivatives. This feature is particularly relevant for artistic authoring. As scale increases, the surface becomes rougher. Conversely, as scale decreases, the surface becomes smoother. These ideas are illustrated in Figure 2.9.

Isotropic Roughness If we scale the microsurface isotropically by a factor α > 0, then the distribution of microfacet slopes is also scaled. This means that we can retrieve the value of D by simply offsetting the slope distribution. Thus, we can define a microfacet NDF with explicit control over isotropic roughness as

D(h; α) = P xh α , ỹh α sec 4 θ h α 2 . (2.24) Beckmann GGX F 3 ( ỹh | xh ) = 1 2 + 1 2 erf( xh ) F 3 ( ỹh | xh ) = 1 2 + σ ỹh πσ + π ỹ2 h + arctan( ỹh σ ) π Q 3 (u| xh ) = erfinv(2u -1) Q 3 (u| xh ) = σ ⎧ ⎨ ⎩ - 1 ribinv(2u, 3 2 , 1 2 ) -1 if u < 1 2 1 ribinv(2-2u, 3 2 , 1 2 ) -1 else.
Notation xh ∈ R: x component of the halfway slope u ∈ [0, 1): uniform variate erfinv: inverse of the error function ribinv: inverse of the regularized incomplete beta function

σ 2 = 1 + xh Table 2.
V -Conditional CDF and QF of analytic visible slope distributions.

(1) Smoother microsurface.

(2) Standard microsurface.

(3) Rougher microsurface.

Figure 2.9 -Geometric interpretation of roughness. Roughness can be interpreted as an amount of vertical stretch that is applied to a standard microsurface, whose roughness is α = 1 by convention.

Note that this expression of D still satisfies Equation (2.8) and thus remains physically sound. The visual impact of isotropic roughness on the microfacet BRDF is illustrated in Figure 2.4.

Elliptically Anisotropic Roughness

We can also provide control over the anisotropy of materials. For this, we can stretch and rotate the microsurface along two major directions. These transformations can be controlled by tweaking the parameters of an ellipse, i.e., two radii α 1 , α 2 > 0 and an angle φ α ∈ [-π/2, π/2]. Thus, we can define a microfacet NDF with explicit control over anisotropic roughness as

D(h; α 1 , α 2 , φ α ) = P -tan θ h • cos(φ h + φ α ) α 1 , -tan θ h • sin(φ h + φ α ) α 2 sec 4 θ h α 1 α 2 .
(2.25)

Note that this expression of D still satisfies Equation (2.8) and thus remains physically sound. Alternatively, we may write the transformation in terms of bivariate distribution parameters, which involve two scale parameters α x , α y > 0 and a correlation coefficient ρ ∈ (-1, 1)

D(h; α 1 , α 2 , φ α ) = P xh α x , α x • ỹh -ρα y • xh α x α y 1 -ρ 2 sec 4 θ h α x α y 1 -ρ 2 . (2.26)
Note that in order to simplify notations, we shall use the notation

P( h; α 1 , α 2 , φ α ) = P std xh α x , α x • ỹh -ρα y • xh α x α y 1 -ρ 2 1 α x α y 1 -ρ 2 , ( 2.27)
so that we can write

D(h; α 1 , α 2 , φ α ) = P( h; α 1 , α 2 , φ α ) sec 4 θ h , (2.28)
where P std ( h) = P( h; 1, 1, 0) gives the microfacet slope PDF of the standard surface, i.e., free of transformations. From an implementation standpoint, Equation (2.26) is preferable to Equation (2.25) whenever trigonometric functions induce performance issues. Note that the bivariate distribution parameters are linked to the ellipse parameters through the relations

α 2 x = 1 2 α 2 1 + α 2 2 + (α 2 1 -α 2 2 ) • cos(2φ α ) (2.29) α 2 y = 1 2 α 2 1 + α 2 2 -(α 2 1 -α 2 2 ) • cos(2φ α ) (2.30) ρ = α 2 1 -α 2 2 α x α y • cos φ α • sin φ α , (2.31)
whose inverses are

α 2 1 = 1 2 α 2 x + α 2 y + (α 2 x -α 2 y ) 2 + 4ρ 2 α 2 x α 2 y (2.32) α 2 2 = 1 2 α 2 x + α 2 y -(α 2 x -α 2 y ) 2 + 4ρ 2 α 2 x α 2 y (2.33) φ = arctan ⎛ ⎝ α 2 x -α 2 y -(α 2 x -α 2 y ) 2 + 4ρ 2 α 2 x α 2 y 2ρα x α y ⎞ ⎠ .
(2.34)

The visual impact of anisotropic roughness on a microfacet BRDF is illustrated in Figure 2.5.

Roughness Invariance of the Masking Term

Defining roughness as a slope-space scaling operation makes it possible to define a wide variety of microfacet NDFs from a single slope distribution by using offset evaluations; we refer to this property as roughness invariance. The Smith masking term is also subject to this property [START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF], which we prove here mathematically. We start with the definition of the Smith masking term, i.e., Equation (2.14), using an elliptically stretched NDF. We employ several substitutions and algebric manipulations provided in the parentheses on the right-hand side for convenience

G 1 (k; α 1 , α 2 , φ α ) = z k Ω + D(h; α 1 , α 2 , φ α ) • kh • dω h = z k R 2 P std xh α x , α x • ỹh -ρα y • xh α x α y √ 1-ρ 2 sec θ h α x α y √ 1-ρ 2 • kh • d h ⎛ ⎝ ω h = ω h ( h) ⇒ dω h = cos 3 θ h • d h⎞ ⎠ = z k R 2 P std xh α x , α x • ỹh -ρα y • xh α x α y √ 1-ρ 2 1 α x α y √ 1-ρ 2 • z k -x k xh -y k ỹh • d h ⎛ ⎜ ⎜ ⎝ sec θ h • h = ⎡ ⎢ ⎢ ⎣ -xh -ỹh 1 ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎠ = z k R 2 P std ( xm , ỹm ) • z k -x k xh -y k ỹh • d m ⎛ ⎜ ⎜ ⎝ xh = α x xm , ỹh = α y (ρ xm + 1 -ρ 2 ỹm ) ⇒ d h = α x α y 1 -ρ 2 • d m⎞ ⎟ ⎟ ⎠ = z k R 2 P std ( xm , ỹm ) • z k -a xm -b ỹm • d m ⎛ ⎝ a = α x x k + ρα y y k b = α y y k 1 -ρ 2 ⎞ ⎠ = z k R 2 P std ( xm , ỹm ) sec θ m • k m • d m ⎛ ⎜ ⎜ ⎝ k = 1 a 2 + b 2 + z 2 k ⎡ ⎢ ⎢ ⎣ a b z k ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎠ = z k Ω + D(m; 1, 1, 0) • k m • dω m ⎛ ⎝ m = m(ω m ) ⇒ d m = sec 3 θ m • dω m ⎞ ⎠ = G 1 (k ; 1, 1, 0) = G 1,std (k ).
This demonstration shows that the Smith term of a scaled surface is equal to the Smith term of an unscaled Smith term evaluated with an offset parameter; the relation is illustrated in Figure 2.10. (1) Visible scaled microsurface.

(2) Visible standard microsurface.

(3) Geometric link between directions.

(4) Visible scaled microsurface.

(5) Visible standard microsurface.

(6) Geometric link between directions.

Figure 2.10 -Roughness invariance of the monostatic shadowing term. Shadowing from any direction on (1, 4) a scaled surface can be transposed into (2, 5) shadowing on the standard surface due to another direction.

(3, 6) The incident directions are linked geometrically.

Importance Sampling with Roughness Parameters Since both NDF and GAF are roughness invariant, it follows by construction (see Equation (2.15)) that the distribution of visible slopes is also roughness invariant.

We can also exploit roughness invariance to produce samples that are distributed according to a scaled microsurface from an unscaled microsurface. Indeed, if m ∈ R 2 constitutes a sample of an unscaled visible microfacet slope distribution from direction k , then the slope

h = ⎡ ⎣ α x • xm α y (ρ • xm + 1 -ρ 2 • ỹm ) ⎤ ⎦ ∈ R 2 (2.35)
constitutes a sample of the visible microfacet slope distribution of the scaled surface from direction k; this operation is illustrated in Figure 2.11. We can thus exploit the quantile functions introduced in the previous section to produce visible slope samples of arbitrary roughness.

(1) Standard microsurface sample.

(2) Rougher microsurface sample.

(3) Geometric link between samples.

Figure Listing 2.10 -C++ implementation for radial microfacet slope distributions.

Radial Slope Distributions Some slope distributions such as the Beckmann and GGX slope distributions are radial distributions, i.e., their PDF can be expressed as a 1D function

P( h) = P r (r h ), (2.36)
where r h = tan θ h = x2 h + ỹ2 h is the radius of the slope in polar coordinates, and function

P r > 0 is normalized such that 2π 0 ∞ 0 r h • P r (r h ) • dr h • dφ h = 1.
(2.37)

We can exploit Equation (2.36) to improve the implementation of microfacet BRDFs built from radial slope distributions. Indeed, for such distributions, the evaluation of the microfacet NDF becomes The Smith masking term can also be simplified into a 1D function. To show this, we start from Equation (2.14) and reexpress the integral in slope space with the change of variables

D(h; α 1 , α 2 , φ α ) = P r ⎛ ⎜ ⎝ xh α x 2 + α x • ỹh -ρα y • xh α x α y 1 -ρ 2 2 ⎞ ⎟ ⎠ sec 4 θ h α x α y 1 -ρ 2 . ( 2
ω h = ω h ( h) (dω h = sec 3 θ h • d h) G 1 (k) = cos θ k R 2 z k -x k xh -y k ỹh • P( h) • d h .
Using the fact that

x k = sin θ k • cos φ k , y k = sin θ k • cos φ k ,
G 1 (k) = cot θ k R 2 cot θ k -(cos φ k • xh + sin φ k • ỹh ) • P( h) • d h .
We then reexpress the integral in polar coordinates,

letting xh = r h • cos φ h and ỹh = r h • sin φ h (d h = r h • dr h • dφ h ) G 1 (k) = cot θ k 2π 0 ∞ 0 cot θ k -r h (cos φ k • cos φ h + sin φ k • sin φ h ) • P r (r h ) • r h • dr h • dφ h , which, using the relation cos(a -b) = cos a • cos b + sin a • sin b, yields G 1 (k) = cot θ k 2π 0 ∞ 0 cot θ k -r h • cos(φ h -φ k )) • P r (r h ) • r h • dr h • dφ h . Letting φ h = φ m + φ k (dφ h = dφ m ),
(dr h • dφ m = cot θ h • d m) yields G 1 (k) = cot θ k 2π 0 ∞ 0 cot θ k -r h • cos φ m • P r (r h ) • r h • dr h • dφ m = cot θ k R 2 cot θ k -xh • P( xh , ỹh ) • d h ,
which may be written more concisely as

G 1 (k) = cot θ k cot θ k -∞ (cot θ k -xh ) • P 1 ( xh ) • d xh = G 1_radial (cot θ k ),
where P 1 is the marginal slope PDF

P 1 ( xh ) = +∞ -∞ P( xh , ỹh ) • d ỹh .
(2.39)

float_t radial::g1_smith_std(const vec3& h, const vec3& k) const { if (k.z <= DJB_EPSILON) return 0.0; if (k.z >= 1.0 -DJB_EPSILON) return 1.0; float_t cot_theta_k = k.z * inversesqrt(1.0 -k.z * k.z);
return g1_radial(cot_theta_k); } Listing 2.12 -C++ implementation for radial microfacet slope distributions.

Validating the Evaluation API Two simple tests can be devised to validate the evaluation API of a microfacet BRDF. The first test consists in checking the normalization of the NDF, i.e., making sure that Equation (2.8) holds

Ω + D(h) • cos θ h • dω h = 1.
The second test consists in checking the normalization of the distribution of the visible NDF, i.e., making sure that the following equation holds for any

k ∈ Ω + Ω + D vis (h; k) • dω h = 1.
To compute these integrals, we first express the integrals in spherical coordinates with the change of variable

ω h = ω h (θ h , φ h ) (dω h = sin θ h • dθ h • dφ h ).
We then make the changes of variable θ h = u 2 1 π/2 (dθ h = πu 1 • du 1 ) and φ h = 2πu 2 (dφ h = 2π • du 2 ) so that we can apply Euler's rule. We thus have

Ω + cos θ h • D(h) • dω h = 2π 0 π/2 0 cos θ h • D(h) • sin θ h • dθ h • dφ h = 2π 2 1 0 1 0 cos(u 2 1 π/2) • D(h) • sin(u 2 1 π/2) • u 1 • du 1 • du 2 = 2π 2 1 N i 1 N j N i ∑ i=0 N j ∑ j=0 cos(u 2 i π/2) • D(h) • sin(u 2 i π/2) • u i ,
and similarly

Ω + D vis (h; k) • dω h = 2π 2 1 N i 1 N j N i ∑ i=0 N j ∑ j=0 D vis (h; k) • sin(u 2 i π/2) • u i ,
where u i = i/N i , u j = j/N j , θ h = πu i /2, and φ h = 2πu j . Note that the change of variable θ h = u 2 1 π/2 produces more evaluations near the pole of Ω + , which usually results in faster convergence than with uniform hemispherical sampling, i.e., with the change of variable θ h = u 1 π/2. Chapter 3

Noncentral Microfacet Theory

(1) Microsurface.

(2) Magnified microsurface region.

(3) Sheared microsurface.

(4) Illuminated microfacets.

(5) Noncentral reflection.

(6) Equivalent central reflection. Similarly to, e.g., albedo texture maps, displacement maps must be filtered to match the resolution of the pixels they project into in order to guarantee consistent surface shading across scales [START_REF] Heckbert | Fundamentals of Texture Mapping and Image Warping[END_REF]. The simplest form of physically sound filtering is linear filtering, which consists in averaging the texels that contribute to the same pixel. But for this approach to work, the texels must interact linearly with shading [START_REF] Bruneton | A Survey of Non-linear Pre-filtering Methods for Efficient and Accurate Surface Shading[END_REF]. Put more intuitively, the average shading of a texture mapped surface should be similar to the shading of the average texture mapped surface. Although this restriction holds most of the time for albedo texture maps (see, e.g., the work of Heitz et al. for counterexamples [START_REF] Heitz | Filtering Color Mapped Textures and Surfaces[END_REF]), it does not for displacement maps. The reason why this is so is one of the main justifications for microfacet theory: the average BRDF of a displaced surface is not equal to the BRDF of its mean (hence smooth) surface.

The response of a proper displacement texture map filter should account for the BRDF of the surface, along with the view-dependent effects produced by the displacements. Since this is a complex task, we restrict for now the base surface BRDF to a Fresnel mirror, and postpone the treatment of the general case to Chapter 5.

Such a restriction is useful here because the filtering configuration we are now considering corresponds to what is handled by microfacet BRDFs with only two notable differences, which we address in the remainder of this chapter. First, displacement mapped surfaces are typically defined in a frame where their mean surface differs from the tangent plane; in Section 3.2, we extend microfacet theory to support noncentral microfacet distributions that can account for such configurations. Second, because in a rendering context we expect to be able to visualize arbitrary portions of the displaced surfaces at multiple resolutions, we need an effective way to retrieve the NDF of the filtered texels; in Section 3.3, we leverage the linear representation introduced by Olano and Baker [START_REF] Olano | LEAN Mapping[END_REF] to efficiently (though approximately) retrieve an anisotropic microfacet NDF at all scales. The combination of these tools constitutes the LEADR mapping framework, which we discuss from a more practical point of view in Section 3.4.

Noncentral Microfacet BRDF Equations

Scalar texture map Normal texture map Base surface Displaced surface Figure 3.2 -Overview of the displacement mapping technique. Displacement mapping consists in displacing a surface along its geometric normal using a scalar texture map, and using a normal texture map to compute the final shading. Here, the scalar texture map is visualized as grey-scale values, while the normal texture map is color-coded.

The purpose of a displacement map is to add high-frequency details on top of a coarser surface. The perturbations produced by the displacements result in a new surface, whose shading can differ quite significantly from the base surface; Figure 3.2 shows such an example. In order to derive a microfacet BRDF for such a surface, we wish to account for configurations where the mean displaced surface differs from the plane tangent of the base surface, in which the displacement is defined. That way, we can rely on the same representation at any scale to evaluate shading. We can accomplish this by exploiting another invariance property of microfacet BRDFs, namely shear invariance: the microfacet BRDF of a sheared surface, whose mean surface normal n ∈ Ω + differs from the tangent up direction, can be retrieved from the components of the same unsheared surface, whose mean surface normal corresponds to the tangent up direction. Whenever the components of the mean surface slope ñ = ( xn , ỹn ) ∈ R 2 associated to the mean surface normal n differs from zero, we say that the microfacet BRDF associated to such a surface is noncentral, with shear parameters ñ (the mean slope can also be interpreted as the amount of shear that has been applied to the microsurface). We also introduce the notions of meso-and macrosurfaces to respectively denote the mean microsurface and the tangent plane.

Noncentral Microfacet NDF If h is a slope defined on the sheared microsurface, then it may be written as the offset slope of an unsheared microsurface h = m + ñ, where ñ ∈ R 2 gives the shear applied to the microsurface (recall that shearing offsets the derivative). It follows that the NDF of any sheared surface may be retrieved from the NDF of the unsheared surface since m = hñ and hence

D(h; ñ) = P std ( xh -xn , ỹh -ỹn ) sec 4 θ h .
Note that shearing preserves the area of the surface projected onto the tangent plane, since for any ñ we have

Ω + D(h; ñ) • cos θ h • dω h = 1. (3.1)
Combining shear invariance with roughness invariance leads to the more general microfacet NDF expression

D(h; α 1 , α 2 , φ α ; ñ) = P std xh -xn α x , α x ( ỹh -ỹn ) -ρα y ( xh -xn ) α x α y 1 -ρ 2 sec 4 θ h α x α y 1 -ρ 2 , ( 3.2)
which gives the microfacet NDF of a surface that has been arbitrarily sheared and scaled from a single microfacet slope PDF. Noncentral Microfacet GAF Just like the microfacet NDF, the Smith masking term is also subject to shear invariance. We can prove this algebraically with a few manipulations provided in the parentheses on the right-hand side for convenience

G 1 (k; ñ) = z k Ω + D(h; xn , ỹn ) • kh • dω h = z k R 2 P std ( xh -xn , ỹh -ỹn ) sec θ h • kh • d h ⎛ ⎝ ω h = ω h ( h) ⇒ dω h = cos 3 θ h • d h⎞ ⎠ = z k R 2 P std ( xh -xn , ỹh -ỹn ) • z k -x k xh -y k ỹh • d h ⎛ ⎜ ⎜ ⎝ sec θ h • h = ⎡ ⎢ ⎢ ⎣ -xh -ỹh 1 ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎠ = z k R 2 P std ( xm , ỹm ) • z k -x k xh -y k ỹh • d m ⎛ ⎜ ⎜ ⎝ xh = xm + xn , ỹh = ỹm + ỹn ⇒ d h = d m ⎞ ⎟ ⎟ ⎠ = z k R 2 P std ( xm , ỹm ) • z k -x k xn -y k ỹn -x k xm -y k ỹm • d m = z k R 2 P std ( xm , ỹm ) • c -x k xm -y k ỹm • d m ⎛ ⎝ c = z k -x k xn -y k ỹn = kn • sec θ n ⎞ ⎠ = z k c c R 2 P std ( xm , ỹm ) • c -x k xm -y k ỹm • d m = z k c z k R 2 P std ( xm , ỹm ) sec θ m • k m • d m ⎛ ⎜ ⎜ ⎝ k = 1 x 2 k + y 2 k + c 2 ⎡ ⎢ ⎢ ⎣ x k y k c ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎠ = z k c z k Ω + D(m; 0, 0) • k m • dω m ⎛ ⎝ m = m(ω m ) ⇒ d m = sec 3 θ m • dω m ⎞ ⎠ = z k c G 1 (k ; 0) = z k c G 1,std (k ),
which reveal that the G 1 term of any sheared microsurface can be determined from the unsheared microsurface.

Combining shear invariance with roughness invariance leads to the more general microfacet Smith monostatic shadowing expression

G 1 (k; α 1 , α 2 , φ α ; ñ) = G 1,std ⎛ ⎜ ⎜ ⎝ 1 √ a 2 + b 2 + c 2 ⎡ ⎢ ⎢ ⎣ a b c ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎠ cos θ k kn • sec θ n , ( 3.3)
where

a = α x x k + ρα y y k , b = α y y k 1 -ρ 2 , and c = z k -x k xn -y k ỹn . Since the G 1 term is shear invariant, it
follows that the microfacet GAF term is also shear invariant. Note that both the ∨-groove and normal map

GAFs are also shear invariant. (1) Visible sheared microsurface.

(2) Visible standard microsurface.

(3) Geometric link between directions. Importance Sampling Since the microfacet NDF and the microfacet GAF terms are shear invariant, it follows by construction (see Equation (2.15)) that the distribution of visible slopes is also shear invariant. Thus, we can also adapt the importance sampling scheme of a standard microsurface to that of an arbitrary sheared and scaled microsurface. Indeed, if m constitutes a sample of the untransformed visible slope PDF from direction

(1) Standard microsurface sample.

(2) Sheared microsurface sample.

(3) Geometric link between samples.

Figure 3.4 -Sampling a sheared microsurface from a standard microsurface.

k , then the slope

h = ⎡ ⎣ α x • xm + xn α y (ρ • xm + 1 -ρ 2 • ỹm ) + ỹn ⎤ ⎦ (3.4)
constitutes a sample of the visible slope PDF from direction k of the sheared and scaled microsurface; this operation is illustrated in Figure 3 // return the associated normal return normalize(vec3(-tx_m, -ty_m, 1)); } Listing 3.3 -Importance sampling a microfacet BRDF with arbitrary roughness and shear parameters in C++.

LEAN/LEADR Mapping

Because the constant increase in texture and screen resolutions has been counteracting the growth in computational power over the years, texture filtering remains an expensive operation. This is why prefiltered textures such as MIP mapped texture maps [START_REF] Williams | Pyramidal Parametrics[END_REF][START_REF] Heckbert | Fundamentals of Texture Mapping and Image Warping[END_REF], which were introduced over thirty years ago, are still widely used and will probably remain popular for quite some time. MIP mapping consists in filtering an input texture at smaller resolutions prior to rendering and storing the output in texture maps of decreasing resolutions, called MIP maps. At render time, only the MIP texels that best match the on-screen resolution need to be evaluated, resulting in greater performance (due to the precomputations, and reduced input/output operations) and images with drastically reduced aliasing. Unfortunately, MIP mapping only works with linearly filterable textures, which means that in order to prefilter a displacement map, we first need to convert it into a linear representation.

As a first step towards this direction, we can use the observation that the filtered response of a displaced Fresnel mirror can be approximated by the evaluation of a noncentral microfacet BRDF. Note that if the GAF can accurately predict the occlusion effects produced by the displacements, then the filtering operation is exact. In any case, this means that we can transpose the problem of representing a displacement map into that of representing a noncentral microfacet slope NDF, i.e., Equation (3.2), from which we can derive a microfacet BRDF; the remainder of this section is dedicated to this problem.

Beckmann Surfaces Finding a generic and linear representation for arbitrary noncentral microfacet NDFs is too difficult. Nevertheless, we can make significant progress if we assume that the microfacet NDF encompasses a Beckmann slope PDF at all scales, i.e.,

P std ( h) ≈ 1 π exp(-x2 h -ỹ2 h ).
This assumption essentially restricts the behavior of the filtered surface to a specific distribution (Gaussian in this case), which simplifies the filtering problem to that of recovering roughness and shear parameters. For this purpose, Beckmann slope PDFs are extremely practical because their parameters can be retrieved from five moments:

E 1 := R 2 xh • P( h; α 1 , α 2 , φ α ; ñ) • d h = xn (3.5) E 2 := R 2 ỹh • P( h; α 1 , α 2 , φ α ; ñ) • d h = ỹn (3.6) E 3 := R 2 x2 h • P( h; α 1 , α 2 , φ α ; ñ) • d h = α 2 x 2 + x2 n (3.7) E 4 := R 2 ỹ2 h • P( h; α 1 , α 2 , φ α ; ñ) • d h = α 2 y 2 + ỹ2 n (3.8) E 5 := R 2 xh • ỹh • P( h; α 1 , α 2 , φ α ; ñ) • d h = ρα x α y 2 + xn ỹn . (3.9)
Note that we provide the derivations that lead to these terms in the next paragraph. The reason why this is useful is that, by exploiting the summation form of these moments, we can retrieve the parameters from any set of N > 0 slope texels that project into a pixel

E 1 = 1 N N-1 ∑ j=0 x j (3.10) E 2 = 1 N N-1 ∑ j=0 ỹ j (3.11) E 3 = 1 N N-1 ∑ j=0 x2 j
(3.12)

E 4 = 1 N N-1 ∑ j=0 ỹ2 j (3.13) E 5 = 1 N N-1 ∑ j=0 x j • ỹ j , (3.14)
where ( x j , ỹ j ) denotes the j-th slope sample of the surface. Since the summation form of the moments is linear with respect to slopes, we can store slopes (rather than normals) in a texture map and precompute MIP maps texels as the average of their higher resolution counterparts. This linear representation was first introduced by

Beckmann Moments

We can recover the scale and shear parameters of a Beckmann NDF from its first and second order slope moments. We prove here the results of Equations (3.5, 3.6, 3.7, 3.8, 3.9). First, we prove that the Beckmann slope PDF with expression

P std ( h) = 1 π exp(-x2 h -ỹ2 h )
is properly normalized, as it integrates to one

R 2 P std ( h) • d h = 1 π ∞ -∞ ∞ -∞ exp(-x2 h -ỹ2 h ) • d xh • d ỹh = 1 π 2π 0 ∞ 0 exp(-r 2 ) • r • dr • dφ ⎛ ⎝ xh = r • cos φ, ỹh = r • sin φ ⇒ d h = r • dr • dφ ⎞ ⎠ = 2 ∞ 0 exp(-r 2 ) • r • dr 2π 0 dφ = 2π = ∞ 0 exp(-x) • dx ⎛ ⎝ r = √ x ⇒ dr = 1 2 √ x • dx ⎞ ⎠ = 1. (3.15)
Next, we prove that the first order moments of the standard PDF integrate to zero. We have

R 2 xh • P std ( h) • d h = 1 π ∞ -∞ ∞ -∞ xh • exp(-x2 h -ỹ2 h ) • d xh • d ỹh = 1 √ π ∞ -∞ xh • exp(-x2 h ) • d xh 1 √ π ∞ -∞ exp(-ỹ2 h ) • d ỹh = 1 = 1 √ π ∞ 0 xh • exp(-x2 h ) • d xh - 1 √ π ∞ 0 xh • exp(-x2 h ) • d xh = 0. (3.16)
The same demonstration applies to the second first order moment, i.e., R 2 ỹh

• P std ( h) • d h = 0.
We also prove that the second order moments of the standard PDF integrate to 1 2 . We have

R 2 x2 h • P std ( h) • d h = 1 π ∞ -∞ ∞ -∞ x2 h • exp(-x2 h -ỹ2 h ) • d xh • d ỹh = 1 √ π ∞ -∞ x2 h • exp(-x2 h ) • d xh 1 √ π ∞ -∞ exp(-ỹ2 h ) • d ỹh = 1 = 1 2 √ π ∞ -∞ xh • f ( xh ) • d xh f ( xh ) = -2 xh • exp(-x2 h ) = 1 2 √ π ⎛ ⎜ ⎜ ⎜ ⎝ xh • f ( xh ) ∞ -∞ = 0 + ∞ -∞ exp(-x2 h ) • d xh = √ π ⎞ ⎟ ⎟ ⎟ ⎠ (integration by parts) = 1 2 .
(3.17)

The same demonstration applies to the second second order moment, i.e., R 2 ỹ2

h • P std ( h) • d h = 1 2 .
Now we can proceed to prove that

E 1 = R 2 xh • P( h; α 1 , α 2 , φ α ; ñ) • d h = R 2 (α x xm + xn ) • P std ( m) • d m ⎛ ⎜ ⎜ ⎝ xh = α x xm + xn , ỹh = α y (ρ xm + 1 -ρ 2 ỹm ) + ỹn ⇒ d h = α x α y 1 -ρ 2 • d m ⎞ ⎟ ⎟ ⎠ = R 2 xn • P std ( m) • d m (Equation (3.16)) = xn . (Equation (3.15))
The proof for E 2 = ỹn can be done similarly.

Next, we have

E 3 = R 2 x2 h • P( h; α 1 , α 2 , φ α ; ñ) • d h = R 2 (α x xm + xn ) 2 • P std ( m) • d m ⎛ ⎜ ⎜ ⎝ xh = α x xm + xn , ỹh = α y (ρ xm + 1 -ρ 2 ỹm ) + ỹn ⇒ d h = α x α y 1 -ρ 2 • d m ⎞ ⎟ ⎟ ⎠ = R 2 (α 2 x x2 m + x2 n + 2α x xn xm ) • P std ( m) • d m (Equation (3.16)) = R 2 α 2 x x2 m • P std ( m) • d m + x2 n (Equation (3.15)) = α 2 x 2 + x2 n .
(Equation (3.17))

The proof for E 4 = α 2 y 2 + ỹ2 n can be done similarly. Finally, we have ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( (

E 5 = R 2 xh • ỹh • P( h; α 1 , α 2 , φ α ; ñ) • d h = R 2 xh • ỹh • P std ( m) • d m ⎛ ⎜ ⎜ ⎝ xh = α x xm + xn , ỹh = α y (ρ xm + 1 -ρ 2 ỹm ) + ỹn ⇒ d h = α x α y 1 -ρ 2 • d m ⎞ ⎟ ⎟ ⎠ = R 2 α x α y (ρ x2 m + 1 -ρ 2 xm ỹm ) • P std ( m) • d m (Equation (3.16)) + R 2 α x ỹn xm • P std ( m) • d m (Equation (3.16)) + ( (
R 2 xn α y (ρ xm + 1 -ρ 2 ỹm ) • P std ( m) • d m (Equation (3.16)) + R 2 xn ỹn • P std ( m) • d m = R 2 α x α y ρ x2 m • P std ( m) • d m + xn ỹn (Equation (3.15)) = ρα x α y 2 + xn ỹn . (Equation (3.17))
Experimental Validation Our filtering scheme predicts the appearance of Beckmann microgeometries. In order to validate this property, we can compare the appearance yielded by our multiscale microfacet BRDF against supersampled renderings of Gaussian microgeometries; Figure 3.5 shows the result of such an experiment based on displaced spheres under directional lighting, which clearly demonstrates a good match for various illumination directions and roughness configurations between our filtering scheme (computed with one sample per pixel) and the reference images (computed using 256 samples per pixel). Note that more exhaustive tests are provided in the supplemental document of the original article [DHI + 13].

Practical Considerations

Precomputations The first thing we need to do to use LEADR mapping is to extract a LEAN map from a displacement map. If the displacement map carries a normal map, then the conversion is straightforward [START_REF] Olano | LEAN Mapping[END_REF].

Here, we show how to precompute LEAN maps from a displacement texture map that only consists of scalar displacement values. In this case, the texture can be interpreted as a discrete, scalar valued function ξ : [0, 1] 2 → R, living in the unit parametric space [0, 1] 2 . To find the slopes produced by ξ at any location in (u, v) ∈ [0, 1] 2 , we need to compute the gradient ∇ξ ≡ [∂ξ/∂u, ∂ξ/∂v] t . We can then express the slopes as

x(u, v) = ∂ξ(u, v) ∂u , ỹ(u, v) = ∂ξ(u, v) ∂v .
In order to compute these derivatives from texel space [0, w] × [0, h], where w and h respectively denote the width and height of the texture, the slopes must be scaled by the resolution of the texture x(p, q) = w • ∂ξ(p, q) ∂p , ỹ(p, q) = h • ∂ξ(p, q) ∂q .

We can use this last equation to approximate the gradient with forward or backward differencing-depending on whether the texel coordinate is even or odd-just as GPUs compute screen-space derivatives. This yields

x(p, q) ≈ w(ξ(p ± 1, q) -ξ(p, q)), ỹ(p, q) ≈ h(ξ(p, q ± 1) -ξ(p, q)).

In practice, however, it is better to favor central differencing

x(p, q) ≈ w 2 (ξ(p + 1, q) -ξ(p -1, q)), ỹ(p, q) ≈ h 2 (ξ(p, q + 1) -ξ(p, q -1))
because it results in smoother shading; Figure 3.6 compares a rendering obtained with both approaches. For the sake of completeness, conversions of displacement maps into normal maps and normal maps into the first MIP of a LEAN/LEADR map are provided in Listings 3.6 and 3.7, respectively (the rest of the LEAN/LEADR MIP maps can be computed using a simple GenerateMipmap() routine). Computing the moments in this way allows support for animated geometry without having to update the LEADR MIP map hierarchy during animation: the distribution is adapted on the fly from the distortions, and adequately modulates the resulting appearance; this is illustrated in Figure 3.7 and Figure 3.8.

Height Scaling One may wish to apply a scaling factor β ≥ 0 on displacements ξ to increase or reduce surface displacement. If β also exhibits low spatial variation across the surface, then it suffices to scale the first order moments E 1 and E 2 by β, and the second order moments E 3 , E 4 , and E 5 by β 2 . Note that this can be done Figure 3.9 -Impact of scale on the displacement mapping technique.

at runtime without having to update the LEADR MIP map hierarchy. Figure 3.9 shows an example of scaled displacement map that uses the same LEAN map.

A Note on CLEAN Mapping Some applications may not be able to afford five floating-point numbers per texel for specular antialiasing purposes. For such scenarios, Baker introduces CLEAN maps [START_REF] Baker | Spectacular Specular -LEAN and CLEAN Specular Highlights[END_REF],

which require only three floating-point numbers per texel at the cost of slightly less satisfying filtering quality (anisotropy is lost). When filtering with a CLEAN map, we are assuming that slope roughness is isotropic, i.e., α x = α y , and ρ = 0. The xn and ỹn parameters are determined similarly to a LEAN/LEADR map using the first slope moments in the x and y directions, which constitute the first two floating-point numbers of a CLEAN map. In order to retrieve the isotropic roughness parameter α = α x = α y , Baker suggests to store and filter the squared magnitude of the second moments, i.e., x2 h + ỹ2 h . The roughness is then determined as follows

α 2 = R 2 ( x2 h + ỹ2 h ) • P( h) • d h =CLEAN 3rd component -x2 n -ỹ2 n . (3.18)
Depending on the textures however, this approach may result in aliasing because some directions may require higher roughness values. In order to guarantee antialiased highlights, the lobe of the slope distribution should be as large as the major radius of the ellipse produced by a LEAN/LEADR map, i.e., max(α 1 , α 2 ). If we set the third component of the CLEAN map to yield the roughness value max(α 1 , α 2 ), then the resulting highlights should always be large enough to avoid aliasing. Note that this is only a conjecture, since this approach has not been tested.

GPU Acceleration

The LEADR mapping technique can be implemented on the GPU to produce high-quality renderings in real time with environment maps. Since a typical GPU rendering pipeline can only afford a limited number of shading evaluations per pixel, we can devise a filtered sampling approach that guarantees artifactfree images at the cost of slight bias, as follows. First, we recall from the previous chapter (see Equation (2.23))

that we can rewrite the rendering equation as

L o = 1 N 2 N-1 ∑ j=0 N-1 ∑ k=0 L i • F(θ d ) • G(i, o) G 1 (o) , having i = 2(o • h)h -o,
where h is the microfacet normal sample simulated with the uniform variates u 1 = j/N and u 2 = k/N. By using prefiltered samples of L i , we can emulate a smooth reconstruction filter that is guaranteed to produce antialiased results; this technique is also known as filtered importance sampling [START_REF] Colbert | GPU-based importance sampling[END_REF][START_REF] Křivánek | Real-time Shading with Filtered Importance Sampling[END_REF]. The problem is then to determine the size of the filter. A simple solution is to compute it according to the solid angle ω covered in the

incident direction i ω ≈ cos θ 3 o • A,
where A ∝ (max(α 1 , α 2 )/N) 2 is the surface area of a sample in slope space [DHI + 13]. This guarantees that the solid angles overlap across neighboring samples and avoids banding artifacts for any value of N ≥ 1. This technique was used to produce the images in, e.g., Figure 3.10 and Figure 3.11.

Chapter 4

Inverted Microfacet Theory

(1) BRDF at macro-vs microscale.

(2) Backscattering at macroscale.

(3) Backscattering at microscale. (1) Can we solve the inverse problem, i.e., is it possible to extract a microsurface from a BRDF? We can tackle this problem by considering the BRDF in (2, 3) backscattering configurations (o = i = h). For such configurations, occlusions due to (4) masking and (5) shadowing become equal, and (6) the evaluation of the BRDF is directly proportional to the fraction (with respect to the entire microsurface) of microfacets that are visible from the backscattering direction. This observation is the key to solve the inverse problem.

Motivation

This chapter encompasses the main contributions of the paper entitled "Extracting Microfacet-based BRDF Parameters from Arbitrary Materials with Power Iterations" [DHI + 15], which shows how to convert the reflectance of an arbitrary material into a microfacet BRDF. In order to alleviate the Fresnel mirror restriction that we imposed to displaced surfaces in the previous chapter, we will approximate the material of the undisplaced surface with a microfacet BRDF. But for this approach to be general (and hence practical), we need to design a procedure that can recover the components of a microfacet BRDF by solely evaluating the input material. Traditionally, this is done with a microfacet BRDF fitting procedure.

Existing microfacet-based fitting procedures may be classified in two categories. On the one hand, we find methods that typically introduce a new parametric microfacet normal distribution, and rely on sophisticated optimizations to fit its parameters [START_REF] Ngan | Experimental Analysis of BRDF Models[END_REF][START_REF] Löw | BRDF Models for Accurate and Efficient Rendering of Glossy Surfaces[END_REF][START_REF] Bagher | Accurate Fitting of Measured Reflectances using a Shifted Gamma Micro-facet Distribution[END_REF]. Although this approach has been used successfully to fit measured materials, current methods suffer from two main problems. First, optimizations are not guaranteed to converge to a valid solution. Second, available parametric models are isotropic and single lobed, which limits the range of materials that they can fit accurately. On the other hand, we find methods that extract an arbitrary microfacet normal distribution into a table [AP07, WZT + 08]. While this approach is considerably more general than the ones from the first category, as it naturally supports complex lobes and anisotropy, existing methods are not entirely satisfactory. Indeed, current methods either rely on crude approximations for visibility effects [START_REF] Ashikhmin | Distribution-based BRDFs[END_REF], or require a complex implementation involving sophisticated optimizations [WZT + 08].

In the following sections, we introduce a new tabulation strategy that takes an arbitrary material as input, and extracts a microfacet normal distribution from backscattering configurations, and then a Fresnel function.

Our approach combines the simplicity of the method of Ashikmnin and Premože [AP07], i.e., it is free from sophisticated optimizations, while retaining an accurate model for visibility effects. Specifically, we use the Smith microfacet visibility model, and also provide the derivations for the ∨-groove and normal map models, which we introduced in Chapter 2. In addition, we provide straightforward routines to convert our fits to Beckmann and GGX microfacet roughness parameters [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF]. Finally, we discuss how implementations may take advantage of microfacet slope space to precompute importance sampling tables that support spatially varying roughness for our fitted materials.

Microfacet Terms Extraction

Backscattering Configurations In order to retrieve a microfacet slope distribution from an input material, we are going to use a simplified form of Equation (2.1). Specifically, we focus on backscattering configurations, which reduces the dimensionality of the BRDF from 4D to 2D, and simplifies the Fresnel term to a constant; Figure 4.1 (2, 3) shows the geometry of a backscattering configuration. Indeed, in such configurations, we have i = o = h, as well as θ d = 0; for convenience, we write F 0 = F(0). Equation (2.1) predicts that the BRDF takes the form

f r (o, o) = F 0 • D(o) • G(o, o) 4 • cos 2 θ o . (4.1)
However, this expression is not entirely accurate because the GAF term G overestimates occlusion effects.

When i = o, shadowing and masking are fully correlated, and the GAF should degenerate into the monostatic

form [Hei14] G(o, o) = G 1 (o). (4.2)
The geometry of shadowing and masking in a backscattering configuration is illustrated in Figure 4.1 (4,5).

The accurate form of the GAF term leads to the microfacet backscattering equation

f r (o, o) = F 0 • D(o) • G 1 (o) 4 • cos 2 θ o . ( 4.3)
With this equation at hand, we may now proceed to the extraction of a microfacet slope distribution P, given an input material f r . Note that the inability of Equation (2.13) to degenerate into Equation (4.2) is known as the "hotspot" problem, and remains an open problem. In practice, we use Equation (4.2) for fitting purposes and Equation (2.13) for rendering. Our results thus overestimate occlusion effects in hotspot configurations by the

factor 1/(2 -G 1 ) ∈ [0.5, 1].
Eigensystem Construction In our microfacet BRDF fitting problem, we are given an input material f r and asked to extract a microfacet Fresnel term and a microfacet slope distribution. As a first step towards this direction, we swap the product F 0 • D and f r to the other side of the equality in the microfacet backscattering equation, i.e., Equation (4.3). This yields

F 0 • D(o) = 4 • f r (o, o) • cos 2 θ o G 1 (o) . ( 4.4)
Next, by replacing D and G 1 in Equation (4.4) by their definitions, i.e., Equation (2.11) and Equation (2.14) respectively, we get an equation for the microfacet slope PDF

F 0 • P(õ) = Ω + K(o, h) • P( h) • dω h , ( 4.5)
where

K(o, h) = 4 • f r (o, o) • cos 5 θ o • oh • sec 4 θ h . ( 4.6)
Equation (4.5) is not trivial to solve: in the mathematics literature, it belongs to the family of multivariate Fredholm equations of the second kind with kernel K [START_REF] Andrei | Handbook of Integral Equations[END_REF]. In our case, the form of the kernel is not known in advance due to the input material f r , so we solve this equation numerically by discretizing Equation (4.5)

with a quadrature rule. Letting w j denote the j-th quadrature rule weight and i = 1, • • • , N, we obtain the new relation

F 0 • P(õ i ) = N ∑ j=1 w j • K(o i , h j ) • P( h j ), (4.7)
where õ1 = h1 , • • • , õN = hN are the (slope) quadrature points, located in R 2 . Now, letting p denote the discretized PDF vector p = (P(õ 1 ), • • • , P(õ N )) t , and K the matrix

K = ⎡ ⎢ ⎢ ⎢ ⎣ w 1 K(o 1 , h 1 ) • • • w N K(o 1 , h N ) . . . . . . . . . w 1 K(o N , h 1 ) • • • w N K(o N , h N ) ⎤ ⎥ ⎥ ⎥ ⎦ , ( 4.8)
Equation (4.7) rewrites as the eigenvalue problem

F 0 • p = K • p. (4.9)
This result shows that the problem of retrieving a microfacet slope PDF given an input material BRDF f r translates into finding an eigenvector p whose components are all nonnegative.

Resolution via Power Iterations Since the entries of matrix K are all nonnegative, the Perron-Frobenius theorem [START_REF] Perron | Zur Theorie der Matrices[END_REF][START_REF] Frobenius | Über Matrizen aus nicht Negativen Elementen[END_REF] states that it always has a unique eigenvector p whose values are all nonnegative. It is thus a valid solution for the microfacet slope distribution P. By solving Equation (4.9), we therefore have the guaranteed ability to compute a valid microfacet slope PDF from an input material BRDF. Furthermore, the theorem states that this eigenvector is associated to the largest eigenvalue of the matrix. As such, we can compute it straightforwardly with the power iteration method.

The power iteration method is based on the property that the eigenvector associated to the largest eigenvalue of a matrix emerges after successive multiplications with a vector. In our implementation, we initialize p =

(1, • • • , 1) t , and successively multiply it by K. In practice, we use only four successive multiplications, which has turned out to be sufficient for all our test cases. Once the vector has been determined, we build a continuous PDF by linearly interpolating the values of p, and normalize it to satisfy Equation (2.12). We store the result in a table, which completes the extraction of P. Algorithm 1 provides pseudocode for our method.

Algorithm 1 Extract P function EXTRACT_P( f r , N) for each i, j ∈ [1, N] do Build kernel matrix K i, j ← w j 4 f r (o i , o i ) cos 5 θ o i o i h j sec 4 θ h j p ← (1, • • • , 1) t for 0 ≤ i < M do Power iterations (we set M = 4) p ← K • p P ← normalize(p)
With the extraction process of P complete, we can now evaluate the microfacet NDF D extracted from the input material thanks to Equation (2.11). It follows from Equation (2.14) that we can also compute the Smith term G 1 . At this point, we can already create a fully functional microfacet BRDF whose microfacets act as "ideal" mirrors, i.e.,

f r,id (i, o) = D(h) • G(i, o) 4 • cos θ i • cos θ o . ( 4.10)
Note that this equation is a special case of Equation (2.1), where F(θ d ) = 1 for any θ d . We can thus turn to the problem of retrieving the microfacet Fresnel term in order to complete the fitting process.

Fresnel Extraction In order to extract the Fresnel term, we compute for each color channel the average ratio between the input material BRDF and Equation (4.10) over all possible θ d configurations, i.e.,

F(θ d ) = E f r (i, o) f r,id (i, o) | ih = cos θ d , (4.11)
where we use the notation E[x|y] to denote the expectation of variable x over the domain that satisfies condition y.

Algorithm 2 provides pseudocode for our method.

The main advantage of our approach is that it is fully automatic. It is very accurate when the behavior of

Algorithm 2 Extract F function EXTRACT_F( f r , f r,id ) for θ d ∈ [0, π/2] do F(θ d ) ← 0 N ← 0 for φ d , φ h ∈ [0, 2π], θ h ∈ [0, π/2] do i ← from_half_diff(h, d) o ← reflect(i, h) F(θ d ) ← F(θ d ) + f r (i, o)/ f r,id (i, o) N ← N + 1 F(θ d ) ← F(θ d )/N
the input material accords to microfacet theory, i.e., when f r is roughly equal to Equation (2.1). When this condition does not hold, our algorithm will propagate the fitting errors into the Fresnel term. Such situations will arise for, e.g., materials with directionally dependent albedos, such as certain fabrics or car paints. In such cases, the extracted behavior of F differs significantly from what is predicted by the Fresnel law, and the term should be regarded as a residual function instead of an actual Fresnel function.

Optimization: Eigensystem Construction for Isotropic Materials Although our fitting process is already complete, we introduce here an optimization for the extraction of P that works for isotropic materials. If the input material is isotropic, then the microfacet NDF D depends only on the elevation angle θ h . It follows from Equation (2.11) that the microfacet slope PDF is also isotropic, which implies that it may be expressed as a 1D radial function

P( h) = g(θ h ). (4.12)
In such cases, the problem of retrieving the 2D function P, i.e., Equation (4.5), simplifies to that of finding the 1D function g. This problem also translates into a (univariate) Fredholm equation of the second kind Inversion with the ∨-Groove GAF We give here the procedure to extract a microfacet slope PDF assuming a ∨-groove GAF. In this case, the inversion procedure is much simpler than for the Smith GAF and only requires a few algebric manipulations

F 0 • g(θ o ) = π/2 0 K (θ o , θ h ) • g(θ h ) • dθ h , (4.13) with kernel K K (θ o , θ h ) = 2π 0 K(o, h) • sin θ h • dφ h . ( 4
P( h) = 4 • cos 6 θ h • f r (h, h) F 0 • min (1, 2 • cos 2 θ h ) , ( 4.16)
which may also be rewritten in the numerically stable form

P( h) = 4 F 0 • cos 4 θ h • f r (h, h) • max cos 2 θ h , 1 2 . ( 4
α 2 x = 2 R 2 x2 h • P( h) • d h α 2 y = 2 R 2 ỹ2 h • P( h) • d h ρα x α y = 2 R 2 xh • ỹh • P( h) • d h.
Note that if P is a Beckmann slope PDF, then our conversion is exact. Conversion to GGX Roughness Note that the 2nd order moments diverge with GGX. We propose an alternative estimation to retrieve the parameters of the scale matrix

α x = R 2 | xh | • P( h) • d h α y = R 2 | ỹh | • P( h) • d h β 1 = R 2 xh ỹh x2 h + ỹ2 h • P( h) • d h β 2 = R 2 ỹ2 h x2 h + ỹ2 h • P( h) • d h ρ = α y α x β 1 β 2 1 + β 2 2 .
this error to the quadrature rules that we use to solve integral Equations (4.5) and (4.13), and consider it as the minimum error produced by our algorithms. The anisotropic case produces the same amount of relative error. α 0.01 0.02 0.05 0.15 δ max 0.03 0.004 0.002 0.0005 Table 4.I -Maximum relative error in backscattering between Beckmann BRDFs with varying roughness and their respective fits computed with our algorithm.

Fitting Isotropic Materials We proceeded to an extensive fitting comparison against the state-of-the-art parametric model, referred hereafter as SGD, of Bagher et al. [START_REF] Bagher | Accurate Fitting of Measured Reflectances using a Shifted Gamma Micro-facet Distribution[END_REF] using the MERL material database of Matusik et al. [START_REF] Matusik | A Data-Driven Reflectance Model[END_REF]. To compute our fits, we initialized our algorithm with N = 90 and made sure that the entire backscattering data was sampled only once. At 32-bit floating-point precision, each of our fits requires 2.1 KB of memory. The exhaustive tests are provided in the supplemental document of the original article [DHI + 15], which also includes more detailed numerical analyses as well as delta-E difference images. Figure 4.3 shows some comparative renderings of both methods as well as with the ABC microfacet model [START_REF] Löw | BRDF Models for Accurate and Efficient Rendering of Glossy Surfaces[END_REF] against the reference for a few materials, using 512 samples per pixel. Note that for the twolayer-gold and changing-paint1 materials, the SGD fitting optimization failed and resulted in flawed images. This example emphasizes one of the strengths of our fitting method over optimization techniques, since, as we showed in the previous section, our fits can not result in such failures. Note that these materials were not the only ones affected by this issue in the database. In general, we believe our method is qualitatively superior to SGD and on par with ABC for metallic materials. Differences with SGD and ABC are most visible in Figure 4.3 at grazing angles. For most other materials in the supplemental document of the original article [DHI + 15], our fits are either on par or slightly below the SGD model. While our observation is mainly qualitative, it is also in agreement with the average delta-E difference image of our supplemental document, which is illustrated in Worse Isotropic Fitting For certain isotropic materials of the MERL database, we noticed that our method could produce fits that were qualitatively less satisfying than previous work. Figure 4.3 carries a few such materials, which can also be found in our supplemental document. In this particular figure, our fit of the aluminaoxide material is worse than those of the ABC and SGD models. We see two reasons why our method would produce less satisfying fits. The first reason is due to the input material itself: since our method extracts the microfacet NDF from backscattering data exclusively, it is highly sensitive to the quality of such configurations. Thus, if backscattering is poorly acquired, then our method will fail at reproducing the input BRDF faithfully.

The second reason is due to our BRDF model: our model is based on a microfacet BRDF model alone and, as such, is limited to material behaviors that are close to what is predicted by the equations. Layered and/or composite materials (e.g., alumina-oxide in Figure 4.3) as well as strong Lambertian component (e.g., some paints and acrylics) tend to be qualitatively less satisfying than the SGD fits in our supplemental document.

For the specific case of isotropic materials with poor backscattering data and/or strong Lambertian components, fitting methods based on optimizations [NDM05, LKYU12, BSH12, WZT + 08] should also perform better than our method in general. As for layered and color-changing materials (e.g., changing-paint1 in Figure 4.3), they

remain a challenging open problem.

Fitting Anisotropic Materials Our method also supports anisotropic materials. To review its performance, we tested some highly anisotropic materials from the recent database of Filip and Vavra [START_REF] Filip | Template-Based Sampling of Anisotropic BRDFs[END_REF]. To compute our fits, we initialized our algorithm with θ N = 90 and φ N = 90. At 32-bit floating-point precision, each of our fits requires 6 KB of memory. Our results were computed using 512 samples per pixel and are illustrated in Note however that, as we predicted in the previous section, our fit of the fabric106 material failed at reproducing the directionally dependent albedo exhibited by the reference.

Speed Our fitting algorithms are very fast: it takes us less than 1 second to fit an isotropic material from the MERL database of Matusik et al. [START_REF] Matusik | A Data-Driven Reflectance Model[END_REF], and less than 20 seconds for an anisotropic material of Filip and

Vavra [START_REF] Filip | Template-Based Sampling of Anisotropic BRDFs[END_REF]. Naturally, fitting performance depends on the resolution of the tables, i.e., on N. We measured the impact of such a factor for both isotropic and anisotropic algorithms. Results are plotted in Fresnel term on an Intel 2.4GHz Core i5 CPU. We believe such timings make our method much faster than previous work.

Memory Because our representation is roughness invariant, it allows us to create a multitude of materials at constant memory costs. As an example, we rendered the scene illustrated in Figure 4.6 using 512 samples per pixel and a few KB of memory. Such rendering configurations are only possible with slope-space tables. Otherwise, the per-pixel sampling rate and/or memory consumption to store importance sampling tables should be increased. For roughness mapped models such as ones shown in Figure 4.6, where the number of different materials is very large, such approaches would have been particularly impractical. Alongside analytic microfacet BRDF models, we believe our tabulation strategy is the first to support such rendering configurations trivially.

Conversion to Analytic BRDFs Although our memory footprint is constant per fitted material, some applications might not be able to afford the storage of precomputed tables. This is typically the case in real-time rendering contexts, where analytic microfacet BRDF models such as the Beckmann or GGX models are a necessity. We showed that our tabulated microfacet slope PDF can be converted straightforwardly to either Beckmann or GGX roughness parameters using slope moments. With such analytic microfacet models, only the Fresnel table needs to be stored. In order to quantify the loss in fitting quality compared to our tabu- lated representation, we placed some renderings obtained with the analytic models next to our tabulated fits in Figure 4.3. We noticed that the GGX model performs generally better than the Beckmann model (see also

Implementation Details

Online Computations We implemented BSDF plugins in the Mitsuba renderer [START_REF] Wenzel | Mitsuba Renderer[END_REF] to use our microfacet model. Our plugins implement the functions eval sample and pdf, which are called by Mitsuba's Monte-Carlo integrator. Because our tables are roughness invariant, we produce a wide variety of roughness effects on the fly by scaling the lookup parameters by α x , α y , and ρ. Slope PDF Normalization We store the tabulated radial slope PDF. Once we have retrieved the coefficients with the power iteration method, we normalize the PDF by dividing each coefficient of the table by the integral Smith Masking Term We compute the Smith masking term with the following integral

k = 2π 0 +∞ 0 r h • P r (r h ) • dr h = 2π π/2 0 tan θ h • P r (tan θ h ) • sec 2 θ h • dθ h r h = tan θ h ⇒ dr h = sec 2 θ h • dθ h = 2π 2 1 0 tan(u 2 π/2) • P r (tan(u 2 π/2)) • sec 2 (u 2 π/2) • u • du θ h = u 2 π/2 ⇒ dθ h = πu • du = 2π 2 N N ∑ i=0 tan(u 2 i π/2) • P r (tan(u 2 i π/2)) • sec 2 (u 2 i π/2) • u i . (u i = i/N)
G 1 (k) = cos θ k Ω + kh • D(h) • dω h = cos θ k 2π 0 π/2 0 kh • D(h) • sin θ h • dθ h • dφ h ⎛ ⎝ ω h = ω h (θ h , φ h ) ⇒ dω h = sin θ h • dθ h • dφ h ⎞ ⎠ = cos θ k 2π 2 1 0 1 0 kh • D(h) • sin(u 2 1 π/2) • du 1 • du 2 ⎛ ⎝ θ h = u 2 1 π/2, φ h = 2πu 2 ⇒ dθ h = πu 1 • du 1 , dφ h = 2π • du 2 ⎞ ⎠ = N i N j • cos θ k 2π 2 ∑ N i i=0 ∑ N j j=0 kh • D(h) • sin(u 2 i π/2) . (u i = i/N i , u j = j/N j )
Note that for isotropic materials, the Smith term is independent from the azimuthal component of the incident vector φ k . So in practice, we only store the results for the case φ k = 0.

F 1 ∈ [0, 1] associated to P 1 can also benefit from symmetry. Indeed, we can write it as

F 1 (x) = 1 2 + sign(x) • |x| 0 P 1 (y) • dy.
So in practice, we precompute for x ∈ [0, +∞) the integral

x 0 P 1 (y) • dy = arctan x 0 P 1 (tan θ) • sec 2 θ • dθ y = tan θ ⇒ dy = sec 2 θ • dθ = π 2 π arctan x 0 P 1 (tan(u 2 π/2)) • sec 2 (u 2 π/2) • u • du θ = u 2 π/2 ⇒ dθ = πu • du = π N j ∑ i=0 P 1 (tan(u 2 i π/2)) • sec 2 (u 2 i π/2) • u i . u i = i/N, j = 2N π arctan x
Once we have computed the CDF, we can compute the quantile function

Q 1 = F -1 1 .
Since F 1 is symmetric, we only need to compute and store half of the quantile function.

Marginal Visible Slope Quantile Function We exploit a property of isotropic materials to make the precomputation of the quantile function of the marginal visible slope PDF practical. For isotropic materials, the distribution of visible slopes is independent of the azimuthal angle of the incident vector. Hence, by choosing φ k = 0, we can express the marginal visible slope PDF as

P vis ( h; k) = max(0, 1 -tan θ k • xh ) • P( h) • G 1 (k).
From this equation, it follows that the marginal visible slope PDF P 2 > 0 is

P 2 ( xh ; k) = ∞ -∞ P vis ( xh , ỹ; k) • d ỹ = max(0, 1 -tan θ k • xh ) • P 1 ( xh ) • G 1 (k).
Assuming that xh < cot θ k , the marginal visible slope CDF is defined as

F 2 ( xh ; k) = xh -∞ P 2 ( x; k) • d x = G 1 (k) xh -∞ (1 -tan θ k • x) • P 1 ( x) • d x = G 1 (k) F 1 ( xh ) 0 (1 -tan θ k • Q 1 (p)) • d p x = Q 1 (p) ⇒ d x = 1 P 1 (Q 1 (p)) • d p = F 1 ( xh ) • G 1 (k) 1 -tan θ k • 1 0 Q 1 (F 1 ( xh ) • u) • du (p = F 1 ( xh ) • u ⇒ d p = F 1 ( xh ) • du) = F 1 ( xh ) • G 1 (k) 1 - tan θ k N N ∑ i=0 Q 1 (F 1 ( xh ) • u i ) . (u i = i/N)
In order to invert F 2 , we use the fact that the CDF may be written as a composition

F 2 ( xh ; k) = g(F 1 ( xh ); k) ⇒ Q 2 (u) = Q 1 (g -1 (u; k)), where g(x; k) = x • G 1 (k) 1 -tan θ k • 1 0 Q 1 (x • u) • du .
Conditional Visible Slope Quantile Function The conditional visible slope PDF is defined as

P 3 ( ỹh | xh ; k) = P vis ( xh , ỹh ; k) P 2 ( xh ; k) = P( xh , ỹh ) P 1 ( xh ) .
For isotropic materials, P 3 is a symmetric PDF. Therefore, we can apply some of the treatments we used to compute F 1 and Q 1 for F 3 and Q 3 , respectively. Note that computing F 3 by integrating the ratio of P and P 1 becomes numerically unstable as soon as we consider large xh values. To avoid this issue, we proceed as follows. Let F ∈ [0, 1] denote the bivariate CDF associated to P. We have

F( xh , ỹh ) = xh -∞ ỹh -∞ P( x, ỹ) • d x • d ỹ = F 1 ( xh ) 0 F 1 ( ỹh ) 0 P(Q 1 (u 1 ), Q 1 (u 2 )) P 1 (Q 1 (u 1 )) • P 1 (Q 1 (u 2 )) • du 1 • du 2 ⎛ ⎝ x = Q 1 (u 1 ), ỹ = Q 1 (u 2 ) ⇒ d x = 1 P 1 (Q 1 (u 1 )) • du 1 , d ỹ = 1 P 1 (Q 1 (u 2 )) • du 2 ⎞ ⎠ = C(F 1 ( xh ), F 1 ( ỹh )),
where C ∈ [0, 1] is the copula CDF [Sch07] associated to the microfacet slope distribution. We can express the PDF of the microfacet slope distribution as the derivative of the copula

P( xh , ỹh ) = ∂ 2 C(F 1 ( xh ), F 1 ( ỹh )) ∂ xh • ∂ ỹh = c(F 1 ( xh ), F 1 ( ỹh )) • P 1 ( xh ) • P 1 ( ỹh ),
where c = ∂ 2 C(u 1 ,u 2 ) ∂u 1 •∂u 2 ≥ 0 is the copula PDF associated to the microfacet slope distribution. It follows that we can write

F 3 ( ỹh | xh ) = ỹh -∞ P( xh , ỹ) P 1 ( xh ) • d ỹ = ỹh -∞ c(F 1 ( xh ), F 1 ( ỹ)) • P 1 ( xh ) • P 1 ( ỹ) P 1 ( xh ) • d ỹ = F 1 ( ỹh ) 0 c(F 1 ( xh ), u) • P 1 (Q 1 (u)) P 1 (Q 1 (u)) • du ⎛ ⎝ ỹ = Q 1 (u) ⇒ d ỹ = 1 P 1 (Q 1 (u)) • du ⎞ ⎠ = ∂C(F 1 ( xh ), F 1 ( ỹh )) ∂u 1 . ∂C(F 1 ( xh ), 0) ∂u 1 = 0
Hence, we can compute and store the partial derivative of the copula in practice, which avoids numerical issues since no divisions are involved. The quantile function Q 3 is computed and stored similarly to Q 1 for multiple xh values.

Chapter 5

Conclusion

(1) Subpixel surface.

(2) Microsurface of the subpixel BRDF.

(3) Combined subpixel surfaces.

subpixel subpixel subpixel (3) We show how to combine the microsurface and the displaced surface together in order to produce a noncentral microsurface, which (4) can be handled with noncentral microfacet theory. (5) The filtered response of a displaced surface can thus be computed from a multiscale microfacet BRDF.

Downsampling Operator

In the previous chapters, we saw how to filter displaced Fresnel mirrors at any scale with noncentral microfacet theory, and how to convert any material into a microfacet BRDF with inverted microfacet theory. With both these tools at hand, we can formulate a more general filtering scheme for displaced surfaces of arbitrary reflectance as the problem of deriving a microfacet BRDF of a surface displaced by the sum of two random processes: one that comes from the displacement map itself, and another that comes from the base material of the surface (if the base material of the surface is a microfacet BRDF, then the displacement is already known;

otherwise, we convert it using the algorithms we presented in Chapter 4). It follows that if we can retrieve the microfacet NDF that emerges from the combination of both processes, then we solve our filtering problem; the entire idea is illustrated in Figure 5.1.

Combining Displacements By noting that summing displacement processes also implies summing their derivatives, i.e., their slopes, we can express the problem of retrieving the slope PDF (and hence, the NDF) of the combined displacement process more formally as follows: If we denote by X and Y the random slope processes yielded by the displacement map, and the base material with PDFs P X ≥ 0 and P Y ≥ 0 respectively, then what is the PDF P Z ≥ 0 of the random process Z = X + Y? We can solve this problem by using the characteristic function [START_REF] Lukacs | A Survey of the Theory of Characteristic Functions[END_REF][START_REF] Gray | Characteristic Functions in Radar and Sonar[END_REF]. The characteristic function C Z ≥ 0 of a bivariate random process Z is defined as the Fourier transform of its PDF, i.e.,

C Z (t) = R 2 e i t h • P Z ( h) • d h. (5.1)
Equivalently, Equation (5.1) can be seen as a complex moment function

C Z (t) = E[e i Zt ]. (5.2)
Since Z = X + Y, we can also write C Z (t) = E[e i(X+Y)t ], so that the characteristic function of Z becomes a function of the PDFs of X and Y

C Z (t) = R 2 R 2 e i t (x+ỹ) • P X (x) • P Y (ỹ) • d x • d ỹ. (5.
3)

The PDF of Z is then the inverse Fourier transform of Equation ( 5.3), which, after a few manipulations and using the definition of the Dirac delta function, expresses P Z as the convolution P X * P Y between P X and P

Y P Z ( h) = 1 4π 2 R 2 e -i ht •C Z (t) • dt = 1 4π 2 R 2 e -i ht • R 2 R 2 e i(x+ỹ)t • P X (x) • P Y (ỹ) • d x • d ỹ • dt (Equation (5.3)) = 1 4π 2 R 2 R 2 R 2 e i(x+ỹ-h)t • P X (x) • P Y (ỹ) • d x • d ỹ • dt = R 2 R 2 δ(x + ỹ -h) • P X (x) • P Y (ỹ) • d x • d ỹ δ(x) def = 1 4π 2 R 2 e i xt • dt = R 2 P X (x) • P Y ( h -x) • d x R 2 δ(x + ỹ -h) • P Y (ỹ) • d ỹ = P Y ( h -x) = (P X * P Y )( h).
We can thus retrieve the microfacet NDF of the combined process by computing a convolution of microfacet slope PDFs. Note that this property only works if X and Y are statistically independent, which is a valid assumption in the case of merging a displacement map into a microfacet BRDF process.

Beckmann Distributions In a physically based rendering context, numerical convolutions are too expensive to be practical. We can avoid numerical convolutions by relying on Beckmann slope distributions. Indeed, Beckmann slope distributions belong to the family of stable distributions [START_REF] Nolan | Stable Distributions -Models for Heavy Tailed Data[END_REF] (as they are equivalent to Gaussian distributions, up to a scale factor), which have the convenient particularity of being invariant to convolution, i.e., the convolution of two Beckmann PDFs is also a Beckmann PDF. More formally, if X and Y follow independent Beckmann distributions with respective PDF parameters {α x,X , α y,X , ρ X , ñX } and {α x,Y , α x,Y , ρ Y , ñY }, then the random process Z = X + Y follows a Beckmann distribution with parameters {α x , α y , ρ, ñ}, where

α 2 x = α 2 x,X + α 2 x,Y α 2 y = α 2 y,X + α 2 y,Y ρ = ρ X α x,X α y,X + ρ Y α x,Y α y,Y (α 2 x,X + α 2 y,X )(α 2 x,Y + α 2 y,Y ) ñ = ñX + ñY .

Positioning and Discussions

The problem of appearance modeling and filtering is not new in the computer graphics literature. Here, we position and discuss some properties and aspects of the contributions of this thesis with respect to existing solutions.

Computer Graphics BRDFs Reflectance modeling has a long history in computer graphics, and entire surveys are dedicated to this only topic [START_REF] Schlick | A Survey of Shading and Reflectance Models[END_REF][START_REF] Kurt | A Survey of BRDF Models for Computer Graphics[END_REF][START_REF] Montes | An Overview of BRDF Models[END_REF]. According to Blinn [START_REF] Blinn | Models of Light Reflection for Computer Synthesized Pictures[END_REF], the very first computer graphics BRDF models acted as ideal Lambertian reflectors. While matte surfaces tend to be well approximated by such BRDFs, they quickly turned out to be inadequate to mimic glossy apperances. In or- Note that if we employ the normal map GAF in our BRDF evaluations, then our framework acts as a normal map filtering algorithm that supports arbitrary reflectance, thus generalizing some previous work in this area.

LEAN Mapping Olano and Baker [START_REF] Olano | LEAN Mapping[END_REF][START_REF] Baker | Spectacular Specular -LEAN and CLEAN Specular Highlights[END_REF] introduced the lightweight representation we use in our framework in the context of normal map filtering. We adopted their representation for the two following reasons.

First, their representation allows for linear filtering of the data in a manner that properly captures filtered reflectance. Second, their representation supports anisotropy and multilayer superposition. The main limitation of LEAN mapping in the context of our work comes from the non-physically based BRDF employed by the model: It lacks proper normalization, as the Jacobian of the reflection operation is missing, and is not capable of reproducing the important masking, shadowing, and projection weighting effects. Moreover, it is only designed for specular microfacets under point and directional lighting. This thesis showed how to account for occlusion effects, as well as arbitrary lighting conditions.

BRDFs Despite its effectiveness to faithfully describe many real-world appearances, microfacet BRDFs remain applicable to a limited set of materials because effects such as transmission, multibounce scattering, diffraction, or layers are not supported by the extended theory we have presented. While it has been shown that certain effects could be incorporated, e.g., transmission [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF], diffraction [START_REF] Stam | Diffraction Shaders[END_REF], or layers [START_REF] Wenzel | A Comprehensive Framework for Rendering Layered Materials[END_REF],

most of the remaining effects are still open problems in the context of deriving a multiscale representation for matter, and should be addressed in order to make the theory even more general.

Curvature Our displaced surface filtering model and its derivations are based on the fundamental assumption that the displacement map is applied over a locally planar patch. In theory, this means that the method is not valid when the pixel footprint covers a curved macrosurface since the curvature must be filtered along with the displacement map. For instance, the claws of the T-rex model in Figure 3.11 are small, smooth, curved, and highly specular. As such, they exhibit aliasing that cannot be filtered by our scheme alone. While, in practice, we could deal with this problem by offsetting roughness heuristically due to curvature [DHI + 13], the problem of properly filtering displacement maps along with large pieces of macrosurfaces is complex and remains an open problem [START_REF] Bruneton | A Survey of Non-linear Pre-filtering Methods for Efficient and Accurate Surface Shading[END_REF].

Multiscale Representation for Matter Our filtering framework makes the per-pixel rendering complexity proportional to the number of displaced surfaces that project into a pixel, where previous algorithms would have been proportional to the total number of texels that project into a pixel. While this shift in complexity is undoubtedly a step forward towards more reasonable rendering times, the ability to handle multiple surfaces inside a pixel remains an open problem. Such scenarios arise for, e.g., trees and fur, whose geometry clearly can not be solely described with displaced surfaces at all scales. A first attempt to deal with such configurations was introduced concurrently to this thesis in the paper entitled "The SGGX Microflake Distribution" [START_REF] Heitz | The SGGX Microflake Distribution[END_REF],

which addresses the problem of downsampling directional data represented as volumes (rather than surfaces), using voxels. The paper shows how the GGX distribution can be extended to work in the spherical (rather than hemispherical) domain, and plugged into a microflake phase function (the 3D analogue of a microfacet BRDF [JAM + 10]), to derive an efficient downsampling operator for 3D directional data. While downsampling voxelized directional data is an important contribution, it only constitutes a small brick in the elaboration of a more general multiscale representation for matter, which should be capable of accounting for correlations between the way matter is oriented, structured, and how it interacts with light. Fulfilling these goals requires to extend our current understanding of how light and matter interact, which should be worth the effort.

  part of the fantastic physically based shading course you have been organizing along with Stephen McAuley so tremendously well for several years in a row now. Antoine Bouthors: for inviting me over at Weta Digital, and giving the chance to integrate some of my research in your rendering pipeline. Wenzel Jakob: for providing the amazing Mistuba renderer; the figures and source code snippets I included in my manuscript are inspired from the documentation you wrote for Mitsuba. Fabrice Neyret and Eric Bruneton: for sharing with me your philosophy of level-of-detail and scientific methodology during my master thesis. Brent Burley: for providing me with the T-Rex model and sharing your insights on BRDF authoring. Thiago Da Costa and Arno Zinke: for providing me with the LAGOA robot asset. Cyril Crassin: for perhaps being the only one I can talk to about

Figure 1 . 1 -Figure 1

 111 Figure 1.1 -Shape versus appearance. Left: Objects with different shapes but the same appearance. Right: Objects with the same shape but different appearance. The plasma ball image was created by Soenke Rahn and is available at http://commons.wikimedia.org/wiki/File:Plasmakugel_%28Plasma_Ball%29.JPG

Figure 1 . 4 -

 14 Figure 1.4 -Theoretical background illustrated.

  double nint = 0; const int cnt = 100; for (int i = 0; i < cnt; ++i) { double u = (double)i / (double)cnt; double theta = u * M_PI * 0.5; double x = tan(theta); double cos_term = cos(theta); nint+= my_func(x) / (cos_term * cos_term); } nint*= M_PI * 0.5 / (double)cnt;

Figure 2 5 Figure 2 . 5 -

 2525 Figure 2.2 -Varying the Fresnel term on a microfacet BRDF.

  vec3 microfacet::evalp(const vec3& i, const vec3& o, const void *user_param) const { const microfacet::params params = user_param ? *reinterpret_cast<const microfacet::params *>(user_param) : microfacet::params::standard(); vec3 h = normalize(i + o); float_t G = gaf(h, i, o, params); if (G > 0.0) { float_t cos_theta_d = dot(o, h); vec3 F = fresnel(cos_theta_d); float_t D = ndf(h, params); return (F * D * G) / (4.0 * o.z); } return vec3(0); } vec3 microfacet::eval(const vec3& i, const vec3& o, const void *user_param) const { return (evalp(i, o, user_param) / i.z); } Listing 2.2 -Evaluating a microfacet BRDF in C++.

Figure 2 . 8 -

 28 Figure 2.8 -Typical surface profiles handled by existing microfacet GAFs. (1, 4) The Smith GAF assumes that the microsurface structure is irrelevant for occlusion effects. (2, 5) The ∨-groove GAF assumes that the microsurface consists exclusively of symmetric ∨-groove cavities. (3, 6) The normal map GAF neglects all occlusion effects.

  ) float_t microfacet::vndf( const vec3& h, const vec3& k, const microfacet::params& params ) const { float_t D = ndf(h, params); float_t G1 = g1(h, k, params); float_t kh = max((float_t)0.0, dot(k, h)); return (kh * D * G1 / k.z); } Listing 2.6 -Evaluating the distribution of visible normals in C++.

  float_t microfacet::ndf(const vec3& h, const microfacet::params& params) const { float_t cos_theta_h_sqr = h.z * h.z; float_t cos_theta_h_sqr_sqr = cos_theta_h_sqr * cos_theta_h_sqr; float_t xslope = -h.x / h.z; float_t yslope = -h.y / h.z; return (p22(xslope, yslope, params) / cos_theta_h_sqr_sqr); } float_t microfacet::p22(float_t x, float_t y, const params& params) const { float_t nrm = params.m_a1 * params.m_a2; float_t x_ = x / params.m_ax; float_t tmp1 = params.m_ax * y -params.m_rho * params.m_ay * x; float_t tmp2 = params.m_ax * params.m_ay * params.m_rho; float_t y_ = tmp1 / tmp2; return p22_std(x_, y_) / nrm; } Listing 2.7 -Evaluating the distribution of visible normals in C++.

  float_t microfacet::g1_smith( const vec3& h, const vec3& k, const microfacet::params& params ) const { float_t a = k.x * params.m_ax + k.y * params.m_ay * params.m_rho; float_t b = k.y * params.m_ay * params.m_sqrt_one_minus_rho_sqr; vec3 kprime = normalize(vec3(a, b, k.z)); return g1_smith_std(h, kprime); } Listing 2.8 -Evaluating the Smith masking term with arbitrary roughness parameters in C++.

  .38) float_t radial::p22_std(float_t x, float_t y) const { return p22_radial(sqrt(x * x + y * y)); } Listing 2.11 -C++ implementation for radial microfacet slope distributions.

  and z k = cos θ k , we factorize the denominator by sin θ k , which yields

  and going back to slope space with xm = r h • cos φ m and ỹm = r h • sin φ m

  bool test_ndf(const djb::microfacet& brdf) { int N_i = 270; int N_j = 360; double nint = 0.0; for (int i = 0; i < N_i; ++i) { double u_i = (double)i / (double)N_i; double theta_h = u_i * u_i * M_PI / 2.0; for (int j = 0; j < N_j; ++j) { double u_j = (double)j / (double)N_j; double phi_h = u_j * 2.0 * M_PI; djb::vec3 h = djb::vec3(theta_h, phi_h); nint+= cos(theta_h) * brdf.ndf(h) * sin(theta_h) * u_i; } } nint*= 2.0 * M_PI * M_PI / (double)(N_i * N_j); return (fabs(nint -1.0) < EPSILON); } Listing 2.13 -C++ microfacet NDF validation test. bool test_vndf(const djb::microfacet& brdf, const djb::vec3& k) { int N_i = 270; int N_j = 360; double nint = 0.0; for (int i = 0; i < N_i; ++i) { double u_i = (double)i / (double)N_i; double theta_h = u_i * u_i * M_PI / 2.0; for (int j = 0; j < N_j; ++j) { double u_j = (double)j / (double)N_j; double phi_h = u_j * 2.0 * M_PI; djb::vec3 h = djb::vec3(theta_h, phi_h); nint+= brdf.vndf(h, k) * sin(theta_h) * u_i; } } nint*= 2.0 * M_PI * M_PI / (double)(N_i * N_j); return (fabs(nint -1.0) < EPSILON); } Listing 2.14 -C++ microfacet visible NDF validation test.

Figure 3 . 1 -

 31 Figure3.1 -Geometric overview of noncentral microfacet theory. Can microfacet theory be made applicable to multiple scales? (1) In order to derive the microfacet BRDF of a microsurface at any scale, microfacet theory must be extended to account for (2) configurations where the mean surface differs from the tangent plane. (3) By shearing such surfaces back to a standard configuration, we show that (4) the BRDF of such surfaces can be determined from the standard microsurface by(5, 6) transforming the illumination and observation directions.

3. 1

 1 MotivationThis chapter encompasses the main contributions of the LEADR mapping paper [DHI + 13], which addresses the problem of accelerating displacement texture map filtering with MIP mapping[START_REF] Williams | Pyramidal Parametrics[END_REF]. A displacement map[START_REF] Robert | Shade trees[END_REF][START_REF] Cook | The Reyes Image Rendering Architecture[END_REF] is a texture that displaces the surface it is mapped onto. In a physically based rendering context, the elements (texels) of such textures typically combine a scalar value and a normal. The former gives the magnitude of the displacement along the surface normal, while the latter provides the final normal of the surface, which should be used for shading. Although historically reserved to offline rendering contexts, displacement maps have progressively become ubiquitous in real-time rendering applications thanks to the availability of hardware accelerated subdivision surfaces[START_REF] Niessner | Feature-adaptive GPU Rendering of Catmull-Clark Subdivision Surfaces[END_REF][START_REF] Niessner | Analytic Displacement Mapping Using Hardware Tessellation[END_REF]. As such, they are now a fundamental tool for rendering.

  float_t microfacet::ndf(const vec3& h, const microfacet::params& params) const { float_t cos_theta_h_sqr = h.z * h.z; float_t cos_theta_h_sqr_sqr = cos_theta_h_sqr * cos_theta_h_sqr; float_t xslope = -h.x / h.z -params.m_tx_n; float_t yslope = -h.y / h.z -params.m_ty_n; return (p22(xslope, yslope, params) / cos_theta_h_sqr_sqr); } Listing 3.1 -Evaluating the NDF term with arbitrary roughness and shear parameters in C++.

  float_t microfacet::g1_smith( const vec3& h, const vec3& k, const microfacet::params& params ) const { float_t a = k.x * params.m_ax + k.y * params.m_ay * params.m_rho; float_t b = k.y * params.m_ay * params.m_sqrt_one_minus_rho_sqr; float_t c = k.z -k.x * params.m_tx_n -k.y * params.m_ty_n; vec3 kprime = normalize(vec3(a, b, c)); float_t g1 = g1_smith_std(h, kprime); if (g1 > 0.0) return g1 * (k.z / c); else return 0.0; } Listing 3.2 -Evaluating the Smith masking term with arbitrary roughness and shear parameters in C++.

Figure 3 . 3 -

 33 Figure 3.3 -Shear invariance of the monostatic shadowing term.

  Olano and Baker[START_REF] Olano | LEAN Mapping[END_REF] under the name of LEAN mapping. It requires five floating-point numbers per MIP texel and is capable of capturing the roughness and main direction of anisotropy of a set of slopes, which is usually sufficient for most filtering cases. // Linear Representation class beckmann::lrep { friend class beckmann; public: // Ctor lrep(float_t E1 = 0.0, float_t E2 = 0.0, float_t E3 = 0.5, float_t E4 = 0.5, float_t E5 = 0.0); private: // members float_t m_E1, m_E2; // first order slope moments float_t m_E3, m_E4; // second order slope moments float_t m_E5; // first order joint slope moment }; Listing 3.4 -C++ implementation of a linear representation for Beckmann microfacet distributions.

  void beckmann::params_to_lrep(const microfacet::params& params, lrep *lrep) { float_t ax, ay, rho, tx_n, ty_n; params.get_pdfparams(&ax, &ay, &rho, &tx_n, &ty_n); (*lrep) = beckmann::lrep(tx_n, ty_n, 0.5 * ax * ax + tx_n * tx_n, 0.5 * ay * ay + ty_n * ty_n, 0.5 * rho * ax * ay + tx_n * ty_n); } void beckmann::lrep_to_params(const lrep& lrep, microfacet::params *params) { float_t tx_n = lrep.m_E1; float_t ty_n = lrep.m_E2; float_t tmp1 = max((float_t)0.0, lrep.m_E3 -lrep.m_E1 * lrep.m_E1); float_t tmp2 = max((float_t)0.0, lrep.m_E4 -lrep.m_E2 * lrep.m_E2); /* clamp the parameters to valid values */ float_t ax = max(1e-5, sqrt(2.0 * tmp1)); float_t ay = max(1e-5, sqrt(2.0 * tmp2)); float_t rho = 2.0 * (lrep.m_E5 -lrep.m_E1 * lrep.m_E2) / (ax * ay); rho = min((float_t)0.99, max((float_t)-0.99, rho)); (*params) = microfacet::params::pdfparams(ax, ay, rho, tx_n, ty_n); } Listing 3.5 -C++ conversion routines for Beckmann microfacet distributions.

Figure 3 . 5 -

 35 Figure 3.5 -Experimental validation of our filtering scheme. We compare the appearance of supersampled (256 samples per pixel) Beckmann microgeometries against the response of our multiscale microfacet BRDF for several roughnesses and lighting conditions. The Microsurface column is a preview of the Beckmann microgeometry.

Figure 3 . 8 -

 38 Figure 3.8 -Surface shearing.

Figure 3 .

 3 Figure 3.10 -A displacement mapped ocean rendered at different scales in real time with LEADR mapping.

Figure 3 .

 3 Figure 3.11 -A production asset (Ptex T-rex model c Walt Disney Animation Studios.) rendered in real time with LEADR mapping.

Figure 4 . 1 -

 41 Figure 4.1 -Geometric overview of inverted microfacet theory. A microfacet BRDF extracts a BRDF from a microsurface.(1) Can we solve the inverse problem, i.e., is it possible to extract a microsurface from a BRDF? We can tackle this problem by considering the BRDF in (2, 3) backscattering configurations (o = i = h). For such configurations, occlusions due to (4) masking and (5) shadowing become equal, and (6) the evaluation of the BRDF is directly proportional to the fraction (with respect to the entire microsurface) of microfacets that are visible from the backscattering direction. This observation is the key to solve the inverse problem.

  .14) void tabular::compute_p22_smith(const brdf& brdf, int res) { int cnt = res -1; float_t dtheta = sqrt(M_PI * 0.5) / (float_t)cnt; matrix km(cnt); for (int i = 0; i < cnt; ++i) { float_t tmp = (float_t)i / (float_t)cnt; float_t theta = tmp * sqrt(M_PI * 0.5); float_t theta_o = theta * theta; float_t cos_theta_o = cos(theta_o); float_t tan_theta_o = tan(theta_o); vec3 fr = brdf.eval(vec3(theta_o, 0.0), vec3(theta_o, 0.0)); float_t fr_i = fr.intensity(); float_t kji_tmp = (dtheta * pow(cos_theta_o, 6.0)) * (8.0 * fr_i); for (int j = 0; j < cnt; ++j) { const float_t dphi_h = M_PI / 180.0; float_t tmp = (float_t)j / (float_t)cnt; float_t theta = tmp * sqrt(M_PI * 0.5); float_t theta_h = theta * theta; float_t cos_theta_h = cos(theta_h); float_t tan_theta_h = tan(theta_h); float_t tan_product = tan_theta_h * tan_theta_o; float_t nint = 0.0; for (float_t phi_h = 0.0; phi_h < 2.0 * M_PI; phi_h+= dphi_h) nint+= max(1.0, tan_product * cos(phi_h)); nint*= dphi_h; km(j, i) = theta * kji_tmp * nint * tan_theta_h / (cos_theta_h * cos_theta_h); } } // compute slope pdf const std::vector<double> v = km.eigenvector(4); for (int i = 0; i < (int)v.size(); ++i) m_p22.push_back(1e-2 * v[i]); m_p22.push_back(0); } Listing 4.1 -Extracting a microfacet slope PDF assuming a Smith GAF in C++.

  microfacet::params tabular::fit_beckmann_parameters(const tabular& tab) { const int ntheta = 128; float_t dtheta = M_PI / (float_t)ntheta; float_t nint = 0.0; float_t alpha; for (int i = 0; i < ntheta; ++i) { float_t u = (float_t)i / (float_t)ntheta; // in [0,1) float_t theta_h = u * u * M_PI * 0.5; float_t cos_theta_h = cos(theta_h); float_t r_h = tan(theta_h); float_t r_h_sqr = r_h * r_h; float_t p22_r = tab.p22_radial(r_h_sqr); nint+= (u * r_h_sqr * r_h * p22_r) / (cos_theta_h * cos_theta_h); } nint*= dtheta * M_PI; /* M_PI = int_0^2pi cos^2 phi dphi */ alpha = sqrt(2.0 * nint); return microfacet::params::isotropic(alpha); } Listing 4.4 -Converting a radial tabulated slope PDF into a Beckmann roughness parameter in C++.

Figure 4

 4 Figure 4.2.

Figure 4 . 2 -

 42 Figure 4.2 -Mean delta-E difference image on the entire MERL database [MPBM03].

Figure 4 . 4 .

 44 Figure 4.4. For each material, the lobe of the BRDF (and hence the microfacet NDF) is captured accurately.

Figure 4 . 5 ,

 45 where the timings include the computations of the slope PDF, the Smith term, the quantile functions, and the SGD

Figure 4 .Figure 4 . 4 -

 444 Figure 4.3 -Side-by-side fitting comparisons of a few representative isotropic materials from the MERL database [MPBM03].

Figure 4 . 5 -

 45 Figure 4.5 -Fitting timings (in seconds) of our algorithms as a function of the number of input material evaluations along the elevation direction.

Figure 4 .

 4 Figure 4.2), which is in agreement with previous observations [TR75, WMLT07, Bur12]. Note that the analytic models are also incorporated in the detailed fitting analysis provided in the supplemental document of the original article [DHI + 15].

Figure 4 . 6 -

 46 Figure 4.6 -Multiple material design on a production asset. (a) Measured BRDFs. (b) Our fitted microfacet BRDFs. (c) Our fitted microfacet BRDFs controled by a roughness texture map. Model courtesy of LAGOA.

  void tabular::normalize_p22() { const int ntheta = 128; const float_t dphi = 2.0 * M_PI; const float_t dtheta = M_PI / (float_t)ntheta; float_t k, nint = 0.0; for (int i = 0; i < ntheta; ++i) { float_t u = (float_t)i / (float_t)ntheta; // in [0,1) float_t theta_h = u * u * M_PI * 0.5; float_t r_h = tan(theta_h); float_t cos_theta_h = cos(theta_h); float_t p22_r = p22_radial(r_h * r_h); nint+= (u * p22_r * r_h) / (cos_theta_h * cos_theta_h); } nint*= dtheta * dphi; // normalize the slope pdf k = 1.0 / nint; for (int i = 0; i < (int)m_p22.size(); ++i) m_p22[i]*= k; } Listing 4.6 -Slope PDF normalization in C++.

( 4 )Figure 5 . 1 -

 451 Figure 5.1 -Geometric overview of our displacement map filtering algorithm. In order to handle (1) displaced surfaces of arbitrary reflectance, we first (2) extract a microsurface from the BRDF of the base surface. (3) We show how to combine the microsurface and the displaced surface together in order to produce a noncentral microsurface, which (4) can be handled with noncentral microfacet theory.(5) The filtered response of a displaced surface can thus be computed from a multiscale microfacet BRDF.

Figure 5

 5 Figure 5.2 -GPU pipeline of our displacement map filtering technique.

Figure 5 .Figure 5 . 3 -

 553 Figure 5.3 -Comparison between (top) naive MIP mapping and (bottom) our downsampling operator. The base material is nickel from the MERL database.

Figure 5 . 4 -Figure 5 . 5 -

 5455 Figure 5.4 -Effect of our downsampling operator on scratch-like displaced surfaces with a few representative isotropic materials from the MERL database [MPBM03].

Table 2 .

 2 I -Popular Fresnel terms used in computer graphics.

	/* Fresnel API */
	namespace fresnel {
	/* Interface */
	class impl {
	public:
	virtual ~impl() {}
	virtual vec3 eval(float_t cos_theta_d) const = 0;
	virtual impl *copy() const = 0;
	};
	} // namespace fresnel
	Listing 2.3 -Fresnel interface in C++.

Table 2

 2

	.II. Figures 2.2 and 2.4 show materials generated with different D terms, while F and G are
	fixed. Note that the GGX distribution is known as the Trowbridge-Reitz model in the physics literature [TR75],

Table 2 .

 2 II -Analytic microfacet slope probability density functions.Listing 2.5 -Evaluating the Smith GAF term in C++.

	float_t
	microfacet::gaf_smith(
	const vec3& h, const vec3& i, const vec3& o,
	const microfacet::params& params
	) const {
	float_t G1_o = g1(h, o, params);
	float_t G1_i = g1(h, i, params);
	float_t tmp = G1_i * G1_o;
	if (tmp > 0.0)
	return (tmp / (G1_i + G1_o -tmp));
	return 0.0; // fully attenuated
	}

  2.11 -Sampling a rough microsurface from a standard microsurface.

	2.4 Implementation Details
	/* Radial microfacets */
	class radial : public microfacet {
	public:
	radial(int gaf = microfacet::GAF_SMITH,
	const fresnel::impl& fresnel = fresnel::ideal(),
	bool shadow = true):
	microfacet(gaf, fresnel, shadow) {}
	// queries
	virtual float_t g1_radial(float_t cot_theta_k) const = 0;
	virtual float_t p22_radial(float_t r_sqr) const = 0;
	private:
	float_t p22_std(float_t x, float_t y) const;
	float_t g1_smith_std(const vec3& k) const;
	};
	vec3
	microfacet::sample(
	float_t u1, float_t u2,
	const vec3& o,
	const void *user_param
	) const {
	// warp the variate with the microfacet parameters
	tx_h = params.m_ax * tx_m;
	float_t choleski = params.m_rho * tx_m
	+ params.m_sqrt_one_minus_rho_sqr * ty_m;
	ty_h = params.m_ay * choleski;
	// return the associated normal
	return normalize(vec3(-tx_h, -ty_h, 1));
	}
	Listing 2.9 -Importance sampling a microfacet BRDF with arbitrary roughness parameters in C++.

const microfacet::params params = user_param ? *reinterpret_cast<const microfacet::params *>(user_param) : microfacet::params::standard(); float_t tx_h, ty_h, tx_m, ty_m; // create a standard variate switch (m_gaf) { case GAF_NONE: sample_nmap(u1, u2, o, &tx_m, &ty_m); break; case GAF_VGROOVE: sample_vgroove(u1, u2, o, &tx_m, &ty_m); break; case GAF_SMITH: sample_smith(u1, u2, o, &tx_m, &ty_m); break; default: abort(); // should never happen }

  .4.

	vec3
	microfacet::sample(
	float_t u1, float_t u2,
	const vec3& o,
	const void *user_param
	) const {
	const microfacet::params params =
	user_param ? *reinterpret_cast<const microfacet::params *>(user_param)
	: microfacet::params::standard();
	float_t tx_h, ty_h, tx_m, ty_m;
	// create a standard variate
	switch (m_gaf) {
	case GAF_NONE: sample_nmap(u1, u2, o, &tx_h, &ty_h); break;
	case GAF_VGROOVE: sample_vgroove(u1, u2, o, &tx_h, &ty_h); break;
	case GAF_SMITH: sample_smith(u1, u2, o, &tx_h, &ty_h); break;
	default: abort(); // should never happen
	}
	// warp the variate with the microfacet parameters
	tx_m = params.m_ax * tx_h + params.m_tx_n;
	float_t choleski = params.m_rho * tx_h
	+ params.m_sqrt_one_minus_rho_sqr * ty_h;
	ty_m = params.m_ay * choleski + params.m_ty_n;

  The slopes in world space can be defined by If the distortions occurring on the macrosurface have low spatial variations, then we can reasonably assume that the mapping distortions represented by the terms ∂u ∂x , ∂v ∂x , ∂u ∂y , and ∂v ∂y are locally constant. It follows that the appropriate slope distribution in world space gives

	∂ξ ∂x	,	∂ξ ∂y	=	∂ξ ∂u	∂u ∂x	+	∂ξ ∂v	∂v ∂x	,	∂ξ ∂u	∂u ∂y	+	∂ξ ∂v	∂v ∂y	.
	E 1 =	∂u ∂x	E 1 +		∂v ∂x	E 2
	E 2 =	∂u ∂y	E 1 +		∂v ∂y	E 2
	E 3 =		∂u ∂x	2	E 3 +		∂v ∂x	2	E 4 + 2	∂u ∂x	∂v ∂x	E 5
	E 4 =		∂u ∂y	2	E 3 +		∂v ∂y	2	E 4 + 2	∂u ∂y	∂v ∂y	E 5
	E 5 =	∂u ∂x	∂u ∂y	E 3 +	∂v ∂x	∂v ∂y	E 4 +	∂u ∂x	∂v ∂y	+	∂v ∂x	∂u ∂y	E 5 .

  .17) Conversion to Beckmann Roughness We retrieve the parameters by computing 2nd order moments like in LEAN/LEADR mapping [DHI + 13]

Acknowledgements

This thesis concludes four invaluable years of research; I have learned so much and met so many great people. I first wish to thank my advisors Jean-Claude Iehl, Pierre Poulin, and Victor Ostromoukhov, who gave me the

your continuous support throughout the years has been invaluable. Financial Support Financial support was partly provided by grants from GRAND and the NSERC Discovery of Pierre Poulin, and the ANR AMCQMCSCGA of Victor Ostromoukhov. I also acknowledge an Explo'ra Doc grant from the Rhône-Alpes region.

leanmap_1(i, j, 0, 0) = slope_x; // E1 leanmap_1(i, j, 0, 1) = slope_y; // E2 leanmap_1(i, j, 0, 2) = 1.f; leanmap_1(i, j, 0, 3) = 1.f; leanmap_2(i, j, 0, 0) = slope_x_sqr + base_roughness_sqr; // E3 leanmap_2(i, j, 0, 1) = slope_y_sqr + base_roughness_sqr; // E4 leanmap_2(i, j, 0, 2) = slope_xy; // E5 leanmap_2(i, j, 0, 3) = 1.f; } } Listing 3.7 -C++ implementation of a routine to compute a LEAN/LEADR map from a normal map.

Surface Stretching and Shearing

In practice, it is more convenient to store the moments of the unit parametric space slopes. We denote these moments as E 1 , E 2 , E 3 , E 4 , and E 5 . The distribution must be expressed in the macrosurface tangent frame (x, y, z). Since meshes are rarely perfectly parameterized, the texture parameterization u(x, y), v(x, y) usually produces distortions, which can be further accentuated in the case of animations. To arrive at this particular result, we start from Equation (4.5) and apply Equation (4.12)

We then proceed with a few simple algebric manipulations, and Equation (4.13) naturally emerges

Note that the choice of the azimuthal angle φ o to fully define o is arbitrary in Equation (4.14). As for the general case, Equation (4.13) may be expressed as an eigenvalue problem of the form

where p = (g(

Since the entries of matrix K are all nonnegative, we can also solve Equation (4.15) with the power iteration method. In practice, we also use four successive multiplications to recover the solution p . Algorithm 3 provides pseudocode for this specialized method.

The choice is arbitrary here Inversion with the Normal Map GAF We give here the procedure to extract a microfacet slope PDF assuming a normal map GAF. In this case, the inversion procedure is even simpler than for the ∨-groove GAF, and only requires a few algebric manipulations Note that if P is a GGX slope PDF, then our conversion is exact.

Experiments

A Note on Precomputations In our implementation, we rely exclusively on rectangle quadrature rules to compute our integrals. This includes the microfacet slope PDF, i.e., w 1 = • • • = w N = 1/N. We store one RGB 1D table for the microfacet Fresnel, and three scalar-valued 1D (resp. 2D) tables for each of P, G 1 , and the quantile function for isotropic (resp. anisotropic) materials. The dimensions of the tables depend on the number of samples that we evaluate from the material in the elevation and azimuthal directions. In the isotropic configuration, the number of samples is equal to N. In this case, our representation stores (including Fresnel and the other tables) 3N + 3N = 6N scalar values. In the anisotropic configuration, the number of samples is equal to N = N θ × N φ , where N θ and N φ respectively denote the number of elevation and azimuthal samples used to evaluate the reflectance of the input material. In this case, our representation stores (including Fresnel and the other tables) 3N θ × (1 + N φ ) scalar values.

Unit Testing We can validate our microfacet slope PDF extraction by testing our algorithms against an analytic model of Equation (4.3) based on a Beckmann distribution [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF] with varying roughness. Table 4.I shows the maximum relative errors that we measured during our experiment in the isotropic case, using N = 360 BRDF evaluations. As can be seen from the reported numbers, the error is small (less than 3%). We attribute Importance Sampling Table Our tabulated representation does not use the distribution of visible normals for importance sampling. Rather, it relies on the less effective scheme of Walter et al. [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF], which consists in sampling the distribution of slopes. We describe here the procedure to build the table in the case of radial slope distributions. In this case, we produce a sample in polar coordinates. The angle is taken uniformly, while the radius is sampled from the quantile function associated to the radial component of the distribution P r . First, we compute the radial CDF F r ∈ [0, 1] associated to P r F r (r h ) = 2π

Once the CDF is computed, we invert it numerically to get the quantile function Q r = F -1 r . In the following paragraphs, we describe a way to precompute the tables for importance sampling with a tabulated distribution of visible normals.

Marginal Slope Quantile Function Before building the quantile functions for importance sampling with the distribution of visible slopes, we need to compute the quantile function associated to the marginal slope PDF P 1 > 0, which is

Note that the property P(x, y) = P(x, -y) is only true for isotropic materials, as P 1 is even, i.e., P 1 (x) = P 1 (-x).

We take advantage of this property to precompute the PDF for x ∈ [0, +∞), rather than x ∈ (-∞, +∞). The CDF These parameters can be computed from the linear representation of Beckmann distributions, since we have [OB10]

which implies that the convolution of Beckmann PDFs can be computed by simply summing their linear representation. Implementation Details The implementation of our displacement map downsampling operator is divided into two stages: a precomputation stage and a runtime stage. The precomputation stage consists in converting the input displacement texture maps into LEAN/LEADR maps using Listing 3.6 and Listing 3.7, and computing the associated MIP maps. Furthermore, if the input BRDF of the displaced surface is not based on a Beckmann microfacet BRDF, then it must also be converted into Beckmann roughness parameters using the procedure we described in Listing 4.4. This completes the precomputation stage. At runtime, we incorporate the displacements that are too fine to be perceived at the pixel level into shading by evaluating a noncentral microfacet BRDF, i.e., the equations we introduced in Chapter 3. In order to retrieve the parameters of the BRDF, we evaluate a MIP texture lookup of the LEAN/LEADR map and merge the retrieved linear data with the base roughness of the surface BRDF using Listing 5.1 and Listing 3.5. Figure 5.2 illustrates the GPU pipeline of our implementation.

Multilobe Representations Some previous work uses multiple lobes to represent and store normal distributions [Fou92, TLQ + 05, HSRG07, TLQ + 08]. While multiple lobes can represent more complex distributions than what we can, their precomputations are problematic. Indeed, in order to make interpolation possible, every lobe in a texel must match the same lobe in the neighboring texels [TLQ + 05, HSRG07, TLQ + 08]. Solving this problem requires heavy nonlinear optimizations, and matching failures may result in visual artifacts when two nonmatching lobes are interpolated. In contrast, our single anisotropic lobe representation allows for lightweight memory storage, simple and fast precomputations, as well as no-failure interpolation.

Future Research

The main purpose of this thesis was to introduce a fully functional filtering framework for arbitrary complex displacement mapped surfaces. Our framework builds upon microfacet theory, which we augmented in the previous chapters, as well as upon the convenient properties of Beckmann distributions. While the presented results can offer significant improvements in terms of rendering quality and speed, there is still much room for improvement; the following paragraphs are dedicated to identifying current limitations and avenues for future work.

Beckmann Distributions Despite their very practical properties (linear representation, stability under convolution), Beckmann distributions are probably not the best parametric model to use for filtering for two reasons. First, their short (Gaussian) tails make them highly sensitive to outliers, which tends to yield overblurry highlights as filtering footprints increase. Second, they may fit measured data quite poorly, which introduces significant bias when fitting some materials. This effect is particularly visible for the gold-metallic-paint2 and violet-acrylic materials in Figures 5.4 and 5.5. In order to improve on this second point, a mixture of Beckmann distributions can be employed; an extensive filtering comparison between a single Beckmann lobe and a mixture of eight lobes is available in the supplemental document of this thesis. Note however that mixture models require more arithmetic operations, as well as an additional random number to determine the lobe to sample.

In contrast, GGX distributions [START_REF] Walter | Microfacet Models for Refraction Through Rough Surfaces[END_REF] seem much more promising, because they fit measured data much better, and their heavy tails should (note that this is a conjecture) result in sharper highlights in the presence of outliers. Deriving a complete filtering framework for displaced surfaces of arbitrary reflectance using GGX distributions remains an open problem.