Ensemble metropolis light transport

5 Ensemble Metropolis Light Transport THOMAS BASHFORD-ROGERS, University of the West of England, UK LUÍS PAULO SANTOS, Universidade do Minho/INESC TEC, Portugal DEMETRIS MARNERIDES and KURT DEBATTISTA, University of Warwick, UK Fig. 1. The DOOR AJAR scene showing the reference on the left, Metropolis Light Transport in the centre, and Ensemble Metropolis Light Transport on the right, rendered at an equal sample count. The use of the ensemble to propose anisotropic transition kernels allows sampling to be adapted to the scene geometry and lighting information, leading to variance reduction as can be seen by lower MSE values shown at the top of each image. This article proposes a Markov Chain Monte Carlo (MCMC)rendering algorithm based on a family of guided transition kernels. The kernels exploit properties of ensembles of light transport paths, which are distributed according to the lighting in the scene, and utilize this information to make informed decisions for guiding local path sampling. Critically, our approach does not require caching distributions in world space, saving time and memory, yet it is able to make guided sampling decisions based on whole paths. We show how this can be implemented efficiently by organizing the paths in each ensemble and designing transition kernels for MCMC rendering based on a carefully chosen subset of paths from the ensemble. This algorithm is easy to parallelize and leads to improvements in variance when rendering a variety of scenes. CCS Concepts: • Computing methodologies →Ray tracing; Additional Key Words and Phrases: Light transport, MCMC, ensemble Authors’ addresses: T. B.-Rogers, University of the West of England, Coldharbour Ln, Bristol BS16 1QY, UK; L. P. Santos, Universidade do Minho/INESC TEC, Campus de Gualtar, Braga 4710-057, Portugal; D. Marnerides and K. Debattista, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK. Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from [email protected]. © 2021 Copyright held by the owner/author(s). Publication rights licensed to ACM. 0730-0301/2021/12-ART5 $15.00 https://doi.org/10.1145/3472294 ACM Reference format: Thomas Bashford-Rogers, Luís Paulo Santos, Demetris Marnerides, and Kurt Debattista. 2021. Ensemble Metropolis Light Transport. ACM Trans. Graph. 41, 1, Article 5 (December 2021), 15 pages. https://doi.org/10.1145/3472294 1 INTRODUCTION Accurately rendering photorealistic imagery requires computing extremely large numbers of light paths in a virtual environment. While most research has been applied to traditional Monte Carlo estimators for rendering, Markov Chain Monte Carlo (MCMC) methods, such as Metropolis Light Transport (MLT) [Veach and Guibas 1997], have shown impressive capabilities for computing light transport efficiently, even in complicated scenes. MCMC algorithms such as MLT generate a chain of paths that follow the distribution of the lighting in the scene. Each new path is generated by applying a transition kernel to the previous path, and probabilistically replacing the previous path with the new path. In the original application of MCMC to graphics [Veach and Guibas 1997], transition kernels were designed to either locally explore regions around an existing path, or to globally explore path space. While these transition kernels have been improved to consider local geometric or lighting information [Li et al. 2015;Otsu et al. 2018], the use of non-local information capturing a wider range of lighting can also be used to guide transition kernels. Such non-local information can be captured by path guiding methods, for example [Müller et al. 2017]; however, this comes at a precomputation and memory cost, and the cached distributions of lighting may not clearly map to transition kernels. Another approach is to exploit the fact that multiple light paths generated by MCMC ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. 5:2 • T. Bashford-Rogers et al. algorithms will be distributed proportional to the lighting in the scene, and as such can be used to generate guided transition kernels. This article proposes such a method that uses ensembles of light paths to guide mutations of existing paths. We name this approach Ensemble Metropolis Light Transport (EMLT). Crucially, these guided transition kernels do not need to be based on caching distributions in world space; they only require moderately sized ensembles in the low tens of thousands of paths, and can be combined to build families of mutation strategies. To summarize, the main contributions of this work are as follows: •The introduction of EMLT, a method that guides sampling based on a complementary ensemble of transport paths. •A family of adaptive and anisotropic proposal distributions for path mutations based on ensemble sampling. •The use of a carefully chosen subset of paths from the ensemble to create guided transition kernels. •Results showing improvement of EMLT over traditional approaches in a range of real-world scenes. 2 BACKGROUND AND RELATED WORK This section introduces the relevant background theory of light transport, path guiding methods which exploit information about the radiance or importance distribution in the scene to reduce variance, and MCMC methods which can efficiently compute images in challenging scenes. 2.1 Light Transport The path integral form of the rendering equation [Kajiya 1986]is given by [Hachisuka et al. 2014; Veach and Guibas 1997] Ij=P hj(x)f(x)dμ(x),(1) and states that the intensity Ijat a pixel jconsists of the contribution f(x)of light paths xweighted by a pixel filter hj(x).The domain of integration is the union of all possible path lengths P=∞ k=2P(k)where P(k)are all paths of length k. In this work, path vertices x0..xklie on the scene manifold M, i.e., integration is with respect to the product area measure μ, and are indexed starting from the light source x0. The contribution of a path of length kis defined as f(x)=Le(x0)G(x0↔x1)k−1 j=1fr(xj−1→xj→ xj+1)G(xj↔xj+1),whereLe is the emitted radiance, Gis the geometry term, and fr is the Bidirectional Reflectance Distribution Function (BRDF). There are multiple ways of solving Equation (1), almost all relying on Monte Carlo estimation: Ij≈1 N N  i=1 hj(x(i))f(x(i)) p(x(i)) ,(2) where p(x(i)) denotes the probability density function (pdf) of sampling the ith path and is a product of probability densities for sampling each vertex to build up the path. This typically consists of sampling the sensor, lens, BRDFs and light sources. Ideally p(x(i)) ∝hj(x(i))f(x(i)); however, this is typically not possible in practice. Therefore, distributions which approximate some components of hj(x(i))f(x(i)) are used (see Christensen and Jarosz [2016] for a survey of these methods). 2.2 Path Guiding Most rendering techniques generate light paths incrementally by sampling the next vertex in a path given the previous vertex. Path guiding approaches build on traditional BRDF and cosine sampling to include information about incoming illumination or importance when generating samples. Most techniques cache a distribution that represents the incoming radiance or importance at a sparse set of locations in a scene, and query locations at runtime using a spatial data structure. Examples include 5D spatio-directional Trees [Lafortune and Willems 1995; Müller et al. 2017], 7D distributions [Pantaleoni 2020], use of various basis functions to store radiance at discrete points in the scene [Bashford-Rogers et al. 2012; Diolatzis et al. 2020; Herholz et al. 2016; Hey and Purgathofer 2002; Jensen 1995; Ruppert et al. 2020; Vorba et al. 2014], or using machine learning methods [Bako et al. 2019; Dahm and Keller 2017]. Path guiding has also been applied to sample partial or complete light paths. Approaches such as neural importance sampling [Guo et al. 2018; Müller et al. 2018; Zheng and Zwicker 2019]have learned a warping in Primary Sample Space (PSS)[Kelemenetal. 2002] which encodes the illumination distribution in PSS based on a small set of paths traced before rendering. However, as these approaches are designed to generate full paths, they face the curse of dimensionality and are more effective in lower dimensional scenarios such as importance sampling one bounce indirect lighting. A related approach to whole path importance sampling was proposed by Reibold et al. [2018] that selectively stores and samples distributions for high contribution paths which were unlikely to be sampled through BRDF sampling. EMLT exploits information about the lighting distribution in the scene, and can use any of the distributions commonly used for path guiding to generate samples. However, our method does not require a spatial cache, and builds distributions on the fly from a small set of paths from a complementary ensemble (see Section 3). 2.3 MCMC MCMC [Hastings 1970; Metropolis et al. 1953] techniques provide another approach to sample a space and were initially applied to rendering as MLT [Veach and Guibas 1997]. In MLT, sampling starts from an initial path x, and then proposes a new path xfrom a transition kernel T(x→x). The transition kernels have to satisfy certain properties in order for the chain to explore the state space: ergodicity meaning all states will be visited by the chain in a finite time, and aperiodicity meaning states will not get stuck in a loop. At the limit, these states are distributed according to a target distribution, f b,whereb=Phj(x)f(x)dμ(x)is a normalization constant. A scalar contribution function, f∗:RS→ R, is defined where Sare the spectra or color channels associated with evaluating f(x). Then, based on the detailed balance condition f∗(x)T(x→x)a(x→x)=f∗(x)T(x→x)a(x→x), the new state xis probabilistically chosen to replace the previous ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. Ensemble Metropolis Light Transport • 5:3 state based on calculating an acceptance probability: a(x→x)=min  1, f∗(x)T(x→x) f∗(x)T(x→x) .(3) This leads to a chain of light paths, each dependent only on the previously sampled light path, which explore the path space. In rendering, the scalar contribution function is typically chosen to be the luminance of the contribution of the path, but other functions can be chosen (see Gruson et al. [2016]; Hoberock and Hart [2010]). bis typically estimated by a separate Monte Carlo estimator, and the initial state of the chain is generated through resampling a path from a small set of paths computed at startup. Then, the resulting Monte Carlo estimator is given by Ij≈b N N  i=1 hj(x(i))f(x(i)) f(x(i)) =b N N  i=1 hj(x(i)),(4) meaning that samples will be distributed according to the integrand. Equation (4) is typically evaluated over the image plane which allows for information about light transport to be shared between pixels, resulting in a significantly more efficient estimator. Many strategies can be designed such that many terms in the numerator and denominator in Equation (3) cancel, something which is especially important to remove the weak singularity in the geometry term. The variance reduction properties of this method also depends on the ability of the transition kernel to explore the state space. Veach and Guibas [1997] proposed a series of transition kernels which were chosen to reduce variance for different types of light transport. Bidirectional mutations were designed to ensure ergodicity through deleting a randomly chosen series of vertices from x, and replacing them with vertices generated through sampling the same pdfs used in a standard Monte Carlo estimator, i.e., BRDF and light source sampling. The remaining strategies, known as perturbations, were designed to explore sub spaces of path space given the state of the path. Lens perturbations explored image space by perturbing the position of the path vertex on the image plane, tracing a path through any specular interactions until a nonspecular vertex is reached, then deterministically connecting to the unchanged light subpath. This connection leaves a geometry term associated with the deterministic connection when evaluating Equation (3). The caustic perturbation was designed to explore caustics through perturbing the outgoing direction for the first vertex on the caustic subpath, following the chain of specular interactions, and then deterministically connecting to the camera. The multi-chain perturbation explores specular-diffuse-specular paths through combining a perturbation on the lens with directional perturbations at each non-specular surface before deterministically connecting to the remaining light subpath. There have been several extensions to the original MLT algorithm which have added or improved mutation strategies such as perturbations in participating media [Pauly et al. 2000], improved sampling of specular chains [Jakob and Marschner 2012], and perturbations in half vector space [Kaplanyan et al. 2014]. Kelemen et al. [2002] introduced mutations to light paths in PSS (PSSMLT). These mutations were improved by Hachisuka et al. [2014]who combined PSSMLT with Multiple Importance Sampling [Veach and Guibas 1995], Bitterli and Jarosz [2019] who detected and perturbed high variance paths in PSS, the use of delayed rejection by Rioux-Lavoie et al. [2020], the use of Hamiltonian Monte Carlo applied to rendering by Li et al. [2015] who used anisotropic Gaussian kernels generated from a path gradient, and Luan et al. [2020] who used the Metropolis-adjusted Langevin algorithm also based on the gradients of the path. Integration in both path space and PSS have been proposed [Bitterli et al. 2018;Otsuetal.2017; Pantaleoni 2017] which allows path space mutations to be combined with PSS mutations. For further information, Šik and Křivánek [2018] provide a detailed survey of MCMC methods in rendering. Closer to our work, adaptive perturbation sizes based on scene geometry were proposed by Otsu et al. [2018], which used cone tracing to estimate how large a perturbation could be based on the surrounding geometry of a path. This was applied starting at the camera, followed specular bounces if any, then traced one extra path vertex to form a perturbed path. Other methods mutate a set of paths but do not directly use these to adapt transition kernels. Energy Redistribution Path Tracing [Cline et al. 2005] combined Path Tracing and MLT by creating many short chains whenever a path would be better explored by MCMC methods than standard Monte Carlo. Segovia et al. [2007] used Multiple-Try MCMC to generate paths for Instant Radiosity [Keller 1997], and Nimier-David et al. [2019]alsoproposeda Multiple-Try MCMC method suitable for vectorized instructions. 2.4 Ensemble MCMC Methods The use of multiple paths have been used in rendering to reduce variance and better explore path space. These methods have largely focused on variants of parallel tempering, also known as Replica Exchange Monte Carlo [Swendsen and Wang 1986]. This uses multiple Markov Chains to explore different spaces, and uses a detailed balance preserving transition to swap chains between spaces. This was introduced to graphics by Kitaoka et al. [2009], and improved by Šik and Křivánek [2016] and Otsu et al. [2013]. These approaches have also been applied to progressive photon mapping [Hachisuka and Jensen 2011], the combination VCM/UPS with MCMC [Šik et al. 2016], and with stratified MCMC on the image plane [Gruson et al. 2020]. Hachisuka et al. [2014]alsouseda pool of chains of different lengths to sample path lengths proportional to their contribution. Another related approach is to use Population Monte Carlo (PMC) [Cappé et al. 2004; Fan et al. 2007; Lai et al. 2007]. This iteratively and adaptively samples and resamples a population of paths proportional to their contribution and guides future samples, typically by adapting the parameters of distributions or kernels used to generate samples. While this is related to our approach, it is not trivial to combine PMC with MCMC methods without biasing the result, and it is also not clear how this approach can be applied when computing high dimensional integrals. One approach outside of the graphics literature which is closely related to our work is Affine Invariant Sampling (AIS) [Goodman and Weare 2010]. This work considered a pool or ensemble of walkers ∈RN, and used the states of all other walkers to guide ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. 5:4 • T. Bashford-Rogers et al. perturbations for each walker. The authors proposed three perturbations: stretch moves which shift a walker’s position toward or away from a randomly sampled walker in the ensemble, a walk move which samples a subset of walkers and builds a Gaussian transition kernel, and a replacement move which aims to reconstruct the whole space and sample from the reconstruction. These techniques were shown to efficiently guide the sampling, especially in the case of complicated distributions. This was further extended by Foreman-Mackey et al. [2013] which proposed a parallel approach to using the ensemble of walkers. We also use a similar approach of partitioning the walkers into two pools, and using one pool to guide sampling in the complementary pool. 3 ENSEMBLE METROPOLIS LIGHT TRANSPORT Our approach is based on an ensemble of paths which capture global information of the distribution of lighting in a scene to guide sampling for each path in the ensemble. We first define an ensemble of chains containing Opaths: X={x1 ,x2 ,..,xO}.(5) This ensemble can be considered to be in PO. Similar to the argument in Goodman and Weare [2010], if we consider a product density using the ensemble F(X)=f(x1)f(x2).. f(xO), then any MCMC algorithm which preserves this density is valid. Such a strategy is to update each path in the ensemble conditioned on the other paths in the ensemble, i.e., following partial resampling [Liu 2008], if when updating path xithe remaining paths in the ensemble {x1 ,..,xi−1 ,xi+1 ,..,xO}remain fixed, then the update of the ith path of the ensemble preserves the joint distribution F. This also allows the other paths in the ensemble to guide the sampling of each path of the ensemble. This implies updating each path in series, as each update relies on fixing the states of all other paths in the ensemble. However, paths can be updated in parallel for all xi∈Xby defining a complementary ensemble, Y={y1 ,y2 ,..,yO}, to guide sampling for each path in X[Foreman-Mackey et al. 2013]. Therefore, each path in Xcan be processed in parallel using Yas guidance for sampling, i.e., the transition kernel takes the form T(xi→xi|Y). This transition kernel can be written as the product of multiple sampling events; in the case of light transport this corresponds to progressively sampling a subpath: T(xi→xi|Y)= k  j=1 K(xi j→xi j|Y),(6) where K(xi j→xi j|Y)is a transition kernel for the jth sampling event of kevents. Specifically, this is the transition kernel associated with perturbing the direction of a path vertex conditioned on the set of paths from the complementary ensemble. This transition kernel can be applied to one or more path vertices, producing a perturbation to a light path. The acceptance probability for updating each path in the ensemble is therefore computed as a(xi→xi)=min  1, f∗(xi)T(xi→xi|Y) f∗(xi)T(xi→xi|Y) .(7) ALGORITHM 1: The EMLT algorithm. Two ensembles of paths X and Yare input, and during rendering, paths from the ensemble Xare processed in parallel it times, and the guided transitions kernels based on Yare used to propose new paths. After all paths in Xare processed, Xand Yare swapped. Input:Xand Y 1while rendering do 2ParFor xi∈X 3for it iterations do See Section 3.4 4xi∼T(xi→xi|Y)See Section 3.3 5a←a(xi→xi|Y)Equation (7) 6Accumulate to Image 7if ξ<athen 8xi←xi 9end 10 end 11 end 12 Swap Xand Y 13 end Once all the paths of the ensemble Xhave been updated, this is referred to as an iteration, the ensembles are swapped X↔Yand paths of Yare updated based on using Xas path guidance:T(yi→ yi|X). However, without loss of generality we refer to Yas the complementary ensemble in the remainder of the text. Algorithm 1 summarizes the EMLT algorithm. Firstly, the two ensembles Xand Yare initialized, then during rendering each path is processed in parallelit times (lines 2 and 3) using the proposed guided transition kernels (lines 4–9). When all paths in an ensemble are processed, then ensembles are swapped (line 12), and the process repeats. The use of ensembles for guiding sampling of paths could be applied to either path space or PSS. One possibility is to apply the use of ensembles to PSS through a strategy which directly perturbs a point in PSS based on other points in the ensemble, similar to AIS. However, due to the difference in the number of random numbers required to sample paths, it is not clear how walkers of different dimensionalities could be used to create any of the transition kernels proposed by Goodman and Weare [2010]. Secondly, interpolating between points in high dimensions, which is the result of applying AIS to PSS, is unlikely to lead to usable paths, especially if there are small regions in PSS containing valid light transport paths. However, the alternative of applying this to path space is also not trivial as samples are no longer in RN, and the strategies outlined in Goodman and Weare [2010] are not immediately applicable. Our proposed transition kernels are designed to be suitable for path space, but are also constructed to inherit the advantages of using an ensemble to guide sampling. This then allows scope for a wide range of new guided transition kernels which are conditioned on the complementary ensemble. While the entire complementary ensemble could be used to create transition kernels, this would be prohibitively expensive when the complementary ensemble is large. An alternative, and significantly faster, approach that we propose in this article is to use a carefully chosen subset of the paths in the complementary ensemble. These paths should be similar, both in interaction types and spatial proximity, such that they can still produce valid guided transition ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. Ensemble Metropolis Light Transport • 5:5 kernels. Section 3.1 describes how to efficiently find and weight the subset of paths from the complementary ensemble, then Section 3.2 describes how guided transition kernels can be constructed from this subset of paths. Finally, Section 3.3 describes how these guided transition kernels can be combined into path perturbations. 3.1 Complementary Ensemble Before describing the transition kernels, we first describe two aspects of using the complementary ensemble for sampling. The first, explained in Section 3.1.1, is how to select paths from the complementary ensemble for sampling. This is important as many of the guided perturbations require paths to be sampled that maintain the same number of path vertices with the same interaction types as the original path. The second aspect deals with the similarity of paths sampled from the complementary ensemble to the original path. This is required as although paths with the same length and interaction types may be sampled from the ensemble, paths which are similar to the original path are likely to lead to better proposal distributions than those further away. Effective use of this similarity between paths is what allows our approach to avoid a spatial cache. Section 3.1.2 describes an approach for measuring similarity between paths. 3.1.1 Finding Paths. As discussed previously, the proposed guided perturbations rely on a subset of Mpaths from Y:ϒ= {υ1 ..υM}.TheseMpaths are located in Ybased on similar properties to a base path xi, such as identical length or the same Heckbert notation interaction types [Heckbert 1990]. This is motivated by two observations: (i) perturbations guided by similar paths, rather than all paths, are likely to explore similar regions of path space leading to higher acceptance probabilities and (ii) perturbations are more likely to succeed since they rely on preserving interaction types. Therefore, the set of paths in ϒis deterministically selected from the paths in Ywith similar properties to xi. This is facilitated by a tree data structure over path lengths and interaction types which can be queried in O(1)time to find a subset of Ywhich matches the desired properties. This is built at the start of the rendering process, or at the end of each iteration, and please see the supplementary material for more details about the construction and traversal of this data structure. If the selection of paths forming ϒwas probabilistic and dependent on xi, then the probability of sampling the set ϒgiven xiwould have to be computed taking into account all paths with similar properties in Ywhich would be prohibitively slow. By performing a deterministic selection, in our case based on a counter which is stored with the ensemble and updated each step, this has the effect of having a minimal impact on performance with the additional benefit that the computation of the acceptance probabilities is significantly simplified as the probability of sampling ϒis not required. 3.1.2 Measuring Similarity. The set of paths returned from querying the ensemble, ϒ, may have similar properties to the current light path xi. However, while some of the vertices in the paths returned may have similar positions in world space to xi,others Fig. 2. Measuring similarity between a base path (red circles) and a path from υn∈ϒ(empty circles). This is a product of the similarity between pairs of path vertices with the same index j:S(xi j,υn j). may not. When developing guided transition kernels, it is useful to have a measure of how similar light paths, or vertices within light paths, are to each other. For instance, some transition kernels can benefit from calculating weights for each vertex from ϒas this is likely to provide a good estimate of nearby lighting. The use of entire light paths in MCMC methods widens the range of methods to measure similarity. While Chaitanya et al. [2018] proposed an effective heuristic of total path length, i.e., the sum of distances between path vertices, we typically do not need to consider the whole path. We develop a heuristic based on the world space position of a set of vertices from ϒ, and vertices in the current light path xi. Other attributes, such as normals, albedo, or surface roughness could be considered, but we found that using the world space position was effective for computing similarity. Specifically, given the j’th vertex from xi,xi j, and the previous vertex, xi j−1, the similarity value can be computed for all paths in ϒby computing the distance to υn jand υn j−1,n∈[1..M]. We define the difference between the world position of two vertices as d(xi j,υn j)=max (|xi j−υn j|2 ,ϵ)−1where ϵis a small positive constant (we use ϵ=0.0001). From this, we define a normalized similarity value as Sxi j,υn j=2 1+e−dxi j,υn j−1.(8) This scaled sigmoid leads to a larger value when vertices are similar, and smaller the further apart they become. See Figure 2for an illustration of similarity computation. The similarity of multiple vertices starting at the j’th position in the path to the k’th position can be computed as Sxi ,υn ,j,k= k  l=j Sxi l,υn l.(9) 3.2 Guided Transition Kernels We first describe guided transition kernels for a single vertex, and then describe how full perturbation strategies can be built from these individual strategies in the following section. All of these methods require information gathered from the set returned from querying the tree structure ϒ. ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. 5:6 • T. Bashford-Rogers et al. Fig. 3. Given some 2D domain, vertices in ϒcan be projected onto that domain (empty circles). One of these is selected (the orange circle) and used by Linear Transition Kernels to form a ray from that point to a projection of a vertex from the current path (the green circle). A distance along this ray is sampled λwhich generates a new point in this domain (the red circle). Linear Transition Kernel. The simplest form of guided transition kernels are the linear transition kernels, which is suitable for guided sampling on the lens. These operate in R2in our implementation. Given coordinates Cxy(xi)∈R2of the current path vertex and the coordinates of a path vertex from ϒ:Cxy(υn i)∈R2, this generates a proposal along a ray in R2:Cxy(x i)=Cxy(υn i)+λ· (Cxy(xi)−Cxy(υn i)). The distance along the ray, λis sampled from a distribution centered on Cxy(xi). This is illustrated in Figure 3. Goodman and Weare [2010] proposed the stretch move which samples λfrom a distribution λ∼д(c)=1 √c,where c∈[1 1+α,1+α], where α∈R+is a scaling term. This density is symmetricд(c)=cд(1 c)(see Goodman and Weare [2010]), and this leadstotheratioofK(xi j→xi j|Y) K(xi j→xi j|Y) =λ, simplifying the acceptance probability. However, other 1D distributions can be sampled to generate λ. For example, a uniform λ∼[1 −β,1+β],β∈R+, or truncated Gaussian can be used, and so long as these are symmetric they simplify in the computation of the acceptance probability. Linear Hemispherical Transition Kernels. While linear transition kernels are defined in R2, many transition kernels are required to be defined over the (hemi)sphere S2. Therefore, we extend the linear transition kernels to the (hemi)sphere. This starts by sampling λfrom one of the linear distributions, then mapping this to a perturbation of the original direction on the sphere. Given two directions in the sphere ω1and ω2, these directions may correspond to a direction on the original path, and the other on a path from ϒ, a new direction ωnew can be sampled along the great arc connecting these two directions: ωnew =ω1sin((ω1·ω2)λ) sin((ω1·ω2)) + ω2sin((ω1·ω2)(1−λ)) sin((ω1·ω2)) , i.e., a slerp between ω1and ω2with parameter 1 −λ(see Figure 4). This leads to K(xi j→xi j|Y) K(xi j→xi j|Y) =sin cos−1(ω1·ω2) sin cos−1(ωnew ·ω2) =1−(ω1·ω2)2 1−(ωnew ·ω2)2. (10) Guided Anisotropic Transition Kernels. The linear and hemispherical transition kernels rely on a single path from the complementary ensemble, and domains in R2or S2.However,more Fig. 4. The Linear Hemispherical Transition Kernel uses the outgoing direction of a path from ϒ(the orange circle), and the current path (the green circle) to propose a new direction (the red circle) along the great arc denoted by the dashed line. Other directions from ϒwhich are not considered are shown as empty circles. information can be gained from utilizing all Mpaths in ϒ.Forexample, a distribution in world space can be fit to the vertices at a certain point along the path, recentered at the current path vertex, and this distribution can be used for sampling. This allows lighting information from multiple paths to inform sampling of the current path, similar to Reibold et al. [2018]. There are multiple methods to achieve this; we describe one such approach. We start with the j’thvertexinapathxjand another vertex in the scene x j−1, and then retrieve the set of path vertices from ϒwhich match the index: υn j,υn j−1∈ϒ. For each subpath, we assign a weight: w(n)=S(xi ,υn ,j,j+1) M k=1S(xj,υk,j,j+1).(11) Then each of these points is projected onto the plane defined by xjand the normal at xj:N(xj). Next, an anisotropic Gaussian N(μ,Σ;ϒ,xj)is fitted to these points via weighted maximum likelihood estimation where the weight of each point is that assigned to each path: μ=1 MM n=1w(n)υn jand Σ=1 M−1M n=1w(n)(υn j−μ)2. This is recentered such that μ=xj, leading to N(xj,Σ;ϒ,xj). This recentering is required such that the sampled point is close to the original, and to ensure that the evaluation of the reverse transition kernel returns a value similar to the proposed transition kernel in the computation of the acceptance probability. The steps of this algorithm are shown in Figure 5. Occasionally, all weights can be zero, or a degenerate covariance matrix can be computed. We detect these cases, and revert to a von Mises–Fisher distribution aligned in the direction x j−1→xjwith a high concentration parameter for sampling. Another approach could be to convolve with an isotropic Gaussian similar to Li et al. [2015]; however, the value to use for the variance of the isotropic Gaussian is unclear in our case. Once the anisotropic Gaussian is defined in world space, it is sampled producing a point z j. This point may not be aligned to the scene geometry, so a ray is traced from xj−1in the direction xj−1→z j, producing a new point on the scene manifold x j.The ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. Ensemble Metropolis Light Transport • 5:7 Fig. 5. The procedure to build and sample anisotropic Gaussian transition kernels. Fig. 6. Generating and sampling using the guided directional transition kernels. density w.r.t. area of computing this point is K(xi j→xi j|Y)= N(z j|xj,Σ;ϒ,xj)G(x j↔xj−1) G(z j↔xj−1), where the final term stems from the ratio of geometry terms resulting from the Jacobian from sampling a point on the plane over which N(z j|xj,Σ;ϒ,xj)is defined, to the scene manifold. To compute the acceptance probability, this process has to be computed in reverse; the distribution N(x j,Σ;ϒ,x j)is first computed, then the vertex xjis projected onto the plane defined by x j and N(xj)leading to a point zj. This leads to a resulting ratio: K(xi j→xi j|Y) K(xi j→xi j|Y) =N(zj|x j,Σ;ϒ,x j)G(z j↔xj−1)G(xj↔xj−1) N(z j|xj,Σ;ϒ,xj)G(zj↔xj−1)G(x j↔xj−1). (12) A simpler version of this approach can be used on the image plane. In this case, an anisotropic Gaussian can be fit to the image plane coordinates of the paths in ϒ, each weighted by the similarity measure. Again, this can be centered at the image plane coordinates of the current path, and a new point on the image plane for the proposed path can be sampled from this distribution. Guided Directional Transition Kernels. Guided anisotropic transition kernels form an anisotropic distribution in world space. However, sometimes it is useful to sample perturbations over solid angle. The linear hemispherical transition kernel performs this, but restricted along a great arc. Another approach is to fit a distribution on S2. Various approaches for this exist, for example tabulated, spherical Gaussian or a von Mises–Fisher distribution. Any distribution on the sphere whose parameters can be estimated from a set of directions can be used. Given a set of normalized directions from some base vertex xjto each member of ϒ,ωn=xj→υn j, and weights computed in the same manner as Equation (11), the parameters of a distribution can be estimated. Similar to the guided anisotropic perturbations, this distribution is recentered around the original direction from the vertex. This can then be sampled generating directions which are guided by nearby paths. Figure 6 shows this process. 3.3 Guided Perturbation Strategies The previous section defined a range of guided transition kernels which are designed to update individual path vertices guided by global information from the ensemble. When perturbing a path, these guided transition kernels can be combined into a wide range of guided perturbation strategies designed to explore different lighting effects. Note that these can be combined with the original mutation and perturbation strategies; this simply adds to the strategies available. We always include the bidirectional mutation strategy from Veach and Guibas [1997] as this ensures ergodicity, thereby guaranteeing that the whole space will be explored. ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. 5:8 • T. Bashford-Rogers et al. The following lists the strategies we have implemented, but many more can be built using combinations of the kernels defined in Section 3.2. Linear Lens Perturbation. Fig. 7. Linear Lens Perturbation. The Linear Lens Perturbation uses the linear transition kernel on the image plane, then similar to Veach and Guibas [1997]tracesa subpath over any specular vertices, then connects to the original path. As ϒcontains paths with similar interaction types and lengths to the current path, this strategy aims to explore the image plane around similar interaction types and also typically moves paths toward higher contribution regions for that path type. This strategy samples a path for the perturbation from ϒ, where the weight for the k’th path is given by S(xi,υk,j,k) N p=1S(xi,υp,j,k),wherejis the index of a vertex on the camera, and kis the index of the first non-specular vertex in the path. Figure 7illustrates this strategy. Linear Caustic Perturbation. Fig. 8. Linear Caustic Perturbation. Linear Caustic Perturbation uses the linear hemispherical transition kernel to perturb the sampled direction on the hemisphere to take into account nearby caustic paths. In this case, ϒ will only contain caustic paths of the same number of path vertices, so they are likely to be exploring a similar region of the scene. This first finds the best path from ϒwhich closest matches the starting point and first specular vertex of the original caustic subpath, and sets the directions ω1=xi c→xi c−1,and ω2=xi c→υn c−1. This then perturbs the direction on the hemisphere, and traces the specular subpath to the first diffuse vertex, and connects to the camera. This is visualized in Figure 8. Linear Multi-Chain Perturbation. Fig. 9. Linear Multi-chain Perturbation. Similar to the multichain strategy described in Veach and Guibas [1997], Linear Multi-chain Perturbation uses the linear transition kernel on the image plane in the same way as the Linear Lens Perturbation. This then traces a specular chain until a non-specular vertex is generated. A deterministic connection to the next specular subpath is then made and this process repeats until the path can be reconnected to the light subpath of the original path (see Figure 9). Anisotropic Path Perturbation. Fig. 10. The Anisotropic Path Perturbation uses the guided anisotropic perturbations for the first and second non-specular interactions from the camera. This can either start from the camera (left image), or toward the camera (right image). This perturbation strategy comprises using the guided anisotropic transition kernels to perturb the current path. This can be applied to any number of path vertices, either from the light source or the eye. As the process of fitting an anisotropic Gaussian is relatively expensive, we restrict this perturbation to the first two vertices from the camera, and randomly select whether to sample from or toward the camera. If sampling from the camera is selected, an anisotropic Gaussian is created on the image plane and sampled, and for all other non-specular interactions the guided anisotropic transition kernels in world space are used, then reconnected to the original path. Likewise, if sampling toward the camera is selected, the guided anisotropic transition kernel is used to generate path vertices which are deterministically connected to the camera. This perturbation strategy helps to explore the local region around the path, and sampling more than one vertex from the camera helps to minimize the impact of the weak singularity in the geometry term near edges visible from the camera (see Figure 10). Another option is to sample with respect to solid angle similar to Otsu et al. [2018] or using the guided hemispherical perturbation which would cancel geometry terms in the calculation of the acceptance probability. Environment Perturbation. Fig. 11. Environment Perturbation. This perturbation is designed to explore environment lighting based on nearby paths. This perturbation uses the linear hemispherical transition kernel or the guided directional transition kernel applied toward the environment map to perturb the direction to the environment map, assuming that the first path vertex from the light is non-specular. This is illustrated in Figure 11. Variants of this strategy can also be applied to other area light sources, or for light source selection. ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021. Ensemble Metropolis Light Transport • 5:9 3.4 Implementation Details This method is initialized similar to the resampling approach described in MLT [Veach and Guibas 1997].Alargesetofpaths is computed with bidirectional path tracing, the contributions of these paths stored and used to scale the result, then these paths are resampled into two subsets: one to generate the ensemble X, and the other used to generate the complementary ensemble Y. There is no limit on the size of Xand Y; however, for convenience we choose them to be the same size of |X|=|Y|=16,384 (please see the supplementary material for further analysis). There is also much flexibility about when to update and swap ensembles (see Section 3for details). Although updating the data structures for finding paths sampling is relatively inexpensive, this does come with some computational overhead of clearing out previous values and reinserting new values. Therefore, each path in the ensemble is perturbed or mutated it times (see Algorithm 1) before the ensembles are swapped, which amortizes the overhead of updating the data structure compared to swapping after every mutation or perturbation. We choose it to take a value of W×H |X|+|Y|,whereW and Hare the width and height of the image plane, respectively, which balances the computation cost of rebuilding the data structures and runtime performance. There is significant freedom to choose the value of αused for the linear transition kernels described in Section 3.2.However,the method does not work well if this is set to a constant, as if the paths in ϒare clustered in a small region of space, e.g., in a small area on the image plane, then αshould be large to facilitate exploration of a small space. Conversely, when paths are spread over a large area αshould be small such that the proposed path is able to explore a similar region of the space. We solve this issue by adapting α based on ϒ. For a linear perturbation on the lens, we first compute a ratio of the bounding box of the image plane coordinates of each path in ϒto the image plane resolution: bres.αis then computed by linearly interpolating between two bounds αland αsbased on aweightwα=(1+e−bres−c σ2)−1,wherec∈[0..1]. This uses a generalized sigmoid as a weighting function as it gives control over where and how fast the weights transition from 0 to 1. We use the parameters αl=0.5, αs=0.05, c=0.1, and σ2=0.02, although the algorithm is quite robust to these values. The supplementary material provides further details on the impact of α. For the linear lens, caustic, and multi-chain perturbations, if the number of paths in ϒis less than two, then the original perturbation strategies are used. This is to handle two situations: one is if no nearby paths are found, then the path can still be perturbed, and secondly, if only one path is found then there is too little information about nearby paths to create a useful sampling distribution. We also set the probabilities of sampling each proposed mutation type to be equal. 4RESULTS EMLT, MLT [Veach and Guibas 1997],andGeometryAwareMLT (GAMLT) [Otsu et al. 2018] were implemented into the same rendering framework for comparison. We implemented the Bidirectional, Lens, Caustic, and Multichain perturbations in MLT, and compare to GAMLT as it is the closest method to ours in terms of using adaptive sized perturbations in path space. We tested the methods in a variety of scenes, from those which exhibit challenging light transport where MCMC methods are expected to perform well, to simpler scenes which represent more common use cases for rendering. All results were computed on a laptop with an i78750H and 16GB RAM. Computation was spread over 12 threads using a thread pool to process paths in parallel and the ensemble was the same size per scene (see Section 3.4 for more information). We set a constant probability for the bidirectional mutation of 1 3. All results were rendered at an average of 64 mutations per pixel to allow equal comparison between methods. 4.1 Indirect Lighting Our method is primarily focused on efficiently computing global illumination. Therefore, we first investigate the performance of EMLT in scenes with indirect lighting only, as direct lighting can be efficiently handled by other techniques in these scenes when compared to MLT. We show results for six scenes which exhibit different types of lighting effects. The DOOR AJAR scene in Figure 1is a challenging scenario where light propagates through the ajar door. Similarly, the BEDROOM scene in Figure 15 has thick, diffuse curtains with a light source on the other side leading to a very challenging lighting scenario. The CLASSROOM scene in Figure 15 is lit by an environment map with light entering through the windows. The KITCHEN scene in Figure 15 shows strong indirect lighting on the back wall and glossy reflections. The CORNELL BOX scene in Figure 16 is representative of many real-world scenes with simple lighting configurations. Finally, the STAIRCASE scene exhibits simple indirect lighting above the stairs, and more complicated indirect lighting under the stairs. Insets in the images show details, and the values printed on the top left of full resolution images correspond to Mean Squared Error (MSE) for the whole image. In Figure 17, we show loglog convergence plots for MSE versus average mutations per pixel for the scenes used in this article to show how error decreases. This shows that there is an improvement in convergence using EMLT (blue line) compared to MLT (green line) and GAMLT (red line; see below for further discussion). This can be seen in the rendered images as a reduction in noise compared to MLT. Diffuse and low glossy surfaces, such as the walls in the DOOR AJAR, CLASSROOM, and CORNELL BOX scenes, or behind the cooker in the KITCHEN scene, exhibit significantly reduced variance with EMLT. This is due to the transition kernels adapting to both the illumination and scene as encoded in the paths in the ensemble. However, EMLT also captures higher frequency lighting effects, and can adapt to higher glossy materials, as can be seen in the BEDROOM scene above the curtains, the metal on the chairs in the CLASSROOM scene, and the strong indirect lighting on the back wall in the KITCHEN scene. EMLT leads to improvements in MSE for all tested scenes comparedtoMLT:2.23×for DOOR AJAR, 1.38×for the BEDROOM scene, 2.11×for the CLASSROOM scene, 1.44×for KITCHEN, 1.71×for the STAIRCASE scene, and 2.93×for the CORNELL BOX scene. ACM Transactions on Graphics, Vol. 41, No. 1, Article 5. Publication date: December 2021.