Fu-Jun Luan,Li-Fan Wu,and Kun Xu

Anisotropic density estimation for photon mapping

Fu-Jun Luan1,Li-Fan Wu1,and Kun Xu1

Photon mapping is a widely used technique forglobalillumination rendering.In the density estimation step of photon mapping,the indirect radiance at a shading point is estimated through a fi ltering process using nearby stored photons;an isotropic fi ltering 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 fi ltering kernel for density estimation to handle such anisotropic eye paths. The anisotropic fi ltering kernel is derived from the recently introduced anisotropic spherical Gaussian representation of BRDFs.Compared to conventional photon mapping,ourmethod isabletoreduce rendering errors with negligible additional cost when rendering scenes containing anisotropic BRDFs.

photonmapping;densityestimation; anisotropic; anisotropic spherical Gaussian

1 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 di ff use surface;the indirect radianceat this intersection point can be approximated by averaging contributions from nearby photons within an encompassing disk of a fi xed radius.This step is usually referred to as density estimation.

Density estimation allows the same photons to be reused in di ff erent 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 a ff ects 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 fi ltering kernels can be applied in order to further reduce bias,such as the cone fi lter[1]and the Epanechnikov kernel[5,6].

We focus on using photon mapping to render sceneswith anisotropicBRDFs.Sincedensity estimation can only be applied on di ff use surfaces, eye rays towards surfaces with anisotropic BRDFs need further tracing in the scene until hitting a di ff use surface.We refer to these eye paths,from viewpoint through re fl ections at anisotropic surfaces to density estimation points at di ff use surfaces, as anisotropic eye paths.Intuitively thinking, density estimation through anisotropic eye paths can be bene fi cial from anisotropic fi ltering kernels. However,existing kernels are mainly isotropic.One exception is photon di ff erentials[7,8]. They use an anisotropic fi ltering kernel for density estimation derived from ray di ff erentials stored in photons, and it is e ff ective in reducing bias when rendering caustics.However,they only take into consideration of the information in light paths but not eye paths in constructing fi ltering kernel,hence they will not be bene fi cial to render anisotropic eye paths.

To improve efficiency in rendering anisotropic eye paths,we propose a novel anisotropic fi lteringkernelfordensity estimation,which considers the anisotropic BRDFs on the eye path. Our method worksasfollows: first,werepresent 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 fi ltering kernel in the direction space according to the gradient of the ASG;after that, the final anisotropic fi ltering 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.

2 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 e ff ects 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 ofphotons and eliminates bias gradually by reducing the radius ofthe radiance estimate kernel progressively.Stochastic progressive photon mapping[13]is more robust to complex scenes and is capable of producing more distributed ray tracing e ff ects,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 adaptiveprogressive photon mapping by selecting local kernel bandwidth adaptively and achieved a higher convergence rate.

Many workshave been done in orderto further improve the robustness of photon mapping. Schj?th et al.[7,8]proposed photon di ff erentials, which shape an anisotropic fi ltering kernel for density estimation derived from ray di ff erentials created in the photon tracing passand 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 bidirectionalway and used MIS to combine light paths and eye paths of di ff erent lengths.A uni fied algorithm[20,21]combining bidirectional path tracing (BPT) and progressive photon mapping(PPM)is more robust to handle complex illuminations and specular-di ff use-specular paths. These MIS based algorithms are more efficient than both BPT and PPM.

Anisotropic appearance. Many real world materialsare anisotropic, such asmetaland wood.Anisotropic appearance exhibits changes with respect not only to the azimuthal di ff erence 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 speci fi c 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 e ff ective 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 e ff ectiveness at approximating anisotropic BRDFs,in our work we use it to represent anisotropic BRDFs.

3 Background

Photon mapping [1,2]is a two-passglobal illumination algorithm.It is good at producing illumination e ff ects 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 di ff use surface.Then the indirect radiance at this intersection point can be evaluated through density estimation. Speci fi cally,the indirect radiance is evaluated as the average of fl uxes of itsNnearest photons:

wherekiterates over theNnearest photons,Φkandikare the fl ux and incident direction of thek-th photon respectively,ois the outgoing direction,ρdenotes the BRDF,andris the radius of the disk that encompasses theNnearest 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

