Decorrelating ReSTIR Samplers via MCMC Mutations

The efficiency of rendering algorithms often hinges on their ability to effectively evaluate similar integrals by reusing samples across pixels [Ward et al. 1988; Jensen 1996; Veach and Guibas 1997; Keller 1997]. In real-time path tracing, sample reuse becomes more critical since tracing rays is computationally intensive even on high-end consumer GPUs [Kilgariff et al. 2018]. Moreover, while existing denoisers drastically improve image quality even at low sample counts [Chaitanya et al. 2017; Schied et al. 2017;, 2018; Kozlowski and Cheblokov 2021; NVIDIA 2022], they are unable to reconstruct features missing from their input samples. Thus, sample reuse across pixels is often the only means to improve sampling quality given limited computational budgets. Compared to methods that generate independent samples, reuse is also at times the only practical approach to render challenging scenes with caustics and tricky lighting [Hachisuka and Jensen 2009; Veach and Guibas 1997].

Recent sampling algorithms for real-time ray tracing achieve massive speedups in scenes with complex illumination by sharing samples spatially within an image and temporally across frames [Bitterli et al. 2020; Ouyang et al. 2021; Lin et al. 2021;, 2022]. These so-called ReSTIR¹ based techniques select N high-contribution samples from a larger streamed candidate pool of size M. They do so by reformulating resampled importance sampling (RIS) [Talbot et al. 2005] in terms of weighted reservoir sampling (WRS) [Chao 1982]. While RIS effectively produces samples in proportion to an arbitrary target function (e.g., the integrand of the rendering equation), WRS makes resampling efficient by reducing storage costs from \(O(M)\) to \(O(N)\). Repeated resampling across pixels then helps distribute important samples over several frames for estimation.

Though ReSTIR derives impressive efficiency gains from correlated sampling, the benefits of repeated resampling are not indefinite. When only a few high-contribution samples have been identified, iterative spatial reuse creates blotchy artifacts as several pixels reuse the same sample (Figure 12, top row). Such undersampling artifacts eventually fade away with temporal reuse over several frames, using a user-specified parameter to balance pixel error with correlations from sample reuse (Figure 4). Unfortunately, simply emphasizing error reduction via greater reuse adds lag under camera movement with dynamically changing lighting and geometry (Section 2.4), and introduces distracting low-frequency artifacts (Figures 2 and 3) akin to those in photon mapping [Hachisuka and Jensen 2009], Metropolis Light Transport (MLT) [Veach and Guibas 1997] and Virtual Point Light (VPL) methods [Dachsbacher et al. 2014].

As spatiotemporal correlations are difficult to quantify, resolving artifacts is challenging. For instance, popular denoisers that compute first- and second-order moments (e.g., Schied et al. [2017]) are less effective given imprecise variance estimates with correlated samples. For ReSTIR, trying to reduce such artifacts by increasing the reservoir size N is also ineffective, as resampling with replacement [Chao 1982] produces duplicate samples in the presence of strong correlations (see Wyman and Panteleev [2021, Figure 19]).

Inspired by work on Sequential Monte Carlo (SMC) [Doucet et al. 2001] and Population Monte Carlo (PMC) [Cappé et al. 2004], we demonstrate that interleaving MCMC mutations with reservoir resampling (Section 3) helps alleviate correlations and impoverishment, especially in scenes with glossy materials and difficult lighting. Unlike MLT where mutations drive information sharing across pixels, our mutations instead primarily mitigate artifacts caused by spatiotemporal reuse, with little-to-no visual impact in scenes where artifacts do not arise (Figure 13). Our approach highlights the complementary strengths of resampling and mutations for real-time rendering: resampling identifies samples with large contributions proportional to a pixel’s target distribution, while mutations diversify the resampled population by locally perturbing samples in proportion to the same target distribution. Furthermore, like Veach and Guibas [1997]’s bias elimination strategy for MLT, we show that resampling eliminates the need for any burn-in period with Metropolis–Hastings (MH) mutations [Metropolis et al. 1953; Hastings 1970] (Section 2.5, Appendix A). This drives considerable image quality improvements from even a single mutation per frame for each reservoir sample (Figures 1, 7, 9, and 11).

From an implementation perspective, our approach requires only simple additions to existing ReSTIR algorithms (see Algorithm 3)—we mutate reservoir samples using Metropolis–Hastings and an appropriate target function every frame after temporal reuse. This is immediately followed by an adjustment to each mutated sample’s contribution weight to maintain detailed balance and ensure unbiased estimation. Overall, our contributions include:

