A method for estimating the errors in many-light rendering with supersampling

Hirokazu Sakai, Kosuke Nabata, Shinya Yasuaki, and Kei Iwasaki,2 ()

Abstract In many-light rendering, a variety of visual and illumination effects, including anti-aliasing,depth of field, volumetric scattering, and subsurface scattering, are combined to create a number of virtual point lights (VPLs). This is done in order to simplify computation of the resulting illumination. Naive approaches that sum the direct illumination from many VPLs are computationally expensive; scalable methods can be computed more efficiently by clustering VPLs, and then estimating their sum by sampling a small number of VPLs. Although significant speed-up has been achieved using scalable methods, clustering leads to uncontrollable errors, resulting in noise in the rendered images. In this paper, we propose a method to improve the estimation accuracy of manylight rendering involving such visual and illumination effects. We demonstrate that our method can improve the estimation accuracy by a factor of 2.3 over the previous method.

Keywords anti-aliasing; depth of field; many-light rendering; participating media

1 Introduction

Many-light rendering methods simplify the complex computation of global illumination into a simple summation of the contributions from many virtual point lights (VPLs) [1]. In many-light rendering,the incident radiance at a point to be shaded (the shading point) is calculated using VPLs. Because the number of VPLs used is generally quite large,previous methods [2-4] have clustered VPLs to streamline the computation. These methods, however, cannot control the errors produced by clustering, resulting in noise in the rendered images. To eliminate noise in such systems,users must adjust clustering parameters in tedious trial-and-error processes.

To address this problem,Nabata et al. [5]proposed a method to estimate errors in the pixel values due to VPL clustering. This method,however,estimates each pixel value using a single shading point and cannot be applied to visual effects such as anti-aliasing and depth-of-field (DOF), which require multiple shading points to estimate the pixel value. Walter et al. [6]proposed a method called multidimensional lightcuts(MDLC), which attempts to control the clustering errors in many-light rendering of various visual effects.This method, however, does not estimate the error in each pixel value, which appears as noise in the rendered images. This paper proposes a method to improve the accuracy of pixel values for visual effects such as anti-aliasing, DOF, and volumetric effects,as shown in Fig.1. Our method automatically partitions clusters and continues sampling them until the relative errors in the pixel values are smaller than a user-specified threshold.

2 Related work

2.1 Scalable VPL rendering

Recent advances in many-light rendering have demonstrated that global illumination effects can be more than adequately approximated using many virtual lights [1, 7]. Keller [8] introduced the instant radiosity method, which calculates the indirect illumination from virtual point lights(VPLs). Walter et al. [2,6]proposed a scalable solution to many-light methods using a hierarchical representation of VPLs,called cuts. Haˇsan et al. [9] represented many-light rendering via a large matrix, and explored the matrix structure by using row-column sampling.Ou and Pellacini [3] clustered the shading points into groups called slices and performed matrix rowcolumn sampling for each slice. Georgiev et al.[10] proposed an importance sampling method for VPLs by recording the contributions of the VPLs at each cache point. Yoshida et al. [11] proposed an adaptive cache insertion method to improve VPL sampling. Wu and Chuang [12] proposed the VisibilityCluster algorithm, which approximates the visibilities between each cluster of VPLs and those of the shading points by estimating the averaged visibility. Huo et al. [4] proposed a matrix samplingand-recovery method to efficiently gather contributions from the VPLs by sampling a small number of them. Although these methods can significantly accelerate many-light rendering, the results rely on user-specified parameters, and finding the optimal parameter values remains a challenging task.

Fig. 1 Rendering with anti-aliasing (256 spp): (a) reference, (b) our method, (c) multidimensional lightcuts (MDLC). (b) and (c) are rendered so that the relative errors are less than 2%. (d) and (e) indicate relative errors (error color-code shown at far right). Values in (d) and (e) give the percentages of pixels with relative errors of less than 2%. Our method improves the estimation accuracy of pixel values by approximately 50% compared to MDLC. The low estimation accuracy of MDLC results in noise as shown in (c).

2.2 Participating media

