Oren and Nayar's rough diffuse¶

The sun's light a rough diffuse surface sends toward the eye.

Oren and Nayar's model of a surface of V-cavities whose facets' slopes deviate by sigma (radians): in place of Lambert's diffuse / pi, diffuse / pi (A + B max(0, l . v - (n . l)(n . v)) / max(n . l, n . v)), A = 1 - 0.5 sigma^2 / (sigma^2 + 0.33), B = 0.45 sigma^2 / (sigma^2 + 0.09) (the qualitative model, with cos(phi) sin(alpha) tan(beta) written in the dot products), weighted by 1 - F(v . h) (lusine.render.schlick, the share the surface's interface does not reflect, as lusine.render.pbr weights its own diffuse), times (n . l) E sunlit: a radiance. sigma 0 is Lambert's exactly. n . v is held to 1e-4 at least. Source: Oren and Nayar, Generalization of Lambert's Reflectance Model (SIGGRAPH 1994), equation 30 and its form in the dot products (Fujii's rewriting). Checked against a reference in numpy, which has its own tests.

Source

Oren and Nayar, Generalization of Lambert's Reflectance Model (SIGGRAPH 1994), equation 30.

Proof

reference, invariant, limit.

oren_nayar.trait

namespace lusine.render;

vector<float<W/m2/sr>, 3> oren_nayar(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> sun_direction,
                                     vector<float, 3> f0, vector<float, 3> diffuse, float<rad> sigma,
                                     vector<float<W/m2>, 3> sun, float sunlit)
{
  let pi = 3.141592653589793;
  let s2 = sigma * sigma;
  let a = 1.0 - 0.5 * s2 / (s2 + 0.33);
  let b = 0.45 * s2 / (s2 + 0.09);
  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 vh = clamp(dot(eye, h), 0.0, 1.0);
  let across = dot(sun_direction, eye) - nl * nv;
  let on = a + b * (across > 0.0 ? across : 0.0) / (nl > nv ? nl : nv);
  let weight = vector<float, 3>(1.0, 1.0, 1.0) - schlick(f0 = f0, cosine = vh);
  return weight * diffuse * (on * nl / pi) * sun * sunlit / 1.0 sr;
}

The sun's light a rough diffuse surface sends toward the eye: Lambert's diffuse / pi scaled by A + B max(0, cos phi) sin alpha tan beta (in dot products), A and B from the facets' slope deviation sigma, weighted by 1 - F. 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/oren_nayar.trait'))

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

// lusine.render.oren_nayar: the sun's light a rough diffuse surface sends toward the eye. Oren and Nayar's model of a
// surface of V-cavities whose facets' slopes deviate by sigma (radians): in place of Lambert's diffuse / pi,
//   diffuse / pi (A + B max(0, l . v - (n . l)(n . v)) / max(n . l, n . v)),
//   A = 1 - 0.5 sigma^2 / (sigma^2 + 0.33),  B = 0.45 sigma^2 / (sigma^2 + 0.09)
// (the qualitative model, with cos(phi) sin(alpha) tan(beta) written in the dot products), weighted by 1 - F(v . h)
// (lusine.render.schlick, the share the surface's interface does not reflect, as lusine.render.pbr weights its own
// diffuse), times (n . l) E sunlit: a radiance. sigma 0 is Lambert's exactly. n . v is held to 1e-4 at least.
// Source: Oren and Nayar, Generalization of Lambert's Reflectance Model (SIGGRAPH 1994), equation 30 and its form in
// the dot products (Fujii's rewriting).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

vector<float<W/m2/sr>, 3> oren_nayar(vector<float, 3> normal, vector<float, 3> eye, vector<float, 3> sun_direction,
                                     vector<float, 3> f0, vector<float, 3> diffuse, float<rad> sigma,
                                     vector<float<W/m2>, 3> sun, float sunlit)
{
  let pi = 3.141592653589793;
  let s2 = sigma * sigma;
  let a = 1.0 - 0.5 * s2 / (s2 + 0.33);
  let b = 0.45 * s2 / (s2 + 0.09);
  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 vh = clamp(dot(eye, h), 0.0, 1.0);
  let across = dot(sun_direction, eye) - nl * nv;
  let on = a + b * (across > 0.0 ? across : 0.0) / (nl > nv ? nl : nv);
  let weight = vector<float, 3>(1.0, 1.0, 1.0) - schlick(f0 = f0, cosine = vh);
  return weight * diffuse * (on * nl / pi) * 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/oren_nayar.trait', 'oren_nayar'))
from lusine.kernels.render.oren_nayar import oren_nayar as reference

print(trait_run.law_program('render/oren_nayar.trait', 'oren_nayar')) from lusine.kernels.render.oren_nayar import oren_nayar 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> f0;
  host vector<float, 3> diffuse;
  host float<rad> sigma;
  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.oren_nayar(items.normal, items.eye, items.sun_direction, items.f0, items.diffuse, items.sigma, items.sun, items.sunlit); }

