Reflectance Models

A reflectance model says how much of the light arriving at a surface from one direction leaves it towards another. The quantity it describes is the bidirectional reflectance distribution function (BRDF), and a reflectance model is a formula for the BRDF with a few parameters: an albedo, a highlight exponent, a roughness, an index of refraction. Renderers evaluate one, and vision methods that read shape, material or lighting from pixel values invert one; the methods linked at the end of this article mostly invert the simplest, the Lambertian model.

Image formation uses the BRDF in its reflection equation and describes, in a paragraph each, the split into diffuse and specular reflection, Phong’s exponent and ambient term, and Shafer’s dichromatic model. Radiometry and photometry defines radiance and irradiance and derives the Lambertian BRDF $\rho/\pi$. This article covers the models themselves: the definition of the BRDF and the three properties that make one physically plausible (positivity, reciprocity and energy conservation), the empirical Phong and Blinn–Phong models and what their normalization constants do, and microfacet models, with the distribution of normals $D$, the masking–shadowing term $G$ and the Fresnel term $F$ that Walter et al. (2007) give for the GGX distribution. The example tests six models for reciprocity, for energy and for Heitz’s normalization test of the masking function, and renders their highlights on a sphere.

Definition

The bidirectional reflectance distribution function of a surface element is the ratio of the radiance it reflects in one direction to the irradiance it receives from another. Nicodemus et al. (1977, Eq. 9, p. 5) define it as

$$f_r(\theta_i, \phi_i; \theta_r, \phi_r) = \frac{dL_r(\theta_i, \phi_i; \theta_r, \phi_r; E_i)}{dE_i(\theta_i, \phi_i)} = \frac{dL_r}{L_i(\theta_i, \phi_i)\cos\theta_i\, d\omega_i} \quad [\text{sr}^{-1}],$$

where the polar angles $\theta$ are measured from the surface normal and the azimuths $\phi$ in the surface plane. The BRDF is a function of two directions, four angles in all (and of wavelength), and its unit is the inverse steradian: it is a density over directions, not a fraction, so it can exceed 1. A mirror’s BRDF is a Dirac delta in the mirror direction (Walter et al., 2007, Eq. 9 writes a smooth specular surface this way). Visionbook (Torralba, Isola and Freeman, 2024, §5.2) introduces the same function as the one that turns incoming ray power into outgoing ray power for a surface normal, wavelength and pair of directions.

A reflectance model is a parametric BRDF. Some are empirical, chosen to look right (Phong, Blinn–Phong); some are derived from a model of the surface, such as a random field of tiny mirrors (microfacet models).

Intuition

Shine a narrow beam on a surface and walk around it with a light meter that measures radiance. An ideal matte surface reads the same from every direction: its BRDF is a constant. An ideal glossy surface reads almost nothing except near the mirror direction, where it reads a great deal: its BRDF is a sharp lobe around the mirror direction. The models below combine or interpolate between these, with a lobe whose width is set by an exponent or a roughness.

The BRDF is the light meter’s reading divided by how much light the beam delivers per unit area, which depends on how obliquely the beam strikes. Dividing by the irradiance rather than the beam’s radiance is what makes the BRDF a property of the surface alone, and it is also what makes a BRDF that is constant over directions describe a surface that looks equally bright from everywhere.

Microfacet models explain the lobe. A rough surface is pictured as many tiny flat mirrors tilted at random. Light from $\mathbf{i}$ leaves towards $\mathbf{o}$ only from the facets that would mirror one into the other, those whose normal is the half-vector $\mathbf{h}$, halfway between $\mathbf{i}$ and $\mathbf{o}$. The brightness in a direction is then proportional to how many facets point along $\mathbf{h}$, reduced by the facets hidden from the light or from the viewer by other facets, and scaled by how much light a single smooth facet reflects at that angle.

Formal Definition

The reflection equation and the albedo

From the definition, a surface lit from all directions reflects towards $\mathbf{o}$ the radiance (our statement, integrating Nicodemus et al., Eq. 6, over incoming directions)

$$L_r(\mathbf{o}) = \int_{\Omega_i} f_r(\mathbf{i}, \mathbf{o})\, L_i(\mathbf{i}) \cos\theta_i\, d\omega_i .$$

For point sources the integral becomes the sum in image formation’s reflection equation (Szeliski’s Eq. 2.86, as that article cites it). The fraction of the flux arriving from one direction that is reflected into the whole hemisphere is the directional–hemispherical reflectance (Nicodemus et al., Table 1, item 3, p. 11),

$$\rho(\theta_i, \phi_i; 2\pi) = \int_{2\pi} f_r(\theta_i, \phi_i; \theta_r, \phi_r) \cos\theta_r\, d\omega_r ,$$

the integral of the BRDF over projected solid angle $d\Omega_r = \cos\theta_r\, d\omega_r$ (note to Table 1). This article calls it the directional albedo and writes it $a(\theta_i)$.

Three properties

• Positivity. $f_r \ge 0$.
• Helmholtz reciprocity. Swapping the two directions does not change the BRDF: $f_r(\theta_1, \phi_1; \theta_2, \phi_2) = f_r(\theta_2, \phi_2; \theta_1, \phi_1)$ (Nicodemus et al., Eq. B7, p. 40). They state that it holds in the absence of polarization and magnetic fields, and that it lets either direction of a pair be the incident one.
• Energy conservation. A surface cannot reflect more flux than it receives, so $a(\theta_i) \le 1$ for every incident direction (our statement of what the directional–hemispherical reflectance, a ratio of reflected to incident flux, must satisfy). For a perfect, lossless reflector $a = 1$.

A BRDF that depends on the azimuths only through their difference $\phi_r - \phi_i$, so that rotating the surface about its normal changes nothing, is called isotropic; it is a function of three angles (our statement). Walter et al. note that anisotropic variants of their distributions exist; all models in this article are isotropic.

Reciprocity is a property of the BRDF. A formula for reflected radiance that is not of the form $f_r E_i$ with a reciprocal $f_r$ need not be reciprocal, and the Phong model below is the standard example.

Lambertian reflection

A perfectly diffuse surface reflects the same radiance in every direction whatever its lighting, which requires a constant BRDF; Nicodemus et al. (Appendix C, Eq. C4, p. 43) show that its value is the diffuse reflectance over $\pi$:

$$f_L = \frac{\rho}{\pi}.$$

The $\pi$ is the projected solid angle of the hemisphere, so $a(\theta_i) = \rho$ for every $\theta_i$. Radiometry and photometry derives this and notes that visionbook’s Lambertian law (its Eq. 5.2), written in units of ray power, has no $\pi$. Visionbook names a synthetic material, spectralon, as the most nearly Lambertian. Real matte surfaces are not Lambertian at oblique angles: Heitz (2014, §3.4) writes the BRDF of a surface made of diffuse microfacets and notes that it has no analytical solution, and that Oren and Nayar (1994) fitted it for a spherical Gaussian distribution of normals with V-cavity masking.

Phong and normalized Phong

Visionbook (§5.2.2) gives the specular part of Phong’s model (Phong, 1975, not read here) as outgoing ray power