Though mutations can help produce images with better visual fidelity, we observe that, similar to blue-noise sampling [Mitchell 1987; Georgiev and Fajardo 2016; Heitz and Belcour 2019], they do not necessarily reduce error (Figure 9). This is true both in scenes with difficult-to-sample lighting (Figure 12, top row), as well as in scenes with ample lighting and diffuse materials where resampling easily finds important samples (Figure 13). Moreover, despite the effectiveness of mutations in scenes with strong correlations, we show in the supplemental document that with our current approach, even infinite mutations cannot eliminate correlations entirely.

We start with the key building blocks of our approach in the next section, and postpone discussion about related work to Section 6 for better context when comparing with our method.

The rendering equation [Kajiya 1986] gives the outgoing radiance \(L_{\text{out}}\) leaving a point y in the direction \(\omega\). Expressed as an integral over directions, it isHere \(L_{\text{e}}\) is the emitted radiance, \(L_{\text{in}}(y, \omega _i)\) is the incoming radiance from the direction \(\omega _i\), \(\rho (y, \omega , \omega _i)\) is the BSDF and \(\theta _i\) is the angle between \(\omega _i\) and the surface normal at y. Absent participating media, the incident radiance \(L_{\text{in}}\) is defined recursively as \(L_{\text{in}}(y, \omega _i) = L_{\text{out}}(t(y, \omega _i), -\omega _i)\); the function \(t(y, \omega _i)\) returns the point on the closest surface from y in direction \(\omega _i\). Integrating over the sphere of directions \(S^2\) then gives the total radiance scattered toward \(\omega\); the rendering equation can be estimated with Monte Carlo aswhere \(p(\omega _i)\) is the probability density function (PDF) with respect to solid angle used to sample incident directions \(\omega _i\).

As in Kajiya’s formulation, sometimes it is more convenient to reformulate Equation (1) over surfaces. To keep the discussion independent of the choice of formulation, we use \(\int _{\Omega } f(x)\ \text{d}x\) to generically represent the integral we want to evaluate with \(\Omega\) as its domain. This integral can likewise be estimated usingwhere \(x_i\) are independent random samples drawn from any source PDF p that is non-zero on the support of f. In rendering, one often draws samples proportional to individual terms of the rendering equation to reduce variance (e.g., the BSDF \(\rho\)). To perform even better importance sampling, ReSTIR instead uses RIS to draw samples approximately proportional to the product of multiple terms in the integrand (e.g., \(L_{\text{in}} \cdot \rho \cdot |\text{cos}\ \theta |\)).

We review RIS and generalized RIS next (Sections 2.1 and 2.2); Section 2.3 discusses a streaming RIS implementation via reservoir sampling. Section 2.4 then describes how correlations arise within ReSTIR due to resampling. Section 2.5 reviews the Metropolis–Hastings algorithm we use in Section 3 to resolve correlation artifacts.

Sample selection with RIS from a target distribution improves with larger populations of candidate samples. ReSTIR provides access to a sizable candidate pool for resampling through spatiotemporal reuse, enabling it to quickly identify high-contribution samples via RIS. However, at times ReSTIR extensively reuses a few samples over multiple frames due to imperfect importance sampling and suboptimal parameters, with no mechanism to easily diversify an existing population of high-contribution samples.

Inspired by Sequential and Population Monte Carlo techniques (Section 6), we interleave reservoir resampling with MCMC mutations to mitigate correlations and sample impoverishment caused by spatiotemporal reuse. Our key observation is that mutating reservoir samples with the same per-pixel target function as RIS helps to quickly decorrelate the resampled population, especially when it contains outliers. In Algorithm 3, we use Metropolis-Hastings to locally perturb temporal reservoir samples selected by Algorithm 2; interleaving with resampling then diversifies the samples ReSTIR shares between pixels. We discuss key aspects of our work next, starting with how to modify mutated samples’ contribution weights to guarantee unbiased results.

Modified contribution weights. A contribution weight \(\widehat{W}\) (Equation (8)) estimates the reciprocal value of the target PDF \(\hat{p}/\int _{\Omega } \hat{p}\) that a sample is approximately distributed according to. \(\widehat{W}\) is needed to compute resampling weights for combining reservoirs (Algorithm 2, lines 7 and 11) and to estimate per-pixel shading (Equation (9)).

