Online Neural Path Guiding with Normalized Anisotropic Spherical Gaussians

2.1 Physically based Rendering and Path Guiding

The radiance of a shading point can be calculated using the rendering equation [Kajiya 1986]: (1) \(\begin{equation} L(\mathbf {x},\omega _o) = L_e(\mathbf {x}) + \int _{S^2}f_s(\mathbf {x},\omega _i, \omega _o)L_i(\mathbf {x}, \omega _i) |\cos \theta _i| d\omega _i, \end{equation}\) which consists of the emitted radiance \(L_e(\mathbf {x})\) at the point and the integral over the sphere \(S^2\), which aggregates the contributions of the incident radiance \(L_i(\mathbf {x}, \omega _i)\) from all directions \(\omega _i\). For each \(\omega _i\), the product of incident radiance, the BSDF \(f_s(\mathbf {x},\omega _i, \omega _o)\), and the geometry term \(|\cos \theta _i|\) signify its contribution to the outgoing radiance. In 3D scenes, due to light bouncing off surfaces, the incoming radiance \(L_i\) can be considered as originating from another point on a different surface. This implies that \(L_i\) also follows the integral form described by the rendering Equation (1) and needs to be evaluated at the point from which the scattered light originates. This makes the problem highly complex, and no analytical solution can be formed in general.

Path tracers solve the rendering equation via the Monte Carlo method, which draws random direction samples \(\omega _i^{\prime }\), and the average of samples is the estimation of the integral, (2) \(\begin{equation} L(\mathbf {x},\omega _o) \approx L_e(\mathbf {x}) + \frac{1}{N}\sum _{k}^{N}\frac{f_s(\mathbf {x},\omega _{ik}^{\prime },\omega _o)L_{ik}(\mathbf {x}, \omega _{ik}^{\prime }) |\cos \theta _{ik}|}{p(\omega _{ik}^{\prime })}, \end{equation}\) where \(p(\omega _{ik}^{\prime })\) is the probability density function (PDF) from which the integrator samples at point \(\mathbf {x}\). Due to the complexity of light transport over the 3D space, the variance of samples is high, and many importance sampling methods have been proposed to reduce this variance (e.g., Veach [1998], Hanika et al. [2015], and Conty Estevez and Kulla [2018]). These methods usually focus on importance sampling of single components of the rendering equation. Here, a universal importance sampling technique for the whole product of the rendering equation can help further improve sampling efficiency. Path guiding, which we describe next, is a genre of methods developed along this direction.

Path guiding methods aim to find sufficient approximation of the incident radiance distribution, which serves to achieve better importance sampling that follows the global illumination distribution rather than local BSDF distribution. Jensen [1995] and Lafortune et al. [1995] first proposed to fit the incident light distribution for more efficient indirect lighting sampling. Vorba et al. [2014] fit a Gaussian mixture to model the distribution estimated by an explicit photon tracing pass. Müller et al. (2017, 2019) proposed the Practical Path Guiding (PPG) algorithm to model the distribution using an SD-tree approach. Based on these techniques, full product guiding methods have also been proposed. Herholz et al. [2016] calculated an approximation of full product on top of an earlier work [Vorba et al. 2014], and Diolatzis et al. [2020] calculated an approximation on top of another study [Müller et al. 2017] to achieve product guiding. Both of these methods operate at the cost of higher overhead as the approximated product needs to be calculated multiple times from the corresponding learned incident radiance distribution. The above techniques share the same idea of partitioning the space and progressively learning discrete distributions. Each discrete distribution is shared among points within the spatial partition. A major issue of this approach is parallax error: The above techniques fail to accurately learn the distribution for close-distance incident radiance where the incoming light quickly changes within the same spatial partition. Ruppert et al. [2020] proposed a parallax-aware robust fitting method to address this issue with a discrete approach. We compared their approach with our network-based approach in Section 7.3. Recent research shows that good path guiding requires a proper blending weight between BSDF and product-driven distributions [Müller et al. 2020], and the learned distribution should be variance aware, since the samples are not zero-variance [Rath et al. 2020]. These findings can help to achieve a more robust path guiding framework.

Recently, neural networks have been used to model scattered radiance distribution; both online and offline learning methods have been proposed. Zhu et al. [2021b] used neural networks to estimate a quad-tree representation of incident radiance distribution using nearest photons as input. Currius et al. [2020] used convolutional neural networks to estimate incident radiance represented by spherical Gaussians; this work is for real-time rendering, but the estimated radiance distribution could be used for path guiding. Zhu et al. [2021a] applied an offline-trained neural network to efficiently sample complex scenes with lamps. This technique requires training a U-Net for more than 10 hours for just a single light source, while the estimated distribution is only a \(16\times 16\) 2D map, which is not sufficient for representing general indirect lighting distributions. A much higher resolution is required for robustly guiding over the entire 3D scene; for instance, PPG’s quad-tree approach can represent a resolution of \(2^{16} \times 2^{16}\)). These methods are categorized as offline learning, since they require training a network offline with massive training samples using ground-truth distributions.