$$\ell_{\text{spec}} = k_s\, (\mathbf{r} \cdot \mathbf{q})^{\alpha}\, \ell_{\text{in}}, \qquad \mathbf{r} = 2(\mathbf{p} \cdot \mathbf{n})\mathbf{n} - \mathbf{p},$$

with $\mathbf{p}$ towards the light, $\mathbf{q}$ the outgoing direction and $\mathbf{r}$ the mirror direction. The formula is written for light, not as a BRDF, and it has no $\cos\theta_i$. Reading $\ell_{\text{in}}$ as the irradiance divided by $\cos\theta_i$, the BRDF that reproduces it is (our conversion, used as model P below, with $\mathbf{i}$ for $\mathbf{p}$, $\mathbf{o}$ for $\mathbf{q}$ and $n$ for $\alpha$)

$$f_P(\mathbf{i}, \mathbf{o}) = \frac{k_s \max(0, \mathbf{r}(\mathbf{i}) \cdot \mathbf{o})^n}{\mathbf{n} \cdot \mathbf{i}} .$$

Because $\mathbf{r}(\mathbf{i}) \cdot \mathbf{o} = 2(\mathbf{n}\cdot\mathbf{i})(\mathbf{n}\cdot\mathbf{o}) - \mathbf{i}\cdot\mathbf{o}$ is symmetric in the two directions, $f_P(\mathbf{i}, \mathbf{o})/f_P(\mathbf{o}, \mathbf{i}) = \cos\theta_o/\cos\theta_i$ (our derivation): this form is not reciprocal. At normal incidence its albedo is $k_s \int \cos^{n+1}\theta\, d\omega = 2\pi k_s/(n+2)$ (our derivation), so with $k_s$ fixed a sharper highlight reflects less energy, and for $n < 2\pi - 2$ a surface with $k_s = 1$ reflects more than it receives.

The normalized Phong BRDF drops the division and scales the lobe so that its albedo at normal incidence is $k_s$ (our derivation):

$$f_{Pn}(\mathbf{i}, \mathbf{o}) = k_s\, \frac{n + 2}{2\pi} \max(0, \mathbf{r}(\mathbf{i}) \cdot \mathbf{o})^n .$$

It is reciprocal by the symmetry above. The constant $(n+2)/(2\pi)$ is the same one Walter et al. (2007, Eq. 30) use for the Phong distribution of microfacet normals, $D(\mathbf{m}) = \chi^+(\mathbf{m}\cdot\mathbf{n})\,(\alpha_p + 2)/(2\pi)\,(\cos\theta_m)^{\alpha_p}$, which they normalize so that its projected area is one. “Phong” thus names two different things here, the reflection model and this distribution of normals; Heitz (2014, footnote 6, p. 93) separates them, noting that the Phong BRDF is not a microfacet model.

Blinn–Phong

The Blinn–Phong model replaces the angle between $\mathbf{r}$ and $\mathbf{o}$ with the angle between the normal and the half-vector:

$$f_{BP}(\mathbf{i}, \mathbf{o}) = c\, \max(0, \mathbf{n} \cdot \mathbf{h})^m, \qquad \mathbf{h} = \frac{\mathbf{i} + \mathbf{o}}{\|\mathbf{i} + \mathbf{o}\|}.$$

It is reciprocal because $\mathbf{h}$ is symmetric in $\mathbf{i}$ and $\mathbf{o}$. Heitz (2014, footnote 6) notes that the original Blinn–Phong BRDF uses a non-normalized distribution. None of the sources read for this article gives a constant $c$, so the example prints two: $c = (m+2)/(2\pi)$, which is Walter et al.’s Eq. 30 applied to $\mathbf{h}$, and $c = (m+2)/(8\pi)$, that constant divided by 4 (our derivation). The factor 4 is the Jacobian between the outgoing direction and the half-vector: Walter et al. (2007, Eq. 14) give $\partial\omega_h/\partial\omega_o = 1/(4|\mathbf{o}\cdot\mathbf{h}|)$, which is $1/4$ at normal incidence, because $\mathbf{o}$ moves twice as fast as $\mathbf{h}$.

The two lobes are related. In the plane of incidence a normal tilted by $\delta$ from $\mathbf{h}$ turns $\mathbf{r}$ by $2\delta$ away from $\mathbf{o}$, so $(\mathbf{r}\cdot\mathbf{o})^n = \cos^n 2\delta \approx 1 - 2n\delta^2$, while $(\mathbf{n}\cdot\mathbf{h})^m = \cos^m\delta \approx 1 - m\delta^2/2$; the two match for $m = 4n$. Out of the plane, a tilt by $\beta$ gives $\mathbf{r}\cdot\mathbf{o} \approx 1 - 2\cos^2\!\phi\, \beta^2$, where $2\phi$ is the angle between light and viewer, against $1 - m\beta^2/2$ for Blinn–Phong. With $m = 4n$ the Phong lobe is therefore wider out of the plane by a factor $1/\cos\phi$ (all our derivation). The two models agree when the light is at the camera and differ more as the angle between them grows.

Microfacet models

Microfacet models were introduced to graphics by Cook and Torrance (1982), based on earlier work in optics by Torrance and Sparrow (1967), as Walter et al. (2007, §2) report. Walter et al. review the theory in a form that this section follows. The surface is replaced by a macrosurface with normal $\mathbf{n}$ whose BRDF reproduces the average scattering of a detailed microsurface. Geometric optics is typically assumed and only single scattering is modelled; wave effects and light that hits the surface more than once are left out or handled separately (§3).

Distribution of normals. $D(\mathbf{m})\,d\omega_m\,dA$ is the microsurface area whose normals lie in $d\omega_m$, over a macrosurface area $dA$. Walter et al. call $D$ a density with units of 1/sr; Heitz (2014, footnote 1) states that, as a distribution of normals per square metre of geometric surface, it is measured in m²/sr, not 1/sr. Its projected area must equal the macrosurface’s (Walter et al., Eqs. 3–4):

$$\int D(\mathbf{m})\,(\mathbf{v}\cdot\mathbf{m})\,d\omega_m = \mathbf{v}\cdot\mathbf{n} \quad\text{for any } \mathbf{v}, \qquad \int D(\mathbf{m})\,(\mathbf{n}\cdot\mathbf{m})\,d\omega_m = 1 .$$

Walter et al. give three distributions. The Beckmann distribution (Eq. 25) arises from Gaussian roughness; the Phong distribution (Eq. 30, above) is empirical, and with $\alpha_p = 2\alpha_b^{-2} - 2$ it is very similar to Beckmann’s, especially for narrow widths (their Fig. 8 shows $\alpha_b = 0.2$ against $\alpha_p = 48$). The GGX distribution is the one they introduced, “to better match our measured data for transmission” (§5.1):

$$D(\mathbf{m}) = \frac{\alpha_g^2\, \chi^+(\mathbf{m}\cdot\mathbf{n})}{\pi \cos^4\theta_m\, (\alpha_g^2 + \tan^2\theta_m)^2} \qquad (\text{Eq. 33}),$$

where $\chi^+(a)$ is 1 for $a > 0$ and 0 otherwise. In their Figure 8, GGX has a narrower peak and stronger tails than Beckmann at the same width. The width parameter $\alpha_g$ is GGX’s own; Walter et al. do not relate it to a “roughness” slider, and a renderer’s roughness parameter may be mapped to $\alpha_g$ in other ways.

