**Physically Based Rendering in Filament**  # About This document is part of the [Filament project](https://github.com/google/filament). To report errors in this document please use the [project's issue tracker](https://github.com/google/filament/issues). ## Authors - [Romain Guy](https://github.com/romainguy), [@romainguy](https://twitter.com/romainguy) - [Mathias Agopian](https://github.com/pixelflinger), [@darthmoosious](https://twitter.com/darthmoosious) # Overview Filament is a physically based rendering (PBR) engine for Android. The goal of Filament is to offer a set of tools and APIs for Android developers that will enable them to create high quality 2D and 3D rendering with ease. The goal of this document is to explain the equations and theory behind the material and lighting models used in Filament. This document is intended as a reference for contributors to Filament or developers interested in the inner workings of the engine. We will provide code snippets as needed to make the relationship between theory and practice as clear as possible. This document is not intended as a design document. It focuses solely on algorithms and its content could be used to implement PBR in any engine. However, this document explains why we chose specific algorithms/models over others. Unless noted otherwise, all the 3D renderings present in this document have been generated in-engine (prototype or production). Many of these 3D renderings were captured during the early stages of development of Filament and do not reflect the final quality. ## Principles Real-time rendering is an active area of research and there is a large number of equations, algorithms and implementation to choose from for every single feature that needs to be implemented (the book *Rendering real-time shadows*, for instance, is a 400 pages summary of dozens of shadows rendering techniques). As such, we must first define our goals (or principles, to follow Brent Burley's seminal paper Physically-based shading at Disney [#Burley12]) before we can make informed decisions. Real-time mobile performance : Our primary goal is to design and implement a rendering system able to perform efficiently on mobile platforms. The primary target will be OpenGL ES 3.x class GPUs. Quality : Our rendering system will emphasize overall picture quality. We will however accept quality compromises to support low and medium performance GPUs. Ease of use : Artists need to be able to iterate often and quickly on their assets and our rendering system must allow them to do so intuitively. We must therefore provide parameters that are easy to understand (for instance, no specular power). We also understand that not all developers have the luxury to work with artists. The physically based approach of our system will allow developers to craft visually plausible materials without the need to understand the theory behind our implementation. For both artists and developers, our system will rely on as few parameters as possible to reduce trial and error and allow users to quickly master the material model. In addition, any combination of parameter values should lead to physically plausible results. Physically implausible materials must be hard to create. Familiarity : Our system should use physical units everywhere possible: distances in meters or centimeters, color temperatures in Kelvin, light units in lumens or candelas, etc. Flexibility : A physically based approach must not preclude non-realistic rendering. User interfaces for instance will need unlit materials. Deployment size : While not directly related to the content of this document, it bears emphasizing our desire to keep the rendering library as small as possible so any application can bundle it without increasing the binary to undesirable sizes. ## Physically based rendering We chose to adopt PBR for its benefits from an artistic and production efficient standpoints, and because it is compatible with our goals. Physically based rendering is a rendering method that provides a more accurate representation of materials and how they interact with light when compared to traditional real-time models. The separation of materials and lighting at the core of the PBR method makes it easier to create realistic assets that look accurate in all lighting conditions. # Notation $$
ewcommand{NoL}{n \cdot l}
ewcommand{NoV}{n \cdot v}
ewcommand{NoH}{n \cdot h}
ewcommand{VoH}{v \cdot h}
ewcommand{LoH}{l \cdot h}
ewcommand{fNormal}{f_{0}}
ewcommand{fDiffuse}{f_d}
ewcommand{fSpecular}{f_r}
ewcommand{fX}{f_x}
ewcommand{aa}{\alpha^2}
ewcommand{fGrazing}{f_{90}}
ewcommand{schlick}{F_{Schlick}}
ewcommand{nior}{n_{ior}}
ewcommand{Ed}{E_d}
ewcommand{Lt}{L_{\bot}}
ewcommand{Lout}{L_{out}}
ewcommand{cosTheta}{\left< \cos \theta \right> } $$ The equations found throughout this document use the symbols described in table [symbols]. Symbol | Definition :---------------------------:|:---------------------------| $v$ | View unit vector $l$ | Incident light unit vector $n$ | Surface normal unit vector $h$ | Half unit vector between $l$ and $v$ $f$ | BRDF $\fDiffuse$ | Diffuse component of a BRDF $\fSpecular$ | Specular component of a BRDF $\alpha$ | Roughness, remapped from using input `perceptualRoughness` $\sigma$ | Diffuse reflectance $\Omega$ | Spherical domain $\fNormal$ | Reflectance at normal incidence $\fGrazing$ | Reflectance at grazing angle $\chi^+(a)$ | Heaviside function (1 if $a > 0$ and 0 otherwise) $n_{ior}$ | Index of refraction (IOR) of an interface $\left< \NoL \right>$ | Dot product clamped to [0..1] $\left< a \right>$ | Saturated value (clamped to [0..1]) [Table [symbols]: Symbols definitions] # Material system The sections below describe multiple material models to simplify the description of various surface features such as anisotropy or the clear coat layer. In practice however some of these models are condensed into a single one. For instance, the standard model, the clear coat model and the anisotropic model can be combined to form a single, more flexible and powerful model. Please refer to the [Materials documentation](./Materials.md.html) to get a description of the material models as implemented in Filament. ## Standard model The goal of our model is to represent standard material appearances. A material model is described mathematically by a BSDF (Bidirectional Scattering Distribution Function), which is itself composed of two other functions: the BRDF (Bidirectional Reflectance Distribution Function) and the BTDF (Bidirectional Transmittance Function). Since we aim to model commonly encountered surfaces, our standard material model will focus on the BRDF and ignore the BTDF, or approximate it greatly. Our standard model will therefore only be able to correctly mimic reflective, isotropic, dielectric or conductive surfaces with short mean free paths. The BRDF describes the surface response of a standard material as a function made of two terms: - A diffuse component, or $f_d$ - A specular component, or $f_r$ The relationship between a surface, the surface normal, incident light and these terms is shown in figure [frFd] (we ignore subsurface scattering for now): ![Figure [frFd]: Interaction of the light with a surface using BRDF model with a diffuse term $ f_d $ and a specular term $ f_r $](images/diagram_fr_fd.png) The complete surface response can be expressed as such: $$\begin{equation}\label{brdf} f(v,l)=f_d(v,l)+f_r(v,l) \end{equation}$$ This equation characterizes the surface response for incident light from a single direction. The full rendering equation would require to integrate $l$ over the entire hemisphere. Commonly encountered surfaces are usually not made of a flat interface so we need a model that can characterize the interaction of light with an irregular interface. A microfacet BRDF is a good physically plausible BRDF for that purpose. Such BRDF states that surfaces are not smooth at a micro level, but made of a large number of randomly aligned planar surface fragments, called microfacets. Figure [microfacetVsFlat] shows the difference between a flat interface and an irregular interface at a micro level: ![Figure [microfacetVsFlat]: Irregular interface as modeled by a microfacet model (left) and flat interface (right)](images/diagram_microfacet.png) Only the microfacets whose normal is oriented halfway between the light direction and the view direction will reflect visible light, as shown in figure [microfacets]. ![Figure [microfacets]: Microfacets](images/diagram_macrosurface.png) However, not all microfacets with a properly oriented normal will contribute reflected light as the BRDF takes into account masking and shadowing. This is illustrated in figure [microfacetShadowing]. ![Figure [microfacetShadowing]: Masking and shadowing of microfacets](images/diagram_shadowing_masking.png) A microfacet BRDF is heavily influenced by a _roughness_ parameter which describes how smooth (low roughness) or how rough (high roughness) a surface is at a micro level. The smoother the surface, the more facets are aligned and the more pronounced the reflected light is. The rougher the surface, the fewer facets are oriented towards the camera and incoming light is scattered away from the camera after reflection, giving a blurry aspect to the specular highlights. Figure [roughness] shows surfaces of different roughness and how light interacts with them. ![Figure [roughness]: Varying roughness (from left to right, rough to smooth) and the resulting BRDF specular component lobe](images/diagram_roughness.png) !!! Note: About roughness The roughness parameter as set by the user is called `perceptualRoughness` in the shader snippets throughout this document. The variable called `roughness` is the `perceptualRoughness` with a remapping explained in section [Parameterization]. A microfacet model is described by the following equation (where x stands for the specular or diffuse component): $$\begin{equation} \fX(v,l) = \frac{1}{| \NoV | | \NoL |} \int_\Omega D(m,\alpha) G(v,l,m) f_m(v,l,m) (v \cdot m) (l \cdot m) dm \end{equation}$$ The term $D$ models the distribution of the microfacets (this term is also referred to as the NDF or Normal Distribution Function). This term plays a primordial role in the appearance of surfaces as shown in figure [roughness]. The term $G$ models the visibility (or occlusion or shadow-masking) of the microfacets. Since this equation is valid for both the specular and diffuse components, the difference lies in the microfacet BRDF $f_m$. It is important to note that this equation is used to integrate over the hemisphere at a _micro level_: ![Figure [microLevel]: Modeling the surface response at a single point requires an integration at the micro level](images/diagram_micro_vs_macro.png) The diagram above shows that at a macro level, the surfaces is considered flat. This helps simplify our equations by assuming that a shaded fragment lit from a single direction corresponds to a single point at the surface. At a micro level however, the surface is not flat and we cannot assume a single ray of light anymore (we can however assume that the incident rays are parallel). Since the micro facets will scatter the light in different directions given a bundle of parallel incident rays, we must integrate the surface response over a hemisphere, noted m in the above diagram. It is obviously not practical to compute the full integration over the microfacets hemisphere for each shaded fragment. We will therefore rely on approximations of the integration for both the specular and diffuse components. ## Dielectrics and conductors To better understand some of the equations and behaviors shown below, we must first clearly understand the difference between metallic (conductor) and non-metallic (dielectric) surfaces. We saw earlier that when incident light hits a surface governed by a BRDF, the light is reflected as two separate components: the diffuse reflectance and the specular reflectance. The modelization of this behavior is straightforward as shown in figure [bsdfBrdf]. ![Figure [bsdfBrdf]: Modelization of the BRDF part of a BSDF](images/diagram_fr_fd.png) This modelization is a simplification of how the light actually interacts with the surface. In reality, part of the incident light will penetrate the surface, scatter inside, and exit the surface again as diffuse reflectance. This phenomenon is illustrated in figure [diffuseScattering]. ![Figure [diffuseScattering]: Scattering of diffuse light](images/diagram_scattering.png) Here lies the difference between conductors and dielectrics. There is no subsurface scattering occurring with purely metallic materials, which means there is no diffuse component (and we will see later that this has an influence on the perceived color of the specular component). Scattering happens in dielectrics, which means they have both specular and diffuse components. To properly modelize the BRDF we must therefore distinguish between dielectrics and conductors (scattering not shown for clarity), as shown in figure [dielectricConductor]. ![Figure [dielectricConductor]: BRDF modelization for dielectric and conductor surfaces](images/diagram_brdf_dielectric_conductor.png) ## Energy conservation Energy conservation is one of the key components of a good BRDF for physically based rendering. An energy conservative BRDF states that the total amount of specular and diffuse reflectance energy is less than the total amount of incident energy. Without an energy conservative BRDF, artists must manually ensure that the light reflected off a surface is never more intense than the incident light. ## Specular BRDF For the specular term, $f_r$ is a mirror BRDF that can be modeled with the Fresnel law, noted $F$ in the Cook-Torrance approximation of the microfacet model integration: $$\begin{equation} f_r(v,l) = \frac{D(h, \alpha) G(v, l, \alpha) F(v, h, f0)}{4(\NoV)(\NoL)} \end{equation}$$ Given our real-time constraints, we must use an approximation for the three terms $D$, $G$ and $F$. [#Karis13a] has compiled a great list of formulations for these three terms that can be used with the Cook-Torrance specular BRDF. The sections that follow describe the equations we picked for these terms. ### Normal distribution function (specular D) [#Burley12] observed that long-tailed normal distribution functions (NDF) are a good fit for real-world surfaces. The GGX distribution described in [#Walter07] is a distribution with long-tailed falloff and short peak in the highlights, with a simple formulation suitable for real-time implementations. It is also a popular model, equivalent to the Trowbridge-Reitz distribution, in modern physically based renderers. $$\begin{equation} D_{GGX}(h,\alpha) = \frac{\aa}{\pi ( (\NoH)^2 (\aa - 1) + 1)^2} \end{equation}$$ The GLSL implementation of the NDF, shown in listing [specularD], is simple and efficient. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float D_GGX(float NoH, float roughness) { float a = NoH * roughness; float k = roughness / (1.0 - NoH * NoH + a * a); return k * k * (1.0 / PI); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [specularD]: Implementation of the specular D term in GLSL] We can improve this implementation by using half precision floats. This optimization requires changes to the original equation as there are two problems when computing $1 - (\NoH)^2$ in half-floats. First, this computation suffers from floating point cancellation when $(\NoH)^2$ is close to 1 (highlights). Secondly $\NoH$ does not have enough precision around 1. The solution involves Lagrange's identity: $$\begin{equation} | a \times b |^2 = |a|^2 |b|^2 - (a \cdot b)^2 \end{equation}$$ Since both $n$ and $h$ are unit vectors, $|n \times h|^2 = 1 - (\NoH)^2$. This allows us to compute $1 - (\NoH)^2$ directly with half precision floats by using a simple cross product. Listing [specularDfp16] shows the final optimized implementation. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #define MEDIUMP_FLT_MAX 65504.0 #define saturateMediump(x) min(x, MEDIUMP_FLT_MAX) float D_GGX(float roughness, float NoH, const vec3 n, const vec3 h) { vec3 NxH = cross(n, h); float a = NoH * roughness; float k = roughness / (dot(NxH, NxH) + a * a); float d = k * k * (1.0 / PI); return saturateMediump(d); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [specularDfp16]: Implementation of the specular D term in GLSL optimized for fp16] ### Geometric shadowing (specular G) Eric Heitz showed in [#Heitz14] that the Smith geometric shadowing function is the correct and exact $G$ term to use. The Smith formulation is the following: $$\begin{equation} G(v,l,\alpha) = G_1(l,\alpha) G_1(v,\alpha) \end{equation}$$ $G_1$ can in turn follow several models, and is commonly set to the GGX formulation: $$\begin{equation} G_1(v,\alpha) = G_{GGX}(v,\alpha) = \frac{2 (\NoV)}{\NoV + \sqrt{\aa + (1 - \aa) (\NoV)^2}} \end{equation}$$ The full Smith-GGX formulation thus becomes: $$\begin{equation} G(v,l,\alpha) = \frac{2 (\NoL)}{\NoL + \sqrt{\aa + (1 - \aa) (\NoL)^2}} \frac{2 (\NoV)}{\NoV + \sqrt{\aa + (1 - \aa) (\NoV)^2}} \end{equation}$$ We can observe that the dividends $2 (\NoL)$ and $2 (n \cdot v)$ allow us to simplify the original function $f_r$ by introducing a visibility function $V$: $$\begin{equation} f_r(v,l) = D(h, \alpha) V(v, l, \alpha) F(v, h, f_0) \end{equation}$$ Where: $$\begin{equation} V(v,l,\alpha) = \frac{G(v, l, \alpha)}{4 (\NoV) (\NoL)} = V_1(l,\alpha) V_1(v,\alpha) \end{equation}$$ And: $$\begin{equation} V_1(v,\alpha) = \frac{1}{\NoV + \sqrt{\aa + (1 - \aa) (\NoV)^2}} \end{equation}$$ Heitz notes however that taking the height of the microfacets into account to correlate masking and shadowing leads to more accurate results. He defines the height-correlated Smith function thusly: $$\begin{equation} G(v,l,h,\alpha) = \frac{\chi^+(\VoH) \chi^+(\LoH)}{1 + \Lambda(v) + \Lambda(l)} \end{equation}$$ $$\begin{equation} \Lambda(m) = \frac{-1 + \sqrt{1 + \aa tan^2(\theta_m)}}{2} = \frac{-1 + \sqrt{1 + \aa \frac{(1 - cos^2(\theta_m))}{cos^2(\theta_m)}}}{2} \end{equation}$$ Replacing $cos(\theta_m)$ by $\NoV$, we obtain: $$\begin{equation} \Lambda(v) = \frac{1}{2} \left( \frac{\sqrt{\aa + (1 - \aa)(\NoV)^2}}{\NoV} - 1 \right) \end{equation}$$ From which we can derive the visibility function: $$\begin{equation} V(v,l,\alpha) = \frac{0.5}{\NoL \sqrt{(\NoV)^2 (1 - \aa) + \aa} + \NoV \sqrt{(\NoL)^2 (1 - \aa) + \aa}} \end{equation}$$ The GLSL implementation of the visibility term, shown in listing [specularV], is a bit more expensive than we would like since it requires two `sqrt` operations. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float V_SmithGGXCorrelated(float NoV, float NoL, float roughness) { float a2 = roughness * roughness; float GGXV = NoL * sqrt(NoV * NoV * (1.0 - a2) + a2); float GGXL = NoV * sqrt(NoL * NoL * (1.0 - a2) + a2); return 0.5 / (GGXV + GGXL); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [specularV]: Implementation of the specular V term in GLSL] We can optimize this visibility function by using an approximation after noticing that all the terms under the square roots are squares and that all the terms are in the $[0..1]$ range: $$\begin{equation} V(v,l,\alpha) = \frac{0.5}{\NoL (\NoV (1 - \alpha) + \alpha) + \NoV (\NoL (1 - \alpha) + \alpha)} \end{equation}$$ This approximation is mathematically wrong but saves two square root operations and is good enough for real-time mobile applications, as shown in listing [approximatedSpecularV]. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float V_SmithGGXCorrelatedFast(float NoV, float NoL, float roughness) { float a = roughness; float GGXV = NoL * (NoV * (1.0 - a) + a); float GGXL = NoV * (NoL * (1.0 - a) + a); return 0.5 / (GGXV + GGXL); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [approximatedSpecularV]: Implementation of the approximated specular V term in GLSL] [#Hammon17] proposes the same approximation based on the same observation that the square root can be removed. It does so by rewriting the expressions as _lerps_: $$\begin{equation} V(v,l,\alpha) = \frac{0.5}{lerp(2 (\NoL) (\NoV), \NoL + \NoV, \alpha)} \end{equation}$$ ### Fresnel (specular F) The Fresnel effect plays an important role in the appearance of physically based materials. This effect models the fact that the amount of light the viewer sees reflected from a surface depends on the viewing angle. Large bodies of water are a perfect way to experience this phenomenon, as shown in figure [fresnelLake]. When looking at the water straight down (at normal incidence) you can see through the water. However, when looking further out in the distance (at grazing angle, where perceived light rays are getting parallel to the surface), you will see the specular reflections on the water become more intense. The amount of light reflected depends not only on the viewing angle, but also on the index of refraction (IOR) of the material. At normal incidence (perpendicular to the surface, or 0 degree angle), the amount of light reflected back is noted $\fNormal$ and can be derived from the IOR as we will see in section [Reflectance remapping]. The amount of light reflected back at grazing angle is noted $\fGrazing$ and approaches 100% for smooth materials. ![Figure [fresnelLake]: The Fresnel effect is particularly evident on large bodies of water](images/photo_fresnel_lake.jpg) More formally, the Fresnel term defines how light reflects and refracts at the interface between two different media, or the ratio of reflected and transmitted energy. [#Schlick94] describes an inexpensive approximation of the Fresnel term for the Cook-Torrance specular BRDF: $$\begin{equation} F_{Schlick}(v,h,\fNormal,\fGrazing) = \fNormal + (\fGrazing - \fNormal)(1 - \VoH)^5 \end{equation}$$ The constant $\fNormal$ represents the specular reflectance at normal incidence and is achromatic for dielectrics, and chromatic for metals. The actual value depends on the index of refraction of the interface. The GLSL implementation of this term requires the use of a `pow`, as shown in listing [specularF], which can be replaced by a few multiplications. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vec3 F_Schlick(float u, vec3 f0, float f90) { return f0 + (vec3(f90) - f0) * pow(1.0 - u, 5.0); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [specularF]: Implementation of the specular F term in GLSL] This Fresnel function can be seen as interpolating between the incident specular reflectance and the reflectance at grazing angles, represented here by $\fGrazing$. Observation of real world materials show that both dielectrics and conductors exhibit achromatic specular reflectance at grazing angles and that the Fresnel reflectance is 1.0 at 90 degrees. A more correct $\fGrazing$ is discussed in section [Specular occlusion]. Using $\fGrazing$ set to 1, the Schlick approximation for the Fresnel term can be optimized for scalar operations by refactoring the code slightly. The result is shown in listing [scalarSpecularF]. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vec3 F_Schlick(float u, vec3 f0) { float f = pow(1.0 - u, 5.0); return f + f0 * (1.0 - f); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [scalarSpecularF]: Scalar optimization of the specular F term in GLSL] ## Diffuse BRDF In the diffuse term, $f_m$ is a Lambertian function and the diffuse term of the BRDF becomes: $$\begin{equation} \fDiffuse(v,l) = \frac{\sigma}{\pi} \frac{1}{| \NoV | | \NoL |} \int_\Omega D(m,\alpha) G(v,l,m) (v \cdot m) (l \cdot m) dm \end{equation}$$ Our implementation will instead use a simple Lambertian BRDF that assumes a uniform diffuse response over the microfacets hemisphere: $$\begin{equation} \fDiffuse(v,l) = \frac{\sigma}{\pi} \end{equation}$$ In practice, the diffuse reflectance $\sigma$ is multiplied later, as shown in listing [diffuseBRDF]. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float Fd_Lambert() { return 1.0 / PI; } vec3 Fd = diffuseColor * Fd_Lambert(); ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [diffuseBRDF]: Implementation of the diffuse Lambertian BRDF in GLSL] The Lambertian BRDF is obviously extremely efficient and delivers results close enough to more complex models. However, the diffuse part would ideally be coherent with the specular term and take into account the surface roughness. Both the Disney diffuse BRDF [#Burley12] and Oren-Nayar model [#Oren94] take the roughness into account and create some retro-reflection at grazing angles. Given our constraints we decided that the extra runtime cost does not justify the slight increase in quality. This sophisticated diffuse model also renders image-based and spherical harmonics more difficult to express and implement. For completeness, the Disney diffuse BRDF expressed in [#Burley12] is the following: $$\begin{equation} \fDiffuse(v,l) = \frac{\sigma}{\pi} \schlick(n,l,1,\fGrazing) \schlick(n,v,1,\fGrazing) \end{equation}$$ Where: $$\begin{equation} \fGrazing=0.5 + 2 \cdot \alpha cos^2(\theta_d) \end{equation}$$ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float F_Schlick(float u, float f0, float f90) { return f0 + (f90 - f0) * pow(1.0 - u, 5.0); } float Fd_Burley(float NoV, float NoL, float LoH, float roughness) { float f90 = 0.5 + 2.0 * roughness * LoH * LoH; float lightScatter = F_Schlick(NoL, 1.0, f90); float viewScatter = F_Schlick(NoV, 1.0, f90); return lightScatter * viewScatter * (1.0 / PI); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [diffuseBRDF]: Implementation of the diffuse Disney BRDF in GLSL] Figure [lambert_vs_disney] shows a comparison between a simple Lambertian diffuse BRDF and the higher quality Disney diffuse BRDF, using a fully rough dielectric material. For comparison purposes, the right sphere was mirrored. The surface response is very similar with both BRDFs but the Disney one exhibits some nice retro-reflections at grazing angles (look closely at the left edge of the spheres). ![Figure [lambert_vs_disney]: Comparison between the Lambertian diffuse BRDF (left) and the Disney diffuse BRDF (right)](images/diagram_lambert_vs_disney.png) We could allow artists/developers to choose the Disney diffuse BRDF depending on the quality they desire and the performance of the target device. It is important to note however that the Disney diffuse BRDF is not energy conserving as expressed here. ## Standard model summary **Specular term**: a Cook-Torrance specular microfacet model, with a GGX normal distribution function, a Smith-GGX height-correlated visibility function, and a Schlick Fresnel function. **Diffuse term**: a Lambertian diffuse model. The full GLSL implementation of the standard model is shown in listing [glslBRDF]. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float D_GGX(float NoH, float a) { float a2 = a * a; float f = (NoH * a2 - NoH) * NoH + 1.0; return a2 / (PI * f * f); } vec3 F_Schlick(float u, vec3 f0) { return f0 + (vec3(1.0) - f0) * pow(1.0 - u, 5.0); } float V_SmithGGXCorrelated(float NoV, float NoL, float a) { float a2 = a * a; float GGXL = NoV * sqrt((-NoL * a2 + NoL) * NoL + a2); float GGXV = NoL * sqrt((-NoV * a2 + NoV) * NoV + a2); return 0.5 / (GGXV + GGXL); } float Fd_Lambert() { return 1.0 / PI; } void BRDF(...) { vec3 h = normalize(v + l); float NoV = abs(dot(n, v)) + 1e-5; float NoL = clamp(dot(n, l), 0.0, 1.0); float NoH = clamp(dot(n, h), 0.0, 1.0); float LoH = clamp(dot(l, h), 0.0, 1.0); // perceptually linear roughness to roughness (see parameterization) float roughness = perceptualRoughness * perceptualRoughness; float D = D_GGX(NoH, roughness); vec3 F = F_Schlick(LoH, f0); float V = V_SmithGGXCorrelated(NoV, NoL, roughness); // specular BRDF vec3 Fr = (D * V) * F; // diffuse BRDF vec3 Fd = diffuseColor * Fd_Lambert(); // apply lighting... } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [glslBRDF]: Evaluation of the BRDF in GLSL] ## Improving the BRDFs We mentioned in section [Energy conservation] that energy conservation is one of the key components of a good BRDF. Unfortunately the BRDFs explored previously suffer from two problems that we will examine below. ### Energy gain in diffuse reflectance The Lambert diffuse BRDF does not account for the light that reflects at the surface and that is therefore not able to participate in the diffuse scattering event. [TODO: talk about the issue with fr+fd] ### Energy loss in specular reflectance The Cook-Torrance BRDF we presented earlier attempts to model several events at the microfacet level but does so by accounting for a single bounce of light. This approximation can cause a loss of energy at high roughness, the surface is not energy preserving. Figure [singleVsMultiBounce] shows why this loss of energy occurs. In the single bounce (or single scattering) model, a ray of light hitting the surface can be reflected back onto another microfacet and thus be discarded because of the masking and shadowing term. If we however account for multiple bounces (multiscattering), the same ray of light might escape the microfacet field and be reflected back towards the viewer. ![Figure [singleVsMultiBounce]: Single scattering (left) vs multiscattering](images/diagram_single_vs_multi_scatter.png) Based on this simple explanation, we can intuitively deduce that the rougher a surface is, the higher the chances are that energy gets lost because of the failure to account for multiple scattering events. This loss of energy appears to darken rough materials. Metallic surfaces are particularly affected because all of their reflectance is specular. This darkening effect is illustrated in figure [metallicRoughEnergyLoss]. With multiscattering, energy preservation can be achieved, as shown in figure [metallicRoughEnergyPreservation]. ![Figure [metallicRoughEnergyLoss]: Darkening increases with roughness due to single scattering](images/material_metallic_energy_loss.png) ![Figure [metallicRoughEnergyPreservation]: Energy preservation with multiscattering](images/material_metallic_energy_preservation.png) We can use a white furnace, a uniform lighting environment set to pure white, to validate the energy preservation property of a BRDF. When energy preservation is achieved, a purely reflective metallic surface ($\fNormal = 1$) should be indistinguishable from the background, no matter the roughness of said surface. Figure [whiteFurnaceLoss] shows what such a surface looks like with the specular BRDF presented in the previous sections. The loss of energy as the roughness increases is obvious. In contrast, figure [whiteFurnacePreservation] shows that accounting for multiscattering events addresses the energy loss. ![Figure [whiteFurnaceLoss]: Darkening increases with roughness due to single scattering](images/material_furnace_energy_loss.png) ![Figure [whiteFurnacePreservation]: Energy preservation with multiscattering](images/material_furnace_energy_preservation.png) Multiple-scattering microfacet BRDFs are discussed in depth in [#Heitz16]. Unfortunately this paper only presents a stochastic evaluation of the multiscattering BRDF. This solution is therefore not suitable for real-time rendering. Kulla and Conty present a different approach in [#Kulla17]. Their idea is to add an energy compensation term as an additional BRDF lobe shown in equation $\ref{energyCompensationLobe}$: $$\begin{equation}\label{energyCompensationLobe} f_{ms}(l,v) = \frac{(1 - E(l)) (1 - E(v)) F_{avg}^2 E_{avg}}{\pi (1 - E_{avg}) (1 - F_{avg}(1 - E_{avg}))} \end{equation}$$ Where $E$ is the directional albedo of the specular BRDF $f_r$, with $\fNormal$ set to 1: $$\begin{equation} E(l) = \int_{\Omega} f(l,v) (\NoV) dv \end{equation}$$ The term $E_{avg}$ is the cosine-weighted average of $E$: $$\begin{equation} E_{avg} = 2 \int_0^1 E(\mu) \mu d\mu \end{equation}$$ Similarly, $F_{avg}$ is the cosine-weighted average of the Fresnel term: $$\begin{equation} F_{avg} = 2 \int_0^1 F(\mu) \mu d\mu \end{equation}$$ Both terms $E$ and $E_{avg}$ can be precomputed and stored in lookup tables. while $F_{avg}$ can be greatly simplified when the Schlick approximation is used: $$\begin{equation}\label{averageFresnel} F_{avg} = \frac{1 + 20 \fNormal}{21} \end{equation}$$ This new lobe is combined with the original single scattering lobe, previously noted $f_r$: $$\begin{equation} f_{r}(l,v) = f_{ss}(l,v) + f_{ms}(l,v) \end{equation}$$ In [#Lagarde18], with credit to Emmanuel Turquin, Lagarde and Golubev make the observation that equation $\ref{averageFresnel}$ can be simplified to $\fNormal$. They also propose to apply energy compensation by adding a scaled GGX specular lobe: $$\begin{equation}\label{energyCompensation} f_{ms}(l,v) = \fNormal \frac{1 - E(l)}{E(l)} f_{ss}(l,v) \end{equation}$$ The key insight is that $E(l)$ can not only be precomputed but also shared with image-based lighting pre-integration. The multiscattering energy compensation formula thus becomes: $$\begin{equation}\label{scaledEnergyCompensationLobe} f_r(l,v) = f_{ss}(l,v) + \fNormal \left( \frac{1}{r} - 1 \right) f_{ss}(l,v) \end{equation}$$ Where $r$ is defined as: $$\begin{equation} r = \int_{\Omega} D(l,v) V(l,v) \left< \NoL \right> dl \end{equation}$$ We can implement specular energy compensation at a negligible cost if we store $r$ in the DFG lookup table presented in section [Image based lights]. Listing [energyCompensationImpl] shows that the implementation is a direct conversion of equation $\ref{scaledEnergyCompensationLobe}$. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vec3 energyCompensation = 1.0 + f0 * (1.0 / dfg.y - 1.0); // Scale the specular lobe to account for multiscattering Fr *= pixel.energyCompensation; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [energyCompensationImpl]: Implementation of the energy compensation specular lobe] Please refer to section [Image based lights] and section [Pre-integration for multiscattering] to learn how the DFG lookup table is derived and computed. ## Parameterization Disney's material model described in [#Burley12] is a good starting point but its numerous parameters makes it impractical for real-time implementations. In addition, we would like our standard material model to be easy to understand and easy to use for both artists and developers. ### Standard parameters Table [standardParameters] describes the list of parameters that satisfy our constraints. Parameter | Definition ---------------------:|:--------------------- **BaseColor** | Diffuse albedo for non-metallic surfaces, and specular color for metallic surfaces **Metallic** | Whether a surface appears to be dielectric (0.0) or conductor (1.0). Often used as a binary value (0 or 1) **Roughness** | Perceived smoothness (0.0) or roughness (1.0) of a surface. Smooth surfaces exhibit sharp reflections **Reflectance** | Fresnel reflectance at normal incidence for dielectric surfaces. This replaces an explicit index of refraction **Emissive** | Additional diffuse albedo to simulate emissive surfaces (such as neons, etc.) This parameter is mostly useful in an HDR pipeline with a bloom pass **Ambient occlusion** | Defines how much of the ambient light is accessible to a surface point. It is a per-pixel shadowing factor between 0.0 and 1.0. This parameter will be discussed in more details in the lighting section [Table [standardParameters]: Parameters of the standard model] Figure [material_parameters] shows how the metallic, roughness and reflectance parameters affect the appearance of a surface. ![Figure [material_parameters]: From top to bottom: varying metallic, varying dielectric roughness, varying metallic roughness, varying reflectance](images/material_parameters.png) ### Types and ranges It is important to understand the type and range of the different parameters of our material model, described in table [standardParametersTypes]. Parameter | Type and range ---------------------:|:--------------------- **BaseColor** | Linear RGB [0..1] **Metallic** | Scalar [0..1] **Roughness** | Scalar [0..1] **Reflectance** | Scalar [0..1] **Emissive** | Linear RGB [0..1] + exposure compensation **Ambient occlusion** | Scalar [0..1] [Table [standardParametersTypes]: Range and type of the standard model's parameters] Note that the types and ranges described here are what the shader will expect. The API and/or tools UI could and should allow to specify the parameters using other types and ranges when they are more intuitive for artists. For instance, the base color could be expressed in sRGB space and converted to linear space before being sent off to the shader. It can also be useful for artists to express the metallic, roughness and reflectance parameters as gray values between 0 and 255 (black to white). Another example: the emissive parameter could be expressed as a color temperature and an intensity, to simulate the light emitted by a black body. ### Remapping To make the standard material model easier and more intuitive to use for artists, we must remap the parameters _baseColor_, _roughness_ and _reflectance_. #### Base color remapping The base color of a material is affected by the "metallicness" of said material. Dielectrics have achromatic specular reflectance but retain their base color as the diffuse color. Conductors on the other hand use their base color as the specular color and do not have a diffuse component. The lighting equations must therefore use the diffuse color and $\fNormal$ instead of the base color. The diffuse color can easily be computed from the base color, as show in listing [baseColorToDiffuse]. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vec3 diffuseColor = (1.0 - metallic) * baseColor.rgb; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [baseColorToDiffuse]: Conversion of base color to diffuse in GLSL] #### Reflectance remapping **Dielectrics** The Fresnel term relies on $\fNormal$, the specular reflectance at normal incidence angle, and is achromatic for dielectrics. We will use the remapping for dielectric surfaces described in [#Lagarde14] : $$\begin{equation} \fNormal = 0.16 \cdot reflectance^2 \end{equation}$$ The goal is to map $\fNormal$ onto a range that can represent the Fresnel values of both common dielectric surfaces (4% reflectance) and gemstones (8% to 16%). The mapping function is chosen to yield a 4% Fresnel reflectance value for an input reflectance of 0.5 (or 128 on a linear RGB gray scale). Figure [reflectance] show those common values and how they relate to the mapping function. ![Figure [reflectance]: Common reflectance values](images/diagram_reflectance.png) If the index of refraction is known (for instance, an air-water interface has an IOR of 1.33), the Fresnel reflectance can be calculated as follows: $$\begin{equation}\label{fresnelEquation} \fNormal(n_{ior}) = \frac{(
ior - 1)^2}{(
ior + 1)^2} \end{equation}$$ And if the reflectance value is known, we can compute the corresponding IOR: $$\begin{equation} n_{ior} = \frac{2}{1 - \sqrt{\fNormal}} - 1 \end{equation}$$ Table [commonMatReflectance] describes acceptable Fresnel reflectance values for various types of materials (no real world material has a value under 2%). Material | Reflectance | IOR | Linear value --------------------------:|:-----------------|:-----------------|:---------------- Water | 2% | 1.33 | 0.35 Fabric | 4% to 5.6% | 1.5 to 1.62 | 0.5 to 0.59 Common liquids | 2% to 4% | 1.33 to 1.5 | 0.35 to 0.5 Common gemstones | 5% to 16% | 1.58 to 2.33 | 0.56 to 1.0 Plastics, glass | 4% to 5% | 1.5 to 1.58 | 0.5 to 0.56 Other dielectric materials | 2% to 5% | 1.33 to 1.58 | 0.35 to 0.56 Eyes | 2.5% | 1.38 | 0.39 Skin | 2.8% | 1.4 | 0.42 Hair | 4.6% | 1.55 | 0.54 Teeth | 5.8% | 1.63 | 0.6 Default value | 4% | 1.5 | 0.5 [Table [commonMatReflectance]: Reflectance of common materials (source: Real-Time Rendering 4th Edition)] Table [fNormalMetals] lists the $\fNormal$ values for a few metals. The values are given in sRGB and must be used as the base color in our material model. Please refer to the annex, section [Specular color], for an explanation of how these sRGB colors are computed from measured data. Metal | $\fNormal$ in sRGB | Hexadecimal | Color ----------:|:-------------------:|:------------:|------------------------------------------------------- Silver | 0.97, 0.96, 0.91 | #f7f4e8 |
ior; // complex IOR, n + ik }; static float illuminantD65(float w) { auto i0 = size_t((w - CIE_D65_START) / CIE_D65_INTERVAL); uint2 indexBounds{i0, std::min(i0 + 1, CIE_D65_END)}; float2 wavelengthBounds = CIE_D65_START + float2{indexBounds} * CIE_D65_INTERVAL; float t = (w - wavelengthBounds.x) / (wavelengthBounds.y - wavelengthBounds.x); return lerp(CIE_D65[indexBounds.x], CIE_D65[indexBounds.y], t); } // For std::lower_bound bool operator<(const Sample& lhs, const Sample& rhs) { return lhs.w < rhs.w; } // The wavelength w must be between 360nm and 830nm static std::complex findSample(const std::vector & samples, float w) { auto i1 = std::lower_bound( samples.begin(), samples.end(), Sample{w, 0.0f + 0.0if}); auto i0 = i1 - 1; // Interpolate the complex IORs float t = (w - i0->w) / (i1->w - i0->w); float n = lerp(i0->ior.real(), i1->ior.real(), t); float k = lerp(i0->ior.imag(), i1->ior.imag(), t); return { n, k }; } static float fresnel(const std::complex & sample) { return (((sample - (1.0f + 0if)) * (std::conj(sample) - (1.0f + 0if))) / ((sample + (1.0f + 0if)) * (std::conj(sample) + (1.0f + 0if)))).real(); } static float3 XYZ_to_sRGB(const float3& v) { const mat3f XYZ_sRGB{ 3.2404542f, -0.9692660f, 0.0556434f, -1.5371385f, 1.8760108f, -0.2040259f, -0.4985314f, 0.0415560f, 1.0572252f }; return XYZ_sRGB * v; } // Outputs a linear sRGB color static float3 computeColor(const std::vector & samples) { float3 xyz{0.0f}; float y = 0.0f; for (size_t i = 0; i < CIE_XYZ_COUNT; i++) { // Current wavelength float w = CIE_XYZ_START + i; // Find most appropriate CIE XYZ sample for the wavelength auto sample = findSample(samples, w); // Compute Fresnel reflectance at normal incidence float f0 = fresnel(sample); // We need to multiply by the spectral power distribution of the illuminant float d65 = illuminantD65(w); xyz += f0 * CIE_XYZ[i] * d65; y += CIE_XYZ[i].y * d65; } // Normalize so that 100% reflectance at every wavelength yields Y=1 xyz /= y; float3 linear = XYZ_to_sRGB(xyz); // Normalize out-of-gamut values if (any(greaterThan(linear, float3{1.0f}))) linear *= 1.0f / max(linear); return linear; } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [specularColorImpl]: C++ implementation to compute the base color of a metallic surface from spectral data] Special thanks to Naty Hoffman for his valuable help on this topic. ## Importance sampling for the IBL In the discrete domain, the integral can be approximated with sampling as defined in equation $\ref{iblSampling}$. $$\begin{equation}\label{iblSampling} \Lout(n,v,\Theta) \equiv \frac{1}{N} \sum_{i}^{N} f(l_{i}^{uniform},v,\Theta) L_{\perp}(l_i) \left< n \cdot l_i^{uniform} \right> \end{equation}$$ Unfortunately, we would need too many samples to evaluate this integral. A technique commonly used is to choose samples that are more "important" more often, this is called _importance sampling_. In our case we'll use the distribution of micro-facets normals, $D_{ggx}$, as the distribution of important samples. The evaluation of $ \Lout(n,v,\Theta) $ with importance sampling is presented in equation $\ref{annexIblImportanceSampling}$. $$\begin{equation}\label{annexIblImportanceSampling} \Lout(n,v,\Theta) \equiv \frac{1}{N} \sum_{i}^{N} \frac{f(l_{i},v,\Theta)}{p(l_i,v,\Theta)} L_{\perp}(l_i) \left< n \cdot l_i \right> \end{equation}$$ In equation $\ref{annexIblImportanceSampling}$, $p$ is the probability density function (PDF) of the distribution of _important direction samples_ $l_i$. These samples depend on $h_i$, $v$ and $\alpha$. The definition of the PDF is shown in equation $\ref{iblPDF}$. $h_i$ is given by the distribution we chose, see section [Choosing important directions] for more details. The _important direction samples_ $l_i$ are calculated as the reflection of $v$ around $h_i$, and therefore **do not** have the same PDF as $h_i$. The PDF of a transformed distribution is given by: $$\begin{equation} p(T_r(x)) = p(x) |J(T_r)|^{-1} \end{equation}$$ Where $|J(T_r)|$ is the determinant of the Jacobian of the transform. In our case we're considering the transform from $h_i$ to $l_i$ and the determinant of its Jacobian is given in \ref{iblPDF}. $$\begin{equation}\label{iblPDF} p(l,v,\Theta) = D(h,\alpha) \left< \NoH \right> |J_{h \rightarrow l}|^{-1} \\ |J_{h \rightarrow l}| = 4 \left< \VoH \right> \end{equation}$$ ### Choosing important directions Refer to section [Choosing important directions for sampling the BRDF] for more details. Given a uniform distribution $(\zeta_{\phi},\zeta_{\theta})$ the important direction $l$ is defined by equation $\ref{importantDirection}$. $$\begin{equation}\label{importantDirection} \phi = 2 \pi \zeta_{\phi} \\ \theta = cos^{-1} \sqrt{\frac{1 - \zeta_{\theta}}{(\alpha^2 - 1)\zeta_{\theta}+1}} \\ l = \{ cos \phi sin \theta, sin \phi sin \theta, cos \theta \} \end{equation}$$ Typically, $ (\zeta_{\phi},\zeta_{\theta}) $ are chosen using the Hammersley uniform distribution algorithm described in section [Hammersley sequence]. ### Pre-filtered importance sampling Importance sampling considers only the PDF to generate important directions; in particular, it is oblivious to the actual content of the IBL. If the latter contains high frequencies in areas without a lot of samples, the integration won’t be accurate. This can be somewhat mitigated by using a technique called _pre-filtered importance sampling_, in addition this allows the integral to converge with many fewer samples. Pre-filtered importance sampling uses several images of the environment increasingly low-pass filtered. This is typically implemented very efficiently with mipmaps and a box filter. The LOD is selected based on the sample importance, that is, low probability samples use a higher LOD index (more filtered). This technique is described in details in [#Krivanek08]. The cubemap LOD is determined in the following way: $$\begin{align*} lod &= log_4 \left( K\frac{\Omega_s}{\Omega_p} \right) \\ K &= 4.0 \\ \Omega_s &= \frac{1}{N \cdot p(l_i)} \\ \Omega_p &\approx \frac{4\pi}{6 \cdot width \cdot height} \end{align*}$$ Where $K$ is a constant determined empirically, $p$ the PDF of the BRDF, $ \Omega_{s} $ the solid angle associated to the sample and $\Omega_p$ the solid angle associated with the texel in the cubemap. Cubemap sampling is done using seamless trilinear filtering. It is extremely important to sample the cubemap correctly across faces using OpenGL's seamless sampling feature or any other technique that avoids/reduces seams. Table [importanceSamplingViz] shows a comparison between importance sampling and pre-filtered importance sampling when applied to figure [importanceSamplingRef]. ![Figure [importanceSamplingRef]: Importance sampling image reference](images/image_is_original.png) Samples | Importance sampling | Pre-filtered importance sampling ---------|-------------------------------|--------------------------------------- 4096 |  | 1024 |  |  32 |  |  [Table [importanceSamplingViz]: Importance sampling vs pre-filtered importance sampling with $\alpha = 0.4$] The reference renderer used in the comparison below performs no approximation. In particular, it does not assume $v = n$ and does not perform the split sum approximation. The pre-filtered renderer uses all the techniques discussed in this section: pre-filtered cubemaps, the analytic formulation of the DFG term, and of course the split sum approximation. Left: reference renderer, right: pre-filtered importance sampling.         ## Choosing important directions for sampling the BRDF For simplicity we use the $ D $ term of the BRDF as the PDF, however the PDF must be normalized such that the integral over the hemisphere is 1: $$\begin{equation} \int_{\Omega}p(m)dm = 1 \\ \int_{\Omega}D(m)(n \cdot m)dm = 1 \\ \int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac{\pi}{2}}D(\theta,\phi) cos \theta sin \theta d\theta d\phi = 1 \\ \end{equation}$$ The PDF of the BRDF can therefore be expressed as in equation $\ref{importantPDF}$: $$\begin{equation}\label{importantPDF} p(\theta,\phi) = \frac{\alpha^2}{\pi(cos^2\theta (\alpha^2-1) + 1)^2} cos\theta sin\theta \end{equation}$$ The term $sin\theta$ comes from the differential solid angle $sin\theta d\phi d\theta$ since we integrate over a sphere. We sample $\theta$ and $\phi$ independently: $$\begin{align*} p(\theta) &= \int_0^{2\pi} p(\theta,\phi) d\phi = \frac{2\alpha^2}{(cos^2\theta (\alpha^2-1) + 1)^2} cos\theta sin\theta \\ p(\phi) &= \frac{p(\theta,\phi)}{p(\phi)} = \frac{1}{2\pi} \end{align*}$$ The expression of $ p(\phi) $ is true for an isotropic distribution of normals. We then calculate the cumulative distribution function (CDF) for each variable: $$\begin{align*} P(s_{\phi}) &= \int_{0}^{s_{\phi}} p(\phi) d\phi = \frac{s_{\phi}}{2\pi} \\ P(s_{\theta}) &= \int_{0}^{s_{\theta}} p(\theta) d\theta = 2 \alpha^2 \left( \frac{1}{(2\alpha^4-4\alpha^2+2) cos(s_{\theta})^2 + 2\alpha^2 - 2} - \frac{1}{2\alpha^4-2\alpha^2} \right) \end{align*}$$ We set $ P(s_{\phi}) $ and $ P(s_{\theta}) $ to random variables $ \zeta_{\phi} $ and $ \zeta_{\theta} $ and solve for $ s_{\phi} $ and $ s_{\theta} $ respectively: $$\begin{align*} P(s_{\phi}) &= \zeta_{\phi} \rightarrow s_{\phi} = 2\pi\zeta_{\phi} \\ P(s_{\theta}) &= \zeta_{\theta} \rightarrow s_{\theta} = cos^{-1} \sqrt{\frac{1-\zeta_{\theta}}{(\alpha^2-1)\zeta_{\theta}+1}} \end{align*}$$ So given a uniform distribution $ (\zeta_{\phi},\zeta_{\theta}) $, our important direction $l$ is defined as: $$\begin{align*} \phi &= 2\pi\zeta_{\phi} \\ \theta &= cos^{-1} \sqrt{\frac{1-\zeta_{\theta}}{(\alpha^2-1)\zeta_{\theta}+1}} \\ l &= \{ cos\phi sin\theta,sin\phi sin\theta,cos\theta \} \end{align*}$$ ## Hammersley sequence ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vec2f hammersley(uint i, float numSamples) { uint bits = i; bits = (bits << 16) | (bits >> 16); bits = ((bits & 0x55555555) << 1) | ((bits & 0xAAAAAAAA) >> 1); bits = ((bits & 0x33333333) << 2) | ((bits & 0xCCCCCCCC) >> 2); bits = ((bits & 0x0F0F0F0F) << 4) | ((bits & 0xF0F0F0F0) >> 4); bits = ((bits & 0x00FF00FF) << 8) | ((bits & 0xFF00FF00) >> 8); return vec2f(i / numSamples, bits / exp2(32)); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [C++ implementation of a Hammersley sequence generator] ## Precomputing L for image-based lighting The term $ L_{DFG} $ is only dependent on $ \NoV $. Below, the normal is arbitrarily set to $ n=\left[0, 0, 1\right] $ and $v$ is chosen to satisfy $ \NoV $. The vector $ h_i $ is the $ D_{GGX}(\alpha) $ important direction sample $i$. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ float GDFG(float NoV, float NoL, float a) { float a2 = a * a; float GGXL = NoV * sqrt((-NoL * a2 + NoL) * NoL + a2); float GGXV = NoL * sqrt((-NoV * a2 + NoV) * NoV + a2); return (2 * NoL) / (GGXV + GGXL); } float2 DFG(float NoV, float a) { float3 V; V.x = sqrt(1.0f - NoV*NoV); V.y = 0.0f; V.z = NoV; float2 r = 0.0f; for (uint i = 0; i < sampleCount; i++) { float2 Xi = hammersley(i, sampleCount); float3 H = importanceSampleGGX(Xi, a, N); float3 L = 2.0f * dot(V, H) * H - V; float VoH = saturate(dot(V, H)); float NoL = saturate(L.z); float NoH = saturate(H.z); if (NoL > 0.0f) { float G = GDFG(NoV, NoL, a); float Gv = G * VoH / NoH; float Fc = pow(1 - VoH, 5.0f); r.x += Gv * (1 - Fc); r.y += Gv * Fc; } } return r * (1.0f / sampleCount); } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [C++ implementation of the $ L_{DFG} $ term] ## Spherical Harmonics Symbol | Definition :---------------------------:|:---------------------------| $K^m_l$ | Normalization factors $P^m_l(x)$ | Associated Legendre polynomials $y^m_l$ | Spherical harmonics bases, or SH bases $L^m_l$ | SH coefficients of the $L(s)$ function defined on the unit sphere [Table [shSymbols]: Spherical harmonics symbols definitions] ### Basis functions Spherical parameterization of points on the surface of the unit sphere: $$\begin{equation} \{ x, y, z \} = \{ cos \phi sin \theta, sin \phi sin \theta, cos \theta \} \end{equation}$$ The complex spherical harmonics bases are given by: $$\begin{equation} Y^m_l(\theta, \phi) = K^m_l e^{im\theta} P^{|m|}_l(cos \theta), l \in N, -l <= m <= l \end{equation}$$ However we only need the real bases: $$\begin{align*} y^{m > 0}_l &= \sqrt{2} K^m_l cos(m \phi) P^m_l(cos \theta) \\ y^{m < 0}_l &= \sqrt{2} K^m_l sin(|m| \phi) P^{|m|}_l(cos \theta) \\ y^0_l &= K^0_l P^0_l(cos \theta) \end{align*}$$ The normalization factors are given by: $$\begin{equation} K^m_l = \sqrt{\frac{(2l + 1)(l - |m|)!}{4 \pi (l + |m|)!}} \end{equation}$$ The associated Legendre polynomials $P^{|m|}_l$ can be calculated from the following recursions: $$\begin{equation}\label{shRecursions} P^0_0(x) = 1 \\ P^0_1(x) = x \\ P^l_l(x) = (-1)^l (2l - 1)!! (1 - x^2)^{\frac{l}{2}} \\ P^m_l(x) = \frac{((2l - 1) x P^m_{l - 1} - (l + m - 1) P^m_{l - 2})}{l - m} \\ \end{equation}$$ Computing $y^{|m|}_l$ requires to compute $P^{|m|}_l(z)$ first. This can be accomplished fairly easily using the recursions in equation $\ref{shRecursions}$. The third recursion can be used to "move diagonally" in table [basisFunctions], i.e. calculating $y^0_0$, $y^1_1$, $y^2_2$ etc. Then, the fourth recursion can be used to move vertically. Band index | Basis functions $-l <= m <= l$ :-----------:|:---------------------------------:| $l = 0$ | $y^0_0$ $l = 1$ | $y^{-1}_1$ $y^0_1$ $y^1_1$ $l = 2$ | $y^{-2}_2$ $y^{-1}_2$ $y^0_2$ $y^1_2$ $y^2_2$ [Table [basisFunctions]: Basis functions per band] It’s also fairly easy to compute the trigonometric terms recursively: $$\begin{align*} C_m &\equiv cos(m \phi)sin(\theta)^m \\ S_m &\equiv sin(m \phi)sin(\theta)^m \\ \{ x, y, z \} &= \{ cos \phi sin \theta, sin \phi sin \theta, cos \theta \} \end{align*}$$ Using the angle sum trigonometric identities: $$\begin{align*} cos(m \phi + \phi) &= cos(m \phi) cos(\phi) - sin(m \phi) sin(\phi) \Leftrightarrow C_{m + 1} = x C_m - y S_m \\ sin(m \phi + \phi) &= sin(m \phi) cos(\phi) + cos(m \phi) sin(\phi) \Leftrightarrow S_{m + 1} = x S_m - y C_m \end{align*}$$ Listing [nonNormalizedSHBasis] shows the C++ code to compute the non-normalized SH basis $\frac{y^m_l(s)}{\sqrt{2} K^m_l}$: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ static inline size_t SHindex(ssize_t m, size_t l) { return l * (l + 1) + m; } void computeShBasis( double* const SHb, size_t numBands, const vec3& s) { // handle m=0 separately, since it produces only one coefficient double Pml_2 = 0; double Pml_1 = 1; SHb[0] = Pml_1; for (ssize_t l = 1; l < numBands; l++) { double Pml = ((2 * l - 1) * Pml_1 * s.z - (l - 1) * Pml_2) / l; Pml_2 = Pml_1; Pml_1 = Pml; SHb[SHindex(0, l)] = Pml; } double Pmm = 1; for (ssize_t m = 1; m < numBands ; m++) { Pmm = (1 - 2 * m) * Pmm; double Pml_2 = Pmm; double Pml_1 = (2 * m + 1)*Pmm*s.z; // l == m SHb[SHindex(-m, m)] = Pml_2; SHb[SHindex( m, m)] = Pml_2; if (m + 1 < numBands) { // l == m+1 SHb[SHindex(-m, m + 1)] = Pml_1; SHb[SHindex( m, m + 1)] = Pml_1; for (ssize_t l = m + 2; l < numBands; l++) { double Pml = ((2 * l - 1) * Pml_1 * s.z - (l + m - 1) * Pml_2) / (l - m); Pml_2 = Pml_1; Pml_1 = Pml; SHb[SHindex(-m, l)] = Pml; SHb[SHindex( m, l)] = Pml; } } } double Cm = s.x; double Sm = s.y; for (ssize_t m = 1; m <= numBands ; m++) { for (ssize_t l = m; l < numBands ; l++) { SHb[SHindex(-m, l)] *= Sm; SHb[SHindex( m, l)] *= Cm; } double Cm1 = Cm * s.x - Sm * s.y; double Sm1 = Sm * s.x + Cm * s.y; Cm = Cm1; Sm = Sm1; } } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [nonNormalizedSHBasis]: C++ implementation to compute a non-normalized SH basis] Normalized SH basis functions $y^m_l(s)$ for the first 3 bands: Band | $m = -2$ | $m = -1$ | $m = 0$ | $m = 1$ | $m = 2$ | :-------:|:------------------------------------:|:-------------------------------------:|:---------------------------------------------------:|:-------------------------------------:|:---------------------------------------------:| $l = 0$ | | | $\frac{1}{2}\sqrt{\frac{1}{\pi}}$ | | | $l = 1$ | | $-\frac{1}{2}\sqrt{\frac{3}{\pi}}y$ | $\frac{1}{2}\sqrt{\frac{3}{\pi}}z$ | $-\frac{1}{2}\sqrt{\frac{3}{\pi}}x$ | | $l = 2$ | $\frac{1}{2}\sqrt{\frac{15}{\pi}}xy$ | $-\frac{1}{2}\sqrt{\frac{15}{\pi}}yz$ | $\frac{1}{4}\sqrt{\frac{5}{\pi}}(2z^2 - x^2 - y^2)$ | $-\frac{1}{2}\sqrt{\frac{15}{\pi}}xz$ | $\frac{1}{4}\sqrt{\frac{15}{\pi}}(x^2 - y^2)$ | [Table [basisFunctions]: Normalized basis functions per band] ### Decomposition and reconstruction A function $L(s)$ defined on a sphere is projected to the SH basis as follows: $$\begin{equation} L^m_l = \int_\Omega L(s) y^m_l(s) ds \\ L^m_l = \int_{\theta = 0}^{\pi} \int_{\phi = 0}^{2\pi} L(\theta, \phi) y^m_l(\theta, \phi) sin \theta d\theta d\phi \end{equation}$$ Note that each $L^m_l$ is a vector of 3 values, one for each RGB color channel. The inverse transformation, or reconstruction, or rendering, from the SH coefficients is given by: $$\begin{equation} \hat{L}(s) = \sum_l \sum_{m = -l}^l L^m_l y^m_l(s) \end{equation}$$ ### Decomposition of $\left< cos \theta \right>$ Since $\left< cos \theta \right>$ does not depend on $\phi$ (azimuthal independence), the integral simplifies to: $$\begin{align*} C^0_l &= 2\pi \int_0^{\pi} \left< cos \theta \right> y^0_l(\theta) sin \theta d\theta \\ C^0_l &= 2\pi K^m_l \int_0^{\frac{\pi}{2}} P^0_l(cos \theta) cos \theta sin \theta d\theta \\ C^m_l &= 0, m != 0 \end{align*}$$ In [#Ramamoorthi01] an analytical solution to the integral is described: $$\begin{align*} C_1 &= \sqrt{\frac{\pi}{3}} \\ C_{odd} &= 0 \\ C_{l, even} &= 2\pi \sqrt{\frac{2l + 1}{4\pi}} \frac{(-1)^{\frac{l}{2} - 1}}{(l + 2)(l - 1)} \frac{l!}{2^l (\frac{l!}{2})^2} \end{align*}$$ The first few coefficients are: $$\begin{align*} C_0 &= +0.88623 \\ C_1 &= +1.02333 \\ C_2 &= +0.49542 \\ C_3 &= +0.00000 \\ C_4 &= -0.11078 \end{align*}$$ Very few coefficients are needed to reasonably approximate $\left< cos \theta \right>$, as shown in figure [shCosThetaApprox]. ![Figure [shCosThetaApprox]: Approximation of $cos \theta$ with SH coefficients](images/chart_sh_cos_thera_approx.png) ### Convolution Convolutions by a kernel $h$ that has a circular symmetry can be applied directly and easily in SH space: $$\begin{equation} (h * f)^m_l = \sqrt{\frac{4\pi}{2l + 1}} h^0_l(s) f^m_l(s) \end{equation}$$ Conveniently, $\sqrt{\frac{4\pi}{2l + 1}} = \frac{1}{K^0_l}$, so in practice we pre-multiply $C_l$ by $\frac{1}{K^0_l}$ and we get a simpler expression: $$\begin{equation} \hat{C}_{l, even} = 2\pi \frac{(-1)^{\frac{l}{2} - 1}}{(l + 2)(l - 1)} \frac{l!}{2^l (\frac{l!}{2})^2} \\ \hat{C}_1 = \frac{2\pi}{3} \end{equation}$$ Here is the C++ code to compute $\hat{C}_l$: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ static double factorial(size_t n, size_t d = 1); // < cos(theta) > SH coefficients pre-multiplied by 1 / K(0,l) double computeTruncatedCosSh(size_t l) { if (l == 0) { return M_PI; } else if (l == 1) { return 2 * M_PI / 3; } else if (l & 1) { return 0; } const size_t l_2 = l / 2; double A0 = ((l_2 & 1) ? 1.0 : -1.0) / ((l + 2) * (l - 1)); double A1 = factorial(l, l_2) / (factorial(l_2) * (1 << l)); return 2 * M_PI * A0 * A1; } // returns n! / d! double factorial(size_t n, size_t d ) { d = std::max(size_t(1), d); n = std::max(size_t(1), n); double r = 1.0; if (n == d) { // intentionally left blank } else if (n > d) { for ( ; n>d ; n--) { r *= n; } } else { for ( ; d>n ; d--) { r *= d; } r = 1.0 / r; } return r; } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## Sample validation scene for Mitsuba ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## Light assignment with froxels Assigning lights to froxels can be implemented on the GPU using two compute shaders. The first one, shown in listing [froxelGeneration], creates the froxels data (4 planes + a min Z and max Z per froxel) in an SSBO and needs to be run only once. The shader requires the following uniforms: Projection matrix : The projection matrix used to render the scene (view space to clip space transformation). Inverse projection matrix : The inverse of the projection matrix used to render the scene (clip space to view space transformation). Depth parameters : $-log2(\frac{z_{lighnear}}{z_{far}}) \frac{1}{maxSlices-1}$, maximum number of depth slices, Z near and Z far. Clip space size : $\frac{F_x \times F_r}{w} \times 2$, with $F_x$ the number of tiles on the X axis, $F_r$ the resolution in pixels of a tile and w the width in pixels of the render target. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #version 310 es precision highp float; precision highp int; #define FROXEL_RESOLUTION 80u layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in; layout(location = 0) uniform mat4 projectionMatrix; layout(location = 1) uniform mat4 projectionInverseMatrix; layout(location = 2) uniform vec4 depthParams; // index scale, index bias, near, far layout(location = 3) uniform float clipSpaceSize; struct Froxel { // NOTE: the planes should be stored in vec4[4] but the // Adreno shader compiler has a bug that causes the data // to not be read properly inside the loop vec4 plane0; vec4 plane1; vec4 plane2; vec4 plane3; vec2 minMaxZ; }; layout(binding = 0, std140) writeonly restrict buffer FroxelBuffer { Froxel data[]; } froxels; shared vec4 corners[4]; shared vec2 minMaxZ; vec4 projectionToView(vec4 p) { p = projectionInverseMatrix * p; return p / p.w; } vec4 createPlane(vec4 b, vec4 c) { // standard plane equation, with a at (0, 0, 0) return vec4(normalize(cross(c.xyz, b.xyz)), 1.0); } void main() { uint index = gl_WorkGroupID.x + gl_WorkGroupID.y * gl_NumWorkGroups.x + gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y; if (gl_LocalInvocationIndex == 0u) { // first tile the screen and build the frustum for the current tile vec2 renderTargetSize = vec2(FROXEL_RESOLUTION * gl_NumWorkGroups.xy); vec2 frustumMin = vec2(FROXEL_RESOLUTION * gl_WorkGroupID.xy); vec2 frustumMax = vec2(FROXEL_RESOLUTION * (gl_WorkGroupID.xy + 1u)); corners[0] = vec4( frustumMin.x / renderTargetSize.x * clipSpaceSize - 1.0, (renderTargetSize.y - frustumMin.y) / renderTargetSize.y * clipSpaceSize - 1.0, 1.0, 1.0 ); corners[1] = vec4( frustumMax.x / renderTargetSize.x * clipSpaceSize - 1.0, (renderTargetSize.y - frustumMin.y) / renderTargetSize.y * clipSpaceSize - 1.0, 1.0, 1.0 ); corners[2] = vec4( frustumMax.x / renderTargetSize.x * clipSpaceSize - 1.0, (renderTargetSize.y - frustumMax.y) / renderTargetSize.y * clipSpaceSize - 1.0, 1.0, 1.0 ); corners[3] = vec4( frustumMin.x / renderTargetSize.x * clipSpaceSize - 1.0, (renderTargetSize.y - frustumMax.y) / renderTargetSize.y * clipSpaceSize - 1.0, 1.0, 1.0 ); uint froxelSlice = gl_WorkGroupID.z; minMaxZ = vec2(0.0, 0.0); if (froxelSlice > 0u) { minMaxZ.x = exp2((float(froxelSlice) - depthParams.y) * depthParams.x) * depthParams.w; } minMaxZ.y = exp2((float(froxelSlice + 1u) - depthParams.y) * depthParams.x) * depthParams.w; } if (gl_LocalInvocationIndex == 0u) { vec4 frustum[4]; frustum[0] = projectionToView(corners[0]); frustum[1] = projectionToView(corners[1]); frustum[2] = projectionToView(corners[2]); frustum[3] = projectionToView(corners[3]); froxels.data[index].plane0 = createPlane(frustum[0], frustum[1]); froxels.data[index].plane1 = createPlane(frustum[1], frustum[2]); froxels.data[index].plane2 = createPlane(frustum[2], frustum[3]); froxels.data[index].plane3 = createPlane(frustum[3], frustum[0]); froxels.data[index].minMaxZ = minMaxZ; } } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [froxelGeneration]: GLSL implementation of froxels data generation (compute shader)] The second compute shader, shown in listing [froxelEvaluation], runs every frame (if the camera and/or lights have changed) and assigns all the lights to their respective froxels. This shader relies only on a couple of uniforms (the number of point/spot lights and the view matrix) and four SSBOs: Light index buffer : For each froxel, the index of each light that affects said froxel. The indices for point lights are written first and if there is enough space left, the indices for spot lights are written as well. A sentinel of value 0x7fffffffu separates point and spot lights and/or marks the end of the froxel's list of lights. Each froxel has a maximum number of lights (point + spot). Point lights buffer : Array of structures describing the scene's point lights. Spot lights buffer : Array of structures describing the scene's spot lights. Froxels buffer : The list of froxels represented by planes, created by the previous compute shader. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #version 310 es precision highp float; precision highp int; #define LIGHT_BUFFER_SENTINEL 0x7fffffffu #define MAX_FROXEL_LIGHT_COUNT 32u #define THREADS_PER_FROXEL_X 8u #define THREADS_PER_FROXEL_Y 8u #define THREADS_PER_FROXEL_Z 1u #define THREADS_PER_FROXEL (THREADS_PER_FROXEL_X * \ THREADS_PER_FROXEL_Y * THREADS_PER_FROXEL_Z) layout(local_size_x = THREADS_PER_FROXEL_X, local_size_y = THREADS_PER_FROXEL_Y, local_size_z = THREADS_PER_FROXEL_Z) in; // x = point lights, y = spot lights layout(location = 0) uniform uvec2 totalLightCount; layout(location = 1) uniform mat4 viewMatrix; layout(binding = 0, packed) writeonly restrict buffer LightIndexBuffer { uint index[]; } lightIndexBuffer; struct PointLight { vec4 positionFalloff; // x, y, z, falloff vec4 colorIntensity; // r, g, b, intensity vec4 directionIES; // dir x, dir y, dir z, IES profile index }; layout(binding = 1, std140) readonly restrict buffer PointLightBuffer { PointLight lights[]; } pointLights; struct SpotLight { vec4 positionFalloff; // x, y, z, falloff vec4 colorIntensity; // r, g, b, intensity vec4 directionIES; // dir x, dir y, dir z, IES profile index vec4 angle; // angle scale, angle offset, unused, unused }; layout(binding = 2, std140) readonly restrict buffer SpotLightBuffer { SpotLight lights[]; } spotLights; struct Froxel { // NOTE: the planes should be stored in vec4[4] but the // Adreno shader compiler has a bug that causes the data // to not be read properly inside the loop vec4 plane0; vec4 plane1; vec4 plane2; vec4 plane3; vec2 minMaxZ; }; layout(binding = 3, std140) readonly restrict buffer FroxelBuffer { Froxel data[]; } froxels; shared uint groupLightCounter; shared uint groupLightIndexBuffer[MAX_FROXEL_LIGHT_COUNT]; float signedDistanceFromPlane(vec4 p, vec4 plane) { // plane.w == 0.0, simplify computation return dot(plane.xyz, p.xyz); } void synchronize() { memoryBarrierShared(); barrier(); } void main() { if (gl_LocalInvocationIndex == 0u) { groupLightCounter = 0u; } memoryBarrierShared(); uint froxelIndex = gl_WorkGroupID.x + gl_WorkGroupID.y * gl_NumWorkGroups.x + gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y; Froxel current = froxels.data[froxelIndex]; uint offset = gl_LocalInvocationID.x + gl_LocalInvocationID.y * THREADS_PER_FROXEL_X; for (uint i = 0u; i < totalLightCount.x && groupLightCounter < MAX_FROXEL_LIGHT_COUNT && offset + i < totalLightCount.x; i += THREADS_PER_FROXEL) { uint currentLight = offset + i; vec4 center = pointLights.lights[currentLight].positionFalloff; center.xyz = (viewMatrix * vec4(center.xyz, 1.0)).xyz; float r = inversesqrt(center.w); if (-center.z + r > current.minMaxZ.x && -center.z - r <= current.minMaxZ.y) { if (signedDistanceFromPlane(center, current.plane0) < r && signedDistanceFromPlane(center, current.plane1) < r && signedDistanceFromPlane(center, current.plane2) < r && signedDistanceFromPlane(center, current.plane3) < r) { uint index = atomicAdd(groupLightCounter, 1u); groupLightIndexBuffer[index] = currentLight; } } } synchronize(); uint pointLightCount = groupLightCounter; offset = froxelIndex * MAX_FROXEL_LIGHT_COUNT; for (uint i = gl_LocalInvocationIndex; i < pointLightCount; i += THREADS_PER_FROXEL) { lightIndexBuffer.index[offset + i] = groupLightIndexBuffer[i]; } if (gl_LocalInvocationIndex == 0u) { if (pointLightCount < MAX_FROXEL_LIGHT_COUNT) { lightIndexBuffer.index[offset + pointLightCount] = LIGHT_BUFFER_SENTINEL; } } } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [Listing [froxelEvaluation]: GLSL implementation of assigning lights to froxels (compute shader)] # Revisions February 20, 2019: Cloth shading - Removed Fresnel term from the cloth BRDF - Removed cloth DFG approximations, replaced with a new channel in the DFG LUT August 21, 2018: Multiscattering - Added section [Energy loss in specular reflectance] on how to compensate for energy loss in single scattering BRDFs August 17, 2018: Specular color - Added section [Specular color] to explain how the base color of various metals is computed August 15, 2018: Fresnel - Added a description of the Fresnel effect in section [Fresnel (specular F)] August 9, 2018: Lighting - Added explanation about pre-exposed lights August 7, 2018: Cloth model - Added description of the "Charlie" NDF August 3, 2018: First public version # Bibliography [#Ashdown98]: Ian Ashdown. 1998. Parsing the IESNA LM-63 photometric data file. http://lumen.iee.put.poznan.pl/kw/iesna.txt [#Ashikhmin00]: Michael Ashikhmin, Simon Premoze and Peter Shirley. A Microfacet-based BRDF Generator. *SIGGRAPH '00 Proceedings*, 65-74. [#Ashikhmin07]: Michael Ashikhmin and Simon Premoze. 2007. Distribution-based BRDFs. [#Burley12]: Brent Burley. 2012. Physically Based Shading at Disney. *Physically Based Shading in Film and Game Production, ACM SIGGRAPH 2012 Courses*. [#Estevez17]: Alejandro Conty Estevez and Christopher Kulla. 2017. Production Friendly Microfacet Sheen BRDF. *ACM SIGGRAPH 2017*. [#Hammon17]: Earl Hammon. 217. PBR Diffuse Lighting for GGX+Smith Microsurfaces. *GDC 2017*. [#Heitz14]: Eric Heitz. 2014. Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs. *Journal of Computer Graphics Techniques*, 3 (2). [#Heitz16]: Eric Heitz et al. 2016. Multiple-Scattering Microfacet BSDFs with the Smith Model. *ACM SIGGRAPH 2016*. [#Hill12]: Colin Barré-Brisebois and Stephen Hill. 2012. Blending in Detail. http://blog.selfshadow.com/publications/blending-in-detail/ [#Karis13a]: Brian Karis. 2013. Specular BRDF Reference. http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html [#Karis13b]: Brian Karis, 2013. Real Shading in Unreal Engine 4. https://blog.selfshadow.com/publications/s2013-shading-course/karis/s2013_pbs_epic_notes_v2.pdf [#Karis14]: Brian Karis. 2014. Physically Based Shading on Mobile. https://www.unrealengine.com/blog/physically-based-shading-on-mobile [#Kelemen01]: Csaba Kelemen et al. 2001. A Microfacet Based Coupled Specular-Matte BRDF Model with Importance Sampling. *Eurographics Short Presentations*. [#Krystek85]: M. Krystek. 1985. An algorithm to calculate correlated color temperature. *Color Research & Application*, 10 (1), 38–40. [#Krivanek08]: Jaroslave Krivànek and Mark Colbert. 2008. Real-time Shading with Filtered Importance Sampling. *Eurographics Symposium on Rendering 2008*, Volume 27, Number 4. [#Kulla17]: Christopher Kulla and Alejandro Conty. 2017. Revisiting Physically Based Shading at Imageworks. *ACM SIGGRAPH 2017* [#Lagarde14]: Sébastien Lagarde and Charles de Rousiers. 2014. Moving Frostbite to PBR. *Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2014 Courses*. [#Lagarde18]: Sébastien Lagarde and Evgenii Golubev. 2018. The road toward unified rendering with Unity’s high definition rendering pipeline. *Advances in Real-Time Rendering in Games, ACM SIGGRAPH 2018 Courses*. [#Lazarov13]: Dimitar Lazarov. 2013. Physically-Based Shading in Call of Duty: Black Ops. *Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2013 Courses*. [#McAuley15]: Stephen McAuley. 2015. Rendering the World of Far Cry 4. *GDC 2015*. [#McGuire10]: Morgan McGuire. 2010. Ambient Occlusion Volumes. *High Performance Graphics*. [#Narkowicz14]: Krzysztof Narkowicz. 2014. Analytical DFG Term for IBL. https://knarkowicz.wordpress.com/2014/12/27/analytical-dfg-term-for-ibl [#Neubelt13]: David Neubelt and Matt Pettineo. 2013. Crafting a Next-Gen Material Pipeline for The Order: 1886. *Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2013 Courses*. [#Oren94]: Michael Oren and Shree K. Nayar. 1994. Generalization of lambert's reflectance model. *SIGGRAPH*, 239–246. ACM. [#Pattanaik00]: Sumanta Pattanaik00 et al. 2000. Time-Dependent Visual Adaptation For Fast Realistic Image Display. *SIGGRAPH '00 Proceedings of the 27th annual conference on Computer graphics and interactive techniques*, 47-54. [#Ramamoorthi01]: Ravi Ramamoorthi and Pat Hanrahan. 2001. On the relationship between radiance and irradiance: determining the illumination from images of a convex Lambertian object. *Journal of the Optical Society of America*, Volume 18, Number 10, October 2001. [#Revie12]: Donald Revie. 2012. Implementing Fur in Deferred Shading. *GPU Pro 2*, Chapter 2. [#Russell15]: Jeff Russell. 2015. Horizon Occlusion for Normal Mapped Reflections. http://marmosetco.tumblr.com/post/81245981087 [#Schlick94]: Christophe Schlick. 1994. An Inexpensive BRDF Model for Physically-Based Rendering. *Computer Graphics Forum*, 13 (3), 233–246. [#Walter07]: Bruce Walter et al. 2007. Microfacet Models for Refraction through Rough Surfaces. *Proceedings of the Eurographics Symposium on Rendering*.