Anisotropic density estimation for photon mapping

Photon mapping is a widely used technique for global illumination rendering. In the density estimation step of photon mapping, the indirect radiance at a shading point is estimated through a filtering process using nearby stored photons; an isotropic filtering kernel is usually used. However, using an isotropic kernel is not always the optimal choice, especially for cases when eye paths intersect with surfaces with anisotropic BRDFs. In this paper, we propose an anisotropic filtering kernel for density estimation to handle such anisotropic eye paths. The anisotropic filtering kernel is derived from the recently introduced anisotropic spherical Gaussian representation of BRDFs. Compared to conventional photon mapping, our method is able to reduce rendering errors with negligible additional cost when rendering scenes containing anisotropic BRDFs.


Introduction
Global illumination is a long and important research direction in computer graphics.Photon mapping [1,2] has always been a widely used technique for global illumination due to its high rendering quality and good efficiency.It is a two-pass algorithm.In the first pass, a large number of photons are emitted from light sources, traced through the scene, and stored in a photon map.In the second pass, each eye ray is traced from the viewpoint into the scene until it hits a diffuse surface; the indirect radiance at this intersection point can be approximated by averaging contributions from nearby photons within an encompassing disk of a fixed radius.This step is usually referred to as density estimation.
Density estimation allows the same photons to be reused in different eye paths, and hence makes photon mapping to be more efficient than Monte Carlo ray tracing [3,4].The radius of the disk is an important parameter in density estimation, which largely affects the bias and variance of the results.Larger radius gives lower variance but higher bias, while smaller radius gives higher variance but lower bias.Instead of equally weighting all photons, smooth filtering kernels can be applied in order to further reduce bias, such as the cone filter [1] and the Epanechnikov kernel [5,6].
We focus on using photon mapping to render scenes with anisotropic BRDFs.Since density estimation can only be applied on diffuse surfaces, eye rays towards surfaces with anisotropic BRDFs need further tracing in the scene until hitting a diffuse surface.We refer to these eye paths, from viewpoint through reflections at anisotropic surfaces to density estimation points at diffuse surfaces, as anisotropic eye paths.Intuitively thinking, density estimation through anisotropic eye paths can be beneficial from anisotropic filtering kernels.However, existing kernels are mainly isotropic.One exception is photon differentials [7,8].They use an anisotropic filtering kernel for density estimation derived from ray differentials stored in photons, and it is effective in reducing bias when rendering caustics.However, they only take into consideration of the information in light paths but not eye paths in constructing filtering kernel, hence they will not be beneficial to render anisotropic eye paths.
To improve efficiency in rendering anisotropic eye paths, we propose a novel anisotropic filtering kernel for density estimation, which considers the anisotropic BRDFs on the eye path.Our method works as follows: first, we represent anisotropic BRDFs using the recently proposed anisotropic spherical Gaussians (ASGs) [9] of BRDFs and obtain anisotropic eye paths by importance sampling anisotropic BRDFs; next, we construct the anisotropic filtering kernel in the direction space according to the gradient of the ASG; after that, the final anisotropic filtering kernel is obtained by projecting it from the direction space to the tangent plane of the density estimation point.Compared to existing works, our experiments demonstrate that our anisotropic kernel yields a better accuracy in rendering anisotropic scenes without incurring additional rendering costs.

