1. Introduction

In rendering a scene, often ray tracing is performed. During ray tracing, starting either at the source (a light in the room) or the observer (a pixel on a camera in the scene), a ray is traced until it intersects with an object or medium. Upon intersection, the ray is then scattered in some manner that primarily depends upon material scattering properties. This work modifies a commonly used optical scatter distribution by folding in certain aspects of a more precise wave optics model while maintaining the simplicity of a closed-form expression that has been studied in previous work [1], and expanded upon here. Then, this work presents a computationally efficient mechanism for the most complex element of the novel BRDF computation: calculation of the polarization factor, $Q$. This optimization sped up computation of $Q$ to be nearly 2 times faster than a naive implementation, has the advantage of $Q$ being inherently polarimetric in both incident and scattered directions, and not requiring computation of a geometric attenuation term. This is in contrast to existing polarimetric models such as proposed by Priest [2,3], which are more computationally cumbersome and do not have the benefits of $Q$ at grazing angles, while still requiring a geometric attenuation term. Even for non-polarimetric rendering, which is the case for this work, use of $Q$ improves radiometric accuracy of the BRDF at grazing angles.

Optical scatter descriptions have existed for several decades. The standard definition of a description of optical scatter known as the bidirectional reflectance distribution function (BRDF) was defined by Nicodemus [4] as the ratio of scattered radiance, $L_s$, to incident irradiance, $E_i$,

The BRDF carries units of $sr^{-1}$, and ranges from 0 to $\infty$ in value. However, the BRDF is related to unitless directional reflectance by integrating over all scattered angles (for Directional-Hemispheric Reflectance, DHR) or over all incident angles (for Hemispherical-Directional Reflectance, HDR), taking into account the power projection area. The HDR is given as,

Several different models may be used, but this paper deals with two particular physics-based models: microfacet models, which assume geometric optics, and scalar wave optics models. This work begins by summarizing previous general results. To help motivate the problem, one material in Matusik’s Mitsubishi Electronics Research Laboratory (MERL) dataset [13] is presented. This dataset was used because it contains complete angular BRDF measurements which do not require fitting to a model in order to compute the unitless reflectance, but when one does so, problems at grazing angle are made clear. That work is then extended by relating the scalar wave optics model coordinate system to the microfacet model, summarizing prior work [14,15]. Then, a prior improvement in high fidelity grazing angle BRDF data is presented [1]. The primary downside of previous work, however, is the apparent substantial computational overhead present in using the polarization factor, $Q$, and the limited scope of the prior BRDF fits (primarily to paints). This work illustrates use of the Modified Cook-Torrance (MCT) model applied to a different set of BRDF data, then presents a non-trivial method of computing $Q$ in a more efficient manner than a naive implementation, easily adapted to scene rendering, and easily adapted to other programming languages. This novel computation presented here will assume all variables are real-valued, not complex-valued, which enables simple porting of the algorithm to a wide range of programming languages. The computation is also readily adapted to use dot products of unit vectors in place of raw angles, which makes it easily portable to scene rendering software. It is anticipated that these results will yield more physically accurate BRDF modeling, particularly for scenes with large incident or scattered angles where the microfacet model is more prone to error, without the apparent substantial additional computational overhead suggested by prior work.

2. Background

In this paper, two major background areas are developed. The first will cover general background on microfacet and wave optics BRDFs, while the second will cover specific background on particular models and BRDF fits analyzed in prior work, which are then extended in this paper.

2.1 Microfacet BRDF models and wave optics

Optical scatter is inherently a wave optics problem, where diffraction may be a significant factor in describing the distribution of light when striking a surface. Although wave optics models such as Rayleigh-Rice [16] and Generalized Harvey-Shack (GHS) [17] exist, they are either not general enough to apply to the wide variety of materials encountered in a scene (in the case of Rayleigh-Rice), or unwieldy for certain computations such as remote sensing or efficient scene rendering (in the case of GHS). Recent work in computer graphics has began to incorporate wave optics effects to some degree, focused on specular microgeometry [18]; here, we discuss diffuse reflection improvements. Due to the computational requirements for a wave optics solution geometric optics approximation from the microfacet model is far more widely used in the computer graphics and scene generation communities. Some common examples of popular microfacet models include: Blinn-Phong [7], Cook-Torrance [19], Ward-Duer [20,21], Ashikhman-Shirley [22], and several other models [3,5,23,24]. Some of these models arise from computer graphics research directly, while others are of use in other applications such as remote sensing or other scene generation applications. In prior work [9], several models were rewritten in a common form to observe similarities, differences, and best practices for microfacet models, as well as to develop a framework in which to easily design new microfacet models using $Q$.

