9512.net
甜梦文库
当前位置:首页 >> 林学 >>

Eurographics Symposium on Rendering 2003, pp. 1–12 Per Christensen and Daniel Cohen-Or (Ed


Eurographics Symposium on Rendering 2003, pp. 1–12 Per Christensen and Daniel Cohen-Or (Editors)

Path Integration for Light Transport in Volumes
Simon Premo?e,1 Michael Ashikhmin2 and Peter Shirley1
2 1 Computer Science Department, University of Utah Computer Science Department, SUNY at Stony Brook

Abstract Simulating the transport of light in volumes such as clouds or objects with subsurface scattering is computationally expensive. We describe an approximation to such transport using path integration. Unlike the more commonly used diffusion approximation, the path integration approach does not explicitly rely on the assumption that the material within the volume is dense. Instead, it assumes the phase function of the volume material is strongly forward scattering and uniform throughout the medium, an assumption that is often the case in nature. We show that this approach is useful for simulating subsurface scattering and scattering in clouds.

1. Introduction The appearance of many materials (e.g., skin, fruits, snow, clouds) cannot be described by a simple BRDF-style re?ectance model. The main reason for this is volumetric scattering which manifests itself in important lighting effects for materials4, 5, 19 and scenes36 . The radiative transfer equation and propagation of light in a scattering medium have been both analytically and numerically studied in astrophysics, atmospheric optics, and more recently medical applications. For most problems with non-trivial boundary conditions, phase functions and initial conditions there are no analytic solutions. Most solutions are based on the diffusion approximation which assumes that enough scattering events have occurred for light to be uniformly scattered in all directions. This approximation has proven useful for generating images with subsurface scattering18 . Monte Carlo methods are also often used to compute radiative transport within a medium. Although simple and powerful, these methods suffer from slow convergence. Finite element methods are also used, but they require large amounts of storage to capture discontinuities and strong directional light distributions. Because the diffusion approximation is only appropriate for dense uniform media19 , there is a gap in the computer graphics literature when accurate approximations are desired for sparse or non-uniform media. This paper attempts to ?ll that gap using an alternative to the diffusion approximation based on Feynman’s path integral approach to solving quantum mechanics problems10 . Path integral formulations of physical processes have been used in physics to solve a
c The Eurographics Association 2003.

wide variety of problems including energy propagation in random media38 and transfer equation39 . The radiative transfer equation describing light propagation can be viewed as a collection of paths taken by radiation as it travels through space. A path integral is an integral over all such possible paths traveled by a photon. Radiative transfer is decomposed into a series of smaller problems formulated by the Green function propagator. Scattered optical ?elds are described using the concept of an ensemble of effective optical paths of partial contributions. This physical picture for treating light transport in multiple scattering media as a collection of most probable paths provides insight into the light propagation in a medium. Unlike the randomized approach to using paths of Metropolis Light Transport45 , path integral methods analytically ?nd the most important paths and develop analytical estimates based on them. Once the most important path is found, the multiple scattering contributions are only computed along the most probable paths and the rest of the paths are dealt with implicitly via analytic integration of multiple scattering using wellknown approximations. We therefore avoid computationally expensive direct numerical simulation of multiple scattering in the medium. We provide some intuition behind path integrals and demonstrate that some useful results of the theory can be obtained without any heavy mathematical tools. We discuss solutions to light propagation as a path integral (a formal sum) and how they can be used for rendering arbitrary scattering materials and media. Because many common phenomena cannot be described using only single scattering, we discuss some observable consequences of multiple scat-

Premo?e et al. / Path Integration for Light Transport in Volumes

tering and how they can be exploited for rendering applications. We show results for inhomogeneous media that cannot be achieved using the diffusion approximation. The restrictions on our method are different; while sparse and inhomogeneous media are allowed, we assume that the density and therefore scattering coef?cients vary smoothly, the phase function is arbitrary but constant within the medium. This expands the class of problems that can be attacked ef?ciently. We note that this paper serves mainly as an introduction to path integral methods and many improvements in ef?ciency have not yet been explored in our implementation. 2. Background and Previous Work Light Transport Approximations There has been much work in approximating radiative transfer in arbitrary media in many ?elds including computer graphics. Perez et al.30 survey and classify global illumination algorithms in participating media in detail. Pharr and Hanrahan31 and Premo?e35 also provide an extensive list of existing methods and background. Here we brie?y review recent methods proven practical and robust. Jensen and Christensen16 presented a two pass photon density estimation method. This method is simple, robust and ef?cient but suffers from large memory requirements to store photons if the extent of the scene is large or the lighting con?guration is very complex. The method is very practical and it has been used for many phenomena including smoke9 , ?re26 , stone8 , and wet materials17 proving its generality. On the other hand, the photon map becomes rather inef?cient in highly scattering media18 . Veach and Guibas45 presented a global illumination algorithm that found important paths, and then explored the path space locally because it was likely that other important paths would be nearby. Pauly et al.27 extended the method for participating media and proposed suitable mutation strategies for paths. Although extremely general and robust, as it could handle any lighting condition and con?guration, it still suffers from the classical Monte Carlo problems of noise and slow convergence. Stam37 presented a solution to multiple scattering in nonhomogeneous materials by solving the diffusion equation using a multigrid method. Jensen et al.18 introduced an analytical solution to the diffusion approximation to multiple scattering, which is especially applicable for translucent materials that exhibit considerable subsurface light transport. Their method relies on the assumption that the multiplyscattered light is nearly isotropic and cannot be easily extended to inhomogeneous materials. Lensch et al.22 implemented this method in graphics hardware and Jensen and Buhler15 extended this diffusion approximation to be computationally more ef?cient by precomputing and storing illu-

mination in a hierarchical grid. Narasimhan and Nayar25 described a physicall-based multiple scattering model for simulating weather effects such as fog, haze, mist and rain. An alternative description of light propagation was done by Pharr and Hanrahan31 who described a mathematical framework for solving the scattering equation in the context of a variety of rendering problems and a numerical solution in terms of Monte Carlo sampling. The scattering equation describes all scattering events inside the object and it does not depend on the incoming illumination. Unfortunately, there has been no other work exploiting this interesting paradigm of scattering objects and interactions between objects on larger scales. Lafortune20 also described the global re?ectance distribution function (GRDF) which corresponds to the scattering equation idea. Path Integral Methods Path integral techniques and functional integration have been widely used in statistical and quantum mechanics to solve propagators for Fokker-Planck and Schr?dinger equations. The concept of photon paths has been well-known in the theory of energy propagation in random media38 . Tessendorf39, 40 used the path integral approach to study the propagation of light in weakly scattering media such as water. Perelman et al.28, 29 described energy transport in a turbid medium using a quasi-particle Lagrangian, from which the most probable paths could be found. Wilson and Wang46 constructed a Lagrangian through a turbid medium using local path descriptors. Miller24 also constructed a stochastic Lagrangian path integral representation for Green’s evolution operator. Constantinou and Demetrescu7 showed the equivalence of the path integral formulation and virtual rays. Gross12 studied multiple scattering of a wave in a system of random and uncorrelated scattering particles. The path integral methods have also been studied in optical tomography and medical imaging. Near-diffusive scattering regimes are very important for obtaining diagnostic information about multilayer biological tissues where standard diffusion approximation fails32, 33 . Jacques and Wang14 presented a basic introduction to the path integral description of photon transport and discussed classical path for describing the most probable path of a photon. 3. Mathematical Preliminaries
MEDIUM PROPERTIES AND LIGHT PATHS . In an arbitrary medium, the underlying optical properties depend on bulk material properties such as density ρ(x), temperature T (x), and the particle absorption and scattering cross-sections, σa and σs . Optical properties of the medium are then described in terms of the scattering coef?cient b(x) = σs ρ(x), the absorption coef?cient a(x) = σa ρ(x), the extinction coef?cient c(x) = a(x) + b(x), and the phase function P(x, ω, ω ). The phase function P describes the probability density of light
c The Eurographics Association 2003.

Premo?e et al. / Path Integration for Light Transport in Volumes
incident light x’ temporal spreading
← ←

ω’ path γP


s’ β(s’)

angular spreading spatial spreading x ω viewing direction

?

Figure 1: Scattering in a highly scattering medium. Original radiance undergoes a series of scattering events that result in angular, spatial and temporal spreading of the original radiance distribution. Generic location in 3 Generic direction Absorption coef?cient at a point Scattering coef?cient at a point Extinction coef?cient at a point Mean cosine of the scattering angle Volume source distribution Phase function Propagator (Green function) Light point path Light propagation direction on the path Displacement along the path Distance along the path (arclength) Arclength of the path γ Spatial variability Mean square scattering angle
?

x ω a(x) b(x) c(x) g Q P(x, ω, ω ) G γP (s) β(s) d(s) s S θ2

Figure 2: Transfer geometry in an enclosed scattering medium. A photon originating at point x with direction ω travels along curved path γP of length s until it reaches the ?nal point x with direction ω. β(s ) is the direction of propagation at arclength parameter s on the curved path.

rameterized path γP (s) until it reaches the ?nal point x with direction ω, (see Figure 2). This path results from an accumulated random walk of propagation directions governed by scattering and absorption events along the path. Because the path γP is in general curved, its total length S is greater than the distance between starting and ending points: |x ? x| ≤ S. The direction of propagation along the path is de?ned by d γP (s) β(s) = ds . The path therefore satis?es the “two-sided” boundary conditions:
γP (0) = x γP (S) = x β(0) = ω β(S) = ω. (1)

Table 1: Notations for important terms used in the paper.

The displacement relative to the point x at distance s is obtained by integrating β:
?

d(s) = coming from incident direction ω scattering into direction ω upon scattering event at point x. The phase function is normalized so that 4π P(ω, ω )d ω = 1 and in most practical settings only depends on the phase angle cos θ = ω · ω . The mean cosine g of the scattering angle is de?ned:
? ?

s 0

β(s )ds .

g=



P(ω, ω )(ω · ω )d ω .

These optical properties are inherent, because they depend solely on the medium and not on the structure of the incoming light ?eld. Upon entering the medium, incoming light undergoes a series of scattering and absorption events that modify both the directional structure of the incoming light ?eld and its intensity. As a result of multiple scattering events, the original radiance distribution undergoes angular, spatial and temporal spreading which result in different radiance distribution. Figure 1 shows spreading effects in an arbitrary highly scattering medium. Table 1 summarizes important terms and quantities used in the paper. Consider a photon that originates at point x with direction ω and traveling in a medium along a curved arclength pac The Eurographics Association 2003.

Light undergoes a series of scattering and absorption events along the path. Note that if we ignore exact backscattering which returns photons back into the same path, the intensity of “original” light will be only diminishing because of these events since in-scattering will be due to photons traveling a different path in the medium. Therefore, if we introduce the effective attenuation τ which determines how much will the light intensity be reduced along the length of the path, the radiance L in the medium will be proportional to
L?
all paths



e?τ( path) .

The path integral formulation of light transport is essentially a mathematically rigorous expression of this simple idea.
RADIATIVE TRANSFER . Light transport in arbitrary media is described by the radiative transport equation2, 13 :

(ω · ? + c(x ))L(x, ω) = b(x)


P(ω, ω )L(x, ω )d ω + Q(x, ω),

(2)

where Q(x, ω) is the source term. In computer graphics, the source terms Q(x, ω) is often due to light emitted by the

Premo?e et al. / Path Integration for Light Transport in Volumes

medium itself (Le (x, ω)). It is often convenient to split the total radiance within the medium into components and write it as the sum of unscattered (direct) radiance Lun , the emission Le and the scattered radiance Lsc : L(x, ω) = Lun (x, ω) + Lsc (x, ω) + Le (x, ω). Here Lun is the radiance which intensity has been reduced due to absorption and outscattering along the length S. Lsc is the radiance that has undergone a series of scattering events and ?nally scattered into a small cone around the observation direction ω.
PROPAGATOR FOR RADIATIVE TRANSFER . The solution of equation 2 is the limit of the corresponding solution for the time-dependent problem where radiance L varies over time t . It is convenient to express time in units of length s as t = s/v where v is the speed of light in the medium. With this notation, the time-dependent radiative transfer (TDRT) equation is

by ?rst ?nding solution G of the elementary transfer problems stated above and then forming the complete solution of equation 3 by using the superposition principle, i.e. integrating G with the initial radiance distribution:
?

L(s, x, ω) =

G(s, x ? x , ω, ω )L0 (x , ω )dx d ω ,

where G(s, x ? x , ω, ω ) is the evolution operator and L0 (x , ω ) is the initial distribution. Conceptually, the notion of the evolution operator is equivalent to the idea of radiative process introduced by Preisendorfer34 . Note that the boundary value problem for the Green function is actually adjoint to the radiative transfer problem, but due to reciprocity (time-reversal invariance) we can solve the light transport problem by reversing the direction of light propagation1 .
PATH INTEGRALS . Consider the problem of ?nding the probability that a point particle at the initial position xi and time ti will reach a ?nal position x f and time t f , i.e. will “propagate” (in the quantum mechanical sense) from xi to x f . This quantity can be expressed using quantum mechanical propagator G(t f ? ti , x f , xi ) which for this problem is the solution of the Schr?dinger equation subject to appropriate initial conditions. If the problem is broken down into a series of shorter time steps with propagators G(t , x, x ), the full propagator is expressed as a superposition of “smaller” Green function:

? + ω · ? + c(x) L(s, x, ω) = ?s b(x)
?



P(ω, ω )L(s, x, ω )d ω + Q(s, x, ω).

(3)

The solution of the TDRT equation can be formulated in terms of the Green evolution operator G which is also called the propagator, the Green function, or the point spread function (PSF). It is de?ned as the solution of homogeneous equation
? + ω · ? + c(x) G(s, x, ω, x , ω ) = ?s b(x)
?

G(t , xi , x f ) =
N →∞
?

lim

···
?

G(t /N , xi , x1 ) · · · G(t /N , xN ?1 , x f )d x1 · · · d xN ?1 .



P(ω, ω )G(s, x, ω , x , ω )d ω ,

(4)

with the initial condition
G(s = 0, x, ω, x , ω ) = δ(x ? x )δ(ω ? ω ).

Physically, the Green propagator G(s, x, ω, x , ω ) represents the radiance at point x in direction ω at time s due to light emitted at time zero by a point directional light source located at x shining in direction ω . For example, in the absence of scattering (b = 0), the solution for the propagator G is
G(s, x, ω, x , ω ) = δ(x ? ωs ? x )δ(ω ? ω ) exp ?
?

The object on the right-hand side is called a path integral. Feynman formulated quantum mechanics using path integrals10 and showed that with an appropriate de?nition of differential measure x in the in?nite-dimensional path space one can write for the quantum mechanical amplitude of a propagating point particle
¤ ?

x f , t f |xi , ti =
¤

? x(t )eiA[x(t )]/h

s 0

where the weight factor contains the classical action A(xi , x f , t f ? ti ) for each path. The classical action is the integral of the Lagrangian over the time the trajectory traverses. The path taken by a classical trajectory could be none other than the one that minimizes the classical action. R ADIATIVE T RANSFER A S A S UM OVER PATHS . The path integral (PI) approach provides a particular way to express the propagator G(s, x, ω, x , ω ). It is based on the simple observation that the full process of energy transfer from one point to another can be thought of as a sum over transfer events taking place along many different paths connecting points x and x , each subjected to boundary conditions restricting path directions at these points to ω and ω , respectively. The full propagator is then just an integral of individual path contribution over all such paths. This object is the path integral de?ned above. Note that this is different than the terminology used in Veach and Guibas45 .
c The Eurographics Association 2003.

a(x ? ω(s ? s ))ds

.

Here the light travels in a straight line and is attenuated by the absorption coef?cient a(x). One can see that in this case, the formulation using the propagator is generally equivalent to simple raytracing. The Green propagator G(s, x, ω, x , ω ) represents the angular distribution and density of rays at point x in direction ω generated by point x in direction ω and is therefore equivalent to raytracing. The concept of the Green function has been used in neutron transport theory3 providing an approach to solving radiative transfer problems with arbitrary boundary conditions

Premo?e et al. / Path Integration for Light Transport in Volumes

Because the integration is performed over the in?nitedimensional path space using the not very intuitive differential measure de?ned for it, the mathematics of path integrals is quite complex10, 21 . Tessendorf41 derived a path integral expression for the propagator G in homogeneous materials. Interested readers are referred to his further work42, 43 for detailed derivations of the path integral formulation. We present here a more intuitive approach suf?cient for our purposes and state results from the literature without derivations. 4. Practical Path Integrals 4.1. The Most Probable Path We ?rst attempt to formulate conditions to ?nd the most probable path (MPP) among all possible ones which satisfy the necessary boundary conditions. We will then assume that the full propagator can be suf?ciently approximated by accounting only for contributions from paths “close to” this special one. Formally, this approach corresponds to evaluating the path integral using WKB (WentzelKramers-Brillouin) expansion6 , which is a well-established method from perturbation theory. In practice, this means that once the MPP (or its approximation) is identi?ed, we simply consider its blurred contribution which approximates the effect of surrounding paths. The details of how this is done are presented in Section 5 below, but one can already see some potential advantages of this approach over more traditional methods. Compared with Monte-Carlo techniques which perform statistical sampling of random paths, the path integral approach attempts to ?nd the most important ones directly and can therefore be considered an extreme form of variance reduction. The PI formulation also does not explicitly rely on further assumptions about the character of radiance distribution which are needed, for example, in diffusion approximation11 . However, if warranted, we can take direct computational advantage of the fact that the radiance distribution becomes more and more blurred as one travels along the MPP. Consider an inhomogeneous medium with positiondependent scattering and absorption coef?cients b(x) and a(x). Let γP (s) be some arclength parameterized path from x to x. The probability density for a photon to reach x while traveling exactly along this path and not by any other possible one can be written as a product of two terms: the probability density of experiencing a series of scattering events which results in this particular path taken and the probability that the photon will “stay alive” at the end of its journey (i.e. not be absorbed along the path). We implicitly assume here that all absorption and scattering events are independent. This allows the second term to be expressed directly from the radiative transfer equation with no scattering which is written along the path in a trivial form dL(s)/ds = ?a(x)L(s) with initial conditions L(0) = 1 where we use loose notations L(s) for the radiance along the path. In this case, the
c The Eurographics Association 2003.

fraction of initial radiance which reaches the end of the path is exactly the probability of a photon not being absorbed. Therefore, using the solution of the equation above we can write
?

p(not absorbed) = exp ?

S 0

a(γP (s ))ds

,

where the argument to the integral is the absorption coef?cient along the path and S is total length of the path. To deal with the scattering term, we adopt an approach similar to Wilson and Wang46 which is inspired by physics, rather than attempting to follow the more mathematically rigorous treatment of Tessendorf42 . Homogeneous media have been considered so far in the literature and in many cases the rigorous mathematical procedures will break down if scattering/absorption coef?cients are allowed to vary across the medium.? A notable exception is the path construction using random walks by Pauly et al.27 . We ?rst split the path into a number of straight line segments connecting positions of scattering events that change photon propagation direction. Then the probability density of a particular path is a product of probability densities that individual scattering events will change the propagation direction “just right” to steer the photon all the way along the path. Using the phase function P(?θ) where ?θ is the change of propagation direction, we get for the total probability density by writing : p(path shape) x =
¤

i∈γP (s)



P(?θi )d ωi ,

(5)

where individual factors correspond to the different (ith ) scattering event along the path. Differential solid angles will eventually become a part of the full differential measure x in the path space and are not of interest for ?nding the MPP since they do not affect relative probabilities of different paths. We now make further assumption that the phase function is strongly peaked in the forward direction which is true for many important media5 . In this case, P(θ) can be approximated with the ?rst terms of its Taylor expansion as 1 ? αθ2 (we drop irrelevant constants here). It can be shown that if we want to maintain the mean cosine of the scattering angle g, α has to be set to α = 1/4(1 ? g) = 1/(2 θ2 ) where θ2 is the mean square scattering angle. Note that although the phase function is strongly forward peaked, this does not mean that the path itself has to deviate by only a small angle from its original propagation direction, which is an assumption often used to simplify derivations. We would also like to treat path γP (s) as a continuous object by taking a limit in equation 5 as scattering events occur often enough along the path. Each scattering event changes the propagation direction by a small amount δθ. In the case of forward peaked phase functions, only such scattering events generally have
¤

? For example, the Fourier transformed RT equation will have a much more complex form in this case, containing convolution over frequencies.

Premo?e et al. / Path Integration for Light Transport in Volumes

signi?cant probability. The expression of interest is
i∈γP (s)



P(?θi ) ≈

i∈γP (s)



1 ? α(δθi )2 =

i∈γP (s)



1 ? αδsi

δθi δsi

2

δsi

,

(6)

will tend to curve more in regions with high scattering coef?cient b(x). Finally, for a homogeneous medium, the MPP will be completely determined by the applied boundary conditions. Explicit expressions can be obtained in some cases using the standard Euler-Lagrange minimization procedure applied to the integral in equation 742 and we will use such results below. 4.2. Spatial Variability Optical properties in a scattering medium can vary arbitrarily spatially. A spatial variability in a medium can be measured by the number of scattering events that occur along the path γP :
?

where we introduced the lengths of path elements between scattering events δsi . Note that by construction these segments are physically constrained to be of ?nite length of the order of 1/b(γP (si )) and can not be simply treated as in?nitely short in the limit since no scattering can physically occur on an in?nitesimal path interval. We also intentionally re-arranged the last expression to highlight its part in square brackets which indeed can be treated as a full differential of some function in our approximation. Taking a logarithm of equation 6, using Taylor series, and replacing one of the δsi with its physical value we obtain

(s) =

s 0

b γP (s ) ds .

i∈γP (s)



ln 1 ? αδsi

δθi δsi

2

δsi

≈ ?α b(γP (si )) δθi δsi
2

i∈γP (s)



δsi .

In uniform media, the spatial variability is just a constant multiple of the distance s: = bs. Given the spatial variability of a path in inhomogeneous media, we can “invert” this equation and write displacement of the ray from its origin x with respect to as:
?

If the scattering scale δs = 1/b is much smaller than the macroscopic scale of the path, the expression in square brackets can be taken to give |d θ/ds|2 ds in this limit and the discrete sum will become an integral along the path. Taking the exponent, we get:
?

d( ) =

0

1 β( )d . b(x + d( ))

(8)

p(path shape) ? exp ?

s 0

α dθ b(γP (s )) ds

2

ds

.

This expression suggests a practical way of constructing actual path in inhomogeneous medium by stitching together straight line segments with lengths given by the local scattering coef?cient. The only extra information we need is local propagation direction which is discussed in the next subsection. 4.3. Finding the MPP Tessendorf43 described the propagation direction β( ) with Euler rotation angles and satisfying boundary conditions (equation 2) using the Fourier series expansion of the angles. We follow a simpli?ed version of his formulation to construct the stationary path β0 which is the path that minimizes attenuation along its length. We mentioned before that due to reciprocity we can construct the path by reversing the direction of light propagation1 . Through the rest of the paper we take advantage of this property and construct the MPP starting from the initial (viewing) direction ω and ending in the ?nal (light) direction ω , although the light is actually moving in the opposite direction. Let R be a rotation matrix that rotates initial direction ω to the z-axis vector z = (0 0 1)T . If the ?nal direction ω is written in spherical coordinates as ? ? sin (θ) cos (φ) ω = ? sin (θ) sin (φ) ? , cos (θ) then the stationary path β0 (s, z, ω ) that uniformly rotates ω

We can write this because of a general property of function limits: if the limit exists, its value does not depend on the speci?c way to take the limit and therefore the particular subdivision of the path we use does not affect the ?nal result. It is comforting to note this expression gives exactly the result of rigorous treatment of Tessendorf39 when applied to a homogeneous medium. The full expression for path probability density is now
p(path) ? exp ?
?

s 0

a(γP (s )) +

α dθ b(γP (s )) ds

2

ds

= (7)

exp (?A( path)),

where A is analogous to classical action along the path. To ?nd the MPP, we need to determine the path which minimizes the effective attenuation in equation 7. It is not possible to write an analytic expression for arbitrary functions a(x) and b(x), but some important general trends can nevertheless be established. Since the expression under the integral is non-negative, equation 7 favors shortest paths with low curvature d θ/ds. For example, if we are interested in the MPP connecting two points without specifying any extra conditions, this will be the straight line connecting these points. If a path has to turn to satisfy boundary conditions, it

c The Eurographics Association 2003.

PSfrag replacements

PSfrag replacements PSfrag replacements x x d d x x x x ω x x ω x θ( ) = 0 θ ( ) = π / 8 ω ω PSfrag replacements ω θ( ω ) = π/4 ω ω ω θ( ) = 0 d θ( ) = 0 θ( ) = π/8 x θ( ) = π/8 ω θ( ) = π/4 θ( ) = π/4 x ω ω x ω

Premo?e et al. / Path Integration for Light Transport in Volumes

ω

d

direction ω , the MPP is constructed by integrating the velocity function β0 . The locations on the MPP are found in terms of displacements d( ) (equation 8) from the starting point x. The displacements d( ) are de?ned implicitly and at ?rst it appears that the spatial variability along the entire path is needed. Note, however, that each displacement d( ) along the path only depends on the value of up to this point and not on the spatial variability of the entire path, allowing an “incremental” construction of the path. The path is constructed by stitching together piecewise linear path segments. We march along the path in steps of size d (step size in spatial variability not in distance). This step size d is arbitrary and is analogous to selecting a step size ds for direct lighting computation when marching along straight line. A sensible value for d can also be estimated from the optical properties and density of the volume. At every step we ?rst obtain the propagation direction β( ) along the existing portion of the path using equations 11 and 10 with the accumulated total spatial variability acc on the path so far substituted for . We then use equation 8 to compute the displacement point d( acc ) along the path.Finding the next point on the path involves increasing the spatial variability of the path so far by d ( acc = acc + d ) and ?nding corresponding displacement d( acc ) from the initial point x. Note that the full path is reconstructed from scratch at each step. So, for every sampling point on the curved path, the MPP is reconstructed from the initial point x and not just from the previous sampling point on the path. Note that at the starting point x the spatial variability is zero (no scattering events encountered so far) which therefore causes the ?rst path segment to be in the initial direction ω. Similarly, the very last segment of the path will by construction be in the direction ω , matching the boundary conditions. The propagation direction of path segments between starting point x and ending point x uniformly (in ) rotates from the initial direction ω to the ?nal direction ω . The segment length in physical space depends on the scattering coef?cient b at previous displacement point d( ). The total spatial variability along the path is just the sum of spatial variabilities d along each segment on the path. This is again analogous to computing the physical length s by summing segments ds along the straight line. As expected, the stationary path β0 is relatively ?at in regions where scattering coef?cient b(x) is small and highly curved where density is high. Figure 3 illustrates construction of the most important path using the described method. Quaternions could provide an alternative and more rigorous approach to uniform rotations. 4.4. Multiple Scattering Phase Function Tessendorf and Wasson44 observe that the width of the phase function grows with the number of scattering events . When the number of scattering events grows large, the probabil-

ω ω θ( ) = π/4 θ( ) = π/8 θ( ) = 0

Figure 3: Construction of the Most Probable Path in a homogeneous medium. The path is constructed by stitching together piecewise linear path segments and integrating the propagation direction β. The propagation direction of path segments between starting point x and ending point x uniformly rotates from the initial direction ω to the ?nal direction ω , therefore matching the boundary conditions.

where

into z over the physical path length S is ? ? sin (θ(s)) cos (φ) β0 (s, z, ω ) = ? sin (θ(s)) sin (φ) ? , cos (θ(s)) θ(s) = θ ? s θ. S

(9)

Note that we still need to apply rotation R?1 that rotates z back to ω such that: β0 (s, ω, ω ) = R?1 β0 (s, z, Rω ). (10)

One can show that such a “uniformly turning” path is exactly the MPP among all paths of ?xed length for the homogeneous medium, i.e. it minimizes effective attenuation given by equation 7. If the stationary path β0 needs to be constructed in inhomogeneous medium with respect to path spatial variability then we will simply replace θ(s) in equation 9 by θ( ) where the light path γP is now parameterized according to the number of scattering events and not the physical path distance: θ( ) = θ ? θ. (11)

Given the starting point x, initial direction ω, and ?nal
c The Eurographics Association 2003.

Premo?e et al. / Path Integration for Light Transport in Volumes

ity of scattering in any direction is equal and the phase function essentially becomes isotropic. It follows from the WKB approximation that the average scattering angle ΘMS after scattering events is: ΘMS = Wasson44 1 ? exp (? ) .

form. In practice, the direct and indirect radiance components are computed separately and if the phase function allows backscattering, it will be computed explicitly in the direct radiance component computation. 5. Algorithm As discussed in section 3, the radiance L received from direction ω at the observation point x is composed of three components: L(x, ω) = Lun (x, ω) + Lsc (x, ω) + Le (x, ω). The unscattered component Li un (x, ω) represents the amount of unscattered light due to the ith light source:
?

Tessendorf and also introduce the idea of the multiply-scattered phase function which is de?ned as the probability of light scattering through an angle θ after scattering events: 1 PMS (θ, ) = P θ N , (12)

1 ? exp (? )

where N is the normalization constant such that 4π PMS (θ, )d ω = 1. Intuitively, equation 12 says that the probability of scattering into an arbitrary angle increases with the number of scattering events. When the number of scattering events is large, the phase function PMS becomes isotropic. Note that equation 12 holds for arbitrary phase function.
?

i Li un (x, ω) = Llight (x, ω) exp ?

∞ 0

c(x ? ωs)ds ,

(15)

where Li light (x, ω) is the radiant exitance in direction ω from th the i light source. In practice, the unscattered radiance Lun and the emitted radiance Le are also the source for the scattered radiance13 :
QS (x, ω) = a(x)Le (x, ω) + b(x)
allNlights

4.5. WKB Approximation For Multiple Scattering Tessendorf41 derived a path integral expression for the propagator G in homogeneous materials:
G(s, x, ω, x , ω ) =
?


i

P(ω, ωi )Li un (x, ωi ). (16)

[d β ][d p]δ(β (0) ? ω )δ(β(s) ? ω)
s
?

δ x? exp b
?

0

ds β(s ) exp (?cs)
s
?

s 0

ds Z ( p(s)) exp i

0

ds p(s ) ·

d β(s ) , ds

(13)

To compute the total radiance L in the medium, all external and internal sources of radiation need to be propagated through the volume to the point x on the volume that the camera is looking at. The evolution operator G from equation 4 propagates all energy to a given observation point and direction. We use the propagator G from equation 14 for the rendering algorithm described in the next subsection. We also use results from Tessendorf and Wasson44 to develop a rendering algorithm using the path integration formulation. Monte Carlo Ray Tracing Monte Carlo ray tracing is an accurate algorithm for solving the radiative transfer equation in arbitrary media. We use it for comparison to evaluate our approximations. We march through the medium in direction ω sampling points along the ray. The light from previous step is attenuated and the light that is inscattered into the viewing direction ω is gathered. The inscattered light is collected recursively for each inscattered ray: Ln+1 (x, ω) =
allNlights

where Z ( p) is the Fourier transform of the normalized Gaussian phase function PG and p is the Fourier transform variable. Interested readers are referred to Tessendorf 42, 43 for detailed derivations of the path integral formulation. If the phase function P is not spatially varying and the single scattering albedo ω0 = b c is smoothly varying in space then the propagator G can be extended to inhomogeneous materials44 , but no details are provided. The WKB approximation done by Tessendorf42, 43 consists of approximating the path integral in equation 13 by ?nding the most probable (stationary) path and integrating all paths ?uctuating around the MPP. The path integration approximation they obtained can be expressed in terms of the Green propagator G:
G(s, x, ω, x , ω ) ? exp ?
?

s 0

ds c(γP (s ))

e ? 1 PMS (θ, ) (14)



The propagator in equation 14 is valid for spatially varying materials as long as the density of the material is smoothly varying and the phase function is uniform within the volume. While the WKB approximation was obtained using the assumption of strongly forward peaked phase functions, there is no other restriction on the speci?c phase function

4π Lsc (x, ωi )P(x, ωi , ω)b(x)?x + N i∑ =1 exp (?c(x)?x)Ln (x + ω?x, ω) where M is the number of directional samples taken. While Monte Carlo ray tracing is robust and powerful, it is also slow because of the large number of rays needed.
c The Eurographics Association 2003.

l M

Lun (x, ωl )P(x, ωl , ω)b(x)?x +

Premo?e et al. / Path Integration for Light Transport in Volumes

5.1. Path Integration Approximation Our algorithm exploits the WKB approximation from previous section. The WKB approximation computes the multiply-scattered light by ?nding the most probable path and then analytically integrates scattered radiance along this path and some neighborhood around this path including quadratic ?uctuations around this path (equation 14). The approximate radiance is then the sum of the singly scattered (“direct”) and multiply scattered components44 :
Lssc (x, ω) =

?

Algorithm 1 Rendering Algorithm Using WKB Approximation
// Preprocessing for ef?ciency purposes for each light source S Compute a light attenuation through the media by ray marching Store light attenuation in a spatial data structure (e.g. octree, deep shadow map) end for // Unscattered (direct) radiance Lssc computation Given starting point x and viewing direction ω for each sample point along the ray in direction ω Compute the unscattered radiance contribution from equation 17 end for // Scattered (indirect) radiance Lmsc computation // using WKB approximation Choose number of important paths nw Choose the sampling length d for each important path βnw for each sample point along the path βnw Find ds and a point x on the path from equation 8 for each set of sampling directions ? Compute the scattered radiance contribution from equation 18 end for end for Add diffuse light source contributions (equation 20) end for

0

QS (x ? sω, ω)) exp ?
?