Masking and shadowing. $G(\mathbf{i}, \mathbf{o}, \mathbf{m})$ is the fraction of facets with normal $\mathbf{m}$ visible from both directions. Walter et al. require it to lie between 0 and 1, to be symmetric in $\mathbf{i}$ and $\mathbf{o}$, and to be zero for facets seen from behind (Eqs. 5–7). They use Smith’s approximation, separable into two one-directional terms (Eq. 23, after Smith, 1967), $G \approx G_1(\mathbf{i}, \mathbf{m})\,G_1(\mathbf{o}, \mathbf{m})$, with, for GGX,

$$G_1(\mathbf{v}, \mathbf{m}) = \chi^+\!\left(\frac{\mathbf{v}\cdot\mathbf{m}}{\mathbf{v}\cdot\mathbf{n}}\right) \frac{2}{1 + \sqrt{1 + \alpha_g^2 \tan^2\theta_v}} \qquad (\text{Eq. 34}).$$

They recommend against the Cook–Torrance $G$, built on a one-dimensional model of parallel grooves, because of its first-derivative discontinuities (§5.1).

Fresnel. A smooth facet reflects the fraction $F$ given by the Fresnel equations. For dielectrics and unpolarized light Walter et al. (Eq. 22) write

$$F = \frac{1}{2}\frac{(g - c)^2}{(g + c)^2}\left(1 + \frac{(c(g + c) - 1)^2}{(c(g - c) + 1)^2}\right), \qquad g = \sqrt{\eta_t^2/\eta_i^2 - 1 + c^2}, \quad c = |\mathbf{i}\cdot\mathbf{m}| .$$

It is small at normal incidence, 0.04 for glass with $\eta_t = 1.5$, and rises to 1 at grazing angles (§5.1). Walter et al. (§2) report that Schlick (1994) gave a cheaper approximation to the Fresnel formula; they do not use it.

The BRDF. Integrating over the facets and keeping only the one whose normal is $\mathbf{h}$ gives (Walter et al., Eq. 20)

$$f_r(\mathbf{i}, \mathbf{o}) = \frac{F(\mathbf{i}, \mathbf{h})\, G(\mathbf{i}, \mathbf{o}, \mathbf{h})\, D(\mathbf{h})}{4\,|\mathbf{i}\cdot\mathbf{n}|\,|\mathbf{o}\cdot\mathbf{n}|} .$$

The 4 is the Jacobian above. Walter et al. note that the Cook–Torrance BSDF has $\pi$ in its place, with a different normalization of $D$, and that later papers agree with 4 (§5). The same paper extends the model to transmission through rough glass, which is its main subject.

The weak white furnace test

A BRDF that models only the first bounce of light off the microsurface cannot reflect all the light, even for a perfect mirror material, because the shadowing term discards rays that would hit another facet and bounce again. Heitz (2014, §3.5) states the consequence: the White Furnace Test, $\int \rho(\mathbf{o}, \mathbf{i})\,|\mathbf{n}\cdot\mathbf{i}|\,d\omega_i = 1$ with $F = 1$, is not met by common microfacet BRDFs, which model a single scattering event. He proposes a weaker test that they must meet: with Fresnel set to 1 and shadowing removed, so that only the masking of the outgoing direction remains, the distribution of reflected rays after the first bounce must be normalized,

$$\int_{\Omega_i} \frac{G_1(\mathbf{o}, \mathbf{h})\, D(\mathbf{h})}{4\,|\mathbf{n}\cdot\mathbf{o}|}\, d\omega_i = 1 \qquad (\text{Heitz, Eq. 36}).$$

It holds only for a masking function that satisfies his Eq. 14, $\cos\theta_o = \int G_1(\mathbf{o}, \mathbf{m})\,\max(0, \mathbf{o}\cdot\mathbf{m})\,D(\mathbf{m})\,d\omega_m$: the visible facets, projected onto the view direction, must cover exactly the projected macrosurface. Heitz derives Smith’s function and the V-cavity function from this condition and finds both correct for their own microsurface profiles (§4). The integral runs over the whole sphere of $\mathbf{i}$, as in the code of his Appendix C, with $D = 0$ where $\mathbf{h}$ points below the surface.

Without masking, $G_1 = 1$, the same integral is $\int \max(0, \mathbf{o}\cdot\mathbf{m})\,D(\mathbf{m})\,d\omega_m/\cos\theta_o$. Using Heitz’s Eq. 14 with Walter et al.’s Eq. 34, this equals $1/G_1(\mathbf{o}) = 1 + \Lambda$, with $\Lambda = (\sqrt{1 + \alpha_g^2\tan^2\theta_o} - 1)/2$ for GGX (our derivation). It exceeds 1 wherever any facet faces away from the viewer, and it depends on $\alpha_g$ and $\theta_o$ only through $\alpha_g \tan\theta_o$.

Properties

• A BRDF is a density, not a fraction. Its values can be large; what is bounded by 1 is its integral over projected solid angle, the albedo.
• Reciprocity constrains the form. Phong’s model, written as reflected light, is not reciprocal; normalized Phong, Blinn–Phong and microfacet BRDFs are.
• Normalization decides energy. An unnormalized lobe reflects energy that depends on its exponent; normalized Phong has the constant $k_s$ in front as its albedo at normal incidence, so that sharpness and brightness can be set separately. Neither of the example’s Blinn–Phong constants gives exactly that.
• Microfacet BRDFs are built to be plausible. $D$ is normalized by projected area (Walter et al., Eq. 4), $G_1$ by the visible projected area (Heitz, Eq. 14), and $F$ comes from the Fresnel equations. In our reading the form is reciprocal whenever $G$ is symmetric, since $\mathbf{i}\cdot\mathbf{h} = \mathbf{o}\cdot\mathbf{h}$ makes $F$ symmetric too.
• Single scattering loses energy. With $F = 1$ a microfacet BRDF with shadowing reflects less than all the light (Heitz, 2014, §3.5; in §7 he proposes a multiple-scattering term as future work). In the example the loss grows with roughness at every printed angle.
• Highlights move with the viewer. A specular lobe’s peak and shape depend on the light–view geometry; a Lambertian surface’s radiance does not depend on the view at all. Brightness constancy lists specular highlights among its failures, because they move with the viewpoint and the light, not with the surface.

Examples

Six reflectance models under three tests and one image

The script implements six model families, all isotropic and in sr⁻¹, with the normal $\mathbf{n} = +z$ for the flat-surface tests. The directions $\mathbf{i}$ (towards the light) and $\mathbf{o}$ (towards the viewer) both point away from the surface, the convention of Walter et al. (Fig. 5); $\mathbf{r}(\mathbf{i}) = 2(\mathbf{i}\cdot\mathbf{n})\mathbf{n} - \mathbf{i}$; $\mathbf{h} = (\mathbf{i} + \mathbf{o})/\|\mathbf{i} + \mathbf{o}\|$.

