Jonathan Dupuy

Eric Heitz

Jean-Claude Iehl

Pierre Poulin

Fabrice Neyret

Victor Ostromoukhov
email: victor.ostromoukhov@liris.cnrs.freric.heitz|fabrice.neyret@inria.frpoulin@iro.umontreal.ca

To appear in ACM TOG 32(6

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture BRDF, microfacet, filtering, LEAN mapping, GPU

Figure 1: A high-quality animated production model (Ptex T-rex model c Walt Disney Animation Studios.

) rendered in real time under directional and environment lighting using LEADR mapping on an NVidia GTX 480 GPU. The surface appearance is preserved at all scales, using a single shading sample per pixel. Combined with adaptive GPU tessellation, our method provides the fastest, seamless, and antialiased progressive representation for displaced surfaces.

Introduction

Rendering applications such as video games commonly employ bump or normal textures (henceforth interchangeably referred to as normal textures) to enhance surface appearance (e.g., [START_REF] Kilgard | A practical and robust bump-mapping technique for todays GPUs[END_REF]). These textures perturb or modify the normal of a simple underlying surface to emulate geometric variations through shading perturbations. Similarly to albedo textures, normal textures must be filtered for antialiasing purposes. But since shading with linearly filtered normals does not result in proper reflectance filtering, these textures cannot exclusively rely on simple methods such as mipmapping. Recently introduced by [START_REF] Olano | LEAN mapping[END_REF], LEAN mapping is an elegant solution to this problem, that has found widespread adoption because of its effectiveness, efficiency, and simplicity [START_REF] Baker | Spectacular specular -LEAN and CLEAN specular highlights[END_REF]].

Normal mapping is an inherited paradigm from the 1980's. At that time, geometric models were coarse because of computing capabilities and memory constraints. Textures cheaply enhanced visual details without increasing geometric complexity. The discrepancy between the resolutions of geometry and texture was the key assumption for texture filtering: within the same large and flat triangle, visibility, curvature, and orientation were considered constant. Consequently, filtering in texture space or in geometry space could reasonably be considered equivalent. This assumption does not hold any more because mesh and texture resolutions can be matched with negligible overhead on modern GPUs.

Filtering the appearance of small-scale geometry is thus a critical emerging problem. Indeed, small-scale geometry produces viewand light-dependent effects that include masking, shadowing, and projection weighting (i.e., the cosine term). Filtering this smallscale geometry while neglecting these visual effects violates energy conservation and can result in objectionable aliasing, popping artifacts, and inconsistent appearance throughout scales. Methods for filtering normal maps do not account for these effects, which is why filtering reflectance from normal mapping is not the same problem as filtering reflectance of small-scale geometry [START_REF] Han | Frequency domain normal map filtering[END_REF][START_REF] Bruneton | A survey of non-linear pre-filtering methods for efficient and accurate surface shading[END_REF] (see Figure 2).

We propose a solution to this problem in the important case where the small-scale geometry is generated by displacement mapping. We introduce Linear Efficient Antialiased Displacement and Reflectance (LEADR) mapping, a technique motivated by LEAN mapping but adapted to reflectance of displacement maps. We use the same lightweight memory representation, but with a physically based BRDF that correctly accounts for view-and light-dependent effects under directional, point, and environment lighting. Our model also leverages commodity GPU features, and produces highquality images matching supersampled ground-truth images. When coupled with an adaptive meshing scheme, LEADR mapping provides a seamless, multi-resolution, and hardware-accelerated surface representation. Our contributions can be summarized as follows:

• We model physically based microfacet BRDFs with noncentered Beckmann microfacet distributions.

• We provide anisotropic BRDF equations for specular and diffuse surfaces of arbitrary roughness.

• We derive an analytical solution for the noncentered anisotropic specular microfacet BRDF with a directional/point light.

• We propose a simple, controllable, noise-and artifact-free sampling scheme for our microfacet model under environment lighting.

• We show how surface height scaling, and tangent stretching and shearing can be integrated efficiently in our reflectance model.

Overview

In order to give a general understanding of the workings of our method, and before diving into the mathematical derivations for our solutions, we provide in Figure 3 an overview of its GPU pipeline.

A series of surface statistics are precomputed from a displacement map, and stored in mipmapped textures. In the vertex processing stages, the geometry is adaptively tessellated and displaced for best rasterization performance. The displacements that are too small to be represented by the final geometry are accounted for in our reflectance model during fragment processing, thus guaranteeing consistent and seamless surface shape and appearance at any scale. First, a noncentered Beckmann distribution of the slopes projected onto the footprint of the pixel is computed. This process leverages modern GPU filtering capabilities and only consists of texture fetch operations. The resulting PDF and incident lighting are then combined to compute the diffuse and specular reflectance from the surface using our two new BRDF models (Sections 4 and 5).

Note that a number of additional mathematical derivations, properties, explanations, and pseudocode is provided in the associated supplemental document. The GPU pipeline for our method. The following data is given or precomputed once per scene: LEADR textures and light sources (points and environments). The processing pipeline is divided in two stages. First, the vertex processing computes all the data projected into each pixel. Then, the fragment processing calculates the noncentered Beckmann distribution in the pixel to be combined with the incident lighting in order to compute the diffuse and specular shading contributions to the pixel.

Related Work

