Rendering Participating Media
with Bidirectional Path Tracing
Eric P. Lafortune and Yves D. Willems

Paper presented at the 7th Eurographics Workshop on Rendering

Please send any correspondence to:Eric Lafortune
Cornell University
Program of Computer Graphics
580 Rhodes Hall
Ithaca, NY 14853-3801
Telephone: +1 - 607 - 255 - 4880
Fax:
+1 - 607 - 255 - 0806
E-mail: Eric.Lafortune@cs.kuleuven.ac.beRendering Participating Media
with Bidirectional Path Tracing
Eric P. Lafortune and Yves D. Willems
Department of Computer Science, Katholieke Universiteit Leuven
Celestijnenlaan 200A, 3001 Heverlee, Belgium
Eric.Lafortune@cs.kuleuven.ac.be

Abstract: In this paper we show how bidirectional path tracing can be extended
to handle global illumination effects due to participating media. The resulting
image-based algorithm is computationally expensive but more versatile than previous solutions. It correctly handles multiple scattering in non-homogeneous,anisotropic media in complex illumination situations. We illustrate its specificadvantages by means of examples.

1 Introduction
Most current global illumination algorithms ignore the influence of the medium throughwhich the light travels. In these illumination simulations, light is only emitted and reflected at surfaces; the absorption, scattering and emission of light as a result of smoke,
dust, fog or flames are neglected. While this is a reasonable approximation for mostdaily scenes, the presence of media can sometimes enhance the realistic appearance ofrenderings. A typical example is the beam of light from the sun shining through a window of a dusty interior. In other applications such as visibility studies for traffic or for
fire exits in buildings, the medium even plays a central role in the simulation problem.
Rushmeier [1] and others [2, 3] have presented extensive overviews of research
efforts and packages for the simulation of light in participating media. Dependingon the applications they handle homogeneous or non-homogeneous and isotropic or
anisotropic media. Many techniques only model single scattering, which suffices for
optically thin media. Most approaches can be classified into a few basic categories:
zonal methods for isotropic scattering,
methods, discrete ordinate methods, point
collocation methods and image-based Monte Carlo methods. A zonal method that handles participating media was presented by Rushmeier [4], as an extension to the basic
radiosity algorithm. The algorithm assumes that the participating medium is isotropic
and computes the exchange of 
light between the discretised volume elements and thesurface elements in the scene.
methods also discretise the volume into finite elements, but now represent the radiance function as a series of spherical harmonics. They
then solve the problem in its differential form, e.g. [5, 6]. Discrete ordinate methods
discretise not only the original volume into finite volume elements, but also the directions of the sphere into finite solid angles, e.g. [7]. This lifts the assumption of thezonal methods that the participating medium is isotropic, at the cost of increased memory requirements. Languénou et al. [2] and Max [3] present algorithms to efficiently
distribute power between the discrete elements. Point collocation methods use higher
order basis functions, such as spherical harmonics, to represent the illumination function over the sphere, e.g. [8].

Monte Carlo methods are generally either an extension to the path tracing algorithm
as proposed by Kajiya [9] or to light tracing algorithms, as pioneered by Arvo [10].
Rushmeier [4] for instance presents an extension to the path tracing algorithm. It traces
random walks starting from the eye point, reflecting at surfaces and additionally scattering in the medium. Roysam et al. [11] and Rozé et al. [12] apply light tracing approaches, tracing random walks from the light sources, until they hit a detector. Pattanaik [13] follows this approach in the context of a zonal method, updating the radiosity values at discrete volume and surface elements.
In this contribution we propose an image-based Monte Carlo solution, now applying
bidirectional path tracing rather than just path tracing. The approach thus integrates previous light shooting and light gathering techniques for participating media. Section 2
recalls the mathematical model of global illumination for scenes with a participating
medium. We discuss the case of non-emitting media with possibly anisotropic scattering, starting with homogeneous media. Section 3 explains how a random walk throughthe scene can be sampled on the basis of the mathematical framework. Section 4 thenshows how these random walks integrate in the bidirectional path tracing algorithm.
This image-based approach is very general, like path tracing, but better suited to handle
complex illumination situations. Section 5 shows how the handling of homogeneous
media can be further extended to non-homogeneous media. Section 6 presents someresults to demonstrate the versatility of the algorithms. We indicate specific advantages
and disadvantages of the algorithm and draw some final conclusions in section 7.