Contribution weights are sample dependent. Thus, a sample that undergoes mutation cannot reuse the weight associated with its original state, i.e., a mutated sample’s contribution weight should provide an unbiased estimate for the sample’s reciprocal target PDF. Our key contribution is to show that the unbiased contribution weight for any mutated sample \(x^k\), from a Markov chain \(x^0, \ldots , x^k, \ldots\), can be computed via the relationEquation (16) does not depend on samples between \(x^0\) and \(x^k\) in the Markov chain and imposes no constraints on computing \(\widehat{W}(x^0)\), which can arise from prior resampling, runs of MH, or a mix of the two. This provides flexibility in where and when to mutate samples during ReSTIR (as long as mutations are confined to a given pixel).

One can get an intuitive feel for Equation (16) by substituting in the expression for \(\widehat{W}(x^0)\) from Equation (8):Notice that the estimated normalization factor for \(\hat{p}\), i.e., the sum of weights w, remains unchanged for both the initial and mutated samples \(x^0\) and \(x^k\). This normalization factor arises via RIS prior to performing mutations (e.g., Algorithm 2, lines 14–15). Meanwhile, MH treats the resampling weights as fixed, simply redistributing a reservoir’s sample population proportionally to the per-pixel target function \(\hat{p}\). Equation (16) then encodes any required correction to a sample’s contribution weight to account for the sample mutation—unlike temporal and spatial resampling which reuse samples across different pixels, this equation does not contain any MIS weight or shift mapping as there is no change of integration domain, with samples only mutated within \(\hat{p}\)’s support.

Start-up bias. Algorithm 3 does not require a burn-in period for mutations, even though the samples used to initialize MH are not distributed exactly according to \(\hat{p}\). This is because we use the unbiased contribution weights of mutated samples for subsequent steps in ReSTIR, including when computing shading and resampling weights for further reuse. This approach eliminates start-up bias completely for any mutated sample \(x^k\) and function f by ensuringAppendix A provides a formal proof. Note that avoiding start-up bias does not imply samples generated using MH are well-distributed according to \(\hat{p}\). However, since we initialize MH using reservoir samples that are already distributed roughly proportional to the target function from resampling, our method does not rely on MH to find important samples (see Figure 13)—rather it decorrelates and diversifies outlier samples by mutating them locally in proportion to \(\hat{p}\) and adjusting their estimated contribution weight accordingly.

When to perform mutations?. Temporal reservoirs often contain stale samples, as ReSTIR assigns higher relative importance to existing samples. We therefore mutate samples output by Algorithm 2 within each pixel (Figure 5), using the same per-pixel target function as RIS for the current frame. Mutating samples randomly after temporal resampling diversifies the inputs to spatial resampling, protecting against possibly escalating amounts of sample impoverishment caused by repeated reuse.

Applying Algorithm 3 within each pixel to mutate samples after the initial or spatial resampling steps in ReSTIR (Section 2.4) is possible but not required. Like mutations, initial resampling serves to rejuvenate the sample population every frame (by introducing new independent samples into the population). Samples from spatial resampling are stored for future reuse; mutating them proportional to the current target function would cause them to lag by one frame.

Finally, Algorithm 3 places no restrictions on MH iteration count. To improve runtime performance, one could adaptively specify mutation counts per pixel (including no mutations) using, for instance, local correlation estimates. We leave development of such heuristics to future work and use a fixed, user-specified number of iterations.

We perform mutations for both direct and indirect illumination in ReSTIR using Kelemen et al. [2002]’s primary sample space (PSS) parameterization. This conveniently allows applying mutations directly to random number sequences used to generate light-carrying paths, while constraining path vertices to remain on the scene manifold. Moreover, it simplifies the use of certain shift mappings in ReSTIR PT, e.g., the random replay shift [Lin et al. 2022, Section 7.2].

