The Beam Radiance Estimate for Vo...
The Beam Radiance Estimate for Volumetric Photon Mapping Wojciech Jarosz, Matthias Zwicker, and Henrik Wann Jensen University of California, San Diego Abstract We present a new method for efficiently simulating the scattering of light within participating media. Using a theoretical reformulation of volumetric photon mapping, we develop a novel photon gathering technique for participating media. Traditional volumetric photon mapping samples the in-scattered radiance at numerous points along the length of a single ray by performing costly range queries within the photon map. Our technique replaces these multiple point-queries with a single beam-query, which explicitly gathers all photons along the length of an entire ray. These photons are used to estimate the accumulated in-scattered radiance arriving from a particular direction and need to be gathered only once per ray. Our method handles both fixed and adaptive kernels, is faster than regular volumetric photon mapping, and produces images with less noise. Keywords: participating media, light transport, global illumination, rendering, photon tracing, photon map, ray marching, nearest neighbor, variable kernel method. Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: raytracing color, shading, shadowing, and texture I.6.8 [Simulation and Modeling]: Monte Carlo G.1.9 [Numerical Analysis]: Fredholm equations integro-differential equations. 1. Introduction The appearance of many natural phenomena, such as human skin, clouds, fire, water, or the atmosphere, are strongly in- fluenced by the interaction of light with volumetric media. Therefore, efficiently rendering scenes with participating me- dia has been an area of interest within computer graphics. This problem is challenging, however, because accurately simulating light transport in participating media is compu- tationally very expensive. Cerezo et al. [CPP*05] provide a recent and comprehensive overview of the wealth of research that has been devoted to address this problem. Light transport in arbitrary participating media is modeled by the radiative transfer equation [Cha60]. The pioneering contributions by Kajiya and von Herzen [KH84] and Rush- meier and Torrance [RT87] are at the beginning of a long list of work on rendering participating media in the computer graphics community. Some of the most popular techniques to date are based on stochastic path-tracing and Monte Carlo integration [PM93, LW96, PKK00]. These approaches are attractive because of their sound underlying theoretical frame- work and their generality. They are unbiased and guaranteed to converge to the exact solution. In addition, it is straight- forward to include heterogeneous media, anisotropic phase functions, and scattering from surfaces. The downside of these approaches is that they suffer from noise that can only be overcome with a huge computational effort. One strategy to solve this issue is to make simplifying as- sumptions about the participating media. For example, homo- geneous media with a high scattering albedo can be modeled accurately using a diffusion approximation [Sta95,JMLH01], which leads to very efficient rendering algorithms. Premoze et al. [PARN04], under the assumption that the medium is tenuous and strongly forward scattering, use a path integral formulation to derive efficient rendering algorithms. Sun et al. [SRNN05] render single scattering in real time, but with- out shadowing effects. In contrast, photon mapping [JC98] improves the efficiency of path-tracing without making additional assumptions about the properties of the medium being rendered. Similar to Monte Carlo methods, photon mapping handles isotropic, anisotropic, homogeneous, and heterogeneous media of ar- bitrary albedo. A disadvantage of photon mapping is that it introduces bias to the solution of the radiative transfer equa- tion. In practice, however, this bias is preferable to the noisy solutions of pure Monte Carlo methods. Intuitively, photon mapping works by splitting the energy emitted by each light source into discrete packets, so called photons. In a first pass, the propagation of light is simulated
2 W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate Symbol Units Description ��s(x) m-1 Scattering coefficient at x ��a(x) m-1 Absorption coefficient at x ��t(x) m-1 Extinction coefficient at x p(x,~ �� ,~ �� 0) sr-1 Normalized phase function ��(x���x0) unitless Optical thickness: R x x0 ��t(x) dx Tr(x���x0) unitless Transmittance: e-��(x���x0) A m2 Surface area V m3 Medium volume ��� sr Hemisphere of directions ���4�� sr Sphere of directions L(x���x0) Wm-2sr-1 Outgoing radiance at x towards x0 L(x,~ �� ) Wm-2sr-1 Incident radiance at x from ~ �� We(x,~ �� ) m-3sr-1 Importance at x towards ~ �� We(x���x0) m-3sr-1 Importance at x towards x0 Table 1: Definitions of quantities used throughout the paper. by scattering photons in the scene. At each scattering event (at surfaces or within the medium) the incident energy carried by a photon is stored in a photon map. In a second pass, the photon map is used to evaluate the radiance at discrete points in the scene by locally computing the photon density. To com- pute the radiance carried along each ray towards the eye, the radiance is estimated within the medium at several sample points along the ray [JC98]. At each point the attenuation through the medium to the eye is computed, and the attenu- ated radiance is added to the ray. The main disadvantage of this procedure, however, is that it is difficult to find a good dis- tribution of sample points along the ray. On one hand, if not enough sample points are used, the result is likely to be noisy. On the other hand, increasing the number of sample points is very costly and can slow down rendering significantly. In this paper, we propose a novel approach for computing the contribution of in-scattered radiance. We gather photons along viewing rays and analytically compute their contri- butions, without point sampling. We present the following contributions: ��� We combine the theory from Veach [Vea98] and Pauly et al. [PKK00] to derive a reformulation of volumetric photon mapping as a solution to the measurement equation. This theory allows for arbitrary measurements of radiance to be computed within participating media, where a measure- ment is simply an integral of the radiance multiplied with a weighting function. ��� Using this new theory, we present an improved radiance estimate for volumetric photon mapping based on ���beam gathering.��� This technique eliminates the need for stepping through the medium to find photons. Instead, it gathers all photons along a ray. We show how to efficiently implement this new gathering technique for both fixed and adaptive smoothing kernels and demonstrate that our method pro- duces images with less noise than conventional photon mapping. The rest of this paper is organized as follows. In Section 2, we review the theory of radiance transport within partici- pating media and the volumetric photon mapping method. In Section 3, we reformulate volumetric photon mapping in terms of the measurement equation and show how the photon map can be used to estimate any measurement of radiance within the scene. In Section 4, we present our new beam radiance estimate using this theory and describe the data structures needed to evaluate it efficiently. Finally, we show comparisons of our approach to conventional photon mapping in Section 5 and discuss avenues of future work in Section 6. 2. Photon Mapping in Participating Media Light transport within participating media is described by the radiative transfer equation (RTE) [Cha60], which defines the radiance that reaches a point x from direction ~ �� as a sum of the exitant radiance from the nearest surface from this direction, and the accumulated in-scattered radiance from the medium between the surface and x (see Figure 1). This can be expressed as: L(x,~ �� ) = Tr(x���xs)L(xs,~ �� )+ Z s 0 Tr(x���xt )��s(xt)Li(xt,~ �� ) dt, (1) where Tr is the transmittance, s is the distance through the medium to the nearest surface at xs = x-s~ �� , and xt = x-t~ �� with t ��� (0,s). We define the remaining quantities in Table 1. The surface radiance, L(xs,~ �� ), at the boundary of the medium is governed by the rendering equation [Kaj86]. The in-scattered radiance, Li(xt ,~ �� ), depends on the radiance ar- riving at xt from all directions~ �� t over the sphere of directions ���4�� and is defined as: Li(xt ,~ �� ) = Z ���4�� p(xt ,~ �� ,~ �� t )L(xt,~ �� t ) d~ �� t , (2) where p is the normalized phase function. Volumetric photon mapping [JC98] solves the RTE using a combination of photon tracing, ray-marching, and density estimation. In a preprocess, packets of energy are shot from light sources, scattered at surfaces and within the medium, and their interactions are stored in a global data structure. Dur- ing rendering, ray marching is used to numerically integrate L(xs, ��) medium object ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� x xs xt Tr(x���xt) ��s(xt) Tr(x���xs) Li(xt, ��) Figure 1: The radiance reaching the eye L(x,~ �� ) is the sum of the radiance from the surface L(xs,~ �� ) and the accumulated in-scattered radiance Li(xt ,~ �� ) along a ray.
W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 3 Conventional Gathering Beam Gathering Gathered Other photons Missed Double-counted Figure 2: Conventional gathering (left) searches for photons in a sphere around numerous samples along the ray. Our method (right) assigns a radius to each photon and performs a single range search to find all photons along the length of the entire ray. Equation 1 for radiance seen directly by the observer, L(x,~ �� ) ��� Tr(x���xs)L(xs,~ �� )+ S-1 ��� t=0 Tr(x���xt )��s(xt)Li(xt,~ �� )���t ! , (3) where ���s is the length of each segment along the ray and x0,...,xs are the sample points for each segment (x0 is the point where the ray enters the medium and xs is a point on a surface past the medium). The most expensive part to compute in Equation 3 is the in-scattered radiance Li, because it involves accounting for all light arriving at each point xt along the ray from any other point in the scene. Instead of computing these values indepen- dently for each location, photon mapping gains efficiency by reusing the computation performed during the photon tracing stage. The in-scattered radiance is approximated using den- sity estimation by gathering photons within a small spherical neighborhood of radius r around each sample location xt, Li(xt ,~ �� ) ��� n ��� i=1 p(xt ,~ �� ,~ �� i)�����i 4 3 ��r3 , (4) where �����i is the power of photon i, and ~ �� i is its incident direction [JC98]. Though photon mapping is much more efficient than brute- force techniques like path tracing, the density estimation requires searching for photons within a global data structure, which is quite expensive. This formulation is suboptimal, firstly because it may gather the same photons more than once if the spherical neighborhoods overlap and, secondly, because it can lead to noise if the step size is too large and photons are omitted (see Figure 2). 3. Reformulation of Volumetric Photon Mapping Our technique solves these shortcomings by querying once for photons along the length of an entire ray, instead of mul- tiple times near points along the ray (see Figure 2). More formally, whereas regular photon mapping estimates Li at discrete points using Equation 4, our main contribution is to directly estimate Z s 0 Tr(x���xt )��s(xt)Li(xt,~ �� ) dt (5) along rays. Though the explanation of photon mapping from the pre- vious section is appealing at an intuitive level, it does not rigorously present the algorithm as a numerical solution to the RTE. Furthermore, this explanation is heavily tied to the geometric interpretation of gathering photons within a disc (on surfaces) or within a sphere (in participating media). In order to avoid these limitations and use the photon map to estimate general radiometric quantities in the volume, such as Equation 5, we use a more flexible derivation of parti- cle tracing methods presented by Veach [Vea98]. We extend this derivation to handle participating media by combining it with the generalized path integral formulation of the radiative transport equation [PKK00] (Section 3.1) and show how to represent particle tracing algorithms like volumetric photon mapping in terms of the measurement equation (Sections 3.2 and 3.3). Finally, we show how to use the same photon maps to estimate more general quantities of radiance (Section 3.4). 3.1. Generalized Path Integral Formulation We use the path integral formulation of the RTE, which arises by recursively expanding the right hand side of Equation 1. Instead of expressing the radiance equilibrium recursively, the resulting path integral formulation is a sum over light- carrying paths of different lengths. In order to do this, we define the path space, the corresponding differential measure, and generalized radiometric terms using notation inspired by Pauly et al. [PKK00]. A light path ��k x l is a set of k + 1 vertices xi. Each path is classified according to its path characteristic l ��� N, which determines for each vertex whether it is in the volume or on a surface. This allows us to integrate over different measures for scattering events at surfaces and within the volume. We define the path characteristic l of a path ��k x l such that the ith bit of l, bi(l), equals 1 if vertex xi is on a surface and bi(l) = 0 if it is in the volume. The space of all paths of length k with characteristic l is therefore: Xk l = ( ��k x l = x0,x1,...,xk xi ��� ( V, if bi(l) = 0 A, if bi(l) = 1 ) , (6) for 1 ��� k ��� and 0 ��� l 2k+1, and where V and A are the media volume and surface area respectively (see Figure 3 for an illustration of this notation). We define the corresponding differential measure at a path vertex xi as: d��l(xi) = ( dV(xi), if xi ��� V, i.e. when bi(l) = 0 dA(xi), if xi ��� A, i.e. when bi(l) = 1 . (7) Additionally, in order to express Equation 1 in terms of paths we need to transform the integration domain from solid
4 W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 1 0 1 1 sensor x0 x1 x2 x3 object medium light source fr (x���2) p(x1) ��� ��� ��� ������������ ��� ��� ��� ��� ������ ��� ������x0 ��� ��� ��� ��� ��������� ��� ��� ��� ��� ������������������������������G ������(3 ������������x ������������������������ ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��( G x 1 ��� ) �� G (x��� 2 ��� x��� 1 )��� �� x ��� 2 ) Figure 3: An example path, ��13, x 3 with length k = 3. The path characteristic, l = 1101b = 13, concatenates the type of scattering event from each path vertex. The radiance transported along the path ��(��13) L x 3 is the emitted radiance at the light multiplied by a series of scattering events (blue) and geometry terms (green). angle (��� and ���4��), to area or volume (A and V) depending on the type of scattering event. This transformation is achieved using the generalized geometry term: ��(x���y) G = V (x���y)Dx(y)Dy(x)��(x)Tr(x���y) kx- yk2 (8) where Dx(y) = ( 1, if x ��� V ~ n (x)��~ �� xy, if x ��� A . (9) We define ~ �� xy to be the unit direction vector from x to y, and Dy(x) is defined symmetrically to Dx(y) for both cases. The visibility function, V , is 1 if x and y are mutually visible, and 0 otherwise. The ��(x) function returns the scattering coeffi- cient ��s(x) if x ��� V and 1 otherwise. Note that Equation 8 simplifies to the regular geometry term if no participating media is present. Similarly, we generalize scattering events and emitted ra- diance. We define f�� to be the generalized scattering function combining the phase function and the surface BRDF f��(xi) = ( p(xi+1 ���xi ���xi-1), if xi ��� V fr(xi+1 ���xi ���xi-1), if xi ��� A , (10) and �� L e is the generalized emitted radiance �� L e(xi ���xi-1) = ( ��a(xi) ��s(xi) Le(xi ���xi-1), xi ��� V Le(xi ���xi-1), xi ��� A , (11) where the emitted radiance in a volume needs to be multiplied by the absorption, not the scattering, coefficient. Given this notation, the generalized path integral formu- lation expresses the outgoing radiance Lo at x1 towards x0 as a sum over all possible light paths arriving at x1. This includes light paths of all lengths k, as well as all possible characteristics l for each length Lo(x1 ���x0) = ��� ��� k=1 2k+1-1 ��� l=0 ��(��k). L x l (12) ��(��k) L x l measures the amount of radiance transported along a path ��k x l and is defined as ��(��k) L x l = Z ��l (xk) ������ Z ��l (x2) �� L e(xk ���xk-1) (13) k-1 ��� j=1 f��(x j) ��(x G j+1 ���x j) ! d��l(x2)������d��l(xk). The radiance transported along an example path is shown in Figure 3. 3.2. The Measurement Equation Many global illumination algorithms can be described in terms of the measurement equation. The measurement equa- tion describes an abstract measurement of incident radiance taken over some set of rays in a scene I = hWe,Li = Z V Z ���4�� We(x,~ �� )L(x,~ �� ) d~ �� dV(x). (14) The importance function We represents an abstract measuring sensor and is defined over the whole ray space V �� ���4�� (though typically We is non-zero for only a small subset of this domain). Path tracing, for instance, measures the contribution of ra- diance arriving over the area of a pixel. Radiosity algorithms integrate the contribution of radiance over basis functions de- fined on the scene geometry. Both of these approaches can be described using Equation 14 with an appropriate importance function. In his dissertation, Veach [Vea98] showed how particle trac- ing methods for surface illumination can also be expressed as a solution to the measurement equation by using the path in- tegral form of the rendering equation. We extend on this idea and use the generalized path integral formulation to describe volumetric photon tracing in the same way. 3.3. Volumetric Photon Tracing Photon tracing methods can be thought of as a way of gen- erating samples from the scene���s equilibrium radiance dis- tribution and then using this single collection of samples to render the entire image. The photon tracing stage generates N weighted sample rays, or photons, (��i,xi,~ �� i), where each (xi,~ �� i) is a ray and ��i is a corresponding weight. Our goal is to use these samples to take unbiased estimates of the radiance as a weighted sum, E " 1 N N ��� i=1 We(xi,~ �� i)��i # = hWe,Li, (15) for an arbitrary importance function We. We must therefore determine the proper distribution of samples for Equation 15 to hold. We start by manipulating the measurement equation on the right-hand side. To express the measurement equation using paths, we expand Equation 14 in terms of the outgoing
W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 5 radiance and convert to area form, introducing an additional geometry term, hWe,Li = Z ��l (x1) Z ��l (x0) We(x0 ���x1)Lo(x1 ���x0) (16) ��(x1 G ���x0) d��l(x0) d��l(x1). By combining with Equations 12 and 13 and moving the summations outside the integrals, this becomes ��� ��� k=1 2k+1-1 ��� l=0 Z ��l (xk) ������ Z ��l (x0) We(x0 ���x1)�� L e(xk ���xk-1) (17) k-1 ��� j=1 f��(x j) ��(x G j+1 ���x j) ! ��(x1 G ���x0)d��l(x0)������d��l(xk) # . The Monte Carlo estimator for Equation 17 using N samples is E " 1 N N ��� i=1 We(xi,0 ���xi,1)Ri # , (18) where each Ri is a random-walk path of length ki generated using Russian roulette Ri = �� L e(xi,ki ���xi,ki-1 ) pdf (xi,ki ,xi,ki-1) (19) ki-1 ��� j=1 1 qi, j f��(xi, j) ��(xi, G j+1 ���xi, j) pdf (xi, j-1) ! ��(xi,1 G ���xi,0), and qi, j was the probability of terminating Ri at the jth vertex. Comparing Equation 18 with 15 we see that in order to satisfy the requirements we need to set ��i = Ri. Connection to Conventional Photon Tracing. Though de- rived in a different fashion, Equation 19 is exactly how con- ventional photon mapping distributes photons within the scene. For instance, for a diffuse area light, photons are emit- ted using a cosine distribution with the power of the light source. In Equation 19 photons are emitted with the radiance of the light source divided by the pdf of choosing a position and direction on the light. These quantities are equivalent. Hence the particles generated above represent differential flux. The correspondence between the photon powers [JC98] and the sample weights is �����i = ��i N . 3.4. Radiance Estimation Using the Measurement Equation The main advantage of the reformulation in Section 3.3 is that it naturally accommodates computation of any measurement of radiance within the scene simply by using an appropriately defined importance function We. In this section, we first show how the conventional estimate for in-scattered radiance can be expressed as a measurement. We then go one step further and show how to derive a beam radiance estimate which approximates Equation 5 along rays directly. Conventional Radiance Estimate. The conventional radi- ance estimate approximates the value of the in-scattered ra- diance Li at fixed points within the scene. To express this using the theory from Section 3, we need to transform Equa- tion 2 into the measurement equation. Since the measurement equation is an integral over all of ray space (V �����4��), we artificially expand Li to also integrate over the volume Li(xt ,~ �� ) = Z p(xt ,~ �� ,~ �� t)L(xt ,~ �� t) d~ �� t (20) = Z���4Z�� V ���4�� ��(kx0 - xtk)p(x0,~ �� ,~ �� t) L(x0,~ �� t) d~ �� t dV(x0). In order to keep the expressions equivalent when we add the integration over volume, we also introduce a Dirac delta function ��. The bottom row of the above equation is now the measure- ment equation where We = ��(kx0 - xtk)p(x0,~ �� ,~ �� t ). Hence we can compute an unbiased estimate using the photon map by evaluating Equation 15 with this importance function. However, in order to obtain a useful estimate of radiance at all points in the scene a normalized kernel function is used in place of the delta function. This is where bias is introduced in the photon mapping method. Another interpretation is that by replacing the delta function with a kernel, photon mapping computes an unbiased estimate of blurred radiance. Jensen and Christensen [JC98] use a constant three-dimensional ker- nel with a radius based on the nth nearest neighbor. This results in the following radiance estimate by applying Equa- tion 15 Li(xt ,~ �� ) ��� Z V Z ���4�� Kn(kx0 - xtk)p(x0,~ �� ,~ �� t) (21) L(x0,~ �� t) d~ �� t dV(x0), ��� 1 N N ��� i=1 Kn(kxi - xtk)p(xi,~ �� ,~ �� i)��i (22) where the kernel Kn is defined as Kn(r) = ( 3 4��dn 3 if r ��� [0,dn] 0 otherwise , (23) and dn is the distance to the nth photon. Note that this is equivalent to the conventional volumetric radiance estimate in Equation 4. Beam Radiance Estimate. A similar procedure can be used to derive an estimate for the accumulated in-scattered radi- ance along an entire ray. To accomplish this, we first expand out Li in Equation 5 and then artificially inflate the resulting
6 W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate r x �� x t �� Figure 4: In the beam radiance estimate, x0 is parameterized in cylindrical coordinates, (t,��,r), about the ray (x,~ �� ). An unbiased estimate would only consider points directly on the ray, while a biased version uses a kernel (shown in grey) to blur the radiance within a beam. expression to integrate over the whole volume, Z sZ 0 ���4�� Tr(x���xt)��s(xt)p(xt ,~ �� ,~ �� t)L(xt ,~ �� t) d~ �� t dt = (24) Z R Z 2�� 0 Z R Z ���4�� ��(r)(H(t)- H(t - s))Tr(x���x0)��s(x0) p(x0,~ �� ,~ �� t)L(x0,~ �� t) d~ �� t dr d�� dt. (25) R is the set of real numbers and x0 is expressed in cylindrical coordinates, (t,��,r), about (x,~ �� ), where r is the radius to the ray (see Figure 4). We have added a Dirac delta function �� as before, and the Heaviside step functions (H(x) = 1 when x 0 and 0 otherwise) constrain the computation to t ��� (0,s). Equation 25 is now equivalent to the measurement equation where the integral over volume has been converted into cylindrical coordinates and where We = ��(r)(H(t) - H(t - s))Tr(x���x0)��s(x0)p(x0,~ �� ,~ �� t). Since the probability of photons landing exactly on the ray (x,~ �� ) is zero, we introduce bias by blurring the radiance and replacing the delta and step functions with a smooth kernel, K. This integral can then be estimated with the measurement equation using the photons as: Z R Z 2�� 0 Z R Z ���4�� K(t,��,r)Tr(x���x0)��s(x0)p(x0,~ �� ,~ �� t) L(x0,~ �� t) d~ �� t dr d�� dt = (26) 1 N N ��� i=1 K(ti,��i,ri)Tr(x���xi)��s(xi)p(xi,~ �� ,~ �� i)��i, (27) where (ti,��i,ri) = xi are the cylindrical coordinates of photon i about the ray. The blurring in the conventional radiance estimate is spher- ical and so the kernel needs to be normalized for 3D. However, with the beam radiance estimate, we blur in two dimensions (perpendicular to the ray) since the radiance we are com- puting already includes the integration along the ray itself. Therefore, the kernel in the beam estimate is normalized for 2D. 3.5. Kernel Radiance Estimation For both the conventional and beam radiance estimates the characteristics of the bias and blur are determined by the smoothing function chosen. Several options exist for applying a smoothing kernel to the photon map data. The kernel method uses a fixed-radius smoothing kernel and results in a uniform blur of radiance within the scene. In practice, using a fixed-width circular kernel implies that in order to evaluate the beam radiance estimate (Equation 26) using the photon map (Equation 15) we only need to con- sider photons which are located within a fixed-radius cylinder about the ray (x,~ �� ). Alternatively, the equivalent dual inter- pretation considers each photon as the center of an oriented disc facing the ray and all photon-discs that intersect the ray need to be found. If the density of photons varies significantly it can be diffi- cult to choose a single radius that works well for all regions of the scene. This can be solved by allowing the size and shape of the blurring kernel to vary spatially. In conventional photon mapping, the nearest neighbor method (k-NN) is used to adapt the kernel width to the local density. Generaliz- ing point-based k-NN to a visually comparable range search along rays is challenging. However, spatial variation can eas- ily be applied to the dual photon-discs interpretation using the variable kernel method (VK) [BMP77]. A smoothing kernel is ���attached��� to each photon and the radius of the ker- nel is allowed to vary from one photon to another, based on the local density. In contrast to k-NN estimation, where the kernel widths vary based on the distance from the evalua- tion location to the data points, in the VK method the kernel widths only depend on the data points themselves. In order to facilitate this, the method relies on a pilot estimate of the local density at each data point in order to assign the kernel widths. This is the approach we take. 4. Algorithm In order to use the dual interpretation to evaluate the beam radiance estimate (Equation 26), we need an efficient way of locating all photon-discs that intersect an arbitrary ray. Additionally, to use variable width kernels we need to effi- ciently compute a radius for each photon in the photon map. At a high level, our volumetric photon mapping technique performs the following five steps: 1. Shoot photons from light sources. 2. Construct a balanced kd-tree for the photons. 3. Assign a radius for each photon. 4. Construct a bounding-box hierarchy over the photon-discs. 5. Use the photon BBH to render the image. Steps 1 and 2 are identical to conventional photon mapping while 3���5 are unique to our approach and are explained in more detail below. Photon Radius Computation. We augment the traditional photon map by associating a radius with each photon. For fixed width kernels the radius is a constant for all photons
W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 7 d e a b c g f Construct balanced kd-tree Compute radius per photon Construct BBH over photon-discs d e a b c g f d e a b c g f Figure 5: After photons have been distributed in the scene, our algorithm constructs a balanced kd-tree (left). We assign a valid radius to each photon by querying in the kd-tree (middle). Finally, we rapidly construct a bounding-box hierarchy over the photon-discs (right) by reusing the same hierarchical structure (shown in red) as the kd-tree. Algorithm 1 CONSTRUCTBBH(p) Require: p is a node in a balanced photon map. Ensure: The subtree at p contains a valid BBH. 1: B ��� BOUNDINGBOX(p.position, p.radius) 2: if p.leftChild then 3: B ��� B ��� CONSTRUCTBBH(p.leftChild) 4: end if 5: if p.rightChild then 6: B ��� B ��� CONSTRUCTBBH(p.rightChild) 7: end if 8: p.bbox ��� B 9: return B and does not need to be explicitly stored. For variable width kernels using the VK method, we perform a density estima- tion at each photon to assign a radius. At each photon we compute the local density by estimating the distance to the nth nearest photon and use this as the photon-disc���s radius. This pilot estimate is performed using the photon map kd-tree. The value n plays the same role as in the conventional radiance estimate, it controls the amount of blur. As an optimization, we only search for the nearest m n photons and estimate the necessary radius for n photons. By assuming locally uniform photon density, if di,m is the dis- tance to the mth photon from photon i, we estimate the radius as ri = di,m 3 p n m . The m parameter controls the sensitivity of the computed radius to the local variation in density. Very small values of m, m 5, can produce noisy radii, which change drastically between neighboring photons, while large values are more expensive to compute. In practice, we have found that m = ��� n works well as a default value and this value was used for all our scenes, significantly accelerating the preprocessing step. Bounding Box Hierarchy Construction. In order to effi- ciently locate all photons-discs which intersect a ray, we construct a bounding box hierarchy. Heuristics for construct- ing efficient BBHs have been extensively studied within the context of ray tracing [WMG*07]. However, the performance characteristics of our ray intersections are different than for regular ray tracing since we are interested in all intersections, not just the first intersection along a ray. Furthermore, the best heuristics tend to induce long construction times. We employ a rapid construction scheme by exploiting the information in the photon map kd-tree and re-using that hierarchy for our BBH. For each photon in the photon map, we compute the bound- ing box of all descendent photon-discs. The bounding box of each node encloses the node���s photon radius and the bound- ing boxes of its two child nodes. The computation starts at the leaves and propagates upwards through the tree. This procedure results in a balanced median-split-style BBH, but unlike traditional BBHs our hierarchy contains photons at interior nodes, not just at the leaves. Figure 5 illustrates the relationship between the kd-tree and the BBH. The BBH can be constructed by passing the root of the photon map kd-tree to Algorithm 4. Given a balanced kd-tree, this linear-time construction procedure is extremely fast and produces BBHs which can be efficiently traversed for nearby photons during rendering. After the BBH is constructed the photon map kd-tree is no longer used and can be freed from memory. Using a BBH and a per-photon radius, an additional 7 floating-point values need to be stored, increasing the size of each photon from 20 bytes to 46 bytes. The Beam Radiance Estimate. During rendering we es- timate the accumulated in-scattered radiance (Equation 5) along viewing rays by locating all photons whose bounding spheres intersect the ray. These photons are found using a standard ray-BBH intersection traversal. The contribution from each photon (��i,xi,~ �� i) is accumulated using Equa- tion 26 however, with the variable kernel method, a kernel Ki is defined per photon. This leads to the following radiance estimate 1 N N ��� i=1 Ki(x,~ �� ,s,xi,ri)Tr(x���xi0)��s(xi0)p(xi,~ �� ,~ �� i)��i, (28) where xi0 = x +ti~ �� is the projection of the photon location xi onto the ray���s direction ~ �� , and ti = (xi - x) ��~ �� . We define
8 W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate Scene Cornell Stage Cars Lighthouse Photon Tracing N Shoot (s) Balance (s) Radius (s) 0.4M 1.50 0.30 2.0 1M 3.25 0.76 5.0 2M 19.0 1.50 2.0 1M 2.83 0.78 6.0 Fixed Radius Comparison r ���t C (m) O (m) 0.4 0.40 3:19 3:00 0.3 0.20 4:21 4:15 0.4 1.25 1:30 1:30 0.4 0.25 1:07 0:59 Adaptive Radius Comparison r+ n ���t C (m) O (m) 0.6 1.5K 0.80 4:03 3:35 0.5 0.5K 0.70 6:38 6:22 0.5 1K 1.25 2:02 1:53 0.5 0.4K 1.00 1:12 1:05 Table 2: Rendering parameters and timings (in seconds (s) and minutes (m)) for all example scenes. Statistics relating to the photon tracing preprocess are shown in the first table. The middle and right tables compare our method (O) to conventional photon mapping (C) using a fixed width kernel and adaptive width kernel. The r column represents the fixed-width kernel radius, while r+ is the maximum search radius and the number of nearest neighbors is n. Scene ��s ��a g Cornell 0.225 0.225 0.00 Stage 0.225 0.225 0.75 Cars 0.06 0.015 0.00 Lighthouse 0.24 0.010 0.75 Table 3: Scattering properties of the participating media used in each of the example scenes. the kernel as Ki(x,~ �� ,s,xi,ri) = ( ri-2K2 di ri if di ��� [0,ri] 0 otherwise , (29) where ri is the pre-computed radius for photon i. We use Sil- verman���s two-dimensional biweight kernel [Sil86] along the ray, K2(x) = 3��-1(1- x2)2, where di is the shortest distance from photon i to the ray. We chose this kernel because it is smooth, efficient to evaluate, and has local support. Equa- tion 28 is the beam radiance estimate, and it replaces the ray marching computation from conventional photon mapping (second row of Equation 3). Heterogeneous Media. For homogeneous media, the trans- mission terms, Tr(x ��� xi0), can be computed explicitly for each photon during gathering. Beam gathering in heteroge- neous media can also be handled efficiently by first marching along the ray and saving the transmission values in a 1D lookup table. Then, during gathering, each photon���s trans- mission is computed by interpolating within the lookup table. This preprocess needs to be performed independently for each ray, just prior to gathering, so the lookup table can be reused. Furthermore, if single scattering is simulated sepa- rately by directly sampling light sources the lookup table can be populated during this marching step. 5. Results We compared our beam gathering technique against conven- tional volumetric photon mapping using ray marching. In order to isolate just the performance of the photon gathering methods, we use the photon map for both single and multiple scattering. We compared results on four test scenes: Cars, Lighthouse, Stage, and a Cornell box. For each scene we compare using a fixed gathering radius for both types of esti- mates, and we also compare the conventional estimate using k-NN to the beam estimate using the VK method. The images were all rendered with a maximum dimension of 1024 pixels with up to 4 samples per pixels on an Intel Core 2 Duo 2.4 GHz machine using one core. In our experimental setup, we first choose suitable gather- ing parameters (search radius and number of nearest neigh- bors n) and render the scenes using our method. We then use the same parameters using conventional photon mapping but tune the minimum step-size ���t to obtain approximately equal render times. Note that ���t is the minimum step-size and that exponential stepping is used to sample the ray according to transmission. Finally, we render a high quality result with conventional photon mapping using a very small step size as a ���reference.��� We show visual comparisons of the methods in Figures 6���7. All images of each scene are rendered using the same photon map. The only differences in quality and performance are due to the gathering method used. We report the render times and gathering parameters, as well as timings for constructing the photon maps in Table 2. We used the Henyey-Greenstein phase function for all scenes with the medium parameters specified in Table 3. In all cases, our method produces significantly higher qual- ity images than conventional photon mapping. This is because querying once for all photons along a ray is more efficient than repeatedly querying for photons near numerous samples along the ray. Not surprisingly, we see that the reduced blur- ring of the adaptive kernel gathering methods is essential for scenes like the Stage and Lighthouse, where concentrated beams of light are visible. However, at the same render time this advantage is difficult to discern in the conventional pho- ton mapping images. Though the k-NN and VK methods both adapt the width of the kernel to the local photon density, they are distinct approaches which result in similar, but not identical, density estimates. This is what accounts for the small differences in blurring between our adaptive results and the k-NN ���refer- ence��� images. However, as our results show, the same value of n produces visually comparable renderings using the two methods.
W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 9 Approx. Equal Time k-NN (6:38) Fixed (4:21) Conventional Radiance Estimate Reference Solution k-NN (5:18:33) Fixed (2:38:57) Beam Radiance Estimate Adaptive (6:22) Fixed (4:15) Figure 6: A comparison between the convention radiance estimate and our beam radiance estimate on the Stage scene with render times provided as (hours:minutes:seconds). Our method (right) produces images with much less noise than an equal time rendering using conventional volumetric photon mapping (middle) for both a fixed radius and an adaptive radius gathering approach. Our method does not require stepping but matches the quality of conventional photon mapping if a very small step size is used (left). 6. Discussion and Future Work Specific trade-offs between the k-NN and VK methods have been extensively studied within the density estimation liter- ature [Sil86]. For computer graphics, a visual advantage of the VK method is that the estimated function inherits all the differential properties of the smoothing kernel. In contrast, even with smooth kernels, the k-NN method results in esti- mates with a discontinuous first derivative. The VK method does, unfortunately, require a pre-process to compute the ker- nel width for each photon. However, the amount of time to compute adaptive radii for each photon using our method is relatively inexpensive (typically less than 1���2% of the total render time). On the other hand, k-NN gathering involves maintaining a priority queue and is much more costly than a fixed radius query. With our method, the gathering procedure is identical whether an adaptive or fixed radius kernel is used and no priority queue is needed. Though adaptive kernel widths can be a distinct advantage, in uniformly illuminated scenes a fixed radius may be suffi- cient. While assigning the radii is inexpensive, each resulting BBH node consumes more memory than a kd-tree node due to the additional 7 floating-point values needed to store the bounding box and radius. If memory is limited and a fixed ra- dius search is acceptable, then the regular photon map kd-tree can be used to perform beam gathering by tracing a cylinder through the kd-tree. This approach could be implemented in the spirit of ray-bundle traversals [WMG*07] by bounding the cylinder using 6 planes. This offers the additional bene- fit of re-using the exact same data structure as conventional photon mapping. Our photon-disc approach using a BBH is not inherently tied to participating media. VK density estimation could also be applied to photon mapping at surfaces. It would be inter- esting to see whether a similar technique could be beneficial for surfaces by querying the BBH for all photons that overlap with a surface location. The pilot estimation preprocess may also be beneficial at detecting and reducing boundary bias in photon mapping for surfaces. The advantages of beam gathering are more pronounced in large open environments where ray marching needs to be per- formed over large distances and where potentially interesting lighting is visible far away. One advantage of the conventional radiance estimate, however, is that a fast preview render can be obtained by using a large step-size. Since no step-size is used for beam gathering, all photons intersecting a ray need to be considered. In fact, these two approaches estimate the lighting integral using different pdfs. The ray marching computation employs exponential stepping, which distributes photon queries with a pdf proportional to the transmission term Tr. In contrast, our approach visits every photon that intersects a line and so concentrates effort where the lighting is important. Since the radiance seen by the eye is a product of the lighting and the transmission, optimally we should concentrate effort according to this product. One way of ex- ploiting this could be to prune the BBH ray traversal using Russian roulette based on some upper-bound on the transmis- sion term and the number of photons in each subtree. Such a scheme could further reduce render times by not considering photons which are too distant and faint to contribute much to the image. We leave this optimization as future work. It would be interesting to explore what other useful radiance estimate could be formed using the measure-
10 W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate k-NN (14:00) Fixed (11:00) k-NN (2:02) Fixed (1:30) Fixed k-NN/Adaptive Adaptive (1:53) Fixed (1:30) Approx. Equal Time Conventional Radiance Estimate Reference Solution Beam Radiance Estimate k-NN (52:45) Fixed (18:37) k-NN (1:12) Fixed (1:07) Fixed k-NN/Adaptive Adaptive (1:05) Fixed (0:59) k-NN (16:23) Fixed (14:43) k-NN (4:03) Fixed (3:19) Adaptive (3:35) Fixed (3:00) Figure 7: The Cornell box, Cars, and Lighthouse scenes. Render times are shown as (minutes:seconds). For both the fixed and adaptive gathering approaches our method produces noise-free results while conventional photon mapping suffers from significant noise, especially around distant light sources.
W. Jarosz & M. Zwicker & H. W. Jensen / The Beam Radiance Estimate 11 ment equation formulation. For instance, Cammarano and Jensen [CJ02] considered the problem of estimating the den- sity of photons within four dimensional ���cylinders��� to prop- erly simulate motion blurred caustics. This process could easily be formulated as a measurement by defining the impor- tance function over space-time. Photon splatting for participating media was presented by Boudet et al. [BPPP05], where the conventional photon mapping method needed to be meticulously transformed into a splatting algorithm. An extra benefit of our reformulation is that Equation 28 can immediately be seen as a splatting operation and could therefore be efficiently implemented in graphics hardware. Furthermore, whereas Boudet et al. only considered fixed-size splats, our VK approach could easily be used to adapt the splat size to the local photon density. A splatting approach could also be used in software to accelerate the computation for eye rays however, the more general beam gathering would need to be used for any secondary rays. 7. Conclusion In this paper, we showed how volumetric photon mapping can be expressed as a solution to the measurement equation. This formulation showed that any measurement of radiance within participating can be estimated using the photon map. We applied this formulation by using the photon map to di- rectly estimate accumulated in-scattered radiance along rays. This approach was implemented using an efficient beam gath- ering method which can be used for both fixed and adaptive width kernels. The resulting algorithm produces images with significantly less noise than conventional volumetric photon mapping using the same render time. References [BMP77] BREIMAN L., MEISEL W., PURCELL E.: Variable ker- nel estimates of multivariate densities. Technometrics 19, 2 (May 1977), 135���144. [BPPP05] BOUDET A., PITOT P., PRATMARTY D., PAULIN M.: Photon splatting for participating media. In GRAPHITE ���05: Proceedings of the 3rd international conference on Computer graphics and interactive techniques in Australasia and South East Asia (New York, NY, USA, 2005), ACM Press, pp. 197���204. [Cha60] CHANDRASEKAR S.: Radiative Transfer. Dover Publica- tions, New York, 1960. [CJ02] CAMMARANO M., JENSEN H. W.: Time dependent photon mapping. In EGRW ���02: Proceedings of the 13th Eurographics workshop on Rendering (Aire-la-Ville, Switzerland, Switzerland, 2002), Eurographics Association, pp. 135���144. [CPP*05] CEREZO E., P��REZ F., PUEYO X., SERON F. J., COIS X. SILLION F.: A survey on participating media rendering tech- niques. The Visual Computer 21, 5 (June 2005), 303���328. [JC98] JENSEN H. W., CHRISTENSEN P. H.: Efficient simulation of light transport in scences with participating media using photon maps. In SIGGRAPH ���98: Proceedings of the 25th annual con- ference on Computer graphics and interactive techniques (New York, NY, USA, 1998), ACM Press, pp. 311���320. [JMLH01] JENSEN H. W., MARSCHNER S. R., LEVOY M., HAN- RAHAN P.: A practical model for subsurface light transport. In SIGGRAPH ���01: Proceedings of the 28th annual conference on Computer graphics and interactive techniques (New York, NY, USA, 2001), ACM Press, pp. 511���518. [Kaj86] KAJIYA J. T.: The rendering equation. In SIGGRAPH ���86: Proceedings of the 13th annual conference on Computer graphics and interactive techniques (New York, NY, USA, 1986), ACM Press, pp. 143���150. [KH84] KAJIYA J. T., HERZEN B. P. V.: Ray tracing volume densities. In SIGGRAPH ���84: Proceedings of the 11th annual con- ference on Computer graphics and interactive techniques (New York, NY, USA, 1984), ACM Press, pp. 165���174. [LW96] LAFORTUNE E. P., WILLEMS Y. D.: Rendering Partici- pating Media with Bidirectional Path Tracing. In Proceedings of the eurographics workshop on Rendering techniques ���96 (London, UK, 1996), Springer-Verlag, pp. 91���100. [PARN04] PREMOZE S., ASHIKHMIN M., RAMAMOORTHI R., NAYAR S. K.: Practical rendering of multiple scattering effects in participating media. In Rendering Techniques (2004), pp. 363��� 373. [PKK00] PAULY M., KOLLIG T., KELLER A.: Metropolis light transport for participating media. In Proceedings of the Euro- graphics Workshop on Rendering Techniques 2000 (London, UK, 2000), Springer-Verlag, pp. 11���22. [PM93] PATTANAIK S. N., MUDUR S. P.: Computation of global illumination in a participating medium by Monte Carlo simulation. The Journal of Visualization and Computer Animation 4, 3 (July��� Sept. 1993), 133���152. [RT87] RUSHMEIER H. E., TORRANCE K. E.: The zonal method for calculating light intensities in the presence of a participating medium. Computer Graphics (SIGGRAPH ���87 Proceedings) 21, 4 (July 1987), 293���302. [Sil86] SILVERMAN B.: Density Estimation for Statistics and Data Analysis. Chapman and Hall, New York, NY, 1986. [SRNN05] SUN B., RAMAMOORTHI R., NARASIMHAN S. G., NAYAR S. K.: A practical analytic single scattering model for real-time rendering. ACM Trans. Graph. 24, 3 (2005), 1040���1049. [Sta95] STAM J.: Multiple Scattering as a Diffusion Process. In Rendering Techniques ���95 (Proceedings of the Sixth Eurographics Workshop on Rendering) (New York, NY, 1995), Hanrahan P. M., Purgathofer W., (Eds.), Springer-Verlag, pp. 41���50. [Vea98] VEACH E.: Robust Monte Carlo Methods for Light Trans- port Simulation. PhD thesis, 1998. [WMG*07] WALD I., MARK W. R., G��NTHER J., BOULOS S., IZE T., HUNT W., PARKER S. G., SHIRLEY P.: State of the art in ray tracing animated scenes. In STAR Proceedings of Eurograph- ics 2007 (Prague, Czech Republic, September 2007), Eurograph- ics Association, pp. 0���0.
The Beam Radiance Estimate for Volumetric Photon Mapping Wojciech Jarosz Matthias Zwicker Henrik Wann Jensen University of California, San Diego
http://www.kevinyank.com Motivation http://mev.fopf.mipt.ru 2 Wojciech Jarosz * In this talk, we are interested in rendering scene with participating media, or scenes where the volume or medium participates in the lighting interactions. * These are just a few example photographs of the types of effects that are caused by participating media.
medium object light source Theoretical Background 3 camera (eye)
object x Volume Rendering Equation 4 * The radiance, L, arriving at the eye along a ray can be expressed using the volume rendering equation. * at a high-level the meaning is pretty simple. * the radiance arriving at the eye is the sum of two terms: