Define Oren-Nayar lobe explicitly as EON model (i.e. energy-compensat…

…ed Fujii) (#174)

As discussed extensively on Slack, our specification of the diffuse lobe is incomplete since we don't make it clear what form of the Oren-Nayar model to take, and there are several options.

As argued on Slack, I think the model proposed by Fujii makes the most sense, as it: 
   - fixes the dark ring artifacts in the original "qualitative Oren-Nayar"
   - is less dark than the original "qualitative Oren-Nayar" (energy preserving at grazing)
   - $\sigma$ parameter has more obvious meaning as the $[0,1]$ "weight of non-Lambertian-ness".
   - looks better in renders (subjectively)

It is also easy to add straightforward energy compensation so that:
   - the model is fully energy preserving, so passes a white furnace test
   - exhibits saturation of the albedo, as physically expected for a multiple scattering microfacet model

This model is also very easy to implement (no tabulation) and efficient (even more-so than qualitative Oren-Nayar), via an accurate fit thanks to the analytical solution for the albedo.

index.html

@@ -570,7 +570,8 @@
 
 ### Glossy-diffuse
 
-The glossy-diffuse slab represents a dielectric with rough GGX [microfacet](index.html#model/microfacetmodel) surface BSDF (whose "specular" parameters are described in the Dielectric base section), embedding a semi-infinite bulk of extremely dense scattering material.
+The glossy-diffuse slab represents the base dielectric (parametrized as described in the Dielectric base section), embedding a semi-infinite bulk of extremely dense scattering material. The BRDF of the slab is the combination of a "glossy" specular lobe provided by immediate reflection from the dielectric interface, and a diffuse lobe provided by scattering off the embedded substrate. This models for example the reflection from shiny, totally opaque surfaces such as dense plastic, rock, and concrete.
+
 We choose to model this concretely as a layer of dielectric "gloss" on top of an opaque slab with a diffuse BRDF:
 \begin{eqnarray}
 M_\textrm{glossy-diffuse} = \mathrm{\mathbf{layer}}(S_\mathrm{diffuse}, S_\textrm{gloss})
@@ -592,32 +593,61 @@
 \begin{eqnarray}
 S_\textrm{diffuse}    = \mathrm{Slab}(f_\mathrm{diffuse})
 \end{eqnarray}
-where $f_\mathrm{diffuse}$ is defined to be the BRDF of the classic Oren-Nayar microfacet model ([#Oren1994], [#Pharr2023]) in which the microfacets are individually Lambertian. The resulting BRDF of the glossy-diffuse slab is the combination of a "glossy" specular lobe provided by reflection from the dielectric interface, and a diffuse lobe provided by bounces off the substrate. This models for example the reflection from shiny, totally opaque surfaces such as dense plastic, rock, and concrete.
+where $f_\mathrm{diffuse}$ is based on the classic Oren-Nayar model ([#Oren1994], [#Pharr2023]), which is a V-cavity based microfacet model in which the microfacets are assumed to be individually Lambertian. The Oren-Nayar model is parametrized by roughness $\sigma$ and an overall RGB scale factor $\boldsymbol{\rho}$. The roughness models the "flattening" effect (Figure [rough_diffuse]), which is observed in real rough diffuse materials [#Oren1994].
+
+![](images/diffuse_smooth.png width=98% align=center) ![](images/diffuse_rough.png width=98% align=center)
+<div class="shifted-caption">
+    ![Figure [rough_diffuse]: Zero roughness ($\sigma=0$, left) and maximum roughness ($\sigma=1$, right) diffuse materials with $\textbf{EON}$ BRDF.](dummy)
+</div>
+
+Unfortunately, the original model suffers from artifacts and also does not conserve energy (i.e. is too dark). We thus opt to define the diffuse BRDF to be a popular improved version of Oren-Nayar introduced by [#Fujii2012], augmented with a simple analytical, reciprocal energy compensation term:
+\begin{eqnarray} \label{EON_brdf}
+f_\mathrm{diffuse}(\omega_i, \omega_o) = f_\mathrm{ON}(\omega_i, \omega_o) + f^\mathrm{comp}_\mathrm{ON}(\omega_i, \omega_o) \ .
+\end{eqnarray}
+This form of the Oren-Nayar model is termed "energy-preserving Oren-Nayar" or $\textbf{EON}$.
+The Oren-Nayar term $f_\mathrm{ON}$ is given by the [#Fujii2012] form [^Oren_Nayar_formula]
+\begin{eqnarray} \label{FON_brdf}
+f_\mathrm{ON}(\omega_i, \omega_o) = \frac{\boldsymbol{\rho}}{\pi} \Bigl( A(\sigma) + B(\sigma) \frac{s}{t} \Bigr) \ .
+\end{eqnarray}
+The roughness parameter $\sigma \in [0,1]$ is given by **`base_diffuse_roughness`**. The RGB $\boldsymbol{\rho}$ parameter is determined by the specified **`base_color`**, as described below. The directional albedo $E_\mathrm{ON}(\omega) = \boldsymbol{\rho}\,\hat{E}_\mathrm{ON}(\omega)$ (and corresponding _average albedo_ $\langle\hat{E}_\mathrm{ON}\rangle$) of the Oren-Nayar term can be determined analytically [^Oren_Nayar_albedo].
 
-Physically there will be darkening and saturation of the observed color due to multiple inter-reflections (scattering) within the gloss layer, as described in the coat darkening section. However in this case, we wish to automatically compensate for this effect by adjusting the albedo of $f_\mathrm{diffuse}$ in order to make the observed color match the input color $\mathbf{C} =$ **`base_weight`** * **`base_color`** (in a similar spirit to the albedo remapping described in the Subsurface section).
+The energy compensation term $f^{\mathrm{comp}}_\mathrm{ON}$ is given in terms of the albedo $\langle\hat{E}_\mathrm{ON}\rangle$ by
+\begin{equation} \label{EON_comp}
+f^\mathrm{comp}_\mathrm{ON}(\omega_i, \omega_o) = \frac{\boldsymbol{\rho}_\mathrm{ms}}{\pi}
+\bigl(1 - \hat{E}_\mathrm{ON}(\omega_i)\bigr)
+\bigl(1 - \hat{E}_\mathrm{ON}(\omega_o)\bigr) \ ,
+\end{equation}
+where the factor $\boldsymbol{\rho}_\mathrm{ms}$ accounts approximately for multiple scattering on the microfacet surface:
+\begin{equation}
+\boldsymbol{\rho}_\mathrm{ms} = \frac{\boldsymbol{\rho}^2}{\pi}
+\frac{\langle\hat{E}_\mathrm{ON}\rangle / (1 - \langle \hat{E}_\mathrm{ON}\rangle)}{1 - \boldsymbol{\rho} \bigl(1 - \langle \hat{E}_\mathrm{ON}\rangle\bigr)} \ .
+\end{equation}
+One can verify that as $\boldsymbol{\rho} \rightarrow 1$, the total directional albedo of the BRDF of Equation [EON_brdf], $E_\mathrm{diffuse}(\omega) \rightarrow 1$, thus the compensation term ensures that the white furnace test passes. Note that in the the zero roughness ($\sigma \rightarrow 0$) limit, the energy compensation term vanishes, and the parameter $\boldsymbol{\rho}$ is equal to the albedo of $f_\mathrm{diffuse}$. As roughness increases, the albedo of $f_\mathrm{diffuse}$ becomes slightly more dark and saturated than $\boldsymbol{\rho}$ due to the multiple scattering, which is physically realistic.
 
+Physically there will be additional darkening and saturation of the observed color due to multiple inter-reflections (scattering) within the gloss layer, as described in the coat darkening section. However in this case, we wish to automatically compensate for this effect by adjusting the $\boldsymbol{\rho}$ parameter of $f_\mathrm{diffuse}$ in order to make the observed color match the input color $\mathbf{C} =$ **`base_weight`** * **`base_color`** (in a similar spirit to the albedo remapping described in the Subsurface section).
 To define the meaning of the specified color in terms of the underlying base albedo, we write the normal-direction reflectance of the glossy-diffuse slab, $\mathbf{E}_\textrm{glossy-diffuse}$ in the general form
 \begin{eqnarray}
 \mathbf{E}_\textrm{glossy-diffuse} = \mathbf{E}_\textrm{spec} + \mathbf{E}_\textrm{diffuse}
 \end{eqnarray}
-where $\mathbf{E}_\mathrm{spec}$ is the normal-direction reflectance of all energy reflected from the dielectric interface without transmission, and $\mathbf{E}_\mathrm{diffuse}$ is the normal-direction reflectance of all energy transmitted through the interface, scattered off the diffuse medium, and transmitted back out.
-We then _define_ $\mathbf{C} =$ **`base_weight`** * **`base_color`** to be such that the reflectance of the remaining energy transmitted into the slab, $\mathbf{E}_\textrm{diffuse}$ in the *zero* **`base_diffuse_roughness`** case (i.e. a Lambertian base), is given by
+where $\mathbf{E}_\mathrm{spec}$ is the normal-direction reflectance of all energy reflected from the dielectric interface _without_ macroscopic transmission, while $\mathbf{E}_\mathrm{diffuse}$ is the normal-direction reflectance of all remaining energy due to (macroscopic) transmission through the interface, multiple scattering in the diffuse medium, and transmission back out.
+
+We then _define_ $\mathbf{C} =$ **`base_weight`** * **`base_color`** to be such that the reflectance of the remaining energy transmitted into the slab, $\mathbf{E}_\textrm{diffuse}$ in the Lambertian case of  *zero* **`base_diffuse_roughness`** (i.e. $\sigma = 0$), is given by
 \begin{eqnarray} \label{glossy_diffuse_albedo_constraint}
 \mathbf{E}_\textrm{diffuse} = \bigl( 1 - \mathbf{E}_\textrm{spec} \bigr) \mathbf{C} \ .
 \end{eqnarray}
 In other words, the selected color $\mathbf{C}$ parametrizes the fraction of the energy transmitted into the dielectric layer that is subsequently transmitted back out due to reflection from a Lambertian interface. Thus since $0 \le \mathbf{C} \le 1$, $\mathbf{E}_\textrm{glossy-diffuse} \le 1$ so energy is always conserved, and if $\mathbf{C}=1$ then $\mathbf{E}_\textrm{glossy-diffuse}=1$ which guarantees that a white $\mathbf{C}$ passes the furnace test. According to this definition, $\mathbf{C}$ is the observed reflection color (viewed at normal incidence under uniform illumination) in areas where the Fresnel reflection is negligible, and otherwise the observed color is a blend of $\mathbf{C}$ with the gray Fresnel reflection that conserves total reflected energy.
 
-Given the required diffuse albedo $\mathbf{E}_\textrm{diffuse}$ according to the formula above, then in principle, the albedo of the diffuse Oren-Nayar lobe $f_\mathrm{diffuse}(\omega_i, \omega_o)$ which generates the required $\mathbf{E}_\textrm{diffuse}$ can be determined. A reasonable practical, albeit non-reciprocal, approximation to the resulting BRDF of the glossy-diffuse slab, which satisfies this requirement is obtained via the non-reciprocal albedo-scaling approximation of equation [non-reciprocal-albedo-scaling]:
+Given the required diffuse albedo $\mathbf{E}_\textrm{diffuse}$ according to the formula above, then in principle, the $\boldsymbol{\rho}$ parameter of the Oren-Nayar model $f_\mathrm{diffuse}(\omega_i, \omega_o)$ which generates the required $\mathbf{E}_\textrm{diffuse}$ can be determined. A reasonable practical, albeit non-reciprocal, approximation to the resulting BRDF of the glossy-diffuse slab, which satisfies this requirement is obtained via the non-reciprocal albedo-scaling approximation of equation [non-reciprocal-albedo-scaling]:
 \begin{equation}
-f_\textrm{glossy-diffuse}(\omega_i, \omega_o) \approx f_\mathrm{dielectric}(\omega_i, \omega_o)  +  \bigl(1 - E_\mathrm{dielectric}(\omega_o)\bigr) \,f_\mathrm{diffuse}(\omega_i, \omega_o) \ ,
+f_\textrm{glossy-diffuse}(\omega_i, \omega_o) \approx f_\mathrm{dielectric}(\omega_i, \omega_o)  +  \bigl(1 - E_\mathrm{dielectric}(\omega_o)\bigr) \,f_\mathrm{diffuse}(\omega_i, \omega_o) \ .
 \end{equation}
-where the albedo of the Oren-Nayar lobe $f_\mathrm{diffuse}(\omega_i, \omega_o)$ must be taken to be $\mathbf{C} =$ **`base_weight`** * **`base_color`** in order to satisfy equation [glossy_diffuse_albedo_constraint].
+In this approximation, the $\boldsymbol{\rho}$ parameter of the Oren-Nayar model in Equation [FON_brdf] and Equation [EON_comp] must be simply $\boldsymbol{\rho} = \mathbf{C}$ in order to satisfy equation [glossy_diffuse_albedo_constraint].
 
 If a reciprocal formulation is desired, the classic Ashikhmin-Shirley or Kelemen model ([#Ashikhmin2000], [#Kelemen2001], [#Kulla2017], [#Kutz2021]) form is
 \begin{equation}
 f_\textrm{glossy-diffuse}(\omega_i, \omega_o) \approx f_\mathrm{dielectric}(\omega_i, \omega_o)  + \mathcal{N} \bigl(1 - E_\mathrm{dielectric}(\omega_i)\bigr) \bigl(1 - E_\mathrm{dielectric}(\omega_o)\bigr) \,f_\mathrm{diffuse}(\omega_i, \omega_o) \ ,
 \end{equation}
-where $\mathcal{N}$ is a normalization factor such that equation [glossy_diffuse_albedo_constraint] is satisfied. In the case of a zero roughness (i.e. Lambertian) Oren-Nayar base $\mathcal{N}$ can be tabulated in terms of the dielectric IOR and roughness ([#Kutz2021]), and the required albedo of the Lambertian base is equal to $\mathbf{C}$ as in the non-reciprocal albedo-scaling approximation. Extending this to the more general case of a non-Lambertian rough Oren-Nayar base would require adding the roughness dimension to the tabulation, and the required Oren-Nayar albedo will not simply equal $\mathbf{C}$. We leave the specific choice of model and these details to the implementation.
+where $\mathcal{N}$ is a normalization factor such that equation [glossy_diffuse_albedo_constraint] is satisfied. In the case of a zero roughness (i.e. Lambertian) Oren-Nayar base $\mathcal{N}$ can be tabulated in terms of the dielectric IOR ([#Kutz2021]), and the required $\boldsymbol{\rho}$ of the Lambertian base is equal to $\mathbf{C}$ as in the non-reciprocal albedo-scaling approximation. Extending this to the more general case of a non-Lambertian rough Oren-Nayar base would require adding the roughness dimension to the tabulation, and the required Oren-Nayar $\boldsymbol{\rho}$ will not simply equal $\mathbf{C}$. We leave the specific choice of model and these details to the implementation.
 
 
 Glossy-diffuse params        | Label             | Type     | Range            | Default             | Description
@@ -861,11 +891,11 @@
     ![Figure [coat]: Coat adds a secondary specular highlight and optional absorption tint](dummy)
 </div>
 
-In the full light transport, the observed color of the coated base is darkened and saturated due to multiple internal reflections from the inside of the coat, which causes light to strike the underlying material multiple times and undergo more absorption, and the observed **`coat_color`** tint also darkens as the incidence angle changes due to the change in path length in the medium. Also, the presence of a rough coat will increase the apparent roughness of the BSDF lobes of the underlying base. We assume that in the ground truth appearance, all these physical effects are accounted for. In the following sub-sections, we detail recommendations for implementation of them. [^porosity]
+In the full light transport, the observed color of the coated base is darkened and saturated due to multiple internal reflections from the inside of the coat, which causes light to strike the underlying material multiple times and undergo more absorption, and the observed **`coat_color`** tint also darkens as the incidence angle changes due to the change in path length in the medium. Also, the presence of a rough coat will increase the apparent roughness of the BSDF lobes of the underlying base. We assume that in the ground truth appearance, all these physical effects are accounted for. In the following sub-sections, we detail recommendations for implementation of them.
 
 ### Darkening
 
-Figure [coat_darkening_grid] shows the physically-correct change in appearance (at normal incidence) of a textured diffuse base with a wood texture and smooth clear-coat as the IOR of the clear-coat is varied, exhibiting darkening due to the internal reflections which increases as IOR increases.
+Figure [coat_darkening_grid] shows the physically-correct change in appearance (at normal incidence) of a textured diffuse base with a wood texture and smooth clear-coat as the IOR of the clear-coat is varied, exhibiting darkening due to the internal reflections which increases as IOR increases. [^porosity]
 
 ![Figure [coat_darkening_grid]: Diffuse base with a clear-coat, for coats of different coat IOR $\eta_c$](images/coat_darkening_grid.png width="60%")
 
@@ -1343,13 +1373,13 @@
 
 [#Burley2018]: Brent Burley, David Adler, Matt Jen-Yuan Chiang, Hank Driskill, Ralf Habel, Patrick Kelly, Peter Kutz, Yining Karl Li, Daniel Teece. *The Design and Evolution of Disney’s Hyperion Renderer*, ACM TOG (2018).
 
+[#Dupuy2023]: Jonathan Dupuy, Anis Benyoub. *Sampling Visible GGX Normals with Spherical Caps*, High-Performance Graphics (2023).
+
 [#Elias2001]: Mady Elias, Lionel Simonot, and Michel Menu. “Bidirectional reflectance of a diffuse background covered by a partly absorbing layer”. In: Optics Communications 191.1 (2001).
 
 [#d'Eon2021]: Eugene d’Eon. *A Hitchhiker’s guide to multiple scattering: Exact Analytic, Monte Carlo and Approximate Solutions in Transport Theory* (2022).
 
-[#d'Eon2023]: Eugene d’Eon. *Student-T and Beyond: Practical Tools for Multiple-Scattering BSDFs with General NDFs*, ACM SIGGRAPH Talks (2023).
-
-[#Dupuy2023]: Jonathan Dupuy, Anis Benyoub. *Sampling Visible GGX Normals with Spherical Caps*, High-Performance Graphics (2023).
+[#Fujii2012]: Yasuhiro Fujii. [*A tiny improvement of Oren-Nayar reflectance model*](https://mimosa-pudica.net/improved-oren-nayar.html) (2018).
 
 [#Georgiev2019]: Iliyan Georgiev, Jamie Portsmouth, Zap Andersson, Adrien Herubel, Alan King, Shinji Ogaki, and Frederic Servant. *Autodesk Standard Surface*, Autodesk white paper (2019).
 
@@ -1430,6 +1460,38 @@
 
 [^anisotropy_g]: Technically, $g$ is the _mean cosine of deflection_ of the phase function, which is not specific to the Henyey--Greenstein phase function model [#d'Eon2021].
 
+[^Oren_Nayar_formula]: The $s$ term is given by (with normal $\mathbf{n}$):
+\begin{eqnarray*}
+  s = \omega_i \cdot \omega_o - (\mathbf{n} \cdot \omega_i) (\mathbf{n} \cdot \omega_o) \nonumber = \cos(\phi_i-\phi_o) \,\sin\theta_i \,\sin\theta_o
+\end{eqnarray*}
+and
+\begin{equation}
+  \frac{1}{t} =
+    \begin{cases}
+      1                                                       & \text{if $s \le 0$}\\
+      1 / \max\left(\mathbf{n} \cdot \omega_i, N \cdot \omega_o\right) & \text{if $s > 0$ .}
+    \end{cases}
+\end{equation}
+The $A, B$ coefficients are
+\begin{eqnarray}
+  A &=& \frac{1}{1 + \left(\frac{1}{2} - \frac{2}{3\pi}\right)\sigma}\ , \nonumber \\
+  B &=& \sigma A \ .
+\end{eqnarray}
+
+[^Oren_Nayar_albedo]: The albedo of the [#Fujii2012] form of the Oren-Nayar BRDF is given by $E_\mathrm{ON}(\omega) = \rho\,\hat{E}_\mathrm{ON}(\omega)$, where
+For unit $\rho$,
+\begin{eqnarray}
+\hat{E}_\mathrm{ON}(\omega) = A + \frac{B}{\pi} G(\omega) \ ,
+\end{eqnarray}
+where (in spherical polars)
+\begin{equation}
+G(\theta) = \sin\theta \Bigl( \theta - \sin\theta \cos\theta \Bigr) + \frac{2}{3}\tan\theta \Bigl(1 - \sin^3\theta - \cos\theta\Bigr) \ .
+\end{equation}
+The corresponding albedo averaged over the hemisphere is given by
+\begin{equation}
+\langle\hat{E}_\mathrm{ON}\rangle = A + \left(\frac{2}{3} - \frac{28}{15\pi}\right) B \ .
+\end{equation}
+
 [^porosity]: In reality, coats can also darken the underlying surface due to a different mechanism where the coat modifies the Fresnel factor of the base due to the coat material filling in air gaps between granules or threads of a porous base material, which reduces the relative IORs at the internal interfaces. This occurs e.g. on adding water to sand or fabric, or adding a penetrating wood finish. We assume here that this effect explicitly does _not_ occur, at present, since we do not have enough knowledge about the properties of the underlying substance to model it. We can only safely assume that the first mechanism of darkening, i.e. internal reflections, occurs.
 
 [^avg_fresnel]: The hemispherical (or average) albedo of the Fresnel factor is defined as $E_F(\eta) \equiv 2 \int_0^1 F(\mu, \eta)\,\mu\,\mathrm{d}\mu$. This can be tabulated, or one can use the convenient analytical approximation [#d'Eon2021] (accurate to within 0.2% in the limited range $\eta \in [1,3]$):