Related works
Photon mapping.Photon mapping is a popular technique of global illumination rendering, which is based on photon density estimation [10].Arvo [11] first combined density estimation with light transport simulation by introducing particle tracing algorithm.Jensen [1,2] introduced photon mapping which approximates the radiance by local photon density estimation.Although photon mapping is a biased algorithm, it is good at rendering various lighting effects like caustics.
To improve efficiency and robustness of photon mapping, Hachisuka et al. [12] proposed progressive photon mapping (PPM), which breaks the memory requirement bottleneck for storing a large number of photons and eliminates bias gradually by reducing the radius of the radiance estimate kernel progressively.Stochastic progressive photon mapping [13] is more robust to complex scenes and is capable of producing more distributed ray tracing effects, e.g., depth of field and motion blur.Hachisuka et al. [14] introduced an error estimation framework for progressive photon mapping, while Knaus and Zwicker [15] presented an asymptotic analysis of converging rates of variance and bias.Knaus and Zwicker [15] introduced a probabilistic framework of progressive photon mapping, which does not need to maintain local photon statistics and can be implemented easier.Kaplanyan and Dachsbacher [16] proposed adaptive progressive photon mapping by selecting local kernel bandwidth adaptively and achieved a higher convergence rate.
Many works have been done in order to further improve the robustness of photon mapping.Schjøth et al. [7,8] proposed photon differentials, which shape an anisotropic filtering kernel for density estimation derived from ray differentials created in the photon tracing pass and can produce better caustic quality with less photons.Spencer and Jones [17] introduced photon relaxation which redistributes the photons into a blue noise distribution.By manipulating the underlying points with feature detection and preservation, photon relaxation achieves noise reduction without increasing bias.
Recently, multiple importance sampling (MIS) [18] has been used in photon density estimation.Vorba [19] extended photon mapping into a bidirectional way and used MIS to combine light paths and eye paths of different lengths.A unified algorithm [20,21] combining bidirectional path tracing (BPT) and progressive photon mapping (PPM) is more robust to handle complex illuminations and specular-diffuse-specular paths.These MIS based algorithms are more efficient than both BPT and PPM.
Anisotropic appearance.Many real world materials are anisotropic, such as metal and wood.Anisotropic appearance exhibits changes with respect not only to the azimuthal difference between incoming and viewing directions, but also to the azimuthal angle of incoming direction.
Kajiya [22] introduced the first anisotropic BRDF model.After that, a number of parametric anisotropic BRDF models have been proposed, such as Ward's model [23] and Ashikhmin's model [24].Various models have also been proposed to handle specific anisotropic materials, such as hair [25], wood [26], and cloth [27,28].
Recently, Xu et al. [9] introduced a representation called anisotropic spherical gaussian (ASG), which is efficient and effective in representing anisotropic spherical functions, such as anisotropic BRDFs and visibilities.Furthermore, ASGs have closeform solutions for integral, multiplication, and convolution operators.Due to its effectiveness at approximating anisotropic BRDFs, in our work we use it to represent anisotropic BRDFs.

Background
Photon mapping [1,2] is a two-pass global illumination algorithm.It is good at producing illumination effects like caustics and is efficient in rendering low-variance images.In the first pass, a large number of photons are emitted from light sources, traced through the scene, and stored in a photon map.In the second pass, each eye ray is traced from the viewpoint into the scene until it hits a diffuse surface.Then the indirect radiance at this intersection point can be evaluated through density estimation.Specifically, the indirect radiance is evaluated as the average of fluxes of its N nearest photons: where k iterates over the N nearest photons, Φ k and i k are the flux and incident direction of the k-th photon respectively, o is the outgoing direction, ρ denotes the BRDF, and r is the radius of the disk that encompasses the N nearest photons.The above equation equally weights all the nearest neighboring photons.To reduce bias, we can weight each photon according to the distance from the photon to the density estimation point.Then, Eq. ( 1) can be rewritten as where d k is the distance from the k-th photon to the density estimation point, and the kernel function w(•) is a monotonic decreasing function.Different kernels have been proposed, such as the cone filter [1] and the Epanechnikov kernel [5,6], both of which have been demonstrated to be useful in reducing bias.we construct the anisotropic filtering kernel in the direction space according to the gradient of the warped ASG (Section 4.2); after that, we project the kernel in the direction space to the local coordinate system of the density estimation point P 1 to obtain the final anisotropic filtering kernel (Section 4.3); finally, we will show how this anisotropic kernel can be used in density estimation (Section 4.4).

Spherical warping
We use ASGs [9] to represent the anisotropic BRDF at P 0 .We choose ASGs as our BRDF representation since ASGs are simple to compute and have good scalability in approximating anisotropic signals.As shown in their work, commonly used parametric anisotropic BRDFs, such as Ward's model [23] and Ashikhmin's model [24], can be well approximated by one ASG.ASG based BRDF representation is based on the microfacet model [29,30].Specifically, the normal distribution function (NDF) is approximated using one ASG, and then the BRDF at a specific view slice is obtained through an ASG spherical warping operator.As shown in Fig. 1, denote the view direction as o, the reflected direction from P 0 to P 1 as r, the anisotropic BRDF ρ(r, o) is approximated by a warped ASG: where M is a smooth function that combines the shadowing term and Fresnel term; G is an ASG; x, y, z are the tangent, bi-tangent, lobe axes, respectively; λ and µ are the bandwidths for xand y-axes, respectively.Those parameters of the ASG G are computed from the NDF function (which is also an ASG) and the view direction o through the ASG warping operator (please refer to their paper [9] for more details of ASGs and the ASG warping operator), hence it is also referred to as the warped ASG.For simplicity, we denote the warped ASG G(r; [x, y, z], [λ, µ]) as G(r) for short.

