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

MONTE-CARLO METHODS IN GLOBAL ILLUMINATION Script written by


MONTE-CARLO METHODS IN GLOBAL ILLUMINATION

Script written by

Szirmay-Kalos L? szl? a o
in WS of 1999/2000

Institute of Computer Graphics Vienna University of Technology

i

Contents
1 Introduction 1.1 Global pass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Local pass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Tone mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Global illumination problem 2.1 The rendering equation . . . . . . . . . . . . . . . . . . 2.2 Measuring the radiance . . . . . . . . . . . . . . . . . . 2.3 The potential equation . . . . . . . . . . . . . . . . . . 2.4 Measuring the potential . . . . . . . . . . . . . . . . . . 2.5 The rendering problem . . . . . . . . . . . . . . . . . . 2.5.1 Geometry of the surfaces . . . . . . . . . . . . . 2.5.2 Bi-directional Re?ection Distribution Functions . 2.5.3 Lightsources . . . . . . . . . . . . . . . . . . . 2.5.4 Measuring devices . . . . . . . . . . . . . . . . 2.6 Numerical solution of the rendering equation . . . . . . 2.6.1 Error measures for numeric techniques . . . . . 2.6.2 Properties of the rendering equation . . . . . . . 2.7 Classi?cation of the solution techniques . . . . . . . . . Optical material models 3.1 Diffuse re?ection . . . . . . . . . . . . . . . . . . . 3.2 Ideal, mirror-like re?ection . . . . . . . . . . . . . . 3.3 Ideal refraction . . . . . . . . . . . . . . . . . . . . 3.4 Non-ideal, specular re?ection . . . . . . . . . . . . . 3.4.1 Phong re?ection model and its modi?cations 3.4.2 Cook-Torrance model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 2 2 5 5 8 8 9 10 10 11 11 12 14 14 15 15 18 18 19 21 21 21 24 28 28 28 28 30 31 31 32 33 34 34 35 37 37 37 39

2

3

4

Solution strategies for the global illumination problem 4.1 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Expansion of the rendering equation: gathering walks . 4.2.2 Expansion of the potential equation: shooting walks . 4.2.3 Merits and disadvantages of expansion methods . . . . 4.3 Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Analysis of the iteration . . . . . . . . . . . . . . . . 4.4 Analytical solution of the rendering equation . . . . . . . . . 4.4.1 Scenes with constant radiance . . . . . . . . . . . . . 4.4.2 Scenes with constant re?ected radiance . . . . . . . .

5

Finite-element methods for the Global Illumination Problem 5.1 Galerkin’s method . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Point collocation method . . . . . . . . . . . . . . . . . . . . . . 5.3 Finite element methods for the diffuse global illumination problem 5.3.1 Geometric methods for form factor computation . . . . .

ii

CONTENTS

iii 41 42 42 43 46 47 49 49 50 51 52 53 53 54 54 57 58 60 60 61 61 62 62 63 66 66 67 68 69 70 71 72 73 74 76 76 79 80 81 82 82 82 83 85 86 86 88 88 90 90 91 91 92 92 93 95

6

Numerical quadrature for high dimensional integrals 6.1 Monte-Carlo quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Quasi-Monte Carlo quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Error Analysis for integrands of ?nite variation: Koksma-Hlawka Inequality . 6.2.2 Generation of the sample points . . . . . . . . . . . . . . . . . . . . . . . . 6.2.3 Generation of low-discrepancy sequences . . . . . . . . . . . . . . . . . . . 6.3 Importance sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Generation of a random variable with a prescribed probability density . . . . 6.3.2 Importance sampling in quasi-Monte Carlo integration . . . . . . . . . . . . 6.3.3 Metropolis sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.4 Application of the VEGAS algorithm . . . . . . . . . . . . . . . . . . . . . Random walk solution of the global illumination problem 7.1 Why should we use Monte-Carlo expansion methods? . . . . 7.2 Quasi-Monte Carlo quadrature for the rendering equation . . 7.2.1 Integrating functions of unbounded variation . . . . 7.3 Importance sampling for the rendering equation . . . . . . . 7.3.1 BRDF sampling . . . . . . . . . . . . . . . . . . . 7.3.2 Lightsource sampling . . . . . . . . . . . . . . . . . 7.3.3 Sampling the lightsources in gathering random walks 7.3.4 Importance sampling in colored scenes . . . . . . . 7.3.5 Multiple importance sampling . . . . . . . . . . . . 7.4 Handling in?nite-dimensional integrals . . . . . . . . . . . 7.4.1 Russian roulette . . . . . . . . . . . . . . . . . . . . 7.4.2 Russian roulette in quasi-Monte Carlo quadrature . . Review of random walk algorithms 8.1 Gathering-type random walk algorithms . . . . . . . 8.1.1 Ray-casting . . . . . . . . . . . . . . . . . . 8.1.2 Visibility ray-tracing . . . . . . . . . . . . . 8.1.3 Distributed ray-tracing . . . . . . . . . . . . 8.1.4 Path-tracing . . . . . . . . . . . . . . . . . . 8.2 Shooting-type walks methods . . . . . . . . . . . . . 8.2.1 Photon tracing . . . . . . . . . . . . . . . . 8.2.2 Light-tracing . . . . . . . . . . . . . . . . . 8.2.3 Random walks for the radiosity setting . . . 8.3 Bi-directional random walk algorithms . . . . . . . . 8.3.1 Bi-directional path-tracing . . . . . . . . . . 8.3.2 Metropolis light transport . . . . . . . . . . 8.3.3 Photon-map . . . . . . . . . . . . . . . . . . 8.3.4 Instant radiosity . . . . . . . . . . . . . . . . 8.4 Global methods . . . . . . . . . . . . . . . . . . . . 8.4.1 Multi-path method using global random lines 8.4.2 Global ray-bundle tracing . . . . . . . . . . 8.4.3 Preprocessing the point lightsources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

8

9

Iteration solution of the global illumination problem 9.1 Why should we use Monte-Carlo iteration methods? . . . . . . 9.2 Formal de?nition of stochastic iteration . . . . . . . . . . . . 9.2.1 Other averaging techniques . . . . . . . . . . . . . . . 9.2.2 Can we use quasi-Monte Carlo techniques in iteration?

10 Review of stochastic iteration algorithms 10.1 Stochastic iteration for the diffuse radiosity . . . . . . . . . . . . . . . . . . . . 10.1.1 Stochastic radiosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.1.2 Transillumination radiosity . . . . . . . . . . . . . . . . . . . . . . . . . 10.1.3 Stochastic ray-radiosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.2 De?nition of the random transport operator for the non-diffuse ?nite-element case 10.2.1 Single ray based transport operator . . . . . . . . . . . . . . . . . . . . . 10.2.2 Stochastic iteration using ray-bundles . . . . . . . . . . . . . . . . . . .

CONTENTS

iv 97 99 99 99 99 100 100 101 101 102 103 103 104 104 104 105 105 105 106 106 107 107 107 107 107 109 109 110 111 116

11 Implementation of the path-tracing algorithm 11.1 Vector module . . . . . . . . . . . . . . . . 11.1.1 Point3D class . . . . . . . . . . . . 11.1.2 Transformation class . . . . . . . . 11.2 Container module . . . . . . . . . . . . . . 11.3 Color module . . . . . . . . . . . . . . . . 11.4 Material models . . . . . . . . . . . . . . . 11.4.1 Diffuse material class . . . . . . . . 11.4.2 Ideal mirror class . . . . . . . . . . 11.4.3 Ideal refracting material class . . . 11.4.4 Specular material class . . . . . . . 11.4.5 General material class . . . . . . . 11.5 Light module . . . . . . . . . . . . . . . . 11.5.1 Emitter class . . . . . . . . . . . . 11.5.2 Positional light class . . . . . . . . 11.6 Model module . . . . . . . . . . . . . . . . 11.6.1 Primitive class . . . . . . . . . . . 11.6.2 Object class . . . . . . . . . . . . . 11.6.3 Virtual world class . . . . . . . . . 11.7 Camera module . . . . . . . . . . . . . . . 11.8 Scene module . . . . . . . . . . . . . . . . 11.8.1 Scene class . . . . . . . . . . . . . 11.9 Dynamic model of path tracing . . . . . . . 11.9.1 Finding the visible primitive . . . . 11.9.2 Detecting the visible lightsources . 11.9.3 Direct lightsource computation . . . 11.9.4 Path tracing . . . . . . . . . . . . . 11.9.5 Rendering complete images . . . . BIBLIOGRAPHY SUBJECT INDEX

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 1

Introduction
The ultimate objective of image synthesis or rendering is to provide the user with the illusion of watching real objects on the computer screen (?gure 1.1). The image is generated from an internal model which is called the virtual world. To provide the illusion of watching the real world, the color sensation of an observer looking at the arti?cial image generated by the graphics system must be approximately equivalent to the color perception which would be obtained in the real world. The color perception of humans depends on the light power reaching the eye from a given direction and on the operation of the eye. The power, in turn, is determined from the radiance of the visible points. The radiance depends on the shape and optical properties of the objects and on the intensity of the lightsources. In order to model this complex phenomenon, both the physical-mathematical structure of the light-object interaction and the operation of the eye must be understood.
rendering observer of the computer screen monitor R G B power λ color perception in the nerve cells Tone mapping power λ radiance λ measuring window device virtual world

power λ observer of the real world real world
Figure 1.1: Tasks of rendering

radiance λ

The image synthesis uses an internal model consisting of the geometry of the virtual world, optical material properties and the description of the lighting in the scene (?gure 1.2). From these, applying the laws of physics (e.g. Maxwell equations) the real world optical phenomena can be simulated to ?nd the light distribution in the scene. This step is called the view-independent step or the global pass of rendering. Then a measurement device, called the eye or camera, is placed into the scene and the light distribution is measured from a given location and orientation. This is called the view-dependent step or the local pass. Note that not all rendering algorithms make a clear distinction between the determination of the view-independent light distribution and the measurement of this 1

1.1. GLOBAL PASS

2

distribution by the camera, but simultaneously compute the light distribution and its effect on the camera. Rendering results in a representation of the perceived image, which is usually the collection of pixel colors or some discrete sampling of the radiance function. The exact simulation of the light perceived by the eye is impossible, since it would require endless computational process. On the other hand, it is not even worth doing since the possible distributions which can be produced by computer screens are limited in contrast to the in?nite variety of real world light distributions. Consequently, color perception is approximated instead of having a completely accurate simulation. The accuracy of this approximation is determined by the ability of the eye to make a distinction between two light distributions. Computer screens can produce controllable electromagnetic waves, or colored light, mixed from three separate wavelengths for their observers. Thus in the ?nal step of image synthesis tone mapping is needed which converts the computed color or radiance information to the R, G, B intensities that can be produced by the color monitor.
geometry of the virtual world material properties lighting camera global rendering (global pass) radiance of surface points image calculation (local pass) radiance of pixels tone mapping R,G,B of pixels

Figure 1.2: Data?ow of rendering

1.1 Global pass
The global pass determines the light re?ected off the surface points at different directions. Since light is an electromagnetic wave, light distribution in a point and at a given direction can be represented by a wavelength-dependent ? function [Abr97, K? n85]. Rendering algorithms usually evaluate this functions at a few representative waveo lengths. On a given wavelength the intensity of the light is described by the radiance. In scenes not incorporating participating media it is enough to calculate the radiance at surface points. The radiance re?ected off a surface point is affected by the emission of this point (lighting), the illumination provided by other surface points and the optical properties of the material at this point (material properties). Formally this dependence is characterized by a Fredholm type integral equation of the second kind, which is called the rendering equation. From mathematical point of view, global pass is the solution of this integral equation for the representative wavelengths.

1.2 Local pass
The local pass means the measurement of the global radiance function by a camera. A camera is a collection of light measuring devices which usually correspond to pixels in the image. A certain measuring device is characterized by a sensitivity function that describes which points and directions may affect the device.

1.3 Tone mapping
Light is an electromagnetic wave, and its color is determined by the eye’s perception of its spectral energy distribution. Due to its internal structure, the eye is a very poor spectrometer since it actually samples and integrates the energy in three overlapping frequency ranges by three types of photopigments according to a widely accepted (but also argued) model. As a consequence of this, any color perception can be represented by three scalars (called ? tristimulus values) instead of complete functions [Abr97, K? n85, Nem90]. o A convenient way to de?ne the axes of a coordinate system in the three-dimensional space of color sensations is to select three wavelengths where one type of photopigments is signi?cantly more sensitive than the other two [SK99c]. This method has been devised by Grassmann, who also speci?ed a criterion for separating the three representative wavelengths. Grassmann laws state that the representative wavelengths should be selected such that no one of them can be matched by the mixture of the other two in terms of color sensation (this criterion is

1.3. TONE MAPPING

3

similar to the concept of linear independence.) An appropriate collection of the representative wavelengths is:

red = 645 nm; green = 526 nm; blue = 444 nm:

(1.1)

Now let us suppose that monochromatic light of wavelength  is perceived by the eye. The equivalent portions of red, green and blue light, or (r, g , b) tristimulus values, can be generated by three color matching functions (r(), g () and b()) which are based on physiological measurements. Note the negative section of r() (and to a less extent in g ()) in ?gure 1.3. It means that not all colors can be represented by positive (r, g , b) values.
R=645nm, G=526nm, B=444nm matching functions 3.5 r(lambda) g(lambda) b(lambda)

3

2.5

2 r,g,b

1.5

1

0.5

0

-0.5 400 450 500 550 lambda[nm] 600 650 700

Figure 1.3: Mean 10-deg color matching functions of Stiles and Burch: r(), g (), b().

If the perceived color is not monochromatic, but is described by an L() distribution, the tristimulus coordinates are computed using the assumption that the sensation is produced by an additive mixture of the perceptions of elemental monochromatic components:

r = L()  r() d; g = L()  g() d; b = L()  b() d:
  

Z

Z

Z

(1.2)

For computer generated images, the color sensation of an observer watching a virtual world on the screen must be approximately equivalent to the color sensation obtained in the real world. Since color sensations are represented by (r, g , b), it means that the tristimulus values should be identical. If two energy distributions are associated with the same tristimulus coordinates, they produce the same color sensation, and are called metamers. In computer monitors and in television screens three phosphor layers can be stimulated to produce red, green and blue light. The objective, then, is to ?nd the necessary stimulus to produce a metamer of the real energy distribution of the light [Sch96, BS95]. This stimulus can be controlled by the (R, G, B ) values of the actual pixel. The (r, g , b) matching functions of ?gure 1.3 depend on the wavelength of the selected primaries, which are not necessarily identical to the wavelengths on which our monitor can emit light. This requires the conversion of tristimulus values by a linear transformation. The calculation of pixel R; G; B values thus consists of the following steps. First the spectrum associated with the pixel is computed. Then the spectrum is matched by three standard color matching functions de?ned by three primaries. Finally, the standard color coordinates are transformed to the monitor color coordinates taking into account the monitor properties. In practice, the standard color system is usually the CIE XYZ system [WS82] which uses three hypothetical primaries to allow the de?nition of any color by positive weights. The linear transformation that converts from the XYZ system to the monitor RGB system can be obtained from the X; Y; Z coordinates of the emissions of the three phosphors and of the white point of the monitor. For a monitor with standard NTSC phosphors and white point, the following transformation can be used [Gla95]:

2 R 3 2 1:967 4 G 5 = 4 0:955
B
0:064

0:548 0:297 X 1:938 0:027 5  4 Y 0:130 0:982 Z

3 2

3 5:

(1.3)

1.3. TONE MAPPING

4
X,Y,Z matching functions 3.5 X(lambda) Y(lambda) Z(lambda)

3

2.5

2 X,Y,Z

1.5

1

0.5

0

-0.5 400 450 500 550 lambda[nm] 600 650 700

Figure 1.4: Mean 10-deg color XY Z matching functions of Stiles and Burch: X (), Y (), Z ()

The whole computation of the (R, G, B ) values in order for the monitor color to be a metamer of the calculated spectrum is called tone mapping. The (R, G, B ) values are positive numbers usually in the range of [0...255] if 8 bits are available to represent them. Unfortunately, not all colors can be reproduced on the computer screen, because of the negative sections of the color matching functions and due to the fact that the number of available intensity levels is usually much less than can be perceived in the real world. Thus tone mapping is also responsible for optimally selecting from the available intensity levels for color reproduction. The mapping from the computed levels to the available ones can be either linear or logarithmic. The latter takes advantage of the logarithmic characteristics of the human perception system [PP98].

Chapter 2

Global illumination problem
In this chapter the mathematical model of the light-surface interaction is presented. This mathematical model is an integral equation, which has to be solved to obtain physically accurate images.

2.1 The rendering equation
Hereinafter, monochromatic light of a representative wavelength  will be assumed, since the complete color calculation can be broken down to these representative wavelengths. The parameters of the equations usually depend on the wavelength, but for notational simplicity, we do not always include the  variable in them. In this section, we brie?y review the measures of the light transport and the mathematical formulation that can compute them.

z θ

ω

z

ω

dθ sinθ dφ

y φ x

dθ dφ x

y

Figure 2.1: De?nition of directions in a spherical coordinate system (left) and calculation of differential solid angles (right)

The directional property of the light energy emission is described in a so-called illumination sphere
or in illumination hemisphere
H which contain those solid angles to where the surface point can emit energy. The surface of transparent materials can emit in any directions of a sphere, while the surface of opaque materials can transfer energy only to the hemisphere that is “above” the surface. Setting up a spherical coordinate system (?gure 2.1), a direction ! can be de?ned by two angles ; , where  is the angle between the given direction and the z -axis, and  is the angle between the projection of the given direction onto the x; y plane and the x-axis. Sets of directions are de?ned by solid angles. By de?nition, a solid angle is a cone or a pyramid, with its size determined by its subtended area of a unit sphere centered around the apex (?gure 2.2). A differential (in?nitesimal) solid angle can also be given by a vector d~ , where the vector equals to a direction of the differential set. ! A differential solid angle can also be expressed by the ;  angles. Suppose that  is modi?ed by d and  is by d. During this the directional vector scans a differential rectangle having d vertical and sin   d horizontal sizes (right of ?gure 2.1), thus the size of the solid angle is

d! = sin   dd: (2.1) The solid angle, in which a differential dA surface can be seen from point p, is the projected (visible) area ~ per the square of the distance of the surface (?gure 2.2). If the angle between the surface normal of dA and the
5

2.1. THE RENDERING EQUATION

6

θ dω dA r

Figure 2.2: De?nition of the solid angle

directional vector from dA to ~ is , and the distance from dA to p is r, then this solid angle is: p ~

d! = dA rcos  : 2

(2.2)

The intensity of the energy transfer is characterized by several metrics in computer graphics depending on whether or not the directional and positional properties are taken into account. The light power or ?ux  is the energy radiated through a boundary per unit time over a given range of the spectrum (say [;  + d]). Since a photon has  = energy where  is the Planck-constant, the power can be h h supposed to be propotional to the number of photons that go through the boundary in a unit time. The power is not always a convenient measure since it also needs the de?nition of a boundary. We can get rid of this additional information if the boundary is de?ned in a differential way focusing on a single surface point and a single direction. The resulting measure is called the radiance. The radiance or intensity L(~ ; ! ) is the differential light ?ux (~ ; dA; !; d! ) leaving a surface element dA x x around ~ in a differential solid angle d! around ! per the projected area of the surface element dA and the size of x the solid angle d! . If the angle of the surface normal and the direction of interest is , then the projected area is dA  cos , hence the radiance is:

x dA; L(~ ; !) = d(~; d!  !; d!) : x (2.3) dA cos  ~ Since a differential surface area can also be seen as a vector dA, which is parallel to the surface normal at the given
point, the radiance can also be obtained in a simpler form

x dA; L(~ ; !) = d(~ ; ~ !; d!) : x dA  d~ !
. dA r
Figure 2.3: Energy transfer between two differential surface elements

(2.4)

θ



θ’ . dA’

Having introduced the most important metrics, we turn to their determination in the simplest case, where there are only two differential surface elements in the 3D space, one (dA) emits light energy and the other (dA0 ) absorbs it (?gure 2.3). According to the de?nition of the radiance (equation (2.3)), if dA0 is visible from dA in solid angle d! and the radiance of the surface element dA is L in this direction, then the ?ux leaving dA and reaching dA0 is:

d = L  dA  d!  cos :
Using equation (2.2), the solid angle can be expressed from the projected area of dA0 , thus we get:

(2.5)

0  dA0 d = L  dA  cos  r2  cos  :

(2.6)

2.1. THE RENDERING EQUATION

7

This formula is called the fundamental law of photometry. Note that if equation (2.2) is used again for the emitter, the transferred power can also be written in the following form: 0 dA0 d = L  dA  cos   2  cos  = L  dA0  d!0  cos 0 : (2.7)

r

Thus similar formula applies for the patch that receives the power as for the patch that emits it. In light-surface interaction the surface illuminated by an incident beam may re?ect a portion of the incoming energy in various directions or it may absorb the rest. It has to be emphasized that a physically correct model must maintain energy equilibrium, that is, the re?ected and the transmitted (or absorbed) energy must be equal to the incident energy. Optically perfect or smooth surfaces will re?ect or transmit only coherent components governed by the laws of geometric optics, including the law of re?ection and the Snellius–Descartes law of refraction. The surface irregularities, however, re?ect or refract the incident light incoherently in any direction. Since the exact nature of these irregularities is not known, light-surface interaction is modeled by means of probability theory. Assume that a photon comes from the direction denoted by ! 0 to point ~ . The probability of re?ection or x refraction at ~ into solid angle d! around ! is expressed by the following probability density function, also called x as transfer probability density:

w(!0 ; ~ ; !)  d! = Prfphoton is re?ected or refracted to d! around ! j coming from !0 g: x

(2.8)