In this section, we represent samples with a path vertex notation \(\bar{\boldsymbol {\mathbf {x}}}= [\boldsymbol {\mathbf {x}}_0, \boldsymbol {\mathbf {x}}_1, \ldots , \boldsymbol {\mathbf {x}}_k] \in \Omega ^k(\mathcal {M})\), with \(\Omega ^k(\mathcal {M})\) the space of all paths of length k on the scene manifold \(\mathcal {M}\) (e.g., \(k = 2\) for direct lighting). Each path \(\bar{\boldsymbol {\mathbf {x}}}\) is uniquely determined² by a vector of random numbers \(\bar{\boldsymbol {\mathbf {u}}}= [u_0, u_1 \ldots ] \in [0, 1]^{O(k)}\). We use S to denote a shift mapping from a base path \(\bar{\boldsymbol {\mathbf {x}}}\) in one pixel to an offset path \(\bar{\boldsymbol {\mathbf {y}}}\) in another pixel, i.e., \(S([\boldsymbol {\mathbf {x}}_0, \boldsymbol {\mathbf {x}}_1, \ldots , \boldsymbol {\mathbf {x}}_k]) = [\boldsymbol {\mathbf {y}}_0, \boldsymbol {\mathbf {y}}_1, \ldots , \boldsymbol {\mathbf {y}}_k]\). Mutated paths and random numbers are represented using \(\bar{\boldsymbol {\mathbf {z}}}\) and \(\bar{\boldsymbol {\mathbf {v}}}\) (resp.).

We prototyped our method in the open-source Falcor rendering framework [Kallweit et al. 2022]. All results use a GeForce RTX 3090 GPU at 1920 \(\times\) 1080 resolution. Our direct lighting implementation uses the same settings as Bitterli et al. [2020], i.e., initial candidate samples \(M = 32\), spatial reuse radius of 30 pixels from the current pixel, and \(M_\text{cap} = 20\). For indirect illumination, we set \(M = 32\) and the spatial reuse radius to 20 similar to Lin et al. [2022], but use a longer temporal history with \(M_{\text{cap}} = 50\) (unless noted otherwise). Our supplementary video shows 1 spp results for all our scenes; Table 1 gives single frame timings.

As we show in Figure 1 and Figures 7–12, short-range correlation artifacts are noticeably reduced in scenes with glossy materials and difficult lighting with just 1–5 mutations; further mutations have diminishing returns in improving image quality (Figures 9 and 10). Mutation cost overhead is generally less than simply increasing sample count (Figure 7), and recent denoisers [NVIDIA 2017] provide considerably better results with our decorrelated samples (Figure 1). Figure 8 shows mutations greatly reduce sample impoverishment, with fewer reservoirs sharing the exact same sample realizations.

Compared to standard path tracing, ReSTIR is much faster at achieving equal-error via correlated sampling for real-time direct [Bitterli et al. 2020, Figure 8] and global illumination [Lin et al. 2022, Figure 13]. Mutations however provide only marginal improvements in mean squared error in ReSTIR samplers (Figures 9 and 13), without ever negatively impacting results. Akin to blue-noise dithering [Georgiev and Fajardo 2016; Heitz and Belcour 2019], our image quality improves despite errors having similar magnitudes. The reason is mutating within a pixel leaves the sum of resampling weight unchanged in Equation (17), and these weights ultimately control RIS estimator variance (Equation (9)). Mutations do slightly reduce variance, as they indirectly alter resampling weights of future samples thanks to spatiotemporal reuse of the new, more diverse sample population; the supplementary document has more details. In Figures 11 and 12 we also ablate \(M_{\text{cap}}\) values to show the greater leeway our approach offers for this parameter, allowing the use of larger values to trade noise for correlation.

Since ReSTIR often suffers from correlation artifacts, we quantify improvements in correlation by computing sample covariance between pixels, which naturally generalizes sample variance. This metric measures the joint variability of two random variables (e.g., whether the error in two pixels varies similarly). For pixels i and j in image I, sample covariance \(\text{Cov}(i, j)\) between i and j is given bywhere K is the number of images used to estimate covariance (we use \(K = 100\)), and \(\bar{I}\) is the average of K images. To capture the joint variability of a pixel with its local neighborhood, in our experiments, we average covariance estimates over boxes of a given radius centered at each pixel. We then further average over the entire image to get a single number. Figure 9 (bottom left) shows average radial covariance decreases with increasing spatial radius. This is expected as ReSTIR only reuses samples in local neighborhoods (so small-scale correlation artifacts are more pronounced); mutations reduce covariance in these short ranges.

Table 1 lists the reduction in average covariance (with pixel radius equal to 8) and the FLIP weighted median score [Andersson et al. 2020] we observe on our scenes—FLIP tends to be marginally more sensitive to short range correlations than mean squared error. As correlations are typically localized, the reduction is even larger for the image insets in our figures compared to the results in Table 1. Ineffective shift mappings in ReSTIR often result in increased correlations; mutations compensate for this shortcoming. For instance, mutations typically have fewer correlation artifacts to resolve with a hybrid shift in ReSTIR PT compared to, e.g., random replay (Figures 13 and 8 respectively), which highlights the benefit of using good shift mappings. In contrast, mutations provide greater covariance reduction in the scenes rendered with ReSTIR DI (Figure 7), where higher covariance stems from vertex reconnections failing to preserve path contributions for low roughness surfaces.