Kernel construction
As shown in Fig. 1, imagining an elliptical filtering kernel is applied at a specific direction to the warped ASG.By observation, it is easy to notice that the optimal choice to preserve the structure in filtered results is to set the minor axis of the ellipse as the gradient direction and set the major axis as the tangent direction (i.e., the direction perpendicular to the gradient direction).Based on this observation, we construct our elliptical filtering kernel in the direction space according to the gradient of the warped ASG.We first compute the gradient g of the warped ASG G(r) using the formula below: where ∇G is computed by directly applying gradient operation to the warped ASG function in the 3D vector space.Since r is a direction (i.e., constrained on the unit sphere) and it is not an arbitrary 3D vector, we need to minus (∇G • r)r to obtain the true gradient to make it perpendicular to direction r.
After that, we simply set the direction of the minor axis u d as the gradient direction u d = g/ g , and set the major axis direction as the tangent direction v d = r × u d .Now, we need to determine the lengths of the minor/major axes.We obtain the axis lengths by constraining relative value changes inside the ellipse within a predefined threshold , and in our implementation we set = 0.02.Along the minor axis, we approximate the value change using first order Taylor expansion at direction r, and the minor axis length l u is obtained by where ∂G/∂u d is the directional derivative of the warped ASG function G along the minor axis direction u d .Since u d is the gradient direction, ∂G/∂u d equals to the length of gradient g .Along the major axis, since the directional derivative along it is zero, we use second order Taylor expansion to obtain the major axis length l v : where ∂ 2 G/∂v 2 d is the second order derivative along the major axis direction v d .We now obtain the minor/major axes (i.e., u and v) by combining its directions and lengths:

Planar projection
We now have obtained the ellipse for filtering in the direction space.We need to further project it to the tangent plane of the density estimation point P 1 , since the density estimation is finally performed on this plane.We denote the minor/major axes of the projected ellipse as s and t.Note that these two projected axes are not necessarily perpendicular to each other.As shown in Fig. 2, take the minor axis as an example, the projected minor axis s should satisfy two conditions: where k is a scalar coefficient and n is the normal direction of the density estimation center.
Combining the above two equations leads to We can similarly obtain the projected major axis t as Fig. 2 Illustration of planar projection.

Density estimation
After obtaining the ellipse on the tangent plane of the density estimation point P 1 , we now explain how to use it as a filtering kernel in density estimation.Specifically, we explain how to compute the weight for each photon in Eq. ( 2).First, we rescale the ellipse to make it fit the encompassing disk of the N nearest neighboring photons, i.e., rescale the ellipse to make the length of the major axis equal to the radius of the disk.The scaled minor/major axes become: where r is the radius of the encompassing disk.After that, for each photon, we denote the offset from the density estimation center to the photon position along the tangent plane as d.As shown in Fig. 3, we can express the offset d as a linear combination of the projected minor/major axes: The coefficients a and b can be easily obtained by equations below: The weight used in density estimation is computed as a Gaussian: Note that we need to normalize the sum of the weights of all photons into one when performing density estimation in Eq. ( 2).

Results
Implementation.
We implement our method on a consumer level PC with an Intel Core i7-Fig.3 Illustration of linear combination.3770 3.4 GHz CPU and 8 GB memory.In our implementation, anisotropic filtering is only applied to density estimation for anisotropic eye paths.For other types of eye paths, e.g., eye ray starting from the viewpoint directly hits a diffuse surface, our method degenerates to use isotropic filtering in density estimation instead.Note that our method focuses on improving the effectiveness of density estimation and hence is perpendicular to other photon mapping improvements such as progressive photon mapping [12,15].To obtain consistent results, we use the progressive photon mapping framework to produce final rendering results.In each iteration, we reduce the radius of density estimation with a reduction parameter α = 2/3.
Comparison.To demonstrate the efficiency of our anisotropic filtering kernel, we compare the rendering results of our anisotropic filtering kernel to those of an isotropic kernel and a constant kernel, using several scenes with anisotropic glossy materials.
In our comparison among different methods, only the filtering kernels are different while other parameters such as the image resolution, the number of photons and the number of iterations keep the same.We compare the results visually, as well as measuring root mean square errors (RMSE) to ground truth.
Figure 4 shows an anisotropic glossy buddha with anisotropic ratio of 10 inside a diffuse textured box.As shown in the enlarged images in Fig. 4, our proposed anisotropic filtering kernel produces less noise than an isotropic filtering kernel with equal rendering time.Note that the time required for constructing the anisotropic filter in our method is negligible compared to other steps in photon mapping, hence, our method won't incur additional costs.Figure 5 shows the RMSE of the rendering results of a constant kernel, an isotropic Gaussian kernel, and our anisotropic Gaussian kernel.Our method consistently gives the smallest error.
Figure 6 shows a scene with a teapot, a fork, and a spoon, and the anisotropic ratios are 10, 3.3, and 3.3, respectively.The surrounding objects are diffuse textures.Figure 7 shows a scene with anisotropic glossy pans and pots whose anisotropic ratios are 10.In both scenes, our method achieves better visual quality and lower RMSE than using isotropic kernel with equal rendering time.

Fig. 5
Fig.5RMSE curve.We plot the RMSE of the rendered images with respect to the number of iterations used in progressive photon mapping in rendering the buddha scene.Notice that our method consistently gives the smallest error.