Convolution-based normal map filtering methods [START_REF] Fournier | Normal distribution functions and multiple surfaces[END_REF][START_REF] Toksvig | Mipmapping normal maps[END_REF][START_REF] Han | Frequency domain normal map filtering[END_REF]] use the fact that, at any scale, the BRDF is the convolution of a base BRDF and a Normal Distribution Function. Because masking, shadowing, and projection weighting are nonlinear functions of the view, light, and normal directions, incorporating these effects into the convolution is difficult. [START_REF] Tan | Multiresolution reflectance filtering[END_REF]2008] use several Gaussian lobes with a masking-shadowing term, but omit the important view-dependent projection weighting effect and do not normalize their BRDF. [START_REF] Heitz | Filtering color mapped textures and surfaces[END_REF] derive an analytical approximation to evaluate viewdependent normal distributions on the fly, in order to filter albedo textures correlated with visibility on displaced surfaces. Neither the BRDF nor reflectance are accounted for in their model. Other methods precompute the complete BRDF at any location and for any scale, and store it in a Bidirectional Texture Function [START_REF] Cabral | Bidirectional reflection functions from surface bump maps[END_REF][START_REF] Becker | Smooth transitions between bump rendering algorithms[END_REF][START_REF] Ma | Level-of-detail representation of bidirectional texture functions for real-time rendering[END_REF][START_REF] Wu | Characteristic point maps. Computer Graphics Forum[END_REF]]. These methods are capable of capturing view-dependent effects, and could theoretically work with displacement maps; however, the dissuasive memory requirements to store 6D spatially varying BRDFs strongly reduce their practicality.

Our approach concerns three main areas of previous work: LEAN mapping, physically based BRDFs, and filtered importance sampling. We review their respective related work in the following paragraphs.

LEAN Mapping. [START_REF] Olano | LEAN mapping[END_REF] introduce a lightweight representation capable of leveraging the tangent-plane parameterization to express the moments of the normal distribution in the same reference frame. Their formulation allows for linear filtering of the data in a manner that properly captures filtered reflectance. Furthermore, their representation supports anisotropy and multilayer superposition, making it overall an excellent tool for real-time rendering. 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 does not incorporate the important masking, shadowing, and projection weighting effects. Moreover, it is only designed for specular microfacets under point and directional lighting.

Physically based BRDFs. Microfacet theory has been used extensively by the graphics community to derive physically based BRDF models [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF][START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Oren | Generalization of Lambert's reflectance model[END_REF][START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]]. It is only recently though that both production and real-time rendering communities have begun leveraging their use [START_REF] Mcauley | Practical physically-based shading in film and game production[END_REF].

Despite the fact that microfacet BRDF derivations bear many similarities with the problem of filtering reflectance from displaced surfaces, they cannot be used directly. Indeed, BRDF formulations do not account for multiple scales, and more importantly assume centered microfacet distributions. This implies a constant average microsurface normal pointing in the "up" direction of the local frame of parameterization, and prevents their use with LEAN mapping, where the average normal can deviate significantly from this direction. As another illustration, [START_REF] Bruneton | Realtime realistic ocean lighting using seamless transitions from geometry to BRDF[END_REF] use a physically based BRDF model [START_REF] Ross | Detailed analytical approach to the gaussian surface bidirectional reflectance distribution function specular component applied to the sea surface[END_REF] to compute the reflectance of procedural ocean displacements at any scale. Since their formulation does not account for noncentered microfacet normal distributions, they are forced to use an approximation, which arguably works for small angular deviations, but fails when the average microfacet normal deviates significantly from the vertical.

Filtered Importance Sampling. For lighting from an environment map, the BRDF must be integrated over all incident lighting directions. Monte Carlo integration is typically used for this problem but produces noisy results and is not adequate for real-time rendering. Colbert and Křivánek [2007;2008] propose a sampling scheme that leverages filtered environment map representations. Instead of sampling directions from the unfiltered environment map, they use samples from the filtered version of the environment map. For a fixed budget of N samples, the algorithm always produces smooth and noise-free results, while exhibiting convergence as N increases. Such behavior is well suited for real-time rendering purposes, where artifact-free tradeoffs between quality and performance are important. We extend this idea to our more complex BRDF models.

Preliminaries

In this section, we show how displacement map filtering can be directly transposed into the problem of deriving a microfacet reflection model. This key observation allows us to use a derivative of LEAN mapping to solve our new shading equations efficiently in the next sections.

Let P be a flat geometric patch that projects onto a pixel, and whose area is normalized under the measure of a low-pass filter kP satisfying P kP (p) dp = 1. We refer to the normal of P as the macronormal ωg, which yields ωg = (0, 0, 1) t in its tangent space. At a smaller scale, P is perturbed by displacements. This process generates micronormals ωn(p) indexed by locations p ∈ P. These micronormals average to ωn, referred to as the mesonormal of P.

Note that contrary to previous microfacet models, we do not assume that ωg = ωn. This configuration is illustrated in Figure 4. For any

Figure 4: Macro-, meso-, and micronormals. The displacement is applied to points pi on patch P, resulting in the process into different sets of normals.

observation direction ωo, the projected area of the plane defined by the mesonormal is

ωn • ωo ωn • ωg = P ωn(p) • ωo ωn(p) • ωg kP (p) dp. (1)
Here, the dot products ωn • ωo and ωn(p) • ωo are the projection weights due to the observation direction. The two other dot products ωn • ωg and ωn(p) • ωg are due to parameterization Jacobians, and respectively measure the areas of the meso-and microsurfaces.

Note that if ωg = ωn, then ωn • ωg = 1.

In order to shade P, we define the intensity I of the pixel it covers as the integration of a filtering operation over all incident lighting directions Ω:

I = ωn • ωg ωn • ωo Ω L(ωi) P R(p, ωo, ωi) kP (p) dp dωi. (2)
The term P R(p, ωo, ωi) kP (p) dp is the "filtering integral" and involves a local reflectance function R(p, ωo, ωi), assumed constant over P. Additionally, L(ωi) is the incident radiance from direction ωi, and ωn•ωg ωn•ωo is a renormalization factor measuring projected area in direction ωo. Under the geometric optics approximation, we define our local reflectance function as

R(p, ωo, ωi) = V (p, ωi) V (p, ωo)× ρ (ωn(p), ωo, ωi) ωn(p), ωo ωn(p), ωi ωn(p) • ωg (3)
where at location p, V (p, ωo) (resp. V (p, ωi)) denotes the masking (resp. shadowing) Heaviside function, and ρ(ωn(p), ωo, ωi) the surface's BRDF. The operator -,gives the cosine of the angle formed between two vectors, clamped to zero if it is negative.

Similarly to PDFs, (low-pass) energy-conserving filters integrate to one over their domain. Therefore, the filtering integral in Equation (2) can be interpreted as the first moment of a function of random variable p, with PDF kP (p)1 :

P R(p, ωo, ωi) kP (p) dp = E[R(p, ωo, ωi)]. (4)
This observation directly links our filtering problem with microfacet BRDF models, since they also integrate reflectance using the statistics of the surface. This allows us to follow an important assumption: Previous work has shown that treating the visibility functions V (p, -) independently from the other parameters of R is an accurate approximation [START_REF] Ross | Detailed analytical approach to the gaussian surface bidirectional reflectance distribution function specular component applied to the sea surface[END_REF][START_REF] Bruneton | Realtime realistic ocean lighting using seamless transitions from geometry to BRDF[END_REF][START_REF] Heitz | Representing appearance and pre-filtering subpixel data in sparse voxel octrees[END_REF]. Intuitively, this is because occlusions never occur at position p directly, as they are due to distant geometry, contrary the other parameters in R that depend only on the normal ωn(p). Based on this observation, we can reasonably write

E[R(p, ωo, ωi)] = E[V (p, ωo)V (p, ωi)] E[Rn(ωn, ωo, ωi)],
having

E[Rn(ωn, ωo, ωi)]
= Ω ρ(ωn, ωo, ωi) ωn, ωo ωn, ωi D(ωn) dωn.

(5)

Here, D is the exact normal distribution function over P:

D(ω) = P δω(ωn(p)) ωn(p) • ωg kP (p) dp, (6)
and satisfies Ω (ωn • ωg) D(ωn) dωn = 1. Thus, using Equations ( 2) and ( 4), we can derive our "microfacet equation":

I = ωn • ωg ωn • ωo × Ω L(ωi) E[V (p, ωo)V (p, ωi)] E[Rn(ωn, ωo, ωi)] dωi, (7)
expressing the pixel intensity as an integral over all incident lighting directions of the masking-shadowing moment E[V (p, ωo)V (p, ωi)] and the reflectance moment [Rn(ωn, ωo, ωi)]. Compared to Equation ( 2), we only assume that the statistics of V (p, -) and ωn are uncorrelated.

E

Gaussian Surfaces

The distribution of normals as formulated in Equation ( 6) makes Equation ( 7) intractable under real-time constraints. Therefore, we will assume that D takes the form of a Beckmann distribution at all scales. Given this assumption, we discuss the memory requirements of our method, derive an analytical formulation for E[V (p, ωo)V (p, ωi)], and instantiate Equation ( 7) for specular and diffuse microfacets. Our choice of distribution is further motivated in Section 6.

Properties

Elevation and slopes are usually the two parameters considered when studying the statistics of random rough surfaces [START_REF] Bourlier | Effect of correlation between shadowing and shadowed points on the Wagner and Smith monostatic one-dimensional shadowing functions[END_REF]]. Choosing such features to follow Gaussian distributions is especially interesting for our needs because they are flexible and lightweight, but also quite general. For instance, many classes of real-world surfaces, such as stochastic, fractal, or multi-layered surfaces, produce Gaussian statistics that can exhibit behaviors that range from perfectly smooth to very rough. Denoting surface slopes as ñ = (xñ, yñ) t , we write their 2D PDF as