Why mutations help?. The supplemental document details why mutations reduce covariance, simplifying down to the following, somewhat unintuitive, phenomenon: without mutations, covariance between pixels i and j stems from the mismatch between the distributions of input samples and the target functions at i and j (Equation (16) in the supplemental). However, in the limit of infinite mutations, covariance is determined by each sample’s mismatch with its own pixel’s target function (Equation (10) in the supplemental) due to the ratio \(\hat{p}(x^0)/\hat{p}(x^k)\) in the mutated contribution weight (Equation (16)); this mismatch tends to be smaller. Though our analysis predicts that covariance does not vanish completely even with infinite mutations, our results show covariance is often reduced with just one mutation.

Our method builds directly on the recent ReSTIR family of algorithms for real-time direct [Bitterli et al. 2020] and global illumination [Ouyang et al. 2021; Lin et al. 2021;, 2022]. We augment spatiotemporal reservoir resampling in ReSTIR with sample mutations, and demonstrate the complementary strengths of resampling and mutations in this framework. In graphics, our approach is most closely related to Metropolis Light Transport (MLT) [Veach and Guibas 1997] and associated techniques [Kelemen et al. 2002; Jakob and Marschner 2012; Lehtinen et al. 2013; Hachisuka et al. 2014; Otsu et al. 2018; Cline et al. 2005; Lai et al. 2007;, 2009; Bashford-Rogers et al. 2021]. In the broader Monte Carlo landscape, our approach belongs to the class of algorithms that jointly use resampling and mutations for sampling problems, such as SMC [Doucet et al. 2001] and PMC [Cappé et al. 2004]. We discuss the relation to MLT, SMC and PMC in more detail next; Table 2 provides a summary. We refer the reader to Bitterli et al. [2020, Section 7] and Lin et al. [2022, Section 9.3] for comparisons between ReSTIR and other rendering algorithms that exploit path reuse and spatial correlations.

Metropolis Light Transport. MLT uses statistically correlated samples generated by Metropolis–Hastings to solve the rendering equation. Unlike algorithms using independent samples, MLT is effective at finding difficult light paths by locally exploring the path space. It reuses samples by mutating high-contribution paths over the image. Algorithmically, our method resembles MLT in various ways. Both techniques require secondary estimators, respectively RIS and bidirectional path tracing (BDPT) [Lafortune and Willems 1993; Veach and Guibas 1995a], to normalize the MH target function. Samples used by these estimators are resampled into a smaller set to initialize MH (our Section 3 and Veach [1998, Chapter 11.3.1]), and contributions of mutated samples are effectively weighted by the same weights (Equation (17)) to remain unbiased (our Appendix A and Veach [1998, Appendix 11. A]).

The crucial difference between our work and MLT lies in how samples are reused across pixels. MLT latches onto high-contribution paths and mutates them over the entire image while resampling just to eliminate start-up bias. Thus, MLT results often contain correlation artifacts caused by mutations, applying MH to both find important samples and redistribute them between pixels. In contrast, ReSTIR derives spatiotemporal reuse exclusively from resampling; in this article, we mutate samples within each pixel to mitigate correlations and sample impoverishment from spatiotemporal resampling. As a result, our method does not require numerous MH iterations, as the primary purpose of mutations is not finding important paths (Figures 12 and 13). Further, our approach suits real-time rendering as it integrates seamlessly into ReSTIR. MLT can be adapted to mutate temporally, but unlike our work, the entire animated sequence must be available in advance [Van de Woestijne et al. 2017].

Several features have recently been added to MLT, including sample stratification [Cline et al. 2005; Gruson et al. 2020], MIS [Hachisuka et al. 2014] and enhanced mutation strategies [Jakob and Marschner 2012; Kaplanyan et al. 2014; Bitterli et al. 2017; Otsu et al. 2018; Luan et al. 2020]. Though we mostly employ simple PSS mutations [Kelemen et al. 2002], many of these improvements can also be incorporated into our approach to reduce correlations faster.

