Parameters Sharing Multi-items Non-parametric Factor Microfacet Model for Isotropic and Anisotropic BRDFs

Abstract

The reflection model of object surfaces is an essential part in photorealistic rendering. While analytical models are unable to express all significant effects of all materials, we turn to data-driven models, which, however, cost a large amount of memory and more computational resources. The non-parametric factor microfacet model designed by Bagher et al. [1] is intended to solve these problems. In this paper, we present a new non-parametric factor microfacet model which has triple specular lobes but retains the original number of parameters by sharing factors among three color channels. The fitting method called AWLS is also extended to solve for the G factor, which makes the fitting process more robust. Moreover, we use the D factor of our model for importance sampling as in the case of analytical models and find it effective for the specular materials. Finally we generalize our model and our fitting method to fit the 150 anisotropic materials. With only 2010 parameters (8KB), it can reconstruct the original data (2 MB) well, which further proves the expressiveness of our microfacet model.

1 Introduction

In computer graphics, bidirectional reflection distribution function (BRDF) [2] has been extensively used to represent material appearance. For a specific wavelength $\lambda$, BRDF is a 4D function $f_{\lambda }(\theta _i, \phi _i, \theta _o, \phi _o)$, which returns the ratio of outgoing radiance to the incoming irradiance incident on the surface. All the notations used in this paper are shown in Fig. 1.

Fig. 1.

Notations used in this paper.

BRDF models are simply divided into analytical models and data-driven models. Analytical models give exact analytical forms with only a few parameters to fit different materials. The microfacet model is one kind of analytical models, which is derived from the microfacet theory. The microfacet theory [3] assumes that rough surface consists of adequate microfacets, which have the same reflection properties and whose orientations obey some kind of distribution. For a given wavelength $\lambda$ (or color channel), the basic structure of the microfacet model is simple:

$$\begin{aligned} \rho (\theta _h, \theta _d, \phi _d) = d + s(\frac{D(\theta _h, \phi _h)F(\theta _d)G_2(\theta _i, \theta _o)}{\cos {\theta _i}\cos {\theta _o}}) . \end{aligned}$$

(1)

If both the index of refraction and the extinction coefficient of the material of interest are known, then $F(\theta _d)$ is determined. The indices of refraction and the extinction coefficients of some materials can be checked in the handbook [4]. Additionally, the masking (or shadowing) factor G can be derived from D, which is first proposed by Smith et al. [5] and generalized by Brown et al. [6] and Walter et al. [7]. If a simple kind of masking-shadowing function $G_2$, e.g. $G_2(\theta _i, \theta _o)=G(\theta _i)G(\theta _o)$, is chosen, then the biggest difference among the microfacet models is normal distribution function $D(\theta _h, \phi _h)$.

There are lots of choices for $D(\theta _h, \phi _h)$. Some of the most commonly used distributions can be found in [7,8,9,10,11,12]. It should be noted that the GGX distribution is the same as the Trowbridge-Reitz distribution. The SGD distribution is a mix of the Trowbridge-Reitz distribution and the Beckmann-Spizzichino distribution. Generally speaking, the shape of the Blinn-Phong distribution is very close to the Beckmann-Spizzichino’s. The Trowbridge-Reitz, SGD and ABC distributions have a narrower peak with a stronger tail but the ABC distribution is only suitable for glossy surfaces. By observing measured data, we can design a more accurate normal distribution function. Actually, with the method Ashikmin et al. [13] introduced, we can arbitrarily design normal distribution functions regardless of whether the corresponding materials exist or not.

However, all the factors derived above are subject to some assumptions, which limit their capability to express some special reflection effects as well as the range of materials their corresponding models suitable to be used. In order to express all the reflection effects and not be bounded to the categories of materials, we can directly use the measured reflection data. However, acquiring accurate and dense measured reflection data is a challenging and time-consuming work. Fortunately, there is a large enough isotropic BRDF database called MERL [14], which has been analyzed by a lot of researchers. Moreover, UTIA¹ BRDF database [15] contains measured reflection data for 150 anisotropic materials. Although it is not that dense, it still can be directly used in a renderer with an appropriate interpolation method. However, directly using the measured data not only needs a large amount of memory storage, but also loses the flexibility for users to edit the material appearance and the simplicity for importance sampling.