2 Mathematical model
In its most commonly used form, the rendering equation defines the radiance function,
assuming that it is invariant along the direction in which it is radiated. The equation
can be extended fairly easily to account for participating media. In this case radiancechanges gradually along a line, as light is absorbed, outscattered to other directions andinscattered from other directions. The integral equation defining the radiance function
inside the medium becomes an integro-differential equation:


 "!#%$'&()

(1)

where:** *

is a point in the 3D space,*.+- is a direction,
is a differential distance along direction [, ],** /! is the sphere of directions, 10* )is a differential solid angle around direction
 107 [ ],is the field radiance function [ 243 ,65],* 8* & is the scattering coefficient [ 973:, ], is the extinction coefficient (often denoted by ;=< ) [ 9:37, ],  is the so-called phase function (often denoted by), which
expresses the fraction of radianceradiated
towards
point
fromdirectionthat10
is scattered in direction [ 9:3 ].
We do not considerof media such as flames or plasma. The scatteringemittance
may be anisotropic:
can be any function, as long as it is normalised over theL(x – sΘ, Θ′)L(x – dΘ, Θ′)L o(x – dΘ, Θ)L(x, Θ x)

Fig. 1. The radiance equation expresses how light in a medium is the result of light inscatteredalong the line and light reflected from the surface.

sphere of directions and reciprocal.
&8 If we& assumefor a moment that the participating  , then the original integroand 
medium is homogeneous, i.e.
differential equation can be rewritten as the following integral equation (fig. 1):(2)







   $ :    $':  "!# = 47 $ =  $= 

where is the distance to the nearest surface andis the outgoing radianceat that point. We can now also take into account the boundary conditions. The relation
between incoming and outgoing radiance at surfaces is expressed by the well-known
integral expression:

"8:     8()8  



cos 

*
where:*  :
*. - "8   is the self-emitted radiance at point in direction (/! 


(3)

10/

[ 243 , 5
], 10
is
the
bidirectionalreflectance
distributionfunction(BRDF)[
973 ],*
 is the hemisphere or sphere of incoming directions at point  ,and the surface normal at point .
 is the angle between direction
Equations (2) and (3) define the behaviour of light in the medium and at surfaces respectively. Due to their integral formulations they are well-suited for Monte Carlo evaluation, which we will discuss in the next section.

3 Tracing a random walk
The sampling of the radiance equations is fairly straightforward and has been describedfor instance by Rushmeier [4]. It is intuitive in the sense that it can be regarded as adirect simulation of light particles scattering at dust or fog particles and reflecting  atsurfaces until they are absorbed. In this case a random walk startsfrom the eye point ,in a direction pointing through a pixel. The radiance value
can be estimatedby performing a random walk. We can introduce a minor optimisation by first rewriting
equation (2) as a single integral:












 $ :    $':  "!# = 47 $ = 



!


"



 $': $#"7

$ +  &%"&'  = (4)

where8

is defined as:   8 "!#&Equation (4) shows that the factor
tion for the outer integral:&'  

medium:(5)can be used as a probability density func-&  

(6)

From a physical point of view, this Poisson distribution expresses the probability densityof a light particle interacting with the medium at distance . If the nearest surface along
the
ray,
term
 6$ at: a distance , lies further away than the point at distance 6, $ the= former
has to be estimated. Otherwise, the latter term
has tobe estimated. The ray tracing function, which determines the nearest surface, thus only
has to look for surfaces within a distance . For instance, in a foggy outdoor scene itis stochastically not always necessary to compute ray intersections with a village in the
distance. The a posteriori probability of the light particle reaching the nearest surface isequal to the transmittance:
=   (7)
surfaceBoth terms* expressions:
gral



8

/8andare estimated by sampling their respective inte-

is estimated by sampling the integral of equation (5) (fig. 2). A suitableis:subcritical PDF (i.e. it does not integrate to 9 ) for the scattered direction*scattering

 


 &(8)

"8