Sequential Monte Car- lo. SMC is a family of Monte Carlo methods used for filtering and tracking problems in Bayesian inference and signal processing [Doucet et al. 2001]. As shown in the inset, the goal is maintaining a population of weighted samples distributed roughly proportional to an evolving target distribution (with unknown normalization factor). Sample weights are adjusted every iteration to reflect each sample’s importance to the most recent distribution. Resampling discards samples with low weights and duplicates those with high weights. Mutations ensure the population does not contain identical samples. Unlike ReSTIR, which uses RIS, SMC methods use weighted importance sampling (WIS) to estimate correlated integrals in a chained fashion, i.e., the current step’s sample weights and normalization factors are defined incrementally based on corresponding quantities from earlier steps [Del Moral et al. 2006]. This allows temporally reusing samples for estimation, instead of generating new samples every frame.

SMC methods have found limited use in rendering: Hedman et al. [2016] maintain a temporally coherent distribution of VPLs for indirect illumination via a sampling method that moves as few of the VPLs between frames as possible. This work is only loosely inspired by SMC, in that it does not perform any sample mutations and generates biased results as it discards VPLs heuristically.

Ghosh et al. [2006] instead sample a sequence of per-pixel target functions for direct illumination of dynamic environment maps. Unlike ReSTIR, which derives its samples from spatiotemporally neighboring pixels, Ghosh et al. [2006] maintain a fixed sample population per pixel that is resampled and mutated to be updated for each frame. Large populations are needed for effective importance sampling, as high-contribution samples are not shared spatially between pixels; in contrast, ReSTIR often stores just a single sample per reservoir. SMC methods likely require MIS weights and shift mappings (like ReSTIR) to resolve bias and correctly derive effective spatiotemporal reuse from neighbors. Similar to our work, mutations mitigate sample impoverishment but do not provide reuse.

Population Monte Carlo. PMC methods also couple resampling and mutations to distribute weighted samples in proportion to a sequence of target functions [Cappé et al. 2004]. The main added feature is they sample using parametric mixture models with simple source PDFs. Mixture probabilities are tuned for each target function using previously generated samples and their importance.

In rendering, the PMC framework has been used for direct lighting [Fan et al. 2007; Lai et al. 2015], global illumination [Lai and Dyer 2007; Lai et al. 2007], and animation [Lai et al. 2009]. Lai et al.’s [2009] work is most relevant to ours: they derive sample reuse by mutating samples spatially and temporally across the image plane using Energy Redistribution Path Tracing (ERPT) [Cline et al. 2005]. Resampling serves to select high-contribution samples while discarding those with small weights; it is also used to refresh the sample population (much like initial resampling in ReSTIR) and eliminate start-up bias from mutations. Unlike our method, they require knowing animated sequences in advance, precluding most real-time applications.

In this article, we provide an unbiased mechanism leveraging MCMC mutations to diversify ReSTIR’s sample population. Often, just a single mutation per pixel effectively mitigates correlation artifacts in glossy scenes with complex lighting. However, as in most MCMC schemes, we cannot accurately predict the number of Metropolis–Hastings iterations needed to reduce correlations below a given threshold. Beyond the analysis in the supplemental document, further investigation is also needed to understand how mutations address sample impoverishment in ReSTIR—not just in terms of the number of duplicate samples (Figure 8), but also the discrepancy characteristics of the resulting sample population.

Mutations in ReSTIR have a non-negligible run-time overhead. Though we demonstrate improvements on an equal-time covariance metric with simple mutation strategies in both ReSTIR DI and PT (Figure 7 and Table 1), more sophisticated mutations [Jakob and Marschner 2012; Kaplanyan et al. 2014; Bitterli et al. 2017; Otsu et al. 2018; Luan et al. 2020] could provide further gains. Our decision to mutate only after temporal (but not spatial) resampling is informed in part by run-time considerations. As mentioned in Section 3, applying mutations selectively (i.e., not at each pixel every frame) based on local correlation heuristics could improve performance. As both mutations and ReSTIR’s initial path candidates serve to rejuvenate the sample population, it may also be worthwhile to carefully balance costs of per-pixel mutations and new path candidates.

Our proposed sample mutations reduce the correlation between nearby pixels, leading to an error distribution (likely) closer to white noise. But blue noise error distributions are often superior with respect to human perception [Mitchell 1987]; perhaps our mutations could be modified to more directly optimize for blue noise characteristics. For example, when deciding mutation acceptance, we might consider both the target function and the neighboring pixel samples, preferring mutations that introduce differing sample values. A further improvement might apply the insights of Heitz and Belcour [2019] to optimize the image-space distribution of error rather than solely considering the sample values.

