Grafikprogrammierung
Bildbasierte Beleuchtung

• Bei der bildbasierte Beleuchtung (engl. "Image-based Lighting") wird die emittierte Strahldichte $L_i$ der Umgebung durch ein sphärische Umgebungsbild (engl. "Environment Map") angegeben
• Selbst ohne Verdeckungen zu berücksichtigen müssten laut Rendering-Gleichung für ein Oberflächenpunkt alle einfallenden Strahldichten $L_i$ mit der BRDF multipliziert und aufsummiert werden. Dies ist wäre aber für eine Echtzeitberechnung zu aufwendig.
• Bei gegegebener parametrisierter BRDF, z.B. Phong BRDF, können die Integrale jedoch vorberechnet werden und parameterisiert abgespeichter werden (Pre-Filtered Environment Map)
• So könnte z.B. bei der Phong-BRDF, der diffuse Anteil durch die Normalenrichtung parametrisiert und der spekulare Anteil durch die Reflektionsrichtung parametrisiert in einem sphärischen Umgebungsbild abgelegt werden
• Häufig werden dabei die Ergebnisse für verschiedene Phong-Glanz-Exponenten $n_s$ in den Mipmap-Leveln einer Textur gespeichert

• Ausgehend von der Rendering-Gleichung wird die modifiziere Phong BRDF eingesetzt (beides Gleichungen bekannt aus Teil 10, Kapitel 1):
  $\begin{align}L_o(\mathbf{v}) &= L_e(\mathbf{v}) + \int\limits_\Omega \mathrm{f}_r(\mathbf{v}, \mathbf{l})\, \, L_i(\mathbf{l}) \cos(\theta) \, d\omega\\ &= L_e(\mathbf{v}) + \int\limits_\Omega \left(\rho_d \frac{1}{\pi} + \rho_s \frac{n_s+2}{2 \pi} \,\cos(\alpha)^{n_s} \right)\, \, L_i(\mathbf{l}) \cos(\theta) \, d\omega \end{align}$
• Betrachten wir zunächst den diffusen Anteil:
  $\begin{align}L_{o,d}(\mathbf{v}) &= \int\limits_\Omega \rho_d \frac{1}{\pi} \, \,L_i(\mathbf{l}) \cos(\theta) \, d\omega = \rho_d \frac{1}{\pi} \int\limits_\Omega L_i(\mathbf{l}) \cos(\theta) \, d\omega\\ &= \rho_d \frac{1}{\pi} \int\limits_{0}^{2\pi}\, \int\limits_{0}^{\pi/2} L_i(\mathbf{l}) \cos(\theta) \sin(\theta) \, d\theta \, d\phi \end{align}$

• Die Approximation des Integrals mit der Riemann-Summe ergibt:
  $\begin{align}L_{o,d}(\mathbf{v}) &= \rho_d \frac{1}{\pi} \int\limits_{0}^{2\pi}\, \int\limits_{0}^{\pi/2} L_i(\mathbf{l}) \cos(\theta) \sin(\theta) \, d\theta \, d\phi\\ &\approx \rho_d \frac{1}{\pi} \sum\limits_{1}^{K}\, \sum\limits_{1}^{J} L_i(\mathbf{l}) \cos(\theta) \sin(\theta) \underbrace{\Delta\theta}_{\frac{\pi/2}{J}} \, \underbrace{\Delta\phi}_{\frac{2\pi}{K}}\\ &= \rho_d \frac{\pi}{K\,J} \sum\limits_{1}^{K}\, \sum\limits_{1}^{J} L_i(\mathbf{l}) \cos(\theta) \sin(\theta) \end{align}$

• Die Riemann-Summe ist in der Regel keine sehr effiziente Lösung, da die Funktion gleichformig abgetastet wird
• Mit der so genannten "Importance Sampling" Theorie kann jedes Integral auch durch folgende Summe approximiert werden:
  $\int\limits_a^b \mathrm{f}(x) \,dx \approx \frac{1}{N} \sum\limits_{n=1}^{N} \frac{\mathrm{f}(x_n)}{\mathrm{p}(x_n)}$
  Dabei ist $\mathrm{p}(x)$ eine beliebige Wahrscheinlichkeitsdichtefunktion (WDF), die die Bedingung erfüllen muss
  $\int\limits_a^b \mathrm{p}(x) \, dx = 1$
• Theoretisch wäre das beste WDF (mit dem geringsten Fehler der Approximation)
  $\mathrm{p}(x) = \frac{\mathrm{f}(x)}{ \int_a^b \mathrm{f}(x)}$
  was bedeutet, dass die WDF der Form der Funktion folgen sollte (d.h. die Dichte sollte höher sein, wenn die Funktionswerte höher sind).

• Die Approximation des Integrals für den diffusen Anteil der Phong BRDF kann somit gelöst werden mit:
  $\begin{align}L_{o,d}(\mathbf{v}) &= \rho_d \frac{1}{\pi} \int\limits_{0}^{2\pi}\, \int\limits_{0}^{\pi/2} \underbrace{L_i(\mathbf{l}) \cos(\theta) \sin(\theta)}_{\mathrm{f}(\theta_n, \phi_n)} \, d\theta \, d\phi\\ &\approx \rho_d \frac{1}{\pi} \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{f}(\theta_n, \phi_n)}{\mathrm{p}(\theta_n, \phi_n)}\\ &= \rho_d \frac{1}{\pi} \frac{1}{N} \sum_{n=1}^{N} \frac{L_i(\mathbf{l}) \cos(\theta_n) \sin(\theta_n)}{\mathrm{p}(\theta_n, \phi_n)} \end{align}$

