BSDF Importance Baking: A Lightweight Neural Solution to Importance Sampling General Parametric BSDFs

Parametric BSDFs are explicitly controlled by material property parameters, such as roughness, refractive index, attenuation coefficient, and etc. For example, single-bounce Smith microfacet BRDFs can be analytically written as:

where $F$ is the Fresnel term, often approximated with a Schlick reflectance $R_{0}$ and an incident direction ${\bm{\omega}}_{i}$. $G$ is the shadowing-masking term accounting for self-occlusion, and $D$ is the normal distribution function (NDF), parameterized on the half vector ${\bm{\omega}}_{h}=\frac{{\bm{\omega}}_{i}+{\bm{\omega}}_{o}}{|{\bm{\omega}}_{i}+{\bm{\omega}}_{o}|}$, as well as roughness $\alpha_{x}$ and $\alpha_{y}$. Given these material parameters and an outgoing direction $\omega_{o}$, 2D BSDF slices can be easily generated where each pixel records the BSDF radiance value of an incident direction $\omega_{i}$.

For multiple-bounce Smith microfacet BSDFs where the light may bounce multiple times before exiting the surface as shown in Fig. 2, there is still no analytical expression. Both the BSDF evaluation and importance sampling need heavy random walk simulations. In this work, we mainly focus on multiple-bounce BSDFs rather than single-bounce BSDFs, since the former are more challenging. Without loss of generality, we focus on GGX NDFs and describe conductors in the main text.

We also study layered BSDFs (Guo et al., 2018), as illustrated in Fig. 2. The configuration for layered BRDFs includes a top layer using a rough dielectric (with roughness $\alpha$ and refractive index $\eta$), a bottom layer using a diffuse BRDF and a homogeneous participating media (with attenuation coefficient $\sigma_{T}$ and albedo $(\mathbf{a}_{R},\mathbf{a}_{G},\mathbf{a}_{B})$) in the middle. Also note that there are no constraints about the number of layers or the types of BSDFs in each layer in our method. We merely initiate from a simple layered BSDF model to prove the generalization of our method.

Additionally, we study the Disney principled BSDF  (Burley, 2012; McAuley et al., 2012). This artists-friendly model was summarized to be few zero to one intuitive parameters, such as metallic, roughness, specular, anisotropic, sheen, and etc. In this paper, we mainly focus on manipulating metallic $\mathbf{m}$ and roughness $\alpha$ values also shown in Fig. 2.

The MC method provides an unbiased estimator to any definite integral in the domain $\Omega$ with an arbitrary integrand $f(x)$

The MC estimator approximates the integral by drawing samples according to a probability density function (PDF) $p(x)$. The $f(X_{k})/p(X_{k})$ is known as the sampling weight. To minimize the estimation variance, importance sampling is desired – the closer between the “shapes” of $f(x)$ and $p(x)$, the lower variance does the MC estimator have.

MC integration is the core of the modern rendering pipeline, solving the rendering equation at each shading point:

Given the BSDF parameters $\epsilon$ and the viewing direction ${\bm{\omega}}_{o}$, the BSDF becomes a 2D slice, and BSDF importance sampling seeks on a good PDF similar to this 2D BSDF slice, together with a sampling technique that generates samples according to this PDF.

Though BSDF importance sampling can be easily calculated if the whole BSDF values at all dimensions are known, it would be extraordinarily difficult for parametric BSDFs because of storage constraints. At the same time, the efficiency of BSDF importance sampling is crucial, since it will be evaluated at every ray bounce for every sample ray from the camera. Massive amount of heavy BSDF importance sampling calculation will dramatically slow down the rendering calculation. Therefore, we propose to separate BSDF importance sampling calculation into the precomputation and compression of sampling data, as well as efficient runtime value lookup.

Aside from the BSDF parameters $\epsilon$ and the outgoing direction ${\bm{\omega}}_{o}$, our importance sampling network takes two random numbers and outputs the sampled incident direction together with its sampling weight:

where sw is the 3-channel sampling weight. The network structure is shown in Fig. 8. Note that since we have full control of the sampling process, the training data, especially the sampling weight to our BSDF importance sampling network, is guaranteed to be correct (unbiased), albeit some noise can remain. Therefore, the source of bias in our results can only originate from the learning process of the network itself, as will be analyzed in the next subsection.