Note that this probability distribution is a mixed, discrete-continuous distribution, since the probability that the light is re?ected exactly to the ideal re?ection direction may be non-zero.

’ h(x, - ω)

ω L(x, ω)

θ’

ω’’

’ L(h(x, - ω) , ω’)

x

Figure 2.4: Geometry of the rendering equation

The light ?ux (out ) leaving the surface at solid angle d! around ! consists of the emission and re?ected (or refracted) components. In order to obtain the re?ected/refracted component, let us assume that photons are coming to area dA around ~ from solid angle d!0 around !0 and their total represented power is in (~ ; dA; !0 ; d!0 ). The probability that a x x photon is re?ected to d! is w(! 0 ; ~ ; ! )  d! thus the expected power leaving the surface element after re?ection or x refration is: To obtain the total re?ected/refracted power, all the possible incoming directions ! 0 should be taken into account and their contributions should be integrated:

w(!0 ; ~ ; !) d!  in (~ ; !0; d!0 ): x x

Z




(w(!0 ; ~ ; !) d!)in (~ ; !0 ; d!0 ): x x

(2.9)

If the surface itself emits energy, i.e. if it is a lightsource, then the emission also contributes to the output ?ux:

e (~ ; !) = Le (~ ; !)  dA  cos   d!: x x
Adding the possible contributions we obtain:

(2.10)

out (~ ; !) = e (~ ; !) + x x

Z



(w(!0 ; ~ ; !) d!)in (~ ; !0 ; d!0 ): x x

(2.11)

2.2. MEASURING THE RADIANCE

8

The ?uxes can be expressed by the radiant intensities according to equations (2.5) and (2.7), thus:

in (~ ; !0; d!0 ) = Lin(~ ; !0 )  dA  cos 0  d!0 ; x x out (~ ; !; d! ) = L(~ ; ! )  dA  cos   d!:  x x
Substituting these terms into equation 2.11 and dividing both sides by dA  d!  cos  we get:

(2.12)

x L(~ ; !) = Le (~ ; !) + Lin (~ ; !0 )  cos 0  w(! ; ~; !) d!0 : x x x cos



Z

0

(2.13)

Let us realize that Lin (~ ; ! 0 ) is equal to L(~ ; ! 0 ) if ~ is the point that is visible from ~ at direction ! 0 , usually x y y x expressed by the so called visibility function h ( i.e. ~ = h(~ ; ! 0 )). y x Using these equations and introducing the bi-directional re?ection/refraction function — or BRDF for short — as de?ned by the following formula

x fr (!0 ; ~ ; !) = w(! ; ~; !) ; x cos
we obtain the following fundamental law, called the rendering equation:

0

(2.14)

L(~ ; !) = Le (~ ; !) + L(h(~ ; !0); !0 )  fr (!0 ; ~ ; !)  cos 0 d!0 : x x x x


Introducing the following notation for the integral operator

Z

(2.15)

(T L)(~ ; !) = x

Z




L(h(~ ; !0); !0 )  fr (!0 ; ~ ; !)  cos 0 d!0 ; x x L = Le + T L:

(2.16)

the short form of the rendering equation can also be established: (2.17)

In fact, every color calculation problem consists of several rendering equations, one for each representative wavelength. Optical parameters (Le ; fr ) obviously vary in the different equations, but the geometry represented by function h does not.

2.2 Measuring the radiance
Having solved the rendering equation, the radiance can be determined for any point and direction. To obtain an image, the power that affect different parts of light sensitive material (retina or the ?lm of a camera) must be determined. Mathematically each distinct part is associated with a measuring device, thus the collection of these devices is, in fact, the model of the human eye or cameras. A measuring device is characterized by a sensitivity function W e (~ ; ! 0 ), which gives a scaling value C if a y light-beam of unit power leaving ~ in direction ! 0 contributes to the measured value and 0 otherwise (for example, y this device can measure the power going through a single pixel of the image and landing at the eye, or leaving a surface element at any direction). Since the power leaving surface area d~ around ~ in solid angle d! around ! is L(~; ! )d~ cos d! , the meay y y y sured value by a certain device is:

ZZ

Operator M is called the radiance measuring operator.


S

d(~; !)  W e (~; !) = y y

ZZ

S

L(~; !) cos   W e (~; !) d~ d! = ML: y y y

(2.18)

2.3 The potential equation
Looking at the light propagation as an interaction between an emitter and a receiver, the radiance describes this interaction from the point of view of the emitter. The same phenomenon can also be explained from the point of view of the receiver, when the appropriate measure is called the potential. By de?nition, the potential W (~ ; ! 0 ) y

2.4. MEASURING THE POTENTIAL

9

y

W(y, ω ) ω W(h(y, ω’) , ω ) h(y, ω’) ω’’ θ

Figure 2.5: Geometry of the potential equation

expresses the effect of that portion of the unit power light-beam emitted by point ~ in direction ! 0 , which actually y lands at a given measuring device either directly or indirectly after some re?ections or refractions. Using probabilistic terms the potential is the product of the scaling value C and the probability that a light-beam emitted by a point in a given direction reaches the measuring device. If only direct contribution is considered, then W (~ ; ! 0 ) = W e (~ ; ! 0 ). To take into account light re?ections, y y we can establish the potential equation [PM95]: In this equation integral operator T 0 — which is the adjoint of T — describes the potential transport

W = W e + T 0 W:

(2.19)

(T 0 W )(~; !0 ) = y

Z

To prove this equation, we can use the probabilistic interpretation of the potential. Denoting the “unit power light-beam that is leaving ~ at direction ! ” by (~ ; ! ), we can write: x x

where  is the angle between the surface normal and the outgoing direction ! .




W (h(~; !0 ); !)  fr (!0 ; h(~; !0 ); !)  cos  d!; y y

(2.20)

W (~; !0 ) = C  Prf(~; !0 ) lands at the deviceg = y y 0 ) goes directly to the deviceg + C  Prf(~ ; ! 0 ) lands at the device after at least one re?ectiong: C  Prf(~; ! y y

(2.21)

The ?rst term is equal to the sensitivity function. Considering the second term, we can obtain:

Z



C  Prf(~; !0 ) lands at the device after at least one re?ectiong = y

C  Prf(h(~; !0 ); !) lands at the deviceg  Prf(~; !0 ) is re?ected to d! around ! at h(~; !0 )g = y y y

Z



W (h(~; !0 ); !)  fr (!0; h(~; !0 ); !)  cos  d!: y y

(2.22)

2.4 Measuring the potential
Alternatively to the radiance, the power arriving at the measuring device can also be computed from the potential. Since

takes into account all possible emission points and directions is

de (~; !0 ) = Le (~; !0 )  cos  d~ d!0 y y y 0 and the probability that the photons of this is proportional to the power of the light-beam emitted by d~ in d! y beam really go either directly or indirectly to the measuring device is W (~ ; ! 0 )=C , the total measured value that y C

ZZ

S

y de (~; !0 )  W (~; ! ) = y C

0

ZZ

where M0 is the potential measuring operator. Note that unlike the radiance measuring operator, the potential measuring operator integrates on the lightsource.


S

W (~; !0 )  Le (~; !0)  cos  d~ d!0 = M0 W; y y y

(2.23)

2.5. THE RENDERING PROBLEM

10

2.5 The rendering problem
Generally, the rendering problem is a quadruple [Kel97]:

where S is the geometry of surfaces, fr is the BRDF of surface points, Le is the emitted radiance of surface points at different directions and e is a collection of measuring functions. Rendering algorithms aim at modeling and simulation of light-surface interactions to ?nd out the power emitted by the surfaces and landing at the measuring devices. The light-surface interaction can be formulated by the rendering equation or alternatively by the potential equation. Since according to the de?nition of the radiance the total power of the scene is the integral of the radiance function

hS; fr (!0 ; ~ ; !); Le (~ ; !); We (~ ; !)i x x x

W

=

ZZ

S

L(~; !) d~ d~ ; y y !

we are interested only in those solutions where the total integral of the radiance function is ?nite. Formally the solution is searched in the L1 function space. Recall that the measured value can be computed by the measuring function from the radiance

ZZ

where M is the radiance measuring operator. Let us introduce the scalar product hu; v i


S

L(~; !) cos   W e (~; !) d~ d! = ML; y y y

(2.24)

hu; vi =

ZZ

where  is the angle between the surface normal and direction ! and thus d~ cos  is the in?nitesimal visible area y from direction ! . Using this scalar product, we can obtain an alternative form of the measuring operator:


S

u(~; !)  v(~; !) d~ d~ = y y y !

ZZ


S

u(~; !)  v(~; !) d~ cos  d!; y y y

ML = hL; W e i:

The potential measuring operator can also be given in a scalar product form:

M0 W = hLe ; W i:

(2.25)

Let us examine a single re?ection of the light. The measured value taking into account a single re?ection can be expressed in two ways: hT Le ; W e i = hLe ; T 0 W e i: (2.26) Thus the radiance and potential transport operators are adjoint operators [M? t81]. a

2.5.1 Geometry of the surfaces
A surface is a set of 3D points which are de?ned by an equation [SK99c]. Those points are said to be in this set, which satisfy the de?nition equation. The equation can produce the points on the surface either implicitly, when the general form is

F (x; y; z ) = 0;

or explicitly, which is preferred in rendering algorithms. The general form of an explicit surface de?nition is:

~ = ~(u; v); u 2 [0; 1]; v 2 [0; 1]: r r (2.27) Points on the surface can be generated by selecting u; v parameters from the unit interval and substituting their values into the function ~(u; v ). For example, a sphere of radius R and of center (x0 ; y0 ; z0 ) can be de?ned either r
by the following explicit equation:

x = x0 + R  cos 2u  sin v; y = y0 + R  sin 2u  sin v; z = z0 + R  cos v; u; v 2 [0; 1];
or by an implicit equation:

(x x0 )2 + (y y0 )2 + (z z0 )2 R2 = 0:

Some rendering algorithms do not work with the original geometry but approximate the surfaces by triangle meshes. This approximation — also called the tessellation — can be realized by selecting n  m points in parameter space u 2 [0; 1]; v 2 [0; 1] and adding the

[~(ui ; vj );~(ui+1 ; vj );~(ui+1 ; vj+1 )] and [~(ui ; vj );~(ui+1 ; vj+1 );~(ui ; vj+1 )] r r r r r r triangles for all i = 1 : : : n 1 and j = 1 : : : m 1 indices to the resulting list (?gure 2.6). For the discussion of
surface de?nition methods using forward and reverse engineering and of transformations refer to [SK99c, Chi88, VMC97, RVW98] and to [SKe95, Kra89, Her91, Lan91], respectively.

2.5. THE RENDERING PROBLEM

11

original surface with isolines

selecting points and triangles

result of tessellation

Figure 2.6: Tessellation process

2.5.2 Bi-directional Re?ection Distribution Functions
The Bi-directional Re?ection Distribution Functions (or BRDFs) characterize the optical material properties. Photorealistic rendering algorithms require the BRDFs not to violate physics. Such BRDF models must satisfy both reciprocity and energy balance, and are called physically plausible [Lew93]. Reciprocity that was recognized by Helmholtz is the symmetry property of the BRDF characterizing re?ections, which is de?ned by the following equation [Min41]:

fr (!; ~ ; !0 ) = fr (!0 ; ~ ; !); x x (2.28) where ! 0 is the vector pointing towards the incoming light and vector ! de?nes the viewing direction. Reciprocity
is important because it allows for the backward tracing of the light as happens in visibility ray-tracing algorithms. Note that reciprocity is valid if the incoming and outgoing beams are in the same material, that is, re?ection BRDFs are expected to be reciprocal but refraction BRDFs are not. Suppose that the surface is illuminated by a beam from direction ! 0 . Energy balance means that the number of outgoing photons cannot be more than how many were received. To examine it formally, we can introduce the albedo [Lew93] that is the ratio of the total re?ected power and incoming power, or equivalently, the probability of the re?ection to any direction. Energy balance means that the albedo cannot be greater than 1:

a(~ ; !0 ) = Prfphoton is re?ected j coming from !0 g = x

Z

Z

H


H

w(!0 ; ~ ; !) d! = x
(2.29)

fr (!0 ; ~ ; !)  cos  d!  1: x

If the BRDF is reciprocal, then the albedo can also be interpreted as the optical response of the surface to homogeneous sky-light illumination of unit intensity:

T1=

Z

Using the de?nition of the transfer probability density w in equation (2.8), we can obtain the following probabilistic interpretation of the BRDF:


H

1  fr (!0 ; ~ ; !)  cos 0 d!0 = a(~ ; !): x x

(2.30)

1 fr (!0 ; ~ ; !)  cos  = d!  Prfphoton is re?ected or refracted to d! around ! j it comes from !0 g: x
Different BRDF models are discussed in [SKe95, SK99c, SK99a, Pho75, Bli77, CT81, War92, HTSG91, ON94, Sch93, NNSK99b, NNSK99a, CSK98].

2.5.3 Lightsources
The lightsources are surfaces that have non-zero Le emission. Rendering algorithms also use abstract lightsources that are not associated with surfaces [SKe95]. The most important types of these abstract lightsource models are the following:



a point-lightsource is located at a given point of the virtual world and its size is zero. The direction of the illumination provided by this lightsource at a point is the vector from the point to the location of the lightsource and the intensity of the illumination decreases with the square of the distance. An example of this category is an electric bulb.

2.5. THE RENDERING PROBLEM

12

 

a directional-lightsource is located at in?nity in a given direction. The direction and the intensity of its illumination are constant everywhere. A directional-lightsource is, in fact, equivalent to an in?nitely distant plane emitter. The sun can be considered as a directional-lightsource. sky-light illumination provides constant illumination at any point and at any direction.

2.5.4 Measuring devices
In order to establish models for the camera, the operation of the human eye can be analyzed when the human is looking at the real world and when he is watching the computer screen.

?e e

ω

θ Φ

?p

y ?e

? p Lp pixel p Φp

|y - e |
watching the real world watching the computer screen

Figure 2.7: A model of the human eye

The human eye contains a lens called the pupil of size e (?gure 2.7). We shall assume that the pupil is small, which is often referred as the pinhole camera model. When the human is watching the computer screen, a pixel p is visible in a small solid angle
p . To provide the illusion that the screen is a replacement of the real world, the p power emitted by the pixel towards the eye should be similar to that  power which would come from the real world and would land at the pupil from solid angle
p . If the emission intensity of the pixel is Lp , then the power landing at the eye from this pixel is where e is the angle between the surface normal on the pupil and the direction of the pixel. We have to ?nd a camera model that computes a measured value P which is stored in the image-buffer to control the monitor. The response of monitors can be characterized by a function B  R, where B represents a scaling according to the brightness settings and R might be non-linear. In order to compensate the monitor’s non-linearity, the look-up table is set to distort the values of the image-buffer with the inverse of this non-linear mapping, that is by R 1 . This distortion is called gamma-correction. Using gamma-correction, the radiant intensity of the pixel is: Since we require that p

p = Lp  e  cos e 
p ;

= , our model camera is expected to provide the following measured value:     P = R R 1 Lp = Lp = e  cos  
 B : B B e p

Lp = B  R(R 1 (P )) = B  P:

(2.31)

Let us assign a measuring device for this pixel. This device is sensitive to those surface points that are visible in this pixel — i.e. in
p — and to those directions that point from the surface points to the pupil. Mathematically, this can be represented by the following measuring function

W e (~; !) = y
where

(C
0

if ~ is visible in
p and ! points from ~ through e; y y otherwise,

(2.32)

1 C = e  cos  
 B : e p

The measuring device provides the following value:

P = ML =

ZZ


S

L(~; !)  W e (~; !)  cos  d~d!: y y y

(2.33)

2.5. THE RENDERING PROBLEM

13

Applying equation (2.2), we can use the following substitutions:

where d~ is a differential area on the pupil and j~ ~j is the distance between the pupil and the visible surface. e y e Substituting these and taking advantage that the pupil e is small, we obtain

2 cos y y e d! = d~  j~ ej2 ; d~ = j~cos~j  d!p ; e y ~ e 

P=

ZZ


p e

y e2 cos L(h(~; !p); !p )  C  j~ ej2  j~cos ~j  cos  d~d!p = e e y ~ e 

ZZ

Z

where eye is the position of the very small pupil, and !p is the direction which points from ~ to the eye position. ~ y Note that the measured value is independent of both the distance and the orientation of the visible surface! This is explained by the fact that as we move farther from a surface, although the power coming from a unit surface area decreases with the square of the distance, the area that can be visible in a given solid angle increases with the same speed. The two factors compensate each other. Similarly, when we look at a surface from a non perpendicular direction, the radiance is weighted by the cosine of this angle in the power, but the surface area that can be seen is inversely proportional to the cosine of the same angle.
window


p

L(h(eye; !p); !p )  C  cos e  e d!p = L(h(eye; !p ); !p ) 
1 B d!p ; ~ ~

p

Z


p e

L(h(~; !p ); !p )  C  cos e d~d!p  e e
p
(2.34)

θp dωp eye f
: focal distance

?p p
pixel

?p θ dωp eye | y - eye | y

Figure 2.8: Evaluation of the measured value as an integral on the pixel (left) and on the surface (right)

Since
p is the collection of those directions that go through the pixel, the measured value can also be expressed as an integral over the pixel area Sp (?gure 2.8). Let ~ be the running point on the pixel, p be the angle between p the pixel normal and vector p eye, and pix be equal to p in the center of the pixel. Since jp eyej = f= cos p ~ ~ ~ ~ where f is the focal distance, i.e. the distance between the eye and the plane of the window, we can establish the following correspondance:

p p 3 p d!p = d~  cos 2 = d~  cos p : jp eyej ~ ~ f2
p

(2.35)

Substituting this to equation (2.34), the measured value becomes

P = L(h(p; !p ); !p ) 
1 B  cos 2p d~: ~ ~ ~ p f
Sp
Let us introduce camera parameter c(p) by ~

Z

3

(2.36)

Since


p =

Z

p

c(p) =
Sp B  cos 2p : ~ f

3

d!p =

Z

p

(2.37)

Sp

cos3 p d~  Sp  cos3 pix ; f2 p f2

(2.38)

the camera parameter can also be expressed in the following exact and approximating forms:

3 1 p cos3 c(p) = B SR cos3 p d~  B cos p3  B : ~  cos pix p p

(2.39)

Sp

2.6. NUMERICAL SOLUTION OF THE RENDERING EQUATION

14

The approximations are accurate if the pixel is small compared to the focal distance. If image space ?ltering is also applied, then Sp may be greater than the pixel and the camera parameter is not necessarily constant but follows the shape of a reconstruction ?lter [SKe95, SK95]. Summarizing the measured value is an integral on the pixel:

P = L(h(p; !p ); !p )  cSp) d~: ~ ~ ~ (~ p
Sp p

Z

(2.40)

Equation (2.34) can also be evaluated on object surfaces. Using the notations of ?gure 2.8, we can obtain:

y d!p = jd~  cos j2 ~ eye y ~  P = L(~; !~!eye)  V
p (~) 
1 B  j~ cos~ j2 d~: y y ~ y y eye y
S p
Let us introduce surface dependent camera parameter g (~) by y

(2.41)

y where ~ is a surface point visible in the pixel. Let the indicator function of those points that are visible in
p be V
p (~). The measured value is then: y

Z

(2.42)

f g(~) =
 B  j1 eyej2 = B  j~ eyej2  R cos3  d~  B  j~ eyej2f  S  cos3  y ~ ~ y y ~ y ~ p p p p pix
Sp
using formula (2.38) that expresses
p . Thus the ?nal form of the measured value is:

2

2

(2.43)

P = L(~; !~!eye)  V
p (~)  g(~)  cos  d~: y y ~ y y y
S
Comparing this formula with equation (2.33), the following sensitivity function can be derived:

Z

(2.44)

W e (~; !) = y

( (! !~!eye)  g(~) y y ~
0
otherwise.

if ~ is visible through the pixel; y

(2.45)

2.6 Numerical solution of the rendering equation
2.6.1 Error measures for numeric techniques
When solving the rendering problem numerically, the required radiance, potential or measured power can only be approximated. To characterize the accuracy of the approximation, error measures are needed. The errors of radiance or potential are called object-space error measures. The error of the measured power is called imagespace error measure. A good error measure should be invariant to changes that modify neither the object space nor the camera. Particularly, the error should not change when a surface is subdivided to two surfaces or the complete scene is scaled — i.e. it should be tessellation and scaling independent —, or when a pixel is subdivided to increase the image resolution — i.e. it should be resolution independent. Tessellation and scaling independent measures can be de?ned if the norm of the error function is divided by the total surface area. The norm [M? t81, ST93] is selected from the following possibilities: a

jjf jj1 =

ZZ

S

vZ Z u jf (~ ; !)j d~ d~ ; jjf jj2 = u (f (~ ; !))2 d~ d~ ; jjf jj1 = max jf (~ ; !)j: x x ! x x ! x t x ~ ;!

S

(2.46)

~ Using any of these, if L is the calculated radiance and L is the real radiance, then an absolute object-space error a and a relative object-space error r are ~ ~ a = jjL S Ljj ; r = jjL ~ Ljj : jjLjj
(2.47)

2.7. CLASSIFICATION OF THE SOLUTION TECHNIQUES

15

Image space norms are based on the powers that go through different pixels p = 1; 2; : : : ; NP .

jjjj1 =

NP X p=1

vN uX u P jp j; jjjj2 = t jp j2 ; jjjj1 = max jp j: p
p=1

(2.48)

Resolution independent measures divide the total error of pixels by the number of pixels (P ). If the calculated and ~ real powers that go through pixel p are p and p , respectively, then an absolute and a relative image space error measures are

~ ~ a = jjN jj ; r = jj ~ jj :
P

jjjj

(2.49)

2.6.2 Properties of the rendering equation
The integral operators of both the rendering and the potential equations are contractions. This statement can be proven from the fact that for physically plausible models, the albedo is less than 1. Here, the 1-norm is used:

kT Lk1 = max jT Lj  max jLj  max jT 1j = max jLj  max ja(~ ; !)j = kLk1  max ja(~ ; !)j: x x kT 0 W k1 = max jT 0 W j  max jW j  max jT 0 1j = max jW j  max ja(h(~; !0 ); !0 )j = y kW k1  max ja(~ ; !)j: x

For the potential, similar results can be obtained:

2.7 Classi?cation of the solution techniques
In order to ?nd the color of a pixel, the radiance has to be evaluated for that surface which is visible through the given pixel. The determination of this surface is called the hidden surface problem or visibility calculation. Having solved the visibility problem, the surface visible in the given pixel is known, and the radiance may be calculated on the representative wavelengths by the rendering equation.

surface 1 surface 2 pixel

eye surface 3

Figure 2.9: Multiple re?ections of light

Due to multiple re?ections of light beams, the calculation of the radiance of the light leaving a point on a surface in a given direction requires the radiances of other surfaces visible from this point, which generates new visibility and shading problems to solve (?gure 2.9). To calculate those radiances, other surfaces should be evaluated, and our original point on the given surface might contribute to those radiances. As a consequence of that, the rendering equation has complicated coupling between its left and right sides, making the evaluation dif?cult. There are three general and widely accepted approaches to attack this coupling: 1. Local illumination methods Local illumination methods take a drastic approach and simplify or disregard completely all the terms causing problems. The unknown radiance inside the integral of the rendering equation is replaced by some approximation of the emission function. Formally, these methods evaluate the following simpli?ed rendering equation to obtain the radiance:

L(~ ; !) = Le (~ ; !) + x x

Z




Llightsource(h(~ ; !0); !0 )  fr (!0 ; ~ ; !)  cos 0 d!0 ; x x

(2.50)

2.7. CLASSIFICATION OF THE SOLUTION TECHNIQUES

16

where Llightsource may be a simpli?cation of the emission function Le . Abstract lightsources, such as point or directional lightsources are preferred here, since their radiance is a Dirac-delta like function, which simpli?es the integral of equation (2.50) to a sum. These methods take into account only a single re?ection of the light coming from the abstract lightsources. Ideal mirrors and refracting objects cannot be rendered with these methods. 2. Recursive ray-tracing Another alternative is to eliminate from the rendering equation those energy contributions which cause the dif?culties, and thus give ourselves a simpler problem to solve. For example, if limited level, say n, coupling caused by ideal re?ection and refraction were allowed, and we were to ignore the other non-ideal components coming from non-abstract lightsources, then the number of surface points which would need to be evaluated to calculate a pixel color can be kept under control. Since the illumination formula contains two terms regarding the coherent components (re?ective and refracting lights), the maximum number of surfaces involved in the color calculation of a pixel is two to the power of the given depth, i.e. 2n . An implementation of this approach is called visibility ray-tracing which uses the following illumination formula:

L(~ ; !) = Le (~ ; !) + Llightsource(h(~ ; !0 ); !0 )  fr (!0 ; ~ ; !)  cos 0 d!0 + x x x x kr (!r ; ~ ; !)  L(h(~ ; !r ); !r ) + kt (!t ; ~ ; !)  L(h(~ ; !t); !t ); x x x x (2.51) where !r and !t are the ideal re?ection and refraction directions, and kr and kt are the integrated BRDF
components representing the ideal re?ection and refraction, respectively. 3. Global illumination solution Local illumination and recursive ray-tracing apply brutal simpli?cations and thus provide physically inaccurate images. Thus these methods cannot be reliably used in engineering applications, such as in lighting design. These application areas require the rendering equation to be solved without relevant simpli?cations, that leads to the family of global illumination methods. Global illumination methods rely on numerical techniques to resolve the coupling in the rendering or potential equations and to solve the integral equation without unacceptable simpli?cations. The three different approaches represent three levels of the compromise between image generation speed and quality. By ignoring more and more terms in the illumination formula, its calculation can be speeded up, but the result inevitably becomes more and more arti?cial. The ignoration of the terms, however, violates the energy equilibrium, and causes portions of objects to come out extremely dark, sometimes unexpectedly so. These artifacts can be reduced by reintroducing the ignored terms in simpli?ed form, called ambient light. The ambient light represents the ignored energy contribution in such a way as to satisfy the energy equilibrium. Since this ignored part is not calculated, nothing can be said of its positional and directional variation, hence it is supposed to be constant in all directions and everywhere in the 3D space. From this aspect, the role of ambient light also shows the quality of the shading algorithm. The more important a role it has, the poorer quality picture it will generate. In order to formally express the capabilities of a certain algorithm, Heckbert has developed a notation that is based on the regular expressions originally proposed for the syntax of formal languages [Hec91]. The elements of the notation are as follows:

Z




A global illumination algorithm is expected to model all types of lightpaths, that is, it must have L[DjS ] E type. Visibility ray-tracing allows multiple steps only for the ideal re?ection or refraction, thus it can model only L[DS ]E path. Finally, local illumination models simulate only a single non-ideal re?ection and fall into the L[D]E category.

      

E is the eye, L is the lightsource, D is a non-ideal — i.e. incoherent — re?ection or refraction, S is an ideal re?ection or refraction,  is the sign of iteration, that is, it can mean 0; 1; 2; : : : applications of the given operator,
[ ] represents optionality,

j means selection.

2.7. CLASSIFICATION OF THE SOLUTION TECHNIQUES

17

local illumination method

local illumination method with shadow computation

recursive ray-tracing

global illumination method

Figure 2.10: Comparison of local illumination method, recursive ray-tracing and global illumination method. Ambient light was also added when the local illumination method and recursive ray-tracing were computed. The images have been generated by a path tracing program [SK99c]. Note that local illumination methods cannot render ideal re?ection or refraction, thus there are no transparent and mirror like objects. Recursive ray-tracing, on the other hand, is unable to follow multiple non-ideal re?ections, thus the illumination diffusely re?ected from the back wall is missing on the spheres. The running times are 90 seconds, 95 seconds, 135 seconds, 9 hours, respectively, which clearly shows the price we should pay for a physically accurate simulation.

Chapter 3

Optical material models
In order to render realistic images, we have to use realistic material models, also called BRDFs. A realistic BRDF model is expected not to violate physics laws including the Helmholtz-symmetry, and the law of energy conservation, and to mimic the features of real materials. The Helmholtz principle states that the incoming and outgoing directions can be exchanged in the BRDF:

fr (!; ~ ; !0 ) = fr (!0 ; ~ ; !): x x

(3.1)

According to energy conservation, a surface — provided that it is not a lightsource — cannot re?ect more photons (more power) than what have been received, thus the albedo de?ned by

a(~ ; !0 ) = x

Z


H

fr (!0 ; ~ ; !)  cos  d! x

cannot be greater than 1. For isotropic models, re?ections of illuminations from incident directions having the same 0 have rotational symmetry, thus the albedo of a given point depends only on incident angle 0 . To emphasize this, the albedo will be denoted by a(0 ). ~ When introducing the BRDF models, the following notations are used: N is the unit normal vector of the ~ ~ ~ surface, L is a unit vector pointing towards the lightsource, V is a unit vector pointing towards the camera, R is ~ ~ ~ ~ ~ ~ the mirror vector of L onto N , H is a unit vector that is halfway between L and V , 0 is the angle between L and ~ ~ ~ ~ ~ N ,  is the angle between V and N , and is the angle between R and V .

3.1 Diffuse re?ection
First of all, consider diffuse — optically very rough — surfaces re?ecting a portion of the incoming light with radiance uniformly distributed in all directions. Looking at the wall, sand, etc. the perception is the same regardless of the viewing direction.

V

N

L

in

θ θ’

L

Figure 3.1: Diffuse re?ection

If the BRDF is independent of the viewing direction, it must also be independent of the light direction because of the Helmholtz-symmetry, thus the BRDF of these diffuse surfaces is constant:

~ ~ fr;di use(L; V ) = kd :
18

(3.2)

3.2. IDEAL, MIRROR-LIKE REFLECTION

19

The albedo of diffuse re?ection is

adi use(0 ) =
Thus for energy conservation should hold.

Z


H

kd  cos  d! =

= Z2 Z 2 =0 =0

kd  cos   sin  dd = kd  :

(3.3)

kd  1=

3.2 Ideal, mirror-like re?ection

~ An ideal mirror complies the re?ection law of geometric optics, which states that incoming direction L, surface ~ ~ normal N and outgoing direction R are in the same plane, and incoming angle 0 equals to outgoing angle  ~ (0 = ). An ideal mirror re?ects only into the ideal re?ection direction R. The BRDF can thus be de?ned by a Dirac-delta function: ~ ~ ~ ~ fr;mirror(L; V ) = kr (0 )  (R V ) : cos 0
(3.4) The albedo of ideal re?ection is

amirror(0 ) =

Z

H

~ ~ kr (0 )  (R V )  cos  d! = kr (0 ): cos 0

(3.5)

Thus for energy conservation kr cannot exceed 1. Even perfect mirrors absorb some portion of the incident light, as is described by the Fresnel equations of physical optics, expressing the re?ection (F ) of a perfectly smooth mirror in terms of the refractive index of the material n, the extinction coef?cient  which describes the conductivity of the material (for nonmetals  = 0), and the angle of the incidence of the light beam. Denoting the incident angle by 0 and the angle of refraction by , the
2.5 aluminium ezust arany rez 2 7 9 aluminium ezust arany rez

8

1.5 kappa n

6

5

1

4

3 0.5 2

0 400

450

500

550 lambda

600

650

700

750

1 400

450

500

550 lambda

600

650

700

750

Figure 3.2: Refraction index of the gold, copper and silver

Fresnel equation expressing the ratio of the energy of the re?ected beam and the energy of the incident beam for the directions parallel and perpendicular to the electric ?eld is:

where | = 1. These equations can be derived from Maxwell’s fundamental formulae describing the basic laws ~ ~ of electric waves. If the light is unpolarized, that is, the parallel (Ek ) and the perpendicular (E? ) electric ?elds have the same amplitude, then the total re?ectivity is:

p

0 2 cos 0 (n |) 2 F? (; 0 ) = cos  + (n + |)  cos 0 ; Fk (; 0 ) = cos 0 + (n + |)  cos  ; cos  (n + |) cos  + cos 









(3.6)

1=2 ~ 1=2 ~ 2 0 ) = F (; 0 ) = jFk  Ek + F?  E? j = Fk + F? : kr ( ~ ~ 2 jEk + E? j2

(3.7)

3.2. IDEAL, MIRROR-LIKE REFLECTION

20

Note that F is wavelength dependent, since n and  are functions of the wavelength. Parameters n and  are often not available, so they should be estimated from measurements, or from the value of the normal re?ection if the extinction is small. At normal incidence (0 = 0), the re?ection is:

2 (n 1)2 + 2 1 2 F0 () = 1 + (n + |) = (n + 1)2 + 2  n + 1 : 1 (n + |) n
n()  1 + pF0 () : 1 F0 ()









(3.8)

Solving for n gives the following equation:

p

(3.9)

F0 can easily be measured, thus this simple formula is used to compute the values of the index of refraction n. Values of n() can then be substituted into the Fresnel equations to obtain re?ection parameter F for other angles of incidence.
"alufresnel" "goldfresnel"

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 700 650 0 10 600 20 30 550 40 50 500 60 70 450 80 90 400 0 750

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 700 650 10 600 20 30 550 40 50 500 60 70 450 80 90 400 750

Aluminum

Gold

"copperfresnel"

"silverfresnel"

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 700 650 0 10 600 20 30 550 40 50 500 60 70 450 80 90 400 0 750

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 700 650 10 600 20 30 550 40 50 500 60 70 450 80 90 400 750

Copper

Silver

Figure 3.3: The Fresnel function of aluminum, gold, copper and silver for different wavelengths (400..700 nm) and incident angles (0..90 degrees)

3.3. IDEAL REFRACTION

21

3.3 Ideal refraction
When ideal refraction happens, the light goes into the material according to the Snellius-Descartes-law, which ~ ~ ~ states that the incident direction (L), surface normal (N ) and the refraction direction (V ) are all in the same plane, and 0 where n is the relative refraction index of the material. Thus the BRDF of ideal refraction is also a Dirac-delta like function:

sin  = n; sin 

~ where T is the ideal refraction direction.

~ ~ ~ ~ fr;refractor(L; V ) = kt  (T V ) ; cos 0

(3.10)

3.4 Non-ideal, specular re?ection
Incoherent — i.e. not ideal, mirror-like — re?ection is usually broken down into diffuse and specular components. The BRDF of the diffuse re?ection has already been discussed. For the specular re?ection, new models are proposed.

3.4.1 Phong re?ection model and its modi?cations
The Phong BRDF [Pho75] was the ?rst model proposed for specular materials. Specular surfaces re?ect most of the incoming light around the ideal re?ection direction, thus the BRDF should be maximum at this direction and should decrease sharply.

H

N L

in

R ψ V
Figure 3.4: Specular re?ection

L

Phong proposed the following function for this purpose:

~ ~ ~ ~ fr;Phong (L; V ) = ks  cos 0 = ks  (R~ V~) (3.11) cos (N  L) ~ ~ where R is the mirror direction of L onto the surface normal. The ks factor is proportional to the Fresnel function, but is less than that, because now the surface is not an ideal mirror. Factor ks can be assumed to be independent of the wavelength and the incident angle for non-conducting
n n
materials (the highlight of a white lightsource on plastics is also white). The original Phong model is not physically plausible since it is not symmetric. Due to this, in photorealistic image synthesis the following reciprocal Phong BRDF is preferred [ICG86]:

~ ~ ~ ~ fr;reciprocalPhong(L; V ) = ks  cosn = ks  (R  V )n

(3.12)

To compute the albedo of the reciprocal Phong model, an appropriate parameterization of the directional sphere ~ is needed (?gure 3.5). Now the north pole of the spherical coordinate system is the ideal re?ection direction R. Let ~ , and by another angle  between its projection onto a plane perpendicular us identify a direction by an angle from R ~ to R and an arbitrary vector on this plane. When the viewing direction is parameterized in this way, the vertical angle is just .

3.4. NON-IDEAL, SPECULAR REFLECTION

22

N R V plane perpendicular to R

ψ

φ

surface

reference direction on the plane perpendicular to R
Figure 3.5: Parameterization for the calculation of the albedo

Unfortunately, the domain of ( ; ) is rather complicated since those pairs that would identify a direction pointing into the object should be excluded from re?ection directions. Without this exclusion the albedo is overestimated. Note that this overestimation can only double the albedo at the worst case which happens at grazing angles when half of the re?ection lobe is cut away by the surface. Representing this effect by a function Cut(0 ), the integral computing the albedo is:

areciprocalPhong(0 ) = Cut(0 ) 

= Z2 Z 2 =0 =0

ks  cosn  cos   sin d d:

Function Cut(0 ) is 1 at perpendicular illumination and decreases to a half at horizontal illumination. When the illumination is perpendicular, then = , thus the albedo for this case is:

areciprocalPhong(0) =

= Z2 Z 2

=0 =0

ks  cosn

 cos  sin d d = 2ks 
ks  n2+ 2 : 

= Z2
=0

cosn+1  sin d = ks  n2 2 : +

Since the albedo is less than this for other illumination angles, energy conservation requires the following relation to hold [LW94]: (3.13) Considering the estimation of the albedo for other illumination angles, we have to realize that even if function

Cut(0 ) is approximated by its maximum, the cos  factor still forbids the symbolic integration. If the surface is very shiny — n is large — then the re?ection lobe is narrow, thus cos  is approximately constant and is equal to cos 0 where cosn is not negligible, thus we can further obtain:

areciprocalPhong(0 )  cos 0 

= Z2 Z 2

=0 =0

 ks  cosn  sin d d = ks  n2+ 1  cos 0 :

The re?ected radiance of this model goes to zero for larger incident angles, which is not realistic when rendering metals (left of ?gure 3.6). This artifact is eliminated in the following max-Phong BRDF [NNSK98b]:

we obtain

~ ~ (R  V )n ~  V ); (N  L)) max ((N ~ ~ ~ The albedo of this max-Phong model is maximum at perpendicular illumination when 0 = 0 and  =
cos ~ ~ fr;maxPhong(L; V ) = ks  max (cos ; cos 0 ) = ks 
= Z2 Z 2 n

(3.14) , thus

amaxPhong (0) =

=0 =0

 ks  cosn  sin  maxcos ; 1) d d = ks  n2+ 2 : (cos

3.4. NON-IDEAL, SPECULAR REFLECTION

23

Figure 3.6: Comparison of the reciprocal Phong (left) and the max-Phong (right) models (n = 20 in the upper row, and n = 5000 in the lower row). The reciprocal Phong model becomes dark for grazing angles.

Figure 3.7: Metallic objects de?ned by the max-Phong model

3.4. NON-IDEAL, SPECULAR REFLECTION

24

As for the reciprocal Phong model, energy conservation requires that

ks  n2+ 2 : 
At non-perpendicular illumination angles, only the cutting of the re?ection lobe attenuates signi?cantly the albedo, for shiny materials cos =max (cos ; cos 0 )  1 in the re?ection lobe. Let us denote the factor representing the cutting of the re?ection lobe by Cut (0 ). Note that Cut (0 ) is similar, but not identical to Cut(0 ).

a(0 ) = Cut (0 )

= Z2 Z 2

=0 =0

cos   ks  cosn  sin  max (cos ; cos 0 ) d d  Cut (0 )  n2+ 1 :

If parameter ks is de?ned from the Fresnel function, we have to specify which incident angle must be used. The angle between the surface normal and the light vector is not appropriate, because the BRDF will not be symmetric and the real normal vector where the re?ection actually happens is not equal to the normal of the mean surface due to the surface irregularities (?gure 3.8). This real normal vector is, in fact, a random variable. If the surface is a ~ ~ collection of randomly oriented ideal mirrors, then those surface elements which re?ect from L to V has normal ~ ~ ~ ~ ~ vector H = (L + V )=2. Thus the cosine of the incoming angle is computed from the (H  L) scalar product for the Fresnel function.

3.4.2 Cook-Torrance model
Specular re?ection can be more rigorously analyzed by modeling the surface irregularities by probability distributions, as has been proposed by Torrance, Sparrow, Cook and Blinn. In their model, the surface is assumed to consist of randomly oriented perfect mirrors, so-called microfacets. As in the previous section, the re?ected light is broken down into diffuse and specular components. The diffuse component is believed to be generated by multiple re?ections on the microfacets and also by emission of the absorbed light by the material of the surface. The diffuse component is well described by Lambert’s law. The specular component, on the other hand, is produced by the direct re?ections of the microfacets. This section discusses the derivation of the bi-directional transfer proba~ ~ bility density w(L; V ) for this specular part. Recall that the transfer function is the following probability density (equation (2.8)):

~ ~ ~ ~ w(L; V )  d! = Prfphoton is re?ected directly to d! around V j coming from Lg: (3.15) ~ ~ Concerning this type of re?ection from direction L to d! around direction V , only those facets can contribute ~ . If re?ection is to happen, the facet must obviously be whose normal is in d!H around the halfway unit vector H

facing in the right direction. It should not be hidden by other facets, nor should its re?ection run into other facets, and it should not absorb the photon for the possible contribution.

N

Figure 3.8: Microfacet model of the re?ecting surface

~ Considering these facts, the event that “a photon is re?ected directly to d! around V ” can be expressed as the logical AND connection of the following events: ~ 1. Orientation: In the path of the photon there is a microfacet having its normal in d!H around H .
2. No shadowing or masking: The given microfacet is not hidden by other microfacets from the photon coming from the lightsource, and the re?ected photon does not run into another microfacet. 3. Re?ection: The photon is not absorbed by the microfacet that is supposed to be a perfect mirror.

3.4. NON-IDEAL, SPECULAR REFLECTION

25

Since these events are believed to be stochastically independent, their probability can be calculated independently, and the probability of the composed event will be their product. Concerning the probability of the microfacet normal being in d!H , we can suppose that all facets have equal ~ area f . Let the probability density of the microfacet normal be P (H ). Blinn [Bli77] proposed Gaussian distribu~ tion for P (H ), since it seemed reasonable due to the central limit value theorem of probability theory:

~ where is the angle of the microfacet with respect to the normal of the mean surface, that is the angle between N ~ , and m is the root mean square of the slope, i.e. a measure of the roughness. and H Later Torrance and Sparrow showed that the results of the early work of Beckmann [BS63] and Davies [Dav54], who discussed the scattering of electromagnetic waves theoretically, can also be used here and thus Torrance proposed the Beckmann distribution function instead of the Gaussian: ~ Let us examine a surface element dA. If a photon arrives from direction L to a surface element dA, the visible ~ ~ area of the surface element will be dA  (N  L), while the total visible area of the microfacets having their normal ~ ~ ~ in the direction around H will be the product of the visible size of a single microfacet (f  (H  L)), and the number of properly oriented microfacets. The number of properly oriented microfacets, in turn, is the product ~ of the probability that a microfacet is properly oriented (P (H )  d!H ) and the number of all microfacets in dA (dA=f ). The probability of ?nding an appropriate microfacet aligned with the photon is the visible size of properly oriented microfacets divided by the visible surface area, thus we obtain: ~ ~ ~ d! ~ ~ ~ Prforientationg = f  (H  L)  P (H )  ~ H  dA=f = P (H )  d!H~ (H  L) : ~ ~ dA  (N  L) (N  L)
L H dω H dHvert dH hor dφ dω dVvert V dV hor
Figure 3.9: Calculation of d!H =d!

~ P (H ) = const  e ( =m)2 ;

(3.16)

1 ~ P (H ) = m2 cos4  e

tan2 m2



:

(3.17)

(3.18)

θH θV

Unfortunately, in the BRDF the transfer probability density uses d! instead of d!H , thus we have to determine ~ De?ning a spherical coordinate system (; ), with the north pole in the direction of L (?gure 3.9), the solid angles are expressed by the product of vertical and horizontal arcs:

d!H =d! [JGMHe88].

d! = dVhor  dVvert ; d!H = dHhor  dHvert:
By using geometric considerations and applying the law of re?ection, we get:

(3.19)

dVhor = d  sin V ; dHhor = d  sin H ; dVvert = 2dHvert :
This in turn yields:

(3.20)

d!H = sin H = sin H = 1 = 1 : ~ ~ d! 2 sin V 2 sin 2H 4 cos H 4(L  H )

(3.21)

3.4. NON-IDEAL, SPECULAR REFLECTION

26

since V

= 2  H . Thus the orientation probability is

~) ~ ~ ~ Prforientationg = P (H )~ (H  L)  d!H  d! = P~(H~  d!: ~ d! (N  L) 4(N  L)

(3.22)

β+2α?π/2

H β V

N α

L

l2 π/2?β

l1

l2

Figure 3.10: Geometry of masking

~ The visibility of the microfacets from direction V means that the re?ected photon does not run into another microfacet. The collision is often referred to as masking. Looking at ?gure 3.10, we can easily recognize that the probability of masking is l1 =l2 , where l2 is the one-dimensional length of the microfacet, and l1 describes the boundary case from where the beam is masked. The angles of the triangle formed by the bottom of the ~ ~ microfacet wedge and the beam in the boundary case can be expressed by the angles = angle(N; H ) and ~ ~ ~ ~ = angle(V ; H ) = angle(L; H ) by geometric considerations and by using the law of re?ection. Applying the sine law for this triangle, and some trigonometric formulae: ~ ~ ~ ~ ~ ~ According to the de?nitions of the angles cos = N  H , cos( + ) = N  V and cos = V  H . If the angle of incident light and the facet normal do not allow the triangle to be formed, the probability of no masking taking place is obviously 1. This situation can be recognized by evaluating the formula without any previous considerations and checking whether the result is greater than 1, then limiting the result to 1. The ?nal result is: ~ ~ The probability of shadowing can be derived in exactly the same way, only L should be substituted for V :
( Prfnot shadowingg = minf2  (N  H )  ~N  L) ; 1g: ~  H) (L ( ~ ~ ~ Prfno shadow and maskg  minf2  (N  H )  (N  V ) ; 2  (N  H )  ~N  L) ; 1g = G(N; L; V ): ~  H) ~ ~  H) (V (L Prfnot maskingg = minf2  (N  H )  (N  V ) ; 1g: ~ ~ (V  H )

= cos( Prfnot maskingg = 1 l1 = 1 sin( + 22 ) 2) = 2  cos cos + ) : l2 sin(=

(3.23)

~ ~

~ ~

(3.24)

~ ~

~ ~

(3.25)

The probability of neither shadowing nor masking taking place can be approximated by the minimum of the two probabilities:

~ ~

~ ~

~ ~

~ ~

(3.26)

~ where variable 0 has been replaced by H ~. equal to H

Note that this approximation upperbounds the real probability, particularly at grazing angles, which can result in the violation of the energy balance. When analyzing the perfect mirrors, we concluded that even a perfect mirror absorb some portion of the incident light, as is described by the Fresnel equations. Since F (; 0 ) is the fraction of the re?ected energy, it also describes the probability of a photon being re?ected at a microfacet that is assumed to be an ideal mirror, giving: ~ ~ Prfre?ectiong = F (; H  L)

~ ~ ~  L since those microfacets that re?ect from L to V

have normal vector

3.4. NON-IDEAL, SPECULAR REFLECTION

27

~ ~ w(L; V )d!:

Now we can summarize the results by multiplying the probabilities of the independent events to express

~ ~ w(L; V )d! = Prforientationg  Prfno mask and shadowg  Prfre?ectiong = ~ P (H )  G(N; L; V )  F (; H  L) d!: ~ ~ ~ ~ ~ (3.27) ~ ~ 4(N  L) ~ ~ The BRDF is the transfer probability divided by d! and the cosine of the outgoing angle (N  V ) , thus we obtain ~) ~ ~ ~ ~ ~ ~ ~ fr;Cook Torrance(L; V ) = ~ P~(H~ ~  G(N; L; V )  F (; H  L): (3.28) 4(N  L)(N  V )

Chapter 4

Solution strategies for the global illumination problem
Since the rendering or the potential equations contain the unknown radiance function both inside and outside the integral, in order to express the solution, this coupling should be resolved. The possible solution techniques fall into one of the following three categories: inversion, expansion and iteration. Operator T represents light-surface interaction, thus each of its application generates a higher-bounce estimate of the light transport (or alternatively T 0 represents potential-surface interaction). For physically plausible optical material models, a re?ection or refraction always decreases the total energy, thus the integral operator is always a contraction. However, when the transport is evaluated numerically, computation errors may pose instability problems if the scene is highly re?ective. As we shall see, expansion and iteration exploit the contractive property of the transport operator, but inversion does not.

4.1 Inversion
Inversion groups the terms that contain the unknown function on the same side of the equation and applies formally an inversion operation: (1 T )L = Le =) L = (1 T ) 1 Le : (4.1) Thus the measured power is

However, since T is in?nite dimensional, it cannot be inverted in closed form. Thus it should be approximated by a ?nite dimensional mapping, that is usually given as a matrix. This kind of approximation is provided by ?nite-element methods that project the problem into a ?nite dimensional function space, and approximate the solution here. This projection converts the original integral equation to a system of linear equations, which can be inverted, for example, by Gaussian elimination method. This approach was used in early radiosity methods, but have been dropped due to the cubic time complexity and the numerical instability of the matrix inversion. Inversion has a unique property that is missing in the other two methods. Its ef?ciency does not depend on the contractivity of the integral operator, neither does it even require the integral operator to be a contraction.

ML = M(1 T ) 1 Le :

(4.2)

4.2 Expansion
Expansion techniques eliminate the coupling by obtaining the solution in the form of an in?nite Neumann series.

4.2.1 Expansion of the rendering equation: gathering walks Substituting the right side’s L by Le + T L, which is obviously L according to the equation, we get:

L = Le + T L = Le + T (Le + T L) = Le + T Le + T 2 L:
Repeating this step n times, the original equation can be expanded into a Neumann series:

(4.3)

L=

n X i=0

T i Le + T n+1 L:
28

(4.4)

4.2. EXPANSION

29

If integral operator T is a contraction, then limn!1 T n+1 L = 0, thus

L=
The measured power is

1 X
i=0

T i Le : MT i Le :

(4.5)

ML =

1 X
i=0

(4.6)

The terms of this in?nite Neumann series have intuitive meaning as well: MT 0 Le emission, MT 1 Le comes from a single re?ection, MT 2 Le from two re?ections, etc.

= MLe comes from the

x2 ’ θ2 L(x, ωp ) ωp ’ θ1 x1
Figure 4.1: The integrand of

ω’ 2

’ ω1 x3

p

MT 2 Le is a two-step gathering walk

In order to understand how this can be used to determine the power going through a single pixel, let us examine the structure of MT i Le as a single multi-dimensional integral for the i = 2 case:

M(T 2 Le ) =

0 0 where Sp is the pixel area, c(p) is the camera parameter, p is a running point on this pixel, !1 and !2 are the ~ ~ directions of the ?rst and second re?ections, respectively (?gure 4.1) and

Sp
01
02

Z Z Z c(p) ~  w (~ )  w (~ )  Le (~ ; !0 ) d!0 d!0 d~: x3 2 2 1 p Sp 1 x1 2 x2

(4.7)

~ 1 = h(p; !p); x ~ ~ 0 ~ 2 = h(~ 1 ; !1 ); x x 0 0 0 ~ 3 = h(~ 2 ; !2 ) = h(h(~ 1 ; !1); !2 ); x x x
and the weights are

(4.8)

0 x ~ 0 w1 = fr (!1 ; ~ 1 ; !p )  cos 1 ; 0 ; ~ 2 ; !0 )  cos 0 : w2 = fr (!2 x 1 (4.9) 2 Thus to evaluate the integrand at point (p; !1 ; !2 ), the following algorithm must be executed: ~ 0 0 1. Point ~ 1 = h(p; !p ) that is visible through the point ~ of the pixel from the eye should be found. This x ~ ~ p can be done by sending a ray from the eye into the direction of p and identifying the surface that is ?rst ~
intersected.

0 0 2. Surface point ~ 2 = h(~ 1 ; !1 ) — that is the point which is visible from ~ 1 at direction !1 — must be x x x 0 and identifying the surface determined. This can be done by sending another ray from ~ 1 into direction !1 x that is ?rst intersected. 0 0 3. Surface point ~ 3 = h(h(~ 1 ; !1 ); !2 ) — that is the point visible from ~ 2 at direction x x x 0 This means the continuation of the ray tracing at direction !2 . 0 !2 — is identi?ed.

0 4. The emission intensity of the surface at ~ 3 in the direction of !2 is obtained and is multiplied with the cosine x terms and the BRDFs of the two re?ections.

4.2. EXPANSION

30

This algorithm can easily be generalized for arbitrary number of re?ections. A ray is emanated recursively 0 0 0 from the visible point at direction !1 then from the found surface at !2 , etc. until !n . The emission intensity at the end of the walk is read and multiplied by the BRDFs and the cosine terms of the stages of the walk. These 0 walks provide the value of the integrand at “point” p; !1 ; !2 ; : : : ; !n . Note that a single walk of length n can be ~ 0 0 used to estimate the 1-bounce, 2-bounce, etc. n-bounce transfer simultaneously, if the emission is transferred not only from the last visited point but from all visited points. The presented walking technique starts at the eye and gathers the illumination encountered during the walk. The gathered illumination is attenuated according to the cosine weighted BRDFs of the path. So far, we have examined the structure of the terms of the Neumann series as a single multi-dimensional integral. Alternatively, it can also be considered as recursively evaluating many directional integrals. For example, the two-bounce transfer is:

MT 2 Le =

Z

~ ~ In order to estimate the outer integral of variable p, the integrand is needed in sample point p. This, in turn, 0 p requires the estimation of the integral of variable !1 at ~, which recursively needs again the approximation of the 0 ~ 0 integral of variable !2 at (p; !1 ). If the same number — say m — of sample points are used for each integral quadrature, then this recursive approach will use m points for the 1-bounce transfer, m2 for the two-bounces, m3 for the three-bounces, etc. This kind of subdivision of paths is called splitting [AK90]. Splitting becomes prohibitive for high-order re?ections and is not even worth doing because of the contractive property of the integral operator. Due to the contraction, the contribution of higher-order bounces is less, thus it is not very wise to compute them as accurately as low-order bounces.
4.2.2 Expansion of the potential equation: shooting walks
The potential equation can also be expanded into a Neumann series similarly to the rendering equation:

Sp

2 2 3 3 Z 6Z 05 05 p w0  6 w1  4 w2  Le d!2 7 d!1 7 d~: 4

01
02

(4.10)

W=
which results in the following measured power

1 X
i=0

T 0i W e ; M0 T 0i W e :

(4.11)

M0 W =

1 X
i=0

(4.12)

M0 W e is the power measured by the device from direct emission, M0 T 0 W e is the power after a single re?ection, M0 T 02 W e is after two re?ections, etc.
y2 θ2 Φ (dy1 , d ω1 ) p ωp θ3 y3
Figure 4.2: The integrand of

ω1 θ1 y1

ω2

M0 T 02 W e is a two-step shooting walk

Let us again consider the structure of M0 T 02 W e :

M0 T 02 W e =

ZZZZ

S
1
2
3

Le (~1 ; !1 )  w0 (~1 )  w1 (~2 )  w2 (~3 )  W e (~3 ; !3 ) d!3 d!2 d!1 d~1 ; y y y y y y

(4.13)

4.3. ITERATION

31

where

~2 = h(~1 ; !1 ); y y ~3 = h(~2 ; !2 ) = h(h(~1 ; !1 ); !2 ) y y y
and the weights are

(4.14)

w0 = cos 1 ; w1 = fr (!1 ; ~2 ; !2 )  cos 2 ; y w2 = fr (!2 ; ~3 ; !3 )  cos 3 : y

(4.15)

Substituting the measuring function of the pin-hole camera (equation (2.45)), we can conclude that M0 T 02 W e can be non-zero only if ~3 is visible through the given pixel, and the integral over
3 is simpli?ed to a single value by y the Dirac-delta function, thus the ?nal form of the two-bounce re?ection is:

ZZZ

S
1
2

Le (~1 ; !1 )  w0 (~1 )  w1 (~2 )  w2 (~3 )  g(~3 ) d!2 d!1 d~1 : y y y y y y

(4.16)

if ~3 is visible though the given pixel and 0 otherwise, where g (~3 ) is the surface dependent camera parameter. y y Thus to evaluate the integrand at point (~1 ; !1 ; !2 ), the following algorithm must be executed: y

y y y 1. The cosine weighted emission of point ~1 in direction !1 is computed. Surface point ~2 = h(~1 ; !1 ) — that is the point which is visible from ~1 at direction !1 — must be determined. This can be done by sending a y ray from ~1 into direction !1 and identifying the surface that is ?rst intersected. This point “receives” the y computed cosine weighted emission. y y y 2. Surface point ~3 = h(h(~1 ; !1 ); !2 ) — that is the point visible from ~2 at direction !2 — is identi?ed. This means the continuation of the ray tracing at direction !2 . The emission is weighted again by the local BRDF and the cosine of the normal and the outgoing direction. y 3. It is determined whether or not this point ~3 is visible from the eye, and through which pixel. The contribution to this pixel is obtained by weighting the transferred emission by the local BRDF, cosine angle and the surface dependent camera parameter.
This type of walk, called shooting, starts at a known point ~1 of a lightsource and simulates the photon re?ecy tion for a few times and ?nally arrives at a pixel whose radiance this walk contributes to. Note that in gathering walks the BRDF is multiplied with the cosine of the angle between the normal and the incoming direction, while in shooting walks with the cosine of the angle between the normal and the outgoing direction.

4.2.3 Merits and disadvantages of expansion methods
The main problem of expansion techniques is that they require the evaluation of very high dimensional integrals that appear as terms in the in?nite series. Practical implementations usually truncate the in?nite Neumann series, which introduces some bias, or stop the walks randomly, which signi?cantly reduces the samples of higher order interre?ections. These can result in visible artifacts for highly re?ective scenes. On the other hand, expansion methods also have an important advantage. Namely, they do not require temporary representations of the complete radiance function, thus do not necessitate ?nite-element approximations. Consequently, these algorithms can work with the original geometry without tessellating the surfaces to planar polygons. Expansion techniques generate random walks independently. It can be an advantage, since these algorithms are suitable for parallel computing. However, it also means that these methods “forget” the previous history of walks, and they cannot reuse the visibility information gathered when computing the previous walks, thus they are not as fast as they could be.

4.3 Iteration
Iteration techniques realize that the solution of the rendering equation is the ?xed point of the following iteration scheme Ln = Le + T Ln 1 ; (4.17)

4.3. ITERATION

32

thus if operator T is a contraction, then this scheme will converge to the solution from any initial function L0 . The measured power can be obtained as a limiting value

ML = nlim MLn: !1

(4.18)

In order to store the approximating functions Ln , usually ?nite-element methods are applied, as for example, in diffuse radiosity [SC94], or in non-diffuse radiosity using partitioned hemisphere [ICG86], directional distributions [SAWG91] or illumination networks [BF89]. There are two critical problems here. On the one hand, since the domain of Ln is 4 dimensional and Ln has usually high variation, an accurate ?nite-element approximation requires very many basis functions, which, in turn, need a lot of storage space. Although hierarchical methods [HSA91, AH93], wavelet or multiresolution methods [CSSD96, SGCH94] and clustering [SDS95, CLSS97, SPS98] can help, the memory requirements are still prohibitive for complex scenes. This problem is less painful for the diffuse case since here the domain of the radiance is only 2 dimensional. On the other hand, when ?nite element techniques are applied, operator T is only approximated, which introduces some non-negligible error in each step. If the contraction ratio of the operator is , then the total accumulated error will be approximately 1=(1 ) times the error of a single step [SKFNC97]. For highly re?ective scenes — i.e. when  is close to 1 — the iteration is slow and the result is inaccurate if the approximation of the operator is not very precise. Very accurate approximations of the transport operator, however, require a lot of computation time and storage space. In the diffuse radiosity setting several methods have been proposed to improve the quality of the iteration. For example, we can use Chebyshev iteration instead of the Jacobi or the Gauss-Seidel method for such ill conditioned systems [BBS96]. On the other hand, realizing that the crucial part of designing such an the algorithm is ?nding a good and “small” approximation of the transport operator, the method called well-distributed ray-sets [NNB97, BNN+ 98] proposes the adaptive approximation of the transport operator. This approximation is a set of rays selected carefully taking into account the important patches and directions. In [BNN+ 98], the adaptation of the discrete transport operator is extended to include patch subdivision as well, to incorporate the concepts of hierarchical radiosity [HSA91]. The adaptation strategy is to re?ne the discrete approximation (by adding more rays to the set), when the iteration with the coarse approximation is already stabilized. Since the discrete approximation of the transport operator is not constant but gets ?ner in subsequent phases, the error accumulation problem can be controlled but is not eliminated. This thesis proposes a new method called the stochastic iteration to attack both the problem of prohibitive memory requirements and the problem of error accumulation. Compared to expansion techniques, iteration has both advantages and disadvantages. Its important advantage is that it can potentially reuse all the information gained in previous computation steps and can exploit the coherence of the radiance function, thus iteration is expected to be faster than expansion. Iteration can also be seen as a single in?nite length random walk. If implemented carefully, iteration does not reduce the number of estimates for higher order interre?ections, thus it is more robust when rendering highly re?ective scenes than expansion. The property that iteration requires tessellation and ?nite-element representation is usually considered as a disadvantage. And indeed, sharp shadows and highlights on highly specular materials can be incorrectly rendered and light-leaks may appear, not to mention the unnecessary increase of the complexity of the scene description (think about, for example, the de?nition of an original and a tessellated sphere). However, ?nite-element representation can also provide smoothing during all stages of rendering, which results in more visually pleasing and dot-noise free images. Summarizing, iteration is the better option if the scene is not highly specular.

4.3.1 Analysis of the iteration
In order to ?nd necessary conditions for the convergence, let us examine two subsequent iteration steps:

Ln = Le + T Ln 1; Ln 1 = Le + T Ln 2: Substracting the two equations and assuming that L0 = 0, we can obtain: Ln Ln 1 = T (Ln
If operator T is a contraction, that is if

(4.19)

1

Ln 2 ) = T n 1 (L1 L0 ) = T n 1 Le :

(4.20)

jjT Ljj <   jjLjj;  < 1;

(4.21)

4.4. ANALYTICAL SOLUTION OF THE RENDERING EQUATION

33

with some function norm, then

jjLn Ln 1 jj = jjT n 1 Le jj < n 1  jjLe jj:

(4.22)

Thus iteration converges to the real solution at least with the speed of a geometric series. The contraction ratio  depends on both the geometry and the average re?ectivity of the scene. From the point of view of the geometry, the contraction ratio is maximum if the environment is closed, when  corresponds to the average albedo. Error caused by the approximation of the transport operator In practice operator T cannot be evaluated exactly, which introduces some error in each step of the iteration. The errors may accumulate in the ?nal solution. This section examines this phenomenon. Assume that an operator T  which is only an approximation of T is used in the iteration. Suppose that both operators are contractions, thus both iterations will converge from any initial function. Let the radiance after n iterations of operator T  be Ln and let us assume that the iteration starts at the solution of the exact equation, that is L0 = L (since the iteration converges to the same limit from any initial function, this assumption can be made). The solution of the approximated equation is L = L1 . The error at step n is

jjLn Ljj = jjLn Ln 1 + Ln
Since if i > 1, then

1

::: + L1 Ljj 

n X i=1

jjLi Li 1 jj:

(4.23)

jjLi Li 1 jj = jjT  Li
we have

1

T  Li 2 jj = jjT  (Li

1

Li 2 )jj    jjLi

1

Li 2 jj  i 1  jjL1 L0 jj
(4.24)

n X i=1

jjLi Li 1 jj  jjL1 L0 jj  (1 +  + 2 + : : : n 1 ):

Letting equations

n go to in?nity, we obtain the error between the ?xed points of the original and the approximated
1 L jjL Ljj  jjL1 L0 jj  (1 +  + 2 + : : :) = jjL1  0 jj :
(4.25)

On the other hand, subtracting the real solution de?ned as the ?xed point of the exact equation L from the ?rst iteration of the approximated operator

= Le + T L

L1 = Le + T  L;
we obtain Thus the ?nal error formula is

L1 L = T  L T L:
L jjL Ljj  jjT 1 T Ljj :

(4.26) (4.27)

4.4 Analytical solution of the rendering equation
In order to test the error and convergence of global illumination algorithms, we need scenes for which the exact solution is known. Generally, there are two alternatives. Either we use simple scenes for which analytical solution is possible, or numerical methods are applied using so many samples that the approximate solution can be accepted as a reference. In order to ?nd analytically solvable scenes, we use a reverse approach. Normally, a scene is analyzed to ?nd the radiance distribution. Now, we start with a prescribed radiance distribution and search for scenes where the radiance would be identical to the given radiance. Two special cases are examined. In the ?rst case we are looking for scenes where the radiance is constant everywhere and at every direction. In the second case we establish criteria for making only the re?ected radiance constant.

4.4. ANALYTICAL SOLUTION OF THE RENDERING EQUATION

34

4.4.1 Scenes with constant radiance

~ Suppose that the constant radiance is L and also that the scene is a closed environment, i.e. looking at any direction we can see a surface, thus the incoming radiance is also constant. Substituting this to the rendering equation we get: ~ ~ x ~ ~ x L = Le (~ ; !) + (T L)(~ ; !) = Le (~ ; !) + L  (T 1)(~ ; !) = Le (~ ; !) + L  a(~ ; !); x x x x (4.28)
since the albedo is the response to homogeneous unit illumination. According to this equation, the radiance will be constant for scenes of arbitrary geometry and of arbitrary emission function if the albedo function is:

If Le is constant, then the required albedo will also be constant. In this special case — called the homogeneous diffuse environment — the geometry is arbitrary, but all surfaces have the same diffuse re?ection and emission [Shi91b, Kel95].

ex a(~ ; !) = 1 L (~ ; !) : x ~ L

(4.29)

4.4.2 Scenes with constant re?ected radiance

~ Let us assume that the re?ected radiance — i.e. the radiance without the emission — is Lr and assume again that the scene is closed. Substituting this into the rendering equation, we obtain: ~ ~ x ~ L(~ ; !) = Le (~ ; !) + Lr = Le (~ ; !) + (T (Le + Lr ))(~ ; !) = Le (~ ; !) + (T Le )(~ ; !) + Lr  a(~ ; !): x x x x x x
This imposes the following constraint on the albedo function: (4.30)

e )(x a(~ ; !) = 1 (T L ~ ~ ; !) : x Lr

(4.31)

Note that if the re?ected radiance is constant, then the albedo is also constant, which is guaranteed if the BRDF is diffuse and has the same fr value everywhere. An appropriate geometry, which makes T Le constant for arbitrary diffuse emission if the BRDF is diffuse and constant, is the internal surface of the sphere as proposed originally by [HMF98]. Formally, let the scene be an inner surface S of a sphere of radius R, the BRDF be constant fr = a= and the emission be diffuse and de?ned by Le (~ ). Using the x

y y d!0 = cos ~ ~ jd~ j~ x 2 y

substitution for the solid angle (equation (2.2)), we obtain the following form of the light transport operator:

(T Le )(~ ; !) = x

Z
S

cos fr  cos ~  Le (~)  j~ ~j2 d~: y y ~y y x x

(4.32)

R θx x

R θy y

Figure 4.3: Geometry of the reference scene

Looking at ?gure 4.3, we can see that inside a sphere cos ~ x ?nal form:

(T Le )(~ ; !) = x

Z
S

R Le(~) d~ y y Z x y y e (~ )  cos ~  cos ~ d~ = fr  Le (~ ) d~ = a  S fr  L y y y j~ ~ j2 y x 4R2 4R2  :
S

y x = cos ~ = j~2R~ j ; thus we can obtain the following y
(4.33)

The response is equal to the product of the albedo and the average emission of the total spherical surface.

Chapter 5

Finite-element methods for the Global Illumination Problem
Iteration requires the representation of the temporary radiance function Ln . So does expansion if view-independent solution is needed since the ?nal radiance distribution must be represented in a continuous domain. The exact representation of such functions might need in?nite data, which is not available for practical algorithms. Thus ?nite approximations are required. To represent a function over a continuous domain, ?nite-element methods can be used which approximate the function in a ?nite function series form:

L(~ ; !)  L(n)(~ ; !) = x x

n X j =1

Lj  bj (~ ; !) = bT (~ ; !)  L x x L

(5.1)

where bj (~ ; ! ) is a system of prede?ned basis functions, and j factors are unknown coef?cients. x This representation can also be seen as projecting the in?nite dimensional space of the possible radiance functions into a ?nite-dimensional function space de?ned by the basis functions.
L b2 L
n

L L b1 b1 b2 b3 b4
n

Figure 5.1: Finite element approximation

L + Ln f cosθ’d ω r L
n

e

~ b2 ~ b1

Figure 5.2: Projection to the adjoint base

Substituting this approximation into the rendering equation we can obtain:

n X j =1

Lj  bj (~ ; !)  x

n X j =1

Le  bj (~ ; !) + T x j

n X j =1

Lj  bj (~ ; !) = x
35

n X j =1

Le  bj (~ ; !) + x j

n X j =1

Lj  T bj (~ ; !): x

(5.2)

5. FINITE-ELEMENT METHODS FOR THE GLOBAL ILLUMINATION PROBLEM

36

Note that equality cannot be guaranteed, since even if j =1 j  bj (~ ; ! ) is in the subspace de?ned by the basis x functions, the integral operator T may result in a function that is out of this space. Instead, equality is required in an appropriate subspace de?ned by adjoint basis functions ~1 (~ ); ~2 (~ ); : : : ~n (~ ) (?gure 5.2). This set is required b x b x b x to be orthogonal to the original base b1 (~ ); b2 (~ ); : : : bn (~ ) in the following sense: x x x

Pn L

hbi (~ ); ~j (~ )i = x b x

( 1 if i = j ,

b Having projected equation 5.2 to the adjoint base — i.e. multiplying it by each adjoint basis functions ~i — we obtain the following system of linear equations:

0 otherwise.

(5.3)

X Li = Le + hT bj ; bii  Lj : i
j

(5.4)

This system of linear equations can also be presented in a vector form:

L = Le + R  L; Rij = hT bj ; ~ii: b P

(5.5)

An adjoint of this linear equation can be derived by selecting the adjoint base as the normalization of the measurement functions. Suppose that each basis function bi is associated with a measurement device Wie that measures the power i leaving the support of the basis function, thus the appropriate adjoint basis function is ~i = Wie =hbi ; Wie i. Thus the measured radiance is b

h
Similarly, the measured emission is

n X j =1

Lj  bj ; Wiei = hbi ; Wiei  Li = Pi : Le  bj ; Wiei = hbi ; Wiei  Le = Pe : j i i
n X j =1

h

n X j =1

Applying measurement operator Wie for equation (5.2), we can obtain the following equation:

hbi ; Wie i  Li = hbi ; Wie i  Le + i
Replacing the radiances by measured values, we get

hT bj ; Wie i  Lj :

(5.6)

Pi = Pe + i
This can also be presented in a matrix form

n X hT bj ; Wie i j =1

hbj ; Wje i  Pj :

(5.7)

P = Pe + H  P;
where

(5.8)

W hb Hij = hTbbj; ;W eiii = hT bj ; ~ii  hhbji ;; Wieii = Rij  hbji ;; Wieii : b b W hj W
e e e j j j

(5.9)

When ?nite-element techniques are used together with expansion, ?nite-element representation can either be used to represent the ?nal result [Kel95], or even be involved in the random walk [PM95]. The latter case may correspond either to the random-walk solution of the linear equation derived by projecting the integral equation, or to the Monte-Carlo evaluation of the multi-dimensional integral containing both the transport and the projection operators. The second case is preferred, because it does not require matrix to be explicitly computed and stored. The main problem of ?nite-element representations is that they require a lot of basis functions to accurately approximate high-variation, high-dimensional functions. Not surprisingly, ?nite-element methods have become really popular only for the diffuse case, where the radiance depends on 2 scalars and is relatively smooth. For solving the non-diffuse case, they are good only if the surfaces are not very specular.

R

5.1. GALERKIN’S METHOD

37

5.1 Galerkin’s method
Galerkin’s method ?nds an approximation of the solution by making the error orthogonal to the set of basis functions. It means that the projection of error to the original base is zero. Formally, in Galerkin’s method the set of basis functions is the same — except for normalization constants — as the set of adjoint basis functions. A particularly important case is the piece-wise constant approximation, when the surfaces and the directions are partitioned into patches A1 ; A2 ; : : : ; An and solid angles
1 ;
2 ; : : : ;
m , respectively. The basis functions are then: 1 if ~ 2 Ai ^ ! 2
j ; x bij (~ ; !) = x (5.10) 0 otherwise. The adjoint basis functions are:

(

~ij (~ ; !) = b x

( 1=(jAi j  j
j j) if ~ 2 Ai ^ ! 2
j ; x
0 otherwise.
n m XX

(5.11)

The system of linear equations determining the unknown Lij values is

Lij = Le + ij
where

k=1 l=1

Lkl  Rijkl ;

(5.12)

Rijkl = hT bkl (~ ; !); ~ij (~ ; !)i = jAi j 1 j
j j  x b x 
5.2 Point collocation method

ZZZ


j Ai


0 bkl (h(~ ; !0 ); !0 )  fr (!0 ; ~ ; !) cos ~ d!0 d~ d!: x x x x

(5.13)

The point collocation method ?nds the unknown coef?cients by requiring the ?nite-element approximation to be exact at the prede?ned knot points only. Formally, it uses Dirac-delta type adjoint basis functions which emphasize these knot-points: ~ij (~ ; !) = ij (~ ~ i ; ! !j ): b x x x (5.14) The coef?cients of the linear equation can be expressed as follows:

Z 0 Rijkl = hT bkl (~ ; !); ~ij (~ ; !)i = bkl (h(~ i ; !0); !0)  fr (!0; ~ i ; !j ) cos ~ i d!0 : x b x x x x



(5.15)

5.3 Finite element methods for the diffuse global illumination problem
If the surfaces have only diffuse re?ection and emission — which is a general assumption of the radiosity method [CG85] — then the rendering (or the potential) equation has a simpli?ed form:

L(~ ) = Le (~ ) + x x

Z


H

0 L(h(~ ; !0))  fr (~ )  cos ~ d!0 : x x x

(5.16)

In this case, the BRDF and the radiance depend on the surface point, but are independent of the direction, which reduces the inherent dimensionality of the problem and simpli?es the ?nite-element representation:

L(~ ; !)  x

n X j =1

Lj  bj (~ ): x

(5.17)

A widely used approach is the application of piece-wise constant basis functions for which bj (~ ) = 1 if ~ is x x on surface element Aj and 0 otherwise. An appropriate adjoint basis is ~j (~ ) = 1=Aj if ~ is on surface element b x x Aj and 0 otherwise. Using this set of basis functions, the original rendering equation is projected to the following linear equation: = e+  (5.18) where

1 Rij = hT bj ; ~i i = Ai  b

L L R L ZZ

Ai


0 bj (h(~ ; !0))  fr (~ )  cos ~ d!0 d~ : x x x x

(5.19)

5.3. FINITE ELEMENT METHODS FOR THE DIFFUSE GLOBAL ILLUMINATION PROBLEM

38

Let us replace the directional integral by a surface integral using equation (2.2):

y y d!0 = d~x cosj~ : j~ ~ 2 y
This is true only when ~ and ~ are visible from each other, otherwise the solid angle of the visible surface is x y obviously zero. In order to extent the formula to this case, a visibility indicator v (~ ; ~ ) is introduced, which is 1 xy if the two points are visible from each other and zero otherwise. Using this substitution we obtain

1 Rij = Ai 

ZZ

Ai S

x y bj (~)  fr (~ )  cosj~~  cos ~  v(~ ; ~) d~ d~ : y x xy y x x ~j2 y

0

(5.20)

Taking advantage that the base functions are zero except for their speci?c domain, cos ~ y cos i are constant on these patches, and assuming that the BRDF on patch i is fi , we get

0 = cos j and cos ~ = x
(5.21)

f Rij = Aii 

ZZ

Note that ij is a product of two factors, the albedo of patch i — that is ai ij which describes the geometry of the scene:

F

R

Ai Aj

cos i  cos j  v(~ ; ~) d~ d~ : xy y x j~ ~j2 x y

= fi   —, and a so called form factor
(5.22)

1 Fij = Ai  1 Z Rij = A
i

ZZ

Ai Aj

cos i  cos j  v(~ ; ~) d~ d~ :   j~ ~j2 x y y x x y

So far we have discussed the light propagation problem from the point of view of gathering. For shooting, similar formulae can be produced if incoming direction ! 0 is replaced by the outgoing direction ! = ! 0 in equation (5.19):

Z

Ai
An adjoint equation can also be derived as a special case of equation (5.8). Let Wie be 1 in points of Ai and at the directions of the hemisphere of Ai . The power equation is then

bj (h(~ ; !))  fr (~ )  cos ~ d! d~ : x x x x

(5.23)

where the different forms of

Hij are taken from equation (5.9): f Z Z Hij = Rij  Ai =Aj = Rji  fi=fj = Aij  bi(h(~; !))  cos j d! d~: y y
Aj


P = Pe + H  P;

(5.24)

(5.25)

In order to solve the projected integral equation or the power equation, basically the same techniques can be applied as for the original integral equation: inversion, expansion and iteration. Both the number of unknown variables and the number of equations are equal to the number of surfaces (n). The calculated i radiances represent the light density of the surface on a given wavelength. To generate color pictures at least three independent wavelengths should be selected (say red, green and blue), and the color information will come from the results of the three different calculations. Thus, to sum up, the basic steps are these:

L

1.

2. Describe the light emission ( e ) on the representative wavelengths, or in the simpli?ed case on the wavei length of red, green and blue colors. 3. Solve the linear equation for representative wavelengths. 4. Generate the picture taking into account the camera parameters by any known hidden surface algorithm. In practical circumstances the number of elemental surface patches can be very large, making the form factor computation and the solution of the linear equation rather time consuming.

Fij form factor calculation.

L

5.3. FINITE ELEMENT METHODS FOR THE DIFFUSE GLOBAL ILLUMINATION PROBLEM

39

5.3.1 Geometric methods for form factor computation
Geometric form factor calculation methods approximate the outer integral of the double integral of the form factor from a single value, thus they apply the following approximation:

1 Fij = Ai 

ZZ

Ai Aj

cos i  cos j  v(~ ; ~) d~ d~  Z cos i  cos j  v(~ ; ~) d~:   j~ ~j2 x y y x x y   j~ i ~j2 xi y y x y
Aj

(5.26)

where ~ i is the center of patch i. x

Aj y dy θj dy cos θ j xi - y
2

Aj

N θi 1 Ai j dy cos θ j cos θ i xi - y
2

N

Ai

xi
Fij

=

nj P2

n pixels j P x P resolution

Figure 5.3: Geometric interpretation of hemisphere form factor algorithm and the discrete algorithm for form factor computation

Nusselt [SH81] has realized that this formula can be interpreted as projecting the visible parts of Aj onto the unit hemisphere centered above ~ i , then projecting the result orthographically onto the base circle of this x hemisphere in the plane of ~ i , and ?nally calculating the ratio of the doubly projected area and the area of the unit x circle ( ). Due to the central role of the unit hemisphere, this method is called the hemisphere algorithm. This means that a single row of the form factor matrix can be determined by solving a visibility problem, where the “window” is a half sphere and the “eye” is in the center of surface i. Continuous hemisphere algorithm First an exact algorithm is presented for the calculation of the hemispherical projection [SKe95].

N

Rl

R l +1

+

+ -

~ ~ Figure 5.4: Hemispherical projection of a planar polygon (Rl and Rl1 are two consecutive vertices of the polygon)

5.3. FINITE ELEMENT METHODS FOR THE DIFFUSE GLOBAL ILLUMINATION PROBLEM

40

~ To simplify the problem, consider only one edge line of the polygon ?rst, and two consecutive vertices, Rl ~ l1 , on it (?gure 5.4). Operator  stands for modulo L addition which can handle the problem that the and R next of vertex l is usually vertex l + 1, except for vertex L 1 which is followed by vertex 0. The hemispherical projection of this line is a half great circle. Since the radius of this great circle is 1, the area of the sector formed ~ ~ ~ ~ by the projections of Rl and Rl1 and the center of the hemisphere is simply half the angle of Rl and Rl1 . Projecting this sector orthographically onto the equatorial plane, an ellipse sector is generated, having the area of ~ the great circle sector multiplied by the cosine of the angle of the surface normal N and the normal of the segment ~ l  Rl1 ). ~ (R The area of the doubly projected polygon can be obtained by adding and subtracting the areas of the ellipse sectors of the different edges, as is demonstrated in ?gure 5.4, depending on whether the projections of vectors ~ ~ Rl and Rl1 follow each other clockwise. This sign value can also be represented by a signed angle of the two vectors, expressing the area of the double projected polygon as a summation:
L 11 X

~ ~ (Rl  Rl1 ) ~ ~ ~ 2  angle(Rl ; Rl1 ) jRl  Rl1 j  N: l=0

~

~

(5.27)

Having divided this by  to calculate the ratio of the area of the double projected polygon and the area of the equatorial circle, the form factor can be generated. This method has supposed that surface Aj is above the plane of Ai and is totally visible. Surfaces below the equatorial plane do not pose any problems, since we can get rid of them by the application of a clipping algorithm. When partial occlusion occurs, then either a continuous (object precision) visibility algorithm is used to select the visible parts of the surfaces, or the visibility term is estimated by ?ring several rays to surface element j and averaging their 0/1 associated visibilities [Tam92]. Discrete hemisphere algorithm and its variations Discrete methods subdivide the hemisphere into ?nite number of areas called “pixels” and assume that what can be seen through these small areas is homogeneous (?gure 5.3). Thus it is enough to determine the visibility through a single point in each pixel. The complicated form of the “hemispherical window” can be simpli?ed if the hemisphere is replaced by other immediate surfaces, such as a hemicube [CG85] or a cubic tetrahedron [BKP91]. In these cases the modi?cation of the geometry must be compensated by appropriate weighting functions in area computation. For hemicube and hemishpere algorithms, the window surface consists of normal rectangles, which can be exploited by the built in transformation and scan-conversion hardware of graphics workstations.

Figure 5.5: A complex scene rendered by the hemicube algorithm (University of Cornell)

Chapter 6

Numerical quadrature for high dimensional integrals
The solution of the rendering equation requires the numerical evaluation of high-dimensional integrals. Numerical quadrature formulae take ?nite samples from the domain, evaluate the integrand for these samples and generate the integral as a weighted sum of the integrand values. The general form for the solution is

I = f (z) dz 

Z

N X i=1

where i is a sample point from the s-dimensional domain V , and w is the weighting function of the given quadrature rule. Classical quadrature rules usually trace back the computation of multi-dimensional integrals to a series of one-dimensional integrals. Thus they select the coordinates of the sample points independently, and the highdimensional points are de?ned as the Cartesian product of the 1-dimensional sample sets. The simplest techniques, including the brick-rule, the trapezoidal-rule, Simpson-rule, etc. space the abscissas equally, which results in a regular grid in the domain. More sophisticated techniques, such as the Gaussian quadrature, selects the sample points non uniforly along each abscissa to decrease the error of the integral.

z

V

f (zi )  w(zi )

(6.1)

?f N ?f 1/N 0 1 N points

n points

?f

n points

Figure 6.1: Error of the brick rule in one and two dimensional spaces

Let us evaluate the error of the brick rule in one and two dimensional spaces (?gures 6.1). Suppose that N sample points are used in the domains of [0; 1] and [0; 1]2 , respectively. In the one dimensional space the error is the sum of the areas of triangles that are between the function and the series of bricks. The average height of a triangle is f=(2N ) where f is the variation of the integrand and N is the number of bricks. Since the base of a triangle is 1=N and the number of triangles is N , the error is:

1 1 1 = f  N  N = f = O N : (6.2) 2N 2N In thep dimensional space the N sample points should be distributed along two coordinates, thus we have two just n = N points in a single row and column. Now the error is the sum of the volume of the roof-like objects that are between the function and the series of bricks. The average height of a volume element is f=(2n), its
41

 

6.1. MONTE-CARLO QUADRATURE

42

base has 1=n2 area and there are n2 elements. Thus the error is

f 1 p 2 = n  n2  n2 = f = O p1 : 2 2 N N





(6.3)

This idea can simply generalized to any dimensions and we can conclude that the integration error of the brick rule in an s-dimensional space is proportional to f=N 1=s (if is the total change, i.e. the variation, of the integrand). From another point of view, it means that the required number of sample points to ?nd an estimate with error  is proportional to (f=)s , thus the computational complexity is exponential with regard to the dimension of the domain. This phenomenon is called as the dimensional explosion or dimensional core. Concerning other classical quadrature rules, f measures higher order changes, such as the distance from the piece-wise linear or polynomial functions, but they also exhibit dimensional core. The dimensional explosion can be avoided by Monte-Carlo [Sob91] or quasi-Monte Carlo [Nie92] integration methods. Monte-Carlo methods trace back the estimation of an integral to the calculation of an expected value, which is estimated by averaging random samples. Quasi-Monte Carlo techniques, on the other hand, use deterministic samples that are uniformly distributed in the integration domain.

6.1 Monte-Carlo quadrature
Monte-Carlo integration converts the calculation of an integral to an equivalent expected value problem. Assume that a random vector variable is uniformly distributed in V , thus its probability density is p( ) = 1=V . The expected value of random variable f ( ) is

z

1 1 E [f (z)] = f (z)  p(z) dz = f (z)  V dz = V  I;
V V

z

z

Z

Z

(6.4)

thus the required integral can easily be found if this expected value is available. According to the theorems of large numbers, if independent random samples 1 ; 2 ; : : : ; N are generated using probability density p, then the ^ expected value can be estimated by the average f :

z z

z

1 E [f (z)]  f^ = N

N X i=1

f (zi ):

(6.5)

^ Estimator f is also a random variable whose expected value is equal to E [f ( )] = I=V . Suppose that the variance ^ of f ( ) is  2 . If the samples are independent random variables, then the variance of estimator f is:

z

z

1 D2 f^ = N 2 2 =  : N i=1
Thus the standard deviation of the estimator is

hi

N X

2

(6.6)

hi D f^ = p : (6.7) N ^ According to the central limit theorem, estimator f will have normal distribution asymptotically, with mean p E [f (z)] = I=V and standard deviation = N . Examining the shape of the probability density of the normal Z N VX p j f (z) dz N f (zi )j < 3V  :
V i=1

distribution, we can conclude that the probability that the distance between the variable and the mean is less than 3 times the standard deviation is 0.997. Thus with 0.997 con?dence level we can say that the (probabilistic) error bound of the Monte-Carlo quadrature is

N

(6.8)

Let us realize that this bound is independent of the dimension of the domain! This property means that by MonteCarlo quadrature the dimensional explosion can be avoided.

6.2 Quasi-Monte Carlo quadrature
In the previous section we concluded that randomly distributing the sample points solves the dimensional core problem. We show that it is also possible with deterministically selected sample points. The resulting method is

6.2. QUASI-MONTE CARLO QUADRATURE

43

called quasi-Monte Carlo quadrature. For the sake of simplicity, assume that an s-dimensional funcion f ( ) needs to be integrated over the domain [0; 1]s . In order to guarantee that f is Riemann-integrable, f is assumed to be piece-wise continuous. This integral is approximated by a ?nite sum as follows:

z

Z

The question is how the sequence of sample points 1 ; 2 ; : : : ; N should be selected to minimize the error of the integral quadrature. A minimal requirement is the stability of the sequence, which means that in asymptotic sense the error is zero for any Riemann integrable function:

z2[0;1]s

1 f (z) dz  N

N X i=1

f (zi ):

(6.9)

z z

z

Z

z2[0;1]s

1 f (z)dz = Nlim N !1

N X i=1

f (zi ):

(6.10)

Sequences meeting this requirement are called uniform or equidistribution (or more precisely 1-equidistribution) sequences [Nie92]. In order to ?nd other necessary requirements for uniform sequences, let us consider the integration of a very simple function which is 1 inside a d-dimensional “brick” originating at the center of the coordinate system and 0 elsewhere: 1 if 0  j1  v1 ; 0  j2  v2 ; : : : ; 0  js  vs ; L( ) = (6.11) 0 otherwise. s Let us denote the volume of this brick by V (A) = j =1 vj : Integrating this function we have:

z

(

z

z

z

Z

Q

z2[0;1]s

L dz =

s Y

j =1

vj = V (A):

(6.12)

If the number of sample points that are inside the s-dimensional “brick” A is m(A), then the ?nite approximation sum is

Since the exact value of the integral is V (A) now, for uniform sequences, the average number of sample points that are located inside a volume should be proportional to the size of this volume:

N 1 X f (z ) = m(A) : i N N i=1

(6.13)

N !1

(A lim mN ) = V (A):

(6.14)

If the size of the sequence of the sample points is not in?nite, then the proportionality requirement can only be approximately met. The maximal error in this approximation is called the discrepancy (or the star-discrepancy) of the sequence:

( D (z1 ; z2 ; : : : zN ) = sup j mNA) V (A)j:
A

(6.15)

If a sequence is appropriate for integral-quadrature, then the approximation sum is expected to converge to the real value of the integral, which requires the discrepancy to converge to 0. This is another necessary requirement which is derived from the examination of a very simple function. Note that this requirement also emphasizes the uniformity of those sequences that can be used for numerical integration. However, this requirement is not only necessary, but also suf?cient, since any Riemann-integrable function can be approximated by piece-wise constant step-functions with arbitrary precision. To show how step-functions can do this, an example of a 2-dimensional function is shown in ?gure 6.2.

6.2.1 Error Analysis for integrands of ?nite variation: Koksma-Hlawka Inequality
In this section, an upper-bound is given to the error of the quasi-Monte Carlo quadrature for functions that have bounded and piece-wise continuous mixed derivatives. The error of the quadrature is shown below:

j

1 f (z) dz N f (zi )j: i=1 z2[0;1]s

Z

N X

(6.16)

6.2. QUASI-MONTE CARLO QUADRATURE

44

f1
b1=f1-f2-f3+f4

f3 f2 f4

b2=f2-f4

b3=f3-f4

=

+

+

+

b4=f4

Figure 6.2: De?nition of functions from bricks originating at the center of the coordinate system

Intuitively this error must depend on two independent factors. On the one hand, if the discrepancy of the sequence of sample points is high, then there are large regions where there are no sample point at all, which increases the error. This means that the error is expected to be proportional to the discrepancy of the sequence of sample locations. On the other hand, the error also depends on how quickly the function changes between the sample points. If the function can change signi?cantly in small domain, then the error can be quite large. However, if the slope of the function is small, then nothing dramatic happens between the sample points thus the error will be small. Measures describing how rapidly a function can change are called variations. For a 1-dimensional function the variation in the sense of Vitali is de?ned as:

VV (f (x)) = lim sup
N M XX i=1 j =1

N X i=1

jf (xi+1 ) f (xi )j:

(6.17)

For a 2-dimensional function, the de?nition is analogous:

VV (f (x; y)) = lim sup

jf (xi+1 ; yj+1 ) f (xi+1 ; yi ) f (xi ; yi+1 ) + f (xi ; yi )j:

(6.18)

Note that for higher dimensions, the variation of Vitali does not provide a good measure: if the function is constant in x, for instance, then the variation is zero, regardless how the function changes depending on y . Thus, it is worth introducing a somehow more stronger variation type, called the Hardy-Krause variation. The variation in the sense of Hardy-Krause is the sum of the variations of the function and of its restrictions to the end of the domain. For dimension 2, the new variation is:

VHK (f (x; y)) = VV f (x; y) + VV f (x; 1) + VV f (1; y):

(6.19)

If a function has bounded and piece-wise continuous mixed derivatives, then its variation is ?nite. For a 2dimensional function meeting this requirement the variation can be given by the following formula:

Z1 Z1 @ 2f (u; v) Z1 @f (u; 1) Z1 @f (1; v) VHK (f (u; v)) = @u@v du dv + @u du + @v dv:
0 0 0 0

(6.20)

The property that a function is not continuous does not necessarily mean that the variation is in?nite. If at most ?nite or countable in?nite discontinuities occur at hyper-planes parallel to the coordinate axes, then the variation is still ?nite. An example of a discontinuous function that have ?nite variation is

f (x; y) =

( 1 if x > x0;

0 otherwise.

(6.21)

However, when the discontinuity is not parallel to the coordinate axes, then the variation is in?nite. A simple function of in?nite variation is [De? 89]: a

f (x; y) =

( 1 if x > y;

0 otherwise.

(6.22)

6.2. QUASI-MONTE CARLO QUADRATURE

45

Now, let us turn to the error of quasi-Monte Carlo integration. The following formula, which expresses the previous intuition that the error depends on the uniformness of the sample points and on the variation of the integrand, is called the Koksma-Hlawka inequality:

j

Z

z2[0;1]s

f (z) dz

N 1 X f (z )j  V  D (z ; z ; : : : z ): i HK 1 2 N N i=1

(6.23)

For the sake of notational simplicity, the Koksma-Hlawka inequality is proven here only for the two-dimensional case (s = 2). The outline of the proof is as follows. We start with the assumption that the function has piece-wise continuous and bounded mixed derivatives, thus the mixed derivatives are Riemann-integrable and their integration gives back the original function. First the function is expressed as the integral of its derivatives, and the function values in the approximation sum are converted accordingly. On the other hand, the integral of f is converted to a similar form using partial integration. Finally the difference of the approximating sum and the integral is examined and upperbounded. First, the value of function f (u; v ) is expressed by its derivatives at point x; y :

f (x; y) = f (1; 1) + [f (1; 1) f (x; 1) f (1; y) + f (x; y)] [f (1; 1) f (x; 1)] [f (1; 1) f (1; y)] = f (1; 1) +
where

Z1 Z1
x y

fuv (u; v) du dv

Z1
x

fu (u; 1) du

Z1
y

fv (1; v) dv:

2 f (u; fuv (u; v) = @ @u@vv) ; fu (u; v) = @f (u; v) ; fv (u; v) = @f (u; v) : @u @v Let us introduce the step function : ( 1 if u  0; v  0; (u; v) = 0 otherwise. Using this step function, the domains of the integrals can be extended to [0; 1], as follows:

f (x; y) = f (1; 1) +
Substituting

Z1 Z1
0 0

fuv (u; v)  (u x; v y) du dv

Z1
0

fu (u; 1)  (u x; 1) du

Z1
0

fv (1; v)  (1; v y) dv:

zi = (xi; yi) into this formula we have:

f (zi ) = f (1; 1) +

Z1 Z1
0 0

fuv (u; v)  (u xi ; v yi ) du dv

Z1
0

fu (u; 1)  (xi ; 1) du

Z1
0

fv (1; v)  (1; yi ) dv:

Averaging these formulae for i = 1; 2; : : : N , the approximating sum has the following form:

N 1 X f (z ) = f (1; 1) + 1 N i=1 i N
where

Z1 Z1
0 0

1 fuv (u; v)  m(u; v) du dv N

Z1
0

1 fu (u; 1)  m(u; 1) du N

Z1
0

fv (1; v)  m(1; v) dv:

m(u; v) =

N X i=1

(u xi ; v yi );

which is the number of points located in the rectangle [(0; 0); (u; v )]. The integral

Z

z2[0;1]2

f (z) dz =

Z1 Z1

v=0 u=0

f (u; v) du dv

can also be converted to a similar form, if partial integration is applied. First the inner integral is considered:

Z1

u=0

f (u; v) du =

Z1

u=0

f (u; v)  1 du = f (1; v)

Z1

u=0

fu (u; v)  u du = F (v)

Then the outer integral is processed in a similar way:

Z1 Z1

v=0 u=0

f (u; v) du dv =

Z1

v=0

F (v) dv =

Z1
v=0

F (v)  1 dv = F (1)

Z1
v=0

Fv (v)  v dv

6.2. QUASI-MONTE CARLO QUADRATURE

46

Substituting the de?nition of F (v ) we get:

Z1 Z1

v=0 u=0

f (u; v) du dv = f (1; 1) +

Z1 Z1
0 0

fuv (u; v)  uv du dv

Z1
0

1 fu (u; 1)  u du N

Z1
0

fv (1; v)  v dv:

Thus the error of quadrature is then:

j j

Z1 Z1

0 0 0 0 0Z1 Z1 1  Z1 Z1 u; @ jfuv (u; v)j  du dv + jfu (u; 1)j du + jfv (1; v)j dvA  sup m(N v) 0 0 0 0
u;v
This is exactly what we wanted to prove.

u; fuv (u; v)  m(N v) uv du dv





N 1X f (u; v) du dv N f (zi )j = i=1 v=0 u=0

Z1 Z1

Z1

u; fu (u; 1)  m(N 1) u du





Z1

fv (1; v)  m(1; v) v dvj  N uv





 = VHK  D(z1 ; z2 ; : : : zN ):

According to this inequality, the error can be upperbounded by the product of two independent factors, the variation of the integrand and the discrepancy of the used sample set. The discrepancy shows how uniformly the set is distributed [Shi91a]. This immediately presents two orthogonal strategies to improve the quality of quadratures. Either we try to make the function ?at by appropriate variable transformations, or use very uniformly distributed sample sets. The ?rst technique is called importance sampling [Sob91], while the second involves the strati?cation [Sob91, Mit96, Arv95] of random points or the application of low-discrepancy series [Nie92, War95, PFTV92, Knu81, Sob91]. If the integrand has ?nite variation, then the error is proportional to the discrepancy of the sequence of sample locations. For carefully selected sample points the discrepancy can converge to zero with almost linear speed. Thus, p quadratures having almost linear convergence can be obtained in this way, which is better than the O(1= N ) speed of Monte-Carlo quadratures. Note, on the other hand, that functions having in?nite variation can also be integrated by quasi-Monte Carlo quadrature. The quadrature will be asymptotically exact for any uniform sequence and for any Riemann integrable function. The fact that the Koksma-Hlawka inequality cannot be applied means that it cannot be used to provide a qualitative measure for the speed of the convergence. Practical experience shows that quasi-Monte Carlo integration outperforms the classical Monte-Carlo integration even for discontinuous functions [SKDP99]. However, the difference in the effectiveness becomes signi?cantly less when the integrand has in?nite variation. This phenomenon will be investigated in the subsequent sections. Since the integrand of the rendering and potential equations is usually discontinuous, this case is very important for computer graphics applications.

6.2.2 Generation of the sample points
As a conclusion of error analysis we can state that we need very uniform sequences for quasi-Monte Carlo quadrature. Regular grids should be rejected because of their O(1=N 1=s ) discrepancy which results in dimensional explosion.
Regular grid 1 1 Random points 1 First 100 Halton points of base (2, 3)

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0 0 0.2 0.4 0.6 0.8 1

0 0 0.2 0.4 0.6 0.8 1

0 0 0.2 0.4 0.6 0.8 1

Figure 6.3: 100 points distributed by a regular grid (left), random distribution (middle) and Halton low-discrepancy sequence (right)

6.2. QUASI-MONTE CARLO QUADRATURE

47

Monte-Carlo method proposed the application of random samples to avoid dimensional explosion. This can also be justi?ed by analyzing the discrepancy. In fact, the discrepancy of a uniformly distributed random series is

O
assymptotically with probability 1.

r

log log N 2N

!

In order to prove this, let us de?ne a new random variable i from uniformly distributed random sample way: 1 if i is in an s-dimensional brick A that originates at the center,

i =

( z

zi in the following

Note that if

zi is uniformly distributed in [0; 1], then the expected value and the variance of i are
E [i ] = V (A); D2 [i ] = V (A) V 2 (A)  1 ; 4

0 otherwise.

where V (A) is the volume of brick A. If the samples are independent random variables, the generated 1 ; 2 ; : : : ; N random variables will also be independent, and of the same distribution having mean E [ ] and variance D2 [ ]. According to the theorem of iterated logarithm [R? n81], the difference between the average of independent random e variables 1 ; 2 ; : : : ; N of the same distribution and their mean E [ ] can be upperbounded in the following way:

Pr

(

X i lim sup N

r 2 2 log log N ) E []  D []  N = 1;

In our case the standard deviation D2 [ ] cannot exceed 1=4 and

X i sup N

m(A) E [] = sup N A

V (A) = D (z1 ; z2 ; : : : zN ); r
log log N = 1: 2N

thus the theorem of iterated logarithm becomes what we want to prove:

Pr lim D (z1 ; z2 ; : : : zN ) 

(

)

6.2.3 Generation of low-discrepancy sequences
Beyond random sequences, however, there are deterministic sequences that have even better discrepancy. The discrepancy of the best sequences known is in the order of O(logs N=N ) or even in the order of O(logs 1 N=N ) if N is known before starting the sequence. These sequences are called low-discrepancy sequences. There are many sequences published in the literature [Nie92, War95, De? 89, Knu81]. The most famous one is probably the a Halton-sequence (its one-dimensional version is also called Van der Corput sequence). The element i of the one-dimensional Halton sequence of base b is de?ned as the radical inverse of the expansion of i in base b. This means that number i is expanded in radix b, then the number is mirrored onto the “radix” point. The ?rst few points in base 2 are shown in table 6.1.

i
1 2 3 4 5 6 7

binary form of i 1 10 11 100 101 110 111

radical inverse 0.1 0.01 0.11 0.001 0.101 0.011 0.111

Hi
0.5 0.25 0.75 0.125 0.625 0.375 0.875

Table 6.1: The calculation of the ith Halton point Hi in base 2

Why is this sequence uniform? Note that the construction algorithm generates as binary form of i all binary combinations of length k before producing a combination of length k + 1. This means that after the radical inverse the sequence Hi will visit all intervals of length 2 k before putting a new point in an interval already visited. Thus the sequence is really uniform.

6.2. QUASI-MONTE CARLO QUADRATURE

48

First 10 Halton points of base (2, 3) 1 1

First 100 Halton points of base (2, 3) 1

First 1000 Halton points of base (2, 3)

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0 0 0.2 0.4 0.6 0.8 1

0 0 0.2 0.4 0.6 0.8 1

0 0 0.2 0.4 0.6 0.8 1

Figure 6.4: The distribution of the ?rst 10,100 and 1000 Halton points in 2 dimensions

On the other hand, as k increases, the algorithm produces a single point in each interval of length 1=2, then in each interval of length 1=4, etc. thus the sequence is not only asymptotically uniform, but also the ?rst N points are fairly uniformly distributed (this is guaranteed by the property that the radical inverse makes the most signi?cant bit the most rapidly changing bit). This is also true in other radix systems as well. If the base is b, then the Halton sequence will place a sample point in all intervals of size b k before putting a new point into any intervals. A Halton sequence is able to generate points that are uniformly distributed in the 1-dimensional [0; 1] interval. If higher dimensional regions, such as rectangles, cubes, etc. should be ?lled uniformly, then different coordinates of the sample vectors can be generated from Halton sequences of different base numbers. In order for these vectors to uniformly ?ll the whole region, the “interdependence” of different coordinates should be as little as possible. To examine this, let us assume that a two-dimensional Halton sequence is generated with base numbers b1 and b2 . According to the behavior of the Halton sequence, the generation algorithm would visit all columns of width b1 k1 before visiting a column again, and similarly it would visit all rows of height b2 k2 before putting a new point into a row. The columns and rows form bk1  bk2 cells. Uniformness in the two-dimensional space means that the 1 2 algorithm is expected to put a point in each cell before putting a second point in any cells. Since the periodicity of the columns and rows are bk1 and bk2 , respectively, the periodicity of the cells is the smallest common multiple of 1 2 bk1 and bk2 . This equals to the their product, that is total number of cells, if b1 and b2 are relative primes. 1 2 This can also be stated in a general way. If a multi-dimensional Halton sequence is to be constructed, then the base numbers of different coordinates must be relative primes. A C++ class that can initialize an arbitrary Halton point and then it can generate incrementally all subsequent points using a very fast algorithm [Kel96b] is presented in the following:
class Halton { double value, inv_base; Number( long i, int base ) { double f = inv_base = 1.0/base; value = 0.0; while ( i > 0 ) { value += f * (double)(i % base); i /= base; f *= inv_base; } } void Next( ) { double r = 1.0 - value - 0.0000000001; if (inv_base < r) value += inv_base; else { double h = inv_base, hh; do { hh = h; h *= inv_base; } while ( h >= r ); value += hh + h - 1.0; } } operator double() { return value; } };

6.3. IMPORTANCE SAMPLING

49

6.3 Importance sampling
Importance sampling is a well-known method of Monte-Carlo integration to reduce variance. The basic idea is to use non-uniform distribution to ?nd sample points, which places more samples where the function is large. More speci?cally it means that

 f (z)  1 X f (zi ) 3V  Z f (z) N I = f (z) dz = p(z)  p(z) dz = E p(z)  N  p(z )  p ; (6.24) N i i=1 V V where p(z) is a probability density in V , the zi points are selected according to this probability density, and the variance  2 is de?ned by  f (z)  " f (z) 2 # Z  f (z) 2 2 = D2 p(z) = E p(z) I = (6.25) p(z) I  p(z) dz: V The probability density p(z) should be selected to minimize the variance. As can be shown easily, the variance can be minimized if p(z) is proportional to the integrand f (z). In order to demonstrate this, let us express the ratio of Z
the integrand and the probability density in the following way:

(z) where = E [ f (z) ] and p

2 = E [( +  (z) E [ +  (z)])2 ] = 2  E [((z))2 ] = 2 : This is obviously minimal if = 0, when the variance is also zero.

V

R ((z))2  p(z) dz = 1. Using this, the variance of the integral quadrature is
(6.27)

f (z) = +  (z); p(z)

(6.26)

Thus in Monte-Carlo integration it is worth applying probability distributions that are large where the integrand is large and small where the integrand is negligible.

6.3.1 Generation of a random variable with a prescribed probability density
We concluded that importance sampling requires random samples generated from a probability density which is proportional — at least approximately — to the integrand. This section examines how such samples can be found. First, let us consider the 1-dimensional case and assume that the domain of the integration is an interval [a; b]. Suppose that we want samples with a probability density that is proportional to a function g (z ). This function is an approximation of the integrand f . If this function is different from the integrand, then the importance sampling is not optimal. The needed probability density p(z ) can be obtained by scaling function g to integrate to 1 as probability densities do:

p(z ) = Rb g(z ) : g(z )dz
a

(6.28)

Note that this gives us a clear explanation why we should use non-optimal densities. If g were equal to f , then the construction of the probability density would require the integral of f . From the probability density p, probability distribution function P is obtained:

P (z ) = g(Z ) dZ:
a

Zz

(6.29)

Prf < z g = PrfP 1(r) < z g = Prfr < P (z )g = P (z ): (6.30) since P (z ) is not decreasing and the probability that r is in an interval of [0; 1] equals to the size of this interval.

A random variable  having probability distribution P can be constructed by transforming another random variable r, which is uniformly distributed in the [0; 1] interval, with the  = P 1 (r) transformation. To prove this, the probability distribution of  is examined:

The multi-dimensional case can be traced back to a series of 1-dimensional constructions. After normalization we have p( ) = g( ) : (6.31)

z

V

z R g(z)dz

6.3. IMPORTANCE SAMPLING

50

The probability density is expressed as a product of 1-dimensional conditional densities:

p(z1 ; z2 ; : : : ; zs ) = p1 (z1 jz2 ; : : : ; zs )  p2 (z2 jz3 ; : : : ; zs )  : : :  ps (zs ): p(z1 ; z2 ; : : : ; zs ) = p1 (z1 )  p2 (z2 )  : : :  ps (zs ):
6.3.2 Importance sampling in quasi-Monte Carlo integration

(6.32)

For these conditional densities the same method can be used recursively as for the 1-dimensional case. If the different coordinates are independent random variables, then we obtain: (6.33)

Classical Monte-Carlo integration places more samples where the integrand is large. The same basic idea also makes sense in quasi-Monte Carlo integration. However, for formal analysis, we have to ?nd another approach since terms like probability density or variance cannot be applied in the deterministic quasi-Monte Carlo framework. The alternative formulation is the integration using variable transformation. Suppose that a function f needs to be integrated in domain V and we have a monotonously increasing mapping T that maps this domain onto V 0 . The integral can also be evaluated in the new domain using the following formula:

Z

1 where @T @ y(y) is the Jacobi determinant of the inverse transformation.

V

f (z) dz = f (T 1(y)) @T @ y(y) dy;
1

Z





(6.34)

V0

If quasi-Monte Carlo integration is used, then domain V 0 is [0; 1]s . In order to quantify the error of the quadrature, we can use the Koksma-Hlawka inequality which expresses the error-bound as a product of the discrepancy of the sample points and the variation of the function. Since the discrepancy of the sample points is independent of the function to be integrated, the error-bound can be controlled by the appropriate transformation of the integrand to reduce variation. In order to reduce the variation, the function should be ?attened. In the ideal case when the integrand is constant, the variation is 0. To make the transformed integrand constant, the Jacobi determinant should be inversely proportional to f . Since the intuitive meaning of the Jacobi determinant is the compression or expansion ratio of the two corresponding sub-cubes in V and V 0 respectively, this criterion states that if the sample points are uniformly distributed in V 0 , then the transformation will concentrate them around regions where f is high. For the sake of simplicity, the details are discussed only for the 1-dimensional case, when the variable transformation has the following form

zmax Z zmin

f (z ) dz = f (T 1(y))  dT dy(y) dy:
0

Z1

1

(6.35)

In the ideal case mapping T makes the integrand have zero variation, that is constant C :

f (T 1(y))  dT dy(y) = C :
From this we can have

1

1 T (z ) = C 

Zz
zmin

f (Z ) dZ: y should be transformed by the inverse of the
= 1. Thus the constant C should be equal

Since mapping T is expected to map to [0; 1], we require that T (zmax ) to following function

zmin

zmax R

f (Z ) dZ .

Summarizing, the uniformly distributed point

Rz f (Z ) dZ T (z ) = zzmin max R f (Z ) dZ :
zmin

(6.36)

Note that this is the same transformation as has been derived for the random samples. It means that the method of importance sampling is independent of whether random samples are transformed in Monte-Carlo quadrature, or deterministic samples are transformed in quasi-Monte Carlo integration.

6.3. IMPORTANCE SAMPLING

51

6.3.3 Metropolis sampling
Importance sampling requires a probability density that is proportional to a function g which, in turn, should mimic the integrand. In the previous section we discussed a random variable transformation method that can generate samples with the required probability density. This method expects function g to be symbolically integrable and its primitive function to be invertible. This requirement often contradicts the requirement of the good approximation of the original integrand. This section discusses a different sampling strategy, proposed by Metropolis et. al [MRR+ 53], that has less expectations towards function g . In fact, it only supposes that the integral of g can be determined either symbolically or numerically. The Metropolis method carries out sampling by establishing a discrete time Markov process f 1 ; 2 ; : : : ; i ; : : :g in the space of samples, whose limiting distribution is proportional to the selected function. A discrete time Markov process visits states which correspond to samples. The Metropolis algorithm constructs this process such that the probability of being in a given state converges to a limiting value and in the limiting case g ( ) = b  p( ), where b = V g( ) d . A Markov process can be de?ned by the state transition probabilities, that is by the conditional probability of the next state provided that the current state is known. In Metropolis method, the next state i+1 of this process is found by letting an almost arbitrary tentative transition function T ( i ! t ) generate a tentative sample t which is either accepted as the real next state or rejected making the next state equal to the actual state using an acceptance probability a( i ! t ). Thus the state transition probability density from state to a different state is: P ( ! ) d = T ( ! ) d  a( ! ): (6.37)

z z

z

R

z z

z

z

z

z

z

z

y

z

z

x

x y y

x y y x y

The event that the process remains in the same state is the complement of moving to any other state, thus the probability of no state transition happening is:

1

Z

x2V;x6=y

P (x ! y) dy = 1

Z

x2V;x6=y

T (x ! y)  a(x ! y) dy:

(6.38)

Let us denote the probability density of being in state at step n by pn ( ). Using the total probability theorem and taking advantage that in Markov processes the future depends only on the present state and is independent of the past, the following recursion can be established for these state probabilities:

x

x

pn+1 (y) =

Z
x2V;x6=y

0 1 Z pn (x)  P (x ! y) dx + B1 P (y ! x) dxC  pn (y): @ A
x6=y

(6.39)

If the limiting probability p(

y) = nlim pn(y) exists, then it should be the ?xed point of this recursion, that is: !1 Z p(y) = p(y) + p(x)  P (x ! y) p(y)  P (y ! x) dx (6.40)
x2V;x6=y

If the Markov process is ergodic, then the limiting distribution is unambigous and is independent of the initial state of the process. The process is ergodic if after a given number of steps any state can be reached from any other state with non-zero probability. The core idea of Metropolis is to construct acceptance probability a( ! ) in such a way that the limiting probability of the process will be p( ) = g ( )=b. Substituting this goal into equation (6.40) we get:

z

g(y) = g(y) +

z Z

x

y

This holds when the total incoming and outgoing ?ows of state require that:

x are balanced. One way of ensuring this is to g(x)  P (x ! y) = g(y)  P (y ! x): (6.42)
(6.43)

V;x6=y

g(x)  P (x ! y) g(y)  P (y ! x) dx:

(6.41)

This condition — that is also called as the detailed balance — means that the transitions between any two states are balanced. Using the formulae (6.37) for state transition probabilities, we can further obtain:

g(x)  T (x ! y)  a(x ! y) = g(y)  T (y ! x)  a(y ! x):

6.3. IMPORTANCE SAMPLING

52

Thus the required ratio of the two acceptance probabilities is:

Any acceptance probability satisfying this requirement makes the limiting distribution proportional to g . Considering the speed of convergence, on the other hand, the state transition probabilities and consequently the acceptance probabilities should be large. Since a probability cannot exceed 1, the optimal acceptance probability is: The algorithm generating a trajectory of the Markov process to obtain samples f 1 ; 2 ; : : : ; N g is as follows:

a(y ! x) = g(x)  T (x ! y) : a(x ! y) g(y)  T (y ! x)

(6.44)

z z for i = 1 to N do Based on the actual state zi , choose another random, tentative state zt using T (zi  zt ) a(zi  zt ) = (g(zt)  T (zt  zi ))=(g(zi )  T (zi  zt )) if a(zi  zt )  1 then zi+1 = zt
Generate uniformly distributed random number r in [0; 1]. if r < a( i t ) then i+1 = t else i+1 = i endif endfor else

) T (y ! x) a(x ! y) = min g(y)  T (x ! y) ; 1 : g(x





z

// accept with probability a(

zi  zt )

z z

z z

z z

6.3.4 Application of the VEGAS algorithm
If neither the integrand nor its approximation are available explicitly, no probability density can be constructed before starting the generation of the samples. However, as the sampling process goes on, the history of previous samples can provide information which the important regions are, thus a probability density can be built adaptively using the information gathered previously. This means that as the samples are evaluated not only the quadrature is computed but a probability density is also approximated which can be proportional to the current estimate of the integrand. This is quite straightforward, but also generates a problem. The estimated probability density is also a highdimensional function, thus its storage would need a lot of storage space. The VEGAS method [Lep80] is an adaptive Monte-Carlo technique that generates a probability density for importance sampling automatically and in a separable form. The reason of the requirement of separability is that probability densities are computed from and stored in histogram tables and s number of 1-dimensional tables need much less storage space than a single s-dimensional table. Formally, let us assume that the probability density can be de?ned in the following product form:

p(z1 ; z2; : : : ; zs ) = p1 (z1 )  p2 (z2 )  : : :  ps (!s ):
It can be shown [Lep80] that the optimal selection of p1 is

(6.45)

p1 (z1 ) /

sZ

Z f 2(z1; : : : ; zD ) : : : p (z ) : : : p (z ) dz2 : : : dzs ; 2 2 s s

(6.46)

and similar formulae apply to p2 ; : : : ; ps . These p1 ; : : : ps functions can be tabulated as 1-dimensional tables. The j +1) j th element of the ith table represents the importance of the directional region where zi 2 [ N ; (jN ]: This immediately presents a recursive importance sampling strategy. The algorithm is decomposed into phases consisting of a number of samples. At the end of each phase weights p1 ; : : : ; ps are re?ned, to provide a better probability density for the subsequent phase. Assuming that p1 ; : : : ; ps are initially constant, a standard MonteCarlo method is initiated, but in addition to accumulating to compute the integral, p1 ; : : : ; ps are also estimated using equation (6.46). Then for the following phase, the samples are selected according to the new pi functions. In order to calculate a sample for zi , for instance, a single random value is generated in the range of 0 and the sum of all elements in the table de?ning pi . Then the elements of the table is retained one by one and summed to a running variable. When this running variable exceeds the random sample, then the searched region is found. The value in this region is then found by uniformly selecting a single point from the region.

Chapter 7

Random walk solution of the global illumination problem
Recall that expansion obtains the measured power in the form of an in?nite Neumann series: ML = 1 MT i Le . i=0 The terms of this in?nite Neumann series have intuitive meaning as well: MT 0 Le = MLe comes from the emission, MT 1 Le comes from a single re?ection, MT 2 Le from two re?ections, etc.

P

7.1 Why should we use Monte-Carlo expansion methods?
Expansion techniques require the evaluation of very high-dimensional — in fact, in?nite dimensional — integrals. When using classical quadrature rules for multi-dimensional integrals [PFTV92], such as for example the trapezoidal rule, in order to provide a result with a given accuracy, the number of sample points is in the order of O(N s ), where s is the dimension of the domain. This phenomenon is called the dimensional core or dimensional explosion and makes classical quadrature rules prohibitively expensive for higher dimensions. However, Monte-Carlo or quasi-Monte Carlo techniques distribute the sample points simultaneously in all dimensions, thus they can avoid dimensional explosion. For example, the probabilistic error bound of MonteCarlo integration is O(N 0:5 ), independently of the dimension of the domain. s-dimensional low discrepancy series [Nie92] can even achieve O(logs N=N ) = O(N (1  )) convergence rates for ?nite variation integrands. Furthermore, classical quadrature cannot be used for in?nite dimensional integrals, thus the Neumann series should be truncated after D terms. This truncation introduces a bias of order D+1  jjLe jj=(1 ), where  is the contraction of the light transport operator. Using a Russian roulette based technique, on the other hand, Monte-Carlo methods are appropriate for even in?nite dimensional integrals. Thus we can conclude that the stochastic approach is indispensable for expansion methods. In computer graphics the ?rst Monte-Carlo random walk algorithm — called distributed ray-tracing — was proposed by Cook et al. [CPC84], which spawned to a set of variations, including path-tracing [Kaj86], lighttracing [DLW93], bi-directional path-tracing [LW93, VG95], Monte-Carlo radiosity [Shi91b, Neu95, PM95], and two-pass methods which combine radiosity and ray-tracing [Shi90, ZS95, WCG87]. The problem of naive generation of walks is that the probability that a shooting path ?nds the eye is zero for a pin-hole camera or very small if a non-zero aperture camera model is used, while the probability that a gathering random path ends in a lightsource may be very little if the lightsources are small, thus the majority of the paths do not contribute to the image at all, and their computation is simply waste of time. Note that shooting is always superior for view-independent algorithms since they do not have to face the problem of small aperture. Thus, on the one hand, random walk must be combined with a deterministic step that forces the walk to go to the eye and to ?nd a lightsource. On the other hand, importance sampling [Sob91] should be incorporated to prefer useful paths along which signi?cant radiance is transferred. Steps of the walk transfer the radiance or the potential in the scene. The source and destination of the transfer can be points in the case of continuous methods or patches in the case of ?nite-element methods. If the algorithm is such that it always selects a single source for shooting or single destination for gathering, then the method is called local method. On the other hand, if many sources and destinations are taken into consideration simultaneously in each transfer, then the method is called global method or multi-path method [Sbe96].

53

7.2. QUASI-MONTE CARLO QUADRATURE FOR THE RENDERING EQUATION

54

7.2 Quasi-Monte Carlo quadrature for the rendering equation
Quasi-Monte Carlo walk techniques mean that instead of generating the next direction randomly, the direction is sampled from a low-discrepancy point set. Since the low-discrepancy sequences have better asymptotic disrepancy than that of the random sequences, quasi-Monte Carlo methods are expected to provide more accurate results. However, the integrand of the rendering equation is discontinuous where the discontinuity is not aligned with the coordinate axes, thus its variation is in?nite. These discontinuities are usually produced by the projected object boundaries. This property makes the Koksma-Hlawka inequality not appropriate for the error analysis and for the prediction of the convergence rates.

7.2.1 Integrating functions of unbounded variation
In this section the convergence speed is examined for functions which are generally smooth but have general discontinuities of ?nite “length”. First the domain of the integration is assumed to be 2-dimensional, then the results will be generalized to arbitrary dimensions.

discontinuity line
1/ N

domain of discontinuity one sample point in each cell
1/ N

grid lines

Figure 7.1: A typical integrand of the rendering equation

Suppose that N samples have been generated to estimate the integral of a function such as in ?gure 7.1 using a low-discrepancy sequence. In order to overcome the dif?culty that the integrand f has in?nite variation, the ~ ^ function is decomposed into two functions, one is smooth f having continuous mixed derivatives and the other f inherits the discontinuity of f (?gure 7.2).
f ~ f ^ f

=

+

~ ^ Figure 7.2: Decomposition of f into a smooth (f ) and a discontinuous (f ) function
Low-discrepancy sequences generate points in a cellular grid in a way that the difference of the number of points in two cells is at most 1. If there are already N number of points, the size of a cell on the ?nest, ?lled level p p ^ is approximately 1= N  1= N . Let us de?ne the domain of f as the set of those cells that are intersected by the discontinuity. This domain is called the domain of discontinuity. The number of such cells is in the order of the “digital length” of the discontinuity curve, which is the product of the maximum extent l and the resolution of the p p grid N . Since each cell has at least 1 (and at most 2) points, the number of points in this domain is at least l N . The error of quadrature is as follows:

j

Z

z2[0;1]2

f (z) dz

N 1 X f (z )j  j i N i=1

Z

z2[0;1]2

Z N N ~(z) dz 1 X f~(zi )j + j ^(z) dz 1 X f^(zi )j: f f N i=1 N i=1 z2[0;1]2
z z z

(7.1)

~ ~ Since f has ?nite variation, the ?rst term in the error is bounded by VHK (f )  D ( 1 ; 2 ; : : : N ).

7.2. QUASI-MONTE CARLO QUADRATURE FOR THE RENDERING EQUATION

55

^ Concerning the second term, the integration of f is estimated taking l N uniformly distributed samples and averaging the result. Since the samples and the discontinuity are not related in any way, we can suppose that this is a normal Monte-Carlo integration [PFTV92]. The uniform property of low-discrepancy sequence guarantees that this pseudo-random set can be assumed to have uniform distribution. If f is the difference between the maximum and minimum values in the domain of discontinuity, then  2  (f )2 . In our case the number of sample points p p N 0 is l N and the size of the domain V is l= N , thus we obtain with 0.997 con?dence level:

p

Z

V

N N p 1 X V X  f^(z) dz = N 0  f^(zi )  3  V  p f 0 = N 0  f^(zi )  3  f  l  N 3=4 : N
0 0

i=1

i=1

(7.2)

^ Taking into account that f is zero outside the domain of discontinuity, equality 7.1 can be expressed as: 1 f (z) dz N f (zi )j  VHK (f~)  D (z1 ; z2 ; : : : zN ) + 3  f  l  N 3=4 : (7.3) i=1 2 z2[0;1] For large N values the second term will be dominant, which results in O(N 3=4 ) error bound. This is poorer than the O(log2 N=N ) bound suggested by the Koksma-Hlawka inequality assuming, for example, the application

j

Z

N X

p

of the Halton sequence. Note that the pointp from where the second term dominates the ?rst one depends on “intensity” f and size of the discontinuity l. The same analysis can be carried out in higher dimensions as well. In s dimensions a discontinuity of size l would intersect V = l  N 1=s volume of cells which would contain N 0 = l  N (s 1)=s sample points. Thus the general error bound is:

j

Z

z2[0;1]s

f (z) dz

N 1 X f (z )j  V (f )  D (z ; z ; : : : z ) + 3  f  pl  N i HK ~ 1 2 N N i=1

(s+1) 2s :

(7.4)

Thus, for quasi-Monte Carlo integration of discontinuous functions, the order of the error bound will be in between the O(N (1 ) ) bound of ?nite variation functions and the O(N 0:5 ) bound of Monte-Carlo quadrature. The higher the dimension of the integral, the closer the Monte-Carlo and quasi-Monte Carlo techniques get in convergence speed. Thus it is still worth using quasi-Monte Carlo quadrature if the size of the discontinuity l is not very large, since in this case the error could be signi?cantly less then for Monte-Carlo quadrature. Numerical evidence using simple functions In order to demonstrate the previous results, the convergences of a 2-dimensional and a 3-dimensional functions are examined, that are simple enough to analytically compute their integrals. The 2-dimensional function is:

1 2a f2 (x; y) = (x + y)  a +otherwise ; if x + y > 1; (7.5) (x + y) a where a is a free parameter in the range of [0; 0:5]. Note that by setting a appropriately, the intensity of the
discontinuity can be controlled without altering either the value of the integral or the variation of the continuous part. If a = 0:5, then the function has ?nite variation, otherwise it has in?nite variation. The results of the simulation are shown in ?gure 7.3. This ?gure shows the maximum error after a given number of samples. The 3-dimensional function is: (7.6) (x + y + z )  a otherwise ; where a is a free parameter in the range of [0; 1=3]. If a = 1=3, then f3 has ?nite variation, otherwise it has not. The error of integration of f3 is summarized in ?gure 7.3. Numerical evidence for the rendering equation The ef?ciency of Monte-Carlo and quasi-Monte Carlo quadratures have been tested for the presented spherical scene (section 4.4.2) assuming a single pixel camera. The error has been measured separately for the different bounces. Looking at the error measurements of ?gure 7.4, we can see that even for integrands of in?nite variation, quasiMonte Carlo methods are still better but they lose their advantage when computing higher bounces as predicted by



 (x + y + z)  a + 0:6 1:8  a if x + y + z > 1; f3 (x; y; z ) =

7.2. QUASI-MONTE CARLO QUADRATURE FOR THE RENDERING EQUATION

56

Errors of integrating function f2 1 MC: a=0 (infinite variation) QMC: a=0 (infinite variation) QMC: a=0.2 (infinite variation) QMC: a=0.5 (finite variation) 0.1 0.1 1

Errors of integrating function f3 MC: a=1/3 (finite variation) QMC: a=0 (infinite variation) QMC: a=0.1 (infinite variation) QMC: a=1/3 (finite variation)

0.01

0.01

0.001

0.001

0.0001

0.0001

1e-05

1e-05

1e-06 100

1000

10000

100000

1e+06

1e-06 100

1000

10000

100000

1e+06

Figure 7.3: Error of integrating f2 (left) and f3 (right)

Error of single-ray based random walk in the reference sphere (D=1, light=25%) 1 Halton random 1

Error of single-ray based random walk in the reference sphere (D=10, light=25%) Halton random

L1 error

0.1

L1 error 1 10 samples 100 1000

0.1

0.01

0.01 1 10 samples 100 1000

Figure 7.4: Error measurements for 1 and 10 bounces in the spherical reference scene (section 4.4.2) where the BRDF is diffuse, the albedo is 0.5, and 25 percents of the area is a diffuse lightsource

7.3. IMPORTANCE SAMPLING FOR THE RENDERING EQUATION

57

the theoretical results. The other important problem in higher dimensions is that although a low-discrepancy series has almost linearly decreasing discrepancy in the asymptotic sense, this discrepancy can still be high for not very many points (in the solution of the rendering equation we rarely use more than 1000 samples for the estimation of a single pixel). In the case of the Halton series, for example, the base of the series strongly affects the initial behavior of the discrepancy. These base numbers are different prime numbers for different dimensions, thus for high-dimensional integrals the base numbers can be quite high, which results in degraded performance.

7.3 Importance sampling for the rendering equation
When solving the rendering equation, usually directional integrals (or surface integrals in other formulation) should be evaluated. These directional integrals have the following form:

T Lin(~ ; !) = Lin(~ ; !0 )  fr (!0 ; ~ ; !)  cos 0 d!0 : x x x



Z

(7.7)

To allow the application of random or low-discrepancy point sets, the integration domain should be transformed to the unit cube or square. To establish such a mapping, ?rst direction ! 0 is expressed by spherical coordinates ; 0 , which converts the directional integral to the following form:

Z



Lin(~ ; !0 )  f x

x r (! 0 ; ~ ; ! )  cos 0 d! 0 =

Z2 Z
=0 0 =0

fr (; 0 )  cos 0  sin 0  Lin d0 d

(7.8)

since d! = sin 0 d0 d. Let us denote fr (; 0 )  cos 0  sin 0 by w(; 0 ), which is the transfer probability density. Now we ?nd a mapping T (; 0 ) = that emphasizes those directions where the integrand is large and projects the domain of the spherical coordinates onto the unit square:

z

Z2 Z
where

=0 0 =0

w(; 0 )Lin (; 0 ) d0 d =

Z

[0;1]2

1 w(T 1 (z))Lin (T 1(z)) dT dz(z) dz = dT 1(z) 1 dz = t(z)





Z

[0;1]2

w(T 1 (z)) Lin(T 1 (z)) dz; t(z)
(7.9)

is the Jacobi determinant of the inverse mapping. If the Jacobi determinant is large, then a small portion of the unit square is mapped onto a large region. Thus sample points that are uniformly distributed in the unit square will be quite rare in these regions. Alternatively, where the Jacobi determinant is small, the sample points will be dense. Considering this, the meaning of t( ) is the density of sample points in the neighborhood of ! = (; 0 ) = T 1( ). This has an illustrative content for the random case. If is uniformly distributed random variable, then the probability density of ! = T 1 ( ) will be t( ). The solution of the rendering equation for a given point (~ ; ! ) requires the evaluation of the following multix dimensional integral (equation (4.5)):

z

z

z

z

z

L(~ ; !) = Le + T Le + T 2 Le + : : : = x

Z

[0;1]2

:::

Z

[0;1]2

Le + w1  Le + w1  w2  Le + : : : dz1 dz2 : : : t t t
1 1 2

(7.10)

This can be estimated by Monte-Carlo or quasi-Monte Carlo quadratures which evaluate the integrand in sample points and average the results. A crucial design decision of such an algorithm is the selection of mappings Ti to have good importance sampling. Using probabilistic approach, it means that the probability of selecting a walk is proportional to its contribution. Following the directions concluded from the Koksma-Hlawka inequality, the mappings should make the integrand ?at — that is of low variation, or constant in the ideal case. Looking at formula (7.10), which is the single multi-dimensional solution of the rendering equation, this decision seems to be hard to made, since there are too many free parameters to control simultaneously. Fortunately, the solution can also be presented in the following recursive form:

L(~ ; !) = Le + x

[0;1]2

Z w1 Z w2 e+ e t1  [L t2  [L + : : :] : : :] dz1 dz2 : : :
[0;1]2

(7.11)

7.3. IMPORTANCE SAMPLING FOR THE RENDERING EQUATION

58

If we could ensure that each of the integrands of the form

[0;1]2

Z wi Z e+ : : :] dzi ti  [L
[0;1]2

is constant (at least approximately), then the integrand of the single multi-dimensional integral will also be constant [SKCP99]. An optimal importance sampling strategy thus requires density ti to be proportional to the product of the incoming illumination Le + : : : and the cosine weighted BRDF wi . Unfortunately, during random walks the incoming non-direct illumination is not known (the random walk is just being done to estimate it). Thus we have to use approximations for which we have three alternatives. Firstly, information about the illumination in the space can be gathered in a preprocessing phase, then this information can be used to obtain probability densities for importance sampling. This is called the global importance sampling. These methods can be classi?ed according to the data structure built in the preprocessing phase. Since the ray-space is 5-dimensional, it is straightforward to apply a 5D adaptive tree [LW96] that is similar to the well-known octree to store radiance information. Jensen proposed the application of the photon-map as the basis of importance sampling [Jen95]. We assigned the power computed in the preprocessing phase to links established between two interacting patches [SKCP98]. The second alternative is using the information gained during previous walks to approximate the illumination. This strategy is called adaptive importance sampling. Adaptive importance sampling methods neither require the non-uniform probability densities to be constructed in advance, nor simplify them to take into account only local properties, but converge to a desired probability density using the knowledge of previous samples. Three techniques are particularly important, which have also been used in rendering: genetic algorithms [LB94] the Metropolis sampling [MRR+ 53, VG97] and the VEGAS method [Lep80, SK98a]. The ?rst use of Metropolis sampling in rendering aimed at speeding up bi-directional path tracing [VG97]. In the third alternative, the problem is simpli?ed and the indirect illumination is not considered in importance sampling. When the directions are generated, we use only transfer probability density wi and Le representing the direct illumination of the actual point. This is called the local importance sampling. It turns out that we have to encounter severe problems when we have to ?nd a mapping which has a density that is proportional to the product of the effects of the transfer probability density and the direct lighting. Consequently, local importance sampling strategies usually use only either wi or Le to identify important directions. The ?rst alternative is called the BRDF sampling, while the second is called the lightsource sampling.

R

7.3.1 BRDF sampling

ti / wi = fr (!in ; ~ ; !out)  cos  sin : x (7.12) In gathering algorithms !out is known,  is the angle between !in and the surface normal, and !in should be determined. In shooting algorithms, on the other hand, !in is known,  is the angle between !out and the surface normal, and !out should be determined. Due to the fact that ti represents density (probability density for Monte-Carlo methods), its integral is 1. Thus
for gathering walks and for non-transparent materials, the ratio of proportionality in equation (7.12) is

BRDF based importance sampling means that at step density wi , that is

i the density ti is proportional to the transfer probability

Z


H

w d!in =

Z


H

fr (!in; ~ ; !out)  cos in d!in = a(~ ; !out ) x x

where a(~ ; !out ) is the albedo of the surface at point ~ in the outgoing direction. Similarly, the proportionality x x ratio for shooting walks is

Z

Thus the weights w


H

w d!out =

Z


H

fr (!in ; ~ ; !out)  cos out d!out = a(~ ; !in ): x x

= wi =ti are the albedos at the visited points.

BRDF sampling for diffuse materials Diffuse materials have constant BRDF, that is

ti (; ) / wi = fr  cos  sin :

7.3. IMPORTANCE SAMPLING FOR THE RENDERING EQUATION

59

The proportionality ratio is found to make ti to integrate to 1:

ti (; ) =

w R wi d! = 2 =2fr  cos  sin  = cos sin  : R R f  cos  sin  dd i
H r
=0 =0

ti (; ) = 21  [2 cos  sin ] ;  where 1=(2 ) is the probability density of  and 2 cos  sin  = sin 2 is the probability density of .
The corresponding probability distribution functions are:

Assuming that the random variables used to produce the two coordinate samples are independent, this density is obtained in a product form: (7.13)

P () =

Z
0

1  2 d = 2 ;

P () = sin 2 d = sin2 :
0

Z

Thus the required  and  random variables can be found by the following transformations of u; v variables that are uniformly distributed in [0; 1] (section 6.3.1):

 = 2  u;  = arcsin v:
The transformed weight after importance sampling is the albedo

p

wi = wi = fr   = a: t
i
BRDF sampling for specular materials Specular materials can be characterized by the reciprocal version of the Phong’s BRDF model, that is

(7.14)

fr (!in ; ~ ; !out ) = ks  cosn  (=2 ); x where is the angle between !out and the mirror direction of !in onto the surface normal, which will be referred to as !r , and (=2 ) indicates that the outgoing direction cannot point into the object, i.e. the angle  between
the surface normal and the outgoing direction should be less than 90 degrees.
N R V plane perpendicular to R

ψ

φ

surface

reference direction on the plane perpendicular to R
Figure 7.5: Parameterization for the calculation of the albedo

In order to appropriately parameterize the directional sphere, now the north pole is selected by the re?ection direction !r (?gure 7.5). Let us identify a direction by an angle from !r , that is by , and by another angle  between its projection onto a plane perpendicular to !r and an arbitrary vector on this plane. BRDF sampling requires a density which satis?es the following criterion:

ti (; ) / wi = ks  cosn  cos ( ; )  (=2 ( ; ))  sin :

7.3. IMPORTANCE SAMPLING FOR THE RENDERING EQUATION

60

Unfortunately, the cos   (=2 density that is proportional only to wi ~ to 1:

) factor forbids the symbolic integration of this formula, thus we will use a = ks  cosn sin . The proportionality ratio is found to make ti to integrate n + ti (; ) = 2 =2ks  cos sin = n2 1 cosn sin : R R k  cosn sin d d s
=0 =0

Assuming that the random variables used for producing the two coordinate samples are independent, this density is obtained in a product form:

where 1=(2 ) is the probability density of  and (n + 1) cosn The corresponding probability distribution functions are:

ti (; ) = 21  [(n + 1) cosn sin ]; 
sin
is the probability density of .

(7.15)

P () = 2 ; 

P ( ) = (n + 1) cosn sin d = 1 cosn+1 :
0

Z

Thus the required  and  random variables can be found by the following transformations of u; v variables that are uniformly distributed in [0; 1]:

 = 2  u;

= arccos(1 v)1=(n+1) :

The transformed weight after importance sampling is

2 wi = wi = nks  cos ( ; )  (=2 ( ; )): ti +1

(7.16)

Different other specular BRDF models are presented and their BRDF sampling is discussed in [War92, NNSK98b, NNSK98a, NNSK99a].

7.3.2 Lightsource sampling
Lightsource sampling is used in direct lightsource calculations [SWZ96] and as a complementary sampling strategy to BRDF sampling in random walks. Since in this case, the samples are selected from a lightsource instead of the directional sphere, the surface integral form of the transport operator is needed:

(T Le )(~ ; !) = x

Z



Le (h(~ ; x

!0 ); !0 )fr (!0; ~ ; !)cos 0 d!0 = x

Z

Se

x y y x y Le (~; !~~ )fr (!~~ ; ~ ; !) cosj~~  cos ~ v(~; ~ ) d~; y y x y x x x ~j2 y
(7.17)

0

where v (~ ; ~ ) is 1 if points ~ and ~ are not occluded from each other and 0 otherwise, and Se is the surface of yx x y non-zero emission. To obtain a Monte-Carlo estimate for this integral, N points ~1 ; : : : ~N are sampled uniformly y y on the lightsource and the following formula is used:

(T Le )(~ ; !)  jSe j  x N

N X i=1

 y Le (~i ; !~i !~ )  v(~i ; ~ )  fr (!~i !~ ; ~ ; !)  cosj~ i  cos2~i : y y x y x y x x x ~j y
i

0

(7.18)

If the scene has a single homogeneous lightsource which is relatively small and is far from the considered point, then the integrand will be approximately constant on the lightsource surface, thus this estimator has low variance.

7.3.3 Sampling the lightsources in gathering random walks
Since lightsource sampling generates samples only on the direct lightsources, it completely ignores indirect illumination. Thus it cannot be used alone in global illumination algorithms, but only as a complementary part of, for example, BRDF sampling. The simplest way to combine the two strategies is to generate all but the last directions of the gathering walk by sampling the BRDF and to compute the last direction by sampling the lightsource. Note that when stopping the walk, the indirect illumination is assumed to be zero, thus following the directions of the lightsources is a reasonable approach.

7.3. IMPORTANCE SAMPLING FOR THE RENDERING EQUATION

61

Another combination strategy is to trace one or more shadow rays from each visited point of the walk towards the lightsources, not only from the last of them. Formally, this approach can be presented as a restructuring of the Neumann series

L = Le + T Le + T 2 Le + T 3 Le : : : = Le + (T Le ) + T (T Le ) + T 2 (T Le ) : : : (7.19) and using lightsource sampling for the Le = (T Le ) integral while sampling the BRDFs when evaluating the T i Le integrals. Practically it means that having hit a surface, one or more shadow rays are traced towards the

lightsources and the re?ection of the illumination of this point is estimated. This re?ection is used as if it were the emission of the surface. This method is particularly ef?cient if the scene consists of point lightsources. Tracing a single ray to each point lightsource, the illumination due to the point lightsources can be determined exactly (with zero variance). If the scene contains small and large lightsources, then lightsource sampling can be used only for the smaller lightsources. Formally, the unknown radiance L is decomposed into two terms:

L = Lep + Lnp Lep + Lnp = Le + T (Lep + Lnp ):
Expressing Lnp we obtain: Introducing the new lightsource term

(7.20)

where Lep is the emission of the small, point-like lightsources, Lnp is the emission of the large area lightsources and the re?ected radiance. Substituting this into the rendering equation we have: (7.21) (7.22)

Lnp = (Le Lep + T Lep ) + T Lnp :

Le = Le Lep + T Lep (7.23) which just replaces the point lightsources (Lep ) by their effect (T Lep ), the equation for Lnp is similar to the
original rendering equation:

Lnp = Le + T Lnp:

(7.24)

When this is solved, small lightsources are replaced by their single re?ection, then the solution is added to the radiance of the small lightsources (equation 7.20).

7.3.4 Importance sampling in colored scenes
So far, we have assumed that the weights containing the BRDFs and the emissions are scalars thus the densities can be made proportional to them. This is only true if the rendering equation is solved on a single wavelength. However, if color images are needed, the rendering equation should be solved on several (at least 3) different wavelengths. If the different wavelengths are handled completely independently, then the proposed importance sampling strategy can be used without any modi?cations. However, this approach carries out geometric calculations, such as tracing rays, independently and redundantly for different wavelengths, thus it cannot be recommended. A better approach is using rays that transport light on all wavelengths simultaneously. In this case the emission and the BRDF can be represented by vectors, thus to allow importance sampling, we need a scalar importance function I that is large when the elements in the vector are large and small when the elements are small. The importance is a functional of the spectrum. A straightforward way is using the luminance of the spectrum since it emphasizes those wavelengths to which the eye is more sensitive.

7.3.5 Multiple importance sampling
So far, we mentioned two basic importance sampling strategies, the BRDF sampling and the lightsource sampling, which are local in the sense that they focus on a single re?ection. It is easy to imagine that if the sampling considers simultaneously many re?ections, then the number of possible strategies increases dramatically. Obviously, we desire to use the best sampling strategy. Unfortunately the performance of a sampling strategy depends on the properties of the scene, which is usually not known a-priori, thus the best strategy cannot be selected. Instead of selecting the best, Veach and Guibas [VG95] proposed to combine several strategies in a way that the strengths of the individual sampling methods are preserved. Suppose that we can use n different sampling techniques for generating random paths, where the distribution of the samples is constructed from several p1 ; :::; pn importance sampling distributions. The number of samples

7.4. HANDLING INFINITE-DIMENSIONAL INTEGRALS

62

taken from pi is denoted by Mi , and the total number of samples by M = i Mi . The Mi values are ?xed in advance before any samples are taken. The “average probability density” of selecting the sample is then

P

p(z) = ^

n X Mi i=1

z

M  pi (z):

(7.25)

Thus the integral quadrature using these samples is

where i;j is the j th sample taken from the ith distribution, and the weights are

z

Z f (z) Mi n Mi n 1 X X f (zi;j ) = X 1 X w (z )  f (zi;j ) f (z) dz = ^ p(z)  p(z) dz  M i=1 j=1 p(zi;j ) i=1 Mi j=1 i i;j pi (zi;j ) ^ ^ [0;1]s [0;1]s
 wi (z) = PnMiMpi (z) (z) : p
k=1 k k

Z

(7.26)

(7.27)

Let us interpret this result when all methods use the same number of samples. pi ( ) is the probability that a sample is generated by method i. The samples are combined with this weight which guarantees that no sample will be accounted for twices. In order to have an unbiased estimation, i wi ( ) = 1 should hold for all .

z

P

z

z

z

7.4 Handling in?nite-dimensional integrals
Expansion, or walk methods require the evaluation of a series of integrals where the dimension goes to in?nity. One way of attacking the problem is truncating the Neumann series, but this introduces some bias, which can be quite high if the scene is highly re?ective. Fortunately, there is another approach that solves the in?nite-dimensional integration problem through randomization. In the context of Monte-Carlo integration, this approach is called the Russian roulette [AK90], but here a somewhat more general treatment is also given that can also justify this approach for quasi-Monte Carlo quadratures.

7.4.1 Russian roulette
Note that the Neumann series contains a sequence of the following integrals:

L=

[0;1]2

Z w Z    [Le + : : :] dz = w (z)  Lin(z) dz = E w  Lin : t
[0;1]2

A Monte-Carlo quadrature would generate random samples in the domain and estimate the integral as an average of the integrand at these samples. Let us further randomize this computation and before each sample let us decide randomly with probability s whether we really evaluate the integrand at the sample point or simply assume that the integrand is zero without any calculations. In order to compensate the not computed terms, when the integrand is really computed, it is divided by probability s. This randomization introduces a new random variable Lref which is equal to w  Lin =s if the integrand is evaluated and zero otherwise. The Monte-Carlo quadrature, which provides the estimate as an expected value will still be correct:

E [Lref ] = s  E Lref j sample is used + (1 s)  E Lref j sample is not used = sE









 w  Lin 
s

+ (1 s)  0 = E w  Lin = L:





(7.28)

The variance of the new estimator, on the other hand, is increased:

D2 [Lref ] = E [(Lref )2 ]

E 2 [Lref ] = s  E

"

1

 in 2 2  in s 1  E [(w  L ) ] + D [w  L ]:



Lref s

2#

+ (1 s)  0 E 2 [Lref ] =
(7.29)

7.4. HANDLING INFINITE-DIMENSIONAL INTEGRALS

63

7.4.2 Russian roulette in quasi-Monte Carlo quadrature
In the context of quasi-Monte Carlo integration Russian roulette should be explained differently because there is no “randomization” in deterministic techniques. This discussion is also used to generalize the basic concept and to show that termination decision can be made using the complete previous path not only the actual state. Instead of randomization, we can suppose that the domain of integration is extended by additional variables on which a contribution indicator function is de?ned that will determine whether or not a higher order term is included in the quadrature (interestingly, in order to get rid of high-dimensional integrals, we increase the dimension of the integration). In order to compensate the missing terms in the integral quadrature, the really computed terms are multiplied by an appropriate factor. If the used contribution indicator is such that the domain where it is nonzero shrinks quickly, then the possibility of requiring very high dimensional integrals is rather low, which saves computation time. However, the integral quadrature will still be correct asymptotically. A term of the Neumann series has generally the following form

L = : : : W (z1 ; : : : ; zn )  Le (z1 ; : : : zn ) dz1 : : : dzn ;

Z

Z

(7.30)

where W ( 1 ; : : : ; n ) is the product of the weights w1 ( 1 )  : : :  wn ( n ) including the cosine functions of the   angles and the BRDFs, or the product of the ratios of weights and densities w1  : : :  wn = w1 =t1  : : :  wn =tn , depending whether or not BRDF sampling has already been applied. Let us extend this by a contribution indicator function C ( 1 ; r1 ; : : : n ; rn ) that is constant 1 if a sample 1 ; : : : n should be taken into account in the integral quadrature and 0 if it should not and its contribution is assumed to be zero. Usually this function is in separable form

z

z

z

z

z

z

z

z

C (z1 ; r1 ; : : : zn ; rn ) =

n Y

where ci ( i ; ri ) = 1 means that the walk must be continued at step i and ci ( i ; ri ) = 0 forces the walk to stop. Function ci ( i ; ri ) can be de?ned, for instance, by a new weight function  ( i ) in the following way:

z

i=1

ci (zi ; ri );

z

ci (zi ; ri ) =

( 1 if (zi ) > ri ,
0
otherwise.

z

z

The “possibility” of really computing a walk 1 ; : : : n is

P (z1 ; : : : zn ) =

Z1

z

We can de?ne the following function of variables r1 ; : : : ; rn ,

r1 =0

:::

z Z1

rn =0

C (z1 ; r1 ; : : : zn ; rn ) dr1 : : : drn :

~ ~ Lr (r1 ; : : : ; rn ) = : : : C (z1 ; r1 ; : : : zn ; rn )  W  Le dz1 : : : dzn ; ~ ~ where W and Le are appropriate modi?cations of W and Le , which can compensate the missing terms. The integral of this function is

Z

Z

(7.31)

Z1

r1 =0

:::

Z1

rn =0

In (r1 ; : : : ; rn ) dr1 : : : drn =

Z1

~ ~ : : : P (z1 ; : : : ; zn )  W  Le dz1 : : : dzn : (7.32) A suf?cient requirement for this integral to be equal to the original integral L is ~ ~ P (z1 ; : : : ; zn )  W  Le = W  Le : (7.33) ~ ~ There are many possible selections of the contribution indicator and the W and Le functions, that can satisfy
this requirement, thus there are many different unbiased estimators. A widely used selection is letting

Z

r1 =0

:::

Z1 Z

Z

rn =0

~ ~ : : : C (z1 ; r1 ; : : : zn ; rn )  W  Le dz1 : : : dzn dr1 : : : drn =

Z

with the probability of the albedo.

~ ~ W = 1; Le = Le and i (zi ) = wi (zi ): which corresponds to continuing the walk after step i only if w(zi ) > ri . If importance sampling has already been applied, then the walk is continued if w > ri . Since w approximates the albedo, this strategy continues the walk

7.4. HANDLING INFINITE-DIMENSIONAL INTEGRALS

64

BRDF sampling for materials of multiple re?ection type Practical re?ection models incorporate different simple BRDFs. For example, a lot of materials can be well modeled by a sum of diffuse and specular re?ections. So far, methods have been presented that are good for either the diffuse or the specular re?ection, but not for the sum of them. Fortunately, Russian-roulette can also be extended to handle these cases. If the re?ection model is a sum of different BRDFs, then a random selection can be made from the different components. Suppose that the transfer probability density is available in the form of a sum of the weights corresponding to elementary BRDFs:

w = w1 + w2 + : : : + wn :
Thus the radiance of a single re?ection is:

L = w  Lin d! = w1  Lin d! + : : : + wn  Lin d!:




Suppose that mappings ti can be found for each integral, that also mimics important directions:

Z

Z

Z

L=

[0;1]2

Z wn Z w1        Lin dz + : : : +  Lin dz = E w1  Lin + : : : + E wn  Lin : t1 tn
[0;1]2

Let us select the ith BRDF with probability pi and weight the resulting radiance by 1=pi or stop the walk with ^  probability p0 = 1 p1 : : : pn . Thus the new random variable L is wi  Lin =pi if the ith model is used, and 0 ^ if no model is selected. The expected value of L will still be correct:

 in  in     ^ E [L] = p1 E w1 p L +: : :+pn E wnp L +(1 p1 : : : pn)0 = E (w1 + : : : + wn )Lin = L: 1 n









(7.34)

This is a Monte-Carlo estimation of the sum. According to importance sampling, the variance will be small   if wi Lin =pi can be made nearly constant. Since we usually do not have a-priori information about Lin , wi =pi can be made a constant number. Thus to obtain a low-variance estimator, an elementary BRDF should be selected   with the probability of its transformed weight wi . Note that the weight wi may either be equal or approximate the albedo, thus a low-variance estimator selects an elementary BRDF with the probability of its albedo. In order to generate an “out” direction from the “in” direction and the surface normal, the following general BRDF sampling algor


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

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

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