Our framework, however, is categorized as online learning, which learns the distribution on the fly, without a ground-truth distribution as reference. Previously, Müller et al. [2020] adopted normalizing flow [Kobyzev et al. 2021] to model the full product of incident radiance distribution. However, with an implicit density model like normalizing flow, each sample/density evaluation requires a full forward pass of multiple neural networks, which introduces heavy computation costs. In a modern path tracer with multiple importance sampling (typically, BSDF sampling and next event estimation (NEE)), this means full forward pass needs to be executed for each surface sample at least two times. Moreover, the training process requires dense usage of differentiable transforms, which makes the training slower than that for regular neural networks. Indeed, Müller et al. [2020] used two GPUs specifically for the normalizing flow’s neural network computation along with the CPU-based path tracing implementation; nevertheless, the sampling speed was still only 1/4 of PPG. Our work, in contrast, proposes the use of an explicitly parameterized density model in closed form that can be learned using a small MLP. The neural network is used to generate the closed-form distribution model rather than actual samples; therefore, we are able to freely generate samples or evaluate the density after a single forward pass. With careful implementation, we show that our method can reach a similar sampling rate to that of PPG, while the result has lower variance. In research that is parallel with our work, Dong et al. [2023] used a small neural network to estimate per-shading-point distribution, which is similar to our approach. However, we also propose the use of NASG, a novel anisotropic model, in place of classic von Mises-Fisher (vMF) distribution. NASG helps to further improve the fitting accuracy and guiding efficiency. We highlight the benefits through experiments in Section 7.

2.2 GPU Path Tracing

GPU path tracing has gained significant attention in recent years due to its high performance, supported by the GPU’s ability to concurrently render many pixels. While many recent research works have focused on constructing real-time path tracers [Lin et al. 2022; Bitterli et al. 2020; Ouyang et al. 2021; Müller et al. 2021; Chaitanya et al. 2017]), GPU-based production renderers leverage GPU to render fully converged, noise-free results faster [Maxon 2023; OTOY 2023]. Compared to CPU architecture, GPU requires a sophisticated programming design for better concurrency, which is needed for overcoming the large memory latency of GPU. One of the general idea of this is ray re-ordering [Meister et al. 2020; Lü et al. 2017]. Laine et al. [2013] proposed a new architecture that splits full path tracing computation into stages (small kernels) and processes threads in the same stage to maximize performance. Recently, Zheng et al. [2022] proposed a more sophisticated framework for GPU path tracing, which leverages runtime compilation for high performance and flexible implementation. While several GPU-based renderer products are available and it has been proven that path guiding can improve sampling efficiency, few products have provided GPU-based path guiding implementation. This is because existing path guiding approaches are designed for CPUs, and specialization for GPUs only emerged recently at the research level [Dittebrandt et al. 2020; Pantaleoni 2020]. Our work provides a low-cost continuous path guiding option for GPU-based path tracers.

2.3 Density Models

Parametric density models have been extensively explored in statistics. Exponential distribution is a representative family, among which Gaussian mixture is the most widely used density model. Vorba et al. [2014] used a 2D Gaussian mixture to model incident radiance distribution; however, since the domain of a 2D Gaussian is the entire 2D plane, when mapping it to unit sphere, it is necessary to discard samples outside the domain, which affects computation efficiency. Dodik et al. [2022] further proposed using 5D Gaussians to model incident radiance distribution over the space. By using a tangent space formulation, it greatly reduced the number of discarded samples.

In the context of physically based rendering, spherical models have been widely leveraged. The spherical Gaussian (SG) is widely used for radiance representation and density modelling (e.g., Wang et al. [2009]). Actually, a normalized SG is equivalent to the vMF distribution in 3D. It has several desirable properties, such as computational efficiency and analytical tractability for integrals. However, its expressiveness is limited due to its isotropic nature, which restricts the shape of the distribution on the sphere. The Kent distribution [1982] presents a possible remedy by generalizing the 3D vMF model to a five-parameter anisotopic spherical exponential model, providing higher expressiveness. However, its application within our framework is impeded by three significant limitations: high computational cost, the absence of a direct sampling method, and issues with numerical precision, particularly when the parameters for controlling concentration are high. These limitations are elaborated upon in Section 4.1.

Other recent works proposed spherical models specifically for graphics. Xu et al. [2013] proposed an anisotropic spherical Gaussian model to achieve a larger variety of shapes; however, no analytic solution was found for its integral, which is problematic for normalization.

Heitz et al. [2016] proposed another anisotripic model with a closed-form expression called linearly transformed cosines (LTC). However, the integral of LTC requires computation of the inverse matrix, and one component requires in total 12 scalars to parameterize a component.

Another family of density models widely adopted in the rendering context is the polynomial family, among which spherical harmonics is the most commonly used model. Polynomial models have been extensively used in graphics [Sloan et al. 2002; Moon et al. 2016]. They can be easily mapped to unit spheres to represent spherical distributions. However, polynomial models are generally limited to capturing low-frequency distributions. To represent high-frequency distributions, a very high degree is required, leading to a significant increase in the number of parameters. Furthermore, there is currently no efficient approach to directly sample a polynomial distribution, although there does exist some relatively expensive importance sampling schemes (e.g., Jarosz et al. [2009]). These were the main challenges in our pilot study that made us give up the idea of adopting them for our path guiding scheme.

