Hao Qin* Xin Sun † Qiming Hou* f Baining Guo† Kun Zhou* † Microsoft Research Asia * State Key Lab of CAD&CG, Zhejiang University Figure 1: A comparison between stochastic progressive photon mapping (SPPM), unified path sampling/vertex connection and merging (UPS/VCM), and our unbiased photon gathering with bidirectional path tracing (UPG+BDPT) after I hour of rendering. SPPM utilizes biased photon mapping to produce a low-variance result, at the cost of over-blurring sharp features.  UPSNCM gains extra benefits from BDPT but the vertex merging part is still biased.  Our method is both unbiased and robust, producing a result that most resembles the reference. Note that the left inset is set to an exposure of 1/64 to make the HDR shadow details visible.
Abstract
Photon mapping (PM) has been widely regarded as an efficient so- lution for light transport simulation, including challenging caustics paths and many-bounce indirect lighting.  The efficiency of PM comes from reusing traced photons. However, the handling of pho- ton gathering in existing PM algorithms is universally biased – the expected value of their results does not necessarily agree with the true solution of the rendering equation.  We present a novel pho- ton gathering method to efficiently achieve unbiased rendering with photon mapping.  Instead of aggregating the gathered photons into an estimated density as in classical photon mapping, we process each photon individually and connect the corresponding light sub- path with the eye sub-path that generates the gather point, creating an unbiased path sample. The Monte Carlo estimate for such a path sample is calculated by evaluating all relevant terms in a strict and unbiased way, leading to a self-contained unbiased sampling tech- nique.  We further develop a set of multiple importance sampling (MIS) weights that allow our method to be optimally combined with bidirectional path tracing (BDPT), resulting in an unbiased render- ing algorithm that can efficiently handle a wide variety of light paths and that compares favorably with previous algorithms. Experiments demonstrate the efficacy and robustness of our method. CR Categories:   I.3.3 [Computer Graphics]:  Three-Dimensional Graphics and Realism——Raytracing Keywords:  global illumination, photon mapping, BDPT
1
Introduction
Extensive study has been devoted to the robust simulation of light transport in virtual scenes.  Photon mapping (PM) [Jensen 2001] has been widely regarded as one of the most efficient solutions when handling challenging caustics paths, or when significant il- lumination comes from many-bounce indirect lighting.  However, its estimation of light transport is usually biased by photon gather- ing, which means that the expected value of rendered results does not necessarily agree with the true solution of the rendering equa- tion [Kajiya 1986].  While progressive variants [Hachisuka et al. 2008; Hachisuka and Jensen 2009; Knaus and Zwicker 2011; Ka- planyan and Dachsbacher 2013a] eventually converge to the correct results, they remain biased at any finite number of samples. Both the efficiency and the bias in photon mapping come from the photon gathering step, during which eye sub-paths are traced from the sensor, and a range search is performed at their endpoints to query a photon map traced from the light sources in a previous pass. The efficiency comes from the reuse of sub-paths corresponding to the gathered photons.  In classical PM methods, the gathered photons are aggregated into an estimated density that approximates local radiance or irradiance. Such an approximation produces bias, typically in the form of blurred results. As shown in Fig. 1 (a) and (b), the shadow in the zoom-in appears blurred. In this paper, we propose a novel photon gathering method as a re- placement of the density estimation to efficiently achieve unbiased rendering with photon mapping. Specifically, instead of aggregat- ing the gathered photons, we process each photon individually and connect the corresponding light sub-path with the eye sub-path that generates the gather point. This connection can be interpreted as a Russian roulette event decided by the gather point. By removing the gather point itself from the completed path, we obtain a path sam- ple that could contribute to the final radiance estimation. While this path reuse scheme has been formulated in previous work [Bekaert et al. 2003], the roulette probability is approximated in the denom- inator of the final estimate, resulting in a significant source of bias in the rendering result. Therefore, we introduce an unbiased algo- rithm to estimate the reciprocal of the roulette probability using a series of Bernoulli trials, making it possible to calculate the Monte Carlo estimate of the path sample by precisely evaluating all its terms.  The end result is a self-contained unbiased sampling tech- nique, which can be incorporated into a wide variety of Monte Carlo light transport algorithms. We additionally develop a set of multi- ple importance sampling (MIS) weights that allow our method to be optimally combined with traditional techniques. The proposed unbiased photon gathering (UPG) can be combined with traditional unbiased techniques for increased robustness.  We first combine it with bidirectional path tracing (BDPT) using MIS, in a manner analogous to unified path sampling (UPS) [Hachisuka et al. 2012] and vertex connection and merging (VCM) [Georgiev et al. 2012], resulting in an unbiased rendering algorithm that can efficiently handle a wide variety of light paths. The unbiased nature of our method also makes it easy to serve as an additional sam- pling technique in primary sampling space Metropolis light trans- port (PSSMLT) [Kelemen et al. 2002]. For glossy and many-bounce light paths, experimental results demonstrate that our method compares favorably to state-of-the- art alternatives, including multiplexed Metropolis light transport (MMLT) [Hachisuka et al. 2014] and UPS/VCM. While our UPG has a higher per-photon overhead than UPS/VCM, the unbiased na- ture of our method allows it to maintain a stable sampling speed, i.e., the number of samples generated per unit time, throughout the rendering, which is higher than the effective BDPT sampling speed that UPS/VCM eventually converges to.  While UPG loses the ability to render strict specular-diffuse-specular ( SDS ) paths as a cost of bias-free rendering, such effects can be approximated using highly glossy surfaces, producing plausible results similar to those generated by UPS/VCM (see Fig. 13).
2
Related Work
Photon Mapping (PM).  The central idea of PM [Jensen 2001] is to approximate local radiance using photon density estimated in photon gathering. 1  The locality of photon gathering makes it ef- ficient to sample challenging caustics paths and many-bounce indi- rect lighting. In addition, PM can reuse light sub-paths, which in- creases the effective sampling rate with minimal cost. The density estimation suppresses sampling noise, but meanwhile, introduces biases. Consequently, PM is predominantly regarded as biased, al- though state-of-the-art variants like Progressive Photon Mapping (PPM) [Hachisuka et al. 2008] and Stochastic Progressive Pho- ton Mapping (SPPM) are consistent [Hachisuka and Jensen 2009; Knaus and Zwicker 2011; Kaplanyan and Dachsbacher 2013a] and converge to the correct result when the sampling rate ap- proaches infinity.  MIS [Veach 1998] has been incorporated into PM [Tokuyoshi 2009; Vorba 2011] approaches to combine multi- ple techniques, although in a biased manner. As an attempt to reduce bias, several researchers proposed to re- move the gather point and connect its prefix vertex to the photon with a shadow ray [Bekaert et al. 2003]. However, the bias reduc- tion is partial at best since they approximated the essential proba- bility term (corresponding to a normalization factor in their paper). Although the way we generate light paths is technically equivalent to Bekaert et al. [2003], we choose to interpret photon gathering as a stochastic process like UPS/VCM, and provide an unbiased esti- mator for the created path sample. The recently proposed particle-guided BDPT [Vorba et al. 2014] shoots photons and importons to compute a full-scene irradiance distribution, which is then used to guide individual ray-surface in- teractions in a separate unbiased pass. Their method operates on a different principle than PM, as the photons are only used as guides and do not make a direct contribution to the final image.  Our method is orthogonal to their approach – their guides can be directly applied to improve sub-path sampling in our method. Our prelimi- nary implementation of a combined method has demonstrated con- siderable benefits (see Fig. 12). Bidrectional Path Tracing (BDPT).  The rendering equation [Ka- jiya 1986] is a physically-based formulation of general light trans- port, which can be solved by sampling a path integral using Monte Carlo methods.  BDPT [Lafortune and Willems 1993; Veach and Guibas 1994] samples paths by tracing light sub-paths and eye sub- paths separately and connecting the endpoints with shadow rays, which is especially efficient for indoor scenes with finite-sized light sources and sensors. With MIS and the path space formulation of BDPT [Veach 1998], a wide variety of path sampling techniques, including ours, can be combined for increased robustness. The final result is unbiased as long as all involved techniques are unbiased. Unified Sampling.  Both PM and BDPT sample paths by connect- ing light sub-paths and eye sub-paths. Compared to BDPT connec- tion, the locality between gather points and photons in photon gath- ering improves sampling efficiency especially for caustics paths. Accordingly, PM and BDPT methods generally specialize in mutu- ally complementary light transport types [Hašan et al. 2009; Vorba 2011]. Therefore, an effective combination of both would be sig- nificantly more robust than either approach. The key is to integrate BDPT connection and photon gathering with a unified formulation of path integrals. VCM [Georgiev et al. 2012] redefines photon gathering as a prob- abilistic bidirectional connection which rejects sub-path pairs with spatially distant endpoints, whereas UPS [Hachisuka et al. 2012] considers BDPT connections as a special case of PM where a pho- ton overlaps the gather point exactly.  The two formulations result in an equivalent path-space combination of SPPM [Hachisuka and Jensen 2009] and BDPT, where MIS could be utilized for increased robustness.  However, both formulations rely on photon density esti- mation, and thereby inherit the estimation biases from photon gath- ering. The final result would remain biased as in photon mapping. In contrast, our unbiased photon gathering formulation excludes the approximated evaluation in density estimation. It directly evaluates the Monte Carlo estimate of the path samples in an unbiased way. Our method can be combined with BDPT in a manner analogous to UPS and VCM, but the final result is unbiased. Markov Chain Monte Carlo (MCMC). MCMC was intro- duced to the graphics community as Metropolis light transport (MLT) [Veach and Guibas 1997], which uses the Metropolis- Hastings algorithm [Metropolis et al. 1953] to improve bidirec- tional path sampling.  Since then, a wide variety of mathemati- cal and practical improvements [Kelemen et al. 2002; Cline et al. 2005; Kitaoka et al. 2009; Hachisuka and Jensen 2011; Jakob and Marschner 2012; Kaplanyan and Dachsbacher 2013b; Kaplanyan et al. 2014; Hachisuka et al. 2014] have been proposed in the MCMC framework. However, analogous to other approaches based on path sampling, MCMC algorithms would still become biased when combined with biased methodology [Hachisuka and Jensen 2011; Kaplanyan and Dachsbacher 2013b]. Many-light Methods.  Many-light rendering [Keller 1997; Kollig and Keller 2004; Walter et al. 2005; Walter et al. 2006; Hašan et al. 2009; Davidovič et al. 2010; Walter et al. 2012] is another effec- tive alternative to Monte Carlo approaches which is known for its high scalability and low noise level. Drawing parallels between vir- tual lights and light sub-paths, many-light methods resemble BDPT more than PM as the virtual lights are connected to eye sub-paths with shadow rays rather than local gathering [Segovia et al. 2006; Walter et al. 2012].  Consequently the end-to-end behavior resem- bles that of BDPT, to which PM remains complementary.  This is especially true when handling specific light transport types such as caustics.  While this paper focuses on photon mapping, VPL schemes with probabilistic shadow ray tests can potentially employ an unbiased reciprocal estimation scheme analogous to ours.
3
Unbiased Photon Gathering
Light transport simulation is generally formulated as path integra- tion [Veach 1998], and classical BDPT [Lafortune and Willems 1993; Veach and Guibas 1994] is its unbiased estimation by con- necting light sub-path / eye sub-path pairs with a shadow ray, as shown in Fig. 2 (a). Photon mapping [Jensen 2001] takes a different approach. Instead of connecting a single light sub-path, the incident radiance at the end point of an eye sub-path is approximated by pho- ton density estimated from nearby photons. As in Fig. 2 (b), it does not perform any visibility test and the bidirectional scattering dis- tribution function (BSDF) term is evaluated for a wrong direction, which creates bias. The equivalent formulations in UPS [Hachisuka et al. 2012] and VCM [Georgiev et al. 2012] unify BDPT connec- tion and photon gathering within the same path integral, however, in a biased way. As illustrated in Fig. 2 (c), they formulate the photon gathering as a hypothetical Russian roulette process without veri- fying the visibility of accepted connections.  In addition, they also approximate the connection probability assuming the surface point within the neighborhood is a perfect disc.  Before UPS/VCM, ad- ditional visibility tests have been proposed to improve the accuracy of density estimation [Bekaert et al. 2003] (see Fig. 2 (d)), where the UPS/VCM connection probability corresponds to a normaliza- tion term. However, they make the same disc approximation, which biases the final estimation.  UPG follows the Russian roulette for- mulation of UPS/VCM, but handles the probability term without any approximation as illustrated in Fig. 2 (e). In this section, we start with the basic formulations of the path in- tegral and BDPT, and then review the formulation of vertex merg- ing [Georgiev et al. 2012] which interprets photon gathering as a technique for path sampling (Sec. 3.1). Based on the formulation of vertex merging, we analyze the causes of the biases from photon gathering (Sec. 3.2).  After that, we introduce our formulation of explicit Russian roulette connection to remove the bias in photon gathering (Sec. 3.3). Finally, we describe an unbiased algorithm to estimate the reciprocal of the sampling probability in photon gath- ering using a series of Bernoulli trials (Sec. 3.4), and a rendering al- gorithm using our unbiased photon gathering under the framework of photon mapping (Sec. 3.5).
Review of Path Integral and Vertex Merging
3.1
Path Integral and BDPT.  We use the Monte Carlo path integral formulation of light transport:
where X denotes a path between a light source and a sensor, n is the number of samples, f(x) is the energy it carries and p(x) is the probability density function used to sample ¤.  Classi- cal BDPT connects a light sub-path y 1 · · · y s to an eye sub-path zt · · · z 1 by tracing the shadow ray y s zt . The resulting full path Xs,t = Y1 · · · Ys Zt · · · z1 has a length of k = s + t — 1. We can define the corresponding Monte Carlo estimator in Eq. (1) as:
where α ⊥ ( R s,t ) and α E ( R s,t ) are the light and eye sub-path terms respectively. For now we will ignore their definitions and focus on f c ( { $,t ), which is the measure term of the sub-path connection
where f s is the bidirectional scattering distribution function (BSDF) and G is the geometry term that converts the area measure at z to the projected solid angle measure [Veach 1998] centered at y s , which also takes the mutual visibility into account. Photon Gathering as Vertex Merging. In photon mapping, the incident radiance at the end point of an eye sub-path is evaluated through photon density estimation, which is computed by gather- ing photons within a local neighborhood. The neighborhood size is typically specified as a radius d. According to the equivalent for- mulations of UPS/VCM, photon gathering can be considered as a process that pairs each eye sub-path z t , · · · z 1 with all light sub- paths yi · · · ys ′ satisfying || ys ′ − z t ′ || < d. Here we follow the vertex merging formulation in VCM as it explains our method more elegantly. Vertex merging requires one of the sub-path endpoints to be re- moved from the final path.  Without loss of generality, we removed z t ′ , and the resulting complete path is y 1 · · · y s ′ z t ′ − 1 · · · · z 1 , which we denote as
where the connection term f c ( ¤ s,t ) in Eq. (2) becomes αξν M ( ¯ x ′ ,t ′ − 1 ), which includes both the measurement contribution fV M ( Xs',t' –1) and the connection probability pV M ( Xs',t' –1)
Bias of Photon Gathering
3.2
Vertex merging in Eq. (4) explicitly formulates photon gathering as a hypothetical Russian roulette process.  The probability p ν M is the corresponding acceptance probability, where z t r is interpreted as the end point of a random tentative ray.  The connection be- tween y1 · · · y s ′ and z t ′ − 1 · · · z 1 is only accepted if z t ′ falls within the neighborhood of y s ′ , which corresponds to the neighborhood search in photon gathering.  While this interpretation seemingly goes against engineering intuition, it renders photon gathering as a bidirectional path sampling technique. This allows us to explic- itly formulate the bias in classical photon mapping, which is the first step towards unbiased photon gathering. Comparing Eq. (4) with the unbiased BDPT formulation in Eq. (2), we can see that if α ν M ( × s ′ ,t ′ − 1 ) in Eq. (4) could have been an unbiased Monte Carlo estimator of f c ( x s,t ) in Eq. (2), the estima- tion C ¥ M ( × s ′ ,t ′ − 1 ) would be unbiased just like C * ( × s,t ). This requires f V M to be equal to f c and p V M to be the real probability for the z t , to accept the Russian roulette connection. However, in vertex merging, C V M is simply evaluated using an equivalent form of the density estimation in classical photon mapping, and neither of the two requirements holds true. Specifically, in vertex merging, f ν M in Eq. (4) is not directly re- lated to the connection between y sr and z t ′ − 1
where fs(Yo'-1 -> Yo'(zt') -> zt'-1) indicates that the BSDF is still evaluated at y s ′ but the outgoing direction is taken as z t ′ → Zt' _ 1 instead of y ,  →  zt' _ 1 In contrast, according to Eq. (3) the strict evaluation should be
Comparing Eq. (5) with Eq. (6), one can see that both the geometry term and the BSDF term at y s , are approximated, because the edge ys' {} zt' -1 is replaced by zt' {} zt' -1, which is an approximation that originated from classical photon mapping. The key idea introduced by UPS/VCM is that the photon gathering neighborhood determines the probability of making such a connec- tion.  In VCM, p ý M in Eq. (4) also takes an equivalent form of photon density estimation using a constant disc filtering kernel PV M (X s ′ ,t ′ − 1 ) = πd 2 p x (Zt ′ − 2 → Zt ′ − 1 → Y s ′ ), (7) where ¯ x ( z t ′ − 2  →  z t ′ − 1  →  y s ′ ) is the probability density to sam- ple vertex y s ′ from two preceding vertices z t ′ − 2 and z t ′ − 1 . The tilde indicates that for the purpose of computing ¯ p x , y s , is always assumed to be visible from z In contrast, the strict evaluation has to properly integrate p x over every possible z t r in the gathering radius, which is the same as the acceptance probability without approximation in VCM and gives
The integration domain S ( y s ′ , d ) is the gathering neighborhood, which consists of all surface points within a sphere of radius d cen- tered at y s ′ . Unlike ¯ ρ x , the integrand p x has the correct visibility term embedded and p x ( zt ′− 2 →  zt ′− 1 →  z ) =  0 if z is not visible from z t ′ − 1 − Comparing Eq. (7) with Eq. (8), p έ integration using one sample at p x ( zt ′ − 2  →  zt ′ − 1  →  y s ′ ) without verifying the visibility, and the area of the neighborhood S ( y , , d ) is also approximated as the simple disk area π d 2 . Biases in existing methods.  According to the term-wise compar- ison between the BDPT estimator in Eq. (2) and PM estimator in Eq. (4), we can see that the most essential change photon mapping introduced is the connection probability p ¥ M . Existing methods typically use a disc approximation for this term, whether it is in- terpreted as a Russian roulette probability [Hachisuka et al. 2012; Georgiev et al. 2012] (Fig. 2 (c)) or as a normalization factor in a customized density estimation [Bekaert et al. 2003] (Fig. 2 (d)). The disc approximation is an essential source of estimation biases, and our method is the first to remove it. Details will be introduced in Sec. 3.3 and Sec. 3.4. Note that the path reuse formulation of photon gathering as well as the connection probability have been described in [Bekaert et al. 2003], and are thus not the contributions of this paper.  Our main contribution is an unbiased algorithm to estimate the reciprocal of the connection probability using a series of Bernoulli trials (Sec. 3.3 and Sec. 3.4), making it possible to calculate the Monte Carlo esti- mate of the path sample for the first time.
Explicit Russian Roulette Connection
3.3
Having analyzed the bias of classical photon mapping, a natural direction is to replace Eq. (5) and Eq. (7) with the corresponding strict evaluations in Eq. (6) and Eq. (8) respectively, which forms an unbiased estimation
Eq. (9) is a strict formulation of the hypothetical Russian roulette event in photon gathering. z t ′ is interpreted as a rejection sample independent of 𝑠 s ′ , t ′ − 1 , which conditionally rejects the connection between the two sub-paths y 1 · · · y s ⁄  and z t ⁄ − 1 · · · z 1 . The connec- tion is only accepted if z t , falls within the neighborhood S ( y s , d ), the probability of which is p c . The explicit Russian roulette connection in Eq. (9) differs from ver- tex merging in Eq. (4) by utilizing the strict evaluations of f c and p c . To truly eliminate bias, all terms in the strict formulas Eq. (6) and Eq. (8) must be evaluated in an unbiased manner.  As Eq. (6) can be evaluated precisely by tracing the shadow ray between y s ′ s and z t ′ − 1 , as proposed by Bekaert et al. [2003], the main challenge comes to Eq. (8), which is the connection probability p c . Note that p x ( z t ′ − 2  → z t ′ − 1  → z ) in Eq. (8) depends on the geom- etry term G ( z t ′ − 1  ↔  z ), which additionally depends on a visibility term, making an analytical integration practically impossible. Thus we introduce a separate Monte Carlo estimate to handle p c . How- ever, since p c is ultimately used as a denominator in Eq. (9), even an unbiased Monte Carlo estimate of Eq. (8) would still bias the final result. Jensen's inequality [Chandler 1987] leads to:
The equality only holds true when the variance is zero, which clearly does not hold for general sub-path pairs. Therefore, we must develop an unbiased method that directly estimates the probability reciprocal r(x) = 1/pc(x).
3.4    Unbiased Probability Reciprocal Estimation
Our probability reciprocal estimation is inspired by Booth's work [2007], which introduces an unbiased solution to the more general problem of estimating the reciprocal of an arbitrary integral. Under Booth's formulation, our method is a special case where the unknown portion of the integrand is a binary function (see details in Appendix C). As illustrated in Fig. 3(a), suppose we literally follow through with the hypothetical Russian roulette test described in Sec. 3.1. For each path x s ′ . t ′ − 1 , we generate the tentative ray z t ′ − 1  → z fol- lowing the distribution p x ( z t ′ − 2  →  z t ′ − 1  →  z ).  Based on the definition in Eq. (8), there is a pc ( X s ′ ,t ′ − 1 ) probability for each ray to be accepted by Russian roulette. Repeatedly generating tentative rays for the same path X s ′ t ′ − 1 gives a series of Bernoulli trials. Assuming that the first accepted ray comes after N( ¤ s ′ ,t ′ − 1 ) tries, one could see that the random variable N( X s ′ ,t ′ − 1 ) follows a geo- metric distribution
and the expectation of N is exactly what
In summary, an unbiased estimation of the probability reciprocal can be computed as the number of tentative rays tested until the first successful Bernoulli trial, i.e., up to and including the first ray that hits the gathering sphere S ( y , , d ). Angular Bound for Tentative Rays. An obvious issue of the aforementioned algorithm is the potentially unbounded number of Bernoulli trials N . However, note that the Bernoulli trials are only triggered with a probability of p c , and N is an estimation of its re- ciprocal. When amortized over all connections, the expected num- ber of Bernoulli trials per connection is simply:
Thus, on average the total number of required Bernoulli trials is equal to the total number of connections. However, since the Rus- sian roulette formulation in Sec. 3.1 assumes all sub-path pairs are tested for connection, the bound is O ( n 2 ) for n sub-paths, which is still prohibitively expensive. To make the Bernoulli trial overhead practical, we must avoid testing potential connections that would be an obvious rejection. As a solution, we exploit the restricted size of the gathering neigh- borhood S( y s , d ). Specifically, we project the neighborhood onto the unit sphere centered around zיי− 1 , and only trace tentative rays within the resulting angular bound as illustrated in Fig. 3(b). Con- sequently, the expected number of Bernoulli trials per connection is reduced by a factor proportional to the squared distance between ys' and z t '_1, providing d ≪ ||ys' = z t '_1 ||. Ignoring the effects of BSDF importance sampling, we have:
where N b ( ¤ s ′ .t ′ − 1 ) is the number of Bernoulli trials after applying the angular bounding. The bound may degenerate to the entire up- per hemisphere for extremely close pairs of z t ′ − 1 and t s ′ . However, such cases only constitute a tiny portion of the Russian roulette con- nections. Based on our experiments, this bounding reduces the total overhead of our algorithm to approximately O ( n ). With N b Bernoulli trials instead of N , the probability reciprocal estimation in Eq. (12) has to be adjusted accordingly:
where Ω b is the set of surface points visible to z t __ 1 and inside the angular bound. p b ( ¯ x s ′ .t ′ − 1 ) is the probability for a tentative ray traced without the angular bound to fall within Ω b . As de- fined in Eq. (13) and Eq. (15), p b ( ¯ s ′ ,t ′ − 1 ) is also the ratio by which the angular bound reduces the number of Bernoulli trials: p b (x s ,t' -1) = E [N b (x s ,t' -1)] /E [N (x s ,t' -1)] Ω b can be chosen arbitrarily as long as it includes all points in S(y s ′ ,d) that are visible to z t ′ − 1 . Note that p ′ ( × s ′ ,τ ′ − 1 ) in Eq. (16) is not a random variable and must be evaluated analytically. This type of integration has been well researched in the context of BSDF importance sampling. We simply bound the random number used to generate the tentative ray direction. Details can be found in the supplementary material, which includes an example mathemat- ical derivation for a combined Lambertian and Phong BSDF. The angular bound also significantly reduces the estimation variance of the probability reciprocal r , as discussed in Appendix A. To prove that the angular bound does not introduce additional bias, we substitute Eq. (15) into Eq. (9) and verify that the expectation of the resulting estimator is still correct
Notice that p c ( { s ′ , t ′ − 1 ) and p b ( { s ′ , t ′ − 1 ) can be taken out of E [ · ] as their values are constant for any given path X s ′ , t ′ − 1 · We note that our method can be regarded as the first application of Booth's method in global illumination. However, a straightfor- ward application of Booth's formulation is prohibitively expensive. We simplify it as a series of Bernoulli trials and optimize it with angular bounding to make it practical, which can be thought of as the auxiliary function in Booth's framework.  This is analogous to the fact that while most unbiased rendering algorithms are applica- tions of the Monte Carlo estimate method, it is the actual sampling technique that is significant in the graphics context.
3.5 Rendering Algorithm
Algorithm 1 summarizes the rendering algorithm using our unbi- ased photon gathering (UPG) under the framework of photon map- ping. Note that the algorithm generates a photon map and performs gathering just like in classical photon mapping. The key difference is in lines 5-16, which generate an unbiased path sample for each photon instead of performing a density estimation. Compared to PPM-based methods, our algorithm shares the same limitation where the photon gathering in line 4 can only use a fixed- radius search rather than the classical k-nearest neighbor (kNN) search [Jensen 2001]. We require an explicit radius to form a well- defined integration domain in Eq. (8). However, unlike PPM, our method is unbiased, which makes progressive rendering a trivial matter of running Algorithm 1 many times and averaging the results (see the error of UPG+BDPT reducing to zero in Fig. 7). Thus, our neighborhood gathering radius d can be tuned arbitrarily and inde- pendently in each rendering session, as d does not affect conver- gence and only impacts computational efficiency. Compared to the photon gathering in UPS/VCM, our method uses an equivalent interpretation to formulate photon mapping in the Algorithm 1 Pseudo code of unbiased photon gathering 1: Generate a fixed number of light sub-paths { y } 2:  Build a photon map using the endpoints { y s ′ } 3:  for each eye sub-path z = zt' . . . z 1 do Gather photons in the radius- d neighborhood S( z t , d ) 4: 5: for the light sub-path y of each gathered photon do Project S ( y s , d ) to an angular bound Ω b around z t ′ − 1 6: N b ← 0 7: 8: do N b ← N b + 1 9: 10: 11:
16: Sample a tentative ray zיוג→ z from Ω x while z & S(y s , d) Create a complete path × = y 1 ··
end for 17:  end for BDPT path space.  However, those approaches focus on combin- ing PM with BDPT in general, and bias is not a major concern. Hence, they utilize approximations for their equivalent terms of p c ( x ). While such approximations are biased, they are simpler and more efficient than our unbiased estimation. This is not surprising as removing bias would unavoidably incur some additional cost. Compared to BDPT, our method has to process a non-trivial number of light and eye sub-paths in batch to be effective. While technically Algorithm 1 remains correct for a single pair of sub-paths, doing so hardly has any practical value due to the additional Bernoulli trials required.  Our scheme retains the inherent behaviors of PM-based approaches, and the most critical factor in efficiency is the reuse of sub-paths [Hachisuka et al. 2012; Georgiev et al. 2012]. Significant benefits are only possible when each eye sub-path can be connected to multiple light sub-paths and vice versa. Finally, two implementation details should be noted.  First, when our method is used in conjunction with BDPT using formulations like UPS/VCM, it is possible to reuse the same set of light and eye sub-paths in both algorithms. As previously discussed, when con- necting the same pair of sub-paths y 1 · · · y s and z t · · · z 1 , BDPT generates a path y 1 · · · y s ž t · · · z 1 which has s + t — 1 edges, whereas our method generates a path y 1 · · · y s Zt – 1 · · · z 1 which has s + t − 2 edges with z t excluded. The complete paths do not overlap in the path space, and they can be treated as independent in practice. Second, in the previous sections we have only discussed about the path X s ′ t ′ − 1 generated by removing z t ′ . In practice one could also consider the alternative that removes y s , instead and gener- ates we found that doing so is not cost-effective.  The more efficient approach is to evaluate one of the two candidates and discard the other, where we make the decision heuristically by removing the vertex generated from the less glossy reflection.  A more detailed explanation is provided in Appendix B. Our UPG can also be taken as an individual unbiased sampling technique and combined with other sampling techniques under the unified framework of MIS. This is analogous to UPS/VCM, where a biased SPPM is combined instead. We follow the VCM framework for combination, since we choose to exclude z t r from the complete path rather than including both endpoints in an extended path space. In Sec. 4 and Sec. 5, UPG is combined with BDPT and PSSMLT re- spectively. The combination of UPG with another unbiased method is, of course, still unbiased.
MIS Combination with BDPT
4
Now we describe how to combine our UPG with BDPT using MIS. With multiple sampling techniques, MIS formulates the Monte Carlo path integral in Eq. (1) as
where m is the number of sampling techniques, n is the sample count taken for technique i. p i ( x ) and w i ( x ) are the corresponding path probability density values and MIS weights. For the combi- nation to be unbiased, the weights w i ( ) must be consistent when the same x is generated from different techniques. In addition, we must meet the normalization and non-negativity requirements:
We use the widely-accepted power heuris.
where β > 0 is a predefined parameter. We take β = 2 as rec- ommended by Veach [1998].  Since Eq. (20) satisfies Eq. (19) by definition, the main implementation challenge lies in computing the relevant probability density terms p j for paths sampled using a different technique i . Given a path x with k edges, there are k + 1 BDPT sampling techniques and k − 1 UPG techniques, corresponding to different choices of the pre-connection sub-path lengths. The number differ- ence comes from the fact that BDPT allows zero-length sub-paths. Regardless of the technique chosen to generate  å, it is straightfor- ward to compute the probability density of BDPT techniques di- rectly from the complete path:
where p 1 and p E are the classical sub-path probability density terms [Veach 1998] included in Eq. (2)
where f L ( Xs,t ) and f L ( Xs,t ) are the measurement contribution terms corresponding to the light and eye sub-paths respectively. On the other hand, the probability density of unbiased photon gath- ering is defined as
in which the connection probability p c cannot be evaluated ana- lytically.  While an unbiased estimation is technically possible, it would not serve the purpose here, since deterministic approxima- tion of MIS weights would not introduce bias as long as Eq. (19) is still satisfied. Hence we revert to the approximation of Eq. (7) used by UPS/VCM when computing p c for MIS weights. However, Eq. (7) may disagree with our real connection probability in Eq. (8). Therefore, a sufficiently large error can lead to subop- timal combination, which is most pronounced when y sr and z t ′ − 1 are spatially close. We thus propose a heuristic clamping to allevi- ate this problem
The intuition behind Eq. (24) is simple. Note that p c is a unit-less probability of a hypothetical Russian roulette event and its correct value should never exceed 1. However, p x is a probability den- sity per unit area and by definition p x ( zt ′ − 2  →  zt ′ − 1  →  y s ′ ) is inversely proportional to ∥ z t ′ − 1 − Y z ′ ∥ 2 . By assuming p x to be constant within a radius of d , Eq. (7) implicitly assumes that d ≪ ‖ z t ′ − 1 − y x ′ ‖ . When the assumption fails, the resulting approximation could produce a value significantly larger than 1. Clamping the value back to 1 restricts this error, as shown in Fig. 4. The clamped result (b) computed using Eq. (24) is noticeably less noisy than the naive MIS result in (a) computed using Eq. (7). Note that such probability clamping is not required in UPS/VCM. Their approaches use photon density estimation to approximate the path integral and produce a smooth result regardless of how small || Zt ′ − 1 − Y s ′ || is. The "probability" greater than 1 implicitly makes the smoother PM result dominate the MIS combination, which is usually desirable in a biased method. Fig. 5 illustrates the effect of combining BDPT and UPG. Note that UPG (b) converges significantly better for many-bounce indirect paths, while BDPT (a) is more efficient for diffuse direct lighting. Like in UPS/VCM, the MIS combination (c) gets the best of both methods and is robust for both caustics and direct lighting.
MCMC Integration
5
MCMC methods sample light transport by sampling a Markov Chain, the state of which typically represents one or more paths. The Metropolis-Hastings algorithm is applied to ensure that the sampling probability density of each state is proportional to a Monte Carlo estimator computed from the corresponding paths. A path sampling technique can either be added as a new type of state mutation or as an extra term in the per-state estimator. We take the later approach and build upon PSSMLT [Kelemen et al. 2002]. The classical PSSMLT state is a conceptually infinite sequence of random numbers, which is used to sample an eye sub-path with s vertices and a light sub-path with t vertices. s × t complete paths are generated by connecting sub-path prefixes as in classical BDPT. The MIS combination of the respective path integral estimators is taken as the MCMC sampling probability density. To add unbiased photon gathering, we first pre-trace a fixed set of light sub-paths as in lines 1-2 of Algorithm 1. When evaluating each PSSMLT state, we take the original eye sub-path, and execute the gathering process in lines 4-16 on all its t – 1 prefixes. The resulting estimators are inserted to the per-state MIS combination, which is then used in the original PSSMLT mutations. The pre-traced light sub-paths are regenerated periodically to minimize correlation artifacts. The simplicity of combining UPG with MCMC methods stems from its unbiased nature. Adding UPG to the set of sampling tech- niques of PSSMLT will not affect the target distribution of the Markov Chain process. The combination would be no more than combining UPG with BDPT. While the combination of PPM with MCMC methods needs to either design another target distribution [Hachisuka and Jensen 2011] or carefully deal with the gather ra- dius shrinkage rate [Kaplanyan and Dachsbacher 2013b]. Fig. 6 illustrates how UPG improves the PSSMLT result. Classical PSSMLT (a) does not converge very well for the many-bounce indi- rect lighting, which is a behavior inherited from BDPT. Combining PSSMLT with UPG (b) results in a noticeable improvement. How- ever, currently the UPG+PSSMLT (c) combination does not per- form as well as UPG+BDPT, as the fixed per-state estimator does not allow flexible adjustment of relative sampling rates. Therefore, unless otherwise stated, we default to UPG+BDPT in comparisons.
6
Experimental Results
We conduct a series of experiments to validate our method as well as to compare with alternative light transport algorithms on a com- puter with four Intel i7 3.40GHz CPU cores and 16GB memory. All methods are implemented in the Mitsuba renderer [Jakob 2010]. We choose to compare with multiplexed Metropolis light transport (MMLT) [Hachisuka et al. 2014] as an unbiased alternative and UPS/VCM as a biased alternative. Bias Experiment.  We experimentally validate the unbiased nature of our method by averaging images generated from independent runs.  As illustrated in the top-left row of Fig. 7, the estimation error of UPG+BDPT approaches zero as the number of images in- creases. In contrast, the consistent-yet-biased UPS/VCM does not converge to the ground truth in such a setup. The RMSE (root mean squared error) curve of UPG+BDPT shows up as a line of slope −0.5 in a log-log plot (Fig. 7 right), which follows the classical square root convergence of unbiased Monte Carlo estimate. In con- trast, UPS/VCM cannot reduce the error beyond a certain threshold, which is a classical behavior of a biased algorithm. Here we should note that this experiment setup uses independent runs to expose bias behavior, in which UPS/VCM has to restart at a fixed initial radius when generating each new image. This should not be con- fused with the parallel runs used in practical UPS/VCM rendering, where participating machines are not fully independent – they have to communicate initially to give each a different initial radius. Contribution Relative and Overhead. Compared with UPS/VCM, our UPG formulation does not need a progressively shrinking radius d. Therefore, in UPG the relative contribution of photon mapping does not diminish as the number of iterations increases. This is illustrated in Fig. 8. Consequently, we have more freedom when fine-tuning d , as it is purely an efficiency issue and is free of convergence-imposed constraints. In Fig. 14, d is chosen to compare with UPS/VCM at a similar level of variance. Compared with traditional photon mapping, UPG introduces two sources of overhead. The first source is the shadow ray traced be- tween z t ′ − 1 and ty s , plus the follow-up BSDF evaluations. The sec- ond source is the Bernoulli trials. The average number of Bernoulli trials per connection is approximately the same across our test scenes, ranging from 2.60 to 2.87.  As illustrated in Fig. 9, most pixels of Fig. 1 require fewer than 3 trials, indicating that the angu- lar bound reduces the complexity of Bernoulli trials from O(n 2 ) to O(n) in practice. Although there are a few long-tail cases triggered at a very low probability, the corresponding performance overhead is acceptable when amortized over the entire image. We quantitatively evaluate the UPG overhead using the kitchen scene in Fig. 8 as an example.  Compared with UPS/VCM, UPG takes approximately 4× time to process each photon. However, as shown in Fig. 8 (c)-(f), the shrinking radius in UPS/VCM gradu- ally reduces the effectiveness of photon mapping during each iter- ation, gradually shifting the dominance to the less efficient BDPT. Here the initial radius of UPS/VCM is 1.5× greater than the ra- dius of UPG. As shown on the right chart of Fig. 8, this allows UPS/VCM to generate light paths at a significantly higher speed than UPG+BDPT in the first few minutes, which leads to the lower estimation variance in Fig. 8 (a). As the rendering progresses, how- ever, the sampling speed of UPS/VCM drops rapidly. On the other hand, UPG+BDPT maintains a constant speed and takes the lead after approximately 2 hours.  Even well before the turning point, the cumulative effect of having more PM samples already allows UPG+BDPT to provide superior rendering quality, as shown in Fig. 8 (b).  This indicates that photon gathering is a more effec- tive sampling technique than BDPT for the kitchen scene, which is dominated by indirect illumination. Although UPG+BDPT maintains a stable sampling speed, the speed still depends on d and some fine tuning is required for optimal rendering efficiency. A larger d generates more connections while potentially consuming more computation time on low-contribution light paths.  A comparison between different values of d is shown in Fig. 10. While our method is not sensitive to small changes of d , excessively large or small values can still be problematic. Comparisons with Biased Alternatives. Fig. 14 compares our UPG with state-of-the-art biased alternatives. As shown in the fig- ures and zoom-ins, the UPS/VCM results exhibit bias artifacts from Figure 11:  Bias artifacts cannot be completely removed without our unbiased probability reciprocal estimation.  (a)-(b) are ren- dered with UPS/VCM. (b) uses the correct fc in place of fc M but the wall corner in the living room scene.  In the kitchen scene with many glossy reflections, the kettle reflects a highlight that does not exist in the reference result. Unbiased connection probability reciprocal estimation is critical when removing bias artifacts. Simply replacing f ℓ M with f c re- duces the artifacts of incorrect evaluations of visibility and shadow. However, the incorrect probability term causes legitimate light paths to become under-weighted, which manifests as artificial dark- ening at sharp features or increased estimation variance, as illus- trated in Fig. 11 (b). Such artifacts remain in the image even after an hour of rendering, as illustrated in Fig. 11 (c). It requires the un- biased probability reciprocal in UPG to completely eliminate such artifacts while keeping the low variance estimation of photon gath- ering, as illustrated in Fig. 11 (d). Comparisons with Unbiased Alternatives. We also compare UPG+BDPT with the state-of-the-art unbiased alternative of MMLT in Fig. 14. For the same computation time, MMLT can also produce a low variance estimation. However, our method does not produce spike artifacts even with challenging paths, and generates less noise for deeper indirect lighting paths such as the living room counter. That is an advantage inherited from photon mapping. We also compare our method with particle-guided BDPT [Vorba et al. 2014].  As shown in Fig. 12, UPG+BDPT (c) outperforms guided BDPT (b) in general. We further made a preliminary attempt of using guided BDPT for sub-path sampling in UPG+BDPT. Al- though our unoptimized implementation is not efficient enough for a same-time comparison yet, it demonstrates some improvements in (d) when compared with plain UPG+BDPT at the same number of paths per pixel. Although our method does not have the low frequency filtering effect of traditional density estimation techniques, UPG still per- forms well for a variety of light paths including many-bounce dif- fuse inter-reflections, generating a smoother result when compared with other unbiased methods. This indicates that the fundamental efficiency of photon gathering goes beyond its filtering effect, and is not limited to caustics paths. Comparison between particle-guided BDPT and
7
Discussion and Future Work
Limitations. The most important limitation of our unbiased pho- ton gathering is the Bernoulli trial overhead, which puts our method at a disadvantage for shallow paths that are already easy to sample. As shown in Table 1, our method is inefficient for direct lighting paths due to the higher per-sample cost.  This is the low variance problem discussed by Veach [1998]. If a specific type of light trans- port has already been well-sampled by existing techniques, com- bining with more techniques would only introduce additional vari- ance. We currently work around this issue by heuristically setting a lower sampling rate for shallow paths and disabling direct lighting paths entirely in our method. This is consistent with the practical usage of PM methods, which often recommend excluding direct lighting [Jensen 2001]. Table 1: RMSE comparison for different path depths. This experi- ment is performed on the stairs scene in Fig. 14.
UPG cannot handle strict specular-diffuse-specular (SDS) paths, where the BSDF terms are Dirac delta functions. This is an open problem for unbiased path sampling which we do not solve. Exist- ing biased estimations such as UPS/VCM approximate such specu- lar bounces by smoothing out the sampled specular photons without fully generating the relevant light paths in a physically based man- ner, which cannot be easily adapted to unbiased rendering. Another issue is that UPG is inefficient for connections with specu- lar bounces. The problem is more pronounced for refraction effects where typical light paths have more specular bounces than non- specular bounces. If one of the sub-path endpoints to be connected is on a specular surface, the corresponding delta function BSDF evaluates to zero, which prevents the connected path from making any contribution.  On the other hand, non-SDS specular paths can still be generated by making a connection at a pair of non-specular vertices. In addition, UPG does handle highly glossy materials efficiently, which can approximate specular effects to a level of sharpness com- parable with UPS/VCM. Fig. 13 compares specular caustics gen- Figure 13:  Caustics rendering.  All results are from 10 minutes of rendering except for the reference.  (a) specular rings rendered with BDPT. (b)-(c) specular rings rendered with UPS/VCM. The initial radius in (c) is a half of the initial radius in (b). (d) Highly glossy rings rendered with UPG+BDPT, whose gathering radius is the same as the initial radius in (b). (c) and (d) look similar but (d) is less noisy. (b) is more blurry compared with (c) and (d). The full images are provided in the supplementary material. erated using a variety of methods.  (a) is the reference generated using BDPT rendering, which is still noisy but illustrates the true "sharpness" of the specular caustics.  In (b), the biased UPS/VCM produces a smooth-but-blurry result.  A smaller initial radius im- proves the sharpness at the cost of a noisy image, as illustrated in Fig. 13 (c). Replacing the specular input with a highly glossy BSDF, UPG+BDPT in Fig. 13 (d) produces a sharpness similar to Fig. 13 (c) but is much less noisy. Discussion. Accuracy and efficiency are the endless pursuit of light transport simulation.  Recently, plenty of efficient methods have been proposed but accompanied with biases.  The error of Monte Carlo estimate comes from both variance and bias.  When variance is significantly reduced to a low level, the error naturally results from the bias.  Variance and bias are like two sides of a coin, without considering each we cannot have a comprehensive un- derstanding.  Biased-but-consistent methods like PPM/SPPM and UPS/VCM strike a nice balance by providing a low variance ini- tially and progressively diminishing the bias to zero.  UPG takes another approach for which we clarify the sources of bias and re- move them at a mathematical level. The elimination of bias restricts the estimation error to variance only and makes our method natively progressive. This is practically significant as parameters like d be- come significantly simpler to tune than in biased methods. Photon mapping has been widely used for efficient light trans- port simulation for decades.  Extensions like PPM/SPPM and UPS/VCM pushed the state-of-the-art for high quality rendering. However, photon mapping has always been considered a biased al- gorithm from the time of its introduction, with speculation that the method is intrinsically biased. UPG takes a critical step to disprove this speculation. In doing so, we decouple photon mapping from traditional biased evaluations, which fully integrates photon map- ping into the family of path sampling techniques by allowing pho- tons to be evaluated in an unbiased manner. Future Work.  It would be interesting to extend our method to vol- ume rendering where the gathering neighborhood or virtual lights can take more flexible shapes such as rays and beams [Sun et al. 2010; Jarosz et al. 2011; Novák et al. 2012; Křivánek et al. 2014]. The radius can be also chosen adaptively and dynamically for better performance. Another interesting direction is to extend the Russian roulette formulation to light path construction mechanisms other than nearest neighborhood search, such as VPL with the probabilis- tic shadow test. Finally, our current MCMC integration is rather straightforward, and more advanced solutions could be developed.
Acknowledgements
We would like to thank the anonymous reviewers for their helpful comments, and Steve Lin for proofreading this paper. This research was supported in part by NSFC (No. 61272305 and No. 61472352) and the National Program for Special Support of Eminent Profes- sionals of China. Bekaert, P., Cools, R., and Slussalek, P. 2003. A custom designed density estimator for light transport. Tech. Rep. MPI-I- 2003-4-004, Max-Planck Institut fr Informatik (Sarrebruck, DE), Saarbrcken. Booth, T. E.  2007.  Unbiased Monte Carlo estimation of the reciprocal of an integral. Nuclear Science and Engineering 156 , 3,403–407. CHANDLER, D. 1987. Introduction to Modern Statistical Mechan- ics. Oxford University Press. CLINE, D., TALBOT, J., AND EGBERT, P. 2005. Energy redistri- bution path tracing. ACM Trans. Graph. (Proc. SIGGRAPH) 24, 3 (July), 1186–1195. D AVIDOVIČ, T., K ŘIVÁNEK, J., H AŠAN, M., S LUSALLEK, P., AND BALA, K. 2010. Combining global and local virtual lights for detailed glossy illumination. ACM Trans. Graph. (Proc. SIG- GRAPH [Asia]) 29, 6 (Dec.), 143:1–143:8. GEORGIEV, I., K ŘIVÁNEK, J.,  D AVIDOVIČ, T., AND Light transport simulation with SLUSALLEK, P. 2012. vertex connection and merging.  ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 31, 6 (Nov.), 192:1–192:10. H achisuka, T., and Jensen, H. W. 2009. Stochastic progres- sive photon mapping.  ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 28, 5 (Dec.), 141:1–141:8. H ACHISUKA, T., AND J ENSEN, H. W.  2011.  Robust adaptive photon tracing using photon path visibility. ACM Trans. Graph. 30 , 5 (Oct.), 114:1–114:11. H ACHISUKA, T., O GAKI, S., AND J ENSEN, H. W. 2008. Progres- sive photon mapping.  ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 27, 5 (Dec.), 130:1–130:8. H ACHISUKA , T., P ANTALEONI , J., AND J ENSEN , H. W. 2012. A path space extension for robust light transport simulation. ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 31 , 6 (Nov.), 191:1– 191:10. H ACHISUKA , T., K APLANYAN , A. S., AND D ACHSBACHER , C. 2014.  Multiplexed Metropolis light transport.  ACM Trans. Graph. (Proc. SIGGRAPH) 33 , 4 (July), 100:1–100:10. HAŠAN, M., KŘIVÁNEK, J., WALTER, B., AND BALA, K. 2009. Virtual spherical lights for many-light rendering of glossy scenes.  ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 28 , 5 (Dec.), 143:1–143:6. Jakob, W., and Marschner, S.  2012. Manifold exploration: A Markov chain Monte Carlo technique for rendering scenes with difficult specular transport. ACM Trans. Graph. (Proc. SIG- GRAPH) 31 , 4 (July), 58:1–58:13. Jakob, W., 2010. Mitsuba renderer. http://www.mitsuba- renderer.org. Jarosz, W., Nowrouzezahrai, D., Thomas, R., Sloan, P.- P., AND ZWICKER, M. 2011. Progressive photon beams. ACM Trans. Graph. (Proc. SIGGRAPH [Asia]) 30 , 6 (Dec.), 181:1– 181:12. J ENSEN, H. W.  2001.  Realistic Image Synthesis Using Photon Mapping. A. K. Peters, Ltd., Natick, MA, USA. K AJIYA, J. T. 1986. The rendering equation. SIGGRAPH Comput. Graph. 20 , 4 (Aug.), 143–150. K APLANYAN , A. S., AND D ACHSBACHER , C.  2013.  Adaptive progressive photon mapping.  ACM Trans. Graph. 32 , 2 (Apr.), 16:1–16:13. K APLANYAN , A. S., AND D ACHSBACHER , C. 2013. Path space regularization for holistic and robust light transport.  Computer Graphics Forum (Proc. of Eurographics 2013) 32, 2, 63–72. K APLANYAN, A. S., H ANIKA, J., AND D ACHSBACHER, C. 2014. The natural-constraint representation of the path space for effi- cient light transport simulation. ACM Trans. Graph. (Proc. SIG- GRAPH) 33 , 4 (July), 102:1–102:13. K ELEMEN , C., S ZIRMAY - K ALOS , L., A NTAL , G., AND C SONKA , F.  2002.  A simple and robust mutation strategy for the Metropolis light transport algorithm. Comput. Graph. Forum 21, 3, 531–540. Keller, A.  1997.  Instant radiosity.  In Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Tech- niques , ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, SIGGRAPH '97, 49–56. Kitaoka, S., Kitamura, Y., and Kishino, F. 2009. Replica exchange light transport. Comput. Graph. Forum 28 , 8, 2330– 2342. K NAUS , C., AND ZWICKER , M. 2011. Progressive photon map- ping: A probabilistic approach. ACM Trans. Graph. 30 , 3 (May), 25:1–25:13. K OLLIG, T., AND K ELLER , A. 2004. Illumination in the presence of weak singularities.  In Monte Carlo and Quasi-Monte Carlo Methods , 245–257. K ŘIVÁNEK, J., G EORGIEV, I., H ACHISUKA, T., V ÉVODA, P., Šik, M., Nowrouzezahrai, D., and Jarosz, W.  2014. Unifying points, beams, and paths in volumetric light transport simulation. ACM Trans. Graph. (Proc. SIGGRAPH) 33 , 4 (July), 103:1–103:13. L AFORTUNE, E. P., AND W ILLEMS, Y. D.  1993.  Bi-directional path tracing.  In Proceedings of CompuGraphics , vol. 93, 145– 153. M ETROPOLIS, N., R OSENBLUTH, A. W., R OSENBLUTH, M. N., Teller, A. H., and Teller, E. 1953. Equation of state cal- culations by fast computing machines. The Journal of Chemical Physics 21 , 6 (June). Novák, J., Nowrouzezahrai, D., Dachsbacher, C., and J AROSZ, W. 2012. Virtual ray lights for rendering scenes with participating media. ACM Trans. Graph. (Proc. SIGGRAPH) 31 , 4 (July), 60:1–60:11. S EGOVIA, B., IEHL, J. C., M ITANCHEY, R., AND P É ROCHE, B. 2006. Bidirectional instant radiosity. In Proceedings of the 17th Eurographics Conference on Rendering Techniques , Eurograph- ics Association, Aire-la-Ville, Switzerland, Switzerland, EGSR '06, 389–397. S UN , X., Z HOU , K., L IN , S., AND G UO , B.  2010.  Line space gathering for single scattering in large scenes.  ACM Trans. Graph. (Proc. SIGGRAPH) 29, 4 (July), 54:1–54:8. Tokuyoshi, Y.  2009.  Photon density estimation using multiple importance sampling. In ACM SIGGRAPH ASIA 2009 Posters , ACM, New York, NY, USA, SIGGRAPH ASIA '09, 37:1–37:1. V EACH, E., AND G UIBAS, L. J. 1994. Bidirectional estimators for light transport.  In Proceedings of the Fifth Eurographics Work- shop on Rendering , Eurographics, 147–162. V EACH, E., AND G UIBAS, L. J.  1997.  Metropolis light trans- port. In Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques , ACM Press/Addison- Wesley Publishing Co., New York, NY, USA, SIGGRAPH '97, 65–76. V EACH, E. 1998. Robust Monte Carlo Methods for Light Transport Simulation . PhD thesis, Stanford, CA, USA. AAI9837162. Vorba, J., Karlík, O., Šik, M., Ritschel, T., and K ŘIVÁNEK, J.  2014.  On-line learning of parametric mixture models for light transport simulation. ACM Trans. Graph. (Proc. SIGGRAPH) 33 , 4 (July), 101:1–101:11. V orba, J. 2011. Bidirectional photon mapping. In Proc. of the Central European Seminar on Computer Graphics (CESCG 11). W ALTER, B., F ERNANDEZ, S., A RBREE, A., B ALA, K., D ONIKIAN , M., AND G REENBERG , D. P.  2005.  Lightcuts: A scalable approach to illumination. ACM Trans. Graph. (Proc. SIGGRAPH) 24, 3 (July), 1098–1107. W ALTER, B., A RBREE, A., B ALA, K., AND G REENBERG, D. P. 2006.  Multidimensional lightcuts.  ACM Trans. Graph. (Proc. SIGGRAPH) 25, 3 (July), 1081–1088. W ALTER , B., K HUNGURN , P., AND B ALA , K. 2012. Bidirectional lightcuts.  ACM Trans. Graph. (Proc. SIGGRAPH) 31 , 4 (July), 59:1–59:11. A The angular bound introduced in Sec. 3.4 also reduces the estima- tion variance of the probability reciprocal r .  Without the angular bound, the estimation variance can be computed based on Eq. (11) and Eq. (12).  Since all relevant terms are computed for the same path X s ′ , t ′ − 1 , we drop the function notations for brevity.
The angular bound factors p out of N as N b /p b , the variance of which can be also evaluated based on Eq. (2
Note that if we ignore the shared ( p c ) 2 term, the remaining 1 − p c and 1 − p c /p b terms are the respective rejection probabilities of the Russian roulette process.  Applying a tight angular bound usually allows most sampled tentative rays to hit the gathering neighbor- hood successfully. This reduces the rejection probability 1 − p c /p o , and by extension the estimation variance, to nearly zero. In prac- tice, the bound is constrained by the analytical integration require- ment and is not always tight. We found that without bounding, the original rejection probability 1 − p c is almost 1 as p c is usually in required Bernoulli trials. With bounding, the probability 1 − p c /p o can be reduced to less than 0.4, which corresponds to a more than 60% reduction of variance. Given two sub-paths y 1 · · · y 2 ⁄  and z 2 · · · z 1 paired by photon gath- ering, generating a full path requires removing one of the sub-path endpoints (i.e., y s ′ or z t ′ ). There are two candidates for removal, which would generate two different paths Xs'-1,t' and Xs',t'-1. However, the two path candidates are spatially close and they are unlikely to make a significantly different contribution to the final image. Therefore, as a cost-motivated optimization, we only evalu- ate one of them and discard the other. Since the bottleneck in evaluating a path is the Bernoulli trials, we do not want the selection criterion to depend on p c . In addition, we prefer a deterministic criterion, as correlating the final Monte Carlo estimate with dependent random variables may bias the result. Fi- nally, MIS-based theories do not apply here, as fundamentally this is a decision between different paths, not a decision between dif- ferent sampling techniques creating the same path. Therefore, we elect to make the choice heuristically. Note that y s ′ is importance-sampled from f s ( y s ′ − 2 → y s ′ − 1 → y α ′ ) and z ′ , is importance-sampled from f s ( z t ′ − 2  →  z t ′ − 1  → z ν ).  Intuitively, removing one of the two endpoints corresponds to reconnecting the corresponding importance-sampled edge to a slightly shifted location. One could formulate the resulting change as
Since a less glossy BSDF is less likely to result in a significant change in C ˙ RR , we heuristically test the surface properties at y s ′ − 1 and z t ′ − 1 , then we remove the endpoint corresponding to the sur- face with the lower glossiness. If y s ′ − 1 is the light source, we as- sign an effective glossiness based on its type, which is positive in- finity for directional and point sources, or a low value comparable to diffuse surfaces for area light sources. If z t ′ _ 1 is the camera, we assign a high effective glossiness corresponding to a surface that produces mirror-like reflections. As shown in Table 2, our heuris- tic significantly reduced the overall variance level.  Note that the high effective glossiness for the camera is required to exclude paths that directly connect the camera to a photon on a diffuse surface, which is the main source of error for X s ′ t ′ − 1 -type connections. Also note that the probability to sample the abandoned path is still non-zero because it can be generated by connecting y 1 · · · y s ′ − 1 z ι ′ and y ্ ′ z t ′ − 1 Table 2: RMSE comparison for different connection choices at the same computational time. The experiment is performed on the liv- ing room scene in Fig. 14. Note that evaluating both connections produced a worse result due to the computational cost wasted on less effective candidates.
C In this section, we briefly introduce the unbiased algorithm for estimating the reciprocal of a general integral proposed by Booth [2007], and discuss its relationship with our method pro- posed in Sec. 3.4. Given a non-negative function f(x) defined in a domain Ω, Booth [2007] aimed to evaluate the reciprocal of its integral
Assuming f (x)dx ∈ (0, 1), Eq. (29) can be evaluated as a series expansion:
Reorganizing the infinite series and taking an independent Monte Carlo estimation for each occurrence of q , we have
where {x j } is an infinite set of stochastic samples in Ω , and we assume ʃ ς 1d x = 1. The series can be terminated using Russian roulette once Π j = 1 ˆ θ j falls below a pre-specified threshold. However, in our case, f is the binary outcome of a Bernoulli trial, and ĝ j always takes a value of 0 or 1. Therefore, the series in Eq. (31) degenerates into a counting process.  Assuming the series of Bernoulli trials succeeds for the first time after N tries, we have:
which is identical to our algorithm in Sec. 3.4.