P22(ñ) = exp -1 2 (ñ -E[ñ]) t Σ -1 (ñ -E[ñ]) 2π |Σ| , ( 8
)
where Σ is the covariance matrix

Σ = σ 2 x cxy cxy σ 2 y with σ 2 x = E[x 2 ñ] -E 2 [xñ] σ 2 y = E[y 2 ñ] -E 2 [yñ] cxy = E[xñyñ] -E[xñ] E[yñ]
D(ωn) = P22(ñ) (ωn • ωg) 4 .
(10)

The factor 1 (ωn•ωg ) 4 is due to the Jacobian ∂ ñ ∂ωn = 1 (ωn•ωg ) 3 and the inverse projection 1 (ωn•ωg ) that normalizes the PDF. Since the moments in Equation ( 9) are linear, we store them in mipmapped textures. Linear (mipmapped) filtering can then be exploited to compute the adequate normal distribution at any scale [START_REF] Olano | Normal distribution mapping[END_REF][START_REF] Olano | LEAN mapping[END_REF]. Note also that we can retrieve the mesonormal ωn directly from this linear data:

ωn = (-E[xñ], -E[yñ], 1) t 1 + E 2 [xñ] + E 2 [yñ] . (11)
Thus, similarly to LEAN mapping, we store the terms

E[xñ], E[yñ], E[x 2 ñ], E[y 2 ñ],
and E[xñyñ], i.e., five scalar values, in two mipmapped texture maps.

Since we use the same representation as LEAN mapping, our method inherits other important properties from it: the base surface roughness can be modeled by offsetting the parameters of the covariance matrix, the model is compatible with typical Blinn-Phong shading since the Beckmann distribution can be used to approximate it, and the linearity of the parameters allows for the combination of displacement layers. Another important property is its compatibility with animation, which is addressed next.

Macrosurface Animation

In this section, we study Gaussian microsurface statistics in the case where the macrosurface is subject to deformations.

Surface Stretching and Shearing. In practice, we store the moments of the slopes defined in texture space (u, v), denoted as