Learning-based density models have been drawing more attention recently [Müller et al. 2020; Gilboa et al. 2021]. The model based on normalizing flow can successfully learn complex distributions [Müller et al. 2020] but suffers from the implicitness and heavy computation mentioned earlier. Marginalizable Density Model Approximation (MDMA) [Gilboa et al. 2021] has been proposed as a closed-form learning-based density model, which is essentially a linear combination of multiple sets of 1D distributions. MDMA showed that closed-form normalization is an essential factor in learning density models, and inspired by their work, we propose such a model for our application.

4.1 Background

One of the main requirements for progressive learning of a distribution with neural networks is that models be easily normalizable, implying that a model must have a closed-form integral.

Marginalizable Density Model Approximation. Gilboa et al. [2021] proposed MDMA, which can be effectively optimized via deep learning. A bivariate dimensional distribution can be modeled as (3) \(\begin{equation} D(x,y) = \sum _{i,j}A_{ij}\phi _{1i}(x)\phi _{2j}(y), \end{equation}\) where \(\phi _{1i}\) and \(\phi _{2j}\) are 1D normalized distributions and \(A_{ij}\) are normalized coefficients that sum to 1.

In practice, the 1D distributions in MDMA can be piecewise linear or piecewise quadratic spline models; however, with its limited number of components, MDMA can learn only a coarse distribution. Therefore, it always gives poor results for distributions with high-frequency spots. Despite the limited accuracy, MDMA works well in learning distributions from sparse noisy samples, compared with non-normalizable models (e.g., polynomial model, spherical harmonics). During training, samples with high energy can lower the contribution of other areas, and eventually the training converges. This feature is crucial in learning distributions from sparse training samples. Normalized Gaussian Mixtures. Normalized Gaussian mixtures can address the accuracy limitations of MDMA: (4) \(\begin{equation} D(\boldsymbol {\mathbf {x}}) = \sum _{i}^{N} A_i\frac{G_i(\boldsymbol {\mathbf {x}};\theta _i)}{K_i}, \end{equation}\) where \(G_i(\boldsymbol {\mathbf {x}};\theta _i)\) are Gaussian distributions parameterized by \(\theta _i\), \(K_i\) are normalizing factors, and \(A_i\) are normalized weights of Gaussians such that \(\sum _{i}A_i = 1\). Gaussian mixtures are highly expressive and also easily normalizable: We look into two common Gaussian distributions that can represent the spherical distributions needed for modeling the lighting: spherical Gaussian and anisotropic spherical Gaussian. Spherical Gaussian. SG is a variant of the univariate Gaussian function defined in the spherical domain: (5) \(\begin{equation} G(\mathbf {v};\boldsymbol {\mu },\lambda) = \exp (\lambda (\boldsymbol {\mu } \cdot \mathbf {v}- 1)), \end{equation}\) where \(\mathbf {v}\) is a unit vector specifying the direction, \(\boldsymbol {\mu }\) is the lobe axis, and \(\lambda\) is a parameter that controls the “sharpness” of the distribution. The normalizing term of an SG can be computed by (6) \(\begin{equation} K = \int _{S^2} G(\mathbf {v};\boldsymbol {\mu },\lambda) d\omega = \frac{2\pi }{\lambda }(1-e^{-2\lambda }). \end{equation}\) Normalized SG is equivalent to the vMF distribution in 3D. Since SG is defined in the spherical domain and is suitable for representing all-frequency light distribution at each sample point, mixture models based on SG have been applied in representing the light map used for precomputed radiance transfer (e.g., Wang et al. [2009] and Currius et al. [2020]). However, SG is an isotropic univariate model and thus less expressive; for example, it cannot represent anisotropic distributions.

Anisotropic Spherical Gaussian. To model anisotropic distributions, Xu et al. [2013] proposed Anisotropic Spherical Gaussian (ASG) to model the light map. ASG can be written as (7) \(\begin{equation} G(\mathbf {v};[\boldsymbol {x},\boldsymbol {y},\boldsymbol {z}],[a,b],c) = c\cdot S(\mathbf {v}; \boldsymbol {z}) \cdot e^{(-a (\mathbf {v}\cdot \boldsymbol {x})^2 - b(\mathbf {v}\cdot \boldsymbol {y})^2)}, \end{equation}\) where \(\boldsymbol {z}\), \(\boldsymbol {x}\), and \(\boldsymbol {y}\) are the lobe, tangent, and bi-tangent axes, respectively, and \([\boldsymbol {x}, \boldsymbol {y}, \boldsymbol {z}]\) form an orthonormal frame; a and b are the bandwidths for the \(\boldsymbol {x}\) and \(\boldsymbol {y}\) axes, respectively; and c is the lobe amplitude. Xu et al. [2013] successfully modeled complex lighting conditions and rendered anisotropic metal dishes using ASG. However, ASG does not have a closed-form solution for the integral, which inhibits its usage for our purpose.