• L: Lambertian, $\rho/\pi$ with $\rho = 0.8$.
• P1, P20: Phong as a BRDF, $k_s(\mathbf{r}\cdot\mathbf{o})^n/(\mathbf{n}\cdot\mathbf{i})$, with $k_s = 1$ and $n = 1$ or 20 (our conversion of visionbook’s form).
• Pn1, Pn20: normalized Phong, $k_s (n+2)/(2\pi)(\mathbf{r}\cdot\mathbf{o})^n$, $k_s = 1$.
• BP80/2pi, BP80/8pi: Blinn–Phong with $m = 80 = 4 \times 20$, with $c = 82/(2\pi)$ (Walter et al.’s Eq. 30 applied to $\mathbf{h}$) and $c = 82/(8\pi)$ (our derivation).
• M0.1, M0.3, M0.6: Walter et al.’s Eq. 20 with GGX $D$ (Eq. 33), Smith $G = G_1 G_1$ (Eqs. 23, 34) and the exact Fresnel term (Eq. 22) from air into $\eta = 1.5$, at $\alpha_g = 0.1$, 0.3 and 0.6.
• M F=1: the same with $F = 1$, a surface of perfect mirrors.

The tests are:

• H1, reciprocity. The largest relative difference $|f(\mathbf{i},\mathbf{o}) - f(\mathbf{o},\mathbf{i})|/(f(\mathbf{i},\mathbf{o}) + f(\mathbf{o},\mathbf{i}))$ over 1000 pairs of directions drawn uniformly on the upper hemisphere (default_rng(0)).
• H2, energy. The albedo $a(\theta_i) = \int f \cos\theta_o\, d\omega_o$ at $\theta_i = 0$, 30, 60, 75 and 85°, by adaptive double quadrature (scipy.integrate.dblquad, relative tolerance $10^{-10}$), using the symmetry about the plane of incidence.
• H3, Heitz’s weak white furnace test for GGX with Smith $G_1$, and with $G_1 = 1$, at $\alpha_g = 0.1$, 0.3, 0.6 and $\theta_o = 0$, 30, 60, 80 and 89°. The integral over $\mathbf{i}$ is changed to one over $\mathbf{h}$ with $d\omega_i = 4(\mathbf{o}\cdot\mathbf{h})\,d\omega_h$ (Walter et al., Eq. 14), over the upper hemisphere of $\mathbf{h}$ with $\mathbf{o}\cdot\mathbf{h} > 0$. This form of the integral is ours; Heitz’s MATLAB code, published under a no-derivatives licence, was not ported.
• H4, a sphere. An orthographic camera looks along $-z$ at a unit sphere, with view direction $\mathbf{v} = +z$, lit by a distant light of unit irradiance at angle $2\phi$ from the view in the $x$–$z$ plane, with $\phi = 5°$ and 22.5°. The radiance is $f\,\max(0, \mathbf{l}\cdot\mathbf{n})$ for P20, Pn20, BP80 (with $c = 82/(8\pi)$) and M0.1. On exact normals the script finds each highlight’s peak in the plane and its full width at half maximum (FWHM) in normal angle, in the plane and out of it through $\mathbf{h}$. It also renders a 513 × 513 image of the sphere and reports the angle between the brightest pixel’s normal and $\mathbf{h}$, and the highlight area, the number of pixels at or above half the image’s maximum.

Six checks verify the code: the Lambertian albedo equals $\rho$ (C1); the furnace routine without masking at $\theta_o = 0$ gives Walter et al.’s Eq. 4, $\int D(\mathbf{m})(\mathbf{n}\cdot\mathbf{m})\,d\omega_m = 1$ (C2); at normal incidence P’s albedo is $2\pi/(n+2)$ and Pn’s is 1 (C3); $F(0°) = 0.04$ and $F(90°) = 1$ (C4); P obeys $f(\mathbf{i},\mathbf{o})\cos\theta_i = f(\mathbf{o},\mathbf{i})\cos\theta_o$ (C5); and the random directions and the sphere’s normals are unit vectors in the upper hemisphere (C6).

The hypotheses were written before any script for this article ran. (H1) L, Pn, BP and M are reciprocal to within $10^{-12}$, and P is not (relative difference above 0.5). (H2a) P20’s albedo at 30° and 60° is within 5% of its normal-incidence value; (H2b) Pn’s albedo never exceeds 1 (to $10^{-9}$), does not increase with $\theta_i$, and for $n = 20$ lies within 5% of $\cos\theta_i$ at 30° and 60°; (H2c) at normal incidence BP80/2pi lies between 3.5 and 4.0 and BP80/8pi between 0.875 and 1.0; (H2d) M’s albedo stays below 1, and M0.1’s at normal incidence lies within 10% of $F(0) = 0.04$; (H2e) M F=1 never exceeds 1, at 85° falls with roughness, and for each roughness is lower at 85° than at normal incidence. (H3) With $G_1$ the furnace integral is 1 to within $10^{-8}$ in all fifteen cases; without it, it exceeds $1 + 10^{-6}$ at 60° and above and grows with both $\theta_o$ and $\alpha_g$. (H4a) P20’s peak is at $\mathbf{h}$ to within $10^{-4}$°, and the peaks of Pn20, BP80 and M0.1 lie on the light side of $\mathbf{h}$ by less than 1°; (H4b) the in-plane FWHM ratio Pn20/BP80 is within 0.02 of 1; (H4c) the out-of-plane ratio is within 0.01 of $1/\cos\phi$; (H4d) the highlight-area ratio is within 0.03 of $1/\cos\phi$. The numerical tolerances ($10^{-12}$ and 0.5 in H1, $10^{-9}$ for “never exceeds”, $10^{-8}$ and $10^{-6}$ in H3, $10^{-4}$° in H4a, and 0.02, 0.01 and 0.03 in H4b–H4d) were set in the plan, after a planning script had printed the reachable values of checks C1–C6 but no hypothesis measure. The script runs in about 12 s with NumPy 2.3.5 and SciPy 1.16.3 (Python 3.13).

import os
os.environ["OPENBLAS_NUM_THREADS"] = os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy import integrate, optimize

# Conventions: unit vectors; i points towards the light, o towards the viewer, both away from the
# surface; nrm is the surface normal. All BRDFs in sr^-1. Angles from the normal.
dot = lambda a, b: (a*b).sum(-1)

def f_L(i, o, nrm, rho=0.8):                       # Lambertian, rho/pi
    return np.full(np.broadcast_shapes(i.shape, o.shape, nrm.shape)[:-1], rho/np.pi)

def f_P(i, o, nrm, n, ks=1.0):                     # Phong as a BRDF: radiance ks (r.o)^n per unit E/cos
    r = 2*dot(i, nrm)[..., None]*nrm - i
    return ks*np.maximum(0.0, dot(r, o))**n/dot(i, nrm)

def f_Pn(i, o, nrm, n, ks=1.0):                    # normalized Phong
    r = 2*dot(i, nrm)[..., None]*nrm - i
    return ks*(n + 2)/(2*np.pi)*np.maximum(0.0, dot(r, o))**n

def f_BP(i, o, nrm, m, c):                         # Blinn-Phong, constant c given explicitly
    h = (i + o)/np.linalg.norm(i + o, axis=-1, keepdims=True)
    return c*np.maximum(0.0, dot(h, nrm))**m

def fresnel(c, eta=1.5):                           # Walter Eq. 22 (eta_t/eta_i = 1.5), c = |i.m|
    g = np.sqrt(eta*eta - 1 + c*c)
    return 0.5*(g - c)**2/(g + c)**2*(1 + (c*(g + c) - 1)**2/(c*(g - c) + 1)**2)