s 0

ds c(x ? (s ? s )ω) ds,

(17)

and
Lmsc (x, ω) =
i=1 4π nw ∞
?



?



0

dsi {QS (x , ω ) × w[x ] × PMS (θ, )},

(18)

where w[x ] = exp (?c(x ) /b(x )) × (exp ( ) ? 1), (19)

and x is an intermediate point on the curved path: x = x + d( ). The ?rst term in equation 19 is the intensity of the ray at a sampling point along the path while the second term is the number of rays that propagate along this path. Note that the number of rays grows exponentially with the distance. Brute force Monte Carlo ray tracing algorithms would have trouble with this due to enormous number of potentially spawned rays. The path integration formulation does not have problems with this exponential growth, because the rays are dealt with via analytic integration. As we saw earlier, the MPP we construct is not uniquely de?ned by the boundary conditions. The total length of the path (either arclength S or total spatial variability ) remains a free parameter since it is not possible to ?nd this parameter without preprocessing or involved analysis. To ?nd the true MPP one would need to perform a search among paths with different lengths and choose the one with minimal effective attenuation. This process can be very expensive for inhomogeneous media. Fortunately, in practice it is suf?cient to construct a small number nw of “quasi-MPPs” with different total lengths (spatial variabilities) using the process described above and sum their contributions, as shown by equation 18. Based on the density and optical properties of the volume, spatial variability of the path is estimated. For example, a raymarching can be performed along three axes and diagonals to estimate the total length of the path parameter. These estimates are then used as a free parameter for the MPP generation. Further investigation is needed to estimate the total path length more robustly and systematically. However, we found that na?ve approach worked well in practice without spending extra time probing the volume. The contributions of diffuse (or indirect) outside light
c The Eurographics Association 2003.