by equation (3). It consists of two terms. The self-emittedis:expressed8
radiance
can simply be evaluated. The second term is estimated by sampling the integral
expression (fig. 3). A common subcritical PDF for the reflected
directionis:  8  

&

c 
(9)
reflection) 

In both cases the radiance value
is then estimated for the sampled direction.
This recursive process results in a random walk through the scene. Because the latter
PDFs are subcritical the random walk terminates with probability 9 .
Stochastic ray tracing algorithms in general compute the direct illumination by explicitly sampling the light sources and tracing additional shadow rays, while the indirect
illumination is computed as before. The process for generating random walks remains
the same.

p medium(s)xp scattering(Θ′)

x – sΘ


Fig. 2. Scattering on the random walk. First the distance is sampled according to  medium  ,then the scattering direction  according to  scattering   .

p reflection(Θ′)

P surface(d)xFig. 3. Reflection on the random walk. If the distancea reflection instead, with an a posteriori probability ofdirection  is sampled according to  reflection   .

x – dΘ

lies past the nearest surface, there issurface  . In that case the reflection

light pathshadow rays

eye pathscreenlight source
eye point
Fig. 4. A schematic overview of bidirectional path tracing in a participating medium. A pair ofrandom walks is constructed: an eye path starting from the eye point, through a pixel, and a lightpath starting from a light source. The paths can scatter inside the medium or reflect at surfaces.
The points on the respective paths are connected by means of shadow rays, which determine thecontribution to the estimated flux of the pixel.