Kent Distribution. Another candidate is the Kent distribution [Kent 1982], which is also an anisotropic spherical distribution, (8) \(\begin{equation} K(\mathbf {v}; \kappa , \beta) = \frac{1}{C(\kappa , \beta)} \exp \left(\kappa (\mathbf {v}\cdot \mathbf {z}) + \beta ((\mathbf {v}\cdot \mathbf {x})^2-(\mathbf {v}\cdot \mathbf {y})^2) \right), \end{equation}\) where \(\kappa , \beta\) are parameters that control anisotropic concentration with \(0\lt 2\beta \lt \kappa\) and \([\boldsymbol {x}, \boldsymbol {y}, \boldsymbol {z}]\) form an orthonormal frame in the same way as ASG. The Kent distribution has a closed-form integral, i.e., its normalizing constant, (9) \(\begin{equation} C(\kappa , \beta) = 2\pi \sum _{i=0}^{\infty } \frac{\Gamma (i+0.5)}{\Gamma (i+1)} \beta ^{2i} \left(\frac{\kappa }{2}\right)^{-2i-\frac{1}{2}}I_{2i+\frac{1}{2}}(\kappa), \end{equation}\) where \(\Gamma (\cdot)\) is the gamma function and \(I_n(\kappa)\) is the modified Bessel function of the first kind. However, since this integral entails the evaluation of a nested infinite series (since \(I_n(\kappa)\) is also an infinite series), sufficient precision necessitates the computation of a significant number of terms, thus imposing a substantial computational overhead. In the context of computer science, the precision of the Kent distribution is further constrained by the limitations on floating-point number precision. Specifically, when the value of \(\kappa\) is large (e.g., \(\kappa \gt 20\)), both the non-normalized density and the normalizing constant approach the representational boundary of floating-point numbers, resulting in a lack of precision. This leads to either a biased result or disrupted learning (where an NaN loss is generated and corrupts the learning process). Within the realm of Monte Carlo rendering, an additional limitation emerges, since Kent distribution lacks a direct sampling method, which necessitates the use of a more computationally expensive rejection method.

4.2 Normalized Anisotropic Spherical Gaussian

Inspired by the models in Section 4.1, we believe a spherical anisotropic exponential distribution can help to achieve more accurate results in our framework. While the Kent distribution exists as a potential candidate, its inherent limitations—namely the necessity of approximating its integral, the precision issue, and the absence of a direct sampling algorithm—prompt us to develop our own distribution, which we call the Normalized Anisotropic Spherical Gaussian, (10) \(\begin{equation} \begin{aligned}&G(\mathbf {v};[\boldsymbol {x,y,z}], \lambda , a) \\ &\quad = {\left\lbrace \begin{array}{ll}\exp \left(2 \lambda \left(\frac{\mathbf {v}\cdot \boldsymbol {z}+1}{2}\right)^{1+\frac{a(\mathbf {v}\cdot \boldsymbol {x})^2}{1-\left(\mathbf {v}\cdot \boldsymbol {z}\right)^2}}-2 \lambda \right)\left(\frac{\mathbf {v}\cdot \boldsymbol {z}+1}{2}\right)^{\frac{a(\mathbf {v}\cdot \boldsymbol {x})^2}{1-(\mathbf {v}\cdot \boldsymbol {z})^2}} & \text{ if } \mathbf {v}\ne \pm \boldsymbol {z} \\ 1 & \text{ if } \mathbf {v}= \boldsymbol {z} \\ 0 & \text{ if } \mathbf {v}= -\boldsymbol {z}, \end{array}\right.} \end{aligned} \end{equation}\) where \([{\boldsymbol {x},\boldsymbol {y},\boldsymbol {z}}]\) forms an orthogonal frame as above, \(\lambda\) is sharpness, and a controls the eccentricity. Figure 3 shows examples of the NASG component with different parameters.

Fig. 3.

View Figure

The derivation of NASG proceeds as follows. First, we employ the standard spherical coordinate expression of \(\mathbf {v}\ne \pm \boldsymbol {z}\) in terms of the polar angle \(\theta\) and the azimuthal angle \(\phi\) (with respect to the orthonormal frame \([{\boldsymbol {x},\boldsymbol {y},\boldsymbol {z}}]\)) to obtain the equation (11) \(\begin{equation} \left(\frac{\mathbf {v}\cdot \boldsymbol {z}+1}{2}\right)^{\frac{a(\mathbf {v}\cdot \boldsymbol {x})^2}{1-(\mathbf {v}\cdot \boldsymbol {z})^2}}=\left(\frac{\cos \theta +1}{2}\right)^{a\cos ^2\phi }, \end{equation}\) which introduces the anisotropic nature of the distribution through the exponent \(a\cos ^2\phi\). Next, we effect a change of variables in the integral of SG (cf. Equation (6)) so that the Jacobian in Equation (11) emerges as part of the corresponding Jacobian. In this way, we arrive at NASG, which satisfies all requirements as a model for learning the distribution needed for unbiased rendering.

Different from ASG [Xu et al. 2013], NASG has an analytical closed-form solution for its integral (see Section B), which makes it easy to normalize: (12) \(\begin{equation} K = \int _{S^2} G(\boldsymbol {\mathbf {v}};[\boldsymbol {x,y,z}], \lambda , a) d \omega =\frac{2 \pi \left(1-e^{-2 \lambda }\right)}{\lambda \sqrt {1+a}}. \end{equation}\) NASG is also expressive, since it can model anisotropic distributions, as well as numerically stable to represent high- and low-frequency distributions accurately.¹ It is also compact, because it is parameterized by only seven scalars, making it highly efficient in GPU-based computation due to its low bandwidth requirement. Additionally, the sampling algorithm of NASG is remarkably efficient (see Section C). These characteristics make NASG a highly feasible and practical option for our framework.

5.1 Network Architecture

We intend our neural network to be as simple as possible, considering the performance requirements for practical use. In this work, we use a four-layer MLP, where each layer has 128 units without bias.

Network Inputs. The input to the model consists of the location of the shading point \(\mathbf {x}\), outgoing ray direction \(\omega _o\), and surface normal \(\mathbf {n}\). The location of the shading point is first encoded into a 57-dimensional vector by positional encoding. Thus, the input size is 63 (57 + 3 + 3) but is padded to 64 for hardware acceleration purposes.² We adopt the simple one-blob encoding [Müller et al. 2019] instead of complex encoding mechanics, such as Multiresolution Hash Encoding [Müller et al. 2022], because there was very little improvement in our pilot study to justify the overhead such mechanics introduce. Details of the encoding scheme can be found in Table 1.

Parameterization of NASG. To improve the efficiency of learning, we reduce the number of parameters used to represent the NASG model. First, we introduce \(\theta ,\phi ,\tau\), the Euler angles representing the orientation of the orthogonal basis with respect to the global axes. Accordingly, \(\mathbf {z}\) and \(\mathbf {x}\) of the orthogonal basis can be represented by (13) \(\begin{equation} \mathbf {z}= \begin{pmatrix} \cos \phi \sin \theta \\ \sin \phi \sin \theta \\ \cos \theta \end{pmatrix}, \ \mathbf {x}= \begin{pmatrix} \cos \theta \cos \phi \cos \tau - \sin \phi \sin \tau \\ \cos \theta \sin \phi \cos \tau + \cos \phi \sin \tau \\ -\sin \theta \cos \tau \end{pmatrix}. \end{equation}\) Thus, we represent NASG with the following seven parameters: \(cos\theta\), \(sin\phi\), \(cos\phi\), \(sin\tau\), \(cos\tau\), \(\lambda\), and a.

Network Outputs. The output of our neural network is a d dimensional vector \(\boldsymbol {o}\), where each component corresponds to the NASG parameter at the shading point \(\mathbf {x}\) and a selection probability c. Assuming the number of mixture components is N, \(d = 5 \times N + 2 \times N + N + 1 = 8N + 1\) (see Table 2). These parameters are further decoded to represent the actual parameters as summarized in Table 2. The mixture weights of NASG \(\mathbf {A}\) are normalized with a softmax function.

5.2 Training

We now describe our loss function and the training process. As shown in Figure 4, we utilize automatic differentiation to train the network effectively with sparse online rendering samples. Since the training data are noisy online samples, we need a loss function that is more robust than regular L1 or L2 losses.

Fig. 4.

View Figure

Kullback–Leibler Divergence. The design of a loss function for distribution learning in PBR context has been studied in several previous works [Müller et al. 2019; , 2020; Zhu et al. 2021b; , 2021a]. By using \(q(\omega _i;\gamma)\) to denote the NASG distribution (with NASG parameters \(\gamma\)) estimated by the neural network, to achieve importance sampling of the rendering equation, the optimal distribution should be proportional to the product: (14) \(\begin{equation} q(\omega _i;\gamma) \propto f_s(\mathbf {x},\omega _i, \omega _o)L_i(\mathbf {x}, \omega _i) |\cos \theta _i|. \end{equation}\) Therefore, we use \(p(\omega _i) = F f_s(\mathbf {x},\omega _i, \omega _o)L_i(\mathbf {x}, \omega _i) |\cos \theta _i|\) as our target distribution, where F is a normalizing term whose value is unknown. We can still train the system without knowing F as shown below.