sources are included as a radiance source at the boundary:
L(x, ω) = exp {?
?

sboundary 0

ds c(x ? ω(s ? s ))}Ld (ω)+ d ω wb [x ]PMS (θ, b )Ld (ω ), (20)

+∑
?

nw



where Ld is the diffuse radiance in a given direction. Similarly to the weight w in equation 19 for the collimated beam13 , the weight for the diffuse light source contribution is: wb [x ] = exp (?c b /b) × (exp ( b ) ? 1), (21)

and b is the number of scattering events when the path reached the boundary x of the volume. The rendering algorithm can be implemented using an existing ray marching routine with little dif?culty. For every viewing ray ω, we ?rst ?nd an intersection with the volume. This is starting point x. The direct component (equation 17) is computed using standard ray marching. At every sample point, the source term (direct radiance and emission) QS is computed and the intensity is then exponentially attenuated. To speed up computation we precompute light attenuation in the volume and store it at all possible depths. For this we use a simpli?ed version of a deep shadow map23 . From each light source we compute attenuation along the ray and store it at every sampling point. This speeds up rendering, because we

Premo?e et al. / Path Integration for Light Transport in Volumes

only pay the price once and the rest is a table lookup. Therefore, the source term QS is just a lookup in the spatial data structure storing the precomputed attenuation for every light source. The indirect component (equation 18) is computed by constructing the most probable path (or some number of them) using the procedure presented in Section 4.3. We march along the curved path in steps such that the number of scattering events is uniform. At every sampling point along the curved path, we compute the contribution from many directions over the entire sphere. Sampling directions are chosen based on viewing and lighting directions (i.e. sample more densely in directions close to the viewing and lighting directions). We have a ?xed number of directions we sample at every point. One can use a more sophisticated metric to determine how many directions to sample. In our implementation, the number of sampling directions is a parameter to the rendering algorithm. Although many directions are sampled at every sampling point along the path, this consists just of a few operations since no rays are actually spawned. This is unlike in a traditional Monte Carlo algorithm where additional rays are created at every sampling point which causes exponential growth in number of rays. Similarly to computing the direct component, the source term QS can be just looked up in the spatial data structure. Once the boundary x of the volume is reached, the diffuse light contributions are explicitly added (equation 20). The rendering algorithm is summarized in Algorithm 1. 6. Results and Discussion The algorithm from the previous section was implemented in C++ (compiled with g++ compiler) as an extension to the traditional Monte Carlo raytracer. All examples were rendered on a Pentium IV 1.7 Ghz with 512 Mb of memory running Linux. The results of the algorithm for a uniform density volume is shown in Figure 4. These images used the scattering parameters for skim milk from Jensen et al.18 and Henyey-Greenstein phase function with g = 0.9. The images are 512 by 384 and used 9 samples per pixel. The runtimes were 25 minutes for the Monte Carlo method and 8 minutes for the PI method. Note that the path integration solution is slightly darker than the brute-force Monte Carlo solution. This is because path integration assumes that paths near the MPP are where all of the light comes from, ignoring contributions that are far from the MPP. However, the degree of darkening is small, showing that the contributions from those paths is small as well. Figure 5 shows the results for a non-uniform density. The images are 512 by 384 and used 9 samples per pixel. The runtimes were 55 minutes for the PI method and and 580 minutes for the Monte Carlo solution. The cloud densities were generated using a Perlin-style noise function, and a Gaussian phase function was used with ? = 0.11. Again, notice some slight darkening. Figure 6 shows the results for a non-uniform density cloud