• Die WDF $ \mathrm{p}(\theta, \phi)$ kann beliebig gewählt werden
• Um ein möglichst gutes Ergebnis nach wenig Abtastungwerten zu bekommen, sollte die WDF der Form der Funktion folgen
• Aber wie können Abtastpositionen mit einer bestimmten WDF generiert werden?
  • In vielen Programmiersprachen ist es leicht gleichverteilte Abtastwerte zu generieren
  • Die Tabellen auf den nächsten Folien zeigen daher die Formeln zur Berechnung von Azimutwinkel $\phi$ und Polarwinkel $\theta$ eines Kugelkoordinatensystems aus zwei gleichverteilten Zufallsvariablen $u$ und $v$ im Bereich [0.0, 1.0] für die Abtastung einer Halbkugel
  • Die entsprechenden Herleitungen können hier gefunden werden:
    Importance Sampling einer Halbkugel
  • Interaktiver Demonstrator:
    Importance Sampling einer Halbkugel

Abtastpositionen der 2D Halten-Sequenz

$y = h_3$

$x = h_2$

• Um wiederholte Abstastung an gleicher Position zu verhinden, werden in der Praxis gerne Pseudozufallssequenzen, wie beispielsweise die Halton-Sequenz verwendet
• Halton-Sequenz (in mehreren Dimensionen)
  $\left(h_2(n), h_3(n), h_5(n), h_7(n), \dots, h_{p_i}(n)\right)^\top$
  Dabei ist $p_i$ die $i$-te Primzahl und $h_r(n)$ berechnet sich, indem der Zahlwert von $n$ (zur Basis $r$) am Dezimalpunkt gespiegelt wird
• Die erzeugten Abtastwerte sind gleichverteilt
• Beliebige Länge der Sequenz möglich

• Beispiele:
  $h_2((26)_{10})= h_2((11010)_2) = (0.01011)_2 = 11/32$
  $h_3((19)_{10})= h_3((201)_3)= (0.102)_3=11/27$

• Wenn die Anzahl $N$ der Abtastwert im Vorfeld feststeht, kann auch die Hammersley-Sequenz verwendet werden, bei der die erste Dimension schneller zu berechnen ist
• Hammersley-Sequenz (in mehreren Dimensionen)
  $\left(\frac{n}{N}, h_2(n), h_3(n), h_5(n), h_7(n), \dots, h_{p_i}(n)\right)^\top$

$y = h_2$

$x = \frac{n}{N}$

• Zurück zum diffusen Anteil der Phong BRDF:
  $\begin{align}L_{o,d}(\mathbf{v}) &= \rho_d \frac{1}{\pi} \int\limits_{0}^{2\pi}\, \int\limits_{0}^{\pi/2} \underbrace{L_i(\mathbf{l}) \cos(\theta) \sin(\theta)}_{\mathrm{f}(\theta, \phi)} \, d\theta \, d\phi \approx \rho_d \frac{1}{\pi} \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{f}(\theta_n, \phi_n)}{\mathrm{p}(\theta_n, \phi_n)}\end{align}$

$\mathbf{v}$

$\mathbf{n}$

• Wenn nun die WDF gewählt wird zu
  $\mathrm{p}(\theta, \phi) = \frac{1}{\pi} \, \cos(\theta) \,\sin(\theta)$
  ergibt sich aufgrund der Unabhänigkeit des diffusen Anteils von der Blickrichtung:
  $\begin{align}L_{o,d}(\mathbf{v}) = L_{o,d}(\mathbf{n}) &\approx \rho_d \underbrace{\frac{1}{N} \sum\limits_{n=1}^{N} L_i(\mathbf{n}, \theta_n, \phi_n)}_{\tiny \mbox{vorberechneter Wert}}\end{align}$
  wobei $\phi_n = 2 \pi \,h_2(n)$ und $\theta_n = \arcsin\left(\sqrt{h_3(n)}\right)$
  und die abgetastete Halbkugel an der Normalen $\mathbf{n}$ ausgerichtet sein muss

• Der spekulare Anteil der Phong BRDF lautet:
  $\begin{align}L_{o,s}(\mathbf{v}) &= \int\limits_\Omega \left(\rho_s \frac{n_s+2}{2 \pi} \,\cos(\alpha)^{n_s} \right)\, \,L_i(\mathbf{l}) \cos(\theta) \, d\omega \end{align}$

$\theta$

$\alpha$

$\mathbf{v}$

$-\mathbf{l}$

$\mathbf{r}$

$\mathbf{n}$

Richtungslichtquelle

$\theta$

$\alpha$

$\mathbf{v}$

$-\mathbf{l}$

$\mathbf{r}$

$\mathbf{n}$

$\alpha$

$\mathbf{r}_v$

Umgebungsbeleuchtung

• Problem: Das Integral hängt sowohl vom Winkel $\alpha$ als auch von $\theta$ ab
• Da letztendlich nur Beiträge um den reflektierten Blickrichtungsvektor $\mathbf{r}_v$ entstehen, wird der Winkel $\theta$ zwischen Lichteinfallsrichtung und Normale als konstant angenommen und mit dem Winkel zwischen der Oberflächennormale $\mathbf{n}$ und $\mathbf{r}_v$ approximiert (entspricht der mittleren Lichteinfallsrichtung):
  $\begin{align}L_{o,s}(\mathbf{v}) &= \int\limits_\Omega \left(\rho_s \frac{n_s+2}{2 \pi} \,\cos(\alpha)^{n_s} \right)\, \,L_i(\mathbf{l}) \cos(\theta) \, d\omega\\ &\approx \cos(\bar{\theta})\,\rho_s \int\limits_\Omega \frac{n_s+2}{2 \pi} \,\cos(\alpha)^{n_s} \, \,L_i(\mathbf{l}) \, d\omega\\ &=\langle \mathbf{n}\cdot \mathbf{r}_v\rangle \,\rho_s \int\limits_\Omega \frac{n_s+2}{2 \pi} \,\cos(\alpha)^{n_s}\, \,L_i(\mathbf{l})\, d\omega \end{align}$

