Modified filtered importance sampling for virtual spherical Gaussian lights

This paper proposes a modification of the filtered importance sampling method, and improves the quality of virtual spherical Gaussian light (VSGL)-based real-time glossy indirect illumination using this modification. The original filtered importance sampling method produces large overlaps of and gaps between filtering kernels for high-frequency probability density functions (PDFs). This is because the size of the filtering kernel is determined using the PDF at the sampled center of the kernel. To reduce those overlaps and gaps, this paper determines the kernel size using the integral of the PDF within the filtering kernel. Our key insight is that these integrals are approximately constant, if kernel centers are sampled using stratified sampling. Therefore, an appropriate kernel size can be obtained by solving this integral equation. Using the proposed kernel size for filtered importance sampling-based VSGL generation, undesirable artifacts are significantly reduced with a negligibly small overhead.


Introduction
The filtered importance sampling method [1] is a variance reduction technique of Monte Carlo integration often used for real-time or interactive rendering, which uses filtering kernels instead of sample points.This paper proposes a modification of filtered importance sampling, and improves the quality of virtual spherical Gaussian light (VSGL) [2] based real-time glossy indirect illumination using this modification.The original filtered importance sampling first samples the center of each filtering kernel according to a probability density function (PDF), and then determines the size of each filtering kernel using the PDF at the sampled center.However, this kernel size determination produces large overlaps of and gaps between filtering kernels for high-frequency PDFs (Fig. 1).This is because the kernel size can be too large when the sampled center is at a local minimum of a high-frequency PDF.Therefore, this paper introduces an appropriate kernel size of filtered importance sampling to reduce these overlaps and gaps.
One effective application of our method is generation of VSGLs using reflective shadow maps [3].Reflective shadow map-based global illumination is well established for real-time rendering.However, stochastic sampling of virtual point lights (VPLs) [4] (i.e., texels of reflective shadow maps which represent one-bounce light subpaths) produces noticeable variance especially for glossy interreflections.To reduce this variance, VSGLs were introduced recently.A VSGL approximates a cluster of VPLs using a Gaussian-based representation.Thanks to this representation, the distribution of VPLs can be filtered with a simple summation operation (e.g., mipmaping).In addition, this representation has an analytical solution of the rendering integral for each VSGL.Therefore, if VSGLs are generated from reflective shadow maps inexpensively, we are able to render one-bounce glossy indirect illumination at real-time frame rates.
Tokuyoshi [2] sampled VPL clusters as VSGLs from a reflective shadow map based on filtered importance sampling to achieve real-time frame rates.However, while this approach is simpler and faster than k-means-based VPL clustering [5,6], it does induce flickering structured artifacts due to the previously mentioned overlaps and gaps.This problem is noticeable when a bidirectional reflective shadow mapping method [7] is used to build the PDF.This is because the bidirectional reflective shadow mapping method produces a dynamic and high-frequency PDF.Using our kernel size, we are able to reduce flickering artifacts significantly for such a high-frequency PDF.
The contributions of our work are as follows: • An appropriate kernel size of filtered importance sampling is introduced to reduce undesirable overlaps of and gaps between filtering kernels.• For image-based PDFs, the above kernel size is computed using a simple numerical approach with negligibly small overhead.• Using the proposed filtered importance sampling method, flickering artifacts are reduced for VSGLbased real-time glossy indirect illumination.

Sampling with pre-integration
Sampling pre-integrated values is often used for image-based lighting.Structured importance sampling [8] stratifies samples hierarchically, and then the illumination is pre-integrated within each stratum.Debevec [9] subdivided an environment map into regions of equal energy using a median cut algorithm.This method approximates each region with a directional light source whose color is the sum of pixel values within the region.Filtered importance sampling [1] is introduced for glossy materials under environment maps.This technique samples prefiltered values using a mipmap, thus it performs at real-time frame rates.This paper modifies the kernel size of filtered importance sampling to improve the sampling quality.