E[uñ], E[vñ], E[u 2 ñ], E[v 2 ñ],
and E[uñvñ]. They are linearly filtered at runtime and transformed into texture space slope moments using Equation (9). In order to solve Equation ( 7), this distribution must be expressed in the macrosurface tangent frame (x, y, ωg). Since meshes are rarely perfectly parameterized, the texture parameterization u(x, y), v(x, y) usually produces distortions. This phenomenon is further accentuated in the case of animation. If we denote as h(u, v) the amplitude of the displacement map, the texture space slopes are then defined as the displacement gradient

(uñ, vñ) = ∂h ∂u , ∂h ∂v .
Therefore, the slopes in world space can be defined by Computing the moments of D in this way allows us to support animated geometry without having to update the LEADR mipmap hierarchy during animation. The distribution is adapted on the fly from the distortions, and adequately modulates the resulting appearance, as illustrated in Figure 5.

(

Masking-Shadowing

Another benefit of using Gaussian surfaces is the availability of an experimentally validated solution to the masking-shadowing moment [START_REF] Smith | Geometrical shadowing of a random rough surface[END_REF][START_REF] Bourlier | Effect of correlation between shadowing and shadowed points on the Wagner and Smith monostatic one-dimensional shadowing functions[END_REF]]:

E[V (p, ωo)V (p, ωi)] = 1 1 + Λ(ωo) + Λ(ωi) . (13)
Here, the Λ function from [START_REF] Smith | Geometrical shadowing of a random rough surface[END_REF] is defined as follows for any direction ω = (sin θ cos φ, sin θ sin φ, cos θ):

Λ(ω) = tan θ +∞ cot θ (x -cot θ) P2(x) dx.
The PDF P2(x) is the slope density along a fixed incident direction ω.

It is also Gaussian, hence

P2(x) = 1 σ(φ) √ 2π exp -(x-µ(φ)) 2 2 σ 2 (φ)
, with parameters

µ(φ) = cos φ E[xñ] + sin φ E[yñ], (14)
σ 2 (φ) = cos 2 φ σ 2 x + sin 2 φ σ 2 y + 2 cos φ sin φ cxy. ( 15
)
Analytical solutions as well as approximations have been derived for Λ [START_REF] Bourlier | Effect of correlation between shadowing and shadowed points on the Wagner and Smith monostatic one-dimensional shadowing functions[END_REF]], but only under the assumption that the distribution is centered, i.e., E[xñ] = E[yñ] = 0. Fortunately, the formula can easily be extended to Gaussians with arbitrary mean by shifting its parameter by the mean slope along the incident direction ω. Using Equations ( 14) and ( 15), we generalize

ν = cot θ-µ(φ) σ(φ) √ 2
and combine it with the approximation given by [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]:

Λ(ω) = exp -ν 2 2ν √ π - erfc(ν) 2 ≈ 1.0-1.259ν+0.396ν 2 3.535ν+2.181ν 2 if ν < 1.6 0 otherwise.
Note that ν < 0 may occur, implying that the incident direction ω is below the surface. This should not happen for the computation of Λ(ωo) because the face would then be backface culled. However, this may happen for ωi, and then Λ(ωi) = +∞. In this case, the masking-shadowing term is 0 because the surface is completely in shadow.

Now that we have expressions for D and E[V (p, ωo)V (p, ωi)], we turn to E[Rn(ωn, ωo, ωi)]. This is the final step to getting a complete expression for Equation ( 7), and simply requires that we define the BRDF term ρ of the microgeometry.

Specular Microfacets

A (mirror) specular BRDF is defined as a function of the half-vector

ω h = ωo+ω i ωo+ω i as follows ρ(ωn, ωo, ωi) = ∂ω h ∂ωi F (ω h , ωi) δω h (ωn) ω h , ωi = F (ω h , ωi) δω h (ωn) 4 ω h , ωi 2 ,
where ∂ω h

∂ω i = 1 4 |ω h •ω i |
is the Jacobian of the reflection transformation [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]], and F is the Fresnel term. After substitution in Equation ( 5) we get

E[Rn(ωn, ωo, ωi)] = F (ω h , ωi) D(ω h ) 4 ,
thus leading to the following expression for Equation ( 7):

I = ωn • ωg ωn • ωo Ω i L(ωi) F (ω h , ωi) 4 [1 + Λ(ωo) + Λ(ωi)] D(ω h ) dωi. (16)
Equation ( 16) provides a complete shading formulation for rough specular surfaces. Its physical soundness can easily be validated since it integrates to one without the shadowing Λ(ωi) and Fresnel F terms for any possible configuration ωo,

E[xñ], E[yñ], σ 2 x , σ 2 y , cxy : ωn • ωg ωn • ωo Ω i 1 4 [1 + Λ(ωo)] D(ω h ) dωi = 1.
Point and Directional Lighting Under a point or directional light δω i emitting radiance L, Equation ( 16) takes an analytical form

I = ωn • ωg ωn • ωo L F (ω h , ωi) 4 [1 + Λ(ωo) + Λ(ωi)] D(ω h ). ( 17
)
When the mesonormal is centered (i.e., ωn = ωg), Equation ( 17) is equivalent to previous work [START_REF] Ross | Detailed analytical approach to the gaussian surface bidirectional reflectance distribution function specular component applied to the sea surface[END_REF][START_REF] Bruneton | Realtime realistic ocean lighting using seamless transitions from geometry to BRDF[END_REF]. Thus, our result can be interpreted as a generalization to the case of Gaussian distributions with arbitrary mean.

Environment Lighting For environment lighting, we reformulate the integral in Equation ( 16) in normal space Ωn, rather than Ωi. The Jacobian ∂ω h ∂ω i then cancels out, yielding

I = ωn • ωg ωn • ωo Ωn L(ωi) F (ω h , ωi) ωn, ωo 1 + Λ(ωo) + Λ(ωi) D(ωn) dωn.
Note that because of the change of variable, ωi must be computed using the reflection operator ωi = 2 (ωn • ωo) ωn -ωo.

With a change of variable D(ωn) dωn = P 22 (ñ) ωn•ωg dñ, we obtain