Reference: random lights and views¶

Random lights and views over a grass-like surface at sigma 0.5: the kernel's radiance within binary32's roundings of the float64 reference's.

In [4]:

def on(normal, eye, sun_direction, sigma, f0=0.03, diffuse=0.1):
    n = len(eye)
    return trait_run.law('render/oren_nayar.trait', 'oren_nayar', normal=normal, eye=eye, sun_direction=sun_direction,
                         f0=np.full((n, 3), f0), diffuse=np.full((n, 3), diffuse), sigma=sigma,
                         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 = on(up, eye, l, 0.5)
want = reference(up, eye, l, np.full((300, 3), float(F32(0.03))), np.full((300, 3), float(F32(0.1))), 0.5,
                 np.ones((300, 3)), np.ones(300))
rel = np.max(np.abs(got - want) / np.maximum(want, 1e-6))
assert rel <= 1e-5, rel
print(f'within {rel:.2g} of the reference, relatively')

def on(normal, eye, sun_direction, sigma, f0=0.03, diffuse=0.1): n = len(eye) return trait_run.law('render/oren_nayar.trait', 'oren_nayar', normal=normal, eye=eye, sun_direction=sun_direction, f0=np.full((n, 3), f0), diffuse=np.full((n, 3), diffuse), sigma=sigma, 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 = on(up, eye, l, 0.5) want = reference(up, eye, l, np.full((300, 3), float(F32(0.03))), np.full((300, 3), float(F32(0.1))), 0.5, np.ones((300, 3)), np.ones(300)) rel = np.max(np.abs(got - want) / np.maximum(want, 1e-6)) assert rel <= 1e-5, rel print(f'within {rel:.2g} of the reference, relatively')

within 2e-07 of the reference, relatively

Invariant: Lambert's at sigma 0¶

A smooth-faceted surface (sigma 0) is Lambert's diffuse weighted by 1 - F: the radiance is the same from every view where the light is fixed head on, but for the Fresnel weight.

In [5]:

head = np.tile([0.0, 0.0, 1.0], (300, 1))
lam = on(up, eye, head, 0.0)
from lusine.kernels.render.schlick import schlick
h = (eye + head) / np.linalg.norm(eye + head, axis=1, keepdims=True)
w = (1 - schlick(np.full((300, 3), 0.03), np.sum(eye * h, axis=1)))[:, 0]
assert np.allclose(lam[:, 0] / w, 0.1 / np.pi, rtol=1e-5)
print('sigma 0: Lambert, view-independent but for 1 - F')

head = np.tile([0.0, 0.0, 1.0], (300, 1)) lam = on(up, eye, head, 0.0) from lusine.kernels.render.schlick import schlick h = (eye + head) / np.linalg.norm(eye + head, axis=1, keepdims=True) w = (1 - schlick(np.full((300, 3), 0.03), np.sum(eye * h, axis=1)))[:, 0] assert np.allclose(lam[:, 0] / w, 0.1 / np.pi, rtol=1e-5) print('sigma 0: Lambert, view-independent but for 1 - F')

sigma 0: Lambert, view-independent but for 1 - F

Limit: backscatter, and the sun behind¶

Lit and seen from the same grazing angle a rough surface is brighter than Lambert's (the full moon's flat disc); with the sun behind it, it sends nothing.

In [6]:

g = np.array([[np.sin(1.3), 0.0, np.cos(1.3)]])
rough, smooth = on(up[:1], g, g, 0.5), on(up[:1], g, g, 0.0)
assert rough[0, 0] > 1.5 * smooth[0, 0]
assert np.all(on(up[:1], g, np.array([[0.6, 0.0, -0.8]]), 0.5) == 0.0)
print(f'backscatter {rough[0, 0] / smooth[0, 0]:.2f} times Lambert; none with the sun behind')

g = np.array([[np.sin(1.3), 0.0, np.cos(1.3)]]) rough, smooth = on(up[:1], g, g, 0.5), on(up[:1], g, g, 0.0) assert rough[0, 0] > 1.5 * smooth[0, 0] assert np.all(on(up[:1], g, np.array([[0.6, 0.0, -0.8]]), 0.5) == 0.0) print(f'backscatter {rough[0, 0] / smooth[0, 0]:.2f} times Lambert; none with the sun behind')

backscatter 1.93 times Lambert; none with the sun behind