In our framework, we use the Kullback–Leibler (KL) divergence to represent the likelihood between our estimated distribution \(q(\omega _i;\gamma)\) and the target one: (15) \(\begin{equation} D_{KL}(p(\omega _i)||q(\omega _i;\gamma)) = \int _{S^2}p(\omega _i)(\log [p(\omega _i)] - \log [q(\omega _i;\gamma)])d\omega _i. \end{equation}\) Since \(p(\omega _i)\) is \(\gamma\) independent, we get (16) \(\begin{equation} \nabla _\gamma D_{KL}(p(\omega _i)||q(\omega _i;\gamma)) = -\int _{S^2}p(\omega _i)\nabla _{\gamma }\log [q(\omega _i;\gamma)] d\omega _i. \end{equation}\) Then, the integral can be replaced with our one-sample estimation: (17) \(\begin{equation} \nabla _{\gamma } D_{KL}(p(\omega _i)||q(\omega _i;\gamma)) = -\, \mathbb {E} \left[ \frac{p(\omega _i)}{\hat{q}(\omega _i;\gamma)}\nabla _{\gamma }\log [q(\omega _i;\gamma)]\right], \end{equation}\) where \(\hat{q}(\omega _i;\gamma)\) denotes the PDF of the distribution the sample obeys during rendering (see below). Although the right-hand side still involves the global scaling factor F, it is canceled in moment-based optimizers such as Adam [Kingma and Ba 2014]; accordingly, we can train the system even though there is an unknown scaling factor in \(p(\omega _i)\) [Müller et al. 2019].

Selection Probability. Following Müller et al. [2020], we do not directly use a distribution estimated by our neural network. Instead, we regard our scattering sampling process as an MIS process that blends the learned distribution and the BSDF distribution \(p_{f_s}(\omega _i)\). As we also learn a selection probability c, as described in Section 5.1, the MIS PDF at shading point \(\boldsymbol {\mathbf {x}}\) is \(\hat{q}(\omega _i;\gamma) = cq(\omega _i;\gamma) + (1-c)p_{f_s}(\omega _i)\). However, optimizing \(\hat{q}(\omega _i;\gamma)\) naively leads to falling into local optima with degenerate selection probability; therefore, our final loss function blends \(\hat{q}(\omega _i; \gamma)\) with \(q(\omega _i;\gamma)\), (18) \(\begin{equation} loss = e D_{KL}(p(\omega _i)||\hat{q}(\omega _i; \gamma)) + (1-e)D_{KL}(p(\omega _i)||q(\omega _i; \gamma)), \end{equation}\) where e is a fixed fraction that we set to 0.2.

7.1 Variance Reduction

We implemented our framework on top of our in-house GPU production renderer and evaluated the quality of the images by comparing them with those produced by plain path tracing using the same number of samples per pixel (SPP). We implemented both the sampling algorithm and the optimizer on the GPU, leaving the CPU to only perform scheduling tasks. We also compared different guiding schemes, i.e., guiding with incident radiance distribution, cosine-weighted radiance distribution, and full product distribution. In addition, we compared NASG with SG, MDMA, 2D Gaussians, and Kent distribution. Here, we used 8-component NASG (64 parameters). For other density models, we adjusted the number of parameters of each Gaussian to match that of NASG: SG uses 14 components (70 parameters), 2D Gaussians uses 12 components (72 parameters), Kent distribution uses 8 components (64 parameters), and MDMA uses 8 1D distributions for both dimensions (80 parameters). For Kent distribution, even when calculations are performed in the logarithmic space to mitigate the precision issues, we found it necessary to restrict the range of \(\kappa\) to \(0 \lt \kappa \lt 20\) and accordingly constrained \(0 \lt 2\beta \lt \kappa\) to prevent biased results or unstable learning.

We run our GPU path tracer on a conventional PC with a four-core i7 9700K CPU and an Nvidia 2080ti GPU. The implementation of our framework introduces roughly 600 MB memory overhead on the GPU, mainly for buffering batched network execution. Our Pytorch optimizer takes another 600 MB of GPU memory to cache tensors and gradients. The test 3D scenes are shown in Figure 5, and we render the scenes with NEE enabled. All of the rendered images are included in the supplemental material. Table 3 shows the results of the quantitative evaluation.

Fig. 5.

View Figure

Learned distribution. Our framework can easily learn a full product. In all of the test scenes, guiding with full product distribution gives lower variance. Guiding with full product helps to effectively reduce variance when rendering smooth surfaces. This is because learning full product not only improves the accuracy of the learned distribution but also helps to learn an optimal selection probability. The difference is more obvious in scenes such as Ajar and Corridor, where the floors are large surfaces with low roughness. When guiding with incident radiance distribution, the lack of BSDF distribution approximation results in many more outliers, as shown in Figure 6(a).

Fig. 6.

View Figure

Distribution Models. Compared to traditional distribution models such as spherical Gaussians, our NASG model is more expressive and achieves lower variance in most scenes, and this is especially the case for scenes with anisotropic lighting conditions, such as Ajar, Ring, and Glossy Tube. In Figure 6(b) we show an example of a scene where the dominant indirect lighting is anisotropic. Although Kent distribution is also a spherical anisotropic distribution, we have observed two limitations: First, the computationally expensive nature of Kent distribution makes it \(5\times\) slower than NASG. Second, due to the constraints imposed to mitigate the precision issue, Kent distribution can only yield a coarse approximation when the target distribution contains high-energy spots. As illustrated in Figure 6(a), the caustics emanating from the smooth floor are not accurately reconstructed when using Kent distribution. In addition the rendered images, several learned distributions with different models are compared in Figure 7.

Fig. 7.

View Figure

7.2 Performance