Figure 4: Upper: MC raytracing. Lower: PI Approximation.

with different illumination that includes indirect illumination and skylight contribution. The image is 640 by 480 pixels and used 9 samples per pixel. The runtime was 95 minutes. Most of the theory in this paper is known to other scienti?c communities. We have described a more intuitive way of constructing path integrals avoiding the heavy mathematics of functional integration that would ordinarily be required. We have shown that this theory can be applied to produce reasonably accurate images. This was accomplished by ?nding the MPP and explicitly computing contributions along that path. All other paths are dealt with via analytic integration around the MPP. The MPP can be also used in traditional Monte Carlo simulations as a form of variance reduction. Our formulation can be used as a starting point for a sophisticated method such as Metropolis Light Transport, which effectively attempts to ?nd the MPP by random walk. We have not pursued a highly ef?cient implementation; there is much room for the development of algorithms based on path integration. For example, a hierarchical implementation analogous to 15 would surely be bene?cial. Our images should be viewed as a proof-of-concept that path integration is a valuable tool for dealing with volumetric lighting. We have not shown that it has advantages over the dipole approximation18 to the diffusion approxic The Eurographics Association 2003.

Premo?e et al. / Path Integration for Light Transport in Volumes

complement to diffusion, and should become part of the toolbox of dealing with volume-based scatterers. Many classical rendering problems could be reformulated and solved faster using the path integral formulation. Understanding the parallels with the classical formulation and thinking of problems in this framework opens many directions for future research. 7. Acknowledgements Thanks to Milan Ikits, Emil Praun, and David Weinstein for reading drafts of the paper. The reviewers provided valuble suggestions for improving this paper. This work was partially supported by NSF grant 99-78099. References
1. Raphael Aronson. Radiative transfer implies a modi?ed reciprocity relation. J. Opt. Soc. Am. A, 14(2):486–490, February 1997. James Arvo. Transfer equations in global illumination. Global Illumination, SIGGRAPH ‘93 Course Notes, August 1993. G. I. Bell and S. Glasstone. Nuclear Reactor Theory. Van Nostrand Reinhold, New York, 1970. Craig F. Bohren. Clouds in a glass of beer. John Wiley and Sons, 1987. Craig F. Bohren. Multiple scattering of light and some of its observable consequences. American Journal of Physics, 55(6):524–533, June 1987. Masud Chaichian and Andrei Demichev. Path Integrals in Physics Volume I. Stochastic Processes and Quantum Mechanics. Institute of Physics Publishing, 2001. C. C. Constantinou and C. Demetrescu. Physical interpretation of virtual rays. J. Opt. Soc. Am. A, 14(8):1799–1803, 1997. J. Dorsey, A. Edelman, H. W. Jensen, J. Legakis, and H. Pedersen. Modeling and rendering of weathered stone. In Proceedings of SIGGRAPH 1999, pages 225–234, August 1999. Ronald Fedkiw, Jos Stam, and Henrik Wann Jensen. Visual simulation of smoke. In Proceedings of ACM SIGGRAPH 2001, pages 15–22, August 2001.

