A fibrous surface's sheen¶

The sun's light a fibrous surface (grass, cloth) scatters toward the eye at grazing angles, a sheen lobe: its colour `sheen` times the surface's diffuse colour, Estevez and Kulla's "Charlie" distribution D = (2 + 1/a) (1 - (n .

h)^2)^(1/(2a)) / (2 pi), a = roughness^2, and Neubelt and Pettineo's visibility V = 1 / (4 (n . l + n . v - (n . l)(n . v))), times (n . l) E sunlit: a radiance. None where the sun is behind the surface (n . l <= 0). n . v is held to 1e-4 at least; the power is the library's pow. Source: Estevez and Kulla, Production Friendly Microfacet Sheen BRDF (Sony Pictures Imageworks, 2017); Neubelt and Pettineo, Crafting a Next-Gen Material Pipeline for The Order: 1886 (SIGGRAPH 2013). Its alpha, roughness^2, is floored at 0.002 as pbr's is, so a sheen of roughness 0 is a sharp sheen and never a division by zero. Checked against a reference in numpy, which has its own tests.

Source

Estevez and Kulla, Production Friendly Microfacet Sheen BRDF (2017); Neubelt and Pettineo, Crafting a Next-Gen Material Pipeline for The Order: 1886 (SIGGRAPH 2013).

Proof

reference, invariant, limit.

sheen.trait

namespace lusine.render;

vector<float<W/m2/sr>, 3> sheen(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> sun_direction,
                                vector<float, 3> diffuse, float sheen, float roughness, vector<float<W/m2>, 3> sun,
                                float sunlit)
{
  let pi = 3.141592653589793;
  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 a = roughness * roughness;
  let inv = 1.0 / (a < 0.002 ? 0.002 : a);
  let s = 1.0 - nh * nh;
  let d = (2.0 + inv) * pow(s > 0.0 ? s : 0.0, 0.5 * inv) / (2.0 * pi);
  let v = 1.0 / (4.0 * (nl + nv - nl * nv));
  let lit = nl > 0.0 ? d * v * nl : 0.0;
  return diffuse * (sheen * lit) * sun * sunlit / 1.0 sr;
}

The sun's light grass or cloth scatters toward the eye at grazing angles: Estevez and Kulla's Charlie distribution with Neubelt and Pettineo's visibility, its colour a share of the surface's diffuse. One of the kernels lusine.render.table draws a table's shapes with.

In [1]:

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
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)

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 F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24 rng = np.random.default_rng(2026)

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('render/sheen.trait'))

print(trait_run.source('render/sheen.trait'))

// lusine.render.sheen: the sun's light a fibrous surface (grass, cloth) scatters toward the eye at grazing angles, a
// sheen lobe: its colour `sheen` times the surface's diffuse colour, Estevez and Kulla's "Charlie" distribution
//   D = (2 + 1/a) (1 - (n . h)^2)^(1/(2a)) / (2 pi),  a = roughness^2,
// and Neubelt and Pettineo's visibility V = 1 / (4 (n . l + n . v - (n . l)(n . v))), times (n . l) E sunlit: a
// radiance. None where the sun is behind the surface (n . l <= 0). n . v is held to 1e-4 at least; the power is the
// library's pow.
// Source: Estevez and Kulla, Production Friendly Microfacet Sheen BRDF (Sony Pictures Imageworks, 2017); Neubelt and
// Pettineo, Crafting a Next-Gen Material Pipeline for The Order: 1886 (SIGGRAPH 2013).
// Its alpha, roughness^2, is floored at 0.002 as pbr's is, so a sheen of roughness 0 is a sharp sheen and never a
// division by zero.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

vector<float<W/m2/sr>, 3> sheen(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> sun_direction,
                                vector<float, 3> diffuse, float sheen, float roughness, vector<float<W/m2>, 3> sun,
                                float sunlit)
{
  let pi = 3.141592653589793;
  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 a = roughness * roughness;
  let inv = 1.0 / (a < 0.002 ? 0.002 : a);
  let s = 1.0 - nh * nh;
  let d = (2.0 + inv) * pow(s > 0.0 ? s : 0.0, 0.5 * inv) / (2.0 * pi);
  let v = 1.0 / (4.0 * (nl + nv - nl * nv));
  let lit = nl > 0.0 ? d * v * nl : 0.0;
  return diffuse * (sheen * lit) * sun * sunlit / 1.0 sr;
}

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.

In [3]:

print(trait_run.law_program('render/sheen.trait', 'sheen'))
from lusine.kernels.render.sheen import sheen as reference