One possible way to solve the problems aforementioned is using 1D arrays to substitute all the factors in the basic structure of the microfacet model. The model after substitution is called non-parametric factor microfacet model in Bagher et al. [1]. As for $G_2$, they used the simplest but inaccurate one, $G_2=G(\theta _i)G(\theta _o)$. Bagher et al. [1] used two slightly different models, the best of which was called Independent-G model. As its name suggests, the G factor is independent and is not derived from the D factor in the Independent-G model. Following the format of the MERL database, each of the factors $D(\theta _h)$, $F(\theta _d)$ and $G(\theta _i)$ (or $G(\theta _o)$), has 90 elements. The fitting objective is

$$\begin{aligned} E = \sum _j{w_j(\rho (\theta _h, \theta _d, \phi _d)_j-\rho ^*_j)^2} \end{aligned}$$

(2)

where $\rho ^*_j$ is a BRDF measurement and $\omega _j$ is the compressive weight. Those who are interested in the weighting scheme can find detailed formulations in [1].

By fitting the MERL measured data, each element of the 1D array and the other two coefficients can be determined. Then using the Independent-G model, original reflection data of each material is approximated with 816 floating-point parameters. It should also be noticed that the Independent-G model is always better than the compared analytical models. So we only use the Independent-G model to compare with our model in the experiments. In addition, we also compare our model with the so-called Bivariate model, which use 10$\times$ more parameters and each parameter has no physical meaning. Taking the images rendered with the original reflection data as the reference images, the PSNRs between the images rendered with the Independent-G model and the reference images are shown in Fig. 8 in [1]. Although the result seems good enough, there is still much room for improvement.

We change the structure of the Independent-G model from a single specular lobe to three specular lobes. However, by sharing factors among three color channels, our model only uses six extra parameters. Moreover, sharing factors conforms to the microfacet theory, which makes our model more explainable and intuitive. More importantly, it’s more suitable for user-editing and importance sampling. In addition, we still can use the simple and easily implemented fitting method which is called alternating weighted least-squares (AWLS) in Bagher et al. [1] to fit the data. The AWLS repeatedly updates each factor in sequence until convergence. Its basic idea is finding the zero point(s) of the first derivative of the fitting objective 2 in order to get the minimum of it. It is easy to prove that its second derivative is positive so the zero point must be the minimum. We do not introduce the AWLS in detail here, but in the next section, we extend it to solve for the G factor. Bagher et al. [1] thought the AWLS was not suitable for the G factor because of the non-trival factor dependency and used the GSS (golden section search) to guarantee to reduce the fitting objective. After extended, the AWLS is also able to reduce the fitting objective steadily while updating the elements of the G factor without using GSS. Finally, microfacet theory does not confine to isotropic materials. So we also use our modified model to fit the anisotropic data.

In the following sections, we first introduce our models for the isotropic and anisotropic measured data respectively. Then we talk about the extended AWLS for the G factor and the simple importance sampling method which helps dramatically decrease the number of samples needed to render some specular materials. After that, we show the fitting results and further analyses about the resulting parameters. The limitations of our model are concluded in the end of the paper.

2 Parameters Sharing Multi-items Model

2.1 Our Model for the Isotropic Materials

In the Independent-G model, the D, F and G factors are different for each color channel. Although such a model can fit the measured data well, it loses the original physical meaning. For example, the D factor represents normal distribution of object surfaces. Therefore, for different color channels (or wavelengths), the D factor is the same. The same conclusion can also apply to the G factor. Only the F factor, the fresnel factor, is dependent on the wavelength of the incident light so it is different for different color channels. Therefore, we share the D and G factors among different color channels, which helps reduce the number of parameters.

Furthermore, experiments conducted by Ngan et al. [16] indicated that two specular lobes helped reduce the fitting error. Lafortune et al. [10] used three specular lobes but Ngan et al. [16] indicated that it made the fitting process very unstable. However, with the extended fitting method AWLS, it is not a problem. In order to compare with the Independent-G model fairly, our model uses three specular lobes, which only adds 6 unavoidable parameters. Obviously, the total number of parameters our model uses is $3+3*3+90*3*3=822$. The resulting model is