2. 3. 4.

Figure 5: Upper: MC raytracing. Lower: PI Approximation.

5.

6.

7. 8.

9.

10. R. P. Feynman and A. R. Hibbs. Quantum mechanics and path integrals. McGraw-Hill Higher Education, 1965.

Figure 6: A cloud rendered using the Path Integration approximation. Diffuse light source contributions were added at the boundary of the volume.

11. R. G. Giovanelli. Diffusion through non-uniform media. Progress in Optics. North-Holland, Amsterdam, 1963. 12. Eugene P. Gross. Path integral approach to multiple scattering. J. Math. Phys., 24(2):399–405, 1983. 13. Akira Ishimaru. Wave Propagation and Scattering in Random Media, Volume I: Single scattering and transport theory; Volume II: Multiple scattering, turbulence, rough surfaces and remote sensing. Academic Press, New York, 1978. 14. Steven L. Jacques and Xuijing Wang. Path integral description of light transport in tissues. In Proceedings of optical tomography and spectroscopy of tissue: theory, instrumentation, and human studies II, pages 488–499, 1997.

mation for dense uniform media, and in fact there is reason to prefer diffusion as density increases because path integration ignores backscattering. However, we have shown that path integration can be applied to sparse and spatiallyvarying volumes where the diffusion approximation is not appropriate 19 . Thus path integration should be viewed as a
c The Eurographics Association 2003.