I = ωn • ωg ωn • ωo L(ωi) F (ω h , ωi) ωn,ωo ωn•ωg 1 + Λ(ωo) + Λ(ωi) P22(ñ) dñ, (19)
which expresses the outgoing radiance from the microsurface as an integration over its slope statistics. This formulation is more convenient because efficient numerical solutions can be derived, as we will show in Section 5. Before diving into these details, we first provide our shading formulation for rough diffuse surfaces.

Diffuse Microfacets

A diffuse BRDF is constant ρ(ωn, ωo, ωi) = 1 π . Substituting this formulation in Equation ( 5) yields E[Rn(ωn, ωo, ωi)] = 1 π Ωn ωn, ωo ωn, ωi D(ωn) dωn.

As in the specular case, we first rewrite Equation ( 7) with this specific formulation of ρ, obtaining

I = ωn • ωg ωn • ωo Ω i L(ωi)
1 π Ωn ωn, ωo ωn, ωi D(ωn) dωn 1 + Λ(ωo) + Λ(ωi) dωi which, after the same change of variable D(ωn) dωn = P 22 (ñ) ωn•ωg dñ, leads to

I = ωn • ωg ωn • ωo Ω i L(ωi) 1 π ωn, ωo ωn, ωi ωn • ωg P22(ñ) dñ 1 + Λ(ωo) + Λ(ωi) dωi.
Even under a single directional light, this equation has no analytical solution. It can however be reformulated into another integral involving irradiance mapping [START_REF] Ramamoorthi | An efficient representation for irradiance environment maps[END_REF].

Denoting as E(ωn) the irradiance map satisfying

E(ωn) = Ω i L(ωi) ωn, ωi 1 + Λ(ωo) + Λ(ωi) dωi, (20)
we derive our final shading expression for a rough diffuse surface:

I = ωn • ωg ωn • ωo 1 π E(ωn) ωn, ωo ωn • ωg P22(ñ) dñ. ( 21
)

Sampling Algorithms

The expressions we have derived for diffuse surfaces (Equation ( 21)) as well as specular surfaces in the special case of environment lighting (Equation ( 19)) must be computed numerically.

In order to maintain high performance, we next derive an efficient sampling scheme, specifically adapted to solve these integrals as efficiently as possible.

Sampling Pattern

Our numerical solution is based on a sampling scheme whose samples are regularly and uniformly distributed on the R 2 plane, and centered at the origin. For a 2D set of N 2 samples, we thus only need to store N 2 positions and weights. We initialize each sample by uniformly sampling a standard 1D normal distribution N (0, 1) so that the set covers an arbitrary interval [-λ, λ] × [-λ, λ], chosen for best results (Figure 6(a)). For instance, we set λ = 1.8 for N = 5. The normalized weight Wj associated with point pj is given by

Wj = exp -1 2 p 2 j N k=1 exp(-1 2 p 2 k )
.

At runtime, we warp each sample in order to match P22 in Equations ( 19) and ( 21

xñ = pj σx + E[xñ] yñ = rxy pj + 1 -r 2 xy p k σy + E[yñ] Wñ = Wj W k
where rxy = cxy/(σxσy) is the correlation coefficient of Σ. This whole process is illustrated in Figure 6, and provides a lightweight collection of samples to compute Equations ( 19) and ( 21), further detailed in the next section. Our sample set converges to the exact value as N increases, and we further discuss their properties in Section 6.

Sampling Environment Lighting

We numerically compute Equation (19) using our sample set:

I = ωn • ωg ωn • ωo N 2 L(ωi) F (ω h , ωi) ωn,ωo ωn•ωg 1 + Λ(ωo) + Λ(ωi) Wñ. ( 22
)
For each slope sample ñ, we compute the associated normal ωn and the reflected vector ωi, using respectively Equations ( 11) and ( 18), to sample the environment map. Directly sampling the environment map from ωi could produce noise or banding artifacts, especially with the small values for N that we target. To prevent this, we derive a method inspired by filtered importance sampling [Colbert and Křivánek 2007;Křivánek and Colbert 2008]: Instead of sampling with a set of delta directions, we emulate a smooth reconstruction filter using a prefiltered environment map. For each sample, we compute the level of detail according to the solid angle α covered in the incident direction ωi

α ≈ ∂ωi ∂ ñ A,
where A = (2λ max(σx, σy)/N ) 2 is the surface area of a sample in slope space. Because we use isotropic samples, we further overestimate it by squaring the length of the sample in the largest direction of eccentricity of P22 (hence the term max(σx, σy)). This guarantees that the solid angles overlap across neighboring samples and avoids banding artifacts for any N . For the sake of completeness, we provide the Jacobian of the slope to reflected vector:

∂ωi ∂ ñ = ∂ωi ∂ωn ∂ωn ∂ ñ = 4 |ωn • ωo| |ωn • ωg| 3 .
The whole computation is detailed in Algorithm 1 and illustrated in Figure 7, and provides efficient and artifact-free results to Equation (19).

Sampling Irradiance

We compute Equation ( 21) in three steps. The first two steps consist in building a spherical harmonic (SH) representation of the irradiance map [START_REF] Ramamoorthi | An efficient representation for irradiance environment maps[END_REF] that we defined in Equation ( 20), that we then sample.

Irradiance Map. Since our irradiance map formulation from Equation ( 20) depends on the masking-shadowing term, we compute the map on the fly using the local frame of P. In Algorithm 2, we compute the first nine SH coefficients of the irradiance map due

Wñ = Wj W k ωn = (-xñ, -yñ, 1)/ x 2 ñ + y 2 ñ + 1 ωi = 2 (ωn • ωo) ωn -ωo J = 4 × |ωn • ωo| × |ωn • ωg| 3 α = J × A I+ = Wñ× F (ωn,ω i ) ωn ,ωo ωn •ωg 1+Λ(ωo)+Λ(ω i ) ×textureLod(ωi, LOD(α)) end for return I × ωn•ωg ωn•ωo end function
to point/directional lights, accounting for the masking-shadowing term by weighting their respective energies by

1 1+Λ(ωo)+Λ(ω i ) .
Unfortunately, interactively computing SH coefficients L lm (environment) for an environment map is too costly. Therefore, we precompute them once on the CPU; at runtime, we multiply them by the masking term 1 1+Λ(ωo) alone, and add them to the coefficients computed for point/directional lights. With this approach, we correctly weight point/directional light contributions by the masking-shadowing term but only use the masking term for environment lights.

+ = 0.488603 × (xi, zi, yi) × P (L21, L2-1, L2-2) + = 1.092548×(xizi, yizi, xiyi)×P L20 + = 0.315392 × (3z 2 i -1) × P L22 + = 0.546274 × (x 2 i -y 2 i ) × P end for L lm + = L lm (environment)/(1 + Λ(ωo)) end function 2001] E(ωn) = 0.429043L22(x 2 n -y 2 n ) + 0.743125L20z 2 n + 0.886227L00 -0.247708L20 + 0.858086(L2-2xnyn + L21xnzn + L2-1ynzn) + 1.023328(L11xn + L1-1yn + L10zn).
Irradiance Integration. Once we have computed the irradiance map, we discretize Equation ( 21) as

