Demystifying the He, Torrance, Sillion and Greenberg BRDF

By Romain Pacanowski

The BRDF Model

This model was introduced by He, Torrance, Sillion and Greenberg at Siggraph in 1991. The model relies on wave-optics, more precisely on the Beckman-Kirchoff theory and is known to be very difficult to use or to implement. In this post, we provide our implementation as well as some details about the model itself, which is not that complex to implement, at least in its unpolarized version!

In the supplemental material of our paper we compare our new BRDF Model with the unpolarized version of the He et al. model that is presented in the appendix of their paper (Equations 69 to 86). Our implementation omits the dirac part (Equation 74) of the model but this should be easy to add if you need it. Finally, as pointed on Westin’s page note that there is an error in equation (25) in the original paper. I have put the errata here as well.

Using the provided code

The C++ source code of our implementation is available here. You may use this code for research, academic or educational purposes only.

If you use this code for your research, do not forget to cite our paper using the following bibtex entry.

For industrial application, contact [Romain Pacanowski](mailto:romain DOT pacanowski AT inria DOT fr) and [Nicolas Holzschuch](mailto:nicolas DOT holzschuch AT inria DOT fr).

Remarks about our implementation

The He et al. BRDF model has several physical parameters:

  • A diffuse coefficients (_diffuse in the provided code)
  • The Index of Refraction as complex number represented by two values in the provided code:
    • _eta
    • _k
  • The standard deviation (\(\sigma_0 \) ) of the surface heights distribution
    • represented by _sigma_0 in the code
  • The autocorrelation length ($\tau$) of the surface heights distribution
    • represented by _tau in the code.

We tried to make the code explicit and well documented but do not hesitate to contact us if you have any questions. As usual with BRDF implementation some attention must be paid to avoid numerical issues.

In addition to the physical parameters described above, the model has additional parameters that control the accuracy of the computation.

  • Numerical Parameters
    • Number of terms (variable _m in He.hpp) used to evaluate the infinite sum of Equation 78
    • Threshold to solve the implicit Equation 81.

The diffraction term

Equation 78: $$ D=\frac{\pi^{2}\tau^{2}}{4,\lambda^{2}},\sum^{\infty}_{m}, \frac{g^{m},e^{-g}}{m!,m},\exp\Big(\frac{-v^{2}_{xy},\tau^{2}}{4m}\Big) $$ represents the most important part of the model regarding the diffraction contribution. This equation can already be found in the monograph of Beckman and Spizzichino (page 86 Eq.35).

In our code, _m is set to 10 by default and this sum is implemented around lines 204–212 in the file He.cpp.

Theoretically, _m depends on $g$ and more precisely on the ratio between $\sigma$ (surface characteristic) and the wavelength $\lambda$ of the light (cf. Equation 79).

When the surface is optically smooth (i.e., $\frac{\sigma}{\lambda} \ll 1$), a small value for _m can be used because the infinite sum converges quickly (more on that in a future post). In the context of the MERL-MIT BRDF database the surface roughness of each sphere is no available making it difficult to assess a “good” default value for _m.

To speedup the evaluation of the BRDF at runtime we precompute the factorial terms and store them in a vector (cf. std::vector<double> _factoTerm; in He.hpp line 62) initialized inside the class constructor (cf. He.cpp lines 50–57).

Solving Equation 81

The implicit Equation 81 of the original paper $$ \sqrt{\frac{\pi}{2}} z_0 = \frac{\sigma_0}{4} (K_i + K_r),\exp\bigg(-\frac{z_0^{2}}{2\sigma_0^{2}}\bigg) $$ is solved each time the BRDF needs to be evaluated. This is another reason why this BRDF model is computationally expensive.

In the provided code, this is done in the method (cf. lines 110–126 in the file He.hpp)

that uses a threshold value (hard coded 1e-5) to stop the numerical iteration. The threshold can be avoided by fixing the number of iterations to the detriment of accuracy.

Once $z_0$ has been computed, it is used to compute $\sigma$ as defined in equation 80: $$ \sigma= \sigma_0\ \Bigg(1+\Big(\frac{z_0}{\sigma_0}\Big)^{2}\Bigg)^{-0.5} $$

This equation is evaluated in the provided code at the line 190 in He.cpp.

Conclusion

The main difficulty we faced, was not the implementation of the model itself, but rather reproducing the fitting results from Ngan et al. published at EGSR 2005. More details about this issue are giving in our supplemental material.

References