wheredkis the distance from thek-th photon to the density estimation point,and the kernel functionw(·)is a monotonic decreasing function.Di ff erent kernels have been proposed,such as the cone fi lter[1] and the Epanechnikov kernel[5,6],both of which have been demonstrated to be useful in reducing bias.

4 Our method

As shown in Fig.1,considering a typical anisotropic eye path,starting from viewpoint to a pointP0on an anisotropic surface,then re fl ected to the density estimation pointP1on a di ff use surface,the anisotropic fi ltering kernel is computed and applied through 4 steps. First,based on the anisotropic sphericalGaussian (ASG)representation [9]of BRDFs,we obtain the 2D BRDF slice of re fl ected directionrat pointP0,and approximate it as an ASG through ASG warping[9](Section 4.1);next, we construct the anisotropic fi ltering 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 pointP1to obtain the final anisotropic fi ltering kernel(Section 4.3); finally,we will show how this anisotropic kernel can be used in density estimation(Section 4.4).

Fig.1 Illustration of anisotropic eye paths.

4.1 Spherical warping

We use ASGs[9]to represent the anisotropic BRDF atP0.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].Speci fi cally,the normal distribution function(NDF)is approximated using one ASG,and then the BRDF at a speci fi c view slice is obtained through an ASG spherical warping operator.As shown in Fig.1,denote the view direction aso,the re fl ected direction fromP0toP1asr,the anisotropic BRDFρ(r,o)is approximated by a warped ASG:

whereMis a smooth function that combines the shadowing term and Fresnel term;Gis an ASG;x,y,zare the tangent,bi-tangent,lobe axes, respectively;λandμare the bandwidths forx-andy-axes,respectively.Those parameters of the ASGGare computed from the NDF function(which is also an ASG)and the view directionothrough 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 ASGG(r;[x,y,z],[λ,μ])asG(r)for short.

4.2 Kernel construction

As shown in Fig.1,imagining an elliptical fi ltering kernel is applied at a speci fi c direction to the warped ASG.By observation,it is easy to notice that the optimal choice to preserve the structure in fi ltered 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 fi ltering kernel in the direction space according to the gradient of the warped ASG.

We first compute the gradientgof the warped ASGG(r)using the formula below:

where?Gis computed by directly applying gradient operation to the warped ASG function in the 3D vector space.Sinceris a direction(i.e.,constrained on the unit sphere)and it is not an arbitrary 3D vector,we need to minus(?G·r)rto obtain the true gradient to make it perpendicular to directionr. After that,we simply set the direction of the minor axisudas the gradient directionud=g/‖g‖,and set the major axis direction as the tangent directionvd=r×ud.

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 prede fi ned threshold∈,and in our implementation we set∈=0.02.Along the minor axis,we approximate the value change using first order Taylor expansion at directionr,and the minor axis lengthluis obtained by

where?G/?udis the directional derivative of the warped ASG functionGalong the minor axis directionud. Sinceudis the gradient direction,?G/?udequals 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 lengthlv:

whereis the second order derivative along the major axis directionvd.We now obtain the minor/major axes(i.e.,uandv)by combining its directions and lengths:

4.3 Planar projection

We now have obtained the ellipse for fi ltering in the direction space.We need to further project it to the tangent plane of the density estimation pointP1, since the density estimation is finally performed on this plane.

We denote the minor/major axes of the projected ellipse assandt.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 axissshould satisfy two conditions:

wherekisa scalarcoefficientandnisthe normal direction of the density estimation center. Combining the above two equations leads to

We can similarly obtain the projected major axistas

Fig.2 Illustration of planar projection.

4.4 Density estimation

After obtaining the ellipse on the tangent plane of the density estimation pointP1,we now explain how to use it as a fi ltering kernel in density estimation. Speci fi cally,we explain how to compute the weight for each photon in Eq.(2).First,we rescale the ellipse to make it fi t the encompassing disk of theNnearest 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:

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).

5 Results

Fig.3 Illustration of linear combination.