Premo?e et al. / Path Integration for Light Transport in Volumes 15. Henrik Wann Jensen and Juan Buhler. A rapid hierarchical rendering technique for translucent materials. ACM Transactions on Graphics, 21(3):576–581, July 2002. 16. Henrik Wann Jensen and Per H. Christensen. Ef?cient simulation of light transport in scenes with participating media using photon maps. In Proceedings of SIGGRAPH 98, pages 311– 320, July 1998. 17. Henrik Wann Jensen, Justin Legakis, and Julie Dorsey. Rendering of wet materials. In Eurographics Rendering Workshop 1999, Granada, Spain, June 1999. Springer Wein / Eurographics. 18. Henrik Wann Jensen, Stephen R. Marschner, Marc Levoy, and Pat Hanrahan. A practical model for subsurface light transport. In Proceedings of SIGGRAPH 2001, pages 511–518, August 2001. 19. Jan J. Koenderink and Andrea J. van Doorn. Shading in the case of translucent objects. In Bernice E. Rogowitz and Thrasyvoulos N. Pappas, editors, Proc. SPIE Vol. 4299 Human Vision and Electronic Imaging VI, pages 312–320, 2001. 20. Eric Lafortune. Mathematical Models and Monte Carlo Algorithms. PhD thesis, Katholieke Universiteit Leuven, 1996. 21. Flor Langouche, Dirk Roekaerts, and Enrique Tirapegui. Functional Integration and Semiclassical Expansions. D. Reidel Publishing Company, 1982. 22. Hendrik P. A. Lensch, Michael Goesele, Phillipe Bekaert, Jan Kautz, Marcus A. Magnor, Jochen Lang, and Hans Peter Seidel. Interactive rendering of translucent objects. In Proceedings of Paci?c Graphics, 2002. 23. Tom Lokovic and Eric Veach. Deep shadow maps. In Proceedings of ACM SIGGRAPH 2000, pages 385–392, July 2000. 24. Steven D. Miller. Stochastic construction of a feynman path integral representation for green’s functions in radiative transfer. Journal of Mathematical Physics, 39(10):5307–5315, 1998. 25. S. G. Narasimhan and S. K. Nayar. Shedding light on the weather. In Proceedings of CVPR, June 2003. 26. Duc Quang Nguyen, Ronald P. Fedkiw, and Henrik Wann Jensen. Physically based modeling and animation of ?re. ACM Transactions on Graphics, 21(3):721–728, July 2002. 27. Mark Pauly, Thomas Kollig, and Alexander Keller. Metropolis light transport for participating media. In Rendering Techniques 2000: 11th Eurographics Workshop on Rendering, pages 11–22, June 2000. 28. Lev T. Perelman, Jun Wu, Irving Itzkan, and Michael S. Feld. Photon migration in turbid media using path integrals. Physical Review Letters, 72(9):1341–1344, 1994. 29. Lev T. Perelman, Jun Wu, Yang Wang, Irving Itzkan, Ramachandra R. Dasari, and Michael S. Feld. Time-dependent photon migration using path integrals. Physical Review E, 6(51):6134–6141, 1995. 30. Frederic Pérez, Xavier Pueyo, and Fran?ois X. Sillion. Global illumination techniques for the simulation of participating media. In Eurographics Rendering Workshop 1997, pages 309– 320, June 1997.
c The Eurographics Association 2003.

