UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Bidirectional importance sampling for illumination from environment maps Burke, David 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2005-0016.pdf [ 10.06MB ]
JSON: 831-1.0051747.json
JSON-LD: 831-1.0051747-ld.json
RDF/XML (Pretty): 831-1.0051747-rdf.xml
RDF/JSON: 831-1.0051747-rdf.json
Turtle: 831-1.0051747-turtle.txt
N-Triples: 831-1.0051747-rdf-ntriples.txt
Original Record: 831-1.0051747-source.json
Full Text

Full Text

Bidirectional Importance Sampling for Illumination from Environment Maps by David Burke B.Sc, University of Prince Edward Island, 2002 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Department of Computer Science)  We accept this thesis as conforming to the required standard  T H E UNIVERSITY OF BRITISH COLUMBIA October 18, 2004 © David Burke, 2004  Abstract  ii  Abstract Image-based representations for illumination are able to capture complex real-world lighting that is difficult to represent in other forms.  Current  importance sampling strategies for image-based illumination have difficulties in the case where both the environment map and the surface BRDF contain important high-frequency detail, for example, when a specular surface is illuminated by an environment map containing small light sources. We introduce the notion of bidirectional importance sampling, in which samples are drawn from the product distribution of both the surface reflectance and the energy in the environment map. Although this makes the sample selection process more expensive, we show significant quality improvements over traditional importance sampling strategies for the same compute time.  iii  Contents  Contents Abstract  ii  Contents  iii  List of Figures  v  Acknowledgements  vii  1  Introduction  1  2  Related Work  3  2.1  Environment Mapping  4  2.2  Prefiltered Environment Mapping  5  2.3  Alternate Bases and Precomputed Transport .  6  2.4  Importance Sampling and Point Relaxation  7  2.5  Combined Importance Sampling  9  3  . .  Monte Carlo Rendering  11  3.1  Environment Mapping  12  3.2  Monte Carlo Estimators and Sampling  15  3.3  3.2.1  Monte Carlo Integration  3.2.2  Sampling a General Distribution  Sampling the Environment Map  ;  15 19 24  Contents  iv  3.4  Improving Shadows: Area-weighted Lighting  27  3.5  Sampling from the BRDF  30  4 Sampling the Product Distribution  35  4.1  Bidirectional Importance Sampling  35  4.2  The Challenge of Bidirectional Sampling  42  4.3  Sample Generation Through Rejection  43  4.4  Sample Generation Through Sampling-Importance Resampling 47  5 Results  51  6  63  Conclusions  Bibliography  64  A Details on the Alias Method for Sampling  70  List of Figures  v  List of Figures 2.1  Environment mapping  4  3.1  The rendering equation illustrated  13  3.2  The transformation method: sampling via an inverted C D F .  21  3.3  Motivation for the alias method  23  3.4  The alias method  24  3.5  Sampling a 2D intensity image  27  3.6  A quantized intensity map  29  3.7  Two interpretations for the BRDF  32  3.8  The Phong BRDF  33  3.9  Sampling from the cosine term/Phong BRDF product . . . .  34  4.1  Angular maps of the E M , BRDF and product distribution . .  40  4.2  Angular maps contrasting sampling methods  41  4.3  Rejection sampling illustrated  43  4.4  Rejection sampling the product distribution  44  4.5  Our SIR method for sampling from the product distribution .  49  5.1  David in Grace Cathedral. Phong exp 10, k = 1.0, kd = 0.0 .  53  5.2  David in Grace Cathedral. Phong exp 50, k — 1.0, kd = 0.0 .  54  5.3  David in Grace Cathedral. Phong exp 50, k — 0.5, kd = 0.5 .  55  s  s  s  List of Figures  vi  5.4  David in St. Peters. Phong exp 50, k — 0.5,  . . . .  57  5.5  David in Grace Cathedral. Diffuse Phong (k = 0.0, kd = 1.0)  58  5.6  David in Grace Cathedral, unzbomed  59  5.7  Comparison with Veach and Guibas(l)  60  5.8  Comparison with Veach and Guibas(2)  61  5.9  Comparison between SIR and converged image  62  s  s  = 0.5  Acknowledgements  vii  Acknowledgements Firstly, I'd like to thank everyone who has contributed ideas, discussion and proof-reading to this work; in paticular, Wolfgang Heidrich, Abhijeet Ghosh, Hendrik Kueck, Simon Clavet, Nando de Freitas, and James Slack. I also want to thank my loved ones for putting up with me: my mother Elizabeth, my father Wayne, my sister Rebecca, and my dear Amy. Also, a big shout-out to Imager Small and everyone else who've kept me sane via games of Quake 3, UT2004, Panel de Pon, Bomberman, Netris, and general lab craziness. In no paticular order: Vlad "Vlady" Kraevoy, Ben "Tron" Forsyth, Ritchie "spankbot" Argue, Hendrik "Space Cowboy" Kueck, Eric "*nipple*" Brochu, James "evil" Slack, Fred "Cheerleader" Kimberly, David "Toasty" Pritchard, Simon "Frog" Clavet, Kristian "zirkusdirektor" Hildebrand, Ahbijeet Ghosh, Roger "Butt Slammer" Tarn, Lisa Streit, Matt Trentacoste, Lewis "Wurmangst" Johnson, and the ever-snipey "Grumpy old man" (you know who you are). Lastly, I'd like to thank der kommuners, past and present, for their friendship and support: Eric Brochu, Eddy Boxerman, Hendrik Kueck, David Pritchard, and Tyson Brochu.  Chapter 1. Introduction  1  Chapter 1  Introduction Image-based representations for illumination, such as environment maps (EMs) and light fields, have received considerable attention in recent years. The primary reason for this attention is that images can capture complex real-world illumination which is difficult to represent in other forms. When integrating image-based lighting, such as environment maps, into a rendering system, the employment of a good sampling strategy for illumination is an important issue. While several researchers have recently worked on this problem, the basic approach that has been taken by most work is an importance sampling strategy based on the energy distribution in the environment map. Unfortunately, such an approach performs poorly for highly specular surfaces, since samples chosen this way have a low probability of lying within the specular lobe. Similarly, if importance sampling is performed based solely on the reflectance function (BRDF) of the surface, then the sampling will not perform well for high-frequency environment maps. This thesis introduces bidirectional importance sampling, a method that samples visibility according to an importance derived from the product of BRDF and environment map illumination. The challenge of this approach is to develop an efficient means of drawing samples from this product distribution. The task is complicated by the fact that the 2D BRDF slice varies from point to point on the surface. Furthermore, the environment map is  Chapter 1.  Introduction  2  usually represented relative to a global coordinate frame, while the BRDF is expressed in a local frame that changes with surface orientation. For these reasons, precomputation approaches, such as those that employ a table of the product distribution, are infeasible. We present two solutions to this problem. The first is a combined rejection sampling and importance sampling scheme that initially generates samples according to the product distribution with rejection sampling, and then uses the resultant samples to estimate local illumination with standard importance sampling. The second approach uses a technique called sampling-importance resampling. In this method, samples are generated from either the lights or the BRDF and then resampled to distribute them according to the product distribution. While both methods increase the cost of sample generation, we demonstrate significant quality improvements for the same compute time under the assumption of BRDF representations that support efficient evaluation and sampling. Our method creates samples on the fly and does not require expensive precomputation. The remainder of the thesis is structured in the following manner. In Chapter 2 we review some of the relevant work in environment map rendering. In Chapter 3 we give an introduction to sampling and Monte Carlo rendering, before describing our approach in Chapter 4. Chapter 5 discusses the results of our methods, and we conclude in Chapter 6. Part of this work has appeared in a different form. Our rejection sampling technique appeared as a Siggraph technical sketch, published with Abhijeet Ghosh and Wolfgang Heidrich [4].  Chapter 2. Related Work  3  Chapter 2  Related Work All rendering systems, both global and local, must at some point compute the direct illumination in the scene; that is, the amount of light arriving at surfaces from primary light sources.  Unfortunately, this task remains  expensive; as such, much research effort has focused on the development of more efficient techniques for computing the direct lighting. Little work deals directly with sampling from environment maps for global illumination applications.  Debevec [6] describes the use of high dynamic range envi-  ronment maps in global illumination in the context of Ward's RADIANCE package [47]; however, no specialized sampling strategies are employed to reduce variance for this specific light representation. This thesis focuses on direct lighting from environment maps. In this chapter, we present an overview of previous work on computing direct illumination from environment maps. After an introduction to environment mapping in Section 2.1, we address the approach of environment map prefiltering in Section 2.2. Section 2.3 summarizes work in alternate bases and precomputed light transport, and Section 2.4 discusses point set relaxation techniques.  Chapter  2.1  2.  Related  4  Work  Environment Mapping  Environment  mapping  was first proposed by B l i n n and Newell [3].  Envi-  ronment m a p p i n g operates under the assumption of a distant environment: when the distance between an object and its environment is large compared to the object's size, i n c o m i n g i l l u m i n a t i o n from the environment can be considered directional. T h a t is, the light sources i n the environment can be thought of as being located at infinity - for example, t h i n k of the d a y t i m e sky lighting the surface of the earth. U n d e r this assumption, the i n c o m i n g i l l u m i n a t i o n at the object can be efficiently stored i n a 2 D image called an environment  map ( E M ) . D u r i n g rendering, the i l l u m i n a t i o n can be indexed  by, for example, the spherical coordinates of the reflected ray, as i l l u s t r a t e d in F i g u r e 2.1.  Figure 2 . 1 : Environment mapping. Incoming light is assumed directional, and can thus be sampled and stored in an image. Rendering involves computing the amount of light reflected towards the viewer (direction w ) ; this requires r  considering light incoming from all directions (w,'s).  Chapter 2. Related Work  2.2  5  Prefiltered Environment Mapping  When computing shading for a non-mirrorlike object, it is generally necessary to consider multiple E M light directions. To mitigate this, the notion of prefiltering  the environment maps was suggested, first by Greene [12]. The  idea is for each E M entry to represent not just the incoming radiance from a single direction, but also the radiance from a larger region that has been integrated against the BRDF lobe of the surface. With this approach, one need only perform a single E M look-up to compute approximate shading for the surface patch. In the case of Greene, the E M was prefiltered with the surface's diffuse reflectance. Heidrich and Seidel [15] extended the approach to glossy BRDFs by prefiltering specular environment maps with radially symmetric Phong lobes. Kautz and McCool [20] addressed other isotropic BRDFs by developing a representation for more complex BRDFs that can still be used to prefilter environment maps. A drawback of these approaches is the precomputation required. This drawback was partially reduced by the hierarchical prefiltering algorithm of Kautz et al. [21] which significantly accelerates the preprocessing time. However, the situation remains difficult when the scene contains multiple BRDFs or a surface with a spatially-varying BRDF. In these cases, the environment must be prefiltered separately against each BRDF model, resulting in large computational and storage costs. Also, efficient representations do not exist for many types of BRDFs, particularly anisotropic models or measured BRDFs. Regardless of whether or not the environment map has been prefiltered, when rendering with a single E M lookup, one does not take visibility or occlusion into account, either during preprocessing or rendering. Thus, one  Chapter 2. Related Work  6  cannot produce images that contain such effects as shadows and multiple scattering.  2.3  Alternate Bases and Precomputed Transport  In some recent work, the illumination and/or BRDF are projected into finite bases, such as spherical harmonics (SH) and wavelets. These finite bases support efficient compression, as only a small number of terms are necessary for representing the source functions to a reasonable fidelity. For example, Ramamoorthi and Hanrahan [33] express the lighting in a spherical harmonics basis and then use only the first nine coefficients to render diffuse objects. However, because ringing artifacts occur when highfrequency lighting is expressed in the SH basis, they must heavily downsample and blur their environment maps before applying their technique. This also implies their technique is unusable for high dynamic range environments. They extend their work in [34] from Lambertian BRDFs to more complex isotropic BRDFs, Sloan et al. [39] introduces precomputed radiance transfer, where visibility and self-transfer effects are taken into account during an offline process. The precomputed radiance is then expressed in spherical harmonics so that it can be compressed and used for interactive rendering. Like Ramamoorthi, their technique is limited to very low frequency lighting. Ng et al. [29] and Liu et al. [27] factorize with wavelet representations which better preserve frequency detail in the environment, albeit at resolutions that are considerably lower than our technique can handle. Ng et al. [30] presents a formulation for factorizing the rendering equation into separate terms, each of which is expressed in a wavelet basis. Their rendering times are non-interactive.  Chapter 2. Related Work  7  Conceptually, what these techniques do is essentially render the scene from all viewpoints and then store a compressed version of these images for access at render time. As such, all flavours of factorization and precomputed transport require hours of preprocessing and huge amounts of storage, often in the gigabytes. What is more, the preprocessing must be repeated every time new geometry or material properties are introduced. Such massive precomputation is prohibitive in most every setting, especially in animation when geometry is constantly changing. Our approach can handle high resolution all-frequency lighting and requires virtually no precomputation. Geometry, environment maps, and surface materials can be substituted on the fly at essentially no cost. While our rendering times are longer, at perhaps twenty seconds instead of two seconds, our algorithm comes with the flexibility of requiring no precomputation.  2.4  Importance Sampling and Point Relaxation  A common approach for rendering from environment maps is a technique called importance sampling, where the full illumination integral is approximated by considering only a small set of sample directions. These directions are chosen according to their contribution to the reflected radiance. One way to achieve this is to approximate the E M with a small number of point lights representing bright regions in the environment. For example, Kollig and Keller [23] propose a scheme that is based on relaxing a point set until its distribution fits the energy distribution of the environment map. Relaxation schemes attempt to serve as a compromise between light intensity and reasonable spatial coverage. Kollig's algorithm is similar to earlier work by Deussen et al. [8] on distributing stipples for digital halftoning. There is  Chapter 2. Related Work  8  also LightGen [5], which seeks to reduce the diffuse illumination from high dynamic range environments into a set of directional light sources using a similar relaxation technique. The problem with point relaxation methods is that they are not proven to converge in 2D, and they are often time consuming for a sufficiently large number of points. As a consequence, techniques that use relaxation precompute only a single relaxed point set, and then use this set to shade all points in the scene. Ostromoukhov et al. [31] presented a technique for distributing 2D point samples which is much faster than relaxation-based approaches, and also appears to produce a good spatial distribution for the points. Agarwal et al. [1] introduced a sampling method for environment maps where the sampled importance takes into account both the energy in the environment map and the solid angle separating the samples. This way, close clustering of environment map samples is avoided, which reduces redundant shadow tests. Like Argarwal et al., we have included the solid angle in the importance term for the lights; see Section 3.4. To distribute their samples, Agarwal et al. use a point relaxation algorithm which is different from that of Kollig and Keller; however, a downside of their approach is that the environment map has to be quantized to generate the distribution. As an extension to their work, Agarwal et al. sort the samples for each shading operation by the magnitude of their contribution to the final illumination. They sample all point lights deterministically, in order of contribution, until the contrast that the remaining lights can add falls below some threshold. This use of the product of BRDF and environment map value is a step towards our approach of drawing samples according to an  Chapter 2. Related Work  9  importance that is the product of BRDF and light distribution. However, because they precompute a point set via relaxation, they are limited to choosing directions only from a set of precomputed samples, thereby introducing quantization artifacts (aliasing and banding). Also, the sorting introduces bias. In contrast, our methods are asymptotically unbiased, and can create samples efficiently on the fly, so that different sampling patterns can be used throughout the scene. Secord et al. [36] describes a fast algorithm for distributing stipples according to image intensities based on pre-integrating and inverting the probability density function (PDF) derived from image intensities. This same approach applies to drawing samples efficiently from environment maps, which are themselves images. We use a variant of this method in our current work; see Section 3.3 for details.  2.5  Combined Importance Sampling  Also similar in spirit to our approach is the work by Veach and Guibas [44] in which they combine sampling from the light sources and sampling from the BRDF to reduce the variance of the results. Prior to rendering, a decision is made as to how many samples to draw from each distribution. The general strategy is to draw more samples from the distribution with higher variance. Unfortunately, employing combined sampling results only in a blend between the variances of the individual distributions (see Section 4.1). Our method goes further than their work to reduce variance in that we sample  directly  from the product distribution, rather than just mixing samples taken from the individual distributions. The work of Szecsi et al. [42] selects sample weights with correlated  Chapter 2. Related Work  10  sampling, a technique that seeks to separate Monte-Carlo estimators into low and high variance components and combine samples from the components to reduce the overall variance. Unfortunately, their results are not convincing. Furthermore, as with Veach and Guibas, decisions must still be made a priori as to how samples are combined. Our approach is mathematically straightforward and allows for direct sampling of the product distribution without any arbitrary guesswork.  Chapter 3. Monte Carlo Rendering  11  Chapter 3  Monte Carlo Rendering The general goal of global illumination is the computation of all light interactions in a scene to generate realistic images.  Many global illumina-  tion algorithms, such as radiosity [11], irradiance gradients [47], and photon mapping [16] consist of two phases: first distributing radiance from the light sources onto surfaces, and then gathering the local radiances to construct an image for a specific viewpoint. A thorough distribution step will permit for a fast gather step, one which can produce quality images in a handful seconds. To date, computing the direct illumination remains the expensive aspect of global illumination; as such, much research effort has focused on coming up with efficient techniques for direct lighting. The work described in this thesis is an approach for tackling the direct illumination problem in the case where both the light sources and the materials are complex - that is, of arbitrarily high frequency. This chapter builds a foundation for the presentation of our methods. Section 3.1 introduces the mathematics of direct illumination rendering from an environment map. In Section 3.4 extends the method for light sampling to include area considerations. Section 3.2 presents some basic notions of sampling and Monte Carlo integration. From these, we construct methods in Sections 3.3 and 3.5 for drawing samples according to the environment map and BRDF.  12  Chapter 3. Monte Carlo Rendering  3.1  Environment Mapping  While following the discussion in this section, refer to Figure 3.1 for a diagram of the geometry involved. The convention is to denote directions with co. The subscripts i and r denote quantities that are incoming and reflected (outgoing), respectively; for example, u>i and to . T  The integral for direct illumination is given by the following expression:  .  L (x,u ) r  r  = / f (x,Ui-+v )Li(ui)(n(x) r  • Ui)V(ui)dui.  r  (3.1)  Equation 3.1 describes the amount of radiance (light) being reflected at a surface point x towards some outgoing direction cu . This reflected radiance r  is computed by integrating the incoming radiance Li{uJi) over the hemisphere centered at x of incoming directions u>i. The reflectance properties of the surface are described with the bidirectional reflectance  distribution  function (BRDF), denoted f {x,u>i —* u> ). The BRDF is a function indiT  r  cating how much of the light incident to x along the direction w; is reflected away towards to . The n(x) • u>i term is known as the cosine term] it is an r  area foreshortening term that scales the incoming radiance Li{u>i) based on orientation of the surface normal n(x) with respect to the incoming light direction LOi. V(ui) is the visibility term, which in our case is binary; that is, V(tUi) is 1 if the surface point x can see the light along Wj and 0 otherwise. In this work, we express the incoming illumination  with an environ-  ment map (EM), a discrete representation of the continuous hemisphere of lights surrounding x. The environment map in our case is a latitudelongitude map - that is, a 2D image that is projected onto the hemisphere. Recall that environment mapping operates under the assumption of distant illumination; that is, each pixel in the E M is treated as a point light source  Chapter  3. Monte Carlo  13  Rendering  Figure 3.1: Shading for surface point x is computed by considering the light from all incoming directions Ui of the hemisphere f i . Each light L(u>i) is weighted by the B R D F (circled i n red), describing the probability that the surface at x reflects that light towards the viewing direction w . T  corresponding to the incoming light summed over a differential patch on the hemisphere. A n E M pixel can be uniquely indexed with a direction u>i expressed in spherical coordinates. Under illumination from an environment map, we can rewrite the direct illumination integral as an explicit summation over the E M entries: JV  L (x,u) ) r  T  =  ^2f (x,Ui(j) r  - » uj )Li{iOi{j)){n{x) r  • Witf))V{m(j))S(w(j))  f  (3.2) where the incoming light direction u>i is now a function of the E M pixel index j. The additional weighting term S(j) represents the solid angle subtended by the spherical patch corresponding to element j. Computing the direct illumination integral involves the evaluation of the  Chapter 3. Monte Carlo Rendering  14  visibility term V(oij). Determining whether or not the surface point x can see the light at u>i typically involves a ray cast. A ray r with origin x and direction tOi is intersected with the scene geometry. If no intersection occurs, x can see the light Li(u)i), and so V(o>j) is set to 1. Otherwise, some part of the scene is occluding x's line of sight to the distant environment along u>i. V(ui) is thus set to 0, negating the light Lj(u;i)'s contribution to the reflected radiance. Performing the ray-scene intersection test makes up the bulk of the work during rendering. Even with the aid of hierarchical acceleration structures, visibility queries dominate the rendering time. In our implementation, for example, typical timings for ray intersection queries with scenes of modest complexity - half a million triangles - are on the order of 4-5 microseconds on our test machine. If one was to precisely evaluate the reflected radiance as indicated in Equation 3.2 by summing over the entire E M , then for sufficiently interesting environment map resolutions (1024x512 pixels), one would have to wait for a couple of seconds to light a single pixel in the output image. A small output image of size 128x128 would take 12 hours to render; this is obviously far too slow. Intuitively, one realizes that it is not necessary to sum over the entire environment map when computing direct illumination, because there are certain regions of the environment that contribute far less light than others. As an example, think of a night scene lit by a full moon in a clear sky. A vast majority of the light received by the scene will be coming only from the small set of directions corresponding to the location of the moon in the sky. In fact, for most real-world lighting conditions, the contrast ratio between the lightest and darkest points in the environment is very high, typically five to eight orders of magnitude [7].  Chapter 3. Monte Carlo Rendering  15  This leads, us to the notion of Monte Carlo ray tracing. When rendering from environment maps, one need only consider a small set of incoming light directions. If these directions are carefully chosen in such a way that their contributions make up a significant amount of the reflected radiance, then one can expect to have computed a decent approximation to the illumination integral. The challenge is to decide how to distribute this small set of sample rays effectively. Possible approaches are to distribute rays according to the intensity of the lights, as suggested above, or according to the reflectance of the surface. In practice, we do this randomly - hence the name Monte Carlo. These topics are discussed in detail in Sections 3.3 and 3.5. First, we require valid mathematics for approximating integrals; this is the subject of the next section.  3.2  Monte Carlo Estimators and Sampling  We present here an overview of Monte Carlo (MC) integration and basic sampling theory.  Many introductory treatises of probability theory and  Monte Carlo methods exist in the literature. The reader is pointed to Kalos and Whitlock [19] and Hammersley and Handscomb [14] for basic texts on probability theory and Monte Carlo methods, including Monte Carlo integration, rejection sampling and importance sampling. For introductions to Monte Carlo methods in ray tracing and rendering, see [9, 17, 18, 37].  3.2.1  Monte Carlo Integration  The goal in rendering is to compute the reflected radiance for every surface point visible to the viewer. In our case, this amounts to evaluating the direct  1  16  Chapter 3. Monte Carlo Rendering  illumination integral for every pixel in the resulting image. As described in the previous section, this is too computationally intense to solve directly. Instead, we employ Monte Carlo integration, a stochastic technique that allows us -to compute approximations for integrals such as the illumination integral for which an exact evaluation is intractable. Suppose we wish to evaluate the integral  Hf)= [ f(x)p(x)dx, •Ix  (3.3)  where the space x is the domain of definition of the function. Of course, /(/) need not necessarily be a continuous integral. As is the case molecular dynamics, for example, /(/) could be a very large sum and x the perhaps uncountable set of possible configurations of the system. In any case, when x is high-dimensional, /(/) becomes infeasible to compute with deterministic integration methods. One can approximate such integrals or large sums /(/) with tractable sums.ijv^/). Given a set of identically independently distributed samples X = {xi, X2, • • •, XN} drawn from some density p(x) defined over the space X, the tractable sum 1  '  N  W) =iv^  / ( x i )  i=l  (3  -  4)  converges asymptotically as N —> oo to  Hf) = / f(x)p{x)dx. x  Let us evaluate the expectation value of our estimator IAT(/)- The expectation of a function f(x) is defined by  (/(*)) = J f(x)p(x)dx.  (3.5)  17  Chapter 3. Monte Carlo Rendering for a continuous variable and (/(z)) =  £/(s)p(s) X  for a discrete variable. Among the properties satisfied by expectation values are the following: (aX)  =  a(X);  i  (3.6)  i  where X is the random quantity and a is any constant. We wish to evaluate the expectation of our estimator i j \ r ( / ) from Equation 3.4: 1  vW)>  N  = <jy £ / ( * * ) > •  (3-8)  Applying Equations 3.6 and 3.7 to the right side of Equation 3.8 yields 1  N  (W)> = ^E</^)>i=l  By the definition of expectation given in Equation 3.5, this can be rewritten as  (i*(f)) = ^ NE / i=l 1  N  = ]vE'(/); =  /(/).  The expectation of 7JV(/) is therefore /(/). IN(I) is an unbiased estimate of /(/), and by the strong law of large numbers, will almost surely converge to the target integral /(/) [2].  Chapter 3. Monte Carlo Rendering  18  The strength of Monte Carlo integration over deterministic integration schemes is essentially that M C integration seeks to evaluate f(x) only in those areas where probability is high - that is, where the contributions of f(x) to the integral /(/) are most significant. Let us now examine the variance of  The variance of a random  quantity X, written var(X), is the expected squared difference between a realization of the variable and its expectation. That is, var(X)  EE =  ((X-(X)f); (X )  - (X) .  2  2  •  For independent random variables X{, variance observes the following properties: var  (E <) x  var(aX)  = E ( ); var  Xi  = a var(X). 2  With these relations in mind, we see that var(7iv(/))  =  v a r  1  ^E/( *)^; x  N  i=l  =  ^var (/(*)).  (3.9)  In rendering, variance manifests itself as per-pixel noise in the output image. Equation 3.9 shows that the variance in a Monte Carlo estimator of an integral /(/) is inversely proportional to the sample size N.  This  result tells us that we can reduce noise in the output image by increasing the number of sample directions we use to approximate the illumination integral.  Chapter 3. Monte Carlo Rendering  19  However, the error in our estimate actually behaves similarly to the standard deviation of the samples; that is, the square root of the variance. Thus is demonstrated the fundamental problem with Monte Carlo integration: the notion of diminishing returns. Because image quality depends on / V , one 2  must quadruple N to halve the error.  3.2.2  Sampling a General Distribution  So far, we have not mentioned how to choose the density function p(x) for distributing samples during Monte Carlo integration. Let us defer this discussion for a time, assuming for the moment that we know a good distribution p(x) to sample from. This section presents two methods for redistributing a uniformly distributed variate according to some arbitrary probability density p(x).  Cumulative Distribution Functions and the Transformation Method A common technique for sampling from general PDFs is known as the transformation method, which requires the notion of the cumulative distribution function. The cumulative distribution function (CDF) is an alternate way of completely describing the distribution of a real-valued random variable X ~ p(x). For every real number x, the CDF is given by C(x)=p{X  < x).  In other words, C(x) is the probability that the random variable X takes on a value less than or equal to x.  f  Chapter 3. Monte Carlo Rendering  20  The CDF can be computed with C(x) — I  p(x)dx,  Jxo  where xo is the smallest value that the variable X may take on, generally —oo. If xo is finite, the random variable can always be translated so that x = 0. 0  We can now use the transformation method to sample from p(x). Specifically, we choose uniform deviates Ui ~ Z^(o,i) and transform them according to  Xj  =  C  _  1  (tii).  Refer to Figure 3.2 for an illustration. Note that since  p(x) > 0, C(x) is monotone increasing, and so C ( x ) exists everywhere. _1  Discrete Distributions We now examine the case where X is a discrete random variable. Here, the probability distribution p(x) describing X is represented by a table of values p[i] — Pr[X = i], where i = 1,... ,N runs through the set of all possible values of X. Such a PDF representation can be created from any discretely sampled function f[i] in 0(n) time via the expression:  P\i] =  ^ N \ l  u  V  Ei=i/W We can evaluate a discrete CDF C[i] according to N  = £>[<]• i = i  Now, to sample from p(x), we use the transformation method as before. That is, choose u ~ W(o,i)  a n  ^ find the table index i that satisfies C[i] < u >  C[i + 1]. This requires a binary search through C[i], and so sampling this way takes O(logiV) time per generated sample. Practically speaking, one  21  Chapter 3. Monte Carlo Rendering  can often do fairly well to use i = [uN\ as a starting point for the search. A 0(1) sampling can be achieved by pre-inverting C,-which takes 0(N) time. Depending on the scenario, this may not be worth the effort. For example, if C is constantly changing, or if only a small number of samples will be drawn, sampling directly from C in O (log TV) time may well be faster than pre-inverting C and sampling in 0(1) time.  Figure 3.2: Sampling via C D F inversion. Left: a I D P D F , P(x). Middle: the corresponding C D F , C(x).  Right: uniform samples transformed by C~ (x). 1  Note along the x-axis how samples have been redistributed according to P(x). Image courtesy of [36].  Faster Sampling: the Alias M e t h o d Another method for sampling from discrete distributions is the alias method, which was proposed by Walker [46]. This ingeniously simple technique computes samples in 0(1) time, but with an initialization that is less costly than the C D F inversion and pre-normalization of the transformation method. Strangely, the alias method has been largely overlooked by the graphics community, where the transformation method seems to be the standard technique for sampling from discrete distributions. As such, this section  Chapter 3. Monte Carlo Rendering  22  presents the idea behind the method of aliases, and the variant of Vose [45] that we use in this work. Let x be a discrete random variable distributed over the set S = {0,1, 2,3} with corresponding probabilities P = {0.1,0.2,0.4,0.3} (see Figure 3.3, left). We want to generate sample values for x according to the distribution P. Marsaglia [28] introduced the following simple approach. Assuming we have a uniform random number generator, one way to sample from P is to generate a uniform random integer i in [1,10] and use it to index into the set L = {0,1,1,2,2, 2, 2, 3, 3, 3} (see Figure 3.3, right). L has been.chosen such that the number of times each element of S appears in L is proportional to the probability of that element occurring, as described by P. Note that the size of L depends on the probabilities in P.  If, for example, x was  distributed over {0,1} with probabilities {1/2953,1 - 1/2953}, this method would require L to have 2954 elements. This simple example distribution can obviously be sampled more efficiently; however, for arbitrarily complex discrete distributions, the required size of L is bounded only by the precision of the random number generator - generally several orders of magnitude. Marsaglia's approach was how arbitrary discrete distributions were sampled prior to Walker's method of aliases, which introduced storage costs that are only linear in the input size while maintaining constant-time sampling. We now describe the sampling approach behind the alias method. Let S = {0,... ,N — 1} be the domain of distribution of x, as before, with corresponding probabilities P — {po,... ,p;v-i}. The method operates on two lists, prob and alias, to generate an integer x from S as follows: given a real-valued uniform variate u ~ [0, TV), let j = [uj. If (u — j) < probj, then choose x = j; otherwise, choose x = alias j. The particular technique we use to create the lists prob and alias is  Chapter  0 Figure 3.3:  1  2  3.  Monte Carlo  3  S.  Rendering  23  i  Marsaglia's sampling method [28]. Left: the target distribution  P. Right: Sampling is achieved by redistributing a uniform random integer i in [1,10] according to L.  that of Vose [45]. Consider again the d i s t r i b u t i o n S = { 0 , 1 , 2 , 3 } ~ P = {0.1,0.2,0.4,0.3}, as illustrated on the left side of F i g u r e 3.4. Vose's algor i t h m can be visualized b y drawing a horizontal line i n the plot of P at p r o b a b i l i t y 4 . Elements of P that are over this line are "chopped off" and their contributions redistributed to elements w h i c h have a smaller p r o b a b i l ity, creating tables prob and alias.  O n the right side of F i g u r e 3.4, the x-axis  corresponds to table indices. T h e heights of the shaded regions correspond to the real-valued entries of prob, while the numbers i n the lighter regions correspond to the indices stored i n alias. <=  - L , probj = NPj.  appears i n alias.  For those elements Pj that are  For those elements that are > jj., probj = 1 and j  S a m p l i n g amounts to p i c k i n g a r a n d o m index j and deter-  m i n i n g if it falls w i t h i n probj. If so, return j, otherwise r e t u r n the reference alias j to a value whose p r o b a b i l i t y overflowed into b i n j .  Chapter  Figure 3.4:  3. Monte Carlo Rendering  Walker's alias method [46].  and the cutoff line at  Left:  24  the target distribution P  = 0.25. Right: Values with large probabilities have  references distributed to bins where probabilities are small.  Refer to A p p e n d i x A for a discussion and i m p l e m e n t a t i o n of the alias m e t h o d , i n c l u d i n g Vose's linear-time a l g o r i t h m for i n i t i a l i z i n g the lists alias  and probs.  3.3  Sampling the Environment Map  In situations where we need to draw samples according to the i m p o r t a n c e of the environment m a p image, we employ an approach similar to Secord et al. [36] for s a m p l i n g from a 2 D P D F . For an R G B image, our n o t i o n of i m p o r t a n c e is the intensity of each pixel - that is, the s u m of the intensities of each color channel. W e begin w i t h an overview of the m e t h o d , followed b y a detailed description.  25  Chapter 3. Monte Carlo Rendering O v e r v i e w o f 2D S a m p l i n g Preprocessing Given an NxM RGB image, we perform the following precomputation:  1. For each pixel, compute an intensity J(x, y) = red(x, y)+green(x, y) + blue(x, y) 2. For each of the y = [1,..., M] scanlines, create a PDF based on the intensities of the pixels in that scanline.  At the same time, com-  pute the average intensity of the pixels in that scanline: I g(y) aV  =  3. Finally, create a PDF from the distribution of M scanline averages lavg  }  .  Sampling One can now generate samples according to the image intensities: 1. Generate U j ~ ^(0,1)2. 2. Redistribute U i a scanline  X i  i  V  t  V  according to the PDF of scanline averages to select  .  3. Redistribute u^ according to the PDF of the pixel intensities in scanx  line  X i  t  V  ,  to select pixel location x^ . x  2D S a m p l i n g i n D e t a i l Our goal is to redistribute a set of uniform 2D samples U = {u\,... where Ui ~ ZY(o,i) according to a PDF p(x), generating a set X = {x\,..., 2  ,u } n  x} n  26  Chapter 3. Monte Carlo Rendering with Xi ~ p(x). Let x; X{.  To determine  X i  >  y  ,  x  and x^  y  denote the x and y coordinate of sample  we compute the cumulative density function  C(y) = rm(t)dt,  (3.10)  Jo  where m(y) is the marginal density function of p{x). In our case, we have a 2D image, and so m(y) can be thought of as the average intensity of scanline y in the environment map. Note that M(y) is non-zero for those scanlines containing non-zero intensity, and so using the transformation method we obtain  Xi  >  y  =  C (UJ - 1  J 2 /  ).  Given X J , the uniform sample U i iV  t  X  PDF of the respective scanline, p (x}. Xiy  is now redistributed according to the This can be accomplished via the  conditional probability distribution  c ( x | x  -  }  =  ^K^y  ( 3  -  n )  and the corresponding cumulative distribution C(x\x ) ity  = / c(s\x^ )ds Jo y  As before, the x coordinate of the new point can be found through C  (Ui  >x  ( -12) 3  Xi  y  X  =  1^1,2/)  Our 2D environment map is a discrete PDF, and so we use either the transformation method or the alias method for sampling in O(l) time. The precomputation of scanline PDFs and data structures have linear requirements in both space and time. The precomputation need only be done once per environment map. Also, the preprocessing is fast, at only milliseconds even for large images; it can thus performed on-the-fly as the image is loaded. Importance sampling from the E M is now just a simple table lookup, where the table indices are uniformly random points. Practically speaking,  Chapter 3. Monte Carlo Rendering  27  Figure 3.5: Sampling a 2D intensity image. Two I D passes are performed: a scanline is selected according to the distribution of scanline averages, and a pixel is then sampled from within that scanline.  rather than transforming a uniform distribution, one can use points taken from a low-discrepancy sequence such as the Halton sequence [36] or a bluenoise distribution [31]. These choices will result in a sample set which is more evenly distributed in space, with fewer clusters of close points. For further information, the reader is referred to Keller's thesis [22] for an excellent presentation of quasi Monte Carlo methods and their use in rendering.  3.4  Improving Shadows: Area-weighted Lighting  In principle, there is always a practical lower-bound on the number of visibility tests one can perform per pixel and still achieve a low-variance image. For example, consider a high-frequency environment map that is dark everywhere except for i V small, bright area light sources. For a surface patch  Chapter 3. Monte Carlo Rendering  28  which can see all N lights, the best thing to do would be to draw a single sample from each of the light sources. If ever we render with less than N samples, there will always be some amount of variance in the resulting image, because visibility will not have been tested for each light. This noise is particularly noticeable in penumbra regions (soft shadow boundaries), as each pixel is seeing a different random sampling of the lights. In practice, this is often acceptable; as will be shown in Chapter 5, our approaches for sampling from the product distribution still give us nice images even for an extremely small number of samples. Sampling purely from the E M or BRDF takes many more samples (hundreds), and thus implicitly generates shadow boundaries which are generally good. Up to this point, sampling from the environment map has meant selecting directions (pixels) based solely on pixel intensity - that is, based on light brightness. In practice, however, it is desirable to distribute light samples with a good spatial distribution as well as by intensity. The motivation for including an area factor is to address the oversampling of small, bright light sources.  Once a small light source has been  sampled from, there is little additional gain in drawing another sample from this light, because it is behaving essentially as a point source. By drawing a single sample from it, we already have an idea of the intensity and color of the light coming from this source. Furthermore, our shadow quality will not improve or smoothen, because the spatial extent of the light is small, and the visibility function is not likely to vary across the light. As an example, consider an environment map containing several small light sources with equal area, with one light being significantly brighter than the others. If directions are chosen based solely on E M pixel intensity, the majority of the samples will come from the bright light. However, because  Chapter 3. Monte Carlo Rendering  29  this light source is small, only one or two samples from this light would be enough. We're better off casting a few samples to some of the dimmer lights in order to smooth out shadow boundaries. To this end, we take the approach of Agarwal et al. [1] and modulate the light intensity by an area term based on the solid angle of the light. The basic idea is to first identify the light sources in the environment map. Then, the importance of all pixels making up a light source are scaled based on the area of that light. In effect, we want to penalize small light sources more than large ones, so that small sources are not oversampled.  Figure 3.6: The Grace Cathedral environment after quantization based on the logarithm of pixel intensity. Observe how area light sources in the environment, such as the altar and windows, have been localized.  Our light classification scheme proceeds in the following manner. We first compute the histogram of the image and quantize it into N intensity levels. The binning is performed on the logarithm of pixel intensity. Because we're interested in high-frequency H D R environment maps, where the contrast  Chapter 3. Monte Carlo Rendering  30  ratio is 5-7 orders of magnitude, the majority of E M pixels are low intensity. So, this classification is performed logarithmically rather than linearly. Next, connected components are found in the quantized image by running a breadth-first search. The area (solid angle) of each connected component is found by summing the areas (solid angles) of each pixel making up the light source. Importances L; of the pixels, originally taken just from intensity, are now scaled by this solid angle area. The new importance is given by L oc LAco. Specifically, we employ the conclusions of Agarwal et al. [1] and select L = L (min(0.01, Aco)) , where b is in the range [0.1,0.2], b  depending on the average size of the light sources in the environment map. The time complexity of the area-weighting algorithm is linear in the size of the E M , and is sufficiently fast so that it can be performed in negligible time during loading. In practice, the bright areas in measured environments are often reasonably well distributed across the sky, and so this factor does not significantly affect the quality of the sample distribution. Nevertheless, the results presented in Chapter 5 were generated using this augmented light importance.  3.5  Sampling from the BRDF  Another method for choosing sample directions for visibility testing is to distribute rays according to surface reflectances (BRDFs). This is particularly helpful in scenes that have, for example, rather uniform light distributions (think outdoor scenes). Sampling based on lighting intensity in such a scene does not significantly localize important directions. However, sampling from glossy or shiny materials in the scene will do well because of their non-uniform reflection properties. Examples of such surfaces include  31  Chapter 3. Monte Carlo Rendering  metals, plastics, anything with a finish or polish, or surfaces such as asphalt or water when viewed at grazing angles. Recall that the BRDF f (x,u>i —> ui ) states what fraction of the power r  r  arriving at a surface point x from an incoming direction u>i will be reflected along an outgoing direction u> (see Figure 3.7). We assume that the BRDF r  is energy preserving; that is, x, ix>i  — >  uj )duj = 1. r  r  (3.13)  In other words, we assume that the BRDF reflects all energy incident upon it, and that no energy is generated or absorbed. While not strictly necessary for computing single-bounce direct illumination, this restriction allows for an alternative interpretation of the BRDF: the BRDF f can be thought of r  as a probability density that describes the probability that an incoming ray of light will be randomly scattered towards a particular outgoing direction. Under this view, each incoming ray maps to a single outgoing ray with equal intensity, but directed along directions sampled according to the magnitude of the BRDF. Furthermore, the BRDF is reciprocal, and so incoming and outgoing directions can be exchanged. Hence, another interpretation is to view the BRDF as a function which describes which incoming light directions uii will contribute to the light reflecting along co . r  There has been much previous work in efficient sampling from BRDFs (for example,  [24, 25, 40] ). There are several analytical representations  for the BRDF that support sampling directly. In particular, we use the energy-preserving (normalized) form [26] of the Phong model [32]: f (x,Ui r  -^u> ) = - - ( n + l)cos"6'-| i  r  6 is the angle between the view direction u  T  .  (3.14)  and the reflected light ray  R, where R is the incoming light direction u>i reflected across the surface  Chapter  3. Monte Carlo Rendering  32  n(x)  (Oi's  V  (b)  (a)  Figure 3.7: Two interpretations for the B R D F . The B R D F describes: (a) the distribution of scattering directions uj for light which is incoming along <*;»•; or, r  (b) the distribution of incoming light directions uii that send light towards u> . r  n o r m a l n ( x ) . See F i g u r e 3.8 for a n i l l u s t r a t i o n . k (n + l ) c o s # is called n  s  the specular term; it describes the w i d t h of shiny highlights observed i n the material. k  a n d kd are referred to as the specular a n d diffuse coefficients,  respectively.  When k  s  = 0; kd — 1, the B R D F reduces to a constant, a n d  s  the surface is said to be diffuse, meaning that light is reflected u n i f o r m l y i n all directions. W h e n k = l;fc<f = 0, i n c o m i n g light is reflected according to s  the cosine lobe. A s such, some reflected directions are preferred over others, and so the reflected radiance is non-uniform over the surface, resulting i n a glint or highlight. A s the specular exponent n is increased, the highlights become taller (brighter) a n d more concentrated. preserves energy, its terms are normalized b y ^  I n order that the B R D F a n d ^ , a n d the coefficients  constrained to k + kd = 1. s  O n e can sample the P h o n g P D F b y sampling from the P D F n + 1  p(0,<t>) =  2TT  cos  (3.15)  33  Chapter 3. Monte Carlo Rendering  n(x)  n(x)  (a)  (b)  Figure 3.8: The Phong B R D F . (a) The amount of light a viewer sees depends on the angle 6 with respect to the reflected light direction R. (b) The diffuse and specular lobes. The specular lobe shape depends on the exponent n .  The pair (#, <f>) is a direction expressed in spherical coordinates; 8 is the polar angle constrained to the upper hemisphere (6 € [0,7r/2]), and' <f> is the azimuthal angle (4> € [0, 27r]). We can sample from Equation 3.15 by transforming a pair of uniform variates  (1*1,1*2)  according to  (0,0) = ( a r c c o s ( ( l - u i ) ^ ) , 2 7 R i ) . 2  (3.16)  Note that our sample directions (6,4>) are distributed about the polar axis. However, the specular term in the BRDF is parameterized about the reflection direction. So, the final step in sampling is to transform our directions from the global polar axis to the local frame of the BRDF. It is common practice to incorporate the cosine term from the illumination integral (Equation 3.1) into the BRDF, yielding fr(x,Ui  -> U ) = r  ^ ( n + l ) c o s 0 + ^ n{x)-Ui. n  27T  (3.17)  7T  Recall that the cosine term n(x) • u>i serves to scale incoming light based on the orientation of the surface with respect to the incoming light direction.  34  Chapter 3. Monte Carlo Rendering  Observe that the Phong lobe is parameterized about the reflection direction, whereas the cosine term is parameterized about the surface normal. This difference in parameterization makes sampling from the product of the B R D F and cosine term more challenging. In fact, to date, there is no analytical method for directly sampling from Equation 3.17 [38]. Instead, we sample according to f  by combing samples drawn exclusively from either  T  the specular term or the diffuse term. The ratio of samples drawn from each term depends on the values k and kd\ in other words, N B R D F samples s  will be drawn as k N specular samples and k^N diffuse samples. Deciding s  which term to sample from is made randomly in order to remain unbiased.  k N diffuse  k N specular  d  s  N samples  Figure 3 . 9 : Sampling from the product of the Phong B R D F and the area term involves combing samples drawn individually from the diffuse and specular terms.  In this case, k  s  = kj, — 0.5, resulting in (approximately) an equal  number of samples being drawn from each term.  Chapter 4. Sampling the Product Distribution  C h a p t e r  35  4  Sampling the Product Distribution The main contribution of this thesis addresses the problem of direct illumination from environment maps.  We propose a bidirectional sampling  approach in which both the energy distribution in the environment map and the reflectance of the BRDF are taken into account.  4.1  Bidirectional Importance Sampling  We operate on the assumption that creating samples from only the environment map or only the BRDF model is inexpensive, and that the visibility test dominates the cost. This assumption holds for scenes with complex geometry and for BRDF models optimized for sampling. In the case of procedural shaders, for example, importance sampling is always difficult, but we can still operate if the shader provides information to support efficient importance sampling as proposed by Slusallek et al. [40]. In any case, visibility is typically evaluated by casting a ray through the scene and intersecting it with geometry. A ray cast operation is orders of magnitude more expensive than drawing samples in virtually every case. Under these assumptions, one can benefit from extra time spent in at-  Chapter 4. Sampling the Product  36  Distribution  taining a good sample distribution that takes both the BRDF and environment map into account. Such a distribution selects only directions for visibility testing that have significant contribution to the reflected radiance of the surface under evaluation. Consider again the direct illumination integral L (to ) — / f (ui Jn r  r  r  —> u ) cos6iLi(<jJi)V(L0i)dwi,  (4.1)  r  where Li is the illumination from the environment map, f is the BRDF, and r  V is the binary visibility term. Conventional approaches perform importance sampling either from the intensity in the environment map according to the probability density  "  '•= T T T T T '  ( 4  -  2 )  or from the BRDF according to ,  \  frjUj  -> U> ) COS 6j r  J n fr ( i  ) COS t'jdcJj  u  Notice that in Equation 4.3 the cosine term has been included in the PDF, as discussed in Section 3.5. Let us derive estimators for the reflected radiance L .  When choosing  r  sample directions uiij ~ q^tOi), Equation 4.1 can be estimated with LAT,L:  r  /  \  LN,L{u ) r  fr( i,J  1  =  u  :  ^7 >  Ur) COS  9ijLi(Uij)V(uJi j) :  —  s  Lite)  L Li(ui)duji N  ;  Jn  •  .  l  K  l  )  Chapter 4. Sampling the Product  37  Distribution  When sampling according to the BRDF with u>ij ~ qf(u>i), the corresponding estimator L^j  is  1 ^  f (u j—>  1  fri^ij  r  9ijLi(u}ij)V(uJi j)  r  :  ^r)cosf?j jLi(o;jj)V'(a;j j) )  i  f {uiju ) r  1  N  u ) cos  it  cosOij  r  0ij<L;  f fr(u)i-*VT)CQ8 sl  i  JQ  ^ j= l  Evaluating the variance of these estimators yields r .2 L  -> /x  var(Ljv,i) = ^ - ^ ( f r i ,  r  x  J f (uJiu ) r  r  ^ i  -> w .)cosf? V'(wi)); T  i  cos 9i  2  When proposing samples from the environment, variance in the result is proportional to the variance in the BRDF. Similarly, when proposing from the BRDF, variance is proportional to the lights. It follows that the greatest reduction in image noise occurs when samples are drawn from the function with greater variance. This is verified by intuition. If the BRDFs are diffuse but the lighting spotty, then directions should be chosen according to the importance of the lights. On the other hand, if light sources in the environment map are relatively broad but the surfaces glossy or shiny, then proposing from the BRDF will be the better thing to do. Either approach will produce significant noise if both the BRDF and the illumination contain any high frequency information. The solution of Veach and Guibas [44] was to linearly combine samples drawn exclusively from the lights or BRDF. However, a mix of samples still suffers from dependence on  Chapter 4. Sampling the Product  38  Distribution  the variances of the individual techniques, because the variance of a sum of terms equates to the sum of the variances of each term. Our approach is to reduce variance by sampling directly from the product distribution,  J  f (oJi  n  r  —> U> ) COS 9 L (uJi)cLj r  i  i  '  i  Observe that the normalization term in Equation 4.4 is the direct illumination integral with the visibility term V(CJJ) omitted: In other words, this term is the exitant radiance in the absence of shadows. We refer to it as L  ns  ("radiance, no shadows"): L  n S  /  -=  fr(Vi  (4.5)  -> U ) COS 9iLi(u>i)duJi. r  Jn. If we draw sample directions CJJJ ~ p{^i) according to the product distribution in Equation 4.4, we can estimate Equation 4.1 with Ljv , where )P  1 ^  LN,M)  frjWjj  r  ^  _  J _ y ^ £ frjUij N  LNA*»T) =  COs9 L (0Jij)V(uJ j)  -> U )  =  i>j  i  i  ,  fhi  u  fr(Vij  r)  COS  -+UJ ) r  ^X>K>)-  9j Ljjujjj) tj  COS  V(uj,j)  .  O Li(uJi j) itj  t  (-) 4  6  i=i  We refer to LJV,  p  as the bidirectional  estimator  for the direct illumination  integral. The evaluation of Equation 4.6 can be interpreted as taking the unoccluded reflected radiance L  ns  and scaling it by the average result of N  visibility tests performed along directions that contribute most significantly to the radiance.  Chapter 4. Sampling the Product Distribution  39  Evaluating the variance, L = -^var(V(o; )) 2  coij ~ piyji) -> v&r(L ) N>p  i  Observe that the variance of the bidirectional estimator for the reflected radiance depends only on the variance in the visibility function. Variance in the lighting and BRDF has been completely factored away by sampling directly from the product distribution. The remaining question is how these samples can be drawn. Figure 4.1 contains angular plots of the probability densities corresponding to the various proposal distributions. The top image of Figure 4.1 is an image of the importance derived from environment map intensity. The center image is a plot of the BRDF for a surface patch whose normal is pointing at the horizon. The BRDF plotted here is a specular Phong BRDF with exponent 50. The bottom image is the product of the E M and BRDF. It is this product distribution that we'd like to sample from, because it directly corresponds to the reflected radiance of the surface patch. Figure 4.2 illustrates the benefit of sampling from the product distribution, as opposed to sampling from the E M or the BRDF individually. Sampling from the BRDF alone misses the bright lights in the environment (top image), while samples drawn from the environment do not lie in the BRDF (middle image). The bottom image, however, shows samples drawn according to the product distribution using our SIR technique, as described in Section 4.4.  Chapter  4. Sampling the Product  Distribution  40  Figure 4.1: From top to bottom: angular plots of the importance function of the Grace Cathedral E M , a specular Phong B R D F of exp 50, and their product.  Chapter  4. Sampling the Product  Distribution  41  F i g u r e 4.2: Samples drawn solely from the B R D F (top) or the environment (middle) vastly undersample the product distribution (bottom). The sample set in the bottom image was generated with our SIR technique; see Section 4.4.  Chapter 4. Sampling the Product  4.2  Distribution  42  The Challenge of Bidirectional Sampling  The challenge in realizing the idea of bidirectional importance sampling is that the product of the BRDF and environment map is too high dimensional to precompute. The BRDF is a 4D function that maps from incoming directions to outgoing directions. The relevant 2D slice of the BRDF, corresponding to a specific outgoing light direction u> , varies from point to point r  in the scene due to changes in the local surface orientation. The environment map is two dimensional, and thus the BRDF-EM product has six dimensions. Even with a coarse discretization of the BRDF, which might cause high frequency features in the BRDF to be lost, precomputing the product distribution and storing it in a table for sampling is obviously prohibitively expensive. Complicating matters further is the fact that the BRDF and environment map are usually specified relative to different coordinate frames. Typical parameterizations for the BRDF are local, either about the surface normal or the reflected light direction. The E M , on the other hand, is expressed in a global frame. Thus, precomputing the product of the BRDF and the environment for a single BRDF orientation is insufficient. One must either compute a sample rotation or add another two dimensions to the precomputed table, corresponding to a discrete sampling of possible surface orientations. The high-dimensional nature of the problem makes it difficult to use precomputation for sample generation. We simply cannot afford to premultiply the two functions together. We suggest the following process for sampling from the product of the lights and the BRDF. First, create samples according to either the environment map or the BRDF. Then, adjust the sample distribution such that the  Chapter 4. Sampling the Product  43  Distribution  directions chosen for visibility testing will be proportional to the product distribution. We have developed two solutions that realize this approach, one based on rejection sampling and the other on the sampling-importance resampling (SIR) algorithm. Note that the complete algorithm is a two-stage approach. That is, the local illumination integral is always estimated with importance sampling, but individual subproblems are solved with either rejection sampling or SIR. Our two sampling methods are detailed in the following two sections.  4.3  Sample Generation Through Rejection  cq(Xi)  x ~q(x) L  x  F i g u r e 4.3: Rejection sampling. A sample Xi ~ q{x) is accepted as being a valid sample of the target distribution p{x) if a uniform variate i n [0, cq(x)] falls under p(x).  Our first approach for sampling from the product distribution is through rejection sampling. Rejection sampling is a technique which allows us to  Chapter 4. Sampling the Product  Distribution  44  sample from a probability distribution known only up to a proportionality constant.  This allows us to sample according to distributions, where the  normalization constant is too expensive to compute explicitly (as is the case with the product distribution from Equation 4.4). Also, rejection sampling provides a simple procedure for sampling from arbitrary distributions that do not have simple analytical forms for sampling. Rejection sampling is illustrated in Figure 4.3. In order to create samples tOij  ~  c-q(ui)  p(u)i),  we can approximate  P(LO%)  with a PDF  q(wi),  such that  for some constant c. We then generate random samples u>ij ~  and accept them with a probability of  p(tOij)/(c  bound the product distribution  fL.  qi^i),  • q(u>ij)).  Figure 4.4: Rejection sampling the product distribution. Both Lmaxf  p(u>i) <  However,  fmaxL  fmaxL  and  is a tighter fit  for this particular B R D F and E M . Thus, the acceptance rate is higher by proposing samples from the lights and rejecting against the max of the B R D F .  In our case, one simple way of bounding p(u>j) from Equation 4.4 is to  Chapter 4. Sampling the Product  Distribution  45  use the global maximum of the BRDF as a conservative estimate for its contribution as illustrated in Figure 4.4.  If L  and f  m&x  max  refer to the  maximum values of both distributions, i.e. £ m a x := m a x ^ L ^ i ) ,  and  /max := m a x ^ g / ^ i ) ,  then we can bound the target PDF  as p(uJi)  < /max -QL^i)-  Note that  qL is just the usual importance from the environment map alone, and so we can sample from it in constant time using CDF pre-integration and inversion or the alias method, as discussed in Section 3.2.2. In order to accept i V visibility samples, we will thus have to create M « /  • N environment  max  map samples Uij through importance sampling, and then accept each sample individually with probability f (iOij) cos 8ij • j Li{u>i)dlOi r  •  u  /max '  L  ns  Both this formula and the final radiance estimate from Equation 4.6 require the normalization term L  ns  from Equation 4.5. We can estimate  this term using information that has already been computed during the ) rejection sampling. Since we already evaluate both BRDF and environment map for M directions  ~  Lns « ^  L  i  ^  qi,{uji),  (  h  j  i  we can estimate L  ns  J2 f K , r  -  as  U, ) COS 9^. 0  (4.7)  i=i  Another interpretation of this method is that we estimate the unoccluded illumination L  ns  with M samples, using importance sampling from  the environment map. However, we only evaluate the visibility for i V of  Chapter 4. Sampling the Product  Distribution  46  those samples whose light contribution is significant enough to make visibility tests worthwhile. The directions for the visibility tests are chosen in an unbiased fashion. So far, we have bounded the actual target PDF as a constant times the environment map PDF. This is appropriate if the BRDF contains mostly low frequencies, i.e. if f  max  is a close bound of the real BRDF. If this is  not the case, then most samples will be rejected and the rejection sampling will become inefficient. Instead, we can perform the same rejection sampling algorithm by approximating the environment map with a conservative bound and then selecting samples according to the real BRDF. Under this scheme, we now have p(wj) < L  max  • qf(ui),  which amounts to generating  samples from the BRDF alone, and then rejecting according to the product distribution as before. Given these two ways of rejection sampling, we usually want to draw the initial samples in such a way that the bounding constant is minimized. That is, if /  m a x  < L  max  we importance sample from the environment map;  otherwise, we importance sample from the BRDF. In practice, we randomly choose which of the two methods to use, where the method with the smaller bounding constant is chosen with a higher probability. This is similar in spirit to the combination of sampling strategies proposed by Veach and Guibas [44]. Our rejection sampling approach has worked well in the experiments that we have performed (see Chapter 5). However, the inherent downside of using rejection sampling is that one cannot guarantee bounds on the execution time for creating a new sample. That is, if the area between c • q(u)i,j) and p(tOij) is large, the probability of sample acceptance will be low. One way of dealing with this is to choose a maximum number of sample attempts in  Chapter 4. Sampling the Product  47  Distribution  the rejection sampling. So long as at least one sample is accepted, this will still yield an asymptotically unbiased estimate of the reflected radiance. If no samples are accepted, a possible strategy could be to test visibility for a random subset. Another possibility which is less expensive is to use the unoccluded illumination, where visibility has not been tested at all. While this introduces bias, it will only happen in very dark areas, where the product of illumination and BRDF is very small. In these areas, the visibility term will not have significant impact anyway. In practice, even in the presence of high specular BRDFs and environments, the rejection sampling acceptance probability has been sufficiently high so as to not require resorting to these biased methods. See Chapter 5 for a complete discussion. The next section presents our second method for sampling from the product distribution, which does not suffer from the unbounded execution time of the rejection sampling.  4.4  Sample Generation Through Sampling-Importance Resampling  In this section, we describe our second technique for sampling from the product distribution by motivating with a description of the importance sampling algorithm. Recall that the idea behind Monte Carlo integration is to approximate the integral /(/) = J f(x)dx  by sampling from a target density p(x) and  computing an estimate to /(/), N (4.8)  48  Chapter 4. Sampling the Product Distribution In our case, p(x) is the BRDF-EM product distribution of Equation 4.4. The idea behind the importance  sampling  algorithm [35] is to sample  from p[x) by introducing an arbitrary proposal distribution q(x) that is easier to sample from than p(x) and whose support includes p(x). Equation 4.8 can be rewritten as  / N  where W(XJ) oc  is a measure of the importance of sample Xj. Conse-  quently, if one can draw samples according to q(x) and evaluate  W(XJ),  a  possible Monte Carlo estimate of /(/) is  (4.9) This leads to the sampling-importance 41, 43].  resampling  (SIR) algorithm [10,  The goal of SIR is to introduce a resampling step to eliminate  samples with low importance ratios w(x). SIR proceeds in two stages. In the first stage, a set of independent random samples X = {xi,... ,XM} is drawn from a proposal distribution q(x). In the second stage, a smaller set of samples Y = {yi, • • •, J/AT} is drawn from X with sample probabilities w(xi) proportional to their importance ratio j^ry- As the number of first-round samples M approaches infinity, the sample set Y can be shown to have been drawn directly from p. A requirement for the algorithm to be efficient is that q{x) be a good approximation to p(x) in the region of visibility. We can apply SIR as a solution to the problem of drawing samples from the lighting-BRDF product distribution. It is desirable for q(x), the proposal density for the first round of sampling, to be a reasonable approximation of  Chapter 4. Sampling the Product Distribution  49  p{x). Because p(x) in our case is the product distribution, natural choices for q(x) are probability densities based on either the BRDF or the environment map. In the second SIR stage, we resample the first-round draws according to w(xi) oc ^ry- The basic intuition here is that directions are proposed from a single distribution, then resampled (filtered) based on their value to the second distribution.  Stage 2  Stage 1 L(co ) itj  CO  L(co ) i2  CO  i,2  CO  i,N  Figure 4 . 5 : Our resampling method. First, M samples are proposed from q/, the P D F of the B R D F . The candidate directions are then resampled based on the incoming light along those directions, producing N samples for visibility testing. JV is generally much less than M.  This leads to our second method for approximating the illumination integral, which proceeds in two stages. In the first stage, M sample directions are drawn according to qL-(wi) or  qf(u>i).  We store the sample and as-  sociated weight w(uij), chosen to be w{u>ij) = f (^i,j) r  when u>ij ~ qL^i)  Chapter 4. Sampling the Product Distribution and  w(u>ij)  = Li(u>ij) when  Uij  50  ~ qf(u)i). The second stage then resamples  the stored sample directions according to the distribution of their weights w(u>i).  Visibility is tested along each of the iV resampled directions.  The total number of samples generated for each pixel is exactly M + N. This an improvement over rejection sampling for two reasons. Firstly, execution time is tightly bounded. We no longer have to wait an indeterminate amount of time (possibly even indefinitely) for the acceptance criteria to pass and accept a sample as being valid. Using the SIR algorithm, samples can be drawn directly from the product distribution in constant time. The second improvement over rejection sampling is that the sample sizes M and TV can be freely chosen, yielding fine control over the tradeoff between quality and time. For example, the BRDF sample size can be adjusted based on the expense of sampling from the BRDF model. Also, it is possible to directly select N - the target number of visibility rays traced per pixel based on, for example, scene complexity. Because the cost of ray tracing typically dominates rendering time, our general approach for generating results has been to fix N and adjust M so as to increase or decrease the variance in resampled directions. As discussed in the next chapter, typical values of M we use are one to two orders of magnitude larger than N. It is interesting to note that prior strategies for sampling, which draw directly from the E M or BRDF alone, are just special cases of our SIR technique where M = N — 1. For example, proposing samples only from the lights can be achieved by SIR sampling with a single first-round sample drawn from the E M , and then trivially resampled in the second round with probability 1.0 because it is the only sample.  Chapter 5. Results  C h a p t e r  51  5  Results This chapter presents the results of our techniques in comparison with previous sampling strategies for rendering from environment maps. While both the rejection sampling and SIR methods increase the cost of sample generation, we demonstrate significant quality improvements for the same compute time under the assumption of BRDF representations that support efficient evaluation. In all of our tests, the illumination comes from high dynamic range (HDR) environment maps, where the contrast ratio between bright and dark areas of the environment is high (several orders of magnitude). In principle, our techniques would work equally well under low dynamic range (LDR) lighting, where the ratio of bright to dark is at most 255 : 1. However, LDR EMs are not of sufficiently high-frequency to be interesting. Furthermore, the trend in rendering is to use HDR environments, which serve as a more accurate representation of real lighting conditions. The images in this chapter were generated with our ray tracer which uses a voxel grid as the acceleration data structure for intersection queries. Our comparisons examine the output quality of the various rendering algorithms for a fixed amount of computation time. We performed these tests on a single-processor 3.0GHz P4 running Linux.  Chapter 5. Results  52  Figures 5.1, 5.2 and 5.3 contain images of Michaelangelo's David in the Grace Cathedral environment. We use here the 700k-triangle version of the David model acquired from the Digital Michaelangelo Project [13]. In our implementation, intersecting a ray with the David model takes on average 6.1 ps on our test machine. The Grace Cathedral environment is a 1024x512 HDR map with a contrast ratio of 10 : 1. In these tests, each algorithm 7  was given 12.0 seconds to render a 176x248 image. This image resolution was chosen in order to be able to distinguish differences between the images when presented in print form. In these three figures, the images on the top row show previous sampling strategies, while the images on the bottom row present our techniques. In each case, the upper-left image was rendered with 100 samples drawn from the environment map, while the upper-right image was rendered with 75 samples from the BRDF. The images on the bottom row were each rendered with 15 visibility samples, chosen from 800 candidates; on the left is the SIR approach, and on the right, rejection sampling. Figure 5.1 shows the David rendered with a Phong BRDF model with coefficients k = 1.0, kd = 0.0 and a specular exponent of 10. The BRDF s  in this case is a broad-lobed gloss. Observe that sampling from the E M produces a nicer image than sampling from the BRDF. This is because Grace is a highly specular environment map, and the BRDF here is glossy but still broad; drawing samples according to the BRDF is less informative than sampling according to the high-frequency lighting. The images generated by our techniques, however, are much better than either simple sampling strategy.  Observe also that our methods are able to achieve this quality  with several times fewer ray casts - in this case only 15 as opposed to 75 or 100.  Chapter 5. Results  53  Figure 5.1: David in Grace Cathedral. Phong exp 10, k = 1.0, kd = 0.0. a  176x248 image rendered in 12.0 seconds. Clockwise from top left: 100 E M samples; 75 B R D F samples; 15/800 rejection samples; 15/800 SIR samples.  Chapter 5. Results  Figure 5.2:  54  David in Grace Cathedral. Phong exp 50, fc = 1.0, kd = 0.0. s  176x248 image rendered in 12.0 seconds. Clockwise from top left: 100 E M samples; 75 B R D F samples; 15/800 rejection samples; 15/800 SIR samples.  Chapter 5. Results  55  Figure 5.3: David in Grace Cathedral. Phong exp 50, k = 0.5, kd — 0.5. s  176x248 image rendered in 12.0 seconds. Clockwise from top left: 100 E M samples; 75 B R D F samples; 15/800 rejection samples; 15/800 SIR samples.  56  Chapter 5. Results  In Figure 5.2, David's BRDF is still a purely specular Phong (k — s  1.0, kd — 0.0), but here the Phong exponent has been increased to 50. Here, the model's surface is very shiny; accordingly, sampling purely from the BRDF produces a better image than sampling from the environment map. Again, however, our techniques go further. Both the lighting and the BRDF in this scene are high-frequency, and yet our methods produce images that are of much higher quality. In Figure 5.3, David is again given a shiny Phong BRDF with exponent 50. This time, however, we add a diffuse term (k = 0.5, kd = 0.5). Observe s  how noisy the image produced by sampling purely from the BRDF has become. BRDF sampling makes little sense in this case, because the presence of a diffuse term sends samples in all directions, making the variance high. Our methods still outperform the straightforward sampling strategies, because by sampling directly from the product distribution we disregard those directions which the BRDF accepts but the E M does not. Figure 5.4 shows David appearing with the same Phong BRDF as in Figure 5.3, but under a different lighting scheme, the St. Peter's Basilica environment map. Like the Grace Cathedral environment, St. Peter's is a 1024x512 E M with a contrast ratio of 10 : 1. However, the content of 7  St. Peter's is somewhat lower frequency than Grace, and as such sampling solely from the E M does not perform as well for this scene as it did in Grace. Consequent, our technique stands out even more.  Chapter 5. Results  57  Figure 5.4: David i n St. Peters. Phong exp 50, k = 0.5, k = 0.5. 176x248 s  d  image rendered i n 12.0 seconds. Clockwise from top left: 100 E M samples; 75 B R D F samples; 15/800 rejection samples; 15/800 SIR samples.  Chapter 5. Results  58  How do our techniques do in situations where the B R D F is very low frequency? The images in Figure 5.5 shows a purely diffuse David (Phong coefficients k = 1.0, kd = 0.1)) in Grace Cathedral. On the left is the image s  resulting from 100 samples drawn purely E M samples. O n the right is the SIR technique with 25 samples ( M = 800, N = 25). In this case, because the B R D F is purely diffuse and thus smoothly varying, one might think that sampling from the E M alone would be all that is necessary to render a decent image. Nevertheless, our SIR technique produces a comparable image with a quarter of the rays, indicating that sampling from the product distribution is always desirable.  Figure 5.5: David in Grace Cathedral. Diffuse Phong (k — 0.0, kd = 1.0). s  176x248, 12.0 seconds. Left: 100 E M samples; Right: 25/800 SIR samples.  Chapter 5. Results  59  Figure 5.6 contains unzoomed images of David in Grace Cathedral. Figures 5.7 and 5.8 contain a comparison of our algorithm to that of Veach and Guibas [44], in which samples drawn from the individual E M and B R D F distributions are combined via a weighted sum to produce an estimate of the reflected radiance. Our technique is able to achieve superior image quality with far fewer rays. The results obtained from the combined sampling approach are an improvement over straightforward sampling approaches; however, our method goes further to reduce variance than simply blending between the variances in the individual distributions. In summary, the results presented here clearly demonstrate that the approach of sampling directly from the product distribution outperforms previous sampling strategies. What is more, we are able to achieve comparable quality with far fewer rays, meaning that our techniques are particularly beneficial to renderings of complex scenes where ray-scene intersection queries are expensive.  Figure 5.6: David i n Grace Cathedral, viewed unzoomed. Phong exp 50, k — 0.5, kd = 0.5. 228x272 image rendered in 18.0 seconds. From left to s  right: 110 E M samples; 13/900 SIR samples; 80 B R D F samples  Chapter 5. Results  Figure 5.7: Comparison between Veach and Guibas [44] (left column, 100 samples) and our SIR technique (right column, 15/800 samples). 176x248 image rendered in 12.0 seconds. Specular Phong (k = 0.5, kd = 0.5), exponents a  10 (top row) and 50 (bottom row).  60  Chapter  5.  61  Results  Figure 5.8: Comparison between Veach and Guibas [44] (left column, 100 samples) and our SIR technique (right column, 15/800 samples).  176x248  image rendered in 12.0 seconds. Phong B R D F , exp 50, k = 1.0, k<i = 0.0 s  Chapter 5. Results  62  Figure 5.9: Our SIR algorithm compared to the full solution. 268x220 image of the Dragon model (871k triangles) i n Grace Cathedral.  Top: converged  image (hours); Bottom: SIR algorithm using 150/3000 samples (90 sec)  Chapter 6.  C h a p t e r  Conclusions  63  6  Conclusions We have presented two Monte Carlo strategies for sampling the incident illumination from environment maps that take into account the product of the light and the surface reflectance.  By providing a means of sampling  from a more complex target distribution, our methods are able to achieve lower variance in the renderings of scenes with high frequency lighting and specular BRDFs, as compared to traditional importance sampling strategies. Note that although our proposed bidirectional methods take longer to generate samples than simpler approaches, the number of samples required to achieve good quality is considerably less than when sampling according to a simple function. For large datasets with complex structures, the time required to trace shadow rays will dominate the rendering time. In such cases, our methods will provide greater benefit over importance sampling from the E M or BRDF alone. On-going work for us is the examination of other sampling strategies that exist in the literature, such as iterative SIR, Metropolis-Hastings, and particle filtering [2]. The general idea behind these strategies is to use samples that have already been drawn as a basis for proposing further, fitter samples. It would be interesting to explore how these methods of sampling from more complicated distributions could be applied to other problems in computer graphics.  64  Bibliography  Bibliography [1] S. Agarwal, R. Ramamoorthi, S. Belongie, and H. W. Jensen. Structured importance sampling of environment maps. ACM Transactions on Graphics (Proc. Siggraph), 22(3):605-612, July 2003. [2] C. Andrieu, N. de Freitas, A. Doucet, and M. Jordan. An introduction to M C M C for machine learning, 2003. [3] J. Blinn and M. Mewell. Texture and reflection in computer generated images. Communications of the ACM, 19(10):542-547, 1976. [4] D. Burke, A. Ghosh, and W. Heidrich. Bidiectional importance sampling for illumination from environment maps. ACM Siggraph 2004 Sketches and Applications, August 2004. [5] J. Cohen and R Debevec.  Light-gen.  HDRshop plugin, 2001.  http://www.ict.usc.edu/jcohen/lightgen/lightgen.html. [6] P. Debevec.  Rendering synthetic objects into real scenes: Bridging  traditional and image-based graphics with global illumination and high dynamic range photography. In Proc. of ACM Siggraph '98, pages 189198, July 1998. [7] P. Debevec and J. Malik. Recovering high dynamic range radiance maps from photographs. In Proc. of ACM Siggraph '97, pages 369-378, 1997.  65  Bibliography  [8] O. Deussen, S. Hiller, C. van Overveld, and T. Strothotte. Floating points: A method for computing stipple drawings. In Proc. of Eurographics, pages 41-50, August 2000. [9] P. Dutre and H. W. Jensen. Siggraph Course Notes: State of the Art in Monte Carlo Global Illumination. 2004. [10] A. Gelman, J. Carlin, H. Stern, and D. Rubin. Bayesian data analysis. Chapman and Hall, London, 1995. [11] C. Goral, K. Torrance, D. Greenberg, and B. Battaile. Modeling the interaction of light between diffuse surfaces. In Proc. of ACM Siggraph '84, pages 213-222, 1984. [12] N. Greene. Environment mapping and other applications of world projections. IEEE CG&A,  6(ll):21-29, 1986.  [13] Stanford Graphics Group.  Digital Michaelangelo project,  2001.  http://graphics .Stanford .edu/projects/mich / . [14] J. Hammersley and D. Handscomb. Monte Carlo Methods. Methuen, London,1964. v  [15] W. Heidrich and H.-P. Seidel. Realistic, hardware-accelerated shading and lighting. In Proc. of ACM Siggraph '99, pages 171-178, August 1999. [16] H. W. Jensen. Global illumination using photon maps. In Eurographics Workshop on Rendering, pages 21-30, 1996. [17] H. W. Jenseri. Siggraph Course Notes: State of the Art in Monte Carlo Ray Tracing for Realistic Image Synthesis: 2001.  66  Bibliography  [18] H. W. Jensen. Siggraph Course Notes: Monte Carlo Ray Tracing. 2003. [19] M . Kalos and P. Whitlock. Monte Carlo Methods. John Wiley & Sons, New York, 1986. [20] J. Kautz and M. McCool. Approximation of glossy reflection with prefiltered environment maps. In Proc. of Graphics Interface, pages 119— 126, May 2000. [21] J. Kautz, P.-P. Vazquez, W. Heidrich, and H.-P. Seidel. Unified approach to prefiltered environment maps. In Eurographics Workshop on Rendering, pages 185-196, 2000. [22] A. Keller. Quasi-Monte Carlo Methods for Photorealistic Image Synthesis. PhD thesis, Aachen, 1998. [23] T. Kollig and A. Keller. Efficient illumination by high dynamic range images. In Eurographics Symposium on Rendering, pages 45-51, June 2003. [24] E . Lafortune, S.-C. Foo, K. Torrance, and D. Greenberg. Non-linear approximation of reflectance functions. In Proc. of ACM Siggraph '97, pages 117-126, August 1997. [25] J. Lawrence, S. Rusinkiewicz, and R. Ramamoorthi. Efficient BRDF importance sampling using a factored representation.  ACM Transac-  tions on Graphics (Proc. Siggraph), 23(3):496-505, August 2004. [26] R. Lewis. Making shaders more physically plausable. Computer Graphics Forum, 13(2):109-120, 1994.  67  Bibliography  [27] X. Liu, P. Sloan, H.-Y. Shum, and J. Snyder. All-frequency precomputed radiance transfer for glossy objects. In Eurographics Symposium on Rendering, June 2004. [28] G. Marsaglia. Generating discrete random variables in a computer. Communications of the ACM, 6(l):37-38, January 1963. [29] R. Ng, R. Ramamoorthi, and P. Hanrahan. All-frequency shadows using non-linear wavelet lighting approximation. ACM Transactions on Graphics (Proc. Siggraph), 22(3):376-381, July 2003. [30] R. Ng, R. Ramamoorthi, and P. Hanrahan. Triple product wavelet integral for all-frequency relighting. ACM Transactions on Graphics (Proc. Siggraph), 23(3):477-487, August 2004. [31] V. Ostromoukhov, C. Donohue, and P.-M. Jodoin. Fast hierarchical importance sampling with blue noise properties. ACM Transactions on Graphics (Proc. Siggraph), 23(3):488-495, August 2004. [32] B. Phong. Illumination for computer generated images.  Communica-  tions of the ACM, 18(6):311-317, June 1975. [33] R. Ramamoorthi and P. Hanrahan. An efficient representation for irradiance environment maps. In Proc. of ACM Siggraph '01, pages 497500, 2001. [34] R. Ramamoorthi and P. Hanrahan. Frequency space environment map rendering. In Proc. of ACM Siggraph '02, pages 517-526, 2002. [35] R. Rubinstein. Simulation and the Monte Carlo Method. John Wiley and Sons, New York, 1981.  Bibliography  68  [36] A. Secord, W. Heidrich, and L. Streit. Fast primitive distribution for illustration. In Eurographics Workshop on Rendering, pages 215-226, June 2002. [37] P. Shirley. Realistic Ray Tracing. A K Peters, Natick, 2000. [38] P. Shirley. Personal communication. May 2004. [39] P. Sloan, J. Kautz, and J. Snyder. Precomputed radiance transfer for real-time rendering in dynamic, low-frequency environments. In Proc. of ACM Siggraph '02, pages 527-536, 2002. [40] P. Slusallek, T. Pfiaum, and H.-P. Seidel. Using procedural renderman shaders for global illumination. In Proc. of Eurographics, volume 14, pages 311-324, 1995. [41] A. Smith and A. Gelfand. Bayesian statistics without tears: a samplingresampling perspective, volume 46, pages 84-88, 1992. [42] L. Szecsi, M . Sbert, and L. Szirmay-Kalos. Combined correlated and importance sampling in direct light source computation and environment mapping. In Proc. of Eurographics, volume 23, pages 585-593, 2004. [43] M . Tanner. Tools for statistical inference. Springer Verlag, New York, 3rd edition, 1996. [44] E . Veach and L. Guibas. Optimally combining sampling techniques for monte carlo rendering. In Proc. of ACM Siggraph '95, pages 419-428, August 1995.  Bibliography [45] M . Vose.  69  A linear algorithm for generating random numbers with  a given distribution.  IEEE  Transactions on Software Engineering,  17(9):972-975, September 1991. [46] A. Walker. An efficient method for generating discrete random variables with general distributions. ACM Transactions on Mathematical Software, 3(3):253-256, September 1977. [47] G. Ward. The RADIANCE lighting simulation and rendering system. In Proc. of ACM Siggraph '94, pages 459-472, August 1994.  Appendix  A p p e n d i x  A. Details on the Alias Method for Sampling  70  A  Details on the Alias M e t h o d for Sampling As mentioned in Section 3.2.2, we use the alias method of Vose [45] for sampling from discrete random distributions. Surprisingly, Vose's work appears to have been overlooked by the graphics community, where the transformation method is the de facto standard for importance sampling. The alias method offers a number of advantages over CDF inversion, the most notable being a fast 0(iV)-time initialization algorithm. In his paper, Vose's description of the approach contains ambiguous pseudocode; we therefore present a class called D i s c r e t e D i s t r i b u t i o n which is similar to the class used in our implementation. D i s c r e t e D i s t r i b u t i o n contains two methods: a constructor for initialization, and a method for drawing samples called sample(). The constructor operates on an input array of non-negative values which are used as distribution weights: D i s c r e t e D i s t r i b u t i o n : .'DiscreteDistribution(double *probs, unsigned num): m_num( num ) , m_probs( probs ) { // C l a s s i f y the normalized weights as being // e i t h e r l a r g e ( > 1/N ) or small ( <= 1/N ). const double stdWeight = 1.0 / m_num; Vector s m a l l , l a r g e ; for  ( unsigned i = 0 ; i < m_num ; ++i ) if  ( m_probs[i] > stdWeight )  Appendix  A. Details on the Alias Method for Sampling  71  large.addObject( i ); else small.addObject( i ); // A s s i g n a l i a s i n d i c e s . m_alias = new unsigned[m_num]; while ( ! l a r g e . isEmptyO && ! small. isEmptyO ) { const unsigned sm = small.getAndRemoveLastObject(); const unsigned l g = large.getAndRemoveLastObject(); m_probs[lg] -= stdWeight - m_probs[sm]; m_probs[sm] *= static_cast<double>( m_num ); m_alias[sm] = l g ; if  ( m_probs[lg] > stdWeight ) large.addObject( l g );  else small.addObject( l g );  > for  ( unsigned i = 0 ; i < s m a l l . s i z e ( ) ; ++i )  m_probs[ s m a l l [ i ] ] = 1.0; f o r ( unsigned i = 0 ; i < l a r g e . s i z e ( ) ; ++i ) m_probs[ l a r g e [ i ] ] = 1.0;  > One can now draw samples from the distribution by calling the sample () method, which returns an index in [0, JV — 1]. unsigned C D i s c r e t e D i s t r i b u t i o n : : s a m p l e ( ) const { // Generate a random number i n [0, m_num). const double u = uniform() * m_num; const unsigned i = static_cast<unsigned>( u ); r e t u r n ( u - i > m_probs[i] ) ? m _ a l i a s [ i ] }  : i ;  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items