Geometric Sample Reweighting for Monte Carlo Integration

TL;DR

This paper introduces a geometric sample reweighting method for Monte Carlo integration that improves convergence and accuracy by linking weight derivation to function reconstruction, applicable to rendering problems.

Contribution

It proposes a novel reweighting scheme based on geometric insights, enhancing Monte Carlo integration without bias and with simple implementation.

Findings

Better convergence than standard methods

Effective in Monte Carlo rendering tasks

Simple to implement and generalizable

Abstract

We present a general sample reweighting scheme and its underlying theory for the integration of an unknown function with low dimensionality. Our method produces better results than standard weighting schemes for common sampling strategies, while avoiding bias. Our main insight is to link the weight derivation to the function reconstruction process during integration. The implementation of our solution is simple and results in an improved convergence behavior. We illustrate its benefit by applying our method to multiple Monte Carlo rendering problems.

Equations

$${\hat{I}}_p=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)},$$

$$\Delta (x_i)=\frac{1}{Np(x_i)}.$$

$$\hat{I}=\sum _{i=1}^{N}w(x_i)f(x_i).$$

$$p_i(x)=\frac{x^{i-1}(1-x{)}^{N-i}}{{\int }_0^1{t}^{i-1}(1-t{)}^{N-i}\mathrm{d}t}.$$

$$\begin{array}{rl}E[x_i]& ={\int }_0^{1}x\cdot \frac{x^{i-1}(1-x{)}^{N-i}}{{\int }_0^1{t}^{i-1}(1-t{)}^{N-i}\mathrm{d}t}\mathrm{d}x\\ & =N\cdot \left(\genfrac{}{}{0ex}{}{N-1}{i-1}\right)\cdot {\int }_0^{1}x\cdot x^{i-1}(1-x{)}^{N-i}\mathrm{d}x\\ & =N\cdot \left(\genfrac{}{}{0ex}{}{N-1}{i-1}\right)\cdot {\int }_0^1{x}^i\cdot \sum _{k=0}^{N-i}\left(\genfrac{}{}{0ex}{}{N-i}{k}\right)1^{N-i-k}(-x{)}^k\mathrm{d}x\\ & =N\cdot \left(\genfrac{}{}{0ex}{}{N-1}{i-1}\right)\cdot {\int }_0^1\sum _{k=0}^{N-i}\left(\genfrac{}{}{0ex}{}{N-i}{k}\right)(-1{)}^k{x}^{k+i}\mathrm{d}x\\ & =N\cdot \left(\genfrac{}{}{0ex}{}{N-1}{i-1}\right)\cdot \sum _{k=0}^{N-i}\left(\genfrac{}{}{0ex}{}{N-i}{k}\right)\frac{(-1{)}^k}{i+k+1}\\ & =N\cdot \frac{(N-1)!}{(N-i)!(i-1)!}\cdot \frac{\Gamma (i+1)\Gamma (N-i+1)}{\Gamma (N+2)}\\ & =N\cdot \frac{(N-1)!}{(N-i)!(i-1)!}\cdot \frac{i!(N-i)!}{(N+1)!}=\frac{i}{N+1}.\end{array}$$

$${\hat{I}}_{GR}=\sum _{i=1}^N{w}_{GR}(x_i)f(x_i)\text{, where }{w}_{GR}(x_i)=\frac{\left|V_i\right|}{C(x_i)\left|\Omega \right|}.$$

$$E\left[V_i\right]$$

$$C(x_i)={\left(\frac{3}{2}\right)}^{b(x_i)}\frac{N}{n+1},$$

$$E\left[w_{GR}(x_i)\right]$$

$$=\frac{1}{\left|\Omega \right|}\frac{{\left(\frac{3}{2}\right)}^{b(x_i)}\frac{|\Omega |}{n+1}}{{\left(\frac{3}{2}\right)}^{b(x_i)}\frac{N}{n+1}}=\frac{1}{\left|\Omega \right|}\times \frac{\left|\Omega \right|}{N}=\frac{1}{N}.\text{hfill}■$$

$$C(x_i)={\left(\frac{3}{2}\right)}^{b(x_i)}\frac{N}{n+s}.$$