• Wenn die Halbkugel für die Integration an $\mathbf{r}_v$ ausgerichtet wird, gilt:
  $\begin{align}L_{o,s}(\mathbf{v}) &\approx \langle \mathbf{n}\cdot \mathbf{r}_v\rangle \,\rho_s \int\limits_{0}^{2\pi}\, \int\limits_{0}^{\pi/2} \frac{n_s+2}{2 \pi} \,\cos(\theta)^{n_s} \, \,L_i(\mathbf{l})\, \sin(\theta) d\theta \, d\phi \\ &\approx \langle \mathbf{n}\cdot \mathbf{r}_v\rangle \,\rho_s \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{f}(\theta_n, \phi_n)}{\mathrm{p}(\theta_n, \phi_n)}\end{align}$

• Wenn nun die WDF gewählt wird zu
  $\mathrm{p}(\theta, \phi) = \frac{n_s + 1}{2 \pi} \, \cos(\theta)^{n_s} \,\sin(\theta)$
  ergibt sich:
  $\begin{align}L_{o,s}(\mathbf{v}) &\approx \langle \mathbf{n}\cdot \mathbf{r}_v\rangle \rho_s \underbrace{\frac{1}{N} \frac{n_s+2}{n_s+1} \sum\limits_{n=1}^{N} L_i(\mathbf{r}_v, \theta_n, \phi_n)}_{\tiny \mbox{vorberechneter Wert}} \end{align}$
  wobei $\phi_n = 2 \pi \,h_2(n)$ und $\theta_n = \arccos\left((1-h_3(n))^{\frac{1}{n_s+1}}\right)$

• Insgesamt werden damit im Shader nur zwei Zugriffe auf vorberechnete Texturen benötigt:
  $\begin{align}L_{o}(\mathbf{v}) &= L_{o,d}(\mathbf{n}) + L_{o,s}(\mathbf{v})\\ &\approx \rho_d \,T_d(\mathbf{n}) + \langle \mathbf{n}\cdot \mathbf{r}_v\rangle \rho_s \,T_s(\mathbf{r}_v, n_s) \end{align}$
• Vorberechnete Umgebungstextur des diffusen Anteils:
  $T_d(\mathbf{n}) = \frac{1}{N} \sum\limits_{n=1}^{N} L_i(\mathbf{n}, \theta_n, \phi_n)$
  mit $\phi_n = 2 \pi \,h_2(n)$ und $\theta_n = \arcsin\left(\sqrt{h_3(n)}\right)$
• Vorberechnete Umgebungstextur des spekularen Anteils:
  $ T_s(\mathbf{r}_v, n_s) = \frac{1}{N} \frac{n_s+2}{n_s+1} \sum\limits_{n=1}^{N} L_i(\mathbf{r}_v, \theta_n, \phi_n)$
  mit $\phi_n = 2 \pi \,h_2(n)$ und $\theta_n = \arccos\left((1-h_3(n))^{\frac{1}{n_s+1}}\right)$

Environment Map

Pre-filtered diffuser Anteil

Ergebnis

Pre-filtered spekularer Anteil $n_s = 600$

• GSN Composer: EnvmapLighting

Important Sampling des diffusen Anteils:

#version 300 es
precision highp float;
out vec4 outColor;
in vec2 tc; // texture coordinate of the output image in range [0.0, 1.0]
uniform int samples; // number of samples
uniform float envMapLevel; // environment map level
uniform sampler2D envMapImage; // environment image

const float PI = 3.1415926535897932384626433832795;

vec2 directionToSphericalEnvmap(vec3 dir) {
  float s = 1.0 - mod(1.0 / (2.0*PI) * atan(dir.y, dir.x), 1.0);
  float t = 1.0 / (PI) * acos(-dir.z);
  return vec2(s, t);
}

mat3 getNormalSpace(in vec3 normal) {
   vec3 someVec = vec3(1.0, 0.0, 0.0);
   float dd = dot(someVec, normal);
   vec3 tangent = vec3(0.0, 1.0, 0.0);
   if(abs(dd) > 1e-8) {
     tangent = normalize(cross(someVec, normal));
   }
   vec3 bitangent = cross(normal, tangent);
   return mat3(tangent, bitangent, normal);
}

// from http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
// Hacker's Delight, Henry S. Warren, 2001
float radicalInverse(uint bits) {
  bits = (bits << 16u) | (bits >> 16u);
  bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
  bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
  bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
  bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
  return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}

vec2 hammersley(uint n, uint N) {
  return vec2(float(n) / float(N), radicalInverse(n));
}