I = ωn • ωg ωn • ωo 1 π N 2 E(ωn) ωn, ωo ωn • ωg Wñ (23)
and compute it numerically with Algorithm 3. We observed that, contrary to specular sampling, irradiance maps need not be prefiltered since they are already low frequency. This prevents artifacts during sampling.

Algorithm 3 Irradiance map sampling function IRRADIANCEMAPSAMPLING(ωo, E[ñ], σx, σy, cxy)

I = 0 for j, k = 1..N do xñ = pj σx + E[xñ] yñ = rxy pj + 1 -r 2 xy p k σy + E[yñ] Wñ = Wj W k ωn = (-xñ, -yñ, 1)/ x 2 ñ + y 2 ñ + 1 I + = Wñ × ωn,ωo ωn•ωg × E(ωn) end for return I × ωn•ωg ωn•ωo × 1 π end function

Discussion

In this section we review the properties of our models and algorithms, and discuss choices we made to elaborate our representation.

BRDF Models. Our BRDF models are designed to be parameterized with noncentered anisotropic Beckmann distributions, and therefore may be interpreted as generalizations of previous work, thus remaining compatible with them. In the case of centered distributions, our specular model corresponds to the one proposed by [START_REF] Ross | Detailed analytical approach to the gaussian surface bidirectional reflectance distribution function specular component applied to the sea surface[END_REF], which is essentially the same as Walter et al.'s model [2007] used with Beckmann distributions. The only difference lies in the masking-shadowing function: Walter et al. use the separable form 1 1+Λ(ωo) × 1 1+Λ(ω i ) while we use the joint form 1 1+Λ(ωo)+Λ(ω i ) proposed by Ross et al., which better models the correlation between visibility for view and light directions. In the case of centered and isotropic Beckmann distributions, our diffuse BRDF model corresponds to Oren and Nayar's model [1994]. Our diffuse BRDF model is more general because it includes noncentered and anisotropic distributions. Within the scope of this paper we focus on Beckmann distributions because they are convenient for multi-scale representations. However, our general formulation for noncentered distributions presented in Equation ( 7) can be used with other normal distribution functions, such as Phong or GGX [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]].

Mono-vs. Multi-lobe Normal Distributions. Some previous work uses multiple lobes to represent and store normal distributions [START_REF] Fournier | Normal distribution functions and multiple surfaces[END_REF][START_REF] Tan | Multiresolution reflectance filtering[END_REF][START_REF] Han | Frequency domain normal map filtering[END_REF][START_REF] Tan | Filtering and rendering of resolution-dependent reflectance models[END_REF]. While multiple lobes can represent more complex distributions, they are more expensive. Also, their main problem lies in the precomputations. To make interpolation possible, every lobe in a texel must match the same lobe in the neighboring texels [START_REF] Tan | Multiresolution reflectance filtering[END_REF][START_REF] Han | Frequency domain normal map filtering[END_REF][START_REF] Tan | Filtering and rendering of resolution-dependent reflectance models[END_REF]. Solving this problem requires heavy nonlinear optimizations, and matching failures may result in visual artifacts when two nonmatching lobes are interpolated. To allow for lightweight memory storage, simple and fast precomputations, as well as no-failure interpolation, we chose to use a single anisotropic lobe representation.

Shadowing.

Physically based rendering techniques should account for shadowing to ensure energy conservation. Our BRDF model accounts for this effect: shadowing that occurs at the level of the macrosurface (typically using a shadow map or a horizon map) is still represented at the level of the microsurface by the BRDF model; energy is conserved as well. However, in a real-life implementation one may decide to neglect shadowing because of performance or memory constraints. In this case, to ensure consistency across scales, the BRDF masking-shadowing term 1 1+Λ(ωo)+Λ(ω i ) must be replaced by the masking term alone 1 1+Λ(ωo) .