def D_ggx(cm, a):                                  # Walter Eq. 33 with chi+(m.n), cm = m.n
    with np.errstate(divide="ignore", invalid="ignore"):
        t2 = (1 - cm*cm)/(cm*cm)
        return np.where(cm > 0, a*a/(np.pi*cm**4*(a*a + t2)**2), 0.0)

def G1_ggx(cv, vm, a):                             # Walter Eq. 34, cv = v.n, vm = v.m
    t2 = np.maximum(0.0, 1 - cv*cv)/(cv*cv)
    return np.where(vm/cv > 0, 2.0/(1 + np.sqrt(1 + a*a*t2)), 0.0)

def f_M(i, o, nrm, a, F=True):                     # Walter Eq. 20 with Smith G = G1 G1 (Eq. 23)
    h = (i + o)/np.linalg.norm(i + o, axis=-1, keepdims=True)
    ci, co, ch, ih = dot(i, nrm), dot(o, nrm), dot(h, nrm), dot(i, h)
    Fr = fresnel(np.abs(ih)) if F else 1.0
    return Fr*G1_ggx(ci, ih, a)*G1_ggx(co, dot(o, h), a)*D_ggx(ch, a)/(4*np.abs(ci)*np.abs(co))

Z = np.array([0.0, 0.0, 1.0])
MODELS = {"L": lambda i, o, n: f_L(i, o, n),
          "P1": lambda i, o, n: f_P(i, o, n, 1), "P20": lambda i, o, n: f_P(i, o, n, 20),
          "Pn1": lambda i, o, n: f_Pn(i, o, n, 1), "Pn20": lambda i, o, n: f_Pn(i, o, n, 20),
          "BP80/2pi": lambda i, o, n: f_BP(i, o, n, 80, 82/(2*np.pi)),
          "BP80/8pi": lambda i, o, n: f_BP(i, o, n, 80, 82/(8*np.pi))}
for a in [0.1, 0.3, 0.6]:
    MODELS["M%.1f" % a] = (lambda a: lambda i, o, n: f_M(i, o, n, a))(a)
for a in [0.1, 0.3, 0.6]:
    MODELS["M%.1f F=1" % a] = (lambda a: lambda i, o, n: f_M(i, o, n, a, F=False))(a)

EPS = dict(epsabs=0.0, epsrel=1e-10)
def albedo(f, ti):                                 # a(theta_i) = int f cos(theta_o) d omega_o (phi-symmetric)
    i = np.array([np.sin(ti), 0.0, np.cos(ti)])
    def g(ph, th):
        o = np.array([np.sin(th)*np.cos(ph), np.sin(th)*np.sin(ph), np.cos(th)])
        return float(f(i, o, Z))*np.cos(th)*np.sin(th)
    return 2*integrate.dblquad(g, 0, np.pi/2, 0, np.pi, **EPS)[0]

def furnace(a, to, masking=True):                  # Heitz Eq. 36, mapped to h: d omega_i = 4 (o.h) d omega_h
    o = np.array([np.sin(to), 0.0, np.cos(to)])
    def g(ph, th):
        h = np.array([np.sin(th)*np.cos(ph), np.sin(th)*np.sin(ph), np.cos(th)])
        oh = o @ h
        G = G1_ggx(o[2], oh, a) if masking else 1.0
        return float(G*D_ggx(np.cos(th), a)*oh/o[2])*np.sin(th)
    def phmax(th):                                 # keep o.h > 0
        if to == 0.0 or th == 0.0:
            return np.pi
        return np.arccos(np.clip(-1/(np.tan(to)*np.tan(th)), -1.0, 1.0))
    return 2*integrate.dblquad(g, 0, np.pi/2, 0, phmax, **EPS)[0]

print("== Checks")
TI = [0, 30, 60, 75, 85]
alb = {k: [albedo(f, np.radians(t)) for t in TI] for k, f in MODELS.items()}
print("C1 Lambertian albedo - rho (0.8), max over theta_i: %.1e (tol 1e-12)" % np.abs(np.array(alb["L"]) - 0.8).max())
print("C2 furnace without G1 at theta_o = 0 (Walter Eq. 4) - 1, alpha 0.1/0.3/0.6:",
      " ".join("%.1e" % (furnace(a, 0.0, masking=False) - 1) for a in [0.1, 0.3, 0.6]), "(tol 1e-12)")
print("C3 albedo(0): P1 - 2pi/3 %.1e, P20 - 2pi/22 %.1e, Pn1 - 1 %.1e, Pn20 - 1 %.1e (tol 1e-12)" % (
      alb["P1"][0] - 2*np.pi/3, alb["P20"][0] - 2*np.pi/22, alb["Pn1"][0] - 1, alb["Pn20"][0] - 1))
print("C4 Fresnel eta 1.5: F(0 deg) - 0.04 = %.1e, F(90 deg) - 1 = %.1e (tol 1e-15)" % (fresnel(1.0) - 0.04, fresnel(0.0) - 1))
rng = np.random.default_rng(0)
def hemi(k):                                       # uniform on the upper hemisphere
    z = rng.random(k); p = 2*np.pi*rng.random(k); s = np.sqrt(1 - z*z)
    return np.stack([s*np.cos(p), s*np.sin(p), z], -1)
I, O = hemi(1000), hemi(1000)
for n in [1, 20]:
    a, b = f_P(I, O, Z, n)*I[:, 2], f_P(O, I, Z, n)*O[:, 2]
    m = (a + b) > 0
    print("C5 P%d: max |f(i,o)cos_i - f(o,i)cos_o|/sum = %.1e over %d nonzero pairs (tol 1e-12)" % (
          n, (np.abs(a - b)[m]/(a + b)[m]).max(), m.sum()))
print("C6 pairs: %d and %d of 1000 with z > 0, max |norm - 1| %.1e (tol 1e-15)" % ((I[:, 2] > 0).sum(), (O[:, 2] > 0).sum(),
      max(np.abs(np.linalg.norm(I, axis=1) - 1).max(), np.abs(np.linalg.norm(O, axis=1) - 1).max())))

print("\n== H1 reciprocity: max |f(i,o) - f(o,i)| / (f(i,o) + f(o,i)) over 1000 random pairs")
for k, f in MODELS.items():
    a, b = f(I, O, Z), f(O, I, Z)
    m = (a + b) > 0
    print("%-11s %.1e  (%d nonzero pairs)" % (k, (np.abs(a - b)[m]/(a + b)[m]).max(), m.sum()))

print("\n== H2 directional albedo a(theta_i)")
print("theta_i " + " ".join("%10s" % k for k in list(MODELS)[:7]))
for j, t in enumerate(TI):
    print("%5d   " % t + " ".join("%10.6f" % alb[k][j] for k in list(MODELS)[:7]))
print("theta_i " + " ".join("%10s" % k for k in list(MODELS)[7:]))
for j, t in enumerate(TI):
    print("%5d   " % t + " ".join("%10.6f" % alb[k][j] for k in list(MODELS)[7:]))
print("cos(theta_i):", " ".join("%.6f" % np.cos(np.radians(t)) for t in TI))