Real-time global illumination
In this paper, we apply the proposed filtered importance sampling technique to real-time global illumination.Interactive global illumination algorithms were surveyed by Ritschel et al. [10].For a comprehensive survey of VPL-based rendering, we refer the readers to Dachsbacher et al. [11].Here we pay attention only to the most relevant works.
Virtual point lights (VPLs) are often used for representing indirect illumination [4].For realtime rendering, single-bounce VPLs lit from point or directional lights can be generated by using reflective shadow maps [3].Hundreds or thousands of VPLs are often resampled from a reflective shadow map by hierarchical sample warping [12] according to a mipmapped image-based PDF [13].
To generate shadow maps for so many lights, Ritschel et al. [14] proposed imperfect shadow maps.In their follow-up paper [7], a bidirectional reflective shadow mapping method was introduced to estimate a viewdependent importance for the image-based PDF.This method roughly computes the contribution of light paths from the eye to each reflective shadow map texel, and thus generates dynamic and highfrequency PDFs.Therefore, it can increase flickering artifacts, even though the total variance can be reduced.To take such dynamic PDFs into account, Barák et al. [15] introduced a temporally coherent sampling technique based on the Metropolis-Hastings algorithm for static light sources.They also proposed tessellation-based imperfect shadow maps to reduce memory usage.While the original VPL method is theoretically unbiased, variance is visible as spiky artifacts especially for glossy materials [16].To avoid this problem, VPLs are often clustered and then represented using a smaller number of area lights for interactive rendering.
Area light approximation via VPL clustering.Dong et al. [5] clustered VPLs using the k-means algorithm, and then approximated visibilities of VPLs using a soft shadow map for each cluster.Prutkin et al. [6] clustered texels of a reflective shadow map based on k-means similar to Dong et al., while they approximated the clusters with area lights for analytical radiance evaluation.They sampled cluster centers using the bidirectional reflective shadow mapping method to improve the quality.Luksch et al. [17] clustered VPLs using a kd-tree to generate virtual polygon lights to update light maps.These virtual area lights were evaluated using analytical form factors which cannot reduce variance caused by high-frequency bidirectional reflectance distribution functions (BRDFs).To allow for all-frequency BRDFs, recently spherical Gaussians have been used.
Spherical Gaussians [18,19] are often used for approximating the rendering of various types of materials under environment maps or area lights [20][21][22][23][24].This is because they have closed-form solutions for the integral, product, and product integral, which are fundamental operations to evaluate rendering integrals.Hence, all-frequency materials can be rendered efficiently.To represent static environment maps using spherical Gaussians, they have been fitted in preprocessing.For dynamic indirect illumination performed at near-interactive frame rates, Xu et al. [25] approximated the outgoing radiance using spherical Gaussians for each triangle primitive lit from distant light sources.Virtual spherical Gaussian lights (VSGLs) [2] were introduced to approximate a set of VPLs.To generate thousands of VSGLs at real-time frame rates, a filtered importance sampling-based approach was used with mipmapped reflective shadow maps.To generate a few VSGLs for more time-sensitive applications such as video games, the total value of all the reflective shadow map texels were computed instead of filtered importance sampling [26].These methods can render caustics, unlike eye-path tracingbased methods (e.g., cone tracing [27,28]).
However, filtered importance sampling-based VSGL generation induces a flickering error for highfrequency PDFs due to inappropriate kernel sizes.In this paper, we introduce an appropriate kernel size of filtered importance sampling.

Filtered importance sampling
Filtered importance sampling can be used when the integrand has a 2D image f (x), where x ∈ [0, 1] 2 is the image-space position.This method first samples each kernel center x i ∈ [0, 1] 2 according to a PDF p(x), and then a filtered value of f (x) is used as each sample value instead of f (x i ).This filtered value is given by a pre-filtered mipmap as follows: where g((x − x i )/s i ) is the unnormalized filtering kernel which has a fixed maximum, s i is the kernel size, a i is the filtering area (i.e., normalization factor) given by is the mipmapped value of f (x) at mip level l i .Let M be the number of texels of f (x), then the filtering area a i is also written as a function of mip level l i : . Křivánek and Colbert [1] determined mip level l i by representing this filtering area using the inverse of the density at sampled center x i as where N is the number of samples, and K is a user-specified parameter to tweak the kernel size (Křivánek and Colbert used K = 4).However, this mip level determination is sensitive to the sampled center x i .When x i is at a local minimum of a high-frequency PDF, the filtering kernel can be too large.Conversely, the filtering kernel can also be too small when x i is at a local maximum.Therefore, undesirable overlaps of and gaps between filtering kernels can be produced.

Our filtering kernels
This paper introduces an appropriate kernel size to reduce overlaps of and gaps between filtering kernels for filtered importance sampling.Sampling according to a PDF is done by computing the inverse cumulative distribution function (CDF) of the PDF.
As shown in Fig. 2, a sampling interval of the vertical axis is the integral of the PDF within each filtering kernel.Therefore, if kernel centers are sampled using stratified sampling, this integral is almost 1 N .Hence, an appropriate kernel size s i is obtained by solving the following integral equation: Since the left side is monotonically increasing with respect to the kernel size, we can obtain the kernel size using a bisection method.When PDF p(x) is given by a 2D image, we can use the mipmap of this PDF, which is also used for sampling x i via hierarchical sample warping [12].Using this mipmap, Eq. ( 2) is rewritten as where p( the mipmapped value of p(x).In this paper, l i is calculated using the bisection method with an iteration count of 12.

Application to virtual spherical Gaussian lights (VSGLs)
In this paper, we demonstrate generation of VSGLs as an effective application of our filtered sampling.A VSGL represents the positional distribution and total radiant intensity of VPLs using a Gaussian and spherical Gaussians, respectively.Since spherical Gaussians have closed-form solutions to evaluate rendering integrals, all-frequency illumination is computed analytically for each VSGL.The VSGL algorithm is composed of the following five phases: reflective shadow map rendering, PDF building, VSGL generation, shadow map rendering, and shading.This paper improves only on the VSGL generation phase.For the detail of VSGLs, please refer to Appendixes A, B, and C.

Mipmap-based VSGL generation
To generate VSGLs, VPLs are first clustered.Then VPL powers are summed and VPL distributions are averaged for each cluster (Fig. 3).To represent VPL distributions with a Gaussian and spherical Gaussians, weighted averages of emission directions, VPL positions, and squared VPL positions weighted by each VPL power are required (for the detail, please refer to Appendix B).Therefore, a reflective shadow map to store the above VPL power and weighted distribution parameters is generated, and then they are mipmapped to approximately obtain total texel values (Fig. 4).Let f (x) be the reflective shadow map, then the total texel value within the ith VPL cluster centered at x i is approximated using the mipmap f (x i , l i ) as follows: where filtering kernel g((x − x i )/s i ) represents the VPL cluster.To sample the kernel center x i and mip level l i , a filtered importance sampling-based approach can be used.The kernel center x i is sampled according to a dynamic and high-frequency view-dependent PDF p(x) given by the bidirectional reflective shadow mapping method.Tokuyoshi [2] determined l i using Eq. ( 1) with K = 1 according to the previous filtered importance sampling.Using this mip level determination, the total value of reflective shadow map texels within each cluster is given by However, since the numerator is filtered while the denominator is not, this sampling method can induce an intensive error due to overlaps of and gaps between filtering kernels.

VSGL generation using our filtering kernels
To obtain an appropriate mip level l i , this paper employs Eq. ( 3) instead of Eq. ( 1) for filtered importance sampling-based VSGL generation.Using this mip level l i , the total texel value within each cluster is approximated as follows: Unlike Eq. ( 4), both the numerator and denominator are filtered using the same kernel.Hence, temporal coherence is improved for a dynamic high-frequency PDF.Furthermore, the approximation error can be reduced if PDF p(x) is approximately proportional to f (x), similar to standard importance sampling.
Controlling the kernel size.For Eq. ( 5), the mip level l i affects only the filter bandwidth.Therefore, the user-specified parameter K can also be used for calculating l i in our case.This is implemented using ).Using K > 1, the temporal coherence is improved, though an overblurring error is induced.This overblurring error is reduced by increasing the number of samples N , similar to the original filtered importance sampling.