An importance sampling method is more practical when it achieves a lower level of error with the same time budget. We conducted an experiment to compare the equal-time performance of our framework with different configurations, as well as with plain path tracing, using the same set of test scenes as in our previous experiments. There are several meta-parameters that can influence performance and accuracy, including the maximum sample size (S), number of hidden units (HU), and number of components (C). In our default configuration, we set \(S = 2^{16}\), HU = 128, and C = 8. We show the results under different configurations in Table 4. In general, the default configuration achieves a good balance between learning accuracy and sampling speed: It achieves the lowest MAPE in four scenes of eight. Even in the scenes where it does not perform the best, its MAPE is close to the best results. With 32 components (\(C=32\)), the result of Box has the lowest MAPE, while the sample count is only 58% of our default configuration. Actually in a simple scene such as Box, where performance is not a bottleneck and all methods achieve significant number of samples, a more precise distribution can be learned with more components, resulting in a lower error. However, in complex scenes where higher variance is introduced by multiple light sources or complex materials, sufficient samples are required for convergence for Monte Carlo rendering. With \(C=32\), sample count is low under limited computation budget, resulting in larger error.

Neural network overhead. In the experiment we measured the time cost of three parts of computation: forward execution for every shading point, training of the neural network, and path tracing. While changing through scenes, we found that typically 15% of the time is used for forward execution, 35% for training, and 50% for path tracing. The training overhead is almost constant among different image resolutions, which makes our framework relatively costly when rendering small images. For example, in the Bidir scene, which uses a small \(512\times 512\) resolution, the ratio of neural network overhead appears to be larger. From our experience, increasing the training steps (i.e., more training batches in each iteration) leads to higher accuracy in the learned distribution but with significant overhead. We conjecture that using meta-learning to accelerate the training process would be a promising direction, and leave this as a future work. Network scale. Intuitively, a larger network, with either more weight parameters or more components, could achieve higher accuracy. However, as neural network computation is the primary overhead in our framework and complex scenes require a significant sample count to converge, it is still worthwhile to trade off accuracy for computation overhead given our limited time budget.

7.3 State-of-the-Art Comparison

To compare our framework with various state-of-the-art path guiding techniques, we implemented it on the Mitsuba renderer. In our Mitsuba implementation, we set \(S=2^{18}\) and \(\nu = 4\) to better show the potential of our framework. The test scenes are shown in Figure 8, and we rendered all of the images in this experiment with two Intel Xeon E5-2680v2 20-thread CPUs (40 threads in total). We compared the results of our framework with the authors’ implementations of PPG [Müller 2019] and PAVMM [Ruppert et al. 2020], along with our implementation of NIS [Müller et al. 2019]. For all techniques, we enabled NEE. To render with PPG, we disabled inverse-variance-based weighting and selection probability optimizaton, while enabling tree splatting using default configuration. For PAVMM, we used BSDF product sampling when available and cosine weight sampling otherwise (the authors’ implementation requires specific scene structure, which conflicts with scene setups), and we set maxSamplesPerLeafNode to 16k, with all other options left as default. We implemented NIS following the original paper; based on our framework, we replaced the sampling-density evaluation component with two-layer quadratic linear normalizing flow, and we added another MLP to learn selection probability. We used the encoded input of our NASG framework as extra features for the transform layer, and the feature dimensions are encoded with one-blob encoding (we used 64 bins). The training sample size and training step’s configuration are aligned with our framework. The major difference of our implementation from that in the original paper is that we only use lightweight MLPs, at the same scale used for our framework, rather than ResNets. The source code of our framework and our NIS implementation are given in the supplemental material.

Fig. 8.

View Figure

The rendering speed of of a CPU-based neural path guiding implementation is strongly affected by the actual implementation, so we mainly compared equal-sample results rather than raw performance.

7.3.1 Variance Reduction of Guiding Method.

To objectively evaluate the efficiency of different guiding algorithms, we first disabled orthogonal features and rendered scenes with our framework and other methods. To ensure a fair comparison, we used the same sampling parameters in all of the methods. The selection probability is set to the fixed fraction 0.5 for all methods. All of the methods have a 128-sample training budget and 384 samples for rendering (512 samples in total). We achieved this by clearing the rendering accumulation after the 128th sample, and we accumulated the results for the rest of the samples.

Figure 1 and Figure 9 reports statistical results and Figure 10 shows the convergence behavior during the experiment. From the results, we can observe that each previous work has both pros and cons; consequently, one technique may perform very well in one kind of scene but poor in others (details discussed below). Our framework, however, while achieving overall low error also provides a more balanced performance in all test scenes.

Fig. 9.

View Figure

Fig. 10.

View Figure

Comparison with PPG. PPG uses quad-tree for directional representation, which can closely model distributions with multiple high-energy points. Hence, in the Torus scene it achieves low variance, while PAVMM struggles due to the limited components used. However, in the Maze scene, the parallax issue makes the PPG result very noisy, as shown in Figure 9(c). Parallax error occurs, because PPG’s spatial partition algorithm has difficulty learning from narrow and close incident radiance. Our framework learns the continuous spatial-varying distribution and does not suffer from such issues.