$$f(x)=10\times \{\begin{array}{ll}\sqrt{-x^2+0.5x}& 0<x<0.25\\ -\sqrt{-x^2+x-0.1875}+0.25& 0.25<x<0.5\\ 20\times (x-0.5)& 0.5<x<0.55\\ 1.0& 0.55<x<0.65\\ -20\times (x-0.7)& 0.65<x<0.7\\ 0.1\times \sin (10\pi \cdot (x-0.7))& 0.7<x<0.8\\ 0.25\times \sin (10\pi \cdot (x-0.8))& 0.8<x<0.9\\ 0.5\times \sin (10\pi \cdot (x-0.9))& 0.9<x<1.0\end{array}$$

$$L_{(i,j)}={\int }_{t_{\text{open}}}^{t_{\text{close}}}{f}_{(i,j)}(t)\mathrm{d}t,$$

$${\hat{I}}_{(i,j)}=\frac{1}{N}\sum _{k=1}^N\frac{f_{(i,j)}({\lambda }_k)}{p_{(i,j)}({\lambda }_k)}$$

$$L_{(i,j)}={\int }_{\mathbb{A}}{f}_{(i,j)}(s)\mathrm{d}s.$$

$$L_{\text{x}}$$

$$=L_e(\text{x})+{\int }_{\mathbb{L}}{f}_s(\text{x})L_e(l\to \text{x})\mathrm{d}l+{\int }_{\Omega }{f}_s(\text{x},\omega )L_i(\omega )\mathrm{d}\omega ,$$

Full text

Sampling and Reconstruction, Monte Carlo Integration, Sample Reweighting, Rendering

††copyright: rightsretained††ccs: Computing methodologies Ray tracing

1. Introduction

Monte Carlo (MC) techniques form the foundation of realistic image synthesis for decades (Cook et al., 1984). The principle is simple: a function is sampled and the samples are combined to approximate its integral. Standard MC is often referred to as brute-force, as its implementation is simple but the variance of the estimation can be high and convergence slow. One method to usually improve the integral approximation is to reconstruct the underlying function from the samples and much previous work devoted its attention to particular cases (e.g., shadows (Egan et al., 2009) or depth of field (Soler et al., 2009)). In this work, we revisit the reconstruction process. We derive an easy-to-implement algorithm to compute sample weights that generally improves the approximation when compared to standard weights for general MC integration.

Our observation is that standard sample weights are often less accurate for lower sampling rates because they do not properly reflect the integration domain nor the local sample density. Our weighting scheme considers all samples of a given set and defines weights based on a geometric partitioning of a low-dimensional integration domain. It results in a consistent estimator that outperforms standard weighting schemes. A major contribution of our work is the derivation of an unbiased estimator. It builds upon this partitioning and applies to sets of independent and identically distributed uniform random samples or stratified samples. Specifically, we propose a novel weighting scheme that is easy to implement and builds upon a sound theoretical derivation. It integrates well into existing rendering pipelines, can be parallelized in conjunction with the unbiased estimator, and we demonstrate its benefit over existing schemes via several rendering problems.

We will first cover prior work and MC integration. We then give the motivation behind our approach (Sec. 3) and present the core of our solution (Sec. 4). Numerical performance and applications to rendering are presented in Sec. 5.

2. Background

MC methods

Since the 80s (Cook et al., 1984), MC integration plays a major role in rendering complex effects, such as motion blur, depth of field, and soft shadows. The complete light transport is described by the rendering equation (Kajiya, 1986), which can be solved using path tracing as an associated MC solution. Nevertheless, not all samples taken during the evaluation of an integral contribute strongly to the result. One strategy to modify subsequent sample choices is to rely on previous samples, i.e., a Markov process. Metropolis sampling (Veach and Guibas, 1997) can handle complex light path configurations by extensively exploring contributing paths once they are discovered. Multidimensional k-d trees (Hachisuka et al., 2008; Guo et al., 2018) can store samples in a global structure, which can then be used as a means to control future sample placement.

While standard Monte Carlo (MC) methods solve a definite integration $I=\int_{\Omega}f(x)\mathrm{d\,}\!{x}$ of a function $f$ over a finite support $\Omega\subset\mathbb{R}^{d}$ by using a random sample set ($\{x_{i}\in\Omega\}$) with the resulting estimator being $\hat{I}_{MC}=\frac{1}{N}\sum_{i=1}^{N}f(x_{i})$, importance sampling influences the sampling process via a probability distribution function (pdf) $p:\Omega\rightarrow\mathbb{R}$ (Veach and Guibas, 1995). The resulting unbiased estimator is:

[TABLE]

which effectively weighs samples differently. Importance sampling is interesting when having knowledge about the scene. For instance, importance sampling the light source works better in scenes with small or point light sources (Dutre et al., 2006; Debevec, 2008). Sampling according to the BSDF works better with glossy to highly glossy surfaces(Shirley, 1991; Ward, 1992; Lafortune et al., 1997). Multiple importance sampling (MIS) combines different such sampling strategies (Veach and Guibas, 1995).

Reweighting

Our solution focuses on the weighting of samples interpreted as an improved function reconstruction. Different weight definitions have been shown to be beneficial for rendering, e.g., derived in Sobolev spaces (Marques et al., 2018). However, these previous solutions target hemispherical illumination integrals and are not generally applicable to other problems. A reweighting scheme was also proposed for addressing firefly artifacts (Zirr et al., 2018) but the solution is biased and limited to narrow application scenarios.

Other specialized reconstruction techniques exist, including solutions for soft shadows (Egan et al., 2009), defocus blur (Soler et al., 2009), and motion blur (Egan et al., 2011), which lead to significant improvements. More complex reconstructions for light fields (Lehtinen et al., 2011) have proven very successful but are biased (though consistent).

Our method is independent of the application scenario and unbiased. It handles general functions and links the weights to Voronoi cell volumes. The latter has also been studied in the context of anti-aliasing problems (Mitchell, 1990), for which the $2$ dimensional voronoi cell volumes bounded with a pixel are directly used as sample weights and leads to improved anti-aliasing effects, but the theory has not been further developed for unbiased solutions, nor generalized to other contexts. Voronoi cell size has been used as weights for Monte Carlo integration in (Vorechovsk et al., 2016), where two ways of treating boundaries have been proposed. In this work, Voronoi cells of given set of samples within a domain are either bounded and clipped by the domain boundary, or extended by periodically adding auxiliary samples that extend the domain. Both approaches are shown to improve numerical performance of MC integrations. However, as we show in the Sec.4, directly using Voronoi cell size as weight results in an biased estimate. Our solution takes advantage of Voronoi tessellation and remains unbiased.

3. Formulation and Problem Statement

Referring to Eq. 1, the estimator of importance sampling is a sum of function value $f(x_{i})$ times a weight $\Delta(x_{i})$, generally:

[TABLE]

Similarly, Riemann integration approximates an integral using function values $f(x_{i})$ times a weight $w(x_{i})$:

[TABLE]

The Riemann weights stem from a partitioning of the support $\Omega$ into hypervolumes. In 1D, these hypervolumes are intervals. Each hypervolume contains exactly one sample and its volume defines the sample’s weight.

The weights $\Delta(x_{i})$ are typically easy to compute but cannot be considered hypervolumes; they would overlap or introduce gaps and cannot easily be linked to a partitioning of $\Omega$. Only with increasing number, due to the stochastic nature of the process, when the samples densely cover the support $\Omega$, the difference in the weight definitions becomes negligible. See Fig. 1 (a), (b) and (c) for an illustration. In consequence, especially for low sample counts, the weights do not well reflect an approximation of the function.

4. Geometric Sample Reweighting

Our goal is to associate weights to samples that define an improved function reconstruction during the integration. We will first define a consistent solution, inspired by Riemann integration. This solution is independent of the sampling pattern and can be applied on any sample set as a post process to improve the approximation. This reweighting is consistent, but not unbiased for all sampling strategies. We then propose a modification to obtain an unbiased estimator for the cases of uniform random samples that are independent and identically distributed (i.i.d.), and samples with stratification. See Sec section 6 for the possibility to generalize our method for samples generated with an analytically known pdf.

4.1. Consistent Estimator

Riemann integration typically assumes a regular partitioning of the domain. Using a Voronoi diagram of the sample points, it is possible to partition the domain $\Omega$ to take sample density into account. A Voronoi diagram is a partition into regions such that the points in each region share the same closest sample location. It can be shown that the Voronoi cell corresponds to the intersection of half spaces defined by hyperplanes that are equidistant to two sample points. The theory of Voronoi diagrams is beyond the scope of this paper but more details can be found in (Aurenhammer, 1991; De Berg et al., 1997).

In our case of a $D$ dimensional problem setting, the diagram will be bounded by the hypercube $(0,1)^{D}$, the domain from which samples are drawn. The volume of each Voronoi cell determines the corresponding sample weight and given that the cells are intersections of half-spaces, they are convex and their volume can be easily computed. The resulting estimator of our approach is $\hat{I}_{CON}=\sum_{i=1}^{N}w_{CON}(x_{i})f(x_{i})\text{, where }w_{CON}(x_{i})=\frac{\left|V_{i}\right|}{\left|\Omega\right|}$.

Implicitly, this construction approximates the integrand via a piecewise-constant representation. Intuitively, to take the most benefit from this interpretation, samples should be chosen with respect to the gradient of the function. Fig. 1 shows an illustration of this strategy. In principle, even more advanced approximations could be used, yet it turns out that such weight definitions, while consistent, lead to a biased estimate. In the following, we will show the reasoning behind this and derive an unbiased estimator for i.i.d. uniform sampling and stratified sampling.

4.2. Deriving an Unbiased Estimator

i.i.d. Uniform Samples

The reason the direct use of the Voronoi cells’ volume is biased is due to the samples whose cell shares a boundary with the domain boundary. To illustrate this situation, we will first consider the 1D case with a set $X$ of $N$ (with $N>3$) i.i.d. uniform samples $X:=\left\{x_{i}\in(0,1)\right\}$ before generalizing to $D$ dimensions, and then to stratified sampling.

One Dimension

Let us assume that the one-dimensional sample set $X$ is sorted from smallest to largest value. We are interested in the expected extent of each Voronoi cell, for which we need to derive the expected distance between two adjacent samples. For this reason, we first determine the expected position of sample $x_{i}$.

From order statics (David and Nagaraja, 2004), we know that the distribution of the $i$-th i.i.d. sample follows the beta distribution, i.e.,

[TABLE]

The expected position of the ordered $i$-th sample $x_{i}$ is then:

[TABLE]

Consequently, we have $E[|x_{i}-x_{i+1}|]=\frac{1}{N+1}$ for $i=1$ to $N-1$ and a similar condition holds for $E[x_{1}-0]=E[1-x_{N}]=\frac{1}{N+1}$. The expected weight is then $\frac{1}{N+1}$ for samples $x_{i}$ with $i=2$ to $N-1$ and $\frac{3}{2}\frac{1}{(N+1)}$ for samples $x_{1}$ and $x_{N}$. The latter weights are larger due to the intervals containing the two boundaries of the domain. Using these weights directly, leads to a consistent but biased estimator.

To render the estimator unbiased, we introduce a per-sample correction coefficient $C$:

[TABLE]

These factors have to be carefully chosen — for instance, $C=1$ would lead to the previously-derived consistent but biased result. The correction coefficient should indeed modify the expected contribution of a sample $x_{i}$ to equal $\frac{1}{N}f(x_{i})$. Following the weight derivation, an unbiased estimator in 1D, we would then define $C(x_{1})=C(x_{N})=\frac{2(N+1)}{3N}$ and $C(x_{i})=\frac{N+1}{N}$ for all other samples. As most samples still share an identical correction factor, it keeps us close to the interpretation of the Voronoi cell volume. In higher dimensions, the definition is less straightforward.

D Dimensions

To derive the correction coefficient $C$ from Eq.4 in $D$ dimensions, we assume a set of $N$ (with $N\geq 3^{D}$, i.e., intuitively, this results in at least one inner point and two boundary points along each dimension) samples in $\Omega=\left(0,1\right)^{D}$. We define the boundary order $b(x_{i})$ of a sample as the amount of cell boundaries of its Voronoi cell that are part of the domain boundary. For instance, in the above one dimensional example, $b(x_{1})=b(x_{N})=1$, and for all other sample points, we have $b(x_{i})=0$.

The cardinality of samples of order $d$ is defined as: $|X_{d}|=\\ \textup{card}\left\{b(x_{i})=d,\forall i\in\left[1,N\right]\right\}$. For such a sample set of $N$ samples, the expected cardinality of samples of order $d$ is $E[|X_{d}|]=\binom{D}{d}(\sqrt[D]{N}-2)^{d}2^{D-d}$. This formula is the $d$-th term in the bionomial expansion of $[$($\sqrt[D]{N}-2$)$+2]^{D}$. To understand this result, one should recall that the expected position of all samples forms a regular grid. Thus, this grid will have a resolution of $n=\sqrt[D]{N}$ along each axis. Starting with one axis, we would find $n$ samples with two boundary samples of order one and all others samples are inner points of order zero. Repeating these samples $n$ times along a new dimension will increment the order of the first repeated set of samples and the last, as these represent a new boundary along this dimension. For all other samples, their boundary order remains unchanged. This process can be done for all $D$ dimensions, thus implying the binomial expansion.

To achieve an unbiased estimator, we first compute the expected Voronoi volume $E\left[V_{i}\right]$ for a sample $x_{i}$. For $D$ dimensions, we have $D+1$ boundary orders from [math] to $D$. As we are dealing with an i.i.d. uniform distribution, in each dimension, we have $n-1$ intervals between samples and two intervals with the boundary, leading to a total of $n+1$ intervals. Therefore, we have:

[TABLE]

Again, for unbiasedness, we need $E\left[w_{GR}(x_{i})\right]=1/N$, thus each sample should expectedly contribute equally. The following definition of the correction coefficients fulfills this property:

[TABLE]

because

[TABLE]

Stratified Samples

The extension to stratified sampling is relatively straightforward, as each stratum is considered an independent unit. This means that the function is independently integrated in each stratum and its whole range is a composition of these units. In consequence, the boundary observation now applies to the boundary of each stratum. For a sample set $X$ of size $N$ generated with $S$ strata, each stratum is expected to contain $\frac{N}{S}$ samples. Let $n=\sqrt[D]{N}$ and $s=\sqrt[D]{S}$, then the integration problem for a stratum with $\frac{N}{S}$ samples would imply the correction coefficients to be:

[TABLE]

5. Results

In this section, we first show the numerical performance (Sec. 5.1) of our scheme and show its application to a few rendering scenarios (Sec. 5.2). In all tests, we compare four different estimators:

1. (1)

Standard MC (i.i.d. uniform sampling) 2. (2)

Our weighted standard MC (i.i.d. uniform sampling) 3. (3)

Stratified MC (stratified sampling) 4. (4)

Our weighted stratified MC (stratified sampling)

5.1. Numerical Performance

The numerical performance is tested with two examples: one for a $1$D MC integration and the other for a $2$D MC integration, which are plotted in Figure 2. The $1$D function is given as:

[TABLE]

For the $2$D function, we take the Lena image (Munson, 1996). The functions were chosen to include discontinuities, large-scale variations and small scale changes and led to a representative behavior of several tests that we have performed. Generally, the MSE drops as more samples are added (Column $1$). Our solutions outperform standard uniform sampling and even stratified sampling by several orders of magnitude and converges around 1000, 100 times faster respectively in 1D and 100, 10 times faster respectively in 2D.

For the case of stratified sampling, we illustrate different amounts of strata for the same sampling count (Column $2$). Our weighting scheme makes this parameter less important, as it achieves a better function approximation.

We next investigate the impact of distributing samples into batches for which we estimate the function integration separately, before deriving the overall estimate by averaging, which would typically be the case for distributed computations. First, we fixate the amount of samples to $100$K (Column $3$). Notice that the performance of standard uniform sampling remains invariant with respect to the amount of samples per batch, as it is already an averaging process. Our solution results in a better approximation for more samples per batch, as it will approximate the function more faithfully, as expected. Similarly, stratified sampling also benefits from more samples per batch, but shows slower convergence.

We also investigate the effect of using different batch sizes for uniform (Column $4$) and stratified sampling(Column $5$). More batches thus means a higher overall sample count and all methods improve with the addition of batches. In all cases, the graphs stop after reaching $100$K samples. Our solution performs best and the graphs also illustrate the convergence over several batches, due to its unbiasedness.

5.2. Application to Rendering

We implemented our method in Mitsuba (Jakob, 2010), targeting one and two dimensional integration problems, namely motion blur (Sec. 5.2.1), dispersion (Sec. 5.2.2), depth of field (Sec. 5.2.3) and illumination integrals (Sec. 5.2.4). We evaluate MSE and visual appearance, as well as convergence behavior. For all implementations, our reweighting operates at a per-pixel level. We apply our method on the level of primary samples, thus all applicable local importance sampling techniques are utilized throughout the pipeline.

5.2.1. Motion Blur

To simulate motion blur, distribution rendering samples the time domain: For a pixel $(i,j)$, the luminance $L_{(i,j)}$ is given by:

[TABLE]

with $t_{\textup{open}}$ and $t_{\textup{close}}$ being the shutter opening and closing time and $f$ incorporating the shutter function. Since time is $1$ dimensional, building a Voronoi partition means sorting and measuring the distance between samples. We tested our implementation in two scenes with animation (Fig. 3 and 4).

5.2.2. Spectral Rendering

Light dispersion can happen at reflective or refractive dielectric materials, leading to effects such as rainbows, resulting from different wavelengths travelling in different directions. Spectral sampling simulates multiple wavelengths in order to capture such effects. To reduce the complexity of the additional spectral dimension(Bergner et al., 2009), hero wavelength spectral sampling(Wilkie et al., 2014) can be used as an approximation:

[TABLE]

Our implementation of spectral sampling uses $15$-bin wavelengths. Hero wavelength sampling is used with $3$ shifted additional wavelength samples (Wilkie et al., 2014). We tested our method with two scenes configured with dispersive di-electric materials (Fig. 5 and 6).

As shown in the results, our method brings down colour noise significantly and dispersive regions look much smoother at low sample rate.

5.2.3. Defocus Blur

A camera with aperture leads to defocus blur/depth of field effects. The aperture is usually modelled as a $2$D shape, e.g., a square, a circle, or a star, which is sampled to determine the origin of each primary sample ray, which passes through the position on the focal plane corresponding to the current pixel. For lens aperture $\mathbb{A}\subset\mathbb{R}^{2}$, we obtain:

[TABLE]

To determine our weights, we use a $2$D Voronoi diagram based on the aperture samples. We tested a simple glossy sphere illuminated using an environment map (Fig. 7).

5.2.4. Direct Illumination

Leaving out irrelevant terms, the luminance $L_{\textbf{x}}$ at scattering point x with one bounce is given by:

[TABLE]

where $L_{e}$ denotes light emission and $l\in\mathbb{L}$ denotes all light sources. In this application, we use light sampling instead of random rays to ensure that the light source is always sampled. Our unbiased reweighting achieves the best convergence and, as shown in the insets, also the smoothest results (Fig. 8).

Observation

In all cases, our solution leads to smoother visual result and less black holes in the falloff regions. From the MSE plots, we can see that standard MC with uniform sampling has the worst performance, while our weighted stratified sampling generally has the best one. Our method improves both uniform sampling and stratified sampling. We can also see that even with uniform sampling as input, our weighted uniform sampling not only improves over the unweighted version, but also has a performance that is as good as our weighted stratified sampling. Precompute sample weights enables a negligible computation overhead.

6. Conclusion

The reweighting scheme in this paper enables a better approximation than standard MC weights. Our solution is general and does not require any prior knowledge about the integrating function. Implicitly, our method approximates this function via a reconstruction from the samples, but does not introduce a bias in the resulting estimator. We showed its practical benefit for various rendering problems.

While we focus on primary samples that are either i.i.d. uniform or stratified in this work, our method can also handle non-uniform sample sets following a distribution of $p(x)$. The expected position of the $i$-th sample $x_{i}$ is then $N\cdot\binom{N-1}{i-1}\cdot\int_{0}^{1}P^{-1}(x)\cdot x^{i-1}(1-x)^{N-i}\mathrm{d\,}\!{x}$, where $P(u)=\int p(u)\mathrm{d\,}\!{u}$. Unfortunately, it is necessary to integrate the distribution function. Approximate schemes remain an area of future work. Similarly, using the method in higher dimensions requires the computation of cell volumes in high-dimensional Voronoi diagrames, which can be costly. One could precompute these weights but we left such accelerations as future work. Finally, it is an exciting opportunity to exploit the generality of our solution to improve other integration problems.