Implementation. We implement our method on a consumer level PC with an Intel Core i7-3770 3.4GHz CPU and 8GB memory.In our implementation,anisotropic fi ltering 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 di ff use surface, our method degenerates to use isotropic fi ltering in density estimation instead.Note that our method focuses on improving the e ff ectiveness 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 fi ltering kernel,we compare the rendering results of our anisotropic fi ltering kernel to those of an isotropic kernel and a constant kernel,using several scenes with anisotropic glossy materials. In our comparison among di ff erent methods,only the fi ltering kernels are di ff erent 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 di ff use textured box.As shown in the enlarged images in Fig.4, our proposed anisotropic fi ltering kernel produces less noise than an isotropic fi ltering kernel with equal rendering time.Note that the time required for constructing the anisotropic fi lter 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 di ff use 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.4 Anisotropic happy buddha.(a)and(b)Photon mapping results using isotropic and anisotropic density estimation kernels, respectively;(c)reference.

Fig.5 RMSE 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.

6 Conclusions

In this work,we have presented an anisotropic fi ltering kernel for density estimation.The kernel takes into consideration of the anisotropic BRDFs on the eye path.It is derived from the ASG representation ofanisotropic BRDFs.Through experiments and comparisons, our proposed anisotropic kernel is demonstrated to produce lower rendering error without incurring additional costs in rendering anisotropic scenes.

Acknowledgements

We thank the anonymous reviewers for their valuable comments. This work was supported by the National High-tech R&D Program of China(No. 2012AA011802)and the National Natural Science Foundation of China(No.61170153). Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use,distribution,and reproduction in any medium,provided the original author(s)and the source are credited.

[1]Jensen,H.W.Global illumination using photon maps. In: Proceedings of the eurographics workshop on Rendering techniques,21-30,1996.

[2]Jensen,H.W.Realistic Image Synthesis Using Photon Mapping.Natick,MA,USA:A.K.Peters,Ltd.,2001.

[3]Cook,R.L.Stochastic sampling in computer graphics.ACM Transactions on GraphicsVol.5,No.1,51-72, 1986.

[4]Cook,R.L.;Porter,T.;Carpenter,L.Distributed ray tracing.ACM SIGGRAPH Computer GraphicsVol. 18,No.3,137-145,1984.

[5]Roland,S.Bias compensation for photon maps.Computer Graphics ForumVol.22,No.4,729-742, 2003.

[6]Walter,B.J.Density estimation techniques for global illumination.Ph.D.Thesis.Cornell University,1998.

[7]Schj?th,L.Anisotropic density estimation in global illumination.Ph.D.Thesis.University of Copenhagen, 2009.

[8]Schj?th,L.;Frisvad,J.R.;Erleben,K.;Sporring, J.Photon di ff erentials.In:Proceedings of the 5th international conference on Computer graphics and interactive techniques in Australia and Southeast Asia, 179-186,2007.

[9]Xu,K.;Sun,W.-L.;Dong,Z.;Zhao,D.-Y.;Wu,R.-D.;Hu,S.-M.Anisotropic spherical Gaussians.ACM Transactions on GraphicsVol.32,No.6,Article No. 209,2013.

Fig.6 Anisotropic teapot.(a)and(b)Photon mapping results using isotropic and anisotropic density estimation kernels, respectively;(c)reference;(d)and(e)RMSE visualization.

Fig.7 Anisotropic frying pan.(a)and(b)Photon mapping results using isotropic and anisotropic density estimation kernels, respectively;(c)reference;(d)and(e)RMSE visualization.

[10]Hachisuka,T.;Jarosz,W.;Bouchard,G.;Christensen, P.;Frisvad,J.R.;Jakob,W.;Jensen,H.W.;Kaschalk, M.;Knaus,C.;Selle,A.;Spencer,B.State of the art in photon density estimation.In:ACM SIGGRAPH 2012 Courses,Article No.6,2012.

[11]Arvo,J.Backward ray tracing.In:Developments in Ray Tracing,SIGGRAPH’86 Course Notes,Vol.12, 259-263,1986.

[12]Hachisuka,T.;Ogaki,S.;Jensen,H.W.Progressive photon mapping.ACM Transactions on GraphicsVol.27,No.5,Article No.130,2008.

[13]Hachisuka T.;Jensen H.W.Stochastic progressive photon mapping.ACM Transactions on GraphicsVol. 28,No.5,Article No.141,2009