VPL rendering suffers from splotches that stem from the singularity of the geometry term relating shading points and VPLs. Clamping is often used to avoid these artifacts. Engelhardt et al. [13] proposed a method to compensate for energy loss due to clamping and proposed a rendering method for participating media using VPLs. Nov′ak et al. [14]proposed virtual ray lights (VRLs) to alleviate this singularity, and several groups have proposed acceleration methods[15, 16] that use VPLs and VRLs. These methods,however, do not estimate the errors due to clustering of the VPLs and VRLs.

2.3 Subsurface scattering

Arbree et al. [17] proposed a scalable rendering method for translucent materials using VPLs. Walter et al. [18] proposed bidirectional lightcuts that supports complex illumination effects, including subsurface scattering. Wu et al. [19] formulated a radiance computation method for subsurface scattering by sampling light-to-surface and surface-tocamera matrices. Although these methods can render translucent materials efficiently by clustering VPLs,they cannot control the errors due to clustering.

To address this problem, we provide an error estimation method that can handle a wide variety of illumination and visual effects. After specifying the relative error tolerance ∈and the probability α, our method stochastically estimates the relative errors with probability α,i.e.,using our method,the relative errors of a proportion α of the pixels in the rendered image is likely to be less than ∈.

3 Background

In many-light rendering, the outgoing radiance L(x,xv) at shading point x toward the viewpoint xvis calculated by

where L is the set of VPLs, I(y) is the intensity of VPL y, and f(y,x,xv) is the material function that encompasses the bidirectional reflectance distribution function (BRDF) on the surface and the phase function within the volume. V(y,x) and G(y,x)are the visibility and geometry terms [1, 7] relating y and x, respectively. To estimate the pixel value for the rendering of the participating media, antialiasing, and DOF, the outgoing radiances from multiple shading points are required, as shown in Fig. 2. For example, when multiple shading points are generated via supersampling for anti-aliasing, the pixel values are calculated using the weighted sum of the outgoing radiances from the shading points. The pixel value I produced by supersampling is calculated as

where G is the set of shading points, and W is the weighting function. The weighting function W and material function f depend on the visual effects and the rendered objects. Details of W and f are provided in Sections 4.4-4.6.

The computational cost of Eq. (2) is proportional to the product of the number of VPLs |L| and the number of shading points|G|. Since,in general,many VPLs are used, the computational cost of Eq. (2) is prohibitive. MDLC is an efficient scalable many-light rendering method that clusters VPLs and shading points [6]. MDLC implicitly represents the hierarchy of clusters with a product graph that consists of pairs of clusters for VPLs and shading points. MDLC estimates the upper bound for clustering error, but does not estimate the error in the pixel value(i.e.,the sum of the errors due to clustering). We improve the estimation accuracy of the pixel value by combining a method for estimating errors [5] with MDLC [6].

Fig. 2 Pixel value estimation using multiple shading points.

4 Proposed method

4.1 Overview

Figure 3 shows an overview of our method, which estimates image pixel values by clustering VPLs and shading points. Clusters of VPLs and shading points are referred to as VPL clusters and shading clusters, respectively. VPL clusters and shading clusters are represented by binary trees, whose leaf nodes correspond to VPLs (or shading points), and whose inner nodes correspond to clusters of VPLs(or shading points). As shown in Fig. 3(b), we refer to the binary trees for the VPL clusters and the shading clusters as the light tree and the shading tree,respectively. Following MDLC, we denote a cluster pair consisting of a VPL cluster CLand a shading cluster CGby (CL,CG).

We create the VPL clusters and build the light tree, used for all pixels, in a preprocess. For each pixel, we first generate shading points and build the shading tree by clustering the shading points. We prepare a priority queue Q that stores the cluster pairs in descending order of standard deviation of each pair. Q is initialized with the pair (CLr,CGr) of root nodes of the light tree and the shading tree. By using this pair, the estimate ^I of the pixel value and the estimated error Δ^I =|^I-I| are initialized, and the standard deviation σ is calculated (see Section 4.2). Then we repeat the following processes:

Fig. 3 Overview of our method, explained for rendering translucent materials. (a) Subsurface scattering of light at xj is calculated by generating shading points xi, which are importance sampled. (b) A VPL cluster CLr and shading cluster CGr are constructed. (c) Two VPLs and shading points from each cluster are sampled from the pair (CLr ,CGr ), then the pixel value ^I, error Δ^I, and standard deviation are estimated.(d) One of the clusters is subdivided if the estimated error Δ^I is greater than tolerance ∈^I. A cluster is selected based on the length of the diagonal of its bounding box weighted by the numbers of VPLs and shading points in each cluster. Here, the VPL cluster CLr is selected and replaced with two child nodes CL1 and CL2; two new pairs (CL1,CGr ) and (CL2,CGr ) are created. The pixel value is estimated from the two pairs (CL1,CGr ) and (CL2,CGr ).

2. Subdivide one of the two clusters,or,and move down one step in the corresponding binary tree. Cluster selection is described in Section 4.3.Replace the selected cluster with two child nodes and create two new pairs.

3. Update ^I and Δ^I, and calculate the standard deviations for the two new pairs. Push the two new pairs into Q.

4. Terminate the process if the estimated error,Δ^I, is smaller than the tolerance, ∈^I, where ∈is the relative error tolerance. Otherwise, return to Step 1.

4.2 Estimating I and ΔI

Outgoing radiances Iidue to the shading cluster CGiilluminated from VPL clusterare calculated by

As summing over all the VPLs and shading points in clusters CLiand CGiis computationally expensive,our method estimates Iiby sampling a small number of VPLs and shading points as follows:

where xkand ykare the k-th sample of the shading point and VPL, respectively, K is the number of samples, and pGand pLare probability mass functions for sampling, calculated as follows:

Substituting these functions into Eq. (4), ^Iiis calculated from:

Computing the error ΔI = |^I -I| necessitates knowing the true value of I in Eq. (2), but the computational cost of obtaining this value is high.Therefore, our method estimates Δ^I using the following equation [5]:

where tαis the α quantile of the t-distribution, N is the number of pairs, and s2iis the sample variance of Iifor the i-th pair of; Eq. (6) is derived in the Appendix.

To select a pair of clusters to be subdivided,the standard deviation σifor each pairis required. Although the sample variance s2ican be used to estimate σi, the accuracy of this estimate is low, since our method estimates ^Iiwith K = 2 samples, as does the previous method [5]. Instead of using the sample variance, our method calculates the standard deviation σifor the i-th pair using the upper bounds of f and G, as follows:

where fuband Gubare the upper bounds of f and G,respectively, within the VPL cluster CLiand shading cluster CGi.is the standard deviation of the visibility function; a maximum value of 0.5 is used. fuband Gubfor CLiand CGiare calculated in a similar way to those in MDLC.

4.3 Cluster selection

We now determine which cluster of the two clusters in the pair with maximum standard deviation is to be subdivided. In MDLC, the cluster is chosen by a refinement heuristic; it basically selects the cluster with the largest bounding box. Although this works well in some cases, it predominantly selects the VPL cluster CL, since VPLs are usually distributed throughout a scene, whereas the shading points are distributed locally, as shown by Fig. 2. In addition,the numbers of VPLs and shading points are often substantially different (e.g., the number of VPLs in Fig. 1 is 35.1k while the number of shading points is only 256). These facts cause unequal subdivisions between VPL clusters and shading clusters, which results in lower estimation accuracy. To address this problem, we propose a cluster selection method that accounts for the size difference between the bounding boxes of the root nodes of the VPL cluster and the shading cluster. We also consider the number of VPLs and shading points when selecting the cluster to be subdivided. Specifically, the bounding box for each shading cluster is scaled using the following two coefficients:

where CLrand CGrare the root nodes of the VPL and shading clusters, respectively, lCLrand lCGrare the diagonal lengths of the bounding boxes for CLrand CGr, respectively, and |L| and |G| are the numbers of VPLs and shading points, respectively.

4.4 Anti-aliasing and DOF

In many-light rendering with anti-aliasing or DOF,Nsppviewing rays are generated via sampling on the screen pixel or camera lens, and the point of intersection between each viewing ray and the surface in the scene is calculated. The material function f at each shading point xiis represented by a BRDF fr.As we use a simple box filter, the weighting function is calculated by W(xi)=1/Nspp; other filters (e.g.,Gaussian filters) can also be applied.