In previous work, the microfacet BRDF model was written as, [9]

Fresnel reflectance is derived in many optics and electromagnetism physics textbooks, but assumes an infinitely large, perfectly smooth surface. For semi-rough surfaces, this term is often still retained despite the assumption of a perfectly smooth surface as a boundary condition used in its derivation. (For optical wavelengths, which are on the order of $10^{-7}\, m$, the infinitely large approximation is typically a good approximation, but assuming perfectly smooth is suspect, especially for rough surfaces. By contrast, the polarization factor, $Q$, described later, was derived assuming an infinitely large, semi-rough surface [16,27,28].) $D(\theta _h)$ is the microfacet distribution function in terms of $\theta _h$, the half angle vector in Rusinkiewicz coordinates, which specifies the microsurface normal’s angle relative to the overall surface normal. $\sigma (\theta _i, \theta _s)$ is a cross-section conversion term. In physics, the cross section is most often expressed in terms of a spherical scatterer; when converting to planar scatter more appropriate for a geometric optics BRDF model, this $\sigma$ cross section conversion term arises. It is equal to $(4 \cos \theta _i \cos \theta _s)^{-1}$, and has been experimentally validated [29]. This term, though physically valid, is mathematically problematic at grazing angle ($\theta _i\rightarrow 90^\circ$ or $\theta _s\rightarrow 90^\circ$), since $\sigma$ diverges to $\infty$ at grazing angle. Mathematically, $G$ compensates for this deficiency by converging to $0$ more rapidly than the cross section conversion term diverges. Microfacet models which include the cross section conversion term typically, though not always, include $G$ for this reason [9].

The primary issue with the microfacet model at grazing angle is related to the interplay between the cross section conversion term, $\sigma (\theta _i,\theta _s)$, which is physically sound but mathematically problematic at grazing angle, and the geometric attenuation term, $G$, which serves an important mathematical function so the BRDF remains finite at grazing angle but is only physically sound for surfaces where the roughness is much larger than the wavelength of light; such surfaces are very rough scatterers, which do not exhibit a strong preferential scatter direction. To illustrate this problem, consider Matusik’s MERL database [13], which is popular in the computer graphics community for scatter research due to being freely available, densely sampled BRDF data. Using Eq. (2) and linear interpolation between data points, the unitless HDR can be computed; see Fig. 1 for the MERL brass.binary data file red, green, and blue channel HDR results as a typical example of the HDR shape for materials in the MERL database [13]. In physics, it is well-known that the unitless reflectance of a surface approaches 1 as the viewing angle approaches grazing, which does begin to occur for HDR in the MERL database as well; however, the HDR then drops rapidly at grazing angle rather than following the expected physical behavior. This occurs even for metals (as shown in Fig. 1), which have minimal volumetric scatter at visible wavelengths. This problem motivated a comparison between microfacet and wave optics BRDF models to determine whether there may be a better approximation to use in a closed-form model in place of $\sigma$ and $G$ in the microfacet BRDF model; this is the subject of the next section.

 

Fig. 1. HDR from MERL brass file for the red, green, and blue channels

Download Full Size | PDF

2.2 Relating microfacet and wave optics BRDF models

Another class of physics-based BRDF models is wave optics models. Light is an electromagnetic wave, and as a surface feature (such as either the variation in surface height or the correlation length in the surface plane) approaches the wavelength of light, diffraction becomes an increasingly important factor in modeling optical scatter precisely. One of the most general scalar wave optics BRDF models was developed in an approximate sense by Harvey and Shack [30], then eventually extended by Krywonos [17,31]; it is known as the Generalized Harvey-Shack (GHS) model.