To validate the functionality of our BSDF sampling network, we compare the importance maps and the binning results between the ground truth and our network’s prediction, as shown in Fig. 9. Thanks to our use of optimal transport that results in smooth importance maps, our network is able to learn to compress them well.

So far, our solution to general parametric BSDF importance sampling has already completed. However, modern rendering pipeline usually supports multiple importance sampling (MIS). In the MIS framework, a renderer should implement not only BSDF sampling, but also BSDF evaluation (returning the BSDF value given the incident and outgoing directions) and PDF query (returning the PDF value similarly). For some BSDFs, these two tasks can be elegantly computed in closed form, e.g., the Disney Principled BSDFs. But for some other BSDFs, it is still difficult to evaluate the BSDF values, e.g., the layered BSDFs.

Therefore, for completeness, we also compress the BSDFs and PDFs, using separate neural networks.

For BSDF evaluation, our evaluation network $\mathcal{E}$ takes an additional incident direction ${\bm{\omega}}_{i}$ and outputs the BSDF value as a 3-channel RGB value:

Note again that the BSDF evaluation network is optional and should only be used when there is no analytical BSDF evaluation scheme. We choose the simple network structure because it fulfills our need already, producing correct BSDF values, as validated in Fig. 1 and Fig. 15. We are aware of other existing neural BSDF/BTF compression methods (Rainer et al., 2019, 2020; Sztrajman et al., 2021; Fan et al., 2022; Zheng et al., 2021). And we also believe that their success can further strengthen our importance sampling scheme when combined together, replacing our simple evaluation network. But we do not extend analysis on them any further, since BSDF evaluation is not our main contribution.

Our PDF query network $\mathcal{P}$ has a very similar definition compared to the BSDF evaluation network. It also takes in the combinations of BSDF parameters as well as the incident directions and returns the PDF value of sampling that direction under solid angle measurement:

For readers unfamiliar with rendering, we make a special note here that our PDF query networks only serve as providing the PDF values for the computation of MIS weights (not sampling weights). A physically correct PDF is certainly welcomed – as it will reduce the variance of the MIS combined results. However, even completely wrong PDFs in this step only leads to higher variance rather than any bias. We elaborate on this more in the next section with a simple but convincing experiment.

Our BSDF evaluation and PDF query networks share a similar lightweight structure, as illustrated in Fig. 10. Note specifically that the emphasis on lightweight neural network design in core rendering (Zhu et al., 2022; Fan et al., 2022) is different from that in deep learning. We use neural networks only as a general tool for efficiently compressing and querying high-dimensional data.

Bias is not desired, but there is no guarantee that neural networks produce fully unbiased results. In our case, the BSDF sampling network learns the sampling weight values, which is equivalent to $\frac{f_{s}\langle{\bm{\omega}}_{i},{\mathbf{n}}\rangle}{\text{PDF}}$. However, this PDF is the correct PDF (rather than the output of our PDF query network) as part of the unbiased training data. And it will not be exposed to any other stages but only together as the sampling weight. Therefore, the only source of bias originates from the learning process of the network itself, and is on sampling weight only, which is usually low frequency and easy to learn as shown in Fig. 7.

For our PDF query network, we do not include any sum-to-one (and other normalization) constraints. As mentioned earlier, this is because the PDF used in MIS can be any value, and the resulting MIS weights may even be negative (Kondapaneni et al., 2019). To further prove this statement, we show the MIS rendering results using a ball status scene with negative uniform MIS PDF, tenth power PDF, and cos-weighted PDF as reference in Fig. 11.

Our neural BSDF importance sampling solution has significant differences from other neural BSDF methods. For example, Fan et al. (2022) sample different layered BRDF parameters by learning a Gaussian lobe and a Lambertian lobe. Similarly, Sztrajman et al. (2021) apply a shallow network to learn the mapping between autoencoded BRDFs and parameters of fitted analytical Blinn-Phong importance sampling parameters. Though Sztrajman et al. (2021) also use a neural network to calculate BRDF importance sampling, their network only provides fitting parameters to the analytical Blinn-Phong model.

Therefore, the core of these importance sampling methods is still the analytical model fitting of BSDF slices. The insight of these strategies relies on a rough prior knowledge of the shape of the BSDF slices. When the prior breaks, the efficiency of sampling will be abysmal. For example, the Blinn-Phong importance sampling model serves as a good sample function for BSDFs with isotropic single specular lobe. However, anisotropic materials or complex multiple-lobe materials cannot be accurately described with only one specular highlight. We show the Blinn-Phong importance sampling fitting results of anisotropic multiple-bounce microfacet BSDF in Fig. 12 for validation.