void main() {
  
  float thetaN = PI * (1.0 - tc.y);
  float phiN = 2.0 * PI * (1.0 - tc.x);
  vec3 normal = vec3(sin(thetaN) * cos(phiN), 
                     sin(thetaN) * sin(phiN), 
                     cos(thetaN));
  mat3 normalSpace = getNormalSpace(normal);

  vec3 result = vec3(0.0);

  uint N = uint(samples);
  
  float r = random2(tc);
  
  for(uint n = 1u; n <= N; n++) {
      vec2 p = hammersley(n, N);
      float theta = asin(sqrt(p.y));
      float phi = 2.0 * PI * p.x;
      vec3 pos = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
      vec3 posGlob = normalSpace * pos;
      vec2 uv = directionToSphericalEnvmap(posGlob);
      vec3 radiance = textureLod(envMapImage, uv, envMapLevel).rgb;
      result +=  radiance;
  }
  result = result / float(samples);
  outColor.rgb = result;
  outColor.a = 1.0;
}

Important Sampling des spekularen Anteils:

#version 300 es
precision highp float;
out vec4 outColor;
in vec2 tc; // texture coordinate of the output image in range [0.0, 1.0]
uniform int samples; // number of samples
uniform float shininess; // specular shininess exponent
uniform float envMapLevel; // environment map level
uniform sampler2D envMapImage; // environment image
const float PI = 3.1415926535897932384626433832795;

vec2 directionToSphericalEnvmap(vec3 dir) {
  float s = 1.0 - mod(1.0 / (2.0*PI) * atan(dir.y, dir.x), 1.0);
  float t = 1.0 / (PI) * acos(-dir.z);
  return vec2(s, t);
}

mat3 getNormalSpace(in vec3 normal) {
   vec3 someVec = vec3(1.0, 0.0, 0.0);
   float dd = dot(someVec, normal);
   vec3 tangent = vec3(0.0, 1.0, 0.0);
   if(abs(dd) > 1e-8) {
     tangent = normalize(cross(someVec, normal));
   }
   vec3 bitangent = cross(normal, tangent);
   return mat3(tangent, bitangent, normal);
}

// from http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
// Hacker's Delight, Henry S. Warren, 2001
float radicalInverse(uint bits) {
  bits = (bits << 16u) | (bits >> 16u);
  bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
  bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
  bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
  bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
  return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}

vec2 hammersley(uint n, uint N) {
  return vec2(float(n) / float(N), radicalInverse(n));
}

void main() {
  
  float thetaN = PI * (1.0 - tc.y);
  float phiN = 2.0 * PI * (1.0 - tc.x);
  vec3 normal = vec3(sin(thetaN) * cos(phiN), 
                     sin(thetaN) * sin(phiN), 
                     cos(thetaN));
  mat3 normalSpace = getNormalSpace(normal);

  vec3 result = vec3(0.0);

  uint N = uint(samples);
  
  float r = random2(tc);
  
  for(uint n = 1u; n <= N; n++) {
      vec2 p = hammersley(n, N);
      float theta = acos(pow(1.0 - p.y, 1.0/(shininess + 1.0)));
      float phi = 2.0 * PI * p.x;
      vec3 pos = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
      vec3 posGlob = normalSpace * pos;
      vec2 uv = directionToSphericalEnvmap(posGlob);
      vec3 radiance = textureLod(envMapImage, uv, envMapLevel).rgb;
      result +=  radiance;
  }
  result = result / float(samples) * (shininess + 2.0) / (shininess + 1.0);
  outColor.rgb = result;
  outColor.a = 1.0;
}

Vertex Shader für das Ergebnisbild:

#version 300 es
precision highp float;
in vec3 position; // input vertex position from mesh
in vec2 texcoord; // input vertex texture coordinate from mesh
in vec3 normal;   // input vertex normal from mesh

uniform mat4 cameraLookAt; //camera look at matrix
uniform mat4 cameraProjection; //camera projection matrix
uniform mat4 meshTransform0; // mesh0 transformation
uniform mat4 meshTransform1; // mesh1 transformation
uniform mat4 meshTransform0TransposedInverse;
uniform mat4 meshTransform1TransposedInverse;

uniform int gsnMeshGroup;

out vec2 tc; // output texture coordinate of vertex
out vec3 wfn; // output fragment normal of vertex in world space
out vec3 vertPos; // output 3D position in world space

void main(){
  mat4 meshTransform;
  mat4 meshTransformTransposedInverse;
  
  if(gsnMeshGroup == 0) {  // transformation of background sphere
   meshTransform = meshTransform0;
   meshTransformTransposedInverse = meshTransform0TransposedInverse;
  } else { // transformation of mesh
   meshTransform = meshTransform1;
   meshTransformTransposedInverse = meshTransform1TransposedInverse;
  }
  tc = texcoord;
  wfn = vec3(meshTransformTransposedInverse * vec4(normal, 0.0));
  vec4 vertPos4 = meshTransform * vec4(position, 1.0);
  vertPos = vec3(vertPos4) / vertPos4.w;
  gl_Position = cameraProjection * cameraLookAt * vertPos4;
}

Fragment Shader für das Ergebnisbild:

#version 300 es
precision highp float; 
precision highp int;
out vec4 outColor;

#define M_PI 3.1415926535897932384626433832795

in vec2 tc; // texture coordinate of pixel (interpolated)
in vec3 wfn; // fragment normal of pixel in world space (interpolated)
in vec3 vertPos; // fragment vertex position in world space (interpolated)

uniform sampler2D envmapBackground; // min_filter="LINEAR" mag_filter="LINEAR"
uniform bool showBackground; // defaultval="true"
uniform sampler2D envmapDiffuse; // min_filter="LINEAR" mag_filter="LINEAR"
uniform sampler2D envmapSpecular; // min_filter="LINEAR" mag_filter="LINEAR"
uniform float diffuseMix; // weighting factor of diffuse color
uniform vec4 diffuseColor; // diffuse color
uniform float specularMix; // weighting factor of specular color
uniform vec4 specularColor; // pecularColor
uniform vec3 cameraPos; // camera position in global coordinate system
uniform int gsnMeshGroup;