Here, for simplicity, we describe methods for rendering translucent materials and participating media with a single ray per pixel; however, with simple modifications, our method can render these materials using multiple rays, as shown in Fig. 8.

4.5 Rendering translucent materials

The outgoing radiance at xoon the surface of a translucent material in viewing direction ωodue to multiple scattered light within the translucent material is calculated using the diffuse bidirectional scattering surface reflectance distribution function(BSSRDF) Rd[20]: subsurface scattering of light at xjin viewing direction ωjis given by

where Tηis the Fresnel transmittance,xiand niare a point and the normal to the surface of the translucent material at xi, respectively, rij=‖xi-xj‖, A is the set of points on the surface of the translucent material,and ? is a set of directions over the hemisphere.

Our method traces Nspprays through each pixel and defines the point of intersection between the j-th ray and the surface of the translucent material as xj.

The pixel value due to the subsurface scattering of light is calculated by

By comparing Eqs. (11) and (2), the weighting function W and material function f are represented as follows:

4.6 Rendering participating media

Like MDLC, our method assumes homogeneous participating media and an isotropic phase function.To calculate a pixel value in the presence of homogeneous participating media, the scattered radiances are integrated along the viewing ray, as shown in Fig.2(b). Let xobe the point of intersection of the viewing ray and the surface in the scene.Thus, the outgoing radiance from xois calculated using the sum of the reflected radiance Lsat xoand the scattered radiance Lmalong the viewing ray, as follows:

where τ(xo,xv) = exp(-σt‖ xo- xv‖) is the transmittance, and σtand σsare the extinction and scattering coefficients, respectively. x(t) is the point on the viewing ray parameterized by the distance t from viewpoint xv. L(xo,xv) in Eq. (12)is calculated using Eq. (1). The weighting function W and material function f for Lsare represented by W(xo) = τ(xo,xv) and the BRDF fr(y,xo,xv),respectively. The scattered radiance L(x(t),xv) at x(t) is calculated using:

where fpis the (isotropic) phase function.

To compute the integral in Eq.(13),shading points xiare generated by uniformly subdividing the viewing ray with step size Δt; Lmis calculated by summing over the shading points and VPLs as follows:

Comparing this equation with Eq. (2), the material function f is represented by the phase function fp,and the weighting function W(xi) is calculated using W(xi)=σsτ(xi,xv)Δt.

5 Results

Figures 1 and 4-8 show our results. Table 1 shows the numbers of VPLs and shading points, as well as the computational time for our method and MDLC measured using a PC with an Intel Xeon E5-2650 v4 2.20 GHz CPU. In all the calculations, the relative tolerance ∈was set to 2%, and the α quantile for Eq. (6) was set to 95%. To measure the estimation accuracy, we defined R∈as the percentage of pixels satisfying ‖ ^I-I ‖＜∈I; the estimation is accurate when R∈is close to α (95% in our experiments).

Figure 1 shows results of our rendering method with anti-aliasing, where the insets compare (a) the reference solution, (b) our method, and (c) MDLC.The reference solution is rendered by summing up all the contributions from the VPLs at all the shading points to calculate the pixel values. Insets (d) and (e)illustrate the relative errors and R∈for our method and MDLC.This figure demonstrates that our method can improve upon the estimation accuracy achieved using MDLC,with corresponding improvements in noise,as shown in Fig.1(c). Figure 4 shows a San Miguel scene rendered with DOF. Figure 5 shows a Cornell box scene filled with homogeneous participating media.As shown in insets (c) of Figs. 4 and 5, the stochastic noise is perceptible when using MDLC, whereas it is imperceptible using our error estimation method (see insets (b)). Figure 6 shows a Sponza scene filled with homogeneous participating media. Figure 7 shows a Cornell box scene where the two boxes consist of translucent material. Figure 8 shows a human head model rendered with a diffuse BSSRDF and anti-aliasing. In all of these experiments, our method showed improved estimation accuracy R∈over MDLC.