We have differentiated our approach against normalizing flow-based importance sampling methods at the beginning of Sec. 4.2, and more comparisons can be found in Fig. 13. For other sampling methods, such as Dupuy et al. (2018), they perform targeted design for specific kind of BSDFs. We consider them too far from our general parametric BSDF sampling, and do not extend the further discussion.

We use an AMD 64-core 3995WX CPU machine for training data generation. The material parameters $\epsilon$ are sampled uniformly in their own space except the roughness is sampled in squared space. For outgoing directions, we sample them uniformly in square root of the polar angle to emphasize the grazing angle cases. Each pair of $(\epsilon,{\bm{\omega}}_{ox},{\bm{\omega}}_{oy})$ will be used for generating a 2D BSDF slice with a resolution of 128$\times$128, where each pixel represents an incoming direction $({\bm{\omega}}_{ix},{\bm{\omega}}_{iy})$ and stores the three-channel $f_{s}\langle{\bm{\omega}}_{i},{\mathbf{n}}\rangle$ value. In our experiments, we generate 32,768 slices for our training on each material.

The dataset for all the networks are generated from these 2D BSDF slices. These 2D BSDF slices can be directly used as training data for BSDF evaluation networks. For the PDF query network’s dataset, we calculate the luminance of three-channel BSDFs to get single-channel data and then get the PDF by normalizing the luminance with the summed solid angle measurement. Finally, for the BSDF sampling network, we reparameterize the PDF from the hemisphere to a unit disk and then perform marginalized inverse transform sampling on the unit disk. This ensures that we are sampling on a correct distribution domain. The importance sampled sample points function as the target distributions of OT. Additionally, we also store the three-channel $\frac{f_{s}\langle{\bm{\omega}}_{i},{\mathbf{n}}\rangle}{\text{PDF}}$ values as the sampling weight. Note that the sampling and PDF calculation are all based on solid angle-measured BSDF radiance values, and there is no parameter space transformation. Therefore, there is no need for Jacobian calculation. In summary, all training data originate from BSDF evaluation.

In the next step, we use Geomloss (Feydy et al., 2019) to provide an initial mapping state, and then we use SOT (Paulin et al., 2020) for further optimization. We choose to use two OT solutions because combining two methods will combine the advantages and avoid the limitations of these two methods. SOT computes swiftly, but we notice apparent crevices on BSDF slices even after tens of thousands of iterations, especially for small roughness situations. The networks can grasp these crevices and produce incorrect dark areas in the rendering results. Geomloss runs slower than SOT, but it provides an outstanding initial state for SOT to optimize further. Note again, Geomloss and SOT are just ways of giving the mapping for the next learning step of our method. Other specific implementations of OT are all available choices.

Input encoding has been proved to provide sharper and more accurate results in several works (Mildenhall et al., 2020; Müller et al., 2021). Therefore, we adopt the input encoding for our query parameters similar to (Müller et al., 2021). The frequency encoding consists of $2k$ sine and cosine functions, with $k$ different frequencies from $2^{0}$ to $2^{k-1}$, and the one blob encoding (Müller et al., 2019) has $l$ bins for the discrete Gaussian values of the encoded numbers. Details of each encoding for different types of materials are shown in Table 1.

For BSDF evaluation networks, in order to capture both highlight area’s intensity and non-highlight area’s color, we use symmetric mean absolute percentage error (SMAPE) on BSDF values to optimize the networks as following:

where $\varepsilon$ is set to $0.01$ and $\text{sg}(\cdot)$ means there is no gradient back propagation here.

For PDF query networks, the predicted values are transformed into log space and then we apply $\ell_{1}$ loss as follows:

For BSDF sampling network, we directly apply $\ell_{1}$ loss on network output, i.e. the incident direction. Additionally, to reduce the number of queries for networks in rendering, our network also returns the sampling weight. It is a low-frequency function (its visualization is shown in Fig. 7) that is averaged from the corresponding three-channel BSDFs. We use $\ell_{1}$ loss for the sampling weight outputs in $\log$ space. Thus, the full loss of our sample network shows below:

where $\hat{\omega}$ indicates the predicted incident direction, and the $\hat{\text{sw}}$ is predicted sampling weight which is a 3 channel value. $\omega_{\text{gt}}$ and $\text{sw}_{\text{gt}}$ are the corresponding ground truth of them. We set $\lambda$ to $0.4$ in our training.

We use the ADAM (Kingma and Ba, 2014) optimizer with default parameters and cosine annealing scheduler to optimize our BSDF sampling networks, and Ranger (Wright, 2019) optimizer to optimize the BSDF evaluation and PDF query networks. The batch size of training is set to $1048576$ and learning rate is set to $0.0001$. Each networks are trained with $500$ epochs in total. It takes about 48 hours to train BSDF sampling networks on each material, and 12 hours to train BSDF evaluation networks and PDF query networks separately on each material.

In Fig. 13, we show an anisotropic conductor vase with $\alpha_{x}=0.3$ and $\alpha_{y}=0.1$ under an environment light, considering direct lighting only. We use our BSDF sampling network trained on the multiple-bounce dataset and compare it against Heitz et al. (2016) and Xie et al.  (2019). For fairness, we set the hidden dimension of the RealNVP network by Xie et al.  (2019) as 32, keeping a similar parameter size as ours, and train their model at the same time as ours. Since their model only provides the sampled direction and its PDF, we call the GT evaluation function and divide it by their PDF as the sampling weights. For all the methods, we perform BSDF sampling only. Our method use less time to achieve equal quality than Heitz et al. (2016). At the same time, the result by Xie et al.  (2019) has higher variance with $2\times$ time cost.

BRDF sampling is also tricky for multiple-bounce microfacet conductors with low roughness, due to the high frequency of the BRDF lobes. To show the effectiveness of our BSDF sampling network on such cases, we compare our method against Heitz et al. (2016) on the teapot scene (with $\alpha=0.01$) with equal-spp in Fig. 14. In this comparison, we perform BSDF sampling only. Our result preserves the correct highlights and shows much less variance than Heitz et al. (2016), with only $0.5\times$ extra cost.

Besides the BSDF sampling network, we also introduce BSDF evaluation and PDF query networks, which are essential for complex BSDFs with heavy evaluations, like the multiple-bounce Smith microfacet BSDFs. Our full solution enables MIS for Monte Carlo rendering. In Fig. 15, we compare our methods (sample network only and full solution), and Heitz et al. (2016). For our method (sample network only), we use the BSDF evaluation and PDF from Heitz et al. (2016) to enable MIS. By comparison, we find that our full solution outperforms the others in terms of rendering quality, with only a slight time overhead.

In Fig. 1 (top), we compare our full solution (three networks) with Heitz et al. (2016) on a variety of objects with different colors and roughness lit by an area light. The roughness of the four objects (from left to right) are set as $\alpha=0.2$, $\alpha=0.6$, $\alpha=0.03$ and $\alpha_{x}=0.08$ and $\alpha_{y}=0.3$, covering from low roughness to high roughness, from isotropic to anisotropic materials. By comparison, our result shows higher quality than Heitz et al. (2016) both visually and quantitatively. However, the time cost of our method is about 2$\times$ slower, due to the network inference. Since our current network inference is a simple CPU implementation, we believe that further optimization can decrease our time cost significantly.

In Fig. 16, we compare our BSDF sampling model (trained on layered materials) with Guo et al. (2018) with equal spp on a Shoe scene lit by an environment map. The shoes have a dielectric material with roughness $\alpha=0.08$ and refractive index $\eta=1.5$ as the top layer, a diffuse with color as the substrate, together with a medium ($\sigma_{T}=0.8$). Both methods perform BSDF sampling only, without MIS. With an equal sampling rate, the result by Guo et al. (2018) is much noisier while having a much longer time than ours. The expensive time cost of Guo et al. (2018) is due to the random walk in the dense medium.

Then, we validate our full solution (three networks) on layered BSDFs, by comparing against Guo et al. (2018) with equal time in Fig. 1 (bottom). In this scene, we show several objects, with different roughness $\alpha$, refractive index $\eta$, attenuation coefficient $\sigma_{T}$, and albedo color lit by an environment map, considering both direct lighting and indirect lighting. Both methods use MIS to render this scene. To achieve equal time, we use 1024 spp for our method and 850 spp for Guo et al. (2018). Our result has a lower variance and higher computational efficiency than theirs.