vec2 directionToSphericalEnvmap(vec3 dir) {
  float s = 1.0 - mod(1.0 / (2.0*M_PI) * atan(dir.y, dir.x), 1.0);
  float t = 1.0 / (M_PI) * acos(-dir.z);
  return vec2(s, t);
}
  
void main() {
  vec3 normal = normalize(wfn.xyz);
  vec3 viewDir = normalize(cameraPos - vertPos);
  
  vec3 rv = reflect(-viewDir, normal);
  
  if(gsnMeshGroup == 0) {
    if(showBackground) {
      // color of envmap sphere
      outColor.rgb = texture(envmapBackground, vec2(1.0-tc.x, tc.y)).rgb;
      outColor.a = 1.0;
    } else {
      discard;
    }
  } else {

    vec3 diff = texture(envmapDiffuse, directionToSphericalEnvmap(normal)).rgb;
    vec3 rd = diffuseMix * pow(diffuseColor.rgb, vec3(2.2));
    vec3 rs = specularMix * pow(specularColor.rgb, vec3(2.2));
    
    // shading front-facing
    vec3 color = rd * diff;
    float rn = dot(rv, normal);
    if(rn > 0.0) {
      vec3 spec = texture(envmapSpecular, directionToSphericalEnvmap(rv)).rgb;
      color += rs * rn * spec;
    }
    
    // shading back-facing
    if(dot(viewDir, normal) < -0.1) {
      color = 0.1 * rs;
    }

    outColor.rgb = pow(color, vec3(1.0/2.2));
    outColor.a = 1.0;
  }
}

• Die GGX Mikrofacetten BRDF (eingeführt in Teil 10, Kapitel 1) wird von aktuellen physik-basierte Renderverfahren ("Physically-Based Rendering", PBR) eingesetzt [Disney 2012] [Unreal Engine 2013] [Frostbite 2014]
• Der spekulare Anteil der GGX Mikrofacetten BRDF lautet:
  $\begin{align}L_{o,s}(\mathbf{v}) &= \int\limits_\Omega \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{4\,\,\langle\mathbf{n} \cdot \mathbf{l}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \cos(\theta) \, d\omega \end{align}$
• Dies ist ein Integral über alle einfallenden Lichtrichtungen $\mathbf{l}$
• Der Verteilung der Mikrofacetten $\mathrm{D}(\mathbf{h})$ ist eine Funktion des Halbvektors $\mathbf{h}$
  $\mathrm{D}_{\tiny \mbox{GGX}}(\mathbf{h}) = \frac{\alpha^2}{\pi \left(\langle\mathbf{n} \cdot \mathbf{h}\rangle^2 (\alpha^2-1)+1\right)^2}$
• Da diese Verteilung als Funktion bekannt ist, soll das Integral durch Importance Sampling über den Halbvektor $\mathbf{h}$ gelöst werden

$\mathbf{v}$

$-\mathbf{l}$

$\mathbf{h}$

$\mathbf{n}$

• Zusammenhang zwischen $\mathbf{l}$ und $\mathbf{h}$
• Wie aus der Zeichnung ersichtlich, kann $\mathbf{l}$ durch Reflektion von $\mathbf{v}$ an $\mathbf{h}$ berechnet werden
• Damit gilt:
  $\mathbf{l} = 2 \,(\mathbf{v}^\top \mathbf{h}) \, \mathbf{h} - \mathbf{v}$

• Integration durch Substitution:
  $\begin{align}L_{o,s}(\mathbf{v}) &= \int\limits_{\Omega_h} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{4\,\,\langle\mathbf{n} \cdot \mathbf{l}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \cos(\theta) \, d\omega\\ &= \int\limits_{\Omega_h} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{4\,\,\langle\mathbf{n} \cdot \mathbf{l}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \cos(\theta) \, \frac{d\omega}{d\omega_h} \,d\omega_h\\ &= \int\limits_{\Omega_h} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{4\,\,\langle\mathbf{n} \cdot \mathbf{l}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \cos(\theta) \, \frac{\sin(\theta)\, d\theta \,d\phi}{\sin(\theta_h) \,d\theta_h d\phi_h} \,d\omega_h \end{align}$
  

  $d\omega = \sin(\theta)\, d\phi\, \cdot d\theta $

with $\theta^\ast = 2 \,\theta_h^\ast$ and $\phi^\ast = \phi_h^\ast$

• Importance sampling:
  $\begin{align}L_{o,s}(\mathbf{v}) &= \int\limits_{0}^{2\pi}\,\int\limits_{0}^{\frac{\pi}{2}}\underbrace{ \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle \,\sin(\theta_h)}_{\mathrm{f}(\theta_h, \phi_h)} \, d\theta_h \, d\phi_h\\ &\approx \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{f}(\theta_{h}, \phi_{h})}{\mathrm{p}(\theta_{h}, \phi_{h})}\\[1.5ex] &=\frac{1}{N} \sum_{n=1}^{N} \frac{\frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{D}(\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v})}{\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle \,\sin(\theta_h)}{\mathrm{D}(\mathbf{h}) \cos(\theta_h) \sin(\theta_h)}\\[1.5ex] &= \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l})\\ \end{align}$