Like Metropolis Light Transport, mutating samples across pixels potentially unlocks further amortization by sharing samples over the entire image (e.g., using the expected values technique [Veach 1998, Section 11.5]). We leave such “cross-pixel” mutations to future work as the change of integration domain from one pixel to another requires a further adjustment to a mutated sample’s contribution weight beyond Equation (16), i.e., as with temporal and spatial resampling, MIS weights and shift mappings are likely needed to correctly mutate a sample in proportion to the target function of a neighboring pixel whose support does not coincide entirely with its own. Alternatively, it is possible that the replica exchange approach adopted by Gruson et al. [2020] enables swapping of independent Markov chains from separate pixels in our framework as well.

Our mathematical presentation of resampling for ReSTIR in Section 2 is essentially unchanged from that of Lin et al. [2021] for fast volume rendering. Consequently, extending their approach to incorporate mutations from Section 3 likely requires no further modifications to the framework for volumes beyond the specific choice of mutation strategies, such as the scattering and propagation perturbations proposed by Pauly et al. [2000, Section 5.1].

More generally, by augmenting ReSTIR with mutations, our work establishes a closer correspondence between the RIS-based resampling techniques developed in graphics, and those in the broader statistics literature such as SMC and PMC. In particular, our approach stands to benefit from techniques such as annealed importance sampling [Neal 2001] used in SMC to reduce variance in the resampling weights [Ghosh et al. 2006, Section 4], as well as from adaptation strategies for mutation kernels developed in PMC to increase acceptance rates [Lai et al. 2007, Section 4.2]. Moreover, as in these fields, mutations in ReSTIR open the door not just to artifact-free integration (of the rendering equation), but also to tracking and filtering problems—for instance using well-distributed sample populations generated by our approach as training data for path guiding [Müller et al. 2017; Müller et al. 2019].

The authors thank Aaron Lefhon, Bill Dally and Keenan Crane for supporting this work. This work was generously supported by an NVIDIA Graduate Fellowship for the first author during his graduate studies at Carnegie Mellon University.

Equation (16) shows how to update the contribution weight \(\widehat{W}(x^k)\) of a mutated sample \(x^k\) from the Markov chain \(x^0, \ldots , x^k, \ldots\), generated with target function \(\hat{p}\). Here we prove this rule yields an unbiased contribution weight for any mutated sample \(x^k\), i.e., for any f with the same or smaller support,We assume sample \(x^0\) initializing the chain has the same support \(\Omega\) as target function \(\hat{p}\). For us, this is guaranteed by chained applications of RIS with a valid shift map in Algorithm 2. Any \(x^0\) chosen by RIS is not distributed exactly proportional to \(\hat{p}\) (unless we have infinite samples), however, its contribution weight \(\widehat{W}(x^0)\) is unbiased and satisfies Equation (22) [Lin et al. 2022]. Next, we show access to \(\widehat{W}(x^0)\) is sufficient to eliminate any startup bias with MH.

Proof. To show Equation (22) holds, we first express the update rule for contribution weight \(\widehat{W}(x^k)\) (for \(k \gt 0\)) in terms of the previous sample \(x^{k-1}\) in the chain, as follows:This is equivalent to Equation (16), shown by recursively unfolding this relationship for all prior samples \(x^{k-1}\) to \(x^1\) in the chain. As in Section 2.5, we also assume a candidate mutation to \(x^{k-1}\) is generated using the proposal density \(T(x^{k-1} \rightarrow z^{k-1})\), with acceptance probability for candidate \(z^{k-1}\) given by Equation (14):Metropolis–Hastings sets \(x^k = z^{k-1}\) with probability a; otherwise \(x^k = x^{k-1}\). This lets us rewrite the expectation in Equation (22):Rearranging the terms slightly yieldsWe now write each expectation as an integral. First, assume an inductive hypothesis \(\mathbb {E}[f(x^{k-1}) \widehat{W}(x^{k-1})]\!=\! \int _{\Omega } f(x)\ \text{d}x\) for the \(k\!-\!1^{\textrm {st}}\) MH iteration. Base case \(k\!=\!1\) holds trivially, as \(\widehat{W}(x^0)\) is an unbiased contribution weight. Next, note that for any integrable function \(g(x^{k-1}, z^{k-1})\), its expectation \(\mathbb {E}[g(x^{k-1}, z^{k-1}) \widehat{W}(x^{k-1})]\) can be rewritten as a conditional expectation over candidate mutations:where T is the proposal density used for mutations. This lets us expand out Equation (26) as follows:Finally, we show that the two double integrals cancel each other out (resulting in \(\mathbb {E}[f(x^k) \widehat{W}(x^k)] = \int _{\Omega } f(x)\ \text{d}x\)) by invoking the detailed balance condition from Equation (15):and rewriting it asSubstituting for \(T(x \rightarrow z) a(x \rightarrow z)\) in the third line of Equation (28) yields the same integral as the second line, but with integration variables x and z swapped. Renaming x and z and swapping the integration order in the third line allows cancellation, simplifying to \(\int _{\Omega } f(x)\ \text{d}x\), yielding Equation (22), and giving a proof by induction.