4 Bidirectional path tracing
Bidirectional path tracing has been introduced by us [14, 15] and has been presented
independently and improved by Veach and Guibas [16, 17]. It combines the ideas ofpath tracing and light tracing by creating random walks, starting not only from the eye
point, but also from the light sources. Figure 4 shows the basic idea. Eye paths start fromthe eye point, through the pixel that is being computed, as in classical path tracing. Light
paths start from a sampled point on a light source in a sampled direction. Both types ofrandom walks are generated as described in section 3. In the extension presented here,
rays can be scattered or absorbed in the medium, as well as being reflected or absorbed
at the surfaces of the scene.
After tracing a pair of random walks, the intersection points on the respective paths
are connected by means of shadow rays. Each shadow ray that is not intercepted byan object in the scene contributes to the estimated flux or average radiance through the
pixel.
The algorithm has proven to be particularly effective for rendering scenes with indirect illumination, compared to ordinary path tracing [18]. The different illuminationtransport paths present different estimators for the flux, each one with its own potential
strength. A judicious combination of the estimators [17] yields a robust estimator witha generally lower variance. For this purpose the alternative estimators for illuminationtransport with the same number of reflections have tosame pa  beconsidered1(    in  the in the 3D
rameter space, for instanceasthesequenceofpoints      space. Theof the light path, starting from a light source  section 
 is part is part of the eye path, starting from the eye.
Thesectionat point
point  .
The illumination contribution of the transport path in this parameter space is:       7  


 


c

 


   
       $   5 &%' (%)
   $   $ 

 5



  

!"#   $  

5

5

 6$

where we use thefor convenience.
 arrow notation8   +*
between point and point.
Thefunctionor the reflection at point on the random walk:8 

 ,*



-

  / 

 ,
 *
 .*


c


 



c

 


5

(10)

is the Euclidean distanceaccounts for the scattering  if the interaction is scattering,  if the interaction is reflection.
where    and   are the angles of the incoming and outgoing directions with the normal
at point .
density for generating
path with the aforementioned eye path / Theprobability this   and light path      is the product of their corresponding
probability densities, as the random walks are independent:         



 

 



  

 

 

(11)

These probability densities in turn can be defined recursively: 0

 

 




     

     1(  

(12)







where
for
  is the   probability  density  if   is known. The initial probability
 anddensities,,
depend on the sampling strategies for the pixel and the light sources respectively. The other probability densities are afunction of the PDFs used for constructing the randomwalk, as presented in section 3.
They depend on the types of interaction at points   and  :










with

  






and


       1(           $  






(13)

5

  scattering   reflection   
  if the interaction at point    is scattering, 
medium     surface$   
$    if the interaction at point   is scattering,

if the interaction at point

if the interaction at pointis reflection,is reflection.

Together with the appropriate probability densities (11), the illumination contribution (10) can be combined into an unbiased estimator, as discussed in [17]. We haveapplied the balance heuristic, for which the weights of the illumination contribution canbe written as follows:     (       









    1(                     1  (14)

where the summation is over all possible types of transport paths for which contributions are added.

5 Extension to non-homogeneous media
While the previous sections only discuss homogeneous participating media, extension
to non-homogeneous media is fairly straightforward. A common approach is to discretise the medium into volume elements, each having its own scattering characteristics.
This representation particularly suits zonal methods and other techniques that alreadydiscretise the scene to compute the illumination. In a ray tracing context the discretisation can be traversed while tracing a ray, treating each element as a homogeneous
medium, e.g. [13]. An alternative and slightly more general representation allows any
function to describe the scattering characteristics of the medium. Typically, this is anempirical noise function that mimics media such as smoke plumes and clouds of dust.
We have opted for the latter approach, additionally bounding the participating media
in the scene by means of simple volumes. For each ray piercing a volume with a participating medium the section of its traversal is determined. This section is then traversed
at discrete intervals where the optical characteristics of the medium are sampled. Using
this information a possible scattering can be sampled or the opacity can be estimated foreach interval, in which the medium is approximated as homogeneous. The handling of
non-homogeneous media thus becomes similar to the handling of homogeneous media.
The traversal obviously requires additional work, depending on the coarseness of thediscrete steps.

6 Implementation and results
We have implemented the presented extension in an existing bidirectional path tracing
program. The scattering characteristics of the medium are modelled using the phasefunction proposed by Schlick [19], which is energy-conserving and reciprocal, andwhich is computationally inexpensive to evaluate and to use as a PDF.
Figure 5 shows a test scene with a few geometric objects with various optical characteristics. In image (b) the  scene is engulfed in a slightly anisotropic cloud of smoke.
It has been rendered at  pixels, with 100 samples per pixel. The bidirectional
path tracing algorithm correctly renders the brightly lit lobe below the light source andthe indirect illumination resulting from it. The combination of objects and the participating medium is smooth.
Figure 6 shows a church interior. The sun casts a beam of light through the stainedglass window, which is only visible in the dusty atmosphere in image (b). The sun
is modelled as a spot light source outside. The stained-glass window is a translucent
surface, scattering the incident light partly diffusely and partly directionally. The illumination of the interior is entirely indirect as a result. Classical path tracing would bevery inefficient at rendering this scene. Only random walks that reach the outside would
get an illumination contribution, resulting in a high variance. Bidirectional path tracing
traces light paths to the inside, scattering and reflecting and eventually meeting the eyepaths there.
Figure 7 illustrates the possible use of bidirectional path tracing for determining
visibility in traffic scenes. In this case the algorithm also simulates the reflection andfocusing of light by the reflectors of the headlights. In image (b) most of the car isilluminated indirectly by light that is scattered in the fog. Rendering more complex
environments with traffic lights, traffic signs, puddles of water, etc. would be straightforward.
Figure 8 shows an interior scene with a fireplace, demonstrating the rendering of a
non-homogeneous participating medium. The smoke rising from the fireplace is illuminated from below. It thus participates in the illumination of the rest of the room.

7 Conclusions
We have discussed the extension of bidirectional path tracing to include the effects ofparticipating media in the simulation. The resulting algorithm correctly handles multiple scattering in non-emitting, non-homogeneous, anisotropic media. It provides unbiased estimators that are guaranteed to converge to the exact solution. Like most imagebased Monte Carlo approaches, the algorithm is versatile and requires little memory. It
is therefore well suited to render scenes with a combination of participating media andcomplex geometries and optical characteristics.
The required computation times may be taxing, however, even for the most hardenedMonte Carlo user. The example images took several hours to render, with the computations distributed over some workstations. The presence of participating media markedly
increases the complexity of the problem. Finite volume element techniques are not as
flexible, but they are undoubtedly faster and more effective for simple scenes. Compared to unidirectional Monte Carlo approaches, the bidirectional algorithm presents asignificant improvement. It handles scenes that were not tractable with classical pathtracing. It also offers a solution for the light tracing approaches used in some visibility studies with participating media, by not requiring the artificially large detectors forreceiving at least some light particles within a reasonable amount of time.

a)b)
Fig. 5. Test scene: (a) without participating media, (b) with a dusty atmosphere.

a)b)Fig. 6. Church interior: (a) without participating media, (b) with a dusty atmosphere.

a)b)Fig. 7. Traffic scene: (a) without participating media, (b) with a foggy atmosphere.