• Insgesamt ergibt sich somit:
  $\begin{align}L_{o,s}(\mathbf{v}) &\approx \frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l})\\ \end{align}$
  mit $\mathbf{h} = \begin{pmatrix}\sin(\theta_h) cos(\phi_h) \\ \sin(\theta_h) \sin(\phi_h)\\ \cos(\theta_h)\end{pmatrix}$  und   $\mathbf{l} = \operatorname{reflect}(-\mathbf{v}, \mathbf{h})$
  wobei $\phi_h = 2 \pi \,h_2(n)$ und $\theta_h = \arccos\left(\sqrt{\frac{1 - h_3(n)}{h_3(n) (\alpha^2-1) + 1} }\right)$
  und $\operatorname{reflect}$ der GLSL Funktion reflect entspricht:
  $\operatorname{reflect}(\mathbf{i}, \mathbf{n}) = \mathbf{i} - 2 \,(\mathbf{i}^\top \mathbf{n}) \, \mathbf{n}$

• GSN Composer: SpheresIBLRef

Fragment Shader (Vertex Shader, wie bisher):

#version 300 es
precision highp float; 
precision highp int;
out vec4 outColor;

#define PI 3.1415926535897932384626433832795

in vec2 tc; // texture coordinate of pixel (interpolated)
in vec3 wfn; // fragment normal of pixel (interpolated)
in vec3 vertPos; // fragment vertex position (interpolated)

uniform sampler2D envMapImage; // environment images
uniform sampler2D baseColorTexture; // base color 
uniform sampler2D roughnessTexture; // roughness texture
uniform sampler2D metallicTexture; // metallic parameter
uniform sampler2D emissionTexture; // emission texture
uniform float reflectance; // Fresnel reflectance
uniform bool showBackground; 
uniform int samplesSpec; //number of samples
uniform float envLevelSpec; // level for envmap lookup
uniform int samplesDiff; // number of samples
uniform float envLevelDiff; // level for envmap lookup
uniform vec3 cameraPos; // camera position in global coordinate system
uniform int gsnMeshGroup;

vec2 directionToSphericalEnvmap(vec3 dir) {
  float s = 1.0 - mod(1.0 / (2.0*PI) * atan(dir.y, dir.x), 1.0);
  float t = 1.0 / (PI) * acos(-dir.z);
  return vec2(s, t);
}

mat3 getNormalSpace(in vec3 normal) {
   vec3 someVec = vec3(1.0, 0.0, 0.0);
   float dd = dot(someVec, normal);
   vec3 tangent = vec3(0.0, 1.0, 0.0);
   if(abs(dd) > 1e-8) {
     tangent = normalize(cross(someVec, normal));
   }
   vec3 bitangent = cross(normal, tangent);
   return mat3(tangent, bitangent, normal);
}

// from http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
// Hacker's Delight, Henry S. Warren, 2001
float radicalInverse(uint bits) {
  bits = (bits << 16u) | (bits >> 16u);
  bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
  bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
  bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
  bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
  return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}

vec2 hammersley(uint n, uint N) {
  return vec2(float(n) / float(N), radicalInverse(n));
}

float random2(vec2 n) { 
	return fract(sin(dot(n, vec2(12.9898, 4.1414))) * 43758.5453);
}
  
float G1_GGX_Schlick(float NdotV, float roughness) {
  float r = roughness; // original
  //float r = 0.5 + 0.5 * roughness; // Disney remapping
  float k = (r * r) / 2.0;
  float denom = NdotV * (1.0 - k) + k;
  return NdotV / denom;
}

float G_Smith(float NoV, float NoL, float roughness) {
  float g1_l = G1_GGX_Schlick(NoL, roughness);
  float g1_v = G1_GGX_Schlick(NoV, roughness);
  return g1_l * g1_v;
}

vec3 fresnelSchlick(float cosTheta, vec3 F0) {
  return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
} 

// adapted from "Real Shading in Unreal Engine 4", Brian Karis, Epic Games
vec3 specularIBLReference(vec3 F0 , float roughness, vec3 N, vec3 V) {
  mat3 normalSpace = getNormalSpace(N);
  vec3 result = vec3(0.0);
  uint sampleCount = uint(samplesSpec);
  float r = random2(tc);
  for(uint n = 1u; n <= sampleCount; n++) {
    //vec2 p = hammersley(n, N);
    vec2 p = mod(hammersley(n, sampleCount) + r, 1.0);
    float a = roughness * roughness;
    float theta = acos(sqrt((1.0 - p.y) / (1.0 + (a * a - 1.0) * p.y)));
    float phi = 2.0 * PI * p.x;
    // sampled h direction in normal space
    vec3 Hn = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
    // sampled h direction in world space
    vec3 H = normalSpace * Hn;
    vec3 L = 2.0 * dot(V, H) * H - V;

    // all required dot products
    float NoV = clamp(dot(N, V), 0.0, 1.0);
    float NoL = clamp(dot(N, L), 0.0, 1.0);
    float NoH = clamp(dot(N, H), 0.0, 1.0);
    float VoH = clamp(dot(V, H), 0.0, 1.0);
    if(NoL > 0.0 && NoH > 0.0 && NoV > 0.0 && VoH > 0.0) {
      // geometry term
      float G = G_Smith(NoV, NoL, roughness);

      // Fresnel term
      vec3 F = fresnelSchlick(VoH, F0);

      vec2 uv = directionToSphericalEnvmap(L);
      vec3 radiance = textureLod(envMapImage, uv, envLevelSpec).rgb;
      result += radiance * F * G * VoH / (NoH * NoV);
    }
  }
  result = result / float(sampleCount);
  return result;
}