Performance of Our Algorithms. We provide an analytical formulation to evaluate our specular BRDF model with point and directional lights. It requires more computations than LEAN mapping alone because of our masking-shadowing term, but has a similar cost to any other physically based model overall. For environment lights, we propose a sampling scheme. We rely on hardware mipmapping and select the best mipmap level to fetch the environment map by appropriately evaluating the solid angle of each fetch. This ensures noise-and artifact-free evaluations for any number of sample points as opposed to Monte Carlo sampling. The number of samples can then be chosen as a tradeoff between bias and performance. Our diffuse BRDF model has no analytical solution and we resorted to our sampling strategy even for a single light source. We tried to find an analytical fitting with rational polynomials, using Pacanowski et al.'s [2012] fitting library. However, our BRDF is 8D (ωo, ωi, E[ñ], σx, σy) and the rational polynomials required to provide a good fit ended up having more than 70 coefficients, and resulted in lower performance than our integration scheme. To amortize the cost of the integration, we first splat all light directions in an irradiance map, and then make the numerical integration on the resulting irradiance map just once. For L lights and N sample points the complexity of the algorithm is O(L) + O(N ), instead of O(L × N ). The cost of the integration is also amortized as the number of light sources increases. In practice, we use a 5 × 5 grid of points for the diffuse integration, because they capture the main view-dependent effects. Note that the spherical harmonic representation of irradiance maps [START_REF] Ramamoorthi | An efficient representation for irradiance environment maps[END_REF] can result in errors of up to 9%. This is accurate enough for many applications, but one would need to sample the lights independently to produce exact results.

An important property of our algorithm is scalability. The further away the rendered object is, the more details get projected onto the same pixel, and the more samples a typical Monte Carlo integration would require to solve Equation (2). Our algorithm is able to recover complex reflectance with computational complexity independent of the size of the filter extent. This key property comes from the fact that we replaced the nonscalable spatial integration (i.e., Equation ( 2)) by a scalable integration in the space of the surface statistics (i.e., Equation ( 7)).

Uniform Sampling. In Section 5.1 we showed how to sample a Gaussian distribution P22 with uniform 2D grids. We tried other sampling strategies before opting for the grid approach. The importance sampling of P22 with a method such as the one of Colbert and Křivánek [2008] is too sensitive to the random samples configuration for small values of N ; the uniform grid provides better results (see our supplemental document for comparisons). We also tried to use Gauss-Hermite quadratures, which are specifically designed for Gaussian kernels. However, this quadrature scheme only works well with polynomial functions, and Equations ( 22) and ( 23) contain many terms that are poorly approximated by polynomials (the function is not even C 1 ). Again, our simple uniform sampling grid provided better results than Gauss-Hermite quadrature. The sampling pattern may be improved with further analysis [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF][START_REF] Subr | Fourier analysis of stochastic sampling strategies for assessing bias and variance in integration[END_REF], but this is outside of the scope of this paper.

Domain of Validity of Our Method.

Our model and its derivations are based on the fundamental assumption that the displacement map is applied over a locally planar patch. In theory, our 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 1 are small, smooth, curved, and highly specular. As such, they exhibit aliasing that cannot be filtered by LEADR mapping alone. In practice, we deal with this problem by offsetting base roughness heuristically due to curvature, and add it to the slope covariance matrix of the displacement map. This way, we merge the normal distribution of the curved macrosurface with the LEADR map. This "trick" effectively removes the aliasing from small and highly curved macrosurfaces, but is not well defined, as we arbitrarily set the mapping from curvature to roughness. 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].

Results

Renderings. We present a set of very rough spheres rendered under different configurations in Figure 8. The top two rows display ground truth computed with supersampling of geometry synthesized directly from normal and displacement maps, and the bottom row displays the results of our method. The images show that normal mapping fails to capture view-dependent effects due to masking and projection weighting. Since normal map filtering methods target ground truth when only normal maps are used, they cannot reproduce these more realistic effects. On the other hand, LEADR mapping accurately reproduces this more complex appearance. We also applied LEADR mapping on a multiscale ocean displacement map, animated with the method of [START_REF] Bruneton | Realtime realistic ocean lighting using seamless transitions from geometry to BRDF[END_REF] (Figure 10). We use our analytic formula to evaluate shading from the sun (which is approximated as a directional light), and our numerical scheme to compute filtered reflection from the environment map. Figure 5 shows that our technique is robust under extreme deformation and animation.

Finally, we tested our approach on a complex, production-quality geometry asset with several highly detailed displacement maps in Figure 1. This result is part of an animated walking cycle, available in the accompanying video, also computed interactively with our technique. While our model has been designed with real-time rendering in mind, it could prove accurate even for some visual effects productions.

Performance. We report various rendering times in the following tables. These numbers, all expressed in milliseconds, are provided as a mean to quantify the performance impact of our model over more naive, yet commonly used shading methods. All numbers are captured on an Intel Core i7-930 CPU at 2.80 GHz, and an NVidia GTX 480 GPU. To obtain a baseline performance value, we rendered a full-screen quad over a 1024 × 1024 framebuffer with a LEADR map of matching resolution. While no filtering occurs, this configuration guarantees that the fragment shader load is the same across all shading configurations, and gives a lower performance bound since coarser mip levels could benefit from superior cache optimizations.

For one directional light source and specular reflection, the baseline are Constant Blinn-Phong LEAN LEADR 0.055 0.146 0.266 0.267

Here, LEADR mapping is as fast as LEAN mapping despite requiring more arithmetic operations. This is most probably due to the bottleneck texture fetches during rendering from the LEADR map. Compared to a naive Blinn-Phong shader, the performance differences are also quite negligible. The constant color configuration outputs constant fragments, and is provided as a measure of the geometry processing overhead.

In the next table, we use the same configuration, but with environment lighting. We compute our specular (resp. diffuse) shading equations at three grid resolutions. Their timings are compared to a single texture fetch along the mean reflected direction for the specular case, and a single irradiance map evaluation using the mean normal for the diffuse case.

Reflection One sample 3 × 3 5 × 5 7 × 7 specular 0.402 0.433 1.704 4.284 diffuse 0.148 1.067 2.750 5.027