31. Matt Pharr and Patrick M. Hanrahan. Monte carlo evaluation of non-linear scattering equations for subsurface re?ection. In Proceedings of SIGGRAPH 2000, pages 75–84, July 2000. 32. A. Ya. Polishchuk and R. R. Alfano. Fermat photons in turbid media: an exact analytic solution for most favorable paths – a step toward optical tomography. Optics Letters, 20(19):1937– 1939, 1995. 33. A. Ya. Polishchuk, J. Dolne, F. Liu, and R. R. Alfano. Average and most-probable photon paths in random media. Optics Letters, 22(7):430–432, 1997. 34. Rudolph W. Preisendorfer. Light Transport on Discrete Spaces. Pergamon Press, 1965. 35. S. Premo?e. Analytic Light Transport Approximations for Volumetric Materials. In Proceeding of Paci?c Graphics, 2002. 36. Holly Rushmeier. Rendering participating media: Problems and solutions from application areas. In Fifth Eurographics Workshop on Rendering, pages 35–56, Darmstadt, Germany, June 1994. 37. Jos Stam. Multiple scattering as a diffusion process. In Eurographics Rendering Workshop 1995. Eurographics, June 1995. 38. V. I. Tatarskii, A. Ishimaru, and V. U. Zavorotny. Wave Propagation in Random Media. Society of Photo-optical Instrumentation Engineers, Bellingham, Washington, 1993. 39. Jerry Tessendorf. Radiative transfer as a sum over paths. Phys. Rev. A, 35:872–878, 1987. 40. Jerry Tessendorf. Comparison between data and small-angle approximaions for the in-water solar radiance distribution. Journal of Optical Society of America A, 5:1410–1418, 1988. 41. Jerry Tessendorf. Time-dependent radiative transfer and pulse evolution. J. Opt. Soc. Am. A, (6):280–297, 1989. 42. Jerry Tessendorf. The Underwater Solar Light Field: Analytical Model from a WKB Evaluation. In Richard W. Spinrad, editor, SPIE Underwater Imaging, Photography, and Visibility, volume 1537, pages 10–20, 1991. 43. Jerry Tessendorf. Measures of temporal pulse stretching. In Gary D. Gilbert, editor, SPIE Ocean Optics XI, volume 1750, pages 407–418, 1992. 44. Jerry Tessendorf and David Wasson. Impact of multiple scattering on simulated infrared cloud scene images. In Proceedings SPIE. Characterization and propagation of sources and backgrounds, volume 2223, pages 462–473, 1994. 45. Eric Veach and Leonidas J. Guibas. Metropolis light transport. In SIGGRAPH 97 Conference Proceedings, pages 65–76, August 1997. 46. Michael J. Wilson and Ruikang K. Wang. A path-integral model of light scattered by turbid media. J. Phys. B: At. Mol. Opt. Phys, 34:1453–1472, 2001.


赞助商链接

更多相关文章:
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图