Physically based shading¶
The radiance a surface sends toward the eye, lit by the sun and the sky.
The Cook-Torrance BRDF, with GGX's distribution (alpha = roughness^2, at least 0.002), the height-correlated Smith visibility and Schlick's Fresnel from the reflectance at normal incidence f0, and Lambert's diffuse weighted by 1 - F; the sun a direct light whose irradiance E arrives with the share sunlit of its disc in view, the sky a uniform radiance A: L = f(v, l) (n . l) E sunlit + A ((1 - F(n . v)) diffuse + F(n . v)). The sky's term is its integral over a uniform environment: Lambert's (1 - F) diffuse exactly, and the smooth surface's directional albedo F(n . v) for the specular. f0 is a surface's own, per channel: lusine.properties. fresnel_f0 of its n and k (a metal's, its element's; lusine.materials.optical); diffuse its colour, zero for a metal. Where f0 is zero in every channel there is no interface and nothing is reflected (Schlick's grazing term aside). Every square root is IEEE's and every power a product. n, v and l are unit vectors. Source: Walter et al., Microfacet Models for Refraction (EGSR 2007); Heitz, Understanding the Masking-Shadowing Function (JCGT 2014); Schlick, An Inexpensive BRDF Model (Eurographics 1994); Karis, Real Shading in Unreal Engine 4 (SIGGRAPH 2013). Checked against a reference in numpy, which has its own tests.
- Source
- Walter et al., Microfacet Models for Refraction (EGSR 2007); Heitz, Understanding the Masking-Shadowing Function (JCGT 2014); Schlick, An Inexpensive BRDF Model (Eurographics 1994); Karis, Real Shading in Unreal Engine 4 (SIGGRAPH 2013).
- Proof
- invariant, limit, reference.
pbr.trait
namespace lusine.render;
vector<float<W/m2/sr>, 3> pbr(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> f0,
vector<float, 3> diffuse, float roughness, vector<float, 3> sun_direction,
vector<float<W/m2>, 3> sun, float sunlit, vector<float<W/m2/sr>, 3> ambient)
{
let pi = 3.141592653589793;
let ones = vector<float, 3>(1.0, 1.0, 1.0);
let rough2 = roughness * roughness;
let alpha = rough2 < 0.002 ? 0.002 : rough2;
let a2 = alpha * alpha;
let nl = clamp(dot(normal, sun_direction), 0.0, 1.0);
let nv = clamp(dot(normal, eye), 0.0001, 1.0);
let h = normalize(eye + sun_direction);
let nh = clamp(dot(normal, h), 0.0, 1.0);
let vh = clamp(dot(eye, h), 0.0, 1.0);
let k = nh * nh * (a2 - 1.0) + 1.0;
let d = a2 / (pi * k * k);
let gv = nl * sqrt(nv * nv * (1.0 - a2) + a2);
let gl = nv * sqrt(nl * nl * (1.0 - a2) + a2);
let g = gv + gl;
let vis = 0.5 / (g < 0.000000000001 ? 0.000000000001 : g);
let interface = dot(f0, f0) != 0.0;
let mh = clamp(1.0 - vh, 0.0, 1.0);
let mh5 = mh * mh * mh * mh * mh;
let fh = interface ? f0 + (ones - f0) * mh5 : f0;
let mv = clamp(1.0 - nv, 0.0, 1.0);
let mv5 = mv * mv * mv * mv * mv;
let fv = interface ? f0 + (ones - f0) * mv5 : f0;
let brdf_cos = (fh * (d * vis) + (ones - fh) * diffuse / pi) * nl;
// the BRDF is per steradian, so f (n . l) E is a radiance: the irradiance over 1 sr
return brdf_cos * sun * sunlit / 1.0 sr + ambient * ((ones - fv) * diffuse + fv);
}
The radiance a surface sends toward the eye, lit by the sun and the sky. The Cook-Torrance BRDF (GGX's distribution, the height-correlated Smith visibility, Schlick's Fresnel from the surface's reflectance at normal incidence f0) with Lambert's diffuse weighted by 1 − F, the sun a direct light by the share of its disc in view, the sky a uniform radiance. A surface's f0 is the library's Fresnel kernel of its n and k (lusine.properties.fresnel_f0; a metal's, its element's), its diffuse colour the game's. lusine.render.table draws a table's shapes with it.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
The kernel¶
The law as the library states it.
print(trait_run.source('render/pbr.trait'))
// lusine.render.pbr: the radiance a surface sends toward the eye, lit by the sun and the sky. The Cook-Torrance BRDF,
// with GGX's distribution (alpha = roughness^2, at least 0.002), the height-correlated Smith visibility and Schlick's
// Fresnel from the reflectance at normal incidence f0, and Lambert's diffuse weighted by 1 - F; the sun a direct light
// whose irradiance E arrives with the share `sunlit` of its disc in view, the sky a uniform radiance A:
// L = f(v, l) (n. l) E sunlit + A ((1 - F(n. v)) diffuse + F(n. v)).
// The sky's term is its integral over a uniform environment: Lambert's (1 - F) diffuse exactly, and the smooth
// surface's directional albedo F(n. v) for the specular. f0 is a surface's own, per channel: lusine.properties.
// fresnel_f0 of its n and k (a metal's, its element's; lusine.materials.optical); diffuse its colour, zero for a metal.
// Where f0 is zero in every channel there is no interface and nothing is reflected (Schlick's grazing term aside).
// Every square root is IEEE's and every power a product. n, v and l are unit vectors.
// Source: Walter et al., Microfacet Models for Refraction (EGSR 2007); Heitz, Understanding the Masking-Shadowing
// Function (JCGT 2014); Schlick, An Inexpensive BRDF Model (Eurographics 1994); Karis, Real Shading in Unreal Engine 4
// (SIGGRAPH 2013).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
vector<float<W/m2/sr>, 3> pbr(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> f0,
vector<float, 3> diffuse, float roughness, vector<float, 3> sun_direction,
vector<float<W/m2>, 3> sun, float sunlit, vector<float<W/m2/sr>, 3> ambient)
{
let pi = 3.141592653589793;
let ones = vector<float, 3>(1.0, 1.0, 1.0);
let rough2 = roughness * roughness;
let alpha = rough2 < 0.002 ? 0.002: rough2;
let a2 = alpha * alpha;
let nl = clamp(dot(normal, sun_direction), 0.0, 1.0);
let nv = clamp(dot(normal, eye), 0.0001, 1.0);
let h = normalize(eye + sun_direction);
let nh = clamp(dot(normal, h), 0.0, 1.0);
let vh = clamp(dot(eye, h), 0.0, 1.0);
let k = nh * nh * (a2 - 1.0) + 1.0;
let d = a2 / (pi * k * k);
let gv = nl * sqrt(nv * nv * (1.0 - a2) + a2);
let gl = nv * sqrt(nl * nl * (1.0 - a2) + a2);
let g = gv + gl;
let vis = 0.5 /;
let interface = dot(f0, f0) != 0.0;
let mh = clamp(1.0 - vh, 0.0, 1.0);
let mh5 = mh * mh * mh * mh * mh;
let fh = interface ? f0 + (ones - f0) * mh5: f0;
let mv = clamp(1.0 - nv, 0.0, 1.0);
let mv5 = mv * mv * mv * mv * mv;
let fv = interface ? f0 + (ones - f0) * mv5: f0;
let brdf_cos = (fh * (d * vis) + (ones - fh) * diffuse / pi) * nl;
// the BRDF is per steradian, so f (n. l) E is a radiance: the irradiance over 1 sr
return brdf_cos * sun * sunlit / 1.0 sr + ambient * ((ones - fv) * diffuse + fv);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('render/pbr.trait', 'pbr'))
def pbr(**given):
return trait_run.law('render/pbr.trait', 'pbr', **given)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float, 3> normal;
host vector<float, 3> eye;
host vector<float, 3> f0;
host vector<float, 3> diffuse;
host float roughness;
host vector<float, 3> sun_direction;
host vector<float<W/m2>, 3> sun;
host float sunlit;
host vector<float<W/m2/sr>, 3> ambient;
vector<float<W/m2/sr>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.pbr(items.normal, items.eye, items.f0, items.diffuse, items.roughness, items.sun_direction, items.sun, items.sunlit, items.ambient); }
Invariant: the white furnace, and energy conserved¶
A white Lambertian surface with no interface under a uniform sky sends exactly the sky's radiance back, to binary32's rounding. A mirror's colour (f0 one, no diffuse) lit by unit irradiance from every direction of the hemisphere returns at most what it receives, and most of it: its directional albedo, summed over 600 directions of the hemisphere for three roughnesses and three angles of view (a smoother lobe than 0.5 wants a finer sum than a page runs).
up = [0.0, 0.0, 1.0]
white = pbr(normal=up, eye=up, f0=[0.0, 0.0, 0.0], diffuse=[1.0, 1.0, 1.0], roughness=1.0, sun_direction=up,
sun=[0.0, 0.0, 0.0], sunlit=1.0, ambient=[0.7, 0.7, 0.7])
assert np.all(np.abs(white - 0.7) <= 2 * EPS), white
nt, npp = 30, 20 # 600 directions, one run a case (the host's queue holds them)
th = (np.arange(nt) + 0.5) / nt * (np.pi / 2)
ph = (np.arange(npp) + 0.5) / npp * (2 * np.pi)
T, P = np.meshgrid(th, ph, indexing='ij')
dw = (np.pi / 2 / nt) * (2 * np.pi / npp) * np.sin(T).ravel()
l = np.stack([np.sin(T) * np.cos(P), np.sin(T) * np.sin(P), np.cos(T)], axis=-1).reshape(-1, 3)
m = len(l)
cases = [(r, a) for r in (0.5, 0.7, 0.9) for a in (0.1, 0.6, 1.2)]
returned = []
for r, a in cases:
L = pbr(normal=up, eye=[np.sin(a), 0.0, np.cos(a)], f0=np.ones((m, 3)), diffuse=[0.0, 0.0, 0.0], roughness=r,
sun_direction=l, sun=[1.0, 1.0, 1.0], sunlit=1.0, ambient=[0.0, 0.0, 0.0])
returned.append(float(np.sum(L[:, 0] * dw)))
assert all(e <= 1.005 for e in returned) and min(returned) > 0.4, returned
for (r, a), e in zip(cases, returned):
print(f'roughness {r}, view {a:.1f} rad: returns {e:.3f}')
roughness 0.5, view 0.1 rad: returns 0.916 roughness 0.5, view 0.6 rad: returns 0.900 roughness 0.5, view 1.2 rad: returns 0.772 roughness 0.7, view 0.1 rad: returns 0.698 roughness 0.7, view 0.6 rad: returns 0.691 roughness 0.7, view 1.2 rad: returns 0.729 roughness 0.9, view 0.1 rad: returns 0.421 roughness 0.9, view 0.6 rad: returns 0.449 roughness 0.9, view 1.2 rad: returns 0.594
Reference: Fresnel from the library's n¶
A black dielectric reflects its f0 of the sky head on, and nearly all of it at grazing incidence; the curve between is Schlick's. Its f0 is the library's own Fresnel kernel at k = 0, ((n − 1)/(n + 1))², water's 2% and glass's 4%, computed here by lusine.properties.fresnel_f0, the kernel the frame loop's table view gives pbr its reflectance from.
angles = np.linspace(0.0, np.pi / 2 - 1e-3, 91)
fig, ax = plt.subplots(figsize=(7, 3.2))
for n_ in (1.33, 1.5):
f0 = float(trait_run.law('properties/fresnel_f0.trait', 'fresnel_f0', n=[n_], k=[0.0])[0])
flat = [0.0, 0.0, 1.0]
F = pbr(normal=np.tile(flat, (len(angles), 1)), eye=np.stack([np.sin(angles), 0 * angles, np.cos(angles)], axis=1),
f0=[f0, f0, f0], diffuse=[0.0, 0.0, 0.0], roughness=0.0, sun_direction=flat, sun=[0.0, 0.0, 0.0],
sunlit=1.0, ambient=[1.0, 1.0, 1.0])[:, 0]
assert abs(F[0] - f0) <= 2 * EPS and F[-1] > 0.99 and np.all(np.diff(F) >= 0), (F[0], f0, F[-1])
ax.plot(np.degrees(angles), F, label=f'index {n_}: f0 {f0:.4f}')
ax.set_xlabel('angle from the normal (degrees)')
ax.set_ylabel('reflected share of the sky')
ax.legend()
plt.show()
Limit: in shadow under a black sky, and no interface¶
With the sun's disc out of view and nothing in the sky, a surface sends nothing; a black surface with no interface (f0 zero in every channel) reflects nothing even at grazing, where Schlick's term alone would say otherwise.
dark = pbr(normal=up, eye=up, f0=[0.04, 0.04, 0.04], diffuse=[0.5, 0.5, 0.5], roughness=0.5, sun_direction=up,
sun=[2.0, 2.0, 2.0], sunlit=0.0, ambient=[0.0, 0.0, 0.0])
none = pbr(normal=up, eye=[np.sin(1.5), 0.0, np.cos(1.5)], f0=[0.0, 0.0, 0.0], diffuse=[0.0, 0.0, 0.0],
roughness=0.5, sun_direction=up, sun=[1.0, 1.0, 1.0], sunlit=1.0, ambient=[1.0, 1.0, 1.0])
assert np.all(dark == 0.0) and np.all(none == 0.0), (dark, none)
print('dark, and nothing reflected')
dark, and nothing reflected