print("\n== H3 weak white furnace test (Heitz Eq. 36), GGX + Smith G1, and with G1 = 1")
print("alpha theta_o   with G1         W-1      without G1")
for a in [0.1, 0.3, 0.6]:
    for to in [0, 30, 60, 80, 89]:
        w1, w0 = furnace(a, np.radians(to)), furnace(a, np.radians(to), masking=False)
        print("%.1f   %5d    %.10f  %8.1e   %.10f" % (a, to, w1, w1 - 1, w0))

print("\n== H4 sphere seen along +z, light at 2*phi from the view in the x-z plane")
V = Z
H4M = {"P20": MODELS["P20"], "Pn20": MODELS["Pn20"], "BP80": MODELS["BP80/8pi"], "M0.1": MODELS["M0.1"]}
def radiance(f, l, nrm):                           # unit irradiance at normal incidence; 0 where l.n <= 0
    ln = dot(l, nrm)
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(ln > 0, f(l, V, nrm)*ln, 0.0)
def nrm_in(t):                                     # in-plane normal at angle t from +z towards +x
    return np.array([np.sin(t), 0.0, np.cos(t)])
def nrm_out(phi, b):                               # tilted by b out of the plane, through h
    return np.array([np.sin(phi)*np.cos(b), np.sin(b), np.cos(phi)*np.cos(b)])
def fwhm(prof, t0, span):                          # full width at half maximum around a peak at t0
    half = prof(t0)/2
    lo = optimize.brentq(lambda t: prof(t) - half, t0 - span, t0, xtol=1e-13)
    hi = optimize.brentq(lambda t: prof(t) - half, t0, t0 + span, xtol=1e-13)
    return hi - lo
N = 513; c = N//2
yy, xx = np.mgrid[0:N, 0:N]
X, Y = (xx - c)/c, (c - yy)/c
inside = X**2 + Y**2 < 1
NR = np.stack([X, Y, np.sqrt(np.clip(1 - X**2 - Y**2, 0, None))], -1)
print("C6 sphere image %dx%d: %d pixels inside, max |norm - 1| %.1e, min n_z %.4f" % (N, N, inside.sum(),
      np.abs(np.linalg.norm(NR[inside], axis=1) - 1).max(), NR[inside][:, 2].min()))
for phi_d in [5.0, 22.5]:
    phi = np.radians(phi_d)
    l = np.array([np.sin(2*phi), 0.0, np.cos(2*phi)])
    print("phi = %.1f deg (light-view angle %.0f deg); 1/cos(phi) = %.4f" % (phi_d, 2*phi_d, 1/np.cos(phi)))
    print("  model  peak(deg)   peak-phi(deg)  FWHM_in(deg)  FWHM_out(deg)  img peak dev(deg)  area(px)")
    res = {}
    for k, f in H4M.items():
        pin = lambda t: float(radiance(f, l, nrm_in(t)))
        pk = optimize.minimize_scalar(lambda t: -pin(t), bounds=(phi - np.radians(20), phi + np.radians(20)),
                                      method="bounded", options=dict(xatol=1e-10)).x
        w_in = fwhm(pin, pk, np.radians(30))
        w_out = fwhm(lambda b: float(radiance(f, l, nrm_out(phi, b))), 0.0, np.radians(30))
        img = np.where(inside, radiance(f, l, NR), 0.0)
        p = np.unravel_index(np.argmax(img), img.shape)
        hvec = nrm_in(phi)
        dev = np.degrees(np.arccos(np.clip(NR[p] @ hvec, -1, 1)))
        area = int((img >= 0.5*img.max()).sum())
        res[k] = (w_in, w_out, area)
        print("  %-5s  %.6f   %+.6f      %.4f       %.4f        %.3f              %d" % (
              k, np.degrees(pk), np.degrees(pk) - phi_d, np.degrees(w_in), np.degrees(w_out), dev, area))
    print("  Pn20/BP80 ratios: FWHM_in %.4f, FWHM_out %.4f, area %.4f" % (
          res["Pn20"][0]/res["BP80"][0], res["Pn20"][1]/res["BP80"][1], res["Pn20"][2]/res["BP80"][2]))

== Checks
C1 Lambertian albedo - rho (0.8), max over theta_i: 1.1e-16 (tol 1e-12)
C2 furnace without G1 at theta_o = 0 (Walter Eq. 4) - 1, alpha 0.1/0.3/0.6: 4.4e-16 0.0e+00 0.0e+00 (tol 1e-12)
C3 albedo(0): P1 - 2pi/3 0.0e+00, P20 - 2pi/22 -5.6e-17, Pn1 - 1 -1.1e-16, Pn20 - 1 -1.1e-16 (tol 1e-12)
C4 Fresnel eta 1.5: F(0 deg) - 0.04 = 0.0e+00, F(90 deg) - 1 = 0.0e+00 (tol 1e-15)
C5 P1: max |f(i,o)cos_i - f(o,i)cos_o|/sum = 1.3e-16 over 694 nonzero pairs (tol 1e-12)
C5 P20: max |f(i,o)cos_i - f(o,i)cos_o|/sum = 1.4e-16 over 694 nonzero pairs (tol 1e-12)
C6 pairs: 1000 and 1000 of 1000 with z > 0, max |norm - 1| 2.2e-16 (tol 1e-15)

== H1 reciprocity: max |f(i,o) - f(o,i)| / (f(i,o) + f(o,i)) over 1000 random pairs
L           0.0e+00  (1000 nonzero pairs)
P1          1.0e+00  (694 nonzero pairs)
P20         1.0e+00  (694 nonzero pairs)
Pn1         0.0e+00  (694 nonzero pairs)
Pn20        0.0e+00  (694 nonzero pairs)
BP80/2pi    0.0e+00  (1000 nonzero pairs)
BP80/8pi    0.0e+00  (1000 nonzero pairs)
M0.1        2.9e-15  (1000 nonzero pairs)
M0.3        2.7e-15  (1000 nonzero pairs)
M0.6        2.8e-15  (1000 nonzero pairs)
M0.1 F=1    0.0e+00  (1000 nonzero pairs)
M0.3 F=1    0.0e+00  (1000 nonzero pairs)
M0.6 F=1    0.0e+00  (1000 nonzero pairs)

== H2 directional albedo a(theta_i)
theta_i          L         P1        P20        Pn1       Pn20   BP80/2pi   BP80/8pi
    0     0.800000   2.094395   0.285599   1.000000   1.000000   3.809524   0.952381
   30     0.800000   2.130229   0.285599   0.880843   0.866025   2.880952   0.720238
   60     0.800000   2.550964   0.285890   0.608998   0.500509   1.024374   0.256093
   75     0.800000   3.709764   0.299238   0.458442   0.271179   0.346067   0.086517
   85     0.800000   8.725410   0.453171   0.363097   0.138293   0.107752   0.026938
theta_i       M0.1       M0.3       M0.6   M0.1 F=1   M0.3 F=1   M0.6 F=1
    0     0.039637   0.035597   0.024310   0.988304   0.877358   0.591523
   30     0.041663   0.037595   0.025875   0.986085   0.863622   0.594009
   60     0.085870   0.060575   0.036731   0.969115   0.818133   0.624722
   75     0.195510   0.098615   0.052830   0.923240   0.808411   0.676011
   85     0.314500   0.138661   0.070984   0.873130   0.847002   0.735475