$$\begin{aligned} \rho _C(\theta _h, \theta _d, \phi _d) = d_C + \sum _{m=1}^{m<=3}{s_C^m(\frac{D^m(\theta _h)F_C(\theta _d)G^m(\theta _i)G^m(\theta _o)}{\cos {\theta _i}\cos {\theta _o}})} \end{aligned}$$

(3)

where the C stands for R, G or B. In addition, it is necessary to notice that, because the D and G factors are shared among three color channels, solving for each element of these two factors uses all the samples from all three color channels.

2.2 Our Model for the Anisotropic Materials

For anisotropic materials, the D factor is dependent on both $\theta _h$ and $\phi _h$, and the G factor is dependent on both $\theta _i$ (or $\theta _o$) and $\phi _i$ (or $\phi _o$). Moreover, to be aligned with the format of the UTIA database, both the factor D and G should be a 6 $\times$ 48 2D array. Since the F factor is dependent on wavelength $\lambda$ and $\theta _d$, and there are only three color channels, the F factor is still a 90-element 1D array for each color channel. So the resulting model is

$$\begin{aligned} \rho _C(\theta _i, \phi _i, \theta _o, \phi _o) = d_C + \sum _{m=1}^{m<=3}{s_C^m(\frac{D^m(\theta _h, \phi _h)F_C(\theta _d)G^m(\theta _i, \phi _i)G^m(\theta _o, \phi _o)}{\cos {\theta _i}\cos {\theta _o}})} \end{aligned}$$

(4)

where the C stands for R, G or B.

Therefore, our model for the anisotropic materials needs $3+3*3+6*48*3*2+90*3=2010$ parameters (about 8KB). We use the tri-linear interpolation method to get the intermediate value for the D and F factor. A helpful trick in practice is to store the weights and positions for each sample before fitting in order to avoid transforming the parametric space over and over again.

2.3 The Extended AWLS for the G Factor

The AWLS fitting methods for D, F, d and s are introduced in detail in [1]. For the G factor, Bagher et al. [1] applied linear interpolation method in the $\theta _i$ and $\theta _o$ space while evaluating the G factor or the cosine factors. Firstly, they used a standard AWLS update step to solve for the $G(\theta _i)$ and $G(\theta _o)$ separately and then averaged them. After that, they used a simple Gaussian smoothing filter to remove high-frequency oscillations. If the first iteration did not reduce the fitting objective, they applied a GSS update.

In order to make the G factor smoother, we use a Catmull-Rom spline to interpolate the G factor and the cosine factors. We forward readers to the Pharr et al.’s book [17] about how to use a Catmull-Rom spline to interpolate 1D control points. So $G(\theta _i)$ (or $G(\theta _o)$) for arbitrary $\theta _i$ (or $\theta _o$) can be expressed as a weighted sum over four control points, where the weights and the particular control points depend on the $\theta _i$ (or $\theta _o$).

Then, we directly use the basic idea of AWLS to solve for the G factor. Assume we want to update the $k^{th}$ element of $G^1$ (updating $G^2$ and $G^3$ is totally the same), denoted by $G^1_k$. All the other factors are held constant, including d, s, $G^2$, $G^3$ and all the other elements of $G^1$. If evaluating BRDF of a specific parametric location with our model needs to use $G^1_k$, which means $G^1_k$ is one of the four control points when evaluating $G(\theta _i)$ or $G(\theta _o)$, then we call such a parametric location has some relationship with $G^1_k$. We denote such parametric locations of all BRDF samples by $(\theta _h, \theta _d, \phi _d)_{k}$, for which the $j^{th}$ parametric location is $(\theta _h, \theta _d, \phi _d)_{kj}$ and the corresponding measured BRDF is $(\rho ^*_C)_{kj}$ where C stands for R, G or B. It is important to emphasize again that the RGB channels of our models is dependent, so $G^1_k$ is the same for each color channel. Because the structure of our model is the same for each color channel, we only take the red channel as an example.