In the GHS model, as in many other wave optics models, direction cosine coordinates popular in Fourier optics texts such as Goodman [32] ($\alpha$, $\beta$, $\gamma$) are used instead of Rusinkiewicz ($\theta _h$, $\phi _h$, $\theta _d$, $\phi _d$) coordinates [26] used in most microfacet models. In a scalar wave optics model that incorporates diffraction effects, in order to obtain a general solution, the Fourier transform of the surface must be computed to obtain the Power Spectral Density (PSD). This results in a substantially more complex computation for the optical scatter distribution, which is the primary reason it is not traditionally used in computer graphics, remote sensing, or other applications where optical scatter is only one of several computationally demanding steps. However, it is popular in modeling stray light or as a non-contact method to evaluate the details of a surface (such as in manufacturing a lens [33]).

Even if the complexity of the Fourier transform can be overcome, such as by assuming the surface follows Gaussian statistics, the general solution is not a single closed-form expression, but rather an infinite summation. This was also developed by Krywonos [17,31], but to better aid the reader in following the relationship to microfacet (Rusinkiewicz) coordinates, this paper will use the equivalent form in prior work [15], given as,

The remaining term in Eq. (4), $\nu _{xy}$, is related to the direction cosine space used in Fourier optics, and for isotropic materials it is given by a relation developed in our prior work, [15]

To convert the ASF to BRDF, Eq. (4) must be multiplied by either Fresnel reflectance, which assumed a perfectly smooth surface in its derivation, or (as is done in the Rayleigh-Rice smooth surface BRDF model) by the polarization factor, $Q$, which was derived from a perturbation of Fresnel reflectance for a surface that is not perfectly smooth [16,27,28]. The polarization factor, $Q$, is split into four different components relating the polarization of the incident and scattered light; for unpolarized light, the equations are summed. The equation for $Q$ is given as, [16]

It is clear from this description that applying a scalar wave optics model directly would be extraordinarily taxing from a computer graphics, scene generation, or remote sensing perspective, due particularly to the Fourier transform (if the surface statistics are non-Gaussian) and the infinite summation (even if the surface statistics are assumed to be Gaussian). For this reason, a scalar wave optics model is not used directly. Instead, the following relationship between microfacet coordinates and direction cosine coordinates was proven: [15]

Although Ewing frequently fitted BRDFs measured at wavelengths other than visible light (UV and IR in addition to visible light), the fundamental physics of optical scatter is unchanged. Equation (5), derived from wave optics, shows that increasing wavelength is mathematically equivalent to reducing RMS surface height variation; that is, as wavelength increases, it is equivalent to having a surface that appears more highly polished.

The MCT model could be thought of as a hybrid model that incorporates parts of scalar wave optics alongside a microfacet distribution function. The scalar wave optics incorporated is in the term below that relates $Q/2$ to other microfacet terms, and is given as [1,14,15],

3. MCT BRDF fits to novel materials

In our previous work, a precise measurement of BRDF at a variety of wavelengths was made for several materials [36]. As discussed previously and shown in Eq. (5), although that work is primarily in the IR, increasing wavelength is mathematically and physically equivalent to decreasing surface roughness from a scalar wave optics scattering perspective [36]. The original intent of that study was to examine the wavelength-dependent behavior of the BRDF, and predates the development of MCT, but since the data represents another source of high-fidelity measurements of the same surface at different wavelengths, this time on other materials not tested by Ewing, the analysis is performed here. The data was fitted to both the original Cook-Torrance model [19] and to the MCT model shown in Eq. (9), both using a Gaussian microfacet distribution. This was performed for only the semi-rough samples measured, since the microfacet BRDF model is known to be much less accurate for polished samples [14,15]. To summarize key results in the original measurements, these data were taken using a Schmitt Measurement Systems Complete Angle Scatter Instrument (CASI), which is capable of measuring up to 9 orders of magnitude in the BRDF and capable of separating emission from reflection in the measured data. The measured data typically has a relative error of under 8% [36]. These results lend further credence to the use of $Q$ (and thus the MCT Model) for microfacet BRDF fitting; the following section will then present an efficient method to implement $Q$ in code.