cos(theta_i): 1.000000 0.866025 0.500000 0.258819 0.087156

== H3 weak white furnace test (Heitz Eq. 36), GGX + Smith G1, and with G1 = 1
alpha theta_o   with G1         W-1      without G1
0.1       0    1.0000000000   4.4e-16   1.0000000000
0.1      30    1.0000000000   4.4e-11   1.0008326401
0.1      60    1.0000000001   6.6e-11   1.0074445783
0.1      80    1.0000000000  -1.5e-13   1.0748117898
0.1      89    1.0000000000   5.1e-14   3.4078083257
0.3       0    1.0000000000   0.0e+00   1.0000000000
0.3      30    1.0000000000   1.2e-11   1.0074445783
0.3      60    1.0000000000   1.2e-11   1.0634713835
0.3      80    1.0000000000  -1.1e-13   1.4867509023
0.3      89    1.0000000000   2.2e-13   9.1080278422
0.6       0    1.0000000000   0.0e+00   1.0000000000
0.6      30    1.0000000000   3.9e-11   1.0291502623
0.6      60    1.0000000000   1.5e-11   1.2211102551
0.6      80    1.0000000000  -5.6e-13   2.2733328433
0.6      89    1.0000000000   2.3e-12   17.6942598948

== H4 sphere seen along +z, light at 2*phi from the view in the x-z plane
C6 sphere image 513x513: 205857 pixels inside, max |norm - 1| 1.1e-16, min n_z 0.0068
phi = 5.0 deg (light-view angle 10 deg); 1/cos(phi) = 1.0038
  model  peak(deg)   peak-phi(deg)  FWHM_in(deg)  FWHM_out(deg)  img peak dev(deg)  area(px)
  P20    5.000000   +0.000000      14.9977       15.0553        0.070              3513
  Pn20   5.061880   +0.061880      14.9059       14.9632        0.155              3470
  BP80   5.061880   +0.061880      14.9692       14.9699        0.155              3480
  M0.1   5.012690   +0.012690      7.4308       7.4306        0.070              863
  Pn20/BP80 ratios: FWHM_in 0.9958, FWHM_out 0.9996, area 0.9971
phi = 22.5 deg (light-view angle 45 deg); 1/cos(phi) = 1.0824
  model  peak(deg)   peak-phi(deg)  FWHM_in(deg)  FWHM_out(deg)  img peak dev(deg)  area(px)
  P20    22.500001   +0.000001      14.9977       16.2414        0.008              3503
  Pn20   22.792376   +0.292376      14.8907       16.1258        0.250              3449
  BP80   22.792383   +0.292383      14.9542       14.9699        0.250              3214
  M0.1   22.560113   +0.060113      7.4346       7.4306        0.008              801
  Pn20/BP80 ratios: FWHM_in 0.9958, FWHM_out 1.0772, area 1.0731

All six checks held. The Lambertian albedo is 0.8 at every angle, GGX’s projected area is 1 (C2), and the normal-incidence albedos of P and Pn match their closed forms to rounding (C3), as do the Fresnel values (C4), Phong’s cosine-ratio identity (C5) and the unit vectors (C6).

Reciprocity. H1 held. Every model except Phong is reciprocal to rounding (at most $2.9 \times 10^{-15}$, for M with Fresnel). Phong’s relative difference reaches 1.0 because its two values differ by the ratio $\cos\theta_o/\cos\theta_i$, which C5 confirms, and that ratio is far from 1 for pairs with one direction near grazing. Only 694 of the 1000 pairs give a non-zero Phong or normalized-Phong value, because the lobe is cut off where $\mathbf{r}\cdot\mathbf{o} \le 0$.

Energy. H2a held: P20’s albedo is 0.285599 at 0° and 30° and 0.285890 at 60°. H2b held: Pn1 and Pn20 fall from 1 with $\theta_i$ and never exceed it, and Pn20’s albedo is 0.866025 at 30° and 0.500509 at 60°, where $\cos\theta_i$ is 0.866025 and 0.5. Both results follow from one identity (our derivation): over the whole hemisphere around $\mathbf{r}$, $\int (\mathbf{r}\cdot\mathbf{o})^n\, \mathbf{o}\, d\omega = 2\pi\mathbf{r}/(n+2)$, so the lobe’s integral against $\cos\theta_o$ would be $2\pi\cos\theta_i/(n+2)$ if the horizon did not cut it. The horizon removes the part of the lobe below the surface, where $\cos\theta_o < 0$, so the true integral is never smaller. Phong divides by $\cos\theta_i$ and so has an albedo of at least $2\pi k_s/(n+2)$ at every angle, which grows near grazing as more of the lobe is cut off: P1 reaches 8.73 at 85° and P20 0.453. Normalized Phong does not divide, so its albedo is at least $k_s\cos\theta_i$, and every printed Pn value lies at or above $\cos\theta_i$. The unnormalized-by-cosine Phong model therefore reflects more light than it receives for broad lobes (P1, 2.09 at normal incidence) and near grazing, while for $n = 20$ it reflects at most 45% of it at the printed angles.

H2c held. The Blinn–Phong lobe normalized as a distribution of normals, BP80/2pi, has an albedo of 3.81 at normal incidence: it reflects almost four times the light it receives, because the outgoing direction sweeps four times the solid angle of the half-vector (our derivation from Walter et al.’s Eq. 14). Dividing by 4 gives 0.952 at normal incidence; the two columns differ by exactly 4 by construction. With the first constant the model reflects more than it receives up to 60° (1.024 there); with the second it stays below 1 at every angle but falls well below $\cos\theta_i$ away from normal incidence (0.256 at 60°), unlike normalized Phong.

H2d held. The microfacet model with Fresnel reflects 0.0396 at normal incidence for $\alpha_g = 0.1$, close to the 0.04 of a single smooth facet of glass, and less for rougher surfaces (0.0356 and 0.0243). Its albedo rises towards grazing, to 0.3145 at 85° for $\alpha_g = 0.1$; in our reading this is the Fresnel term, which Walter et al. note rises to 1 at grazing angles. H2e failed in one of its three parts. With $F = 1$, no albedo exceeds 1 (the largest is 0.988304), and at 85° it falls with roughness (0.873130, 0.847002, 0.735475), but it is not lower at 85° than at normal incidence for every roughness: it is for $\alpha_g = 0.1$ (0.873130 against 0.988304) and 0.3 (0.847002 against 0.877358), and it failed for $\alpha_g = 0.6$, where the albedo is 0.735475 at 85° against 0.591523 at normal incidence and rises at every printed angle. That all values with $F = 1$ lie below 1 follows for every printed case (our derivation). By Heitz’s Eq. 14 for Smith’s $G_1$, the factor $G_1(\mathbf{i}, \mathbf{h})\,D(\mathbf{h})/(4\cos\theta_i)$ alone integrates to 1 over outgoing directions, which is the weak furnace integral with $\mathbf{i}$ and $\mathbf{o}$ exchanged; $1 - a(\theta_i)$ is then the part removed by the remaining factor $G_1(\mathbf{o}, \mathbf{h}) \le 1$, which is zero below the horizon. Read by reciprocity as Heitz’s White Furnace integral with the viewer at $\theta_i$, this is the energy his shadowing term discards. At normal incidence the loss grows with roughness, from 1.2% to 41%.