Comparison with PAVMM. PAVMM performs generally better than PPG. With its parallax-aware fitting method, it learns distributions in the Maze scene most accurately. However, PAVMM produces a rather noisy result in the Ajar scene, as shown in Figure 9(d). This is due to two reasons: First, the major illumination is a line-shaped lit surface on the wall at the right side of the scene (lit by light source outside the door via a small crack), which is anisotropic and difficult for vMF to fit. Second, the fitting algorithm of PAVMM encounters a challenge in modeling the indirect lighting from the specular checkbox floor. In contrast, our proposed framework exhibits greater robustness by leveraging NASG’s ability to accurately capture the anisotropic distribution with a parsimonious set of components, thus achieving low variance.

Comparison with NIS. With the same configuration of training samples and steps, NIS in general achieves similar variance, except in the Maze scene (as shown in Figure 9(c), in which the learning of the 2D distribution suffers from discontinuity at the edge of the azimuth axis). Compared to the results reported in the original NIS paper [Müller et al. 2019], our NIS implementation has a slightly higher variance. We attribute this to two main factors: (a) Our implementation employs a lightweight two-layer MLP-based transform as opposed to the original, relatively more complex, four-layer ResNet-based version, and (b) we limit our training to a subset of samples at each iteration. However, even with this relatively lightweight implementation, NIS takes roughly \(2\times\) longer than our framework: In Table 5, we compare the rendering time of NIS and our NASG framework. The extra overhead is due to the usage of coupled networks, and the normalizing flow needs to be evaluated two times at each non-delta shading point (one for scattering sampling and the other for NEE). The intrinsic characteristics of normalizing flow as an implicit density model, characterized by the utilization of multi-layer transformations, imply that a larger network architecture may be imperative to effectively capture complex distributions with accuracy. Moreover, it is plausible that a greater number of training samples or iterations might be essential to ensure convergence. These adjustments are likely requisite for attaining the results reported in the original paper. However, they also inherently bring a significant increase in computational overhead. In other words, with limited computation resources, our framework can learn more accurate distributions.

Table 5.

Scene	NASG (ours)	NIS
Ajar	188.8 m	520.4 m
Bathroom	56.2 m	166.8 m
Bidir	35.6 mm	89.6 m
Kitchen	50.1 m	163.5 m
Maze	66.4 m	84.42 m
Box	12.4 m	42.7 m
Staircase	46.3 m	217.2 m
Torus	27.3 m	109.2 m
Wine	32.2 m	128.4 m

• NIS uses a normalizing flow to model distribution, which takes roughly twice as much time to calculate compared to our explicit NASG model.

View Table

Table 5. Rendering Time (Minutes) of Proposed NASG Framework and NIS Implementation in Mitsuba

• NIS uses a normalizing flow to model distribution, which takes roughly twice as much time to calculate compared to our explicit NASG model.

8.1 Performance in GPU Path Tracing

A GPU-based brute-force path tracer can be \(20\times\) faster than a CPU-based one. If an importance sampling method provides low per-sample variance yet high overhead, then it cannot defeat a brute-force path tracer on GPU. As a network-based guiding method, the major overhead is the network itself, and, therefore, we must adopt a sampling strategy that can be implemented in a small network. This motivates us to find a solution for learning an explicit distribution rather than an implicit model such as one based on a normalizing-flow. Despite its powerful computation ability, a GPU’s large memory bandwidth can only be leveraged by carefully designed high-concurrency algorithms. This raises the necessity of reducing code branching as much as possible, which conflicts with the nature of path tracing. Otherwise, it suffers from a much more obvious memory latency than CPU, which causes many CPU algorithms to perform relatively slowly on the GPU (e.g., PPG’s quad-tree approach requires multiple random accesses to GPU memory, thus causing long stalls). Recently, with hardware improvement, we see some acceleration hardware that is specifically designed for neural networks, which our framework could potentially benefit from.

8.2 Learning Rate

In our GPU implementation, we used a relatively large learning rate (0.002). This decision was strategically aimed at enabling our model to swiftly navigate the parameter space, thus facilitating rapid convergence. This feature is particularly crucial in leveraging early samples to their fullest potential in our context. However, a large learning rate limits our model’s ability to reach the deepest minima. That being said, at the latter stage of training, the resulting “jitters” remain close to the ground truth, ensuring that our learned distributions are valid approximations. We leave the improvement of optimization strategy as future work.

8.3 Limitations and Future Work

Although in general our configuration for training works well, sometimes we encounter cases where we need to use a smaller learning rate to avoid degenerated distributions. Such degenerated distributions have been observed only a few times when rendering caustics. They could be due to the complexity of the distribution, where only a small amount of light hits the point from a narrow direction. When the early samples based on BSDF fail to produce a representative initial distribution, the system could fail to sample important directions, which could eventually lead to degeneration. Adopting an adaptive learning strategy could possibly improve robustness; we leave this as a future work.

Another interesting extension of our framework could be path guiding in scenes with participating media. This could be achieved in theory (similar to what was done previously [Herholz et al. 2019]), yet the ability to learn 3D distributions remains an area to be further investigated.