Figure 2 shows the original Cook-Torrance model fitted to the BRDF data of a rough aluminum sample at 6.4 $\mu m$. To normalize to differences in overall unit reflectance from one sample to the next, the BRDFs presented were divided by the overall reflectance, $\rho$. Each different colored data point represents a different incident angle. Since the fits did not differ significantly in the backscatter direction, only $\theta _s > 0$ is plotted. The best fit parameters obtained were as follows: $\rho _d = 0.0318, \rho _s = 0.934, \sigma _g = 0.293, n = 8.68\times 10^3,$ and $\kappa = 1.71\times 10^4$ ($\varepsilon = |\tilde {n}|^2 = 3.67\times 10^8$), where $n$ and $\kappa$ are particularly suspect values obtained, although it is common in the microfacet community to use best fit parameters instead of actual parameters for $n$ and $\kappa$, an approach used here as well. Reported values for Al at this wavelength are $n = 12.7$ and $\kappa = 58.8$, $\varepsilon = |\tilde {n}|^2 = 3.62\times 10^3$. [37] Fitting with actual $n$ and $\kappa$ values used instead is left as future work, once the closed-form modeling is more refined.

 

Fig. 2. BRDF divided by unitless DHR for rough Al sample at 6.4 ${\mu m}$, fitted to the original Cook-Torrance BRDF model, using $ln$ error metric

Download Full Size | PDF

Figure 3 contains identical data and maintains a Gaussian distribution function identical to that used by Cook-Torrance, but revises the model using the replacement suggested by Ewing’s paint sample analysis; that is, the MCT model presented in Eq. (9) is used in place of a traditional Cook-Torrance model. The best fit parameters obtained were as follows: $\rho _d = 0.0318, \rho _s = 6.68, \sigma _g = 0.216, n = 0.910$ and $\kappa = 0.740$ ($\varepsilon = |\tilde {n}|^2 = 1.38$). As briefly mentioned earlier, this modification primarily differs as $\theta _i$ or $\theta _s$ increases. This can be seen from the resulting model plots, where the difference in model performance is most pronounced for large $\theta _s$ and when $\theta _i = 50^\circ$, the largest incident angle for which data was available.

 

Fig. 3. BRDF divided by unitless DHR for rough Al sample at 6.4 ${\mu m}$, fitted to Modified Cook-Torrance BRDF using Q in place of F, G, and $\sigma$, using ${\ln }$ error metric

Download Full Size | PDF

To compare the performance of these two models where they differ from each other the most, and to illustrate the result on a linear plot rather than semi-logarithmic plot, Fig. 4 contains the $\theta _i = 50^\circ$ plot from the previous two figures, to include data points with uncertainty, and illustrating both the original Cook-Torrance and MCT models. From this plot, it is clear that the MCT model is a noticeably better fit to the data.

 

Fig. 4. BRDF divided by unitless DHR for rough Al sample at 6.4 ${\mu m}$ for $\theta _i=50^\circ$ on a linear scale, with Cook-Torrance and MCT best fits.

Download Full Size | PDF

As another example of the improvement in BRDF modeling at grazing angle, an Infragold sample was also measured, this time at a wavelength of 5.2 $\mu m$, with the same process illustrated above for aluminum applied; see Fig. 5. Again, it is clear that the MCT model outperforms the standard Cook-Torrance model, with the primary difference being the use of the polarization factor, $Q$, via Eq. (11).

 

Fig. 5. BRDF divided by unitless DHR for rough Infragold sample at 5.2 ${\mu m}$ for $\theta _i=50^\circ$ on a linear scale, with Cook-Torrance and MCT best fits.

Download Full Size | PDF

There were several other BRDF datasets for which this analysis was performed. As mentioned, the dataset in this paper was originally obtained when exploring the BRDF as a function of wavelength, so the same material was measured multiple times with lasers at different wavelengths as the incident light source. In this work, each of these can be treated as an independent sample, where the effect of increasing wavelength is, from a scalar wave optics perspective, equivalent to increasing the smoothness of the surface for the same general material characteristics. To summarize the results obtained, the norm of the relative difference between the BRDF data and each of the two models was computed; a smaller value indicates a better overall match to the data. The results are summarized in Table 1. This norm of the relative difference is computed as,

Table 1. Comparison of the norm of the relative difference, $d$, when fitting BRDF data to the Cook-Torrance model using the ln error metric, and the relative difference using the MCT model, for multiple diffuse samples, at all wavelengths measured.

View Table

In each case, the norm of the relative error for the MCT model is lower than the norm of the relative error for the original Cook-Torrance model, further validating the better performance of the model presented in Eq. (9) and further bolstering our previous result [1].

4. Efficient implementation of $Q$