print(trait_run.law_program('render/sheen.trait', 'sheen')) from lusine.kernels.render.sheen import sheen as reference

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host vector<float, 3> normal;
  host vector<float, 3> eye;
  host vector<float, 3> sun_direction;
  host vector<float, 3> diffuse;
  host float sheen;
  host float roughness;
  host vector<float<W/m2>, 3> sun;
  host float sunlit;
  vector<float<W/m2/sr>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.render.sheen(items.normal, items.eye, items.sun_direction, items.diffuse, items.sheen, items.roughness, items.sun, items.sunlit); }

Reference: random lights and views¶

Random lights and views at roughness 0.4 (a lawn's): the kernel's radiance within binary32's roundings of the float64 reference's (its power the library's fixed sequence).

In [4]:

def lobe(normal, eye, sun_direction, roughness=0.4, sheen=0.5):
    n = len(eye)
    return trait_run.law('render/sheen.trait', 'sheen', normal=normal, eye=eye, sun_direction=sun_direction,
                         diffuse=np.full((n, 3), 0.1), sheen=sheen, roughness=roughness, sun=np.ones((n, 3)),
                         sunlit=np.ones(n))

up = np.tile([0.0, 0.0, 1.0], (300, 1))
v = rng.normal(size=(300, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True)
v[:, 2] = np.abs(v[:, 2]); eye = F32(v).astype(float)
v = rng.normal(size=(300, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True)
v[:, 2] = np.abs(v[:, 2]); l = F32(v).astype(float)
got = lobe(up, eye, l)
want = reference(up, eye, l, np.full((300, 3), float(F32(0.1))), 0.5, float(F32(0.4)), np.ones((300, 3)), np.ones(300))
big = want > 1e-6
rel = np.max(np.abs(got - want)[big] / want[big])
assert rel <= 1e-4, rel
print(f'within {rel:.2g} of the reference, relatively')

def lobe(normal, eye, sun_direction, roughness=0.4, sheen=0.5): n = len(eye) return trait_run.law('render/sheen.trait', 'sheen', normal=normal, eye=eye, sun_direction=sun_direction, diffuse=np.full((n, 3), 0.1), sheen=sheen, roughness=roughness, sun=np.ones((n, 3)), sunlit=np.ones(n)) up = np.tile([0.0, 0.0, 1.0], (300, 1)) v = rng.normal(size=(300, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True) v[:, 2] = np.abs(v[:, 2]); eye = F32(v).astype(float) v = rng.normal(size=(300, 3)); v /= np.linalg.norm(v, axis=1, keepdims=True) v[:, 2] = np.abs(v[:, 2]); l = F32(v).astype(float) got = lobe(up, eye, l) want = reference(up, eye, l, np.full((300, 3), float(F32(0.1))), 0.5, float(F32(0.4)), np.ones((300, 3)), np.ones(300)) big = want > 1e-6 rel = np.max(np.abs(got - want)[big] / want[big]) assert rel <= 1e-4, rel print(f'within {rel:.2g} of the reference, relatively')

within 5.5e-06 of the reference, relatively

Invariant: strongest at grazing¶

Seen and lit along the same direction, the sheen grows as that direction grazes the surface: the fibres' sides catch the light.

In [5]:

th = np.linspace(0.1, 1.45, 28)
dirs = np.stack([np.sin(th), np.zeros_like(th), np.cos(th)], axis=1)
s = lobe(up[:28], dirs, dirs)[:, 0]
assert np.all(np.diff(s) > 0)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(np.degrees(th), s)
ax.set_xlabel('angle from the normal (degrees)'); ax.set_ylabel('radiance (W/(m2 sr))')
plt.show()

th = np.linspace(0.1, 1.45, 28) dirs = np.stack([np.sin(th), np.zeros_like(th), np.cos(th)], axis=1) s = lobe(up[:28], dirs, dirs)[:, 0] assert np.all(np.diff(s) > 0) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(np.degrees(th), s) ax.set_xlabel('angle from the normal (degrees)'); ax.set_ylabel('radiance (W/(m2 sr))') plt.show()

Limit: head on, behind, and none¶

Lit and seen head on there is no sheen (the distribution is 0 at n . h = 1); nor with the sun behind the surface, nor at sheen 0.

In [6]:

out = np.concatenate([lobe(up[:1], up[:1], up[:1]), lobe(up[:1], up[:1], np.array([[0.6, 0.0, -0.8]])),
                      lobe(up[:1], dirs[-1:], dirs[-1:], sheen=0.0)])
assert np.all(out == 0.0)
print('none head on, behind, or at sheen 0')

out = np.concatenate([lobe(up[:1], up[:1], up[:1]), lobe(up[:1], up[:1], np.array([[0.6, 0.0, -0.8]])), lobe(up[:1], dirs[-1:], dirs[-1:], sheen=0.0)]) assert np.all(out == 0.0) print('none head on, behind, or at sheen 0')

none head on, behind, or at sheen 0