Experimental results
Here we present rendering results using 1024 VSGLs generated using our filtered importance sampling with K = 1 on an NVIDIA R GeForce R GTX TM 970 GPU.The frame buffer and reflective shadow map (RSM) resolutions are 1920×1088 and 512 2 , respectively.A tessellation-based imperfect shadow map [15] of resolution 64 2 is employed to evaluate the visibility of each VSGL.To estimate a viewdependent PDF on the reflective shadow map using the bidirectional reflective shadow mapping method, 2048 VSGLs without shadow maps are generated on the G-buffer.For the PDF on the G-buffer, reflectance is used.To perform stratified sampling, the Fibonacci lattice point set using a golden ratio approximation [29] is employed as a quasi-random number.For comparison, this paper uses k-means clustering using 2D image space and 3D world space.In these k-means-based approaches, once clusters are assigned to all the texels, those texels are sorted by cluster ID.Then, to compute the total value of clustered texels, a thread is dispatched for each cluster similar to Prutkin et al. [6].For implementation detail, please refer to Appendix D.
Quality. Figure 5 shows rendered images using different VSGL generation methods.Using the previous kernel size (a), intense artifacts can be produced with low probability, though this sampling method is faster than the k-means-based approaches (c, d).This is because too large of a filtering kernel is produced when the sampled kernel center is at a local minimum of the PDF.On the other hand, our kernel size (b) does not produce these undesirable filtering kernels nor does it noticeably sacrifice performance.
Performance.Table 1 shows the computation time of VSGL generation both for bidirectional reflective shadow mapping (BRSM) (upper row) and final shading (lower row).Our contribution is written in red.
Although our method is a numerical approach, its overhead is a total of about five microseconds more compared to the previous filtered importance sampling-based generation.In addition, our method is about 7-9 times faster than the k-means-based approaches.The difference is significant especially for the bidirectional reflective shadow mapping method, which uses a higherresolution G-buffer (1920×1088) than the reflective shadow map (512 2 ).Table 2 shows the computation time using different PDFs.For these PDFs, the performance of k-means-based approaches is more expensive than Table 1.This is because the last pass "Sum" (which is the summation of texel values based on Prutkin et al.'s implementation) has a linear complexity with respect to the number of texels within a cluster.Comparatively, the performance of our approach is almost independent of the PDFs, because it uses pre-filtered mipmaps.Hence, the proposed method is suitable for applications which require stable performance.
Code size.Table 3 shows the code size of VSGL generation in our implementation.Filtered importance sampling-based approaches use only two compute shaders.One is for the calculation of the additional reflective shadow map (or G-buffer), and the other is for the filtered sampling of VSGLs.The difference of our method from the previous method [2] is only the mip level determination (10 lines of code).On the other hand, the k-meansbased approaches require more compute shaders Fig. 5 Rendered images using different VSGL generation methods for 331k-triangle scene (upper row) and 75k-triangle scene (lower row).When a sampled kernel center is at a local minimum of the PDF, the previous method (a) produces too bright of a VSGL with low probability.On the other hand, our method (b) does not produce such an error similar to k-means-based approaches (c, d).Aliasing artifacts on the glossy table in the upper row are the shadow acne of imperfect shadow maps.than ours.In addition, some of them are dispatched iteratively for the GPU sort.Our method is about five times fewer lines of code than the k-means-based approaches.Kernel size controlling.As shown in Fig. 6, the kernel size of our method is controllable by using the user-specified parameter K unlike k-means-based approaches.Although some illumination appearance is overblurred by using K > 1, the temporal coherence is improved.The parameter K can be used to balance illumination details and temporal coherence according to the liking of a user.

Overlaps and gaps.
Although our method reduces overlaps of and gaps between filtering kernels, they cannot be removed completely for inhomogeneous sample distributions.This problem is alleviated by using stratified sampling.
PDF.Since our method requires a given PDF, it cannot be applied to sampling strategies without the PDF (e.g., sequential Monte Carlo instant radiosity [30]).
Since our method improves only the kernel size, kernel centers can still be temporally incoherent for dynamic PDFs.Although the Metropolis-Hastings algorithm can be used for temporally coherent sampling [15], it is limited to static light sources and has a lack of stratification.This problem induces noticeable artifacts especially for caustics (Fig. 8).Therefore, this paper employs hierarchical sample warping for stratified sampling.If the temporal coherence is more than detailed illumination, K > 1 can be used for our method to improve the temporal coherence.

Conclusions
This paper improved the kernel size of filtered importance sampling to reduce overlaps of and gaps between filtering kernels.Using this modification for VSGL generation, we are able to render glossy indirect illumination with fewer artifacts than the  8 Caustics rendered using our method with the same PDF (514triangle scene).Kernel centers of the upper row and lower row are generated using hierarchical sample warping [12] and the Metropolis-Hastings-based temporally coherent sampling [15], respectively.Due to a lack of stratification, Metropolis-Hastings produces noticeable artifacts.
previous VSGL generation.The overhead of our method is about five microseconds for thousands of VSGLs on a commodity GPU.Although the filtered importance sampling-based approach cannot take into account the difference of higher-dimensional features (e.g., world position) unlike k-means-based approaches, it is simple, fast, and has stable performance.This paper has demonstrated VSGLbased dynamic glossy indirect illumination, but our method is also usable for spherical Gaussian light generation for dynamic environment maps.Since environment maps are 2D light distribution, it might be more suitable than VSGL generation.We would like to investigate its efficiency in the future.

Appendix A Spherical Gaussians
A spherical Gaussian is a type of spherical function and is represented using a Gaussian function γ with respect to a direction vector ω ∈ S 2 as follows: where ξ ∈ S 2 is the lobe axis, and λ is the lobe sharpness.ξ and 1 λ correspond to the mean and variance for the Gaussian function, respectively.The integral of a spherical Gaussian is given by A normalized spherical Gaussian G(ω, ξ, λ) A(λ) is known as the Von Mises-Fisher distribution.For VSGLs, this distribution is used for representing reflection lobes.

A.1 Spherical Gaussian approximation of reflection lobes
Diffuse lobes.For the Lambert BRDF ρ d , the diffuse reflection lobe can be approximated with a spherical Gaussian taking energy conservation into account as follows: where ω ∈ S 2 is the incoming direction, n ∈ S 2 is the surface normal at the world position y ∈ R 3 , ω, n = max(ω • n, 0) is the clamped dot product, and λ d ≈ 2 is the sharpness of the diffuse lobe which is obtained by using the least square method.
Specular lobes.For the microfacet BRDF ρ s , the specular reflection lobe is fitted with a single spherical Gaussian by using Wang et al.'s analytical approximation [22].The BRDF is separated into two factors: the unnormalized normal distribution function D(ω h ) whose maximum is one, and the rest of the factors C(ω) as follows: where ω h = ω + ω ω + ω is the halfway vector of the incoming direction and outgoing direction.Bellshaped normal distribution functions (e.g., Phong [31], Beckmann [32], and GGX [33,34] distributions) can be approximated with a spherical Gaussian as For Beckmann or GGX normal distribution functions, λ h = 2 α 2 , where α is the roughness parameter.Using spherical warping, this can be approximated with a of ω as where ξ s is the reflection vector given by ξ s = 2(ω • n)n − ω , and λ s = λ h 4|ξ s • n| .Hence, the specular lobe is approximated with the following equation: Moreover, since microfacet BRDFs mostly preserve energy for highly glossy surfaces, the specular lobe can be approximated using a normalized spherical Gaussian as follows: where R s is the specular reflectance.Anisotropic spherical Gaussians [19] are also usable in the same manner.

Appendix B Virtual spherical Gaussian lights (VSGLs)
This paper approximates a cluster of VPLs with a VSGL.For a VSGL, the total radiant intensity and positional distribution of VPLs are represented using a spherical Gaussian and isotropic Gaussian distribution respectively.This representation can be computed using a simple summation operation.

B.1 Radiant intensity
The radiant intensity of the jth VPL is given as where Φ j is the power of the jth photon emitted from the light source, ω j ∈ S 2 is the incoming direction of the photon, and n j ∈ S 2 is the surface normal at the VPL position y j ∈ R 3 , and ρ(y j , ω j , ω) is the BRDF.This paper first divides this BRDF into diffuse and specular components (i.e., ρ d and ρ s ).Then, the total radiant intensity of clustered VPLs is approximated with a single spherical Gaussian for each component by using Toksvig's filtering [35].
For ease of explanation, this subsection hereafter describes only a single BRDF component.The total radiant intensity of a VPL cluster S is represented as To compute spherical Gaussian parameters c v , ξ v , and λ v efficiently, each reflection lobe is approximated using Eq. ( 6) or Eq. ( 7) as follows: j∈S Φ j R j where R j is the reflectance, and ξ j and λ j are the axis and sharpness of the reflection lobe at the jth VPL.Then, the weighted average of the normalized spherical Gaussians weighted by the VPL power Φ j R j is approximated with a single spherical Gaussian as A(λ v ) Using Toksvig's filtering, the jth normalized spherical Gaussian is first approximately converted into its averaged direction as ξj = λ j λ j + 1 ξ j .Next, the weighted average of the directions is computed by ξv = j∈S Φ j R j ξj j∈S Φ j R j Finally, the filtered spherical Gaussian is obtained from the weighted average direction as ξ v = ξv ξv , the VSGL viewed from y p .Using Eq. ( 11), the incoming radiance is approximated with the product of two spherical Gaussians which yields a spherical Gaussian as follows: where and

C.2 Shading via product integrals of spherical Gaussians
Since the reflection lobe ρ(y p , ω p , ω) ω, n p can be approximated using spherical Gaussians and anisotropic spherical Gaussians, Eq. ( 9) can be calculated using the analytical product integral.Diffuse reflection.Using Eq. ( 6) and Eq. ( 12), the rendering integral of the diffuse component is calculated using the analytical product integral of two spherical Gaussians.This approach is efficient for a few VSGLs [26].However, a light leak error caused by the spherical Gaussian approximation of reflection lobes cannot be reduced by increasing the number of VSGLs.Unlike the secondary bounce represented by VSGLs, light leaks are noticeable at the first bounce which is more visually important.Therefore, for thousands of VSGLs, the cosine factor at the first bounce is assumed to be a constant and pulled out of the integral [22] as follows: In addition, when λ in is not small, A(λ in ) ≈ 2π λ in can be assumed [20].Therefore, diffuse reflection is inexpensively calculated using the following equation: Specular reflection.While spherical Gaussians are used for VSGLs, this paper employs an anisotropic spherical Gaussian to approximate a specular lobe at a shading point.This is because a specular lobe can be anisotropic even if it is an isotropic BRDF model, especially for shallow grazing angles.For simplicity, anisotropic spherical Gaussians are used only for the first bounce which is more visually important than the second bounce.In addition, the product integral of a spherical Gaussian and anisotropic spherical Gaussian [19] has a reasonable computation cost.An anisotropic spherical Gaussian is defined as Ǵ(ω, ξ x , ξ y , ξ z , η x , η y ) = ω, ξ z e −η x (ω•ξ x ) 2 −η y (ω•ξ y ) 2 where ξ x , ξ y , ξ z are orthonormal vectors, and η x , η y are the bandwidth parameters.Since a specular lobe is approximated with an anisotropic spherical Gaussian as ρ s (y p , ω p , ω) ω, n p ≈ C(ω) Ǵ(ω, ξ x , ξ y , ξ z , η x , η y ), the rendering integral is calculated as L s (y p , ω p ) = where ν = λ in 2 .
Algorithm 1 Mip level calculation using the bisection method lmin ← 0 lmax ← the top mip level of p for k = 1 to the user-specified iteration count do l ← (lmin + lmax)/2 if 4 l M p(x i , l) < K N then lmin ← l else lmax ← l end if end for l i ← (lmin + lmax)/2 binary radix tree [38].For densely distributed cluster centers, this tree-based search is more efficient than using a 2D uniform grid proposed by Prutkin et al. [6].Once clusters are assigned to all the texels, those texels are sorted by cluster ID.Then, to compute the total value of clustered texels, a thread is dispatched for each cluster similar to Prutkin et al. [6].Unlike Prutkin et al., we use a GPU radix sort [39] instead of bitonic sort for the high-resolution reflective shadow map and G-buffer.Although k-means clustering can be improved by updating cluster centers in an iterative fashion, we do not update iteratively in this paper.

Fig. 1
Fig.1Sampling VPL clusters from a reflective shadow map based on filtered importance sampling.Red points: kernel centers sampled according to the PDF (brightness).Orange squares: filtering kernels (i.e., VPL clusters).The previous filtered importance sampling (a) produces significant overlaps of and gaps between filtering kernels.Our method (b) reduces these overlaps and gaps.

Fig. 3
Fig.3Clustered VPLs.Each cluster is approximated with a VSGL by computing the total VPL power and averaged VPL distributions within the cluster.These operations are done by filtered sampling on the reflective shadow map.

2 Fig. 4
Fig. 4 Mipmapped reflective shadow map for VSGL generation.Average emission directions (c, d) and positions (e, f) are weighted by VPL powers (a, b).VSGLs are sampled from this reflective shadow map based on filtered importance sampling.
Fig.8 Caustics rendered using our method with the same PDF (514triangle scene).Kernel centers of the upper row and lower row are generated using hierarchical sample warping[12] and the Metropolis-Hastings-based temporally coherent sampling[15], respectively.Due to a lack of stratification, Metropolis-Hastings produces noticeable artifacts.