[14]Hachisuka,T.;Jarosz,W.; Jensen,H.W.A progressive error estimation framework for photon density estimation.ACM Transactions on GraphicsVol.29,No.6,Article No.144,2010.

[15]Knaus,C.;Zwicker,M.Progressive photon mapping: A probabilisticapproach.ACMTransactionson GraphicsVol.30,No.3,Article No.25,2011.

[16]Kaplanyan, A.S.; Dachsbacher, C.Adaptive progressive photon mapping.ACM Transactions on GraphicsVol.32,No.2,Article No.16,2013.

[17]Spencer, B.; Jones, M. W.Into the blue: Better caustics through photon relaxation.Computer Graphics ForumVol.28,No.2,319-328,2009.

[18]Veach, E.Robust Monte Carlo methods for light transport simulation.Ph.D.Thesis.Stanford University,1998.

[19]Vorba,J.Bidirectional photon mapping.In:The 15th Central European Seminar on Computer Graphics, Vol.2,2011.Available at http://cgg.m ff.cuni.cz/∽jaroslav/papers/2011-bdpm/vorba2011-bdpm.pdf.

[20]Georgiev,I.;Kriv′anek,J.;Davidovic,T.;Slusallek,P. Light transport simulation with vertex connection and merging.ACM Transactions on GraphicsVol.31,No. 6,Article No.192,2012.

[21]Hachisuka,T.;Pantaleoni,J.;Jensen,H.W.A path space extension for robust light transport simulation.ACM Transactions on GraphicsVol.31,No.6,Article No.191,2012.

[22]Kajiya,J.T.Anisotropic re fl ection models.ACM SIGGRAPH Computer GraphicsVol.19,No.3,15-21 1985.

[23]Ward,G.J.Measuring and modeling anisotropic re fl ection.ACM SIGGRAPH Computer GraphicsVol. 26,No.2,265-272,1992.

[24]Ashikhmin,M.;Shirley,P.An anisotropic phong BRDF model.Journal of Graphics ToolsVol.5,No. 2,25-32,2000.

[25]Marschner,S.R.;Jensen,H.W.;Cammarano,M.; Worley,S.;Hanrahan,P.Light scattering from human hair fi bers.ACM Transactions on GraphicsVol.22, No.3,780-791.2003.

[26]Marschner,S.R.;Westin,S.H.;Arbree,A.;Moon,J. T.Measuring and modeling the appearance of fi nished wood.ACM Transactions on GraphicsVol.24,No.3, 727-734,2005.

[27]Irawan,P.;Marschner,S.Specular re fl ection from woven cloth.ACM Transactions on GraphicsVol.31, No.1,Article No.11,2012.

[28]Sadeghi,I.;Bisker,O.;De Deken,J.;Jensen,H.W. A practical microcylinder appearance model for cloth rendering.ACM Transactions on GraphicsVol.32,No. 2,Article No.14,2013.

[29]Cook,R.L.;Torrance,K.E.A re fl ectance model for computer graphics.ACM Transactions on GraphicsVol.1,No.1,7-24,1982.

[30]Torrance,K.E.;Sparrow,E.M.Theory for o ffspecular re fl ection from roughened surfaces.Journal of the Optical Society of AmericaVol.57,No.9,1105-1112,1967.

Fu-Jun Luan is a Ph.D.student in the Computer Science Department at Cornell University.He received his B.S. degreefrom Tsinghua University in 2015.He works on computer graphics with a focus on physically-based rendering and material appearance modeling.

Li-Fan Wu received his B.Eng.degree from Tsinghua University in 2015.He is a Ph.D.student of University of California,San Diego. His research interests include realistic rendering and image synthesis.

Kun Xu is an associate professor in the Department of Computer Science and Technology,Tsinghua University. Hereceived hisPh.D.degreefrom Tsinghua University in 2009. His research interests include realistic rendering and image/video editing.

Other papers from this open access journal are available at no cost from http://www.springer.com/journal/41095. To submit a manuscript,please go to https://www. editorialmanager.com/cvmj.

1TNList, Department of Computer Science and Technology, Tsinghua University, Beijing 100084, China.E-mail:xukun@tsinghua.edu.cn

Manuscript received:2014-11-21;accepted:2015-03-04