The stochastic nature of modern Monte Carlo (MC) rendering methods inevitably produces noise in rendered images for a practical number of samples per pixel. The problem of denoising these images has been widely studied, with most recent methods relying on data-driven, pretrained neural networks. In contrast, in this paper we propose a statistical approach to the denoising problem, treating each pixel as a random variable and reasoning about its distribution. Considering a pixel of the noisy rendered image, we formulate fast pair-wise statistical tests—based on online estimators—to decide which of the nearby pixels to exclude from the denoising filter. We show that for symmetric pixel weights and normally distributed samples, the classical Welch t-test is optimal in terms of mean squared error. We then show how to extend this result to handle non-normal distributions, using more recent confidence-interval formulations in combination with the Box-Cox transformation. Our results show that our statistical denoising approach matches the performance of state-of-the-art neural image denoising without having to resort to any computation-intensive pretraining. Furthermore, our approach easily generalizes to other quantities besides pixel intensity, which we demonstrate by showing additional applications to Russian roulette path termination and multiple importance sampling.
1 Introduction
Monte Carlo (MC) integration is ubiquitous in modern photorealistic rendering. The basic idea is to estimate an integral (pixel intensity) by a (weighted) average of random samples. However, compared to deterministic quadrature rules, MC integration is stochastic in nature and must therefore contend with uncertainty. This stochastic uncertainty manifests as variance—i.e., noise—affecting every sample taken during rendering.
Fig. 1:
Fig. 1: Our statistical denoising method, using online estimates of sample statistics, achieves image quality comparable to current state-of-the-art methods, without any computation-heavy prior training. We compare our denoiser to the approach by Moon et al. [2014] (“Moon CI”), NVIDIA OptiX AI-Accelerated Denoiser (“OptiX”), Intel Open Image Denoise (OIDN), and progressive denoising [Firmino et al. 2022] (“ProDen”). These results have been generated using 256 samples per pixel (SPP) with the following denoising times; Moon CI: 35.3 ms, OptiX: 85.5 ms, OIDN: 19.5 ms, ProDen: 1834.9 ms, ours: 28.0 ms.
The field of MC denoising originated from the idea of applying traditional image filters from the image-processing community to reduce noise in rendering. These filters essentially average multiple estimates (i.e., pixels) in image space, thereby decreasing variance. This approach can be highly effective in situations where the estimates are the same (or similar) in the limit, i.e., in regions of an image that show smooth color gradients. However, in cases of significant difference, such as along feature edges, filtering potentially introduces unwanted bias, often apparent as excessive blurring or color bleeding.
Denoising therefore must achieve a good bias–variance trade-off. In contrast to traditional image filtering, recent denoising approaches rely on neural networks to learn a mapping from a noisy input image, along with auxiliary information (in the form of so-called G-buffers containing low-noise ground-truth scene information), to a low-error approximation of the desired image.
Departing from these pretrained machine learning approaches, in this paper, we propose a general statistical framework for denoising in MC rendering. Specifically, we establish a theoretical connection between minimizing mean squared error (MSE) for pair-wise symmetric weights and Welch’s t-test for normally distributed samples. Using more recent results from the statistics literature, we then generalize this approach to reduce the assumptions on the distributions. As a practical implementation of this framework, we demonstrate an image-space denoising scheme, building upon the well-known joint bilateral filter. Our approach tracks online statistics of the per-pixel samples generated by a state-of-the-art MC renderer. These statistics describe the distribution of samples and allow us to draw statistically valid conclusions about which estimates among neighboring pixels should be combined in order to achieve a near-optimal trade-off between noise and bias. Our statistics-based filter effectively avoids blurring over legitimate features, including lighting effects (such as shadows or caustics) that are not present in the G-buffers (and would therefore be over-blurred by existing denoisers). All these statistics can be estimated online during rendering. Our filter runs on the GPU, taking around 30 ms for 1280 × 720 resolution images on commodity hardware.
Apart from image denoising, we demonstrate that our approach can generalize to other applications, such as Russian roulette (RR) path termination and multiple importance sampling (MIS). Our general approach has an important advantage relating to these generalizations: We do not rely on preexisting training data, therefore we can produce low-error estimates of any quantity of interest, not only the radiance produced by the renderer. For instance, we denoise the per-bounce incident radiance in our RR example or estimates for which sampling strategy outperforms the others (formalized as win rates) during MIS. For both quantities, obtaining sufficient training data for neural denoising approaches would be quite challenging.
In summary, we present the following main contributions:
•
a general statistical framework that puts MC rendering into a statistical context, alongside a theoretical analysis on minimizing mean squared error in the setting of pair-wise testing,
•
an image-space implementation of this framework that produces low-error estimates for any quantity of interest in the context of MC rendering,
•
and applications of this denoising approach to standard image denoising, RR path termination, and MIS.
The stochastic nature of Monte Carlo (MC) rendering introduces error in the form of variance. To reduce this error, one can either compute more samples or use some form of noise reduction, e.g., adaptive sampling or filtering [Huo and Yoon 2021; Zwicker et al. 2015]. A common approach is a posteriori noise reduction or denoising, which operates on the samples generated by MC rendering. Many methods work with (image-space) pixel estimates and range from classical denoising approaches [Overbeck et al. 2009; Xu and Pattanaik 2005] to more recent neural-network-based methods [Back et al. 2022; Firmino et al. 2022].
Classical approaches usually adapt and apply some form of filter while trying to find a balance between noise reduction and introduced bias (usually apparent in the form of excessive blurring). Many use image-filtering kernels, e.g., Gaussian [Rousselle et al. 2011], (joint) bilateral [Li et al. 2012; Liu et al. 2018b; Mara et al. 2017; Park et al. 2013; Rousselle et al. 2013; Sen and Darabi 2012; Xu and Pattanaik 2005; Zheng and Liu 2018], non-local means [Delbracio et al. 2014; Moon et al. 2013; Rousselle et al. 2012; Vicini et al. 2019], wavelet [Dammertz et al. 2010; Kalantari and Sen 2013; Overbeck et al. 2009; Schied et al. 2017; 2018], or non-local Bayes [Boughida and Boubekeur 2017] filters. Other methods build on diffusion [McCool 1999] or higher-order regression [Bauszat et al. 2011; Bitterli et al. 2016; Liu et al. 2018a; Moon et al. 2014; 2015; 2016; Yuan and Zheng 2017; 2018]. In addition to the noisy input image, many of the mentioned approaches incorporate additional auxiliary features like albedos and normals, as well as additional statistics like the corresponding variances. Among these, joint bilateral filtering with auxiliary features is most closely related to our denoising approach.
In a similar direction to our statistics-based approach, Sen and Darabi [2012] presented a method called random parameter filtering (RPF), which uses histograms of sample vectors to calculate the mutual information between scene features and the random parameters of the MC process in order to adjust the weights of a joint bilateral filter and reduce the noise stemming from these random parameters. Even though their approach achieves good results at low sample counts, its computation becomes prohibitively expensive with an increasing number of samples. This downside was addressed by Park et al. [2013] by interpolating sparsely computed mutual information, which reduces but not fully eliminates the scaling issues for high sample counts and requires several seconds or multiple minutes to denoise a low-SPP image. In contrast, our approach is completely independent of the number of samples and requires only a few milliseconds to compute. A similar idea is to represent sample distributions as histograms [Boughida and Boubekeur 2017; Delbracio et al. 2014]. However, histograms require additional memory and their accuracy depends on the number and sizing of bins, which must be chosen a priori. In contrast, we only use a few statistical moments to compactly store information about sample distributions, which can also be efficiently updated with each new sample.
Li et al. [2012] surpass the denoising quality of RPF by using Stein’s unbiased risk estimator (SURE) to select which one out of a discrete set of joint bilateral filters (of different spatial scales) minimizes the error. Similarly, Rousselle et al. [2013] use SURE to select among filters that differ in robustness to noise and sensitivity to image details. Following this general idea of “meta”-approaches, Zheng et al. [2021] combine multiple denoisers, whereas Firmino et al. [2022] apply denoising only when beneficial for convergence. In contrast, our approach does not have the overhead of computing multiple filters, or corresponding derivatives for error estimation, before the final reconstruction. We instead derive filter kernels directly from the pixel statistics and use statistical testing to ensure that only estimates with similar distributions are combined, leading to a fast one-pass approach, which could also improve the results of meta-denoisers by acting as an additional base input.
Mara et al. [2017] presented a real-time approach using a pipeline consisting of several joint bilateral filter stages, and showed results comparable to non-local means (NLM) [Rousselle et al. 2012] and regression-based approaches [Bitterli et al. 2016; Moon et al. 2014]. Their method is specifically tailored for real time and builds on a custom renderer with special handling of direct and indirect illumination. Our approach uses a more straightforward pipeline, which can be easily integrated into any conventional MC renderer with little overhead while offering similar performance during filtering.
Similar to ours, some methods use confidence intervals [Moon et al. 2013; 2016; Sen and Darabi 2012] to exclude pixels from combination during filtering if their corresponding statistics differ significantly. These approaches use confidence intervals in specific subparts of their fairly complex pipelines, which are tailored specifically to MC image denoising. Similar statistical concepts were used by Back et al. [2023], who first applied the concept of uncorrelated statistics to denoising MC renderings. In contrast, we use statistical tests at the core of our efficient denoising pipeline, including higher-order central moments to relax the requirement of normally distributed samples implied by commonly used confidence intervals.
Neural-network-based approaches represent the current state of the art in MC denoising, e.g., NVIDIA OptiX AI-Accelerated Denoiser (OAADN), which is based on Chaitanya et al.’s work [2017], or Intel Open Image Denoise (OIDN) [Áfra 2024]. Both approaches take a noisy image and so-called G-buffers (such as normals or albedos) as inputs. The neural networks then map those inputs to a low-error approximation of the ground truth. For this purpose, the networks have to be trained on numerous input–reference image pairs, so that an effective mapping can be established. Reconstructing image features that suffer from high variance and are not present in the G-buffers, such as lighting effects, is a difficult problem. In these cases, the networks must rely solely on the noisy image to differentiate between legitimate features and noise. Furthermore, the performance of such approaches hinges on the training and available data. Generally, there are no guarantees concerning their convergence behavior. For instance, OIDN [Áfra 2024] includes images at different degrees of convergence into the training to incorporate a notion of convergence. This stands in contrast to several non-neural methods that provide consistency [Back et al. 2023; Bitterli et al. 2016; Moon et al. 2013; Rousselle et al. 2013]. To alleviate this fundamental problem, Firmino et al. [2022] train an additional neural network that predicts per-pixel mixing factors to reduce the weight of the biased reconstruction in favor of the converging input. They use confidence intervals to ensure convergence by limiting the neural mixing weights such that the biased contribution vanishes as the variance decreases with increased sample counts.
Our approach, in contrast, does not require any training, as it directly uses the statistics gathered during rendering to infer relationships between pixels for filtering. Furthermore, we provide a theoretical guarantee concerning the consistency of our filter. We also show in our results that we can denoise quantities other than the radiance reaching the camera. In this paper, we always keep all filtering parameters constant as the number of samples per pixel (SPP) increases to demonstrate the convergence behavior of our method. In the future, automatic parameter tuning as done in neural denoising or integration of our method into meta-denoising strategies could further improve denoising performance.
3 Background and Notation
In this section, we summarize the statistical concepts and methods on which we build our general statistical framework in §4 and establish the connection to the common denoising problem in Monte Carlo rendering. As a starting point, we consider the well-known rendering integral [Kajiya 1986], where each pixel Ii in an image is formed as
where Wi describes the sensor response and L′ is the radiance incident on the i-th pixel, which covers the solid angle Ωi from the camera location x. This radiance must satisfy Kajiya’s rendering equation throughout a virtual scene. Each pixel intensity is commonly approximated by Monte Carlo (MC) integration:
where fk refers to evaluations of the integrand in Eq. (1) for ni randomly drawn samples from Ωi with probability density pk.
We now move on to a more abstract statistical interpretation of the rendering process outlined above. In particular, we interpret each pixel as a statistical estimator \({\hat{\theta }_i(X_1,\dots ,X_{n_i})}\) for the (unknown) ground-truth value θi (in this case Ii, although the same idea applies to other estimators), where we view each sample Xk as a random variable (i.e., a distribution from which the rendered sample fk/pk is drawn). We generally assume unbiased estimators, i.e., \(\mathop{\mathbb {E}}[\hat{\theta }_i] = \theta _i\).
Using results from descriptive statistics, our core idea is to formulate statistical tests to decide which estimators are sufficiently similar so that they can be combined during denoising. In general, we employ the following descriptive statistics of the sample distribution: mean (μ), and the central moments of orders 2–4 (Ml, l ∈ 2, 3, 4), or their standardized variants, variance (σ2 = M2), skewness (\(M_3 / M_2^{3/2}\)), and kurtosis (\(M_4 / M_2^2\)) [Kenney and Keeping 1951]. The central moments are defined as
where \(\bar{X}\) is the sample mean. Crucially, all these statistics can be computed online [Meng 2015], building on the classical result by Welford [1962], i.e., by updating them one sample at a time during rendering without having to store all samples in memory. In our implementation, we use central moments up to the third order.
3.1 Denoising Estimators
In this work, we consider general denoising filters that construct the denoised estimator \(\tilde{\theta }_j\) as a convex combination of (noisy) input estimators \(\lbrace \hat{\theta }_i\rbrace\) :
where wij is the weight assigned to estimator \(\hat{\theta }_i\); all weights must be non-negative, and \(\sum \nolimits _i w_{ij} = 1\). The key question now is: how to find appropriate (sparse) weights wij? In our approach, we split these weights into two parts, Eq. (10): a base filter, such as the well-known joint bilateral filter, and a novel pair-wise statistical membership function mij, which decides which combinations of estimators are admissible.
Formally, the mean squared error (MSE) of the combined estimator \(\tilde{\theta }_j\) can be decomposed into variance and bias
where θj refers to the estimand, i.e., the unknown ground-truth value for the j-th estimator. Note that even though we assume unbiased estimators \(\hat{\theta }_j\), their filtered counterparts \(\tilde{\theta }_j\) generally contain some bias. The variance and bias are given by
where \(\text{Bias}(\hat{\theta }_i, \theta _j) = \mathop{\mathbb {E}}[\hat{\theta }_i]-\theta _j\).
The overall goal of MC denoising in this context is to find weights that minimize the total MSE, or a similar error metric, i.e., achieve an optimal trade-off between variance (noise) and bias, across all denoised estimators (i.e., the whole image):
However, this optimization problem cannot be addressed directly, as the ground-truth values are of course unknown, and we must therefore work with noisy estimates of variance and bias.
4 Statistical Filtering Framework
In this section, we develop our statistical filtering framework. Considering an abstract set of estimators \(\lbrace \hat{\theta }_i\rbrace\), we ask under which conditions combining a subset of these estimators improves mean squared error. Here, we approach the inherent noise-to-bias trade-off from a statistical perspective, formulating a membership function (m), fulfilling the requirements summarized in §4.1. We show how this approach relates to hypothesis testing in statistics, before describing the specifics of our image-space implementation in §5.
4.1 Problem Statement
Considering the general denoising problem introduced in §3, we formulate membership functions that decide which combinations of estimators are admissible during denoising, i.e., under which conditions the combination is likely to improve image quality. In particular, we require the following important properties:
(a) Pair-wise evaluation: for any pair of estimators \((\hat{\theta }_i, \hat{\theta }_j)\), m must be defined as a function of these estimators’ statistics only, i.e., \({m_{ij} = m(\mathcal {S}({\hat{\theta }_i}), \mathcal {S}({\hat{\theta }_j}))}\), where \(\mathcal {S}({\hat{\theta }_i})\) denotes descriptive sample statistics, such as \(\mu _i, \sigma ^2_i, M_{3,i}, M_{4,i}, \ldots\), of estimator \(\hat{\theta }_i\).
(b) Online statistics: all these statistics \(\mathcal {S}\) must be computable by an online algorithm, i.e., by updating their value one sample at a time during Monte Carlo (MC) rendering. Together with the pair-wise property, this requirement ensures that all components of our filtering pipeline can be implemented efficiently in terms of both parallel execution and memory usage.
(c) Symmetry: we require that mij = mji. In practice, symmetry of weights enforces energy preservation during filtering. Relaxing this requirement sometimes produces visually more pleasing results around bright outliers (“fireflies”) at the cost of losing some overall brightness and possibly slightly higher mean squared error (MSE).
(d) Convergence: In order to guarantee convergence of the denoised result \(\tilde{\theta }_i \rightarrow \theta _i \; \forall i\) with increasing sample size, we require that the membership function satisfies
In other words, as the variance approaches zero, the membership function must exclude estimator j from the combination with estimator i if there is any difference in their estimands, which would introduce bias.
(e)Identity: An estimator does not introduce additional bias to itself; we therefore set mii = 1 by definition.
Finally, we can state our main problem as finding membership functions mij that satisfy properties (a–e) and deliver the best possible variance–bias trade-off within these constraints.
4.2 Our Approach
Here, we describe our general framework for formulating statistics-based membership functions respecting the aforementioned requirements. We then use these functions to address the denoising problem, Eq. (8), by defining the filter weights in Eq. (4), as follows:
The first component of the weights, ρij, provides the option to integrate existing filters based on a priori available information into our system. The main reasons for doing so are: (1) to limit the size of the filter kernel, thereby improving the runtime performance by limiting the amount of required membership-function evaluations; (2) to be able to build upon well-known existing approaches; and (3) to include additional low-noise a priori information available from the renderer (G-buffers, e.g., normals or albedos in image space). Note that while a priori information is often useful to preserve some features due to geometry or textures, it is inherently useless for other features that only manifest themselves through sampling (e.g., shadows or caustics). In this regard, our membership function and this a priori weight complement each other, incorporating both empirical and a priori information. In our implementation (§5), we set ρij to represent a joint bilateral filter.
Existing filters are highly effective at reducing variance. The key task of the membership function mij is therefore to limit bias by excluding estimators from the filter whose estimands differ significantly. Following the requirements stated in §4.1, we now consider the MSE minimization problem, Eq. (8), for a pair of input estimators \(\hat{\theta }_i, \hat{\theta }_j\). These estimators are combined with weight w into denoised estimators \(\tilde{\theta }_i = w\hat{\theta }_i + (1-w) \hat{\theta }_j\) and \(\tilde{\theta }_j = w\hat{\theta }_j + (1-w) \hat{\theta }_i\). Note that we enforce symmetry of the weighting here, wij = wji = (1 − w). Assuming unbiased input estimators, the bias of their combination is \(\text{Bias}(\tilde{\theta }_i,\theta _i) = {\mathop{\mathbb {E}}[\tilde{\theta }_i-\theta _i]} = {(w-1) \theta _i + (1-w) \theta _j}\), and analogously for \(\text{Bias}(\tilde{\theta }_j,\theta _j)\). Minimizing the sum of the error terms, Eq. (5), for both denoised estimators then reads:
For this objective, we can find the optimum by setting the derivative with respect to the weight to zero and solving for w* (assuming at least some uncertainty, i.e., \(\text{Var}(\hat{\theta }_i) + \text{Var}(\hat{\theta }_j) \gt 0\)):
Note, however, that we do not know ground-truth values, or estimator variances, and therefore work with (noisy) estimates of these quantities instead, producing noisy results for w*. Consequently, the trivial choice of setting mij = 1 − w* and ρij = 1 does not yield satisfactory results. Our approach is therefore to use an existing smoothing filter for ρij and enforce a binary membership function such that
Here, γ is a threshold that determines how discriminative the membership function is. In this way, we effectively use the pair-wise optimal weight as a test statistic (similar to those used in statistical hypothesis testing) and prevent the estimators from being filtered together if the introduced bias (i.e., the difference of their means) is too large relative to the sum of their variances. In our supplementary document, we show that this test is equivalent to Welch’s t-test, t < γw with \(\gamma _w = \sqrt {1/(2\gamma) - 1}\). Our results use a critical value from Student’s t-distribution of γw = t1 − α/2, ν, with the significance level α = 0.005 and ν = ni + nj − 2, i.e., the upper bound for the degrees of freedom approximated by the Welch–Satterthwaite equation [1946].
Welch’s t-test is related to the confidence interval for the difference of two normally distributed means, see, for example, Eq. (18) in the paper by Curto [2023]. Consequently, we can relax the normality assumption by choosing a different confidence-interval formulation from the literature. Various formulations are available that consider not only the sample variance but also higher-order statistics, such as skewness, of the sample distribution. In our results, we find that the correction of the means proposed by Curto [2023], Eq. (20) there, yields good denoising behavior when combined with a Box-Cox transformation of the samples as detailed in §5.
As an orthogonal extension, we note that allowing asymmetric membership functions can produce visually more appealing results; while fewer constraints can also lead to improved MSE, the potential energy loss due to asymmetric weights may cause overall darker images and therefore slightly worse quantitative errors (Fig. 6 g). Relaxing the symmetry assumption can be easily done by removing the second error term (last two lines) from Eq. (11). In this case, the optimal weight simplifies to
As before, we then replace the unknown ground-truth values with estimates based on the observed sample statistics, using either normality assumption or Curto’s correction, to evaluate the membership function according to Eq. (13).
5 Application to Image-Space Denoising
This section describes how to apply our general statistical framework to image-space denoising. In image space, the indices i and j each denote an individual pixel (not image-space coordinates). In particular, our implementation contains three major components. First, we track online statistics of the samples from the MC renderer, Alg. 1 ; this part is implemented as part of the renderer itself (we use pbrt-v3 [Pharr et al. 2016] for our results). Second, we select a joint bilateral filter to define the a priori weight ρij in Eq. (10). Finally, we implement our statistical denoising pipeline, which takes the noisy image, sample statistics, and a priori weights as input to compute the filter weights, Alg. 2, and produce the final output image (we implement this part on the GPU using CUDA).
Fig. 2:
Fig. 2: Bathroom scene rendered with 8192 SPP. Denoising time for Moon CI: 39.7 ms, OptiX: 91.2 ms, OIDN: 20.2 ms, ProDen: 2027.5 ms, ours: 29.3 ms.
Fig. 3:
Fig. 3: Salle de Bain scene rendered with 32 SPP. Denoising time for Moon CI: 36.4 ms, OptiX: 141.6 ms, OIDN: 22.1 ms, ProDen: 1979.9 ms, ours: 28.5 ms. In the supplementary document, we provide results for this scene at a higher sample count of 512 SPP.
Fig. 4:
Fig. 4: Classroom scene rendered with 256 SPP. Denoising time for Moon CI: 35.6 ms, OptiX: 89.5 ms, OIDN: 20.8 ms, ProDen: 1840.5 ms, ours: 28.3 ms.
Fig. 5:
Fig. 5: Convergence plots for our denoising method, compared to the approach by Moon et al. [2014] (“Moon CI”), NVIDIA OptiX AI-Accelerated Denoiser (“OptiX”), Intel Open Image Denoise (OIDN), and progressive denoising [Firmino et al. 2022] (“ProDen”). Error metrics top-to-bottom: root mean square error (RMSE), mean absolute error (MAE), and structural dissimilarity (DSSIM) (evaluated after tone mapping to low-dynamic-range (LDR) images, because DSSIM requires a bound value range). Scenes left-to-right: Figs. 1–4.
In MC rendering, we work with (weighted) radiance samples that cannot be negative, while a few paths may result in large contributions to a pixel’s intensity. In statistical terms, the distribution of samples tends to be right-skewed. Given these observations, we employ the widely used Box-Cox transformation [Box and Cox 2018] to “normalize” the rendered samples:
We note that for samples generated by MC rendering, choosing λ = 0 is impractical, as many samples may be zero (e.g., paths that do not reach a light source before termination), where the log function is undefined. In our experiments, we find that λ = 1/2 yields good results in practice, as it effectively “compresses” high-valued outliers while avoiding excessive “stretching” of small values toward − ∞. For each transformed sample \(x^{\prime }_k\) arriving from the renderer, we then update the online statistics following Alg. 1 .
As a base filter, we choose the joint bilateral filter [Eisemann and Durand 2004; Petschnigg et al. 2004] to define the a priori weights ρij for any pair of pixels (i, j). These weights follow a Gaussian falloff with distance in combined image space and G-buffers:
where pi denotes the a priori information for pixel i. We use image-space position, RGB albedo color, and surface normal for each pixel, i.e., \(\mathbf {p}_i = (x_i,y_i, r_i,g_i,b_i, n_{ix},n_{iy},n_{iz})^\textsf {T}\), in our results. Moreover, Σ denotes a covariance matrix that controls the rate of falloff for each dimension of pi. To enforce the sparsity of these weights, we limit the filter to a small neighborhood within a constant radius around each pixel in image space; our default covariance matrix is then Σ = diag(10, 10, 0.02, 0.02, 0.02, 0.1, 0.1, 0.1).
In practice, our statistical filter operates with three color channels per pixel in RGB color space and evaluates the membership function for each channel according to Alg. 2 . We set the final weight wij to zero unless all three channels pass the statistical test before applying the filter. In this way, we avoid color shifts that could occur if a channel is evaluated differently from the others.
In summary, our approach builds on three key insights: First, minimizing pair-wise MSE is closely related to Welch’s t-test; this theoretical contribution forms the basis of our method. Second, all collected statistics are noisy estimates; enforcing binary membership functions prevents most of this noise from propagating downstream. Finally, Box-Cox transformation, as well as correcting the mean [Curto 2023; Johnson 1978], makes our method more robust to non-normality, effectively mitigating the influence of outliers (“fireflies”).
The parameters of our method affect the results as follows: The threshold γ, Eq. (13), adjusts denoising strength. For symmetric weights, γ = 0 reverts to the base filter, γ = 1/2 and non-zero variance yields mij = 0 (when i ≠ j), effectively disabling filtering. The covariance matrix in Eq. (16) specifies the extent of the filtering window, depending on distances in image space and differences in G-buffer values; higher values generally increase the number of membership evaluations and lead to smoother results. The Box-Cox parameter λ determines how samples are transformed: λ = 1 keeps their distribution shape unchanged, smaller values correct right-skewed distributions (reducing positive outliers), larger values correct left skewness.
6 Results I — Image Denoising
In this section, we compare our statistical denoiser to state-of-the-art machine-learning-based denoisers (NVIDIA OptiX AI-Accelerated Denoiser and Intel Open Image Denoise), as well as a meta-denoiser [Firmino et al. 2022], on multiple well-known test scenes [Bitterli 2016] in Figs. 1–4. Figure 5 shows a quantitative convergence analysis for these examples. Note that without useful information from the G-buffers (e.g., objects seen in the mirror in Fig. 2), neural methods often fail to reconstruct sharp features. We also compare our membership formulation to the one-sided confidence-interval test proposed by Moon et al. [2013], using the \(99.8\%\) confidence level as proposed in their paper. Here, we compare the membership functions directly without applying their two-step denoising process.
Fig. 6:
Fig. 6: Ablation study: we compare our default settings (center) to (a) the noisy input, (b) the joint bilateral base filter only (i.e., mij = 1 ∀i, j), (c) a Gaussian base filter (i.e., no G-buffers, but using our membership functions), (d) using a smaller base filter (6-pixel radius; default is 20 pixels), (e) a more discerning significance level (α = 0.1 instead of the default α = 0.005), (f) disabling the Box-Cox transformation, (g) allowing asymmetric weights, Eq. (14), and (h) using the one-sided confidence-interval test by Moon et al. [2013]. Errors are given as RMSE/MAE.
We also analyze the effect of various components and parameters of our method in Fig. 6. Our statistical membership function greatly contributes to preserving the hard edges of shadows, which would be over-blurred by the base filter, as such lighting effects are not present in the G-buffers. Note how the G-buffers and our membership functions complement each other in reconstructing image features (Figs. 6 b and c). The Box-Cox transformation further increases performance by “normalizing” the sample distributions.
All results were rendered and denoised on a desktop PC with an AMD Ryzen 9 5950X CPU and an NVIDIA RTX 3080 Ti GPU. Our filtering runtimes are generally comparable to the duration of the inference step of Intel Open Image Denoise (OIDN) and substantially faster than progressive denoising [Firmino et al. 2022]. We set the base-filter radius to a relatively large size of 20 pixels (mostly for 1280 × 720 resolution images). Reducing the radius would further speed up our method, as the number of membership tests scales quadratically with the radius: while for 20 pixels, we require around 28 ms, for 6 pixels, we only require around 11 ms (Fig. 6 d). Note that using the base filter alone would lead to noticeable over-blurring (Fig. 6 b).
7 Results II — Further Applications
The generality of our denoising framework allows us to go beyond simply denoising an image produced by MC rendering. In particular, we see great potential in leveraging the ability to quickly denoise arbitrary statistics of arbitrary quantities tracked during path tracing. This section presents two such applications: including per-bounce approximate path contributions in Russian roulette (RR), as well as an extension to multiple importance sampling (MIS) using MIS strategy win rates.
7.1 Approximate-Contribution Russian Roulette
Here, we extend classic throughput-based RR for path tracing. Instead of determining the path termination probability based on just the camera-centric throughput, we track and denoise the average incoming radiance for a particular bounce index per pixel. The termination probability at each bounce then considers the product of the path throughput and this average incoming radiance. Our goal is to improve termination probabilities for each path individually; thus, we divide the product by the average incoming radiance at the corresponding pixel of the camera.
Formally, the per-bounce incident radiance is defined as
where Ξ is world space and pj, k(x) the probability of sampling a position x within that space for a given pixel j and bounce index k. In practice, we do not need to explicitly evaluate pj, k(x) as this probability is inherent to the path-tracing process. In fact, to estimate \(L^{\prime }_{j,k}\), we just have to average all incoming-radiance samples (over all paths) for a given pixel and bounce index, which can be easily integrated into common path-tracing code.
We then apply an iterative process during rendering: In each iteration, we track new incoming-radiance samples generated during path tracing, progressively refining our estimate of \(L^{\prime }_{j,k}\), Eq. (17), and the corresponding statistics required for denoising. At the end of each iteration, we denoise these estimates as described in §5. From the second iteration onward, where denoised estimates are available, we use modified path termination probabilities according to these radiance estimates.
Fig. 7:
Fig. 7: Comparison between our approximate-contribution Russian roulette (ACRR) and classic throughput-based RR starting at the first bounce (RR-1) and at the fifth bounce (RR-5). Images show RR-1 at 2048 SPP (417.53 s), RR-5 at 1451 SPP (418.05 s) and ACRR at 987 SPP (417.00 s). Timings include GPU upload, denoising, and download. Here, we also show relative mean squared error (rMSE) [Rousselle et al. 2011], as it allows to gauge the sampling performance for individual pixels, regardless of their absolute value.
Fig. 8:
Fig. 8: Comparison of our selective multiple importance sampling (SMIS) to standard multiple importance sampling [Veach 1997] (“Ref. MIS”) and visualization of first-bounce SMIS win rates for BSDF sampling (\(\eta _\text{BSDF}\)) and direct-light sampling (\(\eta _\text{DL}\)) before and after denoising. We use rMSE [Rousselle et al. 2011] to better assess the sampling performance for individual pixels. Images show Ref. MIS at 1024 SPP (734.80 s) and our SMIS at 1063 SPP (732.89 s).
Figure 7 shows improved image errors for our approximate-contribution version (ACRR) in an equal-time comparison over two variants of classic RR taking only path throughput into account. Classic RR often terminates only from a certain (e.g., fifth) bounce onward to avoid prematurely terminating high-contribution paths. Because of our improved path-contribution approximation, ACRR does not rely on this work-around. Rath et al. [2022] show that optimal weights for RR can be approximated using a spatio-directional data structure. While our method cannot achieve their sampling efficiency, we illustrate how our denoising framework, based on regular image buffers, can be used as an alternative. We provide further details and results in the supplementary document.
7.2 Selective Multiple Importance Sampling
As our second example application, we present an extension to (multi-sample) multiple importance sampling (MIS) [Veach 1997; Veach and Guibas 1995]. The basic idea of MIS is to combine multiple sampling strategies by a weighted average, depending on the probability density function (PDF) of each strategy. The sampling strategies typically cover different features of an integrand, for instance, via bidirectional-scattering-distribution-function (BSDF) or direct-light sampling. However, in situations where one sampling strategy is clearly better than the other(s), MIS can introduce additional variance (noise), especially when using the balance heuristic. Veach [1997] also developed the cut-off, power, and maximum heuristics to mitigate this issue, generally moving large weights closer to one and small weights closer to zero. Choosing MIS weights and sample allocations is non-trivial and an active area of research [Grittmann et al. 2019; Kondapaneni et al. 2019; Szirmay-Kalos and Sbert 2022], and the effectiveness of a strategy depends on the given situation. We illustrate how our framework can be used to selectively perform MIS only in beneficial cases, independently of the chosen strategy.
For this, we extend MIS by identifying strategies that should be disabled in order to save computational effort and reduce noise. The primary metric for this decision is the win rate, \(\eta _m = n_m^* / n_i\), of sampling strategy m, where a “win” (counted in \(n_m^*\)) is a non-zero sample whose PDF value of strategy m exceeds the PDF values of all other strategies. When determining this win rate—as in other MIS approaches—we evaluate all sampling strategies for each sample (not just the strategy that produced the sample), even if a strategy has already been disabled earlier. Similarly to the average radiance in the RR example described above, we estimate and denoise the win rates per pixel and per bounce. In contrast to radiance samples, we do not Box-Cox transform win rates. We again proceed iteratively, starting with all MIS strategies enabled, tracking and denoising win rates, and—from the second iteration onward—disabling inferior (ηm < 10− 3) sampling strategies (per pixel and per bounce). Finally, note that disabling a sampling strategy generally does not introduce bias, because two unbiased estimates (a MIS estimate and an estimate from a single sampling strategy) can be naturally combined without introducing bias. Figure 8 shows example results and equal-time convergence behavior for combining bidirectional-scattering-distribution-function and direct-light sampling.
8 Discussion and Conclusion
We have presented a simple yet effective denoising method for Monte Carlo rendering. Using well-known image filters for variance reduction and statistical tests (membership functions) to prevent bias, we achieve state-of-the-art image quality. Our approach is entirely free of pretrained components, using descriptive statistics of the sample distributions instead that can be estimated online during rendering. The properties of our membership function guarantee convergence with increasing sample counts. At low sample counts, high variance dominates the image error. In this case, the membership function becomes less discriminating, thereby reducing that variance and accepting some bias. As rendering progresses and more samples are added, variance decreases and the bias becomes more relevant to the overall error. If bias is significant (relative to variance), we set the corresponding filter weight to zero, based on statistical tests, which essentially eliminates that source of bias. In the limit, as variance approaches zero, any bias is unacceptable and thus prevented by the membership function.
We have focused on binary membership functions in this work. While we derive these functions from optimal pair-wise weights, recall that these weights are noisy estimates in practice. Introducing a threshold and restricting to binary memberships effectively prevents this noise from propagating through the pipeline. We leave investigating continuous membership functions for future work.
In contrast to neural-network-based approaches, our method does not require any computationally expensive training and does not risk adding “wrong” details that were present in the training data but should not be present in the output image.
Including our method as an additional input to existing meta-denoisers, such as progressive denoising [Firmino et al. 2022] or ensemble denoising [Zheng et al. 2021] could in turn deliver another step up in image quality for these methods. Furthermore, we also see great potential for the applicability of our statistics-based approach to other types of estimators during MC rendering. In particular, variance estimates have been used for adaptive sampling [Rousselle et al. 2012], MIS [Grittmann et al. 2019], or path guiding [Rath et al. 2020] in the past. Efficiently denoising variance estimates with our framework could yield improved performance of these methods in the future.
Another promising avenue for future work is the extension of our method to the temporal domain: Conceptually, estimates at different points in time of an animation can be treated equivalently to estimates at different image-space positions. Including such estimates could improve denoising performance and temporal coherence for animation denoising.
Acknowledgments
We would like to thank Lukas Lipp for fruitful discussions, Károly Zsolnai-Fehér and Jaroslav Křivánek for valuable contributions to early versions of this work, and Bernhard Kerbl for help with our CUDA implementation. Moreover, we would like to thank the creators of the scenes we have used: Wig42 for “Wooden Staircase” (Fig. 1), “Grey and White Room” (Fig. S6), and “Modern Living Room” (Fig. S8); nacimus for “Bathroom” (Fig. 3, S5); NovaZeeke for “Japanese Classroom” (Fig. 4, 6); Beeple for “Zero-Day” (Fig. 8); Jay-Artist for “White Room” (Fig. S7); Mareck for “Contemporary Bathroom” (Fig. 2); Christian Freude for “Glass Caustics” (Fig. S10); and Benedikt Bitterli for “Veach Ajar” (Fig. 7, S2), “Veach MIS” (Fig. S4), and “Fur Ball” (Fig. S11). This work has received funding from the Vienna Science and Technology Fund (WWTF) project ICT22-028 (“Toward Optimal Path Guiding for Photorealistic Rendering”) and the Austrian Science Fund (FWF) project F 77 (SFB “Advanced Computational Design”).
Jonghee Back, Binh-Son Hua, Toshiya Hachisuka, and Bochang Moon. 2022. Self-Supervised Post-Correction for Monte Carlo Denoising. In SIGGRAPH ’22: Special Interest Group on Computer Graphics and Interactive Techniques Conference, Vancouver, BC, Canada, August 7 - 11, 2022, Munkhtsetseg Nandigjav, Niloy J. Mitra, and Aaron Hertzmann (Eds.). ACM, 18:1–18:8.
Jonghee Back, Binh-Son Hua, Toshiya Hachisuka, and Bochang Moon. 2023. Input-Dependent Uncorrelated Weighting for Monte Carlo Denoising. In SIGGRAPH Asia 2023 Conference Papers, SA 2023, Sydney, NSW, Australia, December 12-15, 2023, June Kim, Ming C. Lin, and Bernd Bickel (Eds.). ACM, 9:1–9:10.
Pablo Bauszat, Martin Eisemann, and Marcus A. Magnor. 2011. Guided Image Filtering for Interactive High-quality Global Illumination. Comput. Graph. Forum 30, 4 (2011), 1361–1368.
Benedikt Bitterli, Fabrice Rousselle, Bochang Moon, José Antonio Iglesias Guitián, David Adler, Kenny Mitchell, Wojciech Jarosz, and Jan Novák. 2016. Nonlinearly Weighted First-order Regression for Denoising Monte Carlo Renderings. Comput. Graph. Forum 35, 4 (2016), 107–117.
Chakravarty R. Alla Chaitanya, Anton S. Kaplanyan, Christoph Schied, Marco Salvi, Aaron E. Lefohn, Derek Nowrouzezahrai, and Timo Aila. 2017. Interactive reconstruction of Monte Carlo image sequences using a recurrent denoising autoencoder. ACM Trans. Graph. 36, 4 (2017), 98:1–98:12.
José Dias Curto. 2023. Confidence intervals for means and variances of nonnormal distributions. Communications in Statistics - Simulation and Computation 52, 9 (2023), 4414–4430.
Holger Dammertz, Daniel Sewtz, Johannes Hanika, and Hendrik P. A. Lensch. 2010. Edge-avoiding À-Trous wavelet transform for fast global illumination filtering. In Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on High Performance Graphics 2010, Saarbrücken, Germany, June 25-27, 2010, Justin Hensley, Philipp Slusallek, David K. McAllister, and Christiaan P. Gribble (Eds.). Eurographics Association, 67–75.
Mauricio Delbracio, Pablo Musé, Antoni Buades, Julien Chauvier, Nicholas Phelps, and Jean-Michel Morel. 2014. Boosting monte carlo rendering by ray histogram fusion. ACM Trans. Graph. 33, 1 (2014), 8:1–8:15.
Arthur Firmino, Jeppe Revall Frisvad, and Henrik Wann Jensen. 2022. Progressive Denoising of Monte Carlo Rendered Images. Comput. Graph. Forum 41, 2 (2022), 1–11.
Norman J. Johnson. 1978. Modified t Tests and Confidence Intervals for Asymmetrical Populations. J. Amer. Statist. Assoc. 73, 363 (1978), 536–544. http://www.jstor.org/stable/2286597
James T. Kajiya. 1986. The rendering equation. In Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH 1986, Dallas, Texas, USA, August 18-22, 1986, David C. Evans and Russell J. Athay (Eds.). ACM, 143–150.
Nima Khademi Kalantari and Pradeep Sen. 2013. Removing the Noise in Monte Carlo Rendering with General Image Denoising Algorithms. Comput. Graph. Forum 32, 2 (2013), 93–102.
Ivo Kondapaneni, Petr Vévoda, Pascal Grittmann, Tomás Skrivan, Philipp Slusallek, and Jaroslav Krivánek. 2019. Optimal multiple importance sampling. ACM Trans. Graph. 38, 4 (2019), 37:1–37:14.
Yu Liu, Changwen Zheng, and Hongliang Yuan. 2018a. Denoising Monte Carlo Renderings based on a Robust High-order Function. In Proceedings of the 13th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2018) - Volume 1: GRAPP, Funchal, Madeira, Portugal, January 27-29, 2018, Dominique Bechmann, Ana Paula Cláudio, and José Braz (Eds.). SciTePress, 288–294.
Yu Liu, Changwen Zheng, Quan Zheng, and Hongliang Yuan. 2018b. Removing Monte Carlo noise using a Sobel operator and a guided image filter. Vis. Comput. 34, 4 (2018), 589–601.
Michael Mara, Morgan McGuire, Benedikt Bitterli, and Wojciech Jarosz. 2017. An efficient denoising algorithm for global illumination. In Proceedings of High Performance Graphics, HPG 2017, Los Angeles, CA, USA, July 28 - 30, 2017. ACM, 3:1–3:7.
Bochang Moon, José Antonio Iglesias Guitián, Sung-Eui Yoon, and Kenny Mitchell. 2015. Adaptive rendering with linear predictions. ACM Trans. Graph. 34, 4 (2015), 121:1–121:11.
Bochang Moon, Jong Yun Jun, JongHyeob Lee, Kunho Kim, Toshiya Hachisuka, and Sung-Eui Yoon. 2013. Robust Image Denoising Using a Virtual Flash Image for Monte Carlo Ray Tracing. Comput. Graph. Forum 32, 1 (2013), 139–151.
Hyosub Park, Bochang Moon, Soomin Kim, and Sung-Eui Yoon. 2013. P-RPF: Pixel-Based Random Parameter Filtering for Monte Carlo Rendering. In 2013 International Conference on Computer-Aided Design and Computer Graphics, CAD/Graphics 2013, Guangzhou, China, November 16-18, 2013. IEEE, 123–130.
Georg Petschnigg, Richard Szeliski, Maneesh Agrawala, Michael F. Cohen, Hugues Hoppe, and Kentaro Toyama. 2004. Digital photography with flash and no-flash image pairs. ACM Trans. Graph. 23, 3 (2004), 664–672.
Matt Pharr, Wenzel Jakob, and Greg Humphreys. 2016. Physically Based Rendering: From Theory to Implementation (3rd ed.) (3rd ed.). Morgan Kaufmann Publishers Inc., San Francisco, CA, USA. 1266 pages.
Alexander Rath, Pascal Grittmann, Sebastian Herholz, Petr Vévoda, Philipp Slusallek, and Jaroslav Krivánek. 2020. Variance-aware path guiding. ACM Trans. Graph. 39, 4 (2020), 151.
Alexander Rath, Pascal Grittmann, Sebastian Herholz, Philippe Weier, and Philipp Slusallek. 2022. EARS: efficiency-aware russian roulette and splitting. ACM Trans. Graph. 41, 4 (2022), 81:1–81:14.
Fabrice Rousselle, Marco Manzi, and Matthias Zwicker. 2013. Robust Denoising using Feature and Color Information. Comput. Graph. Forum 32, 7 (2013), 121–130.
F. E. Satterthwaite. 1946. An Approximate Distribution of Estimates of Variance Components. Biometrics Bulletin 2, 6 (1946), 110–114. http://www.jstor.org/stable/3002019
Christoph Schied, Anton Kaplanyan, Chris Wyman, Anjul Patney, Chakravarty R. Alla Chaitanya, John Burgess, Shiqiu Liu, Carsten Dachsbacher, Aaron E. Lefohn, and Marco Salvi. 2017. Spatiotemporal variance-guided filtering: real-time reconstruction for path-traced global illumination. In Proceedings of High Performance Graphics, HPG 2017, Los Angeles, CA, USA, July 28 - 30, 2017. ACM, 2:1–2:12.
Christoph Schied, Christoph Peters, and Carsten Dachsbacher. 2018. Gradient Estimation for Real-time Adaptive Temporal Filtering. Proc. ACM Comput. Graph. Interact. Tech. 1, 2 (2018), 24:1–24:16.
Pradeep Sen and Soheil Darabi. 2012. On filtering the noise from the random parameters in Monte Carlo rendering. ACM Trans. Graph. 31, 3 (2012), 18:1–18:15.
László Szirmay-Kalos and Mateu Sbert. 2022. Robust Sample Budget Allocation for MIS. In 43rd Annual Conference of the European Association for Computer Graphics, Eurographics 2022 - Short Papers, Reims, France, April 25-29, 2022, Nuria Pelechano and David Vanderhaeghe (Eds.). Eurographics Association, 17–20.
Eric Veach and Leonidas J. Guibas. 1995. Optimally combining sampling techniques for Monte Carlo rendering. In Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH 1995, Los Angeles, CA, USA, August 6-11, 1995. ACM, 419–428. https://dl.acm.org/citation.cfm?id=218498
Delio Vicini, David Adler, Jan Novák, Fabrice Rousselle, and Brent Burley. 2019. Denoising Deep Monte Carlo Renderings. Comput. Graph. Forum 38, 1 (2019), 316–327.
Changwen Zheng and Yu Liu. 2018. Removing Monte Carlo Noise with Compressed Sensing and Feature Information. In Proceedings of the 13th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2018) - Volume 1: GRAPP, Funchal, Madeira, Portugal, January 27-29, 2018, Dominique Bechmann, Ana Paula Cláudio, and José Braz (Eds.). SciTePress, 145–153.
Shaokun Zheng, Fengshi Zheng, Kun Xu, and Ling-Qi Yan. 2021. Ensemble denoising for Monte Carlo renderings. ACM Trans. Graph. 40, 6 (2021), 274:1–274:17.
Matthias Zwicker, Wojciech Jarosz, Jaakko Lehtinen, Bochang Moon, Ravi Ramamoorthi, Fabrice Rousselle, Pradeep Sen, Cyril Soler, and Sung-Eui Yoon. 2015. Recent Advances in Adaptive Sampling and Reconstruction for Monte Carlo Rendering. Comput. Graph. Forum 34, 2 (2015), 667–681.
Unbiased Monte Carlo path tracing that is extensively used in realistic rendering produces undesirable noise, especially with low samples per pixel (spp). Recently, several methods have coped with this problem by importing unbiased noisy images and ...
We propose a Deep Residual Learning based method that consistently outperforms both the state-of-the-art handcrafted denoisers and learning-based methods for single-image Monte Carlo denoising. Unlike the indirect nature of existing learning-based ...
Denoising Monte Carlo rendering with a very low sample rate remains a major challenge in the photo-realistic rendering research. Many previous works, including regression-based and learning-based methods, have been explored to achieve better rendering ...