As illustrated in Figure 11, LEADR mapping performance scales linearly according to the number of samples used. We believe that significant gains can be made by optimizing the placement of the samples and/or employing sample pruning, as suggested by Jäckel [2005], in order to obtain similar quality but significantly fewer samples. We leave this to future work. For the complex production-quality T-rex model in Figure 1, our framework achieves performance ranging between 40 and 60 milliseconds at an image resolution of 1024 × 1024. For strong closeup and distant views (≈ 128 × 128 pixel coverage), the rendering times are about 200 and 9 milliseconds, respectively. Note that in the accompanying video, we used an NVidia Quadro 6000 GPU with 6 GB of memory so that we could fit the original, fully detailed 2.8 GB displacement textures. We observed similar performance on this platform.

Overall, our method demonstrates a clear quality improvement over naive shading methods, all while maintaining very high performance. Furthermore, our technique exhibits spatial and temporal coherence and always generates alias-free images closely matching ground truth, across all scales.

Conclusion

We have presented LEADR mapping, a lightweight unified representation for physically based BRDFs and small-scale geometries defined with displacement mapping. LEADR mapping generates smooth and antialiased transitions when zooming in and out of displaced surfaces. The BRDF seamlessly emerges from the microgeometry, and accurately captures complex shading effects. We proposed a simple, lightweight, and efficient sampling scheme to compute shading from complex environment lighting in real time. We validated our model against the ground truth and successfully applied it on static and animating/deforming geometries. Unlike supersampling, our method properly reconstructs shading effects when dense subpixel geometric detail is present, even in the difficult case of highly specular microgeometry.

LEADR mapping is limited to the paradigm of displacement mapping, where the local tangent frame is available. Further filtering would require including the macrosurface itself, with curvature and nonlocality. Along with improving our sampling strategies, as discussed earlier, we believe that our approach will promote new work on this important problem.

Figure 2 :

 2 Figure 2: On real surfaces the distribution of visible normals depends on the view direction because of masking and projection weighting.

  Figure3: The GPU pipeline for our method. The following data is given or precomputed once per scene: LEADR textures and light sources (points and environments). The processing pipeline is divided in two stages. First, the vertex processing computes all the data projected into each pixel. Then, the fragment processing calculates the noncentered Beckmann distribution in the pixel to be combined with the incident lighting in order to compute the diffuse and specular shading contributions to the pixel.

  that slopes are mathematically linked to normals by the relations xñ = -xn/zn and yñ = -yn/zn, our normal distribution function D following an anisotropic Beckmann formulation can be expressed as

Figure 5 :

 5 Figure5: Top row: Our representation allows for efficient rendering under surface stretching. Stretching the surface reduces the amplitude of the slopes. Because our representation is linear, this effect can be included at any scale. Note how the specular reflection becomes anisotropic under significant surface compression. Middle row: A quad with diffuse microfacets is deformed by surface shearing. Shearing the surface changes the slope distribution and thus its appearance. As well as stretching, this effect can be accounted for at any scale. Bottom row: Height downscaling reduces the amplitude of the displacement map. Up close, the bumps disappear, and at a distance, the roughness decreases and the surface becomes specular, leading to a mirror-like reflection. All the insets show closer views of the actual displaced geometry.

  ), which results in a new sample position ñ = (xñ, yñ) and weight Wñ. We accomplish this by multiplying the 2D sample (pj, p k ) by the Cholesky decomposition of the covariance matrix of P22 [Jäckel 2005]:

Figure 6 :

 6 Figure 6: Sampling distribution P22. (a) We precompute uniform 1D samples pj. (b) In the fragment shader we compute the Cartesian product and get 2D samples (pj, p k ). (c) We scale, correlate, and shift the samples to match distribution P22.

Figure 7 :

 7 Figure 7: Filtered environment map sampling. (a) We sample slope distribution P22 and for each micronormal sample ωn we compute reflected vector ωi. (b) The surface covered by the samples in slope space are transformed into solid angles by multiplying by the Jacobian of the transformation. (c) For each reflected direction ωi we fetch the environment map with levels of detail determined by the solid angles. The reconstruction is smooth and noise-free.

  For completeness, we recall the formula to sample from the irradiance map given the SH coefficients[Ramamoorthi and Hanrahan Algorithm 2 Computing the irradiance map function INITIRRADIANCEMAPCOEFS for 0 ≤ l ≤ 2 and -l ≤ m ≤ l do L lm = 0 end for for each directional light of radiance L and direction δω i do P = L/(1 + Λ(ωo) + Λ(ωi)) L00 + = 0.282095 × P (L11, L10, L1-1)

  Figure 9 provides other comparisons with more complex geometry: the pillow (top) is illuminated with an environment light only, and the snake (bottom), with an environment map and four point light sources.

Figure 8 :Figure 9 :

 89 Figure 8: Ground truth generated with procedural normal mapping (top row), and displacement mapping (middle row). 5 × 5 sampled LEADR mapping for diffuse microfacets, and 32 × 32 for specular microfacets (bottom row). More complete tables of rendered spheres are provided in the supplemental document.

Figure 10 :

 10 Figure10: Procedural ocean displacement map rendered with our specular microfacet BRDFs, at multiple scales within an image and between an animated zoom-out sequence. From left to right, the ocean is viewed from progressively increasing distances. The specular highlights are correctly and continuously antialiased while zooming out.

Figure 11 :

 11 Figure 11: LEADR mapping render times under environment light using an increasing number of samples on an NVidia GTX 480.

Some filters (such as kP (p) = sinc(p)) produce negative weights and may not be considered as valid PDFs, but this does not impact the equations.

Acknowledgements

The authors would like to sincerely thank Brent Burley and Disney for giving access to the T-rex model and its detailed textures. Laurent Belcour provided much appreciated feedback on an earlier version of the paper, and Derek Nowrouzezahrai was crucial with his timely suggestions. Dupuy and Heitz acknowledge financial support from Explo'ra Doc from the Rhône-Alpes region. In particular, Poulin and Dupuy acknowledge support from GRAND and NSERC, and Ostromoukhov and Dupuy are supported by the ANR excellence chair AMCQMCSCGA, from which Iehl also benefitted.