The primary barrier to implementation of the MCT model given in Eq. (9), or other variant modified microfacet models that use the substitution in Eq. (11), is the complexity for the equation for the polarization factor, $Q$, presented in Eq. (7), particularly when compared to the Fresnel equation. To overcome that barrier, this section proposes an efficient form of computation for $Q$ that runs more efficiently than a naive implementation. (Note that this paper actually proposes an efficient computation for $Q/2$, since it appears that way in the MCT model, and because $Q/2=F$ when $\theta _i=\theta _s$ and $\Delta \phi =\pi$.) The code is available online in Code 1 (Ref. [38]), and is listed at the end of this document as well for ease of the reader to reference.

When timing implementations of $Q$ and Fresnel reflectance in this section, a single index of refraction was used with $\theta _i$ ranging from $0^{\circ}$ to $89^{\circ}$ in $1^{\circ}$ increments, and $\phi _i=0^{\circ}$. For each incident angle, $\theta _s = 0^{\circ} - 89.9^{\circ}$ in $0.01^{\circ}$ increments, and $\phi _s = 0^{\circ} - 359.9^{\circ}$ in $0.1^{\circ}$ increments. This resulted in $32,367,600$ executions of the computation for each incident angle, or $2,913,084,000$ total executions. The norm of the relative difference between the efficient coding of $Q$ presented in this paper and a naive but easy-to-read coding of $Q$ based on Equation (7) was computed to be less than $10^{-13}$, which was near machine precision, validating the accuracy of the code. The Matlab code was vectorized in all cases (naive Q, fast Q, and Fresnel). The Matlab profiler was used to generate timing information. The execution time for Fresnel reflection in Matlab was $9.512\,s$ (or $32.62\,ns$ per execution). Using a naive implementation of $Q$ with no optimization, the code took $24.76\,s$ to run ($84.90\,ns$ per execution), but the optimized version presented in this paper reduces that computational overhead to only $13.33\,s$ total ($45.71\,ns$ per angle), less than double the amount of time required to compute Fresnel reflection using a non-optimized method. The resulting speed-up is nearly a factor of 2 using this optimized form. Recall that when using $Q$, the geometric attenuation term, $G$, present in most microfacet models, does not need to be computed, nor does the $1/\cos ^4\theta _h$ term that is often included in the microfacet distribution function, $D$, as these terms cancel out with terms presented in Eq. (11). The cross section conversion term also cancels out. This results in further time savings for each angle computed when used in a full BRDF model, beyond what is presented here. Since those modifications are trivial, the focus of this section will be on efficient implementation of $Q$.

The version presented here assumes the index of refraction of the scattering material is complex, but does not assume use of complex math. It only assumes use of floating-point numbers (of the desired level of precision, float or double), basic arithmetic operations, basic trigonometric functions (including a 2-parameter $\arctan$ function defined over the entire $[0, 2\pi )$ interval), and square root, are defined. This code is thus easily translated into a compiled language, such as C++, when written in the form presented here. In C++ or other similar languages, each variable could be defined as a constant value, as their values are not changed elsewhere in this code. It would be trivial to instead adapt the function to use the dot product of unit vectors for $\cos \theta _i$, $\cos \theta _s$, and $\cos \Delta \phi$ in place of the angles being explicitly passed in, if the function were used in scene generation software; other computations performed do not require the angles explicitly, and build on the cosine of these angles instead.

In the code presented in Code 1 (Ref. [38]), recurring intermediate computations are made using $\eta =n_t/n_i$ and $\eta _k=k_t/n_i$ as the baseline; that is, $\tilde {\eta }=\tilde {n}_t/n_i$, where ∼ indicates a complex value, but it is explicitly split into real ($\eta$) and imaginary ($\eta _k$) parts. Some of the optimizations are trivial, such as storing temporary variables when values are used multiple times or multiplying a value by itself instead of raising it to the second power. For example, since the magnitude squared is required to compute $Q$, there are many instances where $\Delta \eta ^2 = \eta ^2 - \eta _k^2$ is used, as well as the cross term $\eta _{cp} = \eta \eta _k$ (or minor variations of this value, such as $2\eta _{cp}$). These optimizations are not discussed in detail because they are trivial.