For convenience, we will denote $G^1_k$ by x, $s_R^1\frac{[D^1(\theta _h)]_{kj}[F_R(\theta _d)]_{kj}}{[\cos {\theta _i}\cos {\theta _o}]_{kj}}$ by $y_j$ and

$$ (\rho ^*_R)_{kj}-d_R-\sum _{m=2}^{m<=3}{s_R^m\frac{[D^m(\theta _h)]_{kj}[F_R(\theta _d)]_{kj}[G^m(\theta _i, \phi _i)]_{kj}[G^m(\theta _o, \phi _o)]_{kj}}{[\cos {\theta _i}\cos {\theta _o}]_{kj}}} $$

by $z_j$.

We need to consider three situations. To distinguish these three situations, we use three different subscripts, r, t and s respectively, to substitude j. The first situation is that $[G^1(\theta _i)]_{kr}$ is dependent on $G^1_k$ but $[G^1(\theta _o)]_{kr}$ is not. So $[G^1(\theta _o)]_{kr}$ is a simple constant and we denote it by $G_o^r$. Since all the other elements of $G^1$ are held constant, $[G^1(\theta _i)]_{kr} = a_r^i + c_r^ix$, where $a_r^i$ and $c_r^i$ are also constants. Then the fitting objective for this situation is

$$\begin{aligned} E_r(x) = w_r[z_r - y_r(a_r^i + c_r^ix)G_o^r]^2 . \end{aligned}$$

(5)

The second situation is reverse. So its fitting objective is similar:

$$\begin{aligned} E_t(x) = w_t[z_t - y_t(a_t^o + c_t^ox)G_i^t]^2 . \end{aligned}$$

(6)

Finally, the third situation is that both $[G^1(\theta _i)]_{kr}$ and $[G^1(\theta _o)]_{kr}$ are dependent on $G^1_k$. The corresponding fitting objective is

$$\begin{aligned} E_s(x) = w_s[z_s - y_s(a_s^i + c_s^ix)(a_s^o + c_s^ox)]^2 . \end{aligned}$$

(7)