vec3 diffuseIBLReference(vec3 normal) {
  mat3 normalSpace = getNormalSpace(normal);

  vec3 result = vec3(0.0);
  uint sampleCount = uint(samplesDiff);
  float r = random2(tc);
  for(uint n = 1u; n <= sampleCount; n++) {
    //vec2 p = hammersley(n, N);
    vec2 p = mod(hammersley(n, sampleCount) + r, 1.0);
    float theta = asin(sqrt(p.y));
    float phi = 2.0 * PI * p.x;
    vec3 pos = vec3(sin(theta) * cos(phi), 
                    sin(theta) * sin(phi), 
                    cos(theta));
    vec3 posGlob = normalSpace * pos;
    vec2 uv = directionToSphericalEnvmap(posGlob);
    vec3 radiance = textureLod(envMapImage, uv, envLevelDiff).rgb;
    result +=  radiance;
  }
  result = result / float(sampleCount);
  return result;
}

void main() {
  vec3 normal = normalize(wfn);
  vec3 viewDir = normalize(cameraPos - vertPos);
  
  if(gsnMeshGroup == 0) {
    if(showBackground) {
      // color of envmap sphere
      outColor.rgb = texture(envMapImage, vec2(1.0-tc.x, tc.y)).rgb;
      outColor.a = 1.0;
    } else {
      discard;
    }
  } else {

    vec3 baseColor = pow(texture(baseColorTexture, tc).rgb, vec3(2.2));
    vec3 emission = pow(texture(emissionTexture, tc).rgb, vec3(2.2));;
    float roughness = texture(roughnessTexture, tc).r;
    float metallic = texture(metallicTexture, tc).r;
    
    // F0 for dielectics in range [0.0, 0.16] 
    // default FO is (0.16 * 0.5^2) = 0.04
    vec3 f0 = vec3(0.16 * (reflectance * reflectance)); 
    // in case of metals, baseColor contains F0
    f0 = mix(f0, baseColor, metallic);
    
    // compute diffuse and specular factors
    vec3 F = fresnelSchlick(max(dot(normal, viewDir), 0.0), f0);
    vec3 kS = F;
    vec3 kD = 1.0 - kS;
    kD *= 1.0 - metallic;    
    
    vec3 specular = specularIBLReference(f0, roughness, normal, viewDir); 
    vec3 diffuse = diffuseIBLReference(normal);
    vec3 color = emission + kD * baseColor * diffuse + specular;
    outColor.rgb = pow(color, vec3(1.0/2.2));
    outColor.a = 1.0;
  }
}

• Trotz Important Sampling ist das bisherige Verfahren noch zu langsam für Echtzeit-Anwendungen. Daher wird eine Lösung benötigt, die möglichst viele Anteile der Lösung vorberechnen kann.
• Dafür wird in [Unreal Engine 2013] die Split-Sum-Approximation vorgeschlagen
• Dieser Ansatz teilt die Summen-Formel des bisherigen Verfahrens in zwei Anteile auf die in zwei vorberechneten Umgebungstexturen $T_1$ und $T_2$ gespeichert werden:
  • (1) Pre-Filtered Environment Map $T_1$
  • (2) BRDF Integration Map $T_2$
  $\begin{align} &\frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \,L_i(\mathbf{l})\\ \approx &\underbrace{\frac{1}{N} \sum_{n=1}^{N} L_i(\mathbf{l})}_{T_1} \quad \cdot \quad \underbrace{\frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} }_{T_2}\\ \end{align}$

• Die WDF für das Important Sampling beider Anteile (1 u. 2) ist wie bisher, also
  $\phi_h = 2 \pi \,h_2(n)$ und $\theta_h = \arccos\left(\sqrt{\frac{1 - h_3(n)}{h_3(n) (\alpha^2-1) + 1} }\right)$
• Die Ergebnisse für verschiedene Werte für die Rauheit $\alpha = r_p^2$ mit $r_p$ im Intervall [0.0, 1.0] werden in den Mipmap-Leveln der Textur gespeichert:
  $r_p = \frac{\mbox{mipLevel}}{\mbox{mipCount}}$
• Die Blickrichtung ist zum Zeitpunkt der Vorberechung nicht bekannt. Daher wird von einer senkrechten Blickrichtung ausgegangen, d.h. $\mathbf{n} = \mathbf{v} = \mathbf{r}$

• Durch Einsetzen der Schlick Approximation für $\mathrm{F}$:
  $\begin{align} \mathrm{F}_{\tiny \mbox{Schlick}}(\mathbf{v}, \mathbf{h}) &= \mathrm{F}_0 + \left(1.0 − \mathrm{F}_0\right) \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5\\ &= \mathrm{F}_0 + \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5 - \mathrm{F}_0 \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5\\ &= \mathrm{F}_0 \left(1.0 - \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5\right)+ \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5 \end{align}$
  kann $\mathrm{F}_0$ aus der Summe gezogen werden:
  $\begin{align} T_2 =&\frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{F}(\mathbf{v},\mathbf{h})\,\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle}\\ =&\mathrm{F}_0 \underbrace{\frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \left(1.0 - \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5\right)}_{T_{2,r}}\\ &+ \underbrace{\frac{1}{N} \sum_{n=1}^{N} \frac{\mathrm{G}(\mathbf{l},\mathbf{v}) \,\langle\mathbf{v} \cdot \mathbf{h}\rangle}{\langle\mathbf{n} \cdot \mathbf{h}\rangle\,\,\langle\mathbf{n} \cdot \mathbf{v}\rangle} \left(1.0 − \langle\mathbf{v} \cdot \mathbf{h}\rangle \right)^5}_{T_{2,g}} \end{align}$