However, several of the optimizations required clever manipulation of complex variables in Euler space to add, multiply, or take square roots of vectors in the complex plane through more efficient means than using complex number software packages. When writing complex numbers, they can be written in Cartesian space as $x + iy$, but can be equivalently rephrased as $r \exp {(i\theta )}$, with the latter representation having the benefit of being easy to raise to a power, such as squaring or taking the square root, as is done frequently in computing $Q$ via Equation (7). When performing such a manipulation, $r^2 = x^2 + y^2$ and $\theta = tan^{-1} (y/x)$. The square root then is the square root of r and half the angle of $\theta$, as opposed to two square roots (real and complex component) required in Cartesian space. If the quantity is squared later prior to being added to another term, a square root can be avoided entirely (as is true in some cases here, particularly with terms in the denominator of $Q$). Additionally, the magnitude squared of the exponential term vanishes for the final result, which also greatly reduces the computational complexity when using this notation.

The building blocks of the computation in Euler space are the computations of the variables $\alpha$ and $\beta _{i,s}$, which are angles in the complex plane. The angle $\alpha$ arises from representing $\tilde {n}^2=\varepsilon$ in the Euler notation of $b \exp (i\alpha )$, which can be algebraically reduced from the more obvious form of $\sqrt {(\eta ^2-\eta _k^2)^2 + (2\eta \eta _k)^2}$ to $b=\eta ^2 + \eta _k^2$. The angle in complex space is given from the original expression as,

From here, the terms $a^4_{i,s} = (\Delta \eta ^2-\sin ^2\theta _{i,s})^2 + \eta _{cp}^2$, $a^2_{i,s}$, and $a_{i,s}$ are computed as intermediate values in the computation of various $Q$ components. Then, these building blocks are used to define $c^2_i$ and $c^2_s$ terms, which are both found in the denominator of the $Q_{ss}$ expression, and once each in the $Q_{sp}$ and $Q_{ps}$ expressions. The related terms that are multiplied by $\tilde {n}^2$ in the denominator of $Q_{sp}$, $Q_{ps}$, and $Q_{pp}$ are the $e^2_i$ and $e^2_s$ terms. The terms $w$, $x$, and $p^2$ in the code build up portions of the numerator present in the $Q_{pp}$ term. When making these substitutions, computing the magnitude squared using floating-point numbers, and noting that $|\tilde {n}^2-1|^2=(\Delta \eta ^2-1)^2 + (2\eta _{cp})^2$, the computation reduces to,

5. Conclusions

In this paper, a MCT BRDF model previously developed was presented from theoretical principles by comparing a geometric optics microfacet model to a more accurate, but computationally cumbersome, scalar wave optics model (GHS) [9,15]. The modified model is general enough to be adaptable to other microfacet models as well, to obtain more accurate results as either $\theta _i$ or $\theta _s$ approaches grazing angles, while retaining the closed-form expression for BRDF that is not present in a full GHS scalar wave optics model. The MCT model utilizes the polarization factor in place of the cross section conversion term, $\sigma$, which diverges at grazing angles. This novel model additionally replaces Fresnel reflection (which assumes an infinitely smooth surface) with a perturbation using the polarization factor $Q$, and replaces the need for a shadowing and masking (geometric attenuation) term, $G$, since optical scatter is fundamentally a wave optics phenomenon and not a geometric optics phenomenon. Previous results in fitting MCT were summarized briefly. Then, this work extended fitting the MCT model to a dataset previously collected for a different purpose [36]. The results presented further substantiating evidence supporting the use of the MCT model, and when taken in conjunction with the prior result [1], suggest the model is improved for a wide variety of materials. Finally, an efficient and inherently polarimetric implementation of $Q$ was presented, which executes in Matlab about 2x faster than a naive implementation and can be readily ported to multiple programming languages. Although this work only examines unpolarized modeling, the extension to polarimetric modeling is straightforward; this is left as future work. Additionally, as noted when $Q$ was introduced, the $Q_{pp}$ term has a known issue of rising beyond the expected value of nearly 1 for metals when not at the specular reflection direction; solutions to this issue are left for future work.

This optimization to the model presented in Eq. (9), and the generic expression presented in Eq. (11) that can be applied to any microfacet BRDF model, will result in substantial increased computational efficiency in ray tracing for computer graphics and scene generation rendering, particularly for scenes where grazing angles are of substantial importance in rendering.