Then the total fitting objective correspond to Eq. 2 is $E_R(x) = E_r(x) + E_t(x) + E_s(x)$. Unsurprisingly, its first derivative is just a cubic function, that is, $E_R^{'}(x) = ax^3 + bx^2 + cx + d$. In practice, we use the GSL library to solve the cubic equation $ax^3 + bx^2 + cx + d = 0$. In most case, there is only one real root. If there are three real roots then we need to check which is the best one. Before updating the $G^1_k$, it needs to be checked whether it actually reduces the fitting objective. Therefore, it is guaranteed to reduce the fitting objective after updating all the elements of $G^1$ with such a method.

2.4 The AWLS Method for the Anisotropic Materials

As for the anisotropic materials, using the AWLS method to solve for the D and F factor is similar but there are some problems needed to be addressed. The UTIA measured data is parameterized by $\{\theta _i, \phi _i, \theta _o, \phi _o\}$. Before evaluating the D and F factor, we need to transform the parameter space to $\{\theta _h, \theta _d, \phi _d\}$. We simply use linear interpolation to evaluate the D and F factor. So each parametric location has relationship with 4 control points of the D factor and 2 control points of the F factor. $\{\theta _h, \theta _d, \phi _d\}_{kj}$ represents the $j^{th}$ parametric locations of all BRDF samples that have relationship with the $k^{th}$ element of $D^1$. In the following, we take updating the $k^{th}$ element of $D^1$ as an example. Updating the elements of the F factor is totally the same. For conciseness, we denote the $k^{th}$ element of $D^1$ by x,

$$ {s_R^1(\frac{[F_R(\theta _d)]_{kj}[G^1(\theta _i, \phi _i)]_{kj}[G^1(\theta _o, \phi _o)]_{kj}}{[\cos {\theta _i}\cos {\theta _o}]_{kj}})} $$

by $y_j$ and

$$ (\rho _R^*)_{kj}-d_R-\sum _{m=2}^{m<=3}{s_R^m\frac{[D^m(\theta _h, \phi _h)]_{kj}[F_R(\theta _d)]_{kj}[G^m(\theta _i, \phi _i)]_{kj}[G^m(\theta _o, \phi _o)]_{kj}}{[\cos {\theta _i}\cos {\theta _o}]_{kj}}} $$

by $z_j$.

So the fitting objective for the red channel is $E_R(x)=\sum _jw_j(z_j-(a_j+c_jx)y_j)^2$, which is a simple quatratic function, whose first derivative also can be easily calculated. Finally, because of the parameter space used in the UTIA database, it is easier to update the elements of the G factor but it is still needed to consider the three situations mentioned in the last subsection. There is no need to repeat those fomulations.

2.5 Importance Sampling

There is no satisfying sampling method if we directly use measured data. The default sampling method in pbrt-v2 samples the unit hemisphere with a cosine-weighted distribution. Instead, Bagher et al. [1] used a kind of sampling method, which needed to construct a 2D CDF for each view-slice. It is time-consuming unless we calculate all the CDFs beforehand and store them, which obviously costs too much memory (at least 4G). We use such a method to render noise-free images for comparing, ignoring that it costs too much memory space.

Fortunately, after fitting the measured data to our model, we can use the D factor for importance sampling as in the case of the analytical models. The Independent-G model is unable to do so because the D factors are different for three color channels but our model shares the same D factor instead. For each material in the scenes, it is only needed to calculate $D(\theta _h)\cos {\theta _h}$ and normalize it before rendering. That is, $D(\theta _h)\cos {\theta _h}$ is multiplied by a constant c in order to make sure that

$$\begin{aligned} \frac{\pi ^2}{90}\sum _{i=1}^{90}{c[D(\theta _h)\cos {\theta _h}\sin {\theta _h}]_i} = 1 . \end{aligned}$$

(8)

The normalization formalation 8 can be derived by starting from

$$\begin{aligned} \int _0^{2\pi }{\int _0^{\frac{\pi }{2}}{D(\theta _h)\cos {\theta _h}\sin {\theta _h}}d\theta _hd\phi _h} = 1 . \end{aligned}$$

(9)

Then the resulting $cD(\theta _h)\cos {\theta _h}$ is the pdf (distributed dense function) to draw sample from. We use the Distribution1D class to implement this importance sampling method in pbrt-v3. The pdf $cD(\theta _h)\cos {\theta _h}$ is calculated before rendering and is stored in memory in order to avoid to be calculated whenever it is needed to sample an incoming direction given the outgoing direction. Because our model has three specular lobes and every D factor has 90 elements, 270 extra floating-point numbers are needed to help with the importance sampling.

Using our method, it is only needed to render red-specular-metallic or black-obsidian with 128 samples per pixel to get a good result. The rendering results are shown in Fig. 2.

Fig. 2.

For each material, the most left is the reference image and to the right are one image rendered with the default method with 4 K samples per pixel and two images rendered with our importance sampling method with 64 and 128 samples per pixel respectively. For clarity, only the lower left corners of the spheres are shown (Color figure online).

3 Experimental Results and Analysis

3.1 Result for the Isotropic Materials

In order to compare our model with the Independent-G model, both models are used to reconstruct the reflection data following the format of the MERL database for all the materials, which are then used to render images with 4096 samples per pixel. We also use the HDR environment map (EM) [18] in rendering. In the experiments, two scenes are rendered. One is a sphere in the Grace EM and the other one is the Stanford Buddha in the Grove EM.

Figure 3 shows the average rendering PSNRs of the two scenes for all the materials in the MERL database. We sort the materials by increasing rendering PSNRs of the Independent-G model. It is showed that our model is at least as good as the Independent-G model and in up to 88 cases our model is better than it. It should be noticed that the Bivariate model uses 10$\times$ parameters so in most cases it gets the highest rendering PSNR among all the models showed in Fig. 3. In terms of the rendering PSNRs, our model is closer to the Bivariate model than the Independent-G model and for some materials, our model are even better than the Bivariate model.

We show two of the best rendering results and two of the worst rendering results of our model in Fig. 3. It is hard to see the differences so the PSNRs are shown in the lower left corner of each image. From [1] we know that there is obvious visual error in their fit for silver-metallic-paint but Fig. 3 shows that our model is much more better than the Independent-G model. However, for some other materials like violet-acrylic or nylon, there is no much improvement.

Fig. 3.

Average rendering PSNRs for the MERL isotropic materials. (a) and (b) are two of the best rendering results of our model. (c) and (d) are two of the worst rendering results of our model (Color figure online).

3.2 Result for the Anisotropic Materials

Just like last section, we render a sphere in the Grove EM and the woven cloth using our anisotropic model and the original measured data. The average rendering PSNR of our model is 106.80 and only for 8 materials the rendering PSNRs are lower than 100. We also calculate the reconstruction fidelity, which is also evaluated by the standard PSNR but is calculated between the original measured data and the one reconstructed by our model. We call such a PSNR numerical PSNR. Both the rendering PSNR and the numerical PSNR are shown in Fig. 4. We can compare the numerial PSNR with Fig. 6 in [15]. Although such a compare seems to be meaningless because the purpose of Filip et al. [15] is just to find a fixed number (e.g. 1438) of the most important parametric locations to reconstruct the measured data instead of using any model to fit the measured data. By directly comparing the results numerically it is convinced that our model is good enough even for the anisotropic materials.

Fig. 4.

The rendering PSNRs and the numerical PSNRs for the UTIA anisotropic materials. (a) is one of the worst rendering results. (b) and (c) are two of the best rendering results. For each pair of images, the image rendered with our model is on the left and the reference image is on the right. There is little visual error even for the worst result.

3.3 Data Analysis

We show the resulting factors for some materials in Fig. 5. Obviously, the curves of the D factor for aluminium and violet-acrylic are similar to the curves of the analytical models. This is the reason why the analytical models, e.g. Löw’s ABC model [12], fit aluminium and violet-acrylic as well as the Independent-G model (see the results in the supplementary material of Bagher et al. [1]). But the curve of the D factor for polyurethane-foam is totally different. So analytical models do not fit it well. Moreover, the curves of the F factor are similar for each color channel. It is convinced that using the same F factor for each color channel leads to little loss. For some specular materials, e.g., aluminium and violet-acrylic, there are high-frequency oscillations in the G factor if we do not use the Gaussian smoothing filter after updating the G factor. Although using the Gaussian smoothing filter results in a higher fitting objective, the rendering PSNR for aluminium or violet-acrylic is almost the same. Additionally, it is not necessary to use the Gaussian smoothing filter for the diffuse materials because there are few high-frequency oscillations in the G factor. Therefore, we only use the Gaussian smoothing filter for the specular materials.

Fig. 5.

The D, F and G factors for aluminium, violet-acrylic and polyurethane-foam respectively.

4 Limitation and Future Work

Both our model and the Independent-G model simply use the Lambertian model as the diffuse lobe. Actually, few materials in the MERL database show the Lambertian respond. So a more complicated diffuse component used in [19] is considered to be added in our model. Such a diffuse component considers the remaining energy, after substracting the total energy from the energy for the specular reflection. However, a new fitting method is needed for such a complicated model.

Furthermore, we only use a simple shadowing-masking function $G_2(\theta _i, \theta _o)=G(\theta _i)G(\theta _o)$, which completely ignores the correlation between the masking and shadowing. As Heitz et al. [20] indicates, since some correlations always exist, such a $G_2$ always overestimates the shadowing. Although we have the other three more complicated forms for $G_2$ (see Heitz et al. [20]), the fitting process is troublesome if we continue to use the AWLS fitting method. Again, a new fitting method is needed to help us deal with this complicated model.

5 Conclusion

By using triple specular lobes and sharing factors among three color channels, we design a new model with 822 parameters to fit the MERL isotropic data and with 2010 parameters to fit the UTIA anisotropic data. We do not change the weighting scheme used in Bagher et al. [1] but we do extend their AWLS method to solve for the G factor, which makes updating the elements of the G factor reduce the fitting objective steadily. We offer a simple method for importance sampling which dramatically reduces the number of samples per pixel needed by the specular materials. Our model is better than the Independent-G model for most materials and even better than the Bivariate model for some materials. Note that the Bivariate model uses 10 times number of parameters our model uses. Finally, our model for the UTIA data is also good. The rendering PSNRs for almost all the materials are above 100 and there is no noticeable visual error.