Table 1 Statistics for our method and MDLC

Fig. 4 San Miguel scene with DOF.

Fig. 5 Cornell box scene with participating media.

Fig. 6 Sponza scene with participating media.

Fig. 7 Cornell box scene with two boxes of translucent materials.

Fig. 8 Human head scene (subsurface scattering with anti-aliasing(8 spp)).

In Figs. 9 and 10, cluster selections without scaling(the second column), scaling by cl(the third column),and scaling by cd(the fourth column) are compared.Scaling by cdyields the best estimation accuracy R∈for five out of six of these scenes. Moreover,comparing results for these scenes without scaling with those with scaling demonstrates that scaling by cdconsistently improves the estimation accuracy. In Fig. 9(a), the first row of the San Miguel scene with anti-aliasing shows relative errors greater than 5%,especially around the ivy-covered wall in the first column, but are reduced in the third column. In Fig. 9(d), the fourth row of the kitchen scene with DOF, relative errors greater than 5% can be seen,especially around the shadow boundaries and outlines.We think that these low estimation accuracies arise from the uneven cluster subdivisions (i.e., shading clusters are not subdivided, whereas VPL clusters are predominantly subdivided). Comparing scalings by cland cdshows that estimation accuracies in the three scenes deteriorate when using cl, although clyields the best estimation accuracy in the San Miguel scene with DOF (c). Based on these experiments, we used the scaling coefficient cdto select the clusters.

Fig. 9 Cluster selection without scaling (the second column), scaling by cl (the third column), and scaling by cd (the fourth column). The left column shows the reference solutions and the other images show the relative errors in false color. The values in each relative error image show R∈. Scaling using cd yields better estimation accuracies than not scaling or scaling by cl.

Fig. 10 Cluster selection without scaling (the second column), with scaling by cl (the third column), and scaling by cd (the fourth column) for subsurface scattering (above) and participating media (below). Although scaling by cl provides better estimation accuracy R∈than not scaling in the human head scene (above), relative errors greater than 5% can be seen in the neck. Scaling using cd provides the best estimation accuracy and relative errors greater than 5% are not seen. In the Cornell-box scene with homogeneous participating media (below), since VPLs and shading points are distributed throughout the scene, the improvements due to the use of cd are subtle, but its use does not impair the estimation accuracy.

While our method outperformed MDLC in estimation accuracy,its computational time was from 1.11 to 8.55 times greater. We attempted to adjust the parameter ∈so that the estimation accuracy of MDLC was similar to our result in Fig. 5. After several trial-and-error processes and re-renderings,R∈of MDLC became 91.2% with ∈= 0.5%. The computational time for MDLC with ∈= 0.5% was 33.5 s, which is comparable to our result (35.6 s).However, our method does not require the tedious trial-and-error processes needed for MDLC.

6 Conclusions and future work

We have presented a scalable many-light rendering method that can improve the estimation accuracy for various visual and illumination effects. Our method automatically partitions VPL and shading clusters so that pixel errors are smaller than the relative tolerance ∈.Currently, our method is limited to homogeneous participating media with isotropic scattering. We would like to lift this limitation in future work.Moreover, we intend to apply our method to motion blur.

Acknowledgements

This work was partially supported by JSPS KAKENHI 15H05924 and 18H03348.

Appendix Derivation of the error estimate ΔI

If the samples in each pair follow a normal distribution and Neyman allocation is used for the number of samples for each pair, their statistics T follow a tdistribution with (n-N) degrees of freedom:

where N is the number of pairs, niis the number of samples for the i-th pair (CLi,CGi), n is the total number of samples (i.e., n=N i=1ni), and siis the sample variance of the i-th pair. The error ΔI = ^I-I is calculated using the α quantile in the t-distribution,tα, as

To estimate the pixel value, it is known to be more efficient to use a large number of pairs with a small number of samples than a small number of pairs with many samples. Since at least two samples for each pair is required to calculate the sample variance si,we sample two VPLs and shading points from the pair (K =2 in Eq. (4)). By substituting ni=2 and n=2N in Eq. (15), the error I can be simplified to give Eq. (6)