Fig. 8. Fireplace scene with a non-homogeneous participating medium.

References
1. H. Rushmeier, “Rendering participating media: Problems and solutions from application areas,” in Proceedings of the Fifth Eurographics Workshop on Rendering,(Darmstadt, Germany), pp. 35–55, June 1994.
2. E. Languénou, K. Bouatouch, and M. Chelle, “Global illumination in presence ofparticipating media with general properties,” in Proceedings of the Fifth Eurographics Workshop on Rendering, (Darmstadt, Germany), pp. 69–85, June 1994.
3. N. Max, “Efficient light propagation for multiple anisotropic volume scattering,”
in Proceedings of the Fifth Eurographics Workshop on Rendering, (Darmstadt,
Germany), pp. 87–104, June 1994.
4. H. Rushmeier, Realistic Image Synthesis for Scenes with Radiatively Participating
Media. PhD thesis, The Sibley School of Mechanical and Aerospace Engineering,
Cornell University, 1988.
5. J. Kajiya and V. Herzen, “Ray tracing volume densities,” Computer Graphics,
vol. 18, pp. 165–174, July 1984.
6. J. Stam, “Multiple scattering as a diffusion process,” in Proceedings of the Sixth
Eurographics Workshop on Rendering, (Dublin, Ireland), pp. 69–79, June 1995.
7. Ch. Patmore, “Simulated multiple scattering for cloud rendering,” in Proceedingsof ICCG, (Bombay, India), pp. 59–70, Feb. 1993.
8. N. Bhate and A. Tokuta, “Photorealistic volume rendering of media with directional scattering,” in Proceedings of the Third Eurographics Workshop on Rendering, (Bristol, UK), pp. 227–245, May 1992.
9. J. Kajiya, “The rendering equation,” Computer Graphics, vol. 20, pp. 143–150,
Aug. 1986.
10. J. Arvo, “Backward ray tracing,” Aug. 1986.
11. B. Roysam, A. Cohen, P. Getto, and P. Boyce, “A numerical approach to the computation of light propagation through turbid media: Application to the evaluationof lighted exit signs,” IEEE Transactions on Industry Applications, pp. 661–669,
May 1993.
12. C. Rozé, B. Maheu, and G. Gréhan, “Evaluations of the sighting distance ina foggy atmosphere by Monte Carlo simulation,” Atmospheric Environment,
vol. 28, no. 5, pp. 769–775, 1994.
13. S. Pattanaik, Computational Methods for Global Illumination and Visualisationof Complex 3D Environments. PhD thesis, NCST Birla Institute of Technology &Science, Pilani, India, Feb. 1993.
14. E. Lafortune and Y. Willems, “Bi-directional path tracing,” in Proceedings ofCompuGraphics, (Alvor, Portugal), pp. 145–153, Dec. 1993.
15. E. Lafortune and Y. Willems, “A theoretical framework for physically based rendering,” Computer Graphics Forum, vol. 13, pp. 97–107, June 1994.
16. E. Veach and L. Guibas, “Bidirectional estimators for light transport,” in Proceedings of the Fifth Eurographics Workshop on Rendering, (Darmstadt, Germany),
pp. 147–162, June 1994.
17. E. Veach and L. Guibas, “Optimally combining sampling techniques for Monte
Carlo rendering,” Computer Graphics, vol. 29, pp. 419–428, Aug. 1995.
18. E. Lafortune, Mathematical Models and Monte Carlo Algorithms for Physically
Based Rendering. PhD thesis, Katholieke Universiteit Leuven, Belgium, Feb.
1996.
19. Ch. Schlick, Divers Eléments pour une Synthèse d’Images Réalistes. PhD thesis,
Université Bordeaux 1, France, Nov. 1992.