ReSTIR uses MIS weights during resampling (Equations (12) and (13)) to mitigate noise and bias from reusing samples across pixels. We provide explicit expressions for the MIS weights used in Algorithm 2 here; Lin et al. [2022, Equations (37)–(38)] provide similar expressions for Pairwise MIS weights needed for spatial resampling.

Let \(S_j\!:\!\Omega _j\!\rightarrow \!\Omega _i\) denote the shift map from pixel j to pixel i. Let \(x_i\) and \(x_j\) further represent the corresponding samples for these pixels, and \(S_j(x_j)\!=\!x_j^{\prime }\) and \(S_j^{-1}(x_i)\!=\!x_i^{\prime }\) the respective shift mapped values. The MIS weights for \(x_i\) and \(x_j^{\prime }\) are then given byWe set \(m_i(x_i) = 1\) and \(m_j(x_j^{\prime }) = 0\) when valid shifts do not exist for \(x_i\) and \(x_j\) (resp.); Lin et al. [2022, Section 5.6] discusses properties of these MIS weights in detail.

A mutation involving a reconnection vertex \(\boldsymbol {\mathbf {y}}_i\) requires modifying random numbers not just for \(\boldsymbol {\mathbf {y}}_i\), but also for the non-mutated vertices \(\boldsymbol {\mathbf {y}}_{i+1}\) and \(\boldsymbol {\mathbf {y}}_{i+2}\). This is because the solid angle PDFs used to sample outgoing directions \(\nu _{i}\) and \(\omega _{i+1}\) in Section 4.3 depend on the mutated incoming directions \(\nu _{i-1}\) and \(\nu _{i}\) (resp.). Here we derive Equation (20) by first noting the joint PDF for connecting the mutated reconnection vertex \(\boldsymbol {\mathbf {z}}_i\) to \(\boldsymbol {\mathbf {z}}_{i+1}\) and \(\boldsymbol {\mathbf {z}}_{i+1}\) to \(\boldsymbol {\mathbf {z}}_{i+2}\) in the surface area measure is:This is a product of delta functions as \(\boldsymbol {\mathbf {z}}_{i+1}=\boldsymbol {\mathbf {y}}_{i+1}\) and \(\boldsymbol {\mathbf {z}}_{i+2}=\boldsymbol {\mathbf {y}}_{i+2}\) are the only valid vertex positions. In the PSS to path space mapping, the joint PDF for the mutated random numbers \(\bar{\boldsymbol {\mathbf {v}}}_{i+1}\) and \(\bar{\boldsymbol {\mathbf {v}}}_{i+2}\) (for vertices \(\boldsymbol {\mathbf {y}}_{i+1}\) and \(\boldsymbol {\mathbf {y}}_{i+2}\)) is related to \(p(\boldsymbol {\mathbf {z}}_{i+2}, \boldsymbol {\mathbf {z}}_{i+1}\ |\ [\boldsymbol {\mathbf {z}}_0, \boldsymbol {\mathbf {z}}_1, \ldots , \boldsymbol {\mathbf {z}}_i])\) via a Jacobian determinant:This PDF serves as our proposal density \(T(\bar{\boldsymbol {\mathbf {u}}}\rightarrow \bar{\boldsymbol {\mathbf {v}}})\) for mutations, which then yields:The delta functions in the first line cancel since they are symmetric. In the third line, we use the fact that the Jacobian determinant of a sampling scheme is the same as its inverse PDF [Kelemen et al. 2002, Section 2]. The final step substitutes in the definition of the Jacobians relating the solid angle and area measures.