• Das Ergebnis der Vorberechnung kann im Rot- u. Grünkanal einer Textur gespeichert werden, die in x-Richtung mit $\langle\mathbf{n} \cdot \mathbf{v}\rangle$ und in y-Richtung mit der Rauheit $r_p$ parametrisiert ist (jeweils im Bereich [0.0, 1.0])

$T_{2,r}$

$T_{2,g}$

kombiniert in Rot- u. Grünkanal

$\langle\mathbf{n} \cdot \mathbf{v}\rangle \, \longrightarrow$

$r_p \, \longrightarrow$

• GSN Composer: SpheresIBL

Fragment Shader (Vertex Shader, wie bisher):

#version 300 es
precision highp float; 
precision highp int;
out vec4 outColor;

#define PI 3.1415926535897932384626433832795

in vec2 tc; // texture coordinate of pixel (interpolated)
in vec3 wfn; // fragment normal of pixel (interpolated)
in vec3 vertPos; // fragment vertex position (interpolated)
uniform sampler2D envmapImage; 
uniform sampler2D prefilteredEnvmap; 
uniform sampler2D brdfIntegrationMap; 
uniform sampler2D diffuseMap; 
uniform sampler2D baseColorTexture; 
uniform sampler2D roughnessTexture; // roughness texture
uniform sampler2D metallicTexture; // metallic texture
uniform sampler2D emissionTexture; // emission texture"
uniform float reflectance; // Fresnel reflectance
uniform bool showBackground;
uniform vec3 cameraPos; // camera position in global coordinate system
uniform int mipCount; // number of usable mipmap levels
uniform int gsnMeshGroup;

vec2 directionToSphericalEnvmap(vec3 dir) {
  float s = 1.0 - mod(1.0 / (2.0*PI) * atan(dir.y, dir.x), 1.0);
  float t = 1.0 / (PI) * acos(-dir.z);
  return vec2(s, t);
}

// adapted from "Real Shading in Unreal Engine 4", Brian Karis, Epic Games
vec3 specularIBL(vec3 F0 , float roughness, vec3 N, vec3 V) {
  float NoV = clamp(dot(N, V), 0.0, 1.0);
  vec3 R = reflect(-V, N);
  vec2 uv = directionToSphericalEnvmap(R);
  vec3 prefilteredColor = textureLod(prefilteredEnvmap, uv, 
                                     roughness*float(mipCount)).rgb;
  vec4 brdfIntegration = texture(brdfIntegrationMap, vec2(NoV, roughness));
  return prefilteredColor * ( F0 * brdfIntegration.x + brdfIntegration.y );
}

vec3 diffuseIBL(vec3 normal) {
  vec2 uv = directionToSphericalEnvmap(normal);
  return texture(diffuseMap, uv).rgb;
}

vec3 fresnelSchlick(float cosTheta, vec3 F0) {
  return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
} 

void main() {
  vec3 normal = normalize(wfn);
  vec3 viewDir = normalize(cameraPos - vertPos);
  
  if(gsnMeshGroup == 0) {
    if(showBackground) {
      // color of envmap sphere
      outColor.rgb = texture(envmapImage, vec2(1.0-tc.x, tc.y)).rgb;
      outColor.a = 1.0;
    } else {
      discard;
    }
  } else {

    vec3 baseColor = pow(texture(baseColorTexture, tc).rgb, vec3(2.2));
    vec3 emission = pow(texture(emissionTexture, tc).rgb, vec3(2.2));;
    float roughness = texture(roughnessTexture, tc).r;
    float metallic = texture(metallicTexture, tc).r;
    
    // F0 for dielectics in range [0.0, 0.16] 
    // default FO is (0.16 * 0.5^2) = 0.04
    vec3 f0 = vec3(0.16 * (reflectance * reflectance)); 
    // in case of metals, baseColor contains F0
    f0 = mix(f0, baseColor, metallic);
    
    // compute diffuse and specular factors
    vec3 F = fresnelSchlick(max(dot(normal, viewDir), 0.0), f0);
    vec3 kS = F;
    vec3 kD = 1.0 - kS;
    kD *= 1.0 - metallic;    
    
    vec3 specular = specularIBL(f0, roughness, normal, viewDir); 
    vec3 diffuse = diffuseIBL(normal);
    
    vec3 color = emission + kD * baseColor * diffuse + specular;
    outColor.rgb = pow(color, vec3(1.0/2.2));
    outColor.a = 1.0;
  }
}

• [Disney 2012] Brent Burley: Physically Based Shading at Disney SIGGRAPH 2012 SIGGRAPH 2012 Course: Practical Physically Based Shading in Film and Game Production
• [Unreal Engine 2013] Brian Karis: Real Shading in Unreal Engine 4, SIGGRAPH 2013 Course: Physically Based Shading in Theory and Practice
• [Frostbite 2014] S. Lagarde, C. de Rousiers: Moving Frostbite to PBR, SIGGRAPH 2014 Course: Physically Based Shading in Theory and Practice
• [Walter et al. 2007] B. Walter, S. R. Marschner, H. Li, K. E. Torrance: Microfacet Models for Refraction through Rough Surfaces. Eurographics Symposium on Rendering, 2007
• [Walter 2005] B. Walter: Notes on the Ward BRDF. Technical report PCG-05-06, Cornell University