The furnace test. H3 held. With Smith masking the weak white furnace integral is 1 in all fifteen cases, to within $6.6 \times 10^{-11}$, so Walter et al.’s $G_1$ for GGX passes Heitz’s normalization test at every roughness and view angle printed. Without masking the integral exceeds 1 at every printed angle away from the normal, by 0.08% at 30° for $\alpha_g = 0.1$ and by a factor of 17.7 at 89° for $\alpha_g = 0.6$: the facets turned away from the viewer are counted as visible. The values match $1 + \Lambda = 1/G_1(\mathbf{o})$ from the formal section (our derivation) in the five cases we checked by hand: 1.0008326 ($\alpha_g = 0.1$, 30°), 1.0074446 (0.1, 60°), 1.0074446 (0.3, 30°), 1.0291503 (0.6, 30°) and 1.2211103 (0.6, 60°). The cases $\alpha_g = 0.1$ at 60° and 0.3 at 30° print the same value, 1.0074445783, because $\alpha_g\tan\theta_o$ is the same, 0.1732, for both.

Highlights on a sphere. H4 held in all four parts. P20’s peak is at $\mathbf{h}$ (deviation 0.000000° and 0.000001°): its radiance $k_s(\mathbf{r}\cdot\mathbf{v})^n$ peaks exactly where the mirror direction meets the viewer. Normalized Phong and Blinn–Phong carry the extra factor $\mathbf{n}\cdot\mathbf{l}$ in their radiance, which pulls the peak towards the light, by 0.062° at a light–view angle of 10° and 0.292° at 45°. To first order the shift is $\tan\phi/(4n)$, 0.063° and 0.297° (our derivation and arithmetic). The microfacet peak moves less, 0.013° and 0.060°; in our reading because its lobe is about half as wide (FWHM 7.43° against 14.9°). In the plane of incidence normalized Phong and Blinn–Phong with $m = 4n$ have the same width to within 0.5% (ratio 0.9958 at both angles). Out of the plane they agree at 10° (0.9996, against $1/\cos\phi = 1.0038$) but not at 45°, where the Phong lobe is 7.7% wider (1.0772 against 1.0824), and the rendered highlight is 7.3% larger in area (1.0731). P20 shows the same widening: its out-of-plane width grows from 15.06° to 16.24° while its in-plane width stays at 15.00°. In the 513 × 513 image the brightest pixel’s normal lies within 0.25° of $\mathbf{h}$ for every model. One pixel is 1/256 of the radius, a step of about 0.22° to 0.24° in normal angle near these peaks (our arithmetic), so the image is too coarse to measure the shifts that the continuous profiles find.

Process note. The plan for this example listed nine departures from the batch plan’s description, before the pilot: (1) Phong as a BRDF is fixed as $k_s(\mathbf{r}\cdot\mathbf{o})^n/(\mathbf{n}\cdot\mathbf{i})$, our conversion, and run at $n = 1$ and 20; (2) Blinn–Phong is printed with two constants, since no read source gives its normalization, which makes six model families rather than five; (3) M is also run with $F = 1$, to show the energy that shadowing removes; (4) the weak furnace integral is computed in half-vector coordinates through Walter et al.’s Jacobian, not on the sphere of $\mathbf{i}$ as in Heitz’s code, which is not ported; (5) the width comparison is normalized Phong against Blinn–Phong, since both carry $\cos\theta_i$ in their radiance, with Phong printed beside them, widths measured on continuous profiles of exact normals, the image giving the peak pixel and highlight area, at two light–view angles; (6) Walter et al.’s Beckmann–Phong relation $\alpha_p = 2\alpha_b^{-2} - 2$ is stated from the source and not computed; (7) checks C1–C6 were added; (8) the planning research’s expectation that Phong’s albedo exceeds 1 for $k_s = 1$ was replaced by the exact normal-incidence value $2\pi k_s/(n+2)$, so $n = 1$ runs beside $n = 20$; (9) the research’s expectations for normalized Phong and the microfacet model near grazing were restated as H2b and H2e. One edit was made to the script before it first ran: the indicator $\chi^+(\mathbf{m}\cdot\mathbf{n})$ of Walter et al.’s Eq. 33 was added to the GGX distribution, and the image radiance was set to zero where $\mathbf{l}\cdot\mathbf{n} \le 0$, since the sphere has normals facing away from the light. The pilot’s output was identical to the first run’s, and the first run gave identical output when run twice; no setting, measure or printed value changed after the pilot.

Common Misconceptions

• “The BRDF is the fraction of light reflected.” It is a density per steradian and can be much larger than 1 near a mirror direction. The fraction is the albedo, its integral over projected solid angle.
• “The Lambertian BRDF is the albedo.” It is the albedo divided by $\pi$ (Nicodemus et al., Eq. C4). Formulas written for reflected light rather than radiance, such as visionbook’s Eq. 5.2, leave the $\pi$ out.
• “Phong is a BRDF.” Phong’s model is a formula for reflected light. Written as a BRDF it is not reciprocal, and its energy depends on its exponent; the example’s P1 reflects more than twice the light it receives at normal incidence.
• “A Blinn–Phong exponent is a Phong exponent.” For the same width in the plane of incidence the Blinn–Phong exponent is about four times the Phong exponent, and the shapes still differ out of that plane when the light and the camera are apart.
• “Microfacet BRDFs conserve energy, so they reflect all the light when Fresnel is 1.” They model one bounce. The shadowing term removes the rest, up to 41% of the light at normal incidence for $\alpha_g = 0.6$ in the example; Heitz’s weak test checks masking alone for this reason.
• “GGX’s α is the roughness.” It is the width parameter of Walter et al.’s Eq. 33. Beckmann’s $\alpha_b$ and the Phong distribution’s exponent $\alpha_p$ are different parameters, and a renderer’s roughness setting may be mapped to $\alpha_g$ in its own way.
• “Specular highlights are where the normal bisects light and view.” Among the example’s models, that holds exactly only for Phong’s radiance; models with a $\cos\theta_i$ factor or a $1/\cos\theta_o$ factor move the peak slightly, by tenths of a degree in the example.

Where It Is Used

Image formation puts the BRDF in its reflection equation, and radiometry and photometry supplies the units. Shape from shading and photometric stereo invert the Lambertian model to recover normals; photometric stereo also surveys methods that fit Torrance–Sparrow, Ward or other non-Lambertian models, or rely only on general BRDF properties such as isotropy, and learned methods trained on images rendered with measured or principled BRDFs. Shape from shading lists non-Lambertian reflectance among the effects that change its image irradiance equation. Brightness constancy, the assumption behind the optical flow constraint equation, fails at specular highlights; in our reading it holds for Lambertian surfaces under unchanged light, whose radiance does not depend on the view. Differentiable rendering differentiates images with respect to material parameters, and its early rasterizers assumed Lambertian materials. Neural radiance fields map a position and a viewing direction to colour, which in our reading models view-dependent reflection without an explicit BRDF. The lidar range equation uses a Lambertian target’s $\Gamma/\pi$, and multispectral satellite imaging the same $\pi$ when it converts radiance into reflectance.