Light-Driven Global Illumination with a Wavelet Representation of Light Transport by Robert R. Lewis B. S. Physics, Haxvey Mudd College, 1974 M. A. Astronomy, University of California at Berkeley, 1979 M. S. Computer Science and Engineering, Oregon Graduate Institute, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doc to r of Ph i losophy in THE FACULTY OF GRADUATE STUDIES (Department of Computer Science) we accept this thesis as conforming to the required standard The University of British Columbia February 1998 © Robert R. Lewis, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date 6, DE-6 (2/88) Abstract This thesis considers the problem of global illumination: the modelling of light as it travels through a scene interacting with the objects contained within the scene. Starting with a description of the problem and a discussion of previous work, we explore a new approach called light-driven global illumination that offers several advantages over its predecessors: a lower asymptotic complexity, a wider range of representable surface interaction phenomena, and an absence of the need for "meshing" - object surface subdivision needed primarily to represent shadows. Light-driven global illumination is intermediate between local and global illumination. Representing light with wavelet basis functions, we are able to treat both the interaction between two surfaces and the interaction of a surface with a radiation field in a source-to-destination model that applies to whole surfaces, not just small elements. We have found this "wavelet radiative transfer" to be a valid way to generate and store complex global light field data as four-dimensional textures for incorpora-tion in local illumination solutions. Wavelets can considerably reduce the otherwise substantial storage and reconstruction problems associated with doing this. We include several examples of this. We also discuss plausible illumination models, which are required to make ii light-driven global illumination work theoretically. Like wavelet radiative transfer, these models have application in other areas of rendering besides global illumination. Finally, we develop the theory behind light-driven global illumination and apply it successfully to some simple examples. While we find the algorithm to be quite slow compared to other well-known rendering algorithms, we analyze what is needed to make it competitive. In conclusion, we find that representing light with wavelets has a set of advantages that are independent of the comparative inefficiency of the light-driven algorithm. . . 111 Contents Abstract ii Contents iv List of Tables xi List of Figures xii Acknowledgements xvi Dedication xvii 1 Introduction 1 1.1 Conventions . . . . . 2 2 Illumination 7 .2.1 What is Illumination? 7 2.2 Relation to Computer Graphics 7 2.3 Local and Global Illumination 8 2.4 Radiative Transport Theory 10 2.4.1 Radiance and the Transport Operator T 10 iv 2.4.2 The Blocking Operator B • 12 2.4.3 Reflection, Transmission, and the Surface Interaction Opera-tor S 13 2.5 Renderers 16 2.5.1 Luminaire Models 17 2.5.1.1 Point Luminaire 17 2.5.1.2 Directional Luminaire 18 2.5.1.3 Other Luminaires . 18 2.5.1.4 What Constitutes a Luminaire? . . 19 2.5.2 Propagation Models : : 19 2.5.3 Object Models 21 2.5.3.1 Polygons 21 2.5.3.2 Parametric Surfaces 21 2.5.3.3 Implicit Surfaces 22 2.5.3.4 Constructive Solid Geometry 23 2.5.4 Illumination Models 23 2.5.4.1 Phong Illumination Models 23 2.5.4.2 Torrance-Sparrow Illumination Models 27 2.5.4.3 Neumann-Neumann Illumination Models 28 2.5.4.4 Minnaert Illumination Models 31 2.5.5 Surface Models . . 32 2.5.6 Shading Models . 33 2.5.7; Rendering Techniques 34 2.5.7.1 The Rendering Equation 34 2.5.7.2 Depth-Buffering 36 v 2.5.7.3 Raytracing 36 2.5.7.4 Radiosity : 37 2.5.7.5 Photon. Maps 39 2.5.7.6 Comparison and Summary 39 3 The Lucifer Algorithm 41 3.1 Isolation 41 3.2 Power Computation for Isolated Volumes 42 3.3 Solution Order and Convergence 43 3.4 Spatial Partitioning 44 3.5 The Lucifer Algorithm 44 3.6 Complexity 47 3.7 Light Through a Window 51 3.8 Radiance Representation 54 3.8.1 Box Discretization . 54 3.8.2 Wavelets 55 4 Making Illumination Models More Physically Plausible 58 4.1 Energy Conservation 59 4.1.1 Making Phong Illumination Models Conserve Energy 62 4.1.2 Do Torrance-Sparrow Illumination Models Conserve Energy? 64 4.1.3 Making Neumann-Neumann and Minnaert Illumination Models Conserve Energy 67 4.2 Making Illumination Models Reciprocal 67 4.2.1 Are Phong Illumination Models Reciprocal? 67 4.2.2 Are Torrance-Sparrow Illumination Models Reciprocal? . . . 68 vi 4.2.3 Are Separable Illumination Models Reciprocal? 68 4.3 An Energy-Conserving, Reciprocal Illumination Model 69 4.4 Summary of Plausibility Results . . 71 5 Wavelet Radiative Transport and Surface Interaction 73 5.1 Radiance in Nusselt Coordinates 74 5.2 Radiance Representation : 76 5.3 Multidimensional Wavelets 76 5.4 Irradiance . 79 5.5 Power Flux : • • 79 5.6 Transport 80 5.7 Surface Interaction •'. 81 6 Wavelet Radiative Transfer Implementation and Practicum 85 6.1 Naming Conventions 86 6.2 Indexing Wavelets . 87 6.3 The Waveletlndex Class 88 6.4 Storing Wavelet Coefficients 89 6.4.1 Flashing Coefficients 89 6.4.2 Multichannel Grouping 89 6.4.3 Hashing Nodes Instead of Coefficients . 90 6.5 The WaveletCoefficientTree Class 91 6.6 Representing Transport Geometry 92 6.6.1 Rectilinear 93 6.6.2 Perspective •". 94 6.6.3 Bilinear 94 vii 6.7 The TransportGeometry Class . 96 6.8 Choice of Wavelet 96 6.9 Problems with Transport Coefficient Computation 97 6.9.1 Some Source Points Do Not Project Into Destination Space . 97 6.9.2 The Projected Hypercube Has Curved Sides 98 6.10 Integration Techniques 98 6.10.1 Reducing the Dimensionality 99 6.10.2 Numerical Quadrature 100 6.10.3 Comparison 101 6.11 Functional Decomposition of WRT . . . . 104 6.11:1 tg.oracleQ . . • 104 6.11.2 tcJntegrateQ .' . 105 6.11.3 wn_pull() and wnjpushQ . . 105 6.11.4 tcevalQ 106 6.11.5 tc-propagate() 106 6.11.6 wet-transport() .' 107 6.12 Example of Transport 108 6.13 Example of Radiance Compression 110 Lucifer Implementation and Practicum 114 7.1 The Cell Class 114 7.2 Bidirectional Reflectance and Transmittance Distribution Functions 116 7.3 Implementation Choices 116 7.3.1 Integration Scheme 117 7.3.2 Storing Transport Coefficients 117 7.3.3 Luminaire Model '. 118 viii 7.3.4 Object Model 119 7.3.5 Illumination Model 119 7.3.6 Shading Model 119 7.3.7 Surface Model . '. . 1 2 0 7.4 Surface Interaction with Haar Wavelets • • • 120 7.5 Functional Decomposition of Lucifer 122 7.5.1 wct-transportConfined() 122 7.5.2 wct_transportRestricted() 123 7.5.3 wLblocking() 123 7.5.4 cLcfcTransblockQ 124 7.5.5 wctJransportWallToObjQ . 126 7.5.6 wct-addInteraction() . . . . . . 126 7.5.7 wctJransportObjToWall() - 127 7.5.8 cLtransport() 128 7.5.9 cLtransblockQ 128 7.5.10 cLbalanceQ 129 7.6 Results . .' 129 7.7 Analysis . .- . 134 8 Conclusions and Future Work 137 8.1 Fitting Plausible Shaders to Physical Data 139 8.2 Fitting Plausible Shaders to Physical Data 139 8.3 Curved Surfaces . 139 8.4 Importance 139 8.5 Unifying Illumination Modelling and Surface Modelling 140 ix Appendix A One-Dimensional Wavelet Properties 141 A . l Scaling Functions and Wavelets 141 A.2 Multiresolution Refinement Equations 143 A.3 Orthogonal Wavelets 144 A.4 Biorthogonal Wavelets • . . 144 A.5 Wavelet Projections and Approximation 145 A.6 The Fast Wavelet Transform 145 A;7 Wavelet Compression . . 146 Appendix B Pseudocode for Haar Surface Interaction 148 Bibliography 155 Index 163 x List of Tables 1.1 Table of Commonly-Used Symbols. 3 1.2 Light Measurement Terminology 6 2.1 Effects Available with Established Rendering Techniques 40 4.1 Summary of Plausibility Results 70 6.1 Major WRT Classes and Their Abbreviations 86 6.2 Number of Bits Required to Store a Nonstandard Multiresolution Index as a Function of the Maximum Resolution Level . . . . . . . . . 87 • 7.1 Lucifer Test Cases 131 xi List of Figures 2.1 A Typ ica l G loba l I l luminat ion Prob lem 9 2.2 Surface Geometry for Reflection and Transmission 14 2.3 I l luminat ion M o d e l Test Conf igurat ion 25 2.4 Sphere Shaded wi th a Phong I l luminat ion M o d e l 26 2.5 Sphere Shaded wi th Neumann-Neumann and Minnaer t I l luminat ion Models 30 2.6 Compar ison of a Square Shaded wi th a Neumann-Neumann I l lumi-nation M o d e l and wi th a Phong Shader Us ing B l inn 's Ff 30 3.1 A n App l ica t ion of Isolation 42 3.2 Octree Cel l Nomenclature . 45 3.3 The Lucifer A lgor i thm 46 3.4 " B a l l s " M o d e l 51 3.5 A r ray of F i xed Direct ion V iews 52 3.6 A r ray of F i xed Posi t ion Views 53 3.7 Box Discret izat ion Examp le 55 3.8 Or ig ina l Image and Its Reconstruct ion from 4-Dirnensional Wavelet Coefficients 56 xii 4.1 Specular Integrals HS(S) for Phong 's F f and B l inn 's F? 61 4.2 Sphere Shaded wi th an Energy-Conserv ing Phong I l luminat ion M o d e l 61 ; 4.3 Direct ional-Hemispher ical Reflectances (Phong, B l i nn , Torrance-Sparrow) 65 4.4 Direct ional-Hemispher ical Reflectances (Neumann-Neumann, Minnaer t ) 66 4.5 Direct ional-Hemispher ical Reflectance (Reciprocal Phong) 70 4.6 Sphere Shaded wi th a Reciprocal Phong I l luminat ion M o d e l 71 • 5.1 Nusselt Coordinates . 74 6.1 The Waveletlndex t y p e d e f 88 6.2 The WaveletCoefficientTree t y p e d e f 91 6.3 Coord inate Systems 92 6.4 The TransportGeometry t y p e d e f 95 6.5 Relat ive Eucl idean Distance vs. T ime for Different Integration Tech-niques and Numbers of Samples 101 6.6 Relat ive M a x i m u m Depar ture vs. T ime for Different Integration Techniques and Numbers of Samples 102 6.7 Simplif ied Funct ional Decomposi t ion for Wavelet Radia t ive Transfer 104 6.8 tcjpropagate() Pseudocode 106 6.9 wet-transport() Pseudocode 108 6.10 Geometry Used for Examp le of Transport 109 6.11 A Stained Glass W i n d o w 110 6.12 Four-Dimensional Resul ts of Wavelet Radiat ive Transport I l l 6.13 Detai led V iew of the Br ightest Par t of F igure 6.12 112 6.14 Examp le of Wavelet Transpor t 112 6.15 Effect of Compressing the Wavelet Radiance Dis t r ibu t ion 113 xi i i 7.1 The Ce/Z typedef 115 7.2 Simplified Functional Decomposition for Lucifer 122 7.3 Unblocked, Partially Blocked, and Completely Blocked Propagation 123 7.4 cLcfcTransblock() Pseudocode 125 7.5 wctJransportWallToObj() Pseudocode 126 7.6 wctJransportObjToWall() Pseudocode . . . 127 7.7 cLtransport() Pseudocode 127 7.8 cLtransblock() Pseudocode 128 7.9 cLbalance() Pseudocode 129 7.10 Configuration for Lucifer Example 130 7.11 Lucifer Sequence 132 7.12 Lucifer Example 133 7.13 Lucifer C P U Time vs. Z^'x • • 134 7.14 Lucifer Peak Memory Usage vs. Z""x 135 A . l Haar Wavelet and Smoothing Functions 142 A.2 Daubechies-4 Wavelet and Smoothing Functions 143 A.3 Linear Spline Wavelet and Smoothing Functions 143 A. 4 The Fast Wavelet Transform 146 B. l Pseudocode for Surface Interaction I: top level 149 B.2 Pseudocode for Surface Interaction II: Delta_ks 149 B;3 Pseudocode for Surface Interaction III: DeltaJs 150 B.4 Pseudocode for Surface Interaction IV: Delta_kr 151 B.5 Pseudocode for Surface Interaction V: DeltaJr 152 B.6 Pseudocode for Surface Interaction VI: Delta_x 153 xiv B.7 Pseudocode for Surface Interaction V I I : De l ta .y 154 xv Acknowledgements The author is grateful to John Buchanan, Pat Hanrahan, Paul Lalonde, Pierre Poulin, Leena-Maija Reissell, and, of course, his supervisor, Alain Fournier, for encouragment and numerous discussions on this thesis and related subjects. ROBERT R. LEWIS The University of British Columbia February 1998 xvi For V i c to r i a , the patience, the long hours, and, of course, the love. A n d for A l a i n , EDGAR BERGEN: Mor t imer , I can' t understand how one person can possibly be so stupid! MORTIMER SNERD: Wel l , yuh see M r . Bergen, I got this fella helpin ' me... xv i i Chapter 1 Introduction Realistic images are images which are indistinguishable from images captured from the physical world. The creation of such images is a fundamental goal in computer graphics. Research in the modelling of shape, of properties of surfaces, and of the interaction of light with matter is directed towards this goal, as evidenced in major texts on the subject such as Glassner [24], Foley, et al. [18], and Watt [80]. For the past decade, some of the most realistic images in computer graphics have been produced by global illumination techniques. These techniques are limited, however, in the range of lighting and surface characteristics that they can represent. This thesis presents a new technique in realistic image synthesis called "light-driven global illumination" that allows a much wider range of characteristics to be modelled. After this introduction, we will examine the basics of local and global illumination in Chapter 2. Chapter 3 will present the new global illumination al-gorithm itself. Chapter 4 will present an analysis of the "plausibility" of current illumination models. (As we will see, energy-conserving illumination models, which are a subset of physically plausible illumination models, are a theoretical prereq-1 uisite for our technique.) Chapter 5 explains the approach we use to perform the local solutions required by the global solution, a formulation of radiative transfer and surface interaction that is based on wavelet techniques. Chapters 6 and 7 de-scribes implementations of ideas from the preceding chapters. Chapter 8 describes our conclusions. 1.1 Conventions Table 1.1 lists symbols that will be commonly used throughout this thesis. Equations taken from other sources will be modified to use these symbols. 2 Table 1.1: Table of Commonly-Used Symbols. The "*" is a wildcard matching character. symbol definition angle between N and H B blocked radiative transport operator scaling constants for various illumination models specular "half-angle" dA element of surface area du>i sin 9id8id<f>i, an element of incident solid angle dur sin 9rddrd(f>r, an element of reflected solid angle £>* facet slope distribution functions E irradiance EN normal irradiance F Fresnel factor fr bidirectional reflectance distribution function (BRDF) P* x s specular shading functions <t>i incident azimuthal angle 4>r reflected azimuthal angle ®L angle between incident and viewing directions G global illumination operator G geometrical attenuation factor H bisector of source and viewing directions 1 identity operator continued on next page 3 Table 1.1: Table of Commonly-Used Symbols. The "*" is a wildcard matching character. continued from previous page symbol definition ambient reflectance coefficient kd diffuse reflectance coefficient ks specular reflectance coefficient h* directional-hemispherical reflectances ka fraction of total energy reflected specularly La ambient radiance Lt incident radiance Lr reflected radiance Lt transmitted radiance M exitance m RMS slope of the surface N the surface normal < specular exponents the hemisphere surrounding N the hemisphere surrounding —N reflected incident direction S surface interaction operator s incident direction T radiative transport operator continued on next page 4 Table 1.1: Table of Commonly-Used Symbols. The "*" is a wildcard matching character. continued from previous page symbol definition Oi incident polar angle Or reflected polar angle transmitted polar angle V direction of viewer 5 Physics Radiometry Radiometric Units (SI) flux angular flux density incident flux density emergent flux density radiant energy radiant power radiance irradiance radiosity (or exitance) radiant intensity joules [ J = kgm2s~2 ] watts [ W = Js'1 ] [ Wm^sr'1 ] [ Wm~2 ] [ Wm'2 ] [ Wsr~2 ] Table 1.2: Light Measurement Terminology We have attempted to be compatible with the ANSI/IES standard [36] wher-ever possible. Table 1.2, taken from Cohen and Wallace [14] with minor alterations, lists the equivalent terminology in the physics and radiometry (the science of the physical measurement of electromagnetic energy) communities of the quantities we will be discussing. We will not be considering photometry (the psychophysical measurement of the visual sensation produced by the electromagnetic spectrum), as determining what viewers actually see when looking at a scene. This may be treated as an or-thogonal postprocessing step after an accurate computation of the light that reaches them. The latter is our goal. 6 Chapter 2 Illumination This chapter will describe the theory of illumination that underlies rendering and then goes on to an overview of previous work in global illumination in computer graphics. 2.1 What is Illumination? The study of illumination is the study of light and how it interacts with matter to produce visible scenes. Its principles are derived from those of physics. As we will see in subsequent chapters, the equations describing illumination are relatively easy to write down, but hard to solve in non-trivial cases. 2.2 Relation to Computer Graphics Computer graphics is the use of computers to produce images. Frequently, it is desirable that those images look realistic. "Realistic" in our context means that viewers get the impression that what they are seeing when they look at a computer-7 generated (i.e., rendered) art i fact (photograph, video screen, etc.) is an image of a scene that exists in the physical world outside the computer. It does not mean that they receive the same st imulus that they would in a physical sett ing. Af ter a l l , it is usually quite obvious when one is looking at a photograph or a video image. W h a t we want is for the art i fact to be indist inguishable f rom a photograph or video image of a scene in the physical wor ld: to be realistic. It is therefore product ive to study i l luminat ion for this purpose. 2.3 Local and Global Illumination There are two broad categories of problems that arise when at tempt ing to solve for i l luminat ion: local and global . Loca l i l luminat ion problems are confined to smal l (often infinitesimal) areas or volumes. They usually involve a single light source (which we wi l l refer to as a luminaire), a surface or volume element, and a viewing direct ion. The basic question of local i l luminat ion is: How does the element interact wi th the incident light to produce the light leaving the element in the viewing direction? G loba l i l luminat ion describes how light is distr ibuted in a scene: a collection of objects, including luminaires, immersed in a given medium. F igure 2.1 shows a typical global i l luminat ion problem. G loba l i l luminat ion solutions must consider mult iple reflections. Par t of that solution may invoke local i l luminat ion solutions. Fournier [22] describes their interrelation further. The goal of a global i l luminat ion scheme may vary wi th the needs of the user. In order of increasing scope, global i l luminat ion problems may require solution at... • a small set of points. For example, the user may only need to know if there 8 directional light source Figure 2.1: A Typical Global Illumination Problem is sufficient light for people seated in a library to read by. (These sorts of problems are more often found in illumination engineering than in computer graphics.) all visible surfaces. Given a scene with non-participating media and a viewer at a fixed location, compute the illumination of every surface that the viewer sees. (This is the "viewer dependent" problem shown in Figure 2.1.) all surfaces. Like the previous class, but producing a solution that represents the illumination on every surface, which, for example, allows the user to do a "walkthrough" as a postprocessing step. (This is the "viewer independent" problem.) all points. In the presence of participating media, it is in general necessary to compute illumination not just at surfaces, but throughout the volume of a scene. 9 2.4 Radiative Transport Theory This section will present the physical basis underlying global illumination as mod-elled by computer graphics. It will start with a discussion of radiance and how it relates to other physical quantities, 2.4.1 Radiance and the Transport Operator T Let us first discuss some of the basics of how light is represented. The fundamental quantity is radiance, which is defined in the ANSI/IES standard [36] to be the amount of power passing in a given direction though a given surface per unit area (perpendicular to the direction of travel) per unit solid angle. Radiance at a point P in a direction S is usually represented by a function L(P, S). We note, as the standard does, a distinction between "radiance" and "spec-tral radiance". Spectral radiance adds additional units of inverse wavelength to radiance to allow for a distribution of radiance over wavelength. This thesis (like most of the rendering literature) takes radiance to be mo-nochromatic and assumes one can construct a polychromatic image by combining "channels" of independently computed monochromatic images. We assume further that these channels are "independent": that they do not interact with each other, as would be the case in the presence of phenomena like fluorescence. Also like most of the literature, we ignore polarization and assume time invariance. One very useful property of radiance, or "scene radiance" as it is called, is its proportionality to "image irradiance", the quantity detected by a camera viewing a scene. This relationship is derived in Chapter 10 of Horn [34]. Radiance's most useful property for rendering, however, is its invariance: In a non-participating medium and in the absence of diffraction phenomena (both 10 of which we assume here), the radiance given off at a point P 0 on a surface in a direction S is constant until reaching another surface1. We can express this as the action of a radiative2 transport operator T acting on a source radiance Ls to produce a destination radiance L&: ^(P,S) = T £ S = L S (P 0 (P,S) ,S) (2.1) Since we also assume that the radiation travels in a straight line3, P = P 0 + St (2.2) for some scalar t > 0. This principle underlies raytracing. Technically, we should allow for the light travel time from source to des-tination, but our assumption of time invariance makes this unnecessary. We can therefore treat light travel as instantaneous. It is evident from Equation (2.1) that for any scaling factor a and radiance distribution L, TaL = oTL. (2.3) T is also linear in the following sense: given two sources A and B with radiances LSA and LSB and LDA = TLSA and Lds = TLSB, then the resulting radiance Ld is Ld(P, S) = T {LsA + LsB) = TLsA + TLsB = LdA + LdB (2.4) if, from the point of view of P, neither A nor B obstructs the other. 'This is derived, for instance, in Siegel and Howell [68]. 2For brevity, we may drop the "radiative" throughout this thesis, since we are only concerned with radiative transport, as opposed to conductive, convective, or any other form of energy trans-port. 3 Note that we are not ruling out refraction, which takes place at a surface, not during propaga-tion. 11 Let us review the assumptions we have made about radiance and T . These are all based on the way light behaves under what we consider "ord inary" rendering circumstances: • straight- l ine (ray) propagation • l inearity (in the absence of blocking) • instantaneous transport • energy conservation • steady-state power transfer • non-interact ing monochromat ic channels • absence of diffraction phenomena • absence of polar izat ion phenomena It is impor tant to note that while all of these assumptions are compat i -ble wi th classical electromagnetic theory (under certain circumstances), it is these empir ical ly-based assumptions that drive our discussion, not the theory itself. Other areas of investigation which make similar assumptions (e.g., neutron transport) may find this work appl icable. 2.4.2 The Blocking Operator B We noted in Section 2.4.1 the compl icat ions that can arise in the presence of more than one destinat ion object. We can account for them by the incorporat ion of a blocking operator B defined for an obstruct ing object O: L 6 (P ,S) = B I , = 6(P ,S)I , (P > S) (2.5) 12 where 6(P, S) is a geometric funct ion which is 0 if a ray star t ing on the source object at point P in direct ion S intersects O and 1 otherwise. If we then have a source object S and a destination object D, and if O does not lie "beh ind" D, we can express the result ing radiance Lj, as Ld = TBLS. (2.6) 2.4.3 Reflection, Transmission, and the Surface Interaction Oper-ator s In this section, we wi l l review the formulat ion of the fundamental equation describing physically-based i l luminat ion model l ing, presented below as Equat ion (2.11). Our presentation follows those given in Cook and Torrance [15], Foley, et a l . [18], and Cohen and Wal lace [14]. F igure 2.2 shows the (usual) geometry for two cases of local i l luminat ion. In both of them, an inf ini tesimal element of surface area dA is being i l luminated by an incident radiance L{ coming from an inf initesimal solid angle du{ surrounding direction S . (A l l vectors presented here are are of unit magnitude.) A n observer (or measuring device) is located in direct ion V . N is the surface normal at dA. On the left, the observer is above the "hor izon" of dA in the hemisphere the hemisphere surrounding N . In this case, the observer receives radiat ion reflected from dA. O n the right, the observer is below the horizon in the hemisphere surrounding — N . In this case, the observer receives radiat ion t ransmit ted by dA. A s treatment of these two cases is very simi lar, let us consider the reflection case first. We assume all reflected radiat ion passes unobstructed into We use N as the z-axis of a polar coordinate system so that we can specify 13 Figure 2.2: Surface Geometry for Reflection and Transmission S by the usual incident polar and azimuthal angles 0,- and <fo and V by the usual reflected polar and azimuthal angles 0r and (j>r4. To consider an illumination model from the standpoint of power balance, we start from the equation of power balance given in Siegel and Howell [68]: dE = Li ( N • S) dui (2.7) where dE is the irradiance resulting from the illumination of the patch in direction S subtending the solid angle du{. Using the definition of / r , the bidirectional re-flectance distribution function (hereafter, BRDF) , this irradiance gives rise to dLr, the resulting radiance of the patch: dLr = fr(S,V)dE (2.8) 4 For the time being, we wil l assume 4>, and <j>r are measured from some locally-defined frame of reference. 14 In general, fr is a function of the incident direction (#,-, and the reflected direction (9r, 4>r)5- For the sake of brevity, we will write it as / r (S , V ) . As we have assumed an opaque, nonemissive surface, the only contributions to Lr can come from 0 N . Therefore, we can integrate du>i over fiN to get Lr = [ fr(S,V)Li(N-S)dut (2.9) N We allow for transmission, as shown on the right of Figure 2.2, by considering that radiation can come from anywhere and can go in any direction. This means that we need to generalize the limits of our directional integration to refer not to fi+ and Q^, but to Jl^, the hemisphere (either fi+ or that contains V and fiN, the hemisphere opposite 0£J. We also need to allow for a transmitted component of radiance using the B T D F ft. Equation (2.9) then becomes U= I ft(S,V)Li | N - S | dui (2.10) In general, then, the total radiance for a non-emissive surface is L = Lr + Lt. If the surface is emissive, we can add an additional surface emissivity term Le when the observer is in fiN6. To summarize, then, the total radiance given off by dA is L r ( V ) = Le(V)+ [ / r ( S + , V ) I 8 - ( S + ) | N - S + | dut + / / t ( S - , V ) L t ( S - ) | N - S " | dui (2.11) where Lr is the total radiance given off of a surface with normal N , Le is the surface emissivity, L\ is the incident radiance, is the reflection hemisphere (contains 5 fr may also vary over the surface, but that variation is usually treated as part of texturing rather than illumination modelling, so we will ignore it here. As mentioned above, we are also ignoring its possible dependence on wavelength. 6 0 ^ is somewhat trickier, since if the surface is emissive, then the interior of the object is also likely to be emissive, which would imply a volume emissivity rather than a surface emissivity. As we are ignoring participating media in this section, we will therefore restrict ourselves to surface emissivity only. 15 N , the "v iewing" direct ion V , and the "posit ive source" direction S + ) , fi^ is the transmission hemisphere (opposite Q^, containing the "negative source" direct ion S - ) , fr is the B R D F , ft is the bidirect ional t ransmit tance distr ibut ion funct ion ( B T D F ) , dui is sin QidQidcpi, 0{ is the incident polar angle, and <fo is the incident az imuthal angle. A local i l luminat ion solution is entirely characterized by its B R D F and B T D F functions. Note that we ignore the spatial variat ion of the radiances and dist r ibut ion functions, since it symbol ical ly factors out of Equat ion (2.11), We represent combined reflection and transmission surface interaction by the integral operator <S: Lr = Le + SLi (2.12) 2.5 Renderers In computer graphics, a realistic rendering package (which we wil l refer to hereafter as a Tenderer) is a (predominantly) software system used to create realistic images. It usually does this by performing a simplif ied solut ion of the equations describing light and matter interact ion. It is useful to consider a rendering package as being composed of seven parts: • a luminaire model: how light sources are represented • a propagation model: how light travels f rom source to surface, surface to sur-face, and surface to viewer • an object model: how objects in the scene are represented • an illumination model: how light interacts wi th the surface of an object 16 • a surface model: how small-scale surface structure is represented • a shading model: how light passing to the viewer is converted into an image • a rendering technique: how all the models that constitute the package interact to produce an image All of these are necessary for a rendering package, although in some cases they can be trivial. In the remainder of this section, we will discuss each of these parts and how they relate to each other. 2.5.1 Luminaire Models Models of luminaires, or light sources, can range from the simple to the sophisticated. In this section, we will cover several popular ones. The most important attribute of a luminaire model is its dimensionality: 0 (a point luminaire), 1 (a linear luminaire), or 2 (an area luminaire). 2.5.1.1 Point Luminaire A point luminaire assumes that a finite amount of light is coming from a single point. If we take $ as the amount of emissive power being given off by the luminaire, the radiance of a point luminaire is: ^ = i ^ ( S - S') (2-13) where r is the distance from the luminaire to the illuminated point and S' is the direction from the illuminated point to the luminaire. The Dirac delta function S(x) is defined by /oo f{x)8(x)dx = f{0) (2.14) -oo 17 for any function f(x). 2.5.1.2 Directional Luminaire A directional luminaire is effectively a point luminaire at infinity: every point on the illuminated surface receives light from the same direction. Unlike a point luminaire, however, a directional luminaire is not characterized by its power (which is infinite) and a position, but by its normal irradiance EN, the amount of power per unit area transferred to a surface with normal incidence7, and a fixed direction S'. We can express the radiance of a directional luminaire as LE = ENS{S - S'). (2.15) 2.5.1.3 Other Luminaires Directional and point luminaires were the only ones available in early Tenderers (with the exception of diffuse area luminaires used in radiosity). Verbeck and Greenberg, for example, describe in [77] how to construct approximation of linear and area light sources from point light sources. They also describe how to model the effects of flaps and focussing. In 1990, Poulin and Amanatides [59] were able to model linear luminaires directly. The next year, Tanaka and Takahashi [74] did the same thing for area luminaires. Both of these groups worked with sources that were isotropic in the sense that the emitted radiance was constant with respect to direction for each differential line or area element. In 1992, Ashdown [2] drew on work from the illumination engineering com-munity to incorporate measurements of real, non-isotropic luminaires using "near-7EN is the value of irradiance E when the light is coming from the normal direction N. 18 field photometry". More recently, Lalonde and Fournier [41] was able to efficiently represent an area luminaire using wavelets. 2.5.1.4 What Constitutes a Luminaire? The distinction between a luminaire and any other object is somewhat artificial. Light interacting with a surface produces reflected and transmitted radiance distri-butions which, as far as a T e n d e r e r is concerned, should be indistinguishable from emitted distributions. The technique we will present in this thesis will make no (essential) distinctions between emitted, reflected, and transmitted radiance distri-butions. 2.5.2 Propagation Models Light travels from surface to surface through a medium. The easiest medium to represent is a non-participating medium - a medium that does not in any way impede radiation travelling through it. A vacuum is one such medium. The propagation model in such a medium is the one we used to describe the invariance of radiance under the T operator in Equation (2.1): the radiance remains constant along a ray from source to destination. In the presence of a participating medium, propagation becomes more diffi-cult to represent. There are three phenomena which must be taken into account: emission, absorption, and scattering. We have already mentioned emission from surfaces, but it can also happen from an incandescent gas such as a fire, a gaseous nebula, or the Sun. When the model of emission cannot be confined to a two-dimensional manifold, it must be treated as part of propagation. 19 Absorption removes energy from the beam entirely, converting it to a non-visible form of energy such as heat. Smoke, for instance, can make a very absorptive medium. Scattering does not eliminate energy from the beam, but redirects it. It may remove energy along the line-of-sight from source to destination (outscattering), but it may also scatter energy into the line-of-sight from other sources (inscattering). So, for instance, a street light in fog, a typical scattering medium, appears surrounded by a halo of light that was redirected towards the viewer by the medium. Allowing for participating media turns Equation (2.1) into an integro-differential equation. Taking Glassner's [24] formulation of it and substituting our own notation, we get: where j is the volume emissivity, 0 is the directional sphere, K is the volume outscat-tering or inscattering probability function, cja is absorption function. Equation (2.16) is easy to interpret. The left hand side is the spatial variation of radiance in direction S . On the right hand side, the first term is the volume emissivity adding power (radiance per unit length, actually) to the beam. The . second is the amount of power scattered into the beam. The third is the sum of the power scattered out of the beam and the amount of power absorbed by the medium. A considerable amount of work has been done in the fields of physics (e.g. Chandrasekhar [10]), astronomy (e.g. Mihalas [51]), oceanography (e.g. Mobley [53]), and elsewhere to understand participating media. Little, however, has been S - V L ( P , S ) = i ( P , S ) (2.16) 20 done in computer graphics. See Rushmeier [63] for an overview of these. • 2.5.3 Object Models The construction of object models is what is normally referred to as "modelling" in computer graphics parlance. The principle role of object models is to define surfaces: places where radiance distributions change dramatically over an interval of path length so small that it may be considered infinitesimal. The surface of a properly-defined object subdivides space into an "inside" and an "outside". 2.5.3.1 Polygons Polygons, especially triangles, are by far the most common object modelling con-struct. Each /V-sided polygon is specified by a list of ./V vertices. Three-dimensional objects may be modeled with polyhedra, i.e., sets of polygons. Each polygon repre-sents a "facet" of the object. By convention, the vertex list is specified in counter-clockwise order when the facet is viewed from "above". Special care must be taken to ensure that the object that results from the union of all the facets has a consistent inside and an outside, Otherwise, most Tenderers produce a confusing result. 2.5.3.2 Parametric Surfaces Parametric surfaces are made up of "patches": sets of points {P} defined by an equation of the form P = P ( « , u ) (2.17) where 0 < u < 1 and 0 < v < 1. 21 Al though there have been at tempts made to render patches directly, such as Shantz and Lien [67] and Lane et a l . [43], it typical ly proves efficient to convert each patch to a set of polygons. One of the most widely-used algori thms for this is due to Lane et a l . [43]. G iven a patch, the algor i thm tests it for planari ty wi th in a specified tolerance. If it is wi th in the tolerance, the patch is replaced by a corresponding polygon. If the tolerance is exceeded, the patch is subdivided using well-known subdivision procedures and the algor i thm is recursively appl ied. 2.5.3.3 Implicit Surfaces A n impl ic i t surface is made up of the set of al l points {P} which satisfy / ( P ) = 0 (2.18) for a funct ion / . L ike parametr ic surfaces, impl ic i t surfaces are generally converted to polygons for the purposes of rendering. One of the most widely-used approaches is the "marching cubes" a lgor i thm, developed independently by W y v i l l , et a l . in 1986 [86] and by Lorensen and Cl ine in 1987 [48]. Th is a lgor i thm voxellizes object space and then uses the value of / at each of the eight vertices of each voxel to generate polygons that span the voxel. One of the shortcomings of marching cubes and simi lar algori thms is a sensi-t iv i ty to sampl ing problems: Structures wi th spat ial frequencies above the Nyquis t l imi t may not be polygonal ized correctly. In part icular, there is no guarantee of pre-serving the topology of the polygonalized object. Recent work, however, by Stander and Har t [73] shows promise of al leviat ing this problem. 22 2.5.3.4 Constructive Solid Geometry A notable extension of impl ic i t surfaces is "construct ive solid geometry" or C S G . It extends the definit ion given in Equat ion (2.18) to use / to classify any point in Eucl idean 3-space as being inside (/.< 0), on ( / = 0) or outside ( / > 0) an object. It is thereby possible to define composite objects in terms of set operations act ing on such pr imit ives. A comprehensive treatment of C S G is the book by M a n t y l a [49]. C S G represents a convenient way to model objects in some but not all ap-pl ication domains. In addi t ion, C S G objects are easily rendered wi th a raytracer (discussed in Section 2.5.7.3). Faster rendering can be achieved either by apply ing the marching cubes algor i thm to polygonalize C S G objects or by performing C S G on polygonalized pr imit ive objects, but just as wi th impl ic i t surfaces, sampl ing prob-lems can occur. 2.5.4 Illumination Models We wil l review several i l luminat ion models in use in computer graphics. Space does not permit us to consider some of the more elaborate ones (such as He, et a l . [33] or Oren and Nayar [57]), but as these are not yet in wide use, perhaps the reader wi l l forgive the omission. Our discussion wi l l center on reflection models, as most Tenderers in computer graphics pay l i t t le attent ion to transmission. 2.5.4.1 Phong Illumination Models Phong i l luminat ion models are generally of the form (see, for example, Foley, et a l . [18], equation [16.15]) Lr = kaLa + [kd{N-S) + ksFs{S,V))EN (2.19) 23 where ka is the "ambient" reflection coefficient, La is an ambient radiance uniformly distr ibuted over ^l^8, kd is the diffuse reflection coefficient, ks is the specular re-flection coefficient, and EN is the normal irradiance parameter iz ing a direct ional luminaire. There are two popular choices for Fs. The original one given in Phong [58] corresponds to ( ( R ; • V ) " ^ if (R j • V ) > 0 F f ( S , V ) = < (2.20) 0 otherwise where R ; = 2 N ( N • S) — S is the "ref lected" source direction and is the Phong specular reflection exponent. The other popular choice for Fs comes from B l inn [6]: F f ( S , V ) = ( N • H ) n ? (2.21) where H = (S + V ) / |S + V | is the unit vector halfway between S and V and nf is the B l inn specular reflection exponent. Since we are not considering emissivi ty or transmission in this section, let us see how these i l luminat ion models correspond to Equat ion (2.9). Let the patch be i l luminated by the combinat ion of a direct ional source of normal irradiance EN in direction S (i.e., Q{ = 6s, 4>i = ^ s ) on £2^ and a radiance La constant over Us ing Equat ion (2.15), we can model the result ing incident radiance as Li = La + EN6(cos6i - cos 0 S ) 6(<t>i - <h). (2.22) If we substi tute this into Equat ion (2.9), we get Lr = La f / r ( S , V ) ( N - S ) d w t - + £ ? J V / r ( S I V ) ( N - S ) (2.23) Jilt For the sake of simplicity, we have assumed that dA has an unobstructed view of 24 directional light source V Figure 2.3: Illumination Model Test Configuration Our goal is to make Equations (2.19) and (2.23) equivalent for all possible La and Ejy. This means that the non-ambient terms of both equations must be equal, so This will allow us to convert a Phong illumination model into its correspond-ing B R D F . For the purposes of comparison of the various models, let us adopt the (com-mon) test configuration whose geometry is shown in Figure 2.3. A single directional luminaire shines at a sphere. Viewed from the center of the sphere, the luminaire is located at angle $L from the viewing direction. Figure 2.4 shows a series of images generated with the test configuration with a Phong illumination model using Ff. In this series, ks varies between 0 and 1 and kd is taken to be 1 — ks. The angular distribution of the specular peak is 25 0.00 0.25 0.50 0.75 1.00 10° 45° 90° 13 20" 5° 45° 90° 135° 30° 45° 90° 135 40° 50° °l 45° 90° 135° 45° 91 § IQji • • ) • i « ill O i l a CI C l l C t i ( [ K [ .... CCIC.CM Figure 2.4: Sphere shaded with a Phong illumination model, using Blinn's , for an assortment of incident angles with respect to the viewer ($/,), specular coefficients (ks), and specular distribution half-angles (3). qualitatively characterized9 by the specular half-angle f3, defined by Ff = \ = cos"? /3 Hence, for a given 3, In 2 (2.25) (2.26) s In cos 3 The ambient terms of Equations (2.19) and (2.23) must also be equal 1 0 . If we equate these and notice that for 7 / -1 2TT J, Oft (N • V ) 7 dur 7 + 1 (2.27) 9 Only in the case of Phong illumination modelling with Ff does 0 have a direct and obvious geometrical interpretation, but, as we will see (and Blinn [6] points out), it is qualitatively useful in other cases. 1 0 We can look at this as a consistency constraint: The same B R D F we use to shade a directional luminaire must also shade an ambient light source. 26 we get ka = kdir + ks[ FS{S,V)dui (2.28) So consistency demands that we have only two degrees of freedom in selecting ka, kd, and ks. In what follows, we wi l l take ka to be dependent upon kd and ks. 2 .5 .4 .2 T o r r a n c e - S p a r r o w I l l u m i n a t i o n M o d e l s The other major class of i l luminat ion models was first proposed in Torrance and Sparrow [75] and applied to computer-generated imagery in B l inn [6]. We shal l , however, follow the development given in Cook [15]. Torrance-Sparrow i l luminat ion models can be formulated direct ly in terms of their B R D F : WG D ^ S - V > = , ( N . S ) ( N . V ) ( 2 ' 2 9 ) where F is the Fresnel coefficient, G is the geometrical at tenuat ion factor, and D is the facet slope distr ibut ion funct ion. The Fresnel coefficient for unpolarized light and zero ext inct ion ([15] ignores extinct ion) is (9-c)2 \ (c (g + c) - 1) 21 F = '— 1 + , ) ' r ^7 (2.30) 1 + ( (g c ) - l ) 2 (c (g-c) + 1)' 2 ( 5 + c ) 2 where c = ( V • H ) , g2 = n2 + c 2 — 1, and n is the index of refract ion. The geometrical at tenuat ion factor is G — min <f 1 2 ( N ' H ) ( N ' V ) 2 ( N • H ) ( N • S) \ G ~ M I N I 1 ' (v^i) ' ( V T H ) J ( 2 ' 3 1 ) There are several choices for the facet slope dist r ibut ion funct ion. B l i nn [6] suggests three of them. The first corresponds to a Phong i l luminat ion model: D i = bi cos C l a (2.32) 27 where c o s a = ( N • H ) . The second is the Gauss ian one original ly used in Torrance and Sparrow [75]: D2 = 6 2 e - ( C 2 ° » 2 (2.33) The th i rd is f rom [76]: D3 (2.34) cos 2 a(c§ — 1) + 1 In al l of these, the 6's are arbi t rary constants analogous to the fc's in Phong i l lumi-nat ion models. The c's (empirically) determine the width of the specular lobe. A s B l inn [6] observes, if we define 3 to be the value of a at which a d ist r ibut ion drops to half its peak value, we have In 2 ci = C2 = In cos 3 cos 2 3 — 1 c 3 = C O S " 3-y/2 Cook [15] considers an addi t ional possibi l i ty or iginat ing wi th Beckmann and Spizzichino [4], which we wi l l include here as 1 _ | l - c o s 2 a \ DA = 5 -r-e W " " ' W (2.35) 4mz cos 4 a where m is the R M S slope of the surface. Unl ike D\ - D3, there is no arbi t rary b constant for this d is t r ibut ion. The relat ionship of m to the corresponding value of 3 is m = (2.36) V l n 2 - 4 1 n cos3 v ' 2 .5 .4 .3 N e u m a n n - N e u m a n n I l l u m i n a t i o n M o d e l s In [55], Neumann and Neumann discuss "separable" i l luminat ion models (i.e., those whose B R D F is of the form a(S) r ( V ) for some functions a and r) and how their use 28 can speed up radiosity computat ion in non-diffuse environments. A s an example, they describe a "lacquer model " of a purely diffuse mater ial covered by a semi-transparent " lacquer" that absorbs but does not scatter l ight that passes through it. The result ing B R D F they derive is / r ( S , V ) = oexp{- 5(^ + p ?i v y)} (2.37) where c and s are constants that characterize the model. We can make this comparable wi th Equat ions (2.24) and (2.29) by defining 6/v as the value of fr at S = V = N . The equation then becomes / , ( S , V ) = * W H [ p { - . ( i j i l 1 . + _ l _ _ 2 ) } (2.38) A s we did before, we can relate s to a more geometrical ly meaningful quant i ty 3 that qual i tat ively measures the width of the specular peak. Keep the i l luminat ion normal (S = N ) and increase the angle between V and N unti l fr drops to half of its max imum (i.e. V = N ) value. We define the result ing angle to be 3. We can relate s to 3: s = - - ^ - n (2.39) 1 + cos 3 K ' Figure 2.5 shows what Neumann-Neumann i l luminat ion models (and the subsequently-discussed Minnaer t i l luminat ion models) look like when applied to a sphere. Incident light for each i l luminat ion model is scaled to produce a peak unsaturated radiance at normal incidence (<£/, = 0) and then held constant as and 3 are varied. Note that as increases, the image radiance decreases, unlike Phong i l lu-minat ion models. A l so notice that the l imb of the sphere ( ( N • V ) = 0) is always dark. 29 p 10° 20° 30° 40° illumination * model l 0 ° 30 ° 60 ° 0 ° 30 ° 60 ° 0" 30 ° 60 ° 0 ° 30 ° 60 ° 50° Neumann-Neumann Minnaert ©•IODIODIODDOC © • I O D I O D & O D D O C Figure 2.5: Sphere shaded with Neumann-Neumann and Minnaert illumination models for an assortment of incident angles with respect to the viewer (4>L) and specular distribution half-angles (3). illumination ^ model L Phong/Blinn 10" 20' 30" 20° 10° 20" 30 I II ir 1 m n Figure 2.6: Comparison of a Square Shaded with a Neumann-Neumann Illumination Model and with a Phong Shader Using Blinn's . Neumann-Neumann illumination models exhibit peculiar behaviour when be-ing applied to a specular surface. As Figure 2.5 shows, they do produce an acceptable specular peak. Nevertheless, for a given incident angle the resulting radiance always peaks when the viewer is looking in the normal direction. For a highly specular surface, we should expect the radiance to peak somewhat closer to the reflected direction. Instead, the peak occurs where the exponent of Equation (2.38) is mini-mized, which may not be the same direction. Figure 2.6 illustrates this, comparing a Neumann-Neumann model with a specular Phong model using Blinn's F s s . To make sure the incident angle is constant, 30 we have replaced the sphere used in the test configuration of Figure 2.3 with a square whose normal is pointing at the observer. As the figure shows, increasing causes the highlight to move to the right for the Phong model, as we would expect of a specular surface, but not for the Neumann-Neumann one. This is not to say that Neumann-Neumann models do not apply to real-world phenomena, only that if they do, those phenomena do not resemble specularity. 2.5.4.4 Minnaert Illumination Models \ In [52], Minnaert describes an illumination model derived from observations of the Moon. His original model is Lr = bM(N-S)k(N -Vf-1^ (2.40) for some constants 6^ and k. This corresponds to the B R D F / r (S , V ) = 6 M ( ( N • S) (N • V ) ) * - 1 (2.41) We can relate k to an angle 8 defined as in the previous section: k = 1 - (2.42) In cos p-Figure 2.5 contrasts Minnaert illumination models with Neumann-Neumann illumination models. It is difficult to tell them apart. Their numerical values in these images differ by no more than 2%. (As a computational aside, since Minnaert illumination models are also separable, this suggests that a Minnaert illumination model should be able to take the place of a Neumann-Neumann illumination model with fewer arithmetic operations in most cases, especially if k is an integer.) Minnaert illumination models exhibit the same peculiar behaviour as Neumann-Neumann illumination models when applied to a specular surface: for a fixed in-31 cident angle, the resulting radiance peaks when the viewer looks in the normal direction. 2.5.5 Surface Mode l s The standard presentation of B R D F given in Section 2.4.3 leaves out an important fact: both of them may vary spatially, hence we should allow for 1 1 : We could deal with this somewhat by replacing the object being rendered with a collection of smaller objects with spatially-invariant surface properties, but this would cause the number of objects to increase considerably. Having to do this for a mottled surface such as a brick or an orange peel would inflate the number of objects to be rendered considerably and undesirably. One way conventional Tenderers deal with spatial variation of the B R D F is by texture mapping, factoring it into spatial and directional components: where the function t represents the surface texture: the spatial variation of reflective properties on a scale smaller than that of the object. u(P) and v(P) map spatial coordinates to texture coordinates. It is convenient to think of this as "wrapping" the surface with paper printed with a pattern that represents the texture. While this is a considerable reduction in flexibility from Equation 2.43, it is often adequate in practice. "Even this is an oversimplification, as it assumes that light encountering the surface at point P is reflected from the surface at the same point. A more thorough treatment takes into account subsurface scattering and is treated by Hanrahan and Krueger [31] and Lalonde [41]. / r ->/ r (P,S,V). (2.43) fr^t(u(F),v(P))fr(S,V) (2.44) 32 Another way Tenderers represent spatial BRDF as well as geometric variation is by allowing the surface normal used to compute fr to vary with texture coordi-nates -u(P) and u(P)..This was first devised by Blinn [7] and is referred to as bump mapping. As the name suggests, bump mapping is particularly effective at repre-senting surfaces whose reflective properties are the result of geometric variations such as roughness or machining as opposed to variations in material composition. Work continues in this vein. Poulin and Fournier [60] extend BRDF mod-elling to include anisotropic reflection phenomena, and Fournier [21] describes how to model a given BRDF with a "normal distribution function". 2.5.6 Shading Models There is a tendency to use the term "shading model" synonymously with "illumi-nation model", but as Foley et al. [18] point out, a shading model may or may not make use of an illumination model. An illumination model is a description of what is actually happening in modelling or "world space". A shading model is a description of what is to be shown in viewing or "screen space". A shading model is what actu-ally generates pixels. It is possible for a shading model to use an illumination model to generate a pixel, but it could also use a screen-oriented interpolation scheme to do so. For instance, the Goraud shading model typically uses an illumination model to compute radiances at polygonal vertices, but then interpolates these in screen space during scan conversion. On the other hand, the Phong shading model (not to be confused with the Phong illumination model) computes normals at vertices, interpolates the normals in screen space during scan conversion, and (typically) applies an illumination model at each pixel. 33 An essential part of a shading model is the number of samples of an illumi-nation model required per pixel. If this number is greater than one, as in the case of jittered or distribution raytracing, the samples must be filtered. If this number is less than one, the samples must be interpolated. The use of a shading model is indicative of a perceived difference between what is necessary to model physically what is going on in the scene and what is nec-essary to generate an acceptable image. We will return to this idea in our discussion of radiosity in Section 2.5.7.4. 2.5.7 Rendering Techniques The most distinctive part of a Tenderer is the rendering technique it uses. In this section, we will examine the predominant rendering techniques, ray tracing and radiosity, along with a new and promising one, photon maps. 2.5.7.1 The Rendering Equation In his landmark work [38], Kajiya showed the unity of rendering techniques by connecting them to his form 1 2 of Equation (2.11). We can do much the same thing in the current context by noting that in Equation (2.12), in a state of thermodynamic equilibrium, Li should be equal to TBLr. In other words, if we want a truly global solution, there should only be a single radiance distribution L which should satisfy L = Le + STBL (2.45) which we rewrite as L = QLe (2.46) 1 2 Like several other authors, Kajiya formulated Equation (2.11) in terms of "point-to-point" transfers, while we have chosen a "direction-to-pbint-to-direction" approach. The two are entirely equivalent. 34 where we have defined the global illumination operator Q to be G = ( I - STB) - l (2.47) where 1 is the identity operator. We can expand G as Kajiya does: G = 1 + STB+(STB)2 + •••. (2.48) He then points out an interesting physical interpretation of this series. The first term represents light that is being given off by the surface on which we are computing L itself, the second term represents light from sources that directly illu-minate the surface (a single "bounce" or surface interaction), and so on, so that the nth term represents light after n — 1 "bounces". He also relates convergence of the series to energy conservation13. As Kajiya points out, rendering techniques approximate G in different ways. They also simplify the integral in Equation (2.11) (i.e., in S) differently. We will discuss these distinctions in the following sections. Over time, these techniques have undergone hybridization 1 4 in order to improve their efficiency or functionality as they get incorporated into real Tenderers, but as we are more interested in their theoretical distinctions here, we will deal with each rendering technique in its "classical" form. Due to space constraints, we will confine our discussion to currently-popular rendering techniques, omitting such work as that of Moravec [54], Weiler-Atherton [81], Warnock [79], and Roberts [62], which are now mainly of historical interest. 1 3 Mathematically, this is equivalent to requiring that the spectral radius of the operator STB be less than unity. 1 4 Tha t ' s the nice name for it. 35 2.5.7.2 Depth-Buffering By far the simplest rendering technique is depth buffering, also known as "z-buffering" First developed by Catmull [9] in 1974, the algorithm requires a buffer that contains a scalar value for each pixel. This value is the "depth" - the distance from the viewpoint to the nearest object that the pixel represents. Initially, all values in the depth buffer are (conceptually) infinite. While primitives are being scan converted, the algorithm computes some measure of the distance from the object to the view plane for each pixel. In the case of opaque objects, this distance is compared with the depth buffer value. If (and only if) the new distance is smaller than the previous distance, the colour buffer is set .to the appropriate value for the new primitive and the new distance is stored in the depth buffer. With some limitations, non-refractive transparency is possible with a depth buffer. Williams [84] presents a way to use depth buffers to generate (hard) shadows. Since it places no constraints on the method used to compute the values that go into the colour buffer, depth buffering is compatible with all local illumination models. Depth buffering can be implemented in hardware or software. The OpenGL graphics library ([40] and [8]) includes built-in support for depth-buffered rendering. Hardware depth buffering provides very fast rendering. OpenGL makes use of it if it is present and simulates it in software if it is not. 2.5.7.3 Ray tracing Classical ray tracing simplifies Equation (2.11) by setting fr and ft to be Dirac 5-functions (which we will describe below). This permits (mirror-like) reflection and refraction. The second way is by setting L{ to include a ^-function, which, as we will 36 see in Sections 2.5.1.1 and 2.5.1.2, permits point and directional luminaires. These are mutually exclusive since we cannot allow both the incident radiance and the B R D F or B T D F to be ^-functions. Since each surface interaction involves samples of the B R D F and/or B T D F and possibly casting additional reflection/refraction rays, ray tracing is compatible with all illumination models. 2.5.7.4 Radiosity Although radiosity methods had been present in the mechanical engineering commu-nity for some time, they were first introduced to the computer graphics community by Goral, et al. [25] in 1984. In its classic form radiosity takes L to be isotropic and constant over a patch and therefore equal to ^ , where B is the radiosity of the patch. As derived in, for example, Cohen and Wallace [14], this leads to the classic radiosity equation: (I - pF) B = E (2.49) where, given N patches, pF is an N X N matrix whose elements are pjFj,-, p; is the (isotropic) reflectivity of patch i, Fji is the purely geometric form factor: the amount of energy leaving patch j that arrives at patch i, B is an A^-element column vector of radiosities, and E is an A^-element column vector of patch emissive exitances. Several conventional numerical linear algebra schemes have been applied to solve Equation (2.49), including Gauss-Seidel [25], Southwell (also known as "progressive refinement") [26], and multigrid [46]. Classical radiosity is compatible only with diffuse illumination models, but more recent research has extended this to glossy reflection. Immel et al. [35] devel-oped a "global cube" algorithm that collected radiances instead of radiosity. This 37 i approach produced realistic images, but was costly in both time and memory, es-pecially for highly specular surfaces. Rushmeier and Torrance [64] developed the notion of "extended form factors" to deal with specular (including mirror) surfaces. Sillion et al. [69] deal with glossy reflection by representing both B R D F s and radi-ances with spherical harmonics. Even more recently, Fournier [19] has shown how to incorporate non-diffuse phenomena by the use of separable B R D F s . An important consideration in viewing the results of radiosity and other global illumination computations is the role of what Christensen, et al. [11] (citing Reichert [61]) refer to as the "final gather" step. It takes place after Equation 2.49 or its equivalent has computed radiance distributions (or their equivalent) on all visible surfaces and amounts to a reevalu-ation of Equation 2.11 at the point in the scene intersected by a ray through the middle of each pixel. As Christensen, et al. [11] put it, Formally, this final gather corresponds to changing to a piecewise-constant basis, where the support of each basis function is the projection of a pixel onto a surface in the scene. This basis is tailored to be visually pleasing. The final gather smoothes the discontinuities in the wavelet representa-tion and makes highlights, textures, and shadows crisper. A final gather is effectively a local illumination calculation performed at the end of a global computation. It is an implicit acknowledgement of the fact that the criteria for physically-correct global illumination and a "visually pleasing" image are not identical. 38 2.5.7.5 Photon M a p s One of the most promising of recent rendering techniques are photon maps, devel-oped by. Jensen [37]. A photon map is a four-dimensional distribution of incident photons collected on a surface. The first step of photon map global illumination is the generation of two photon maps built up by casting photons from all luminaires. The first photon map is a high-resolution caustics photon map made up of photons cast only in the direction of specular objects. The second photon map is made up of photons cast towards all objects. The second step is the rendering step. In it, the photon maps of incident photons are combined with arbitrary B R D F s to produce reflected radiances made up of four terms: direct illumination, specular reflection, caustics, and soft indirect illumination. Each term is computed with appropriate photon maps and B R D F specular and/or diffuse components. As Jensen puts it: The photon map is used to generate optimized sampling directions, to reduce the number of shadow rays, to render caustics, and to limit the number of reflections traced as the scene is rendered with distribution ray tracing. Lucifer, the algorithm we will present in Chapter 3, is coeval with photon maps. We will see that the two have much in common. 2.5.7.6 Comparison and Summary This section will conclude with an analysis of the strengths and weaknesses of each of these models. 39 Chronological Order of Development —> Depth Ray Photon Effect Buffer Tracing Radiosity Maps diffuse surfaces yes yes yes yes specular highlights yes yes no yes transparency yes yes no yes mirror reflection no yes no yes refraction no yes no yes sharp shadows yes yes no yes soft shadows no no yes yes diffuse lighting yes no yes yes colour bleeding no no yes yes caustics no no no yes Table 2.1: Effects Available with Established Rendering Techniques Table 2.1 lists a number of illumination effects and whether or not the vari-ous established rendering techniques are capable of them. One word of caution: If we evaluate this table for specific Tenderers, it becomes clear that it is an oversim-plification. Some Tenderers which refer to themselves as raytracers, for instance, are capable of soft shadows. This is because those Tenderers have incorporated radiosity or other methods in hybrid form and are not "pure" raytracers in that sense. For this reason, Table 2.1 should not be taken to refer to specific Tenderers, only to rendering techniques. Indeed, ray tracing and radiosity techniques can be enhanced (often by Monte Carlo methods) to produce a wider range of effects than Table 2.1 indicates. As in many other areas in which they are applied, however, Monte Carlo techniques exhibit slow convergence. To be able to produce these and other effects, we will elaborate on a global illumination technique that allows a wider range of local illumination models. 40 Chapter 3 The Lucifer Algorithm The algorithm we present here, which we refer to as the Lucifer1 algorithm, is a sequential computational analogue of how energy distributes itself in a scene. It combines the physical concepts of isolation and energy conservation with a refineable spatial partitioning scheme. This work is a outgrowth of work done by Fournier, Fiume, Ouellette, and Chee [20] (hereafter referred to by its project name "FIAT"). 3.1 Isolation For possibly complex scenes, we can take a divide-and-conquer approach based on the principle of isolation, illustrated in Figure 3.1 (right). Isolation is a conceptual tool: we place any part of our scene (including the viewer) within a volume V. Suppose then that somehow we can determine the distribution of radiance on the surface dV of V . Isolation says that for the purposes of solving for global illumination anywhere outside V, we can effectively replace the 'From the Latin for "light bearer". 41 directional Unlit source o Figure 3.1: An Application of Isolation contents of V with dV's radiance distribution. Isolation gives us a way to deal with complex scenes. If we cannot deal with the whole scene, break it into volumes we can deal with as local illumination prob-lems and transfer radiance distributions along the boundaries between the volumes. 3.2 Power Computation for Isolated Volumes For any volume V,-, physics demands energy conservation. In the steady state, this is equivalent to power conservation: + $ e m , t = *o»M + *»b,, i (3.1) where $ , n , = / I L i n \N-S\dudA (3.2) is the flux entering V,-, $ o « t , , - = / / ^ c u t | N • S| dudA (3.3) • 12 is the flux leaving V j , * . m , , - = / / jdudVi-r [ [ Le \N-S\dudA (3.4) is the flux generated within V ; , f i N and fi^, at any point on dV are the unit hemispheres entirely outside and entirely inside V , respectively, L i n is the radiance coming into V ; from fiN, LoM is the radiance passing out of V j from Q^, Q is the unit directional sphere, j is the volume emissivity within V , Le is the surface emissivity on all emissive surfaces 5 e m ( V t ) within V , and <J>ab3,i is the amount of radiant power absorbed and not converted into radiant energy again (at least in Vj ) . 3.3 Solution Order and Convergence As with raytracing and radiosity, we are constructing a sequential analogue of what nature does in parallel. We need a sequential ordering. A reasonable way to proceed is to maintain the volumes in a queue sorted in order of decreasing undistributed power $ i n j + #em)t- and to always concentrate our efforts on the volume V j at the front of the queue. Once the undistributed power in this volume is distributed, the volume moves to the end of the queue. Note that, from Equations (3.3) and (3.4), undistributed power is computed from radiances, emissivities, and geometrical information about V j . If, during the distribution process, we ensure that no more energy leaves V j than was either incident upon or emitted within it, we can guarantee that the total undistributed power in the scene (i.e., in the queue) is monotonically non-increasing, since all the power distributed by V j that does not get absorbed becomes incident on another volume or leaves the scene entirely, possibly reaching the observer. 43 One way to ensure that no excess energy is created is to require an energy-conserving local illumination model. We will present some necessary constraints for such a model in Section 4, where we will consider several illumination models commonly in use in computer graphics from the aspects of both energy conservation and Helmholtz reciprocity2. 3.4 Spatial Partitioning We require a tessellation of the scene that can be easily refined as needed as the solution progresses. Several data structures permit this, but the simplest one for our purposes is the octree. As shown in Figure 3.2, we take octrees to be composed of cubic cells. The root cell contains the entire scene. Each cell is either a parent cell or a leaf cell. A parent cell has eight child cells. A leaf cell has no child cells. Only leaf cells are "volumes" in the sense of our previous discussion. Non-leaf cells are purely structural. Each cell except the root cell has seven sibling cells. We also choose a minimum size for an octree cell. Cells that we cannot deal with that are this size will be solved trivially, but still conserving energy. 3.5 The Lucifer Algorithm Figure 3.3 shows our rendering algorithm pseudocode. We start off with a scene s with a single cubic cell, rootCell(s). We then create a priority queue q with, initially, a single element, rootCell(s). q is always maintained in decreasing order of undistributed power. 2Note that Lucifer as presented here does not require reciprocity. 44 parent cell Figure 3.2: Octree Cell Nomenclature If q is empty or sumUndistPower(q), the sum of the undistributed powers of all cells in the queue, is lower than some prespecified powerTolerance, we consider the scene rendered. Otherwise, we remove the cell with the largest amount of undistributed power, c, from the front of q. All further computation in the main loop is con-cerned only with c. This is why we refer to our algorithm as "light-driven": it is always concerned with the cell with the largest undistributed power. We test c against the first of these cases that it matches: • If it is empty, we call propagatef) to transfer the radiance coming into c to its neighbors. • If it is a cell that we know how to pass power through, we call balance() 45 render(scene s) { queue q; cell c; q = createPriorityQueue(rootCell(s)); while (sumUndistPower(q) > powerTolerance) { c = dequeue(q); if (isEmptyCell(c)) propagate(c); else if (canDealWith(c)) balance(c); else if (size(c) > minimumSize) subdivide(c, q)\ else trivialize(c); } } Figure 3.3: The Lucifer Algorithm to perform the redistribution of the reflected and refracted radiance to c's neighbors. • If it is above a specified minimum size, we call subdivide() to split c into eight child cells and add these back to q in order. If it meets none of these criteria, we call trivialize() to perform an ad hoc solu-tion that may be some combination of the first two cases, but which in any case guarantees that power balance is preserved. propagate() is conceptually straightforward. It makes use of the invariance of radiance in the direction of propagation. Radiance at a given point in a given direction on an incoming cell wall defines a ray which may intersect an outgoing cell wall at a new point and direction. (The direction is, of course, the same, but may be expressed in a different coordinate frame.) We may alternatively project backwards from the outgoing cell wall to the incoming one. 4 6 balance() is somewhat more complicated. The ambiguity of the phrase "know how to pass power through" at this level of description is intentional. Our current interpretation of it is that the cell contains at most a single convex object. This allows us to decompose the procedure into four subtasks: 1. Transport radiance from the incoming walls to the object (or that part of the object which lies within the cell). 2. Interact the incoming radiance with the object's B R D F (and B T D F ) to pro-duce a reflected radiance being given off by the object. 3. Transport reflected and emitted radiance from the surface of the object to the outgoing walls. 4. Transport the unobstructed part of the incoming radiance to the outgoing walls. Subtasks 1 and 3 are variations on propagatef), allowing transport to and from surfaces as well as cell walls. Subtask 2 requires surface interaction, as dictated by Equation (2.11). Subtask 4 is another variation on propagate(), including a geometric blocking function that inhibits propagation if the ray intersects the object within the cell. We call this variation transblockf). 3.6 Complexity Let us compare the complexity of Lucifer and two rendering techniques that can produce comparable effects: raytracing and radiosity. Our measure of complexity is the mean time required per displayed pixel. We assume the scene has Nob-t objects, of which i V l u m are (initially) luminaires. 47 In classical raytracing, computing the intersection point of a ray and an object in the scene requires 0(Noh)) time per ray 3 If, to reduce aliasing effects, we oversample by a factor of iV r a y rays per pixel, computing the primary (starting at the observer) intersections is 0(NI!LyNobi). To determine shadows, we must cast N[um shadow rays to light sources. As Whitted [83] first showed, we can treat reflective, refractive, and transmissive phenomena by constructing a tree at each intersection whose links correspond to reflected, refracted, and transmitted rays and whose nodes represent ray intersections with the scene. If we let D be the mean number of nodes in this tree (minimum of one), then noting that each ray, shadow or reflected/refracted, costs A o b j , the efficiency of classical raytracing is 0(NI!LyNobi(D + A r l u m ) ) . As summarized by Arvo and Kirk in [1], many techniques have been devised to improve on this. With suitable restrictions on such Tenderer components as illumination model, luminaire model, and the spatial distribution of objects, considerable speedups are possible. Fujimoto et al. [23] cite rendering times that are practically independent of the number of objects in the scene (i.e.,'0(Arra.yZ))) but this is quite optimistic. In classical radiosity, the computation of form factors is 0(N2bi) and while the inversion of the transfer matrix is in principle as high as 0(N03bl), it easily can be made 0(N2bi) by iterative and other means. Considerable speedups of radiosity are possible with hierarchical radiosity, as presented by Hanrahan et al. in [32]. Smits et al. in [71] extended this work by introducing clustering. As reported there, hierarchical radiosity has a complexity of 0(N2uif + Np::,tl,h), where Nsuit is the number of initial surfaces (proportional to Nobi) and NpiLich is the final number of patches (which is at least a linear function of Nsur{). Clustered radiosity, on the other hand, 3This is assuming a limiting depth - see Fiume [17] for additional details. 48 has a complexity of 0(NSUIf log N3Ulf + iVpatch). The linear or superlinear dependence of ^ p a t c h ° n -^surf .has given rise to a considerable amount of interest in optimal meshing schemes. Obviously, a clever adversary can defeat these speed-up schemes, but in practice the expected values apply. In Lucifer, interactions are never between objects, but between a cell's walls and either the walls of neighboring cells or the object contained within the cell. Even though a particular cell may need to be balanced more than once, the number of times that needs to happen is dependent on the scene geometry, not on 7V c e l l. If the proportion / b a c k of power $ that comes back to the cell is fixed, then in order for the power to drop by a fractional tolerance ep w r, the cell must be balanced 7V b a l times, where f&?* = ^*- (3-5) Taking the log of both sides, we find that ATbal is proportional to the log of ep w r. Hence, as may be inferred from Figure 3.3 the asymptotic efficiency of Lucifer is 0(T c e l l7V c e l l) where T c e l l is the cost of operation within a cell (a function of the maximum resolution used to express the radiance). Since our spatial partitioning scheme guarantees the number of objects, Nceli is 0(iVobj), Lucifer is intrinsically 0(Nobi) - asymptotically more efficient than other global illumination schemes. It is important to note that in these techniques visibility is an important part (often the most important part) of the cost, and what makes it superlinear. Lucifer in a real sense "clusters" visibility, since it is carried along with the radiance representation from cell to cell. As an example, consider a light source blocked by a large object. Before it reaches the blocker, the light flux is represented normally, and small shadows will be carried along as well. After it reaches the blocker, a large portion 49 of the light flux will be zero, the small shadows included in it will be absorbed. Any scheme capable of compressing the light flux data will take advantage of this, and objects and cells in the shadow will not even interact with the resulting light front. To illustrate the trade-off consider the following scene. We have 7Vobj ~ 107 objects4. Any approach that has to consider pairs (either for explicit light transfer or for visibility) will have to deal with an order of 10 1 4 such interactions. To represent light fluxes at an acceptable level of detail, assume that the space is 10m by 10m by 3m, and we want about 1 cm accuracy. If we impose a maximum cell size of l m 3 , that means we have have 300 cells to balance. For each cell wall, we then require a positional resolution of 100 by 100 elements. If we assume the same resolution for direction (l/100th of a radian is approximately 0.57°) , we need to represent 108 elements on each wall. If we can use wavelet compression by 1000:1 (and we will see in Section 6.13 that this is not unreasonable), radiances on the walls can be represented with 105 coefficients. Even assuming that at each step of the propagation we deal with 0(N2) interactions between elements (wall-to-object, wall-to-wall, and object-to-wall), that means ~ 10 1 0 such interactions. Multiplying by the number of cells, we get ~ 3x 10 1 2 interactions, which is a small fraction of the pairwise number, and this is not taking into account the speedups possible due to the hierarchical nature of the wavelet representation: If two coefficients do not interact, neither do any of their children. The actual number of such computations is a function of how many times and at what level of detail we have to do this, but even at this level we have some hope of winning. An important point here is that the Lucifer numbers do not significantly 4 We should note here that in the Lucifer context, the definition of the term "object" might be more powerful than in some other contexts. 50 Figure 3.4: "Balls" Model increase as a function of the number of primitives in the scene. 3.7 Light Through a Window Let us consider an example of what radiance L on a cell wall looks like. L is a function of four variables (2 positional, 2 directional). Figure 3.4 shows an instance of the "balls" model from Haines's Standard Procedural Database of test models for raytracers (described in [29]) rendered by a typical raytracing T e n d e r e r . We have adapted this Tenderer to perform 4-dimensional ray tracing after a "light through a window" model. In this, we imagine light from a scene going through a unit square window. We further imagine sampling the radiance over the incoming hemisphere of directions on a uniformly spaced grid of positions. (In Lucifer, the "window" corresponds to a cell wall, but we'll discuss the 51 1.0 — Figure 3.5: Array of Fixed Direction Views more general case here.) Before proceeding further, let us transform our definition of direction from a (6, (f>) parameterization to (fix,fj,y), where \ix = sin 9 cos (f> ny = sin 9 sin cf> (3-6) This avoids discrepancies close to 9 = 0 and makes for a more uniform grid in Cartesian 3-space. \j?x + fi2 = 1 defines the unit directional circle (hereafter, UDC) . All real directions must lie within this circle. The modified Tenderer permits us to create an array of images, each image corresponding to a single direction (Figure 3.5) or a single position on the window 52 0.75 — 0.5 0.25 — •••$9 X 0.0 0.25 0.5 0.75 1.0 Figure 3.6: Array of Fixed Position Views (Figure 3.6). In the former case, the individual images are parallel projections and in the latter case they are "fisheye" views. Note that, as in a perspective view, except in the special case of an orthographic (parallel) projection, a parallel projection of a sphere is not round. From these examples, we can see that a 4-dimensional representation of ra-diance exhibits the same combination of discontinuities and relatively smooth areas we find in 2-dimensional images. r>3 3.8 Radiance Representation Radiance on a cell wall is a potentially discontinuous, generally non-analytic function of four variables (2 positional, 2 directional). We can represent radiance L at a point P on a cell wall or a surface and in a direction (JJLx, /j,y) with a finite element expansion with Nj degrees of freedom: Even though the basis functions B{ may take on values for fj,x + /J,2 > 1, we can disregard their behaviour there, since we never evaluate them in that region. Choices for B{ include: box discretization, Fourier, discrete cosine, orthog-onal polynomials, "light fields" (as described by Levoy and Hanrahan in [44]) and wavelets. 3.8.1 B o x Discre t iza t ion For box discretization, as done by FIAT, the P>i are constant within a quantized direction for each quantized position. Figure 3.7 shows the result of Lucifer using a box discretization on a simple model consisting of three squares forming a corner of a (root) cell with two spheres suspended within the cell. The square on the upper right is a diffuse white emitter and the other two squares are diffuse red (bottom) and gray (upper left) reflectors. One sphere is a mirror and the other is a diffuse green reflector. This image illustrates some of the effects possible with Lucifer - diffuse surfaces, colour bleeding, soft shadows, and mirror reflection - all produced with the same rendering technique. There are various ways to perform propagation and object interaction of the discretized radiances. FIAT describes one possible way, equivalent to the one (3.7) 54 j Figure 3.7: Box Discretization Example we used for Figure 3.7. They have a common shortcoming, however: by quantizing angles and positions, box discretization causes various rendering artifacts. These can be reduced by increasing the number of quantized positions or directions, but this causes memory requirements to increase dramatically, even with dynamic allocation of memory. In practice, the CPU time and memory required for box discretization is limiting. Figure 3.7, for example, required about 9.3MB of memory and 3.2 CPU-hours on an SGI Indigo 2 workstation. 3.8.2 Wavelets Wavelets are a promising alternative to other basis functions. They were successfully applied to radiosity solutions by Gortler, et al. [28] and Schroder, et al. [65]. Schroder and Hanrahan [66] and Christensen, et al. [12] extended this work from 55 Figure 3.8: Original Image (left) and Its Reconstruction from 4-Dimensional Wavelet Coefficients (right) radiosity to the representation of radiance itself. In this section, we will consider the general applicability of such a wavelet representation of radiance to Lucifer. Appendix A describes some of the properties of one-dimensional wavelets. There are two properties of particular interest for radiance representation. First is their ability to approximate L2 functions, even those with discontinuities, with a relatively sparse set of coefficients. As an example of this, we apply wavelet compression to the 4-dimensional data shown in Figure 3.5. On the left of Figure 3.8 is a magnified view of one image from Figure 3.5 in a particular direction, on the right is a view in the same direction that was reconstructed from a 4-dimensional wavelet transform of the original data compressed by 96% - only the top 4% (in magnitude) of the original coefficients were retained. It is important to note that this compression was in all four dimensions 56 of the data prior to reconstruction. There is an infinite number of base scaling functions <j>(x) and therefore an infinite number of possible wavelets. For the example above, we used the L = 2 "Coiflet" wavelet described in [16], for these, but other wavelets gave qualitatively similar results. Part of our work will be to find which wavelet basis works best for radiance representation. Most importantly, we need to evaluate the tradeoff between improving approximation by increasing Nv, the number of vanishing moments, and minimizing operation count by reducing Wh, the size ("support") of {hm}, the set of wavelet coefficients. Other considerations in selecting a wavelet include: • Is <j>(x) symmetric about some value of x? • Does 4>(x) need to have an analytical form? (This leads to biorthogonal wavelets.) • Does refinement of the wavelet interpolate the coarser values? • What is the tradeoff between reducing Wh and representing discontinuities compactly? The second property of wavelets relevant to Lucifer is their dyadic nature, which corresponds well with our octree spatial subdivision scheme. Splitting a cell whose walls represent radiance with with wavelets amounts to a partitioning of the wavelet coefficients. We will consider these issues further in Chapter 5. 57 Chapter 4 Making Illumination Models More Physically Plausible In our presentation of the Lucifer algorithm in Chapter 3, particularly in Section 3.3, we described the need for illumination models that conserved energy. We will de-scribe such models in this chapter. This work has been published separately in reference [45] and (more accessibly) in reference [47]. Illumination modelling computation is an essential part of any rendering al-gorithm. Getting an exact physical model of the interaction of light with a surface is, for most surfaces occurring in the real world, a very difficult problem. Conse-quently, much effort has been expended on finding approximations that are both good-looking and quickly computed. An extensive, if somewhat dated, summary of these illumination models is in Hall [30]. Looking good and being quickly computable are sufficient criteria for most raytracing Tenderers (see Section 2.5.7.3). If one of these illumination models is used in a radiosity computation, however, it is necessary to additionally ensure that, since 58 radiosity is based on principles of energy conservation (see Kajiya [38], for example), the illumination model itself conserves energy. The global illumination model we present in Chapter 3 also requires energy conservation. Another recent technique described in Larson (ne' Ward) [78] constructs functions to fit actual BRDFs. Such fits are likely to be more successful if the func-tions themselves are, like the data they fit (as demonstrated empirically in Clarke and Perry [13]), consistent with physics (at least within experimental uncertainty). So a need has arisen for illumination models that are not only good-looking and easy to compute, but are also "physically plausible". 1 In this chapter we will look at the correspondence between the fundamental reflectance equation (2.9) and traditional illumination models. Then we will look at ways to modify those models to make them more physically plausible. 4.1 Energy Conservation The first physical constraint we will examine with respect to illumination models is that of energy conservation. Physically plausible illumination models must obey energy conservation. In a steady-state scene, the input power and the rest of the scene configuration (objects and their positions, surface properties, etc.) are taken not to vary with time. In such a situation, energy conservation is synonymous with power conservation. The total amount of power reflected, i.e., M dA, where M is the exitance, must be less than or equal to the total power incident Ed A, where E 'We use the term plausible here in contrast to that of feasible in Neumann and Neumann [55]. A "feasible" illumination model is one that we can imagine constructing physically. This is not always possible. The weaker definition of a "plausible" illumination model is one whose existence does not violate physics. 59 is the irradiance. This must hold for all areas dA in our scene, hence M < E (4.1) From Siegel and Howell [68], we have an equation similar to Equation (2.7) describing the exitance dM due to a reflected radiance Lr radiated into an infinites-imal solid angle dur around a direction V , so that dM = (N • V ) Lr dujr (4.2) We substitute Equation (2.9) into this and integrate over Q N to get M = / \ / , / r ( S 1 V ) L i ( N . S ) ( N - V ) d w 1 - d w r (4.3) We can also integrate Equation (2.7) to get E= f Li (N • S) dui (4.4) So, if we make the trivial assumption that E > 0, we divide both sides of Equation (4.1) by E to get /«+ /n+ ^(S, V ) L{ (N • S) (N • V ) ^ <Ly / n + L , - ( N • S) dwi - 1 ( 4 - 5 ) Energy conservation does not depend upon the particular Li distribution. Given any Li, Equation (4.5) must hold, so, as we did with the Phong illumination models in Section 2.5.4.1, let us use a directional light source of the form Li = EN 5(cos0i - cosdt) <5((/>,• - <&) (4.6) to represent a directional source of normal irradiance EN in a direction S. Accord-ing to the ANSI/IES standard [36], M/E in this case becomes the "directional-hemispherical reflectance", which we will refer to as kp2. Integrating the ^-functions 2 In Neumann and Neumann [55] , this is referred to as "albedo", but that usage is imprecise as the definition of that term does not require a unidirectional source. In addition, "albedo" is not defined in the standard. 60 Figure 4.1: Specular Integrals HS(S) for Phong's F f (left) and Blinn's F? (right) p 10° 45" 90" 135" 20l 45" 90° 135° 30l 45° 90" 135° 40' 45" 9<>" 135° 50l 45° 91 0.00 • ! | n ) • ) m j I § 1 0.25 E 1 D • : f 0.50 • c 1 I • | • 0.75 t j f c 1 E 3 D NO j • 1.00 D o c 1 • E > 1 Figure 4.2: Sphere shaded with an energy-conserving Phong illumination model, using Blinn's F S B , for an assortment of incident angles with respect to the viewer ( $ L ) I specular fractions (ka), and specular distribution half-angles (/3). and cancelling out common factors, we get kp = I / r ( S , V ) ( N - V ) f l ! u ; r < 1 (4.7) 61 4.1.1 Making Phong Illumination Models Conserve Energy Let us apply these results to a Phong illumination model to see what constraint(s) energy conservation leads to. Recalling Equation (2.27), Equation (4.3) becomes M = EN[kdn(N-S) + ksHs{S)} (4.8) where we have defined HS(S)= f F 8 ( S , V ) ( N - V ) d u v (4.9) Jot, Figure 4.1.shows Hs evaluated numerically usi-rig both Equations (2.20) and (2.21) for a variety of specular exponents. Note that Hs is a function of the incident direction and specular exponent only and that it can be thought of as an integral operator acting on a given Fs. As we might expect, Equation (4.4) becomes E = EN(N-S) (4.10) so that k p = kdn + k s ^ t (4.11) To guarantee energy conservation regardless of reflection geometry, it is nec-essary to guarantee that kp < 1 for all incident directions. But there is a problem here. Given the F s 's in Equations (2.20) and (2.21) and regardless of S, it is always the case that Fs > 0 and, furthermore, there is always some non-vanishing region of £2+ over which Fs > 0. That means that Hs is always > 0, as Figure 4.1 illustrates. So that if ks > 0, it is always possible to choose close enough to 90° that kp will be greater than one. We therefore conclude that the specular terms of Phong illumination models do not conserve energy at sufficiently large incident angles. 62 After Neumann and Neumann [55], let us consider a different formulation of an illumination model. Start from Equation (2.19), but suppose that, instead of being constant, we allowed ks to vary with S in such a way that energy conservation was maintained. (As Equation (4.11) shows, we are not getting any trouble from the diffuse term, so we will leave it alone.) Let ka be the fraction of exitance that is reflected specularly: _ Jn+ dMspec K = M 1 + ^TT(N-S) ksHJS) We can solve Equations (4.11) and (4.12) for kd and ks to get (4.12) kd = K kp (1 - ka (4.13) kp K (N • S) HS(S) so we can rewrite the B R D F for the new illumination model as / R(S,V) = * P l - k c kaFs{S,V) (4.14) HS(S) We can also construct the analogue of Equation (2.19) to express this result in terms of radiance: Lr = kaLa + kp (N • S) i - K kaFs(s,v) 7T HS(S) Ed where h — h 1 +*-(4 ( N-s»^r , f c'-1) (4.15) (4.16) corresponds to Equation (2.28). Figure 4.2 shows what such an illumination model looks like when applied to a sphere with a single directional light source and no ambient radiance. For this 63 figure, we have used Blinn's F°. We have also taken kp = 1, since any other value would just be a uniform reduction by a constant factor in image radiance. Notice that, unlike Figures 2.4 and 2.5, the highly specular parts of the printed images are necessarily clamped in order to show the diffuse parts. 4.1.2 Do Torrance-Sparrow Illumination Models Conserve Energy? Figure 4.3 shows some numerical integrations of Equation (4.7), contrasting Phong illumination models with Torrance-Sparrow illumination models. All Torrance-Sparrow illumination models were computed with a Fresnel factor F = 1 (i. e. a large index of refraction) to show the worst case. One way to produce an energy-conserving Torrance-Sparrow illumination model suggests itself: simply choose any value of bj such that 'kT]' [ b3 \ max (4.17) is the maximum value as shown in Figure 4.3. (The Beckmann-max Spizzichino distribution is not a problem as long as its integral is always less than unity, and it has no 6-coefficient to adjust anyway.) Nevertheless, doing this would probably be a mistake. To see why, look at the plot for kTl/b\, the Torrance-Sparrow illumination model with the Phong microfacet distribution. Notice that it does not diverge as Oi 90°, even though k^/ks, the corresponding Phong illumination model with a Phong specular term, does diverge. The same is true for kj2/b2 compared to kB/ks. Why should this be? The answer lies in the geometrical attenuation factor G. As Oi —> 90°, G is guaranteed to be less than or equal to unity and, if ( V - H ) > 0 (i.e., V and S are not antiparallel), it will vanish in the limit. 64 where 0 1 ' ' 1 0 1 ' • 1 0 30 60 90 0 30 60 90 Figure 4.3: Directional-hemispherical reflectances as functions of incident angle for a Phong illumination model (kp), a Blinn illumination model (kf), and Torrance-Sparrow illumination models with Phong (kj1), original Torrance-Sparrow (kj2), Trowbridge (kj3), and Beckmann-Spizzichino (kj4), microfacet distributions. 65 Figure 4.4: Directional-hemispherical reflectances as functions of incident angle for a Neumann-Neumann illumination model (k^) and a Minnaert illumination model But what does this really mean? If we go back to the derivation of the geo-metrical attenuation factor in Torrance and Sparrow [75], we see that G is designed to compensate for the blocking of light that falls on a facet and the masking of light that the facet reflects. The blocking and masking agents are themselves other facets. This leads to a critical question for Torrance-Sparrow illumination models and energy conservation: What happens to the light that gets blocked or masked? The illumination model does not treat secondary reflection. Instead, it acts as though the blocked or masked light were completely absorbed by the surface. This is unlikely. For this reason, while it may be reasonable to consider the use of Torrance-Sparrow illumination models as ad hoc basis functions to fit empirical data, as was done in Larson [78], we should do so realizing that it is not really "fair" to use Torrance-Sparrow illumination models in an energy-conserving context. Basis functions that properly account for blocked and masked light are needed, but we will not attempt to derive them here. 66 4.1.3 Making Neumann-Neumann and Minnaert Illumination Models Conserve Energy Figure 4.4 shows some numerical integrations of Equation (4.7) for Neumann-Neumann and Minnaert illumination models. Contrast these with those of Figure 4.3. Like the Torrance-Sparrow illumination models, kp is bounded in both cases, so we can put a limit on bpj or &jvf to assure energy conservation. In the case of a Minnaert illumination model, we can go a bit further and note that kp can be determined analytically (using Equation (2.27), as was done in Woodham and Lee [85]). The resulting B R D F can be formulated directly in terms of A;,,: / r (S , V ) = fcpi*±i) ((N • S) (N • V ) ) ^ 1 (4.18) Z7T where, as always, any value of kp between 0 and 1 will guarantee energy conservation. 4.2 Making Illumination Models Reciprocal The second physical constraint we will examine with respect to illumination models is that of Helmholtz reciprocity. A physically plausible illumination model ought to obey Helmholtz reciprocity (see Siegel and Howell [68]). In terms of the B R D F , this means that / r ( S , V ) = / r ( V , S ) (4.19) for all V and S in £2+. 4.2.1 Are Phong Illumination Models Reciprocal? Using the B R D F of a Phong illumination model given in Equation (2.20), and expressing FS in the functional form F S ( S , V ) , we see that such an illumination 67 model will be reciprocal if F S(S,V) F S(V,S) (4.20) (N-V) (N-S) Substitution of both from Equation (2.20) and F S B from Equation (2.21) reveals that neither of these illumination models is reciprocal3. Is our energy-conserving modified Phong illumination model reciprocal? Ap-plying Equation (4.19) to Equation (4.14), we are asking if F S(S,V) F,(V,S) H,(S) . H.(V) { • 1 Again, the answer is no for both F s 's . 4.2.2 Are Torrance-Sparrow Illumination Models Reciprocal? By inspection, it is easy to see that the Torrance-Sparrow illumination models are all reciprocal. This should come as no surprise, as the assumption of reciprocity was part of their derivation in [75]. Unfortunately, the arguments made above about their energy conservation still limits their plausibility. 4.2.3 Are Separable Illumination Models Reciprocal? It is also easy to see by inspection that both Neumann-Neumann and Minnaert illumination models are reciprocal. Again, this is because reciprocity was part of their derivation. As with Torrance-Sparrow illumination models, separable illumination mod-els could be used as ad hoc basis functions. Care needs to be taken, though, to retain reciprocity. Given two separable BRDFs / ri(S, V) = ai(S)ri(V) and / r2(S, V) = «2(S)r 2(V), a simple linear combination of the form / r(S,V) = ci/ r i(S,V) + C2/r2(S, V) for some constants c\ and C2 is not, in general, reciprocal. 'Even though F? (S, V) = Ff (V,S). and Ff(S, V) = Ff(V,S). 68 We can, however, generalize the work of Westin, et al. [82] and note that if we can express the B R D F in a separable matrix product form / r ( S , V ) = [ Q ( S ) ] T M Q ( V ) (4.22) where Q(x) is some (possibly non-linear) vector function of x and M is a constant matrix, then this B R D F is reciprocal for any symmetric M . This is how we could combine multiple Neumann-Neumann, Minnaert, or other similar basis functions: set them to be the elements of Q and find a symmetric M to fit our data, as Westin, et al. did for spherical harmonics. 4.3 An Energy-Conserving, Reciprocal Illumination Model Objections can be raised to all of the illumination models we are presented so far, either on the grounds of implausibility (Phong) or of behaviour that, while plausible, is unlikely to fit a real B R D F (Torrance-Sparrow, Neumann-Neumann, Minnaert). Consider instead a Phong illumination model formulated like Equation (2.24), but using Blinn's and omitting the ( N • S) in the denominator of the specular term, we find fr(S,V) = kd + ksF?(S,V) (4.23) Obviously, since Ff is reciprocal, this B R D F is reciprocal. (We could also have done this with Fp, since it is also reciprocal.) Figure 4.5 shows the resulting kp. It is bounded, so we can always conserve energy by limiting ks and kd. (Unfortunately, we cannot formulate the illumination . model in terms of kp and k„ as we did above, since doing this makes the illumination model non-reciprocal.) 69 beta = 10 deg 0 30 60 90 Figure 4.5: Directional-hemispherical reflectance as a function of incident angle for a reciprocal Phong illumination model, using Blinn's (k^) Plausibility Conserves Other Illumination Model Energy? Reciprocal? Objections ' Phong no no Energy-Conserving Phong yes no Torrance-Sparrow yes yes no secondary reflection Neumann-Neumann yes yes non-specular behaviour Minnaert yes yes non-specular behaviour Reciprocal Phong-Blinn yes yes Table 4.1: Summary of Plausibility Results Figure 4.6 shows some images produced with a reciprocal Phong illumination model. As in Figure 2.4, ks varies between 0 and 1 and kd is taken to be 1 — ks. While resembling Figure 2.4, the images for large are dimmer, as we might expect from the absence of the (N • S) in the specular denominator. Nevertheless, they are not as diminished as those of the separable illumination models in Figure 2.5 (which does not even bother showing > 60°) . 70 p 0.00 0.25 10" 45" 90" 135' 20° 45" 9(1" 135' 30° 45" 90° 135° 40" 45° 90° 135" 50" 45" 91 0.50 0.75 1.00 II 11 ( ; ! ( I t l I ! IO B H , J l II 1 % c i • I'll 1 * up 1 w i i n 1 Figure 4.6: Sphere shaded with a reciprocal Phong illumination model for an assort-ment of incident angles with respect to the viewer specular coefficients (ks), and specular distribution half-angles (3) 4.4 Summary of Plausibility Results We have examined a number of illumination models commonly used in graphics, looking at their plausibility in terms of energy conservation and reciprocity. Our results are summarized in Table 4.1. As originally defined, Phong illumination models fail on both counts. It is possible to modify a Phong illumination model to conserve energy and even, as shown in Equation (4.14), have an energy-based parameterization, but this rules out satisfying reciprocity. Torrance-Sparrow illumination models are reciprocal and appear to conserve energy, but their underlying derivation fails to account for blocked and masked 71 energy. They may still be useful, however, as ad hoc basis functions. Neumann-Neumann and Minnaert illumination models are similar. Both are plausible: they conserve energy and are reciprocal. Minnaert illumination models have been used successfully to fit radiometric data. While it would be worth trying one of them as a basis, we expect that they will prove less useful with highly specular surfaces because both illumination models peak undesirably in the normal direction. A differently-modified Phong illumination model given in Equation (4.23) is reciprocal and can be constrained to conserve energy. 72 Chapter 5 Wavelet Radiative Transport and Surface Interaction In Section 3.8.2, we discussed how wavelets show promise for the compact repre-sentation of radiance needed by Lucifer. In this chapter, we will see that apart from compression, representing radiance in terms of a wavelet basis with direction expressed in "Nusselt coordinates" makes several calculations of relevance to illu-mination computation easier. In particular, we will show how to construct discrete representations of the radiative transport operator T (defined in Section 2.4.1) and the surface interac-tion operator S (defined in Section 2.4.3) in terms of inner products of smoothing functions. The blocked propagation operator B (defined in Section 2.4.2) can be constructed from the realization of T and will be presented in Chapter 7. Notice that these all act directly on wavelet coefficients themselves and do not require a full inverse wavelet transform. 73 0.5 K . U D C Figure 5.1: Nusselt Coordinates 5.1 Radiance in Nusselt Coordinates As we discussed in Section 2.4.1, radiance at a point P in a direction S is usually represented by a function L(P,S). If we confine our discussion to surfaces, we can assume a planar (possibly local) parameterization for P of (u, v). S is then typically represented in polar and azimuthal coordinates (9, 4>) according to the local frame of reference. Consider the x, y, and z direction cosines corresponding to a direction (9, cj>): [ix = sin 9 cos (f> Hy = sin 9 sin 4> fMz = cos 9 (5-1) We take (^x,ny) to be an alternative parameterization of direction. It is convenient in what follows for all variables to vary between the extrema of 0 and 1, so let us make a change of the directional variables from (/j,x, ny), which we used in Section 3.7 to define the U D C , to (K, A): K = M , ± 1 A = ^ + l ( 5 2 ) Figure 5.1 shows the relation between (9,cf>) and (K, A) graphically. It also shows 74 the U D C . To convert integration over (8, <fi) to integration over (K, A), the determinant of the Jacobian is: so, assuming L{ is zero for directions outside the directional limiting circle, Equa-tion (2.11) becomes ft f i That the integral no longer contains trigonometric functions should come as no surprise. We have simply used a differential form of the "Nusselt analog" ([56], but see Cohen and Wallace [14] for a description in English): the amount of power per unit area transferred from a differential solid angle du{ is proportional to dn{d\i, the area of the surface that the projection of dcoi on a unit sphere subtends. For this reason, we refer to n and A as "Nusselt coordinates". We also note that, since pix + fj,21+fi2. = l and since each vector is defined only over a hemisphere, not the whole directional sphere, we can express S+, S~, and V all unambiguously in terms of their respective incident and reflected K'S and A's. Simply put, it is always clear which sign to attach to the square root. Other ways to parameterize the directional component of a radiance distri-bution are possible. Light fields (as in Levoy and Hanrahan [44]) and lumigraphs (as in Gortler, et al. [27]), are very promising approaches for display purposes. Christensen, et al. [12] use a combination of a gnomonic projection and "stretch" to map directions to the unit square. None of these approaches, however, leads to the simplification of surface interaction that Equation (5.4) demonstrates. 0 ( K , A ) cos 8 sin 8 4 (5.3) L = Le + 4 75 5.2 Radiance Representation What are the characteristics of a four-dimensional radiance distribution L(x, y, K, A)? The easiest way to visualize this is as "light through a window" where an observer at position (x, y) on a window casts a ray in direction (K, A). For a fixed direction, the resulting two-dimensional projection is a parallel projection1, as shown in the individual cells of Figure 3.5. For a fixed position, the distribution in (K, A) would be a "fisheye" view, as shown in the individual cells of Figure 3.6. In both cases, the result is an image, so we can deal with those radiance distributions as we deal with images. Radiance at a point on a surface is a potentially discontinuous, generally non-analytic function. We can approximate it with a finite element expansion with Nf degrees of freedom: Nf L(x,y,K,\) = ^2bjBj(x,y,K,\) (5.5) i=i Choices for the basis functions Bi include box discretization (a la FIAT), Fourier, discrete cosine, spherical harmonics, orthogonal polynomials, and wavelets. We are particularly interested in wavelets because, unlike the other bases listed, their basis functions are of limited support and they can represent discontinuities compactly. They are also capable of considerable compression. 5.3 Multidimensional Wavelets Appendix A summarizes the properties of one-dimensional wavelets. In this section, we will describe multidimensional wavelets with the intention of applying them to radiative transport and surface interaction. 'The special case (K, A) = ( i , | ) is an orthographic projection. 76 For D-dimensional coordinates q = (gi, ? 2 , . - • • , ? £ ) ) , (5.6) we can define a set of multidimensional wavelet basis functions indexed by a standard multiresolution index j = {y\l\,m\,V2,m?2,... ,VD,m>D) (5.7) where i>3, which we call the "basis selector", determines the combination of one dimensional smoothing and wavelet functions: £ j ( q ) (?l)<^mj(</2) ^ ^ • ( 9 1 ) ^ ^ ( 9 2 ) • •<t>tj ML {.QD) i>j = 0 D D • •<f>ij mi (QD) V>',= 1 D D <t>i[m[ ( ? i ) ^ i m i (92) • • . ^ m , D (go) (5.8) 1 We refer to the special case of v3 — 0 as the "pure smoothing" component, as the corresponding basis function is made up of only smoothing functions. We can apply Equation (A.10) multi-dimensionally. If we have j as in Equa-tion (5.7) and define k as then D d=l d-d "*d"-d (5.9) (5.10) (Note the use of the subscripted inner product described in Section A.3.) JS; is iden-tical to 5j as defined in Equation (5.8), but with the primal wavelets and smoothing functions replaced by their duals. 77 This is the "standard" multidimensional Cartesian product basis. It is also possible to constrain l\ = V2 = ... = lJD = P, resulting in the so-called "nonstandard" basis. In general multidimensional (especially image-oriented) applications, as cited in Daubechies [16] and in Schroder et al. [65], the nonstandard bases are preferred because of their "square" support. In this paper, we are considering 4-dimensional, nonstandard basis functions, so let us enumerate the coordinates with the 4-vector q = (U,V,K,\) and the basis functions with a nonstandard multiresolution index (5.11) j = (i/,l,mu,mv,mK,mx) (5.12) Using Equations (A.l), (A.3), (A.5), (A.6) and (5.8), all the basis functions at level I of the pyramid can be written in terms of the v = 0 basis functions at level E m hm'uhmivhmiKhm^£0(/+1)(2m+m/)(q) v = 0 l^m' arn!uhm:vhmlKhm^B0^+1^2m+m'){'i) V=\ Jim'9m'uSm'vSm'liSm'xBo(l+l)(2m+m'){ci) V = 15 where m' = (m'u, m'v) m'K, m\). So for any function / , (/ | Buim)q = Y^,Wvm> | S 0 (f + l)(2m+m')) where Wvm> is (4 times) a product of smoothing and wavelet coefficients. 78 (5.13) (5.14) 5.4 Irradiance Irradiance is computed as E(x,y)= [ Li(x,y,e,<P) | N • S | du{. (5.15) Again making use of Equation (5.3), we have E(x,y) = 4f f Li(x,y,K,X)dKidXi. (5.16) Jo Jo The limits of the both integrations are 0 and 1, but if Li is zero outside the U D C (see Figure 5.1), we can safely extend the integration limits to — oo and +oo. This allows us to say that if Li(x,y,K,X) = Y,bjBi(x,y,K, A ) (5-17) j then E{x,y) = 4YJb-i{B-i\l)KX (5.18) j is the wavelet representation of the irradiance. The inner products on the right hand side are usually easy to compute in tab-ular form, if not analytically, making particular use of Equation (A.4) to eliminate many coefficients. 5.5 Power Flux The power flux passing through an area A is defined as $= / / Li(x,y,6,<f>) | N • S | duidA = f E{x,y)dA (5.19) JA 7n£ JA If we have a spatial parameterization (u,v) <-> (x,y) that maps the unit square to A and back, then we have r)ix iA dudv. (5.20) *=. /' C E(x(u,v)Mu,v)) Jo Jo o{u,v) 79 If, as above, we take E(x,y) = 0 for (x,y) outside of A and if we extend Equa-tion (5.17) to include this parameterization: Li(x,y,K,X) = ^2biBi(u(x,y),v(x,y),K,X), (5.21) j then we can incorporate Equation (5.18) to get the wavelet representation of flux: * = 4 £ * i ( S j | • (5-22) 5.6 Transport We represent radiance as where J , (q) = J > k £ k ( q ) (5.23) k bk = (L\ B k ) q (5.24) and k is defined as in Equation (5.9). Radiance travels from a source point q s to a destination point q ^ . If we have a mapping of q s —> q ^ , we can compute Ld(qd) = J24Bk(qd) (5.25) k where bi = (Ls(qs(.))\Bk = z2bIT^ (5-26) j j is defined as in Equation (5.7), and we define geometry-dependent "transport coefficients" T j k = ( f i j (g s ( - ) ) |B k ) , (5.27) 80 Using the multidimensional refinement shown in Equation (5 .13), given T(Q/ .M.) on level lj, we can compute all coefficients on the coarser level above it in the pyra-mid: m' and given Tj(oj f cmfc) on level 4, we can compute T jK(/ f c - i )m f c ) = zZW^m"Ti{0ik(2mk+m")) (5.28) (5.29) where m ? = (m3u, mJv, mj,, m3x) and m/, = (m^, mjj, m^). This means that we can compute all transport coefficients strictly in terms of pure smoothing components: (5.30) 5.7 Surface Interaction Using Equation (5.8), let us define a mixed primal-dual, four-dimensional, nonstan-dard wavelet basis: ( fern' (Ks)<t>Vm> (K)4>ljmi (Kr)<f>lJmJ (K) "J = o A K \ ^ > m ^ ( « » ) ^ > m i ( A » ) ^ / i m > («r)0,>m> (Ar) V = 1 ^ > m j > » ) V W . (A s )^- , (/t r)^m J"(A r) ^ = 2 ' (5.31) F J ( K S , A S , K R , A R ) — < where ^Umu{Ka)lPlJmi {K)^l3m3 («r )^/> m i (Ar) ^ = 15 A K A j = (^,^,ra J K ,ro\ ,r<,m^) (5.32) So Fj (K S , A S , « R , A R ) is Bj (KS, A S , K R , A R ) with dual scaling functions and wavelets substituted for the primal scaling functions and wavelets in the incident directional components only. 81 . If we then represent the B R D F in Nusselt coordinates with this basis: fr ( K S , A S , Kr, A R ) = ^2 / J F j A S , Kr, A R ) (5.33) j and Li(x,y,Ks,Xs) - Y^bkBk(x,y,ns,\s) (5.34) k where k= (vk,lk,mku,mk,mk,mkx) (5.35) then, applying Equation (5.4) and again requiring either fr or L ; (or both) to vanish outside the U D C (thus allowing us to extend our integration to (—00..00)), the reflected radiance is Lr = 4 f1 fr(S+,V)Li(S+)dKsd\a JO JO = 4 E - E J i 6 k < i ? j i B k > K . , A . ( 5 - 3 6 ) j k We seek a wavelet representation of the post-interaction radiance: Lr - ^2brnBn(x,y,Kr, Xr) (5.37) n Again using the basis representation of the reflected radiance as in Equations 5.23 and 5.24, it follows that brn = (Lr I Bn) (5.38) So, substituting Equation 5.36, we find • bn = *z2EfIb*{(F>\B*)Ks,X3\Bn) : (5.39) J k Let us now simplify notation. We can rewrite Fj and Bk compactly by defining a function Cfm(x) that takes on the value of the smoothing function or the 82 (5.40) wavelet depending on a single binary value 8: <t>im{x) 8 = 0 lplm{x) 8 = 1 and similarly for a f with <f> and tp in place of <j> and V> respectively. We also take indexable (from 0 to 3) representations of the arguments to Fj and 5k, respectively p = (KS, X S , K R , Ar) (5.41) and q = ( U , V , K S , \ S ) . (5.42) If we now adopt the notation where refers to the ctth bit of the binary form of i (i.e. i2~a mod 2), we can express the innermost inner product of Equa-tion 5.39 as 1 . j .k 3 1 =il il n (**) n <i C?KS d\s (5.43) and removing factors that do not depend on the variables of integration from the integral, we have f <°> I f <2> c"0> I f ( 3 ) .Sim' 1 Vm* As ^ ( ^ ^ ( ^ e ^ K ) ^ ^ ^ ) If we now define a multiresolution delta tensor (5.44) A C T — C (") I t (T> we can rewrite Equation 5.44 as (5.45) (5.46) 83 and we can continue to apply the A symbol to simplify the whole nested inner product in Equation 5.39 (<*S I * k > » . , A , I B*)XtVtKriXr = A g A ] 3 A™ A i i A g A g . (5.47) The As are very much dependent upon our choice of wavelet, so we will defer discussing their computation to Chapter 7. We can do the same thing with a B T D F and so represent general surface interactions: reflection, refraction, and transmission. 84 Chapter 6 Wavelet Radiative Transfer Implementation and Pract icum In this chapter, we discuss the implementation of the W R T algorithm (hereafter, "WRT'), along with a practicum ("a course of study ... that involves the super-vised practical application ... of previously studied theory"). We apply some of the concepts of Chapter 5 to a classic illumination problem: the transport of radi-ation between two arbitrarily-oriented polygons. Before that, however, we describe naming conventions, design decisions, and classes that apply to both WRT and the Lucifer implementation we discuss in Chapter 7. These implementations take up approximately 30,000 lines of C code. We therefore emphasize that the classes and pseudocode we present here are in most cases simplified for clarity. 85 Class Prefix Meaning Polygon pgn a 2D polygon with an arbitrary number of edges Transform tf a conventional 4 x 4 afhne 3D transform Transform2d tf2d a conventional 3 x 3 afnne 2D transform TransportGeometry tg geometry for W R T WaveletCoefRcientTree wet a W C T (sparse tree of wavelet coefficients) Waveletlndex wi indexes a wavelet coefficient WaveletNode wn a single node of a W C T (16 f l o a t values per channel) Table 6.1: Major WRT Classes and Their Abbreviations. 6.1 Naming Conventions Throughout WRT and Lucifer implementations, we have followed certain program-ming conventions that we feel have improved reliability and flexibility and allowed many of the benefits of object-oriented programming while developing in a highly portable environment1. We describe those conventions here to enhance the read-ability of what follows. The most visible conventions are those we use to name objects (i.e., C data structures) and their associated methods. They are an adapted form of what are known as "Hungarian" conventions2, partly in deference to the nationality of their chief developer, Charles Simonyi [70]. As implemented here, all classes have C t ypede f s associated with them. Their names are mixed case. For example, three dimensional vectors have a t y p e d e f "Vector". Each class also has a lower-case prefix. Every instance of that class (i.e., object) has a name that begins with that abbreviation. The rest of the object name describes the particular object in mixed case. An example of the "Vector" 'indeed, the code moves between IBM A1X, SGI IRIX, and Intel (RedHat) Linux with fewer than 50 system-dependent lines of code, mostly due to differing system header files. 2 See McConnell [50] for an additional discussion of Hungarian naming, particularly as practiced at Microsoft. 86 ^max ' 1 2 3 4 5 6 7 8 Number of Bits 9 1 4 1 8 2 3 2 7 3 1 3 5 4 0 ^raax 9 1 0 11 1 2 1 3 , 1 4 1 5 16 Number of Bits 4 4 4 8 5 2 5 6 6 0 6 4 6 8 7 3 Table 6 . 2 : Number of Bits Required to Store a Nonstandard Multiresolution Index as a Function of the Maximum Resolution Level class might be "vNormal". The prefix, followed by a also prefaces the names of ("member") functions that take a data structure (or a pointer to such a data structure) as a first parameter. This allows us to represent the object-oriented concept of a "method" in C. For instance, we have a function vjmagQ that returns the magnitude of its first (and only) argument, a Vector. Table 6 . 1 lists the major classes of WRT (and Lucifer) with their prefixes. 6.2 Indexing Wavelets Wavelet indices contain the six integer values referred to in Section 5 . 1 : • the basis selector v £ { 0 . . . 1 5 } , as used in Equation ( 5 . 8 ) , • the level / £ { 0 . . ./max}, where / m a x is the maximum depth of resolution, and • a vector of four offset values m with each component ra; € jo .. .I1 — 1 j (for Haar). Theoretically, the total number of bits required to represent a wavelet index is 4 + [log 2(/m a x + 1 ) ] + 4 / m a x . Table 6 . 2 shows this quantity for several possible values of / m a x . If we represented the coefficients as single channel 32-bit floating 8 7 typedef struct { short int nu\ short int /; short int m[4]; } Waveletlndex; Figure 6.1: The Waveletlndex typedef point values, then for any value of lm3uX > 6 we would use more storage for indices than for coefficients. For speed of access, it is advantageous to use data that is at least 8-bit byte-aligned. 6.3 The Waveletlndex Class Figure 6.1 shows the wavelet index typedef we have adopted for our implementation. The names and sequence of the components are consistent with Equation 5.12. For the time being, we have implemented all of these fields as short (16 bit) unsigned integers. This allows us to go as far as level wi.l = 16 without overflowing the offsets. Experience suggests that the maximum level we will use will be much less than 16. We could achieve a minor reduction in memory usage by using unsigned char variables for wi.nu3. The implemented structure requires 96 bits, which is 19 bits (31 %) larger than the theoretical minimum given in Table 6.2, but as we shall see in the next section, it is possible to group coefficients in such a way that the overhead of having to store their indices is minimal. 3but not wi.l, since it needs to represent levels from 0 to the maximum level inclusively. 88 6.4 Storing Wavelet Coefficients Because wavelet coefficients are hierarchical in nature, we refer to the data structures we devise here to store them as "wavelet coefficient trees" (hereafter, "WCT"s). Representing a W C T with a maximum level of resolution / m a x requires I 6 ' m a x + 1 possible wavelet coefficients per channel. Clearly, compression is called for. In Sec-tion A.7, we describe how wavelet properties provide an L 2-optimal compression strategy - thresholding low-magnitude coefficients. In this case, the data was single-channel and single dimensional (i.e., each wavelet coefficient's support was unique), but the results are independent of the wavelet dimensionality. We not only want to guarantee a good approximation of the data with the compressed coefficients, but also efficient use of storage and fast reconstruction. We will give an example of wavelet compression in Section 6.13. 6.4.1 Hashing Coefficients Given a sparse set of wavelet radiance coefficients {&k}> we need to store them in a way that facilitates the mapping of the wavelet index k —r bk needed to perform the transport operation Equation (5.26). The obvious way to do this is with a hash table. Not knowing the set of destination indices {k} in advance prohibits perfect hashing, so it is necessary that the hashing scheme allows for collisions and that the hash table entries contain k as well as 6k values. 6.4.2 Multichannel Grouping As we mentioned in Section 2.4.1, we have been treating data monochromatically throughout this thesis. In practical applications, however, we must evaluate Equa-tion (5.26) for all three (or however many) channels. Having assumed a non-89 participating medium, the transport coefficients Tjk are achromatic. Evaluation of Equation (5.26) simply means multiplying each element of a (now) 3-vector b- by the same value of Tjk and accumulating the result in another 3-vector b^. In the absence of the need for compression, it would be convenient to group all channels of 6| into a single group, rather than create a separate representation for each channel. In the presence of compression, however, each channel has its own threshold. If one component is above its threshold but the other two are below theirs, saving the latter is a waste of storage. This would be mitigated, however, if there were a high degree of correlation between the magnitudes of the coefficients. Certainly, a wavelet representation of white light given off by a luminaire displays a high degree of correlation. When this light is reflected off a surface of varying spectral reflectivity, however, that correlation will be diminished. By how much depends on the nature of the surface. We have two choices here: to group or not to group multichannel data. For the time being, we have chosen the former - believing that there is sufficient correlation and advantage in retrieval speed in most situations to justify grouping. This definitely requires further study. 6.4.3 Hashing Nodes Instead of Coefficients In one dimension, a node in the wavelet pyramid contains a single wavelet coefficient. In four dimensions, such a node has one pure wavelet (y = 15) coefficient and fourteen mixed wavelet/smoothing [y G {1...14}) coefficients. (Recall that the pure smoothing (y = 0) coefficient may be reconstructed from the node's ancestors, if any, and only needs to be kept at the root node.) All coefficients at a given node correspond to basis functions with the same (Haar) or similar (other bases) support. 90 typedef struct { int iWnOfHash[int mxnIWnOfHash]; struct WaveletNode { Waveletlndex wi\ unsigned short mskChild; unsigned short mskNu; int iWclFirst; } wnBase[int nWn]; struct WaveletCoefficientList { float 6 [3]; int iWclNext; } wclBase [int nWc/]; double umrad; } WaveletCoefficientTree; Figure 6.2: The WaveletCoefficientTree typedef We might therefore expect to find a higher degree of correlation in magnitude between coefficients that belong to the same node and coefficients that do not. This suggests that we can reduce the index storage overhead by hashing entire nodes rather than individual coefficients. 6.5 The WaveletCoefficientTree Class Figure 6.2 shows the fundamental structure we use to represent wavelet coefficients: WaveletCoefficientTree. Note that the code in the figure is not legal C code: we have combined several subsidiary classes and moved some component definitions around to indicate dynamically-sized structures in an obvious (we hope) way. Given a WaveletCoefficientTree wet, wet.iWnOfHashf] is the hash table itself. It is indexed by the hash of a (pure smoothing) wavelet index. Its entries are indices into wct.wnBase[], the array of wavelet nodes. Note that the use of integer indices rather than pointers allows the dynamic resizing of wct.wnBase[]. Each node contains the corresponding wavelet index. This is not only use-91 source y, parametric object destination object world 0 H, / parametric Id Figure 6.3: Coordinate Systems. The gray quadrilaterals represent the support of source and destination bases. ful for hash collision detection, but by indexing wct.wnBase[], we can traverse all nonzero entries in wet directly, wet. wnBase[i].mskChild is a 16-bit mask indicating which of node t's 16 children are also present in wet. wct.wnBase[i].mskNu is a 16-bit mask indicating which basis selectors are present in wet for node i. wct.wnBasefij.iWclFirst is the index of the first element of a list of wavelet coefficients in wct.wnBase[i].wclBase[]. The elements are ordered in increasing basis selector value, which is derived from wct.wnBase[i].mskNu. Each element wct.wnBase[i].wclBase[j] contains the index wct.wnBase[i].wclBase[j].iWclNext of the next-higher set of multichannel coefficients belonging to node i. (The last ele-ment in the list has this index set to -1.) 6.6 Representing Transport Geometry To establish the q s f-> mapping we need in order to transport radiance from source to destination, we must deal with several coordinate systems, as shown in Figure 6.3: source parametric, source object, world, destination object, and desti-nation parametric. Obviously, source object to destination object is best done with 92 a conventional affine transform (with projection), but there are several choices pos-sible for the parametric object mappings: rectilinear, perspective, and bilinear. All of these are local to the sending or receiving surface. 6.6.1 Rectilinear This is the simplest possible mapping: x W 0 u (6.1) y 0 H v W and H are the dimensions of a bounding rectangle. If the object is a rectangle, an obvious strategy is to choose coordinates in which W is its width and H is its height. This will ensure that the (inverse) transformed object completely fills the unit square. Otherwise, or in the more general case of an arbitrarily-sided polygon, when computing Tjk we must clip the (square) support of Bj (source) or i?k (destination) against the polygon when integrating. A major advantage of the rectilinear mapping is the ease of computation of the Jacobian determinant: which makes the power flux computation shown in Equation (5.22) very easy: d{x,y) d(u,v) = WE (6.3) j especially in the case of Haar wavelets: $ = AWHbQ. (6.4) 93 6.6.2 Perspective This mapping allows us to represent the more general quadrilaterals without the need to clip: xw A n A12 A13 u yw = A 2 i A 2 2 A 2 3 V w A 3 2 1 1 where the Ajj's are easily-determined functions of the quadrilateral vertices. A straight line in perspective parametric coordinates au + bv + c = 0 trans-forms to a straight line a'x + b'y + c' in object coordinates. This means that a quadrilateral in object coordinates will map to a quadrilateral in parametric coor-dinates and vice versa. This has favourable implications for transport coefficient computation that we will discuss below. One drawback of a perspective mapping is that if the quadrilateral ap-proaches degeneracy (a triangle, for instance), the appearance of a uniform grid in parametric space becomes increasingly nonuniform. The power flux computation • of Equation (5.22) is not as easy as in the rectilinear case, but may still be analytically done for basis functions Bj with closed-form representations, such as splines. 6.6.3 Bilinear As with 2-D perspective, this mapping also allows us to represent quadrilaterals without clipping. If the quadrilateral is defined by four points {po, p i , P 2 , P 3 } (in C C W order), the customary bilinear mapping applies: X 1 V Po Pl 1 - U (6.6) y P3 P2 u 94 t y p e d e f s t r u c t { Transform tfSobjWld; Transform tfDobjWld; Transform tfSobjDobj; Transform tfR; Polygon *pgnSobj; Transformed tf2dSparToSobj\ Transformed tf2dSobjToSpar; Polygon *pgnDobj;. Transformed tf2dDparToDobj; Transformed tfedDobjToDpar; } TransportGeometry; Figure 6.4: The TransportGeometry typedef Unlike the perspective case, we can treat triangles as degenerate quadrilaterals, al-beit with some irregular object space meshing. However, this mapping is not without its own drawbacks. While the parametric-to-object mapping is straightforward, the inverse object-to-parametric mapping is, in fact, double valued: a given (x,y) usu-ally has two solutions (u\v), one of them inside the unit rectangle, one outside. In order to distinguish the two cases, we must clip in object space before inversion. A more serious drawback is that straight lines are not, in general, preserved. The inverse projection of a straight line ax + by + c = 0 into parametric space' is, in general, a hyperbola. This also has implications for transport coefficient computation: it complicates determination of the limits of integration. Again, while the still more complicated Jacobian of this transform makes the power flux computation of Equation (5.22) more difficult, it may still be done for choices of JBJ with closed-form representations, such as splines. 95 6.7 The Transport Geometry Class The TransportGeometry t y p e d e f shown in Figure 6.4 contains all geometric informa-tion necessary to compute wavelet radiative transfer between source and destination polygonal objects. Given a TransportGeometry object tg, tg.tfSobjWld is the affine 3D transform from source object to world positional coordinates, tg.tfWldDobj is the affine 3D transform from world to destination object positional coordinates, tg.tfR is the orthogonal 3D rotational transform from source to destination directional coor-dinates, tg.pgnSobj and tg.pgnDobja.ve the source and destination polygons (in object coordinates), tg.tf2dSparToSobj is the affine 2D transform from source parametric to source object coordinates, and tg.tf2dDobjToDpar is the affine 2D transform from source parametric to source object coordinates. The implemented TransportGe-ometry class contains inverses of most of these transforms as well, since some of the integration schemes we use actually project destination points back into source space. 6.8 Choice of Wavelet Except where noted, our discussions in Appendix A and Chapter 5 did not depend on any particular choice of wavelet. For implementation purposes, we have to choose one. There are several good reasons for choosing Haar wavelets. As we mentioned in Section A.6, the fast wavelet transform can be performed in O(N) time, but if we allow for a varying dimensionality D and wavelet basis, it is easy to see from Equation (5.13) that the complexity is actually 0(W^N) where Wh is the (varying) maximum width of the {hj} and {hj} (and, consequently, {gj} and {gj}) coefficient sets. 96 For this reason, as the dimensionality increases, the rapidly-increasing opera-tion count makes narrower filters more and more desirable, even though wider filters generally have better approximation properties. Since Haar is the narrowest possi-ble wavelet filter (Wh — 2), it seems a wise strategy to make any multidimensional efforts first with Haar and move to wider bases later if Haar proves unsatisfactory. An additional advantage of Haar wavelets over the others is the simplifica-tion of the calculation of the transport coefficients, irradiance and power flux. As Equations (5.28) and (5.29) have shown, these coefficients can be computed entirely in terms of pure smoothing calculations. A four-dimensional Haar pure smoothing basis is a function that is constant (= 4') within a hypercube and zero outside of it. The resulting transport coefficients are volume integrals of the overlap between such a hypercube in destination parametric space and the object which is a projection of a hypercube in source parametric space into the destination space4. 6.9 Problems with Transport Coefficient Computation As the preceding subsection suggests, using Haar wavelets turns transport coefficient computation into volume integral evaluation. There are two practical problems that complicate the computation of that volume. 6.9.1 Some Source Points Do Not Project Into Destination Space The source hypercube defines a range of positional and directional coordinates q 5 . Not all of these coordinates may map to points in the destination plane, much less the destination quadrilateral. This complicates any attempt at direct evaluation of 4 This may not be such a great simplification. After all, any arbitrarily complex 3-dimensional integral may be trivially turned into a 4-dimensional volume integral! 97 the transport integrals. 6.9.2 The Projected Hypercube Has Curved Sides Even if all points in the source hypercube map to the destination plane, the nature of the resulting volume, is not trivial. Needless to say, the projection of a source para-metric hypercube into destination parametric space is not a hypercube. If, however, we could choose coordinate systems such that the source hypercube mapped to a polytope in destination space, we could take advantage of computational geometric techniques to first clip it against the destination hypercube and then compute the volume of the resulting polytope. Unfortunately, this is not possible because the source hypercube does not project to a polytope. Consider the coordinate systems described in Section 6.6. Even if we choose rectilinear parametric object mappings, the source object to destination object transform involves a projection. Hence, at best, the qs —> qd mapping has a nonlinear dependence on the directional components. As a result, the projection of the source hypercube into destination para-metric space has curved sides. Furthermore, the curvature is such that we cannot guarantee that the convex hull of the polytope formed by projecting the 16 corners of the source hypercube into destination space contains the hypervolume5. 6.10 Integration Techniques For these reasons, we must resort to multidimensional numerical integration schemes. Before considering candidate techniques, we first make an observation about the di-5 We have not fully pursued the possibility of an approximation of the projected hypercube by its convex hull here. 98 mensionality required for numerical integration. 6.10.1 Reducing the Dimensionality As Equation (5.27) indicates, the computation of transport coefficients is intrinsi-cally four-dimensional: two directional integrals and two positional integrals. If we take the two outermost integrals over direction and restrict our discussion, to pure Haar smoothing components (VJ = Vk — 0, as Equation (5.30) permits), it is then evident that T j k = 4';+'* JJs Gj(Kd, Xd) Ajk{nd, Xd) dnd dXd (6.7) 5k is the (square) directional support of Bk. Gj(nd, Xd) is a geometric function which is equal to one if a ray from the destination plane projected backwards along the ( « d i Xd) direction reaches the source plane and falls within the directional support of B-}. Otherwise, it is zero. Ajk(nd, Xd) is the area of the intersection of the spatial support of Bk in destination parametric space with the projection (in the Kd, Xd direction) to destination parametric space of the spatial support of Bj. We are now in a position to evaluate the coordinate mappings given in Sec-tion 6.6 to see which of them makes Ajk(nd, Xd) easy to compute. All the mappings transform lines of constant u or v to lines in object space, so any of them would work for the source parametric to source object mapping. Only rectilinear and per-spective mappings, however, transform arbitrary lines in object space to lines in parametric space. If we choose either of them for our destination object to desti-nation parametric mappings, computation ot A-ik(Kd,Xd) amounts to clipping the projected quadrilateral to the spatial support of Bk using a conventional polygon clipping algorithm and computing the resulting area. This allows us to reduce the dimensionality that we need to integrate numerically from four to two. 99 6.10.2 Numerical Quadrature Regardless of the number of dimensions, numerical integration techniques are all based on some form of quadrature: „ N3amp / / ( x ) d x « J2 ™>/(x') (6-8) j=l where Ns^mp is the number of samples. Techniques differ principally in their choices of weights Wi and sample points x;. Zwillinger [87] provides an extensive survey of these. The ones we have chosen to evaluate are: • trapezoidal: a regular grid approximating / linearly in each dimension • Romberg: a multilevel (Richardson) extrapolation (in terms of the grid spac-ing) of trapezoidal results • Monte Carlo: x;'s chosen from a (possibly stratified) pseudo-random sequence6 • • Halton: similar to Monte Carlo, but using quasi-random numbers generated according to number theoretical considerations • Hammersley: an alternative quasi-random method While pseudo-random and quasi-random techniques are generally preferred for multidimensional quadrature, we include trapezoidal and Romberg techniques to explore the two-dimensional case, where Monte Carlo has error of order 0(N~JJ2)7 and the trapezoidal rule has error of order 0 ( iV~ 3 / 2 ) . We must be careful to use these error estimates cautiously, however, since our integrand, G(nd, A<f) Ajk («;<*, ^d) contains discontinuities that error analyses do not account for. 6Glassner [24] is a good overview of pseudo- and quasi-random integration methods in image synthesis. 7This is true regardless of the dimensionality of the problem and is one of the appeals of Monte Carlo integration. 100 o c c 3 w > 04 0.1 0.01 10 2D Monte Carlo, Unstratified ..., Stratified 2x2 ..., Stratified 4x4 ... , Stratified 8x8 4D Monte Carlo, Unstratified ..., Stratified 2x2x2x2 ..., Stratified 4x4x4x4 Quasirandom, Halton Hammersley Romberg, Linear ... , Quadratic ... , Cubic 2D Trapezoidal 4D Trapezoidal o --X--•-X--100 1000 Time (sec.) Figure 6.5: Relative Euclidean Distance vs. Time for Different Integration Tech-niques and Numbers of Samples 6.10.3 Compar i son To evaluate the accuracies of the various methods, we have applied them to a test problem: computing the set of all Tjk for a given j and geometry. To analyze the results we can treat this set as a vector in a A'-dimensional space, where K = 15'max+i a n c j i s maximum level we have coefficients for. (I.e., we have K values of k.) To compare results against a reference, we create a reference vector T-^ using an extremely long (25 CPU-minutes) integration time and then adopt two metrics. The first, which we call the relative Euclidean distance (hereafter, RED) is the 101 3 a, Q S >< as 0.1 0.01 2D Monte Carlo, Unstratified -e-. . . , Stratified 2x2 --©--. . . , Stratified 4x4 ... , Stratified 8x8 o 4D Monte Carlo, Unstratified ... , Stratified 2x2x2x2 ... , Stratified 4x4x4x4 -•+-• Quasirandom, Halton -B-Hammersley -a-Romberg, Linear - * -. . . , Quadratic . . . , Cubic -x-2D Trapezoidal 4D Trapezoidal - * -10 100 1000 Time (sec.) Figure 6.6: Relative Maximum Departure vs. Time for Different Integration Tech-niques and Numbers of Samples Euclidean distance (L2 norm) of Tjk from T£! divided by the magnitude of T£u. Ek fe) Figure 6.5 shows how the R E D varies with integration time8 and i V 3 a m p for the various techniques. Note that we allow the degree of extrapolation for Romberg integration to vary from linear to quadratic to cubic. The times shown are for an IBM Model RS/6000 POWERserver 560 workstation. Rather than NSS,MP, time is the appropriate abscissa here. Since the time required to evaluate the integrand varies with dimensionality, fewer samples do not necessarily imply faster computation. 102 The second metric, which we call relative maximum departure (hereafter, RMD) , is more conservative. It is the maximum absolute difference between Tjk and (LK, norm) divided by the maximum absolute value of T-£: maxk Tjk - Tg maxk (6.10) Figure 6.6 shows how the R M D varies for the same parameters as Figure 6.5. From these plots, we can draw several conclusions: • The two figures are qualitatively similar: a method that does well by one metric generally does well by the other. This gives us some confidence that these metrics are valid. • Stratification improves Monte Carlo results. (This comes as no surprise.) • The time required for Monte Carlo methods being linear in . /V 3 a m p , they all approximate the expected 0(N~^p2) behaviour. • For a given integration time, the 4D trapezoidal method does generally worse than the other methods. • Increasing the degree of extrapolation for Romberg (2D) integration generally makes matters worse. This is presumably a result of the discontinuities in the integrand we discussed above. • Quasi-random methods give similar and comparatively good results. This reinforces the findings of Keller [39]. The most surprising observation, however, is how well the straightforward 2D trape-zoidal rule does. For the R E D metric, it is comparable with the best of the other 103 wct_transport() tg_oracle() tc _propagate() o \ tc_eval() tc_lntegrate() Figure 6.7: Simplified Functional Decomposition for Wavelet Radiative Transfer methods, quasi-random, and for the R M D metric it does noticeably better for short integration times. 6.11 Functional Decomposition of WRT Figure 6.7 shows a simplified functional decomposition for WRT. In this section, we will discuss each block individually and describe how it implements the ideas presented in Chapter 5. We will proceed in a bottom-up sequence. 6.11.1 tg-oracle () This function is a geometric query function. Given a TransportGeometry and source and destination Waveletlndices, it performs several fast query functions. • If either source or destination wavelet directional supports lie entirely outside the U D C , they cannot interact. Projecting the 16 vertices of the destination support hypercube project onto 104 the plane that contains the source, the vertices can be classified by the "fast reject" part of a four-dimensional version of the Cohen-Sutherland line clipping algorithm (see Foley, et al. [18], for example). • If all 16 destination vertices map into points contained within the source sup-port, we assume that the entirely of the destination support does, so an exact computation of the inner product (which is then proportional to the volume of the destination support hypercube) is possible. The latter two items require an assumption that the projected destination support vertices define a convex hull for the entire projected support. As we noted in Section 6.9.2, this is not in general the case. Nevertheless, by disabling this oracle and comparing results, we have found empirically that it is a reasonable approximation, at least for oracular purposes. 6.11.2 tc-integrate () This function is responsible for performing the actual integration. It is capable of using any of the schemes described in Section 6.10. 6.11.3 wri-.pull() and wn-push() These functions act on individual nodes of a W C T to perform, respectively, a single-level, four-dimensional (i.e. 16 value), Haar analysis (i.e., forward transform or "pull"9) or synthesis (i.e., inverse transform or "push"). 9 T h e useful "pull" and "push" terminology is from Hanrahan, et al. [32]. 105 tcjpropagate(WaveletIndex wiS, double bS[], Waveletlndex wiD, TransportGeometry tgSD, WaveletCoefficientTree wctD) if wiD.I == maximum wavelet level + 1 return tc-eval(wiS, wiD, tgSD) else if tg_oracle(wiS, wiD, tgSD) gives exact result return it else for each child wiDChild of wiD tc[wiDChild] = tc-propagate(wiS, bS, wiDChild, tgSD, wctD) for each channel chan for each child wiDChild of wiD wn[chan][wiDChild] = bS[chan] * tc[wiDChild] wn-pull(wn) if wiD.I is not the top destination level wn[chan][wiD.nu = 0] = 0 wctD[wiD] = wn Figure 6.8: tc-propagate() Pseudocode 6.11.4 tc-eval() Given source and destination wavelet indices, this function computes a single trans-port coefficient. It will call itself recursively if either source or destination basis selectors are nonzero, although as it is currently used in WRT, tcevalQ is only called to evaluate pure smoothing coefficients10. At the user's request, it will call tg.oracle() to attempt to find an alternative to numerical quadrature. Otherwise, it will call tgJntegrateQ to perform that quadrature. 6.11.5 tc-propagate() This function propagates a single source wavelet coefficient to a destination wavelet coefficient tree wctD. Figure 6.8 shows its pseudocode. Like tc-eval(), it will work with a nonzero source basis selector, but in the context of WRT is only called upon with pure source smoothing coefficients. It is important to note that wctD will l c T h e ability to work with non-zero basis selectors is used in self-test mode. 106 contain no pure smoothing coefficients except at a specified top level. The source pure smoothing coefficients corresponding to wiS are passed to the function as b0[]. There are three alternatives: • If the destination wavelet is one more than the maximum wavelet level, it returns the smoothing coefficient it gets from tceval(). • If tg_oracle() gives an exact result (often zero), it returns that result. Note that wctD does not need to be updated in this case, since the exact result is either zero or a pure smoothing result below the top level, neither of which requires storage. • Otherwise (and most importantly) the function applies itself recursively to the 16 child indices of wiD, collecting the 16 returned smoothing coefficients. For each channel, it then multiplies each coefficient by the source coefficient for that channel, and uses wn-pull() to apply a "pull" to convert the result to wavelet coefficients. If wiD.l is not the top destination level, the pure smoothing coefficient is redundant and is set to zero. The result is then ready to be stored in wctD. . 6.11.6 wct_transport() Figure 6.9 shows the pseudocode for this function, which transports an entire source wavelet coefficient tree wctS to a destination wavelet coefficient tree wctD.-The function first uses tg.oracle() to reject impossible interactions. If there is a node wn, indexed by wiS in wctS, it will call itself recursively to descend the tree. First, though, to the pure smoothing coefficients of wn, it adds bSOfJ (which is zero when this function is initially called before recursion). It then invokes wn.pushf) to 107 wct-transport(WaveletIndex wiS, double bS0[], WaveletCoefficientTree wctS, TransportGeometry tgSD, Waveletlndex wiDTop, WaveletCoefficientTree wctD): if tg-oracle(wiS, tgSD, wiDTop) says, an interaction is not possible return wn — wctS[wiS] if wn exists, for each channel chan, wn[chan][0] = w«[c/jan][0] + bS0[chan] wn-push(wn) for each child wiSChild of wiS wct-transport(wiSChild, wn[*\[wiS Child], wctS, tgSD, wiDTop, wctD) else tc-propagate(wiS, bSO, wiDTop, tgSD, wctD) Figure 6.9: wet-transport() Pseudocode convert the wavelet coefficients to finer pure smoothing coefficients at the level of the children of wiS. wet-transport () then applies itself recursively to each of these children, passing the elements of wn as new values of bS0[]. If wn does not exist in wctS, there are no nodes in wctS at or below wiS, so the only thing that needs to be propagated is the pure smoothing value bS0[]. wet-transport() calls tc-propagate() to do this. The change of prefix is an indication that below this call level, we are no longer concerned with a source W C T , but with individual pure smoothing (multichannel) coefficient vectors. 6.12 Example of Transport For a test configuration, we imagine light shining through the square stained glass window shown in Figure 6.11. As illustrated in Figure 6.10, the incident light shines down with a distribu-tion peaking at an angle of 45° from the horizontal, diffused by the glass according to a distribution proportional to the 4th power of the cosine of the angle between 108 Side View Camera View floor Figure 6.10: Geometry Used for Example of Transport the propagation direction and the peak direction. The light is transported to the floor and is collected there. Figure 6.12 shows the complete, inversely-transformed 4 D results. Each of the small images represents the spatial variation of radiance in the fixed direction given by the image's position in the matrix. Figure 6.13 is a detailed view of the brightest part of Figure 6.12. This com-putation of 3 2 x 3 2 x 3 2 x 3 2 coefficients (in red, green, and blue channels) compressed by 95% required 12.4 hours of C P U time on an I B M R S / 6 0 0 0 POWERserver 560 workstation. Needless to say, this is impractical for a single frame, but the result is reusable. Figure 6.14 shows several frames generated with an otherwise conventional raytracer modified to treat a wavelet radiance distribution as a "4 -D texture". Each frame is generated from a different camera position. Given the diffuse nature of 109 Figure 6.11: A Stained Glass Window. This is the spatial component of radiance for the test configuration. the problem, the resulting images look adequate and contain none of the "noise" common to the usual Monte Carlo approach to this sort of problem. 6.13 Example of Radiance Compression We can use the test scene to illustrate the effects of compression discussed in Sec-tion 6.4. Figure 6.15 shows the effect of compressing nodes of the wavelet represen-tation of the radiance by differing ratios prior to transport in the above example. It is evident that considerable compression is possible. Although beyond the scope of this thesis, more sophisticated reconstruction techniques could be applied to remove the artifacts visible for high compression. 110 Figure 6.12: Four-Dimensional Results of Wavelet Radiative Transport 111 Figure 6.13: Detailed View of the Brightest Part of Figure 6.12 Figure 6.14: Example of Wavelet Transport. The resulting wavelet representation of radiance shown in Figure 6.12 can be sampled from different directions. 112 Figure 6.15: Effect of Compressing the Wavelet Radiance Distribution by (left-to-right and top-to-bottom) 2:1, 20:1, 200:1, 1000:1, 2000:1, and 10000:1 113 C h a p t e r 7 Lucifer Implementation and Practicum This chapter describes the implementation of the Lucifer algorithm, which we de-scribed in Chapter 3. This brings together all previous threads of discussion. It uses both the theory of W R T described in Chapter 5 and its implementation, described in Chapter 6. In addition, as we discussed in Chapter 3, in order for Lucifer to con-verge (theoretically), it is necessary to make use of the energy-conserving shaders described in Chapter 4. As mentioned in Chapter 6, Lucifer makes use of all the classes that W R T uses. In this section, we will describe those classes that are unique to Lucifer. We will then describe the implementation choices we've made, cover a functional decomposition of it, and conclude with an example. 7.1 The Cell Class Figure 7.1 shows the typedef we use to represent cells. 114 typedef struct -Cell { BoundingBox bx; struct -Cell *parent; int Ivl; int priority; WaveletCoefficientTree *wc£[6]; 600/ isLeaf; union { struct { int nObj; Object **obj\ bool isEmShot; } leaf; struct { struct -Cell *child[8\; } nonleaf; } ujis; } Cell; Figure 7.1: The Cell t y p e d e f If cl is an instance of Cell, cl.bx is its bounding box, cl.parent is its parent cell (or N U L L if it is the root cell), cl.lvlis the cell level (0 for the root cell), cl.priority is the cell's priority used for the priority queue described in Section 3.5, cl.wctf] is an array of six W C T s , one for each wall of the cell, cl.isLeafls T R U E if the cell is a leaf cell and F A L S E if it is not. If cl is a leaf cell, cl.uJs.leaf.nObj is the number of objects it contains and cl.uAs.leaf.obj is an array of (pointers to) those objects. A flag cl.U-is.leaf.isEmShot is set depending on whether or not emissive objects within the leaf cell have had their emissivities transferred to the cell walls (i.e. "shot"). Notice that this flag is associated with an cell rather than an object. This allows emissive objects to span cells. As we shall see below, it is the responsibility of wct-transport() to assure that only that part of an object's emissivity contained within the cell being balanced gets transferred to the walls. 115 If cl is not a leaf cell, cl.uJs.nonleaf.child is an array of (pointers to) its children. 7.2 Bidirectional Reflectance and Transmittance Distri-bution Functions Like a W C T , the wavelet representation of a B R D F or B T D F is a nonstandard multiresolution 4-tensor. After the work of Lalonde and Fournier [42], we can use the same object as described in Section 6.4 with a slight reinterpretation of the meaning of the components of the wavelet indices (see Section 6.2). To represent a. B R D F or B T D F , we take the first two dimensions of the wavelet index (i.e., the ones referred to as wi.m[0] and wi.mflj) to refer to the incident directional components of the distribution and the last two dimensions (wi.m[2] and wi.m[3]) to refer to the reflected or transmitted directional components. As some time in the future, we may wish to allow for positional variation of BRDFs and B T D F s by using a multiresolution 6-tensor that adds the two additional surface positional components. It is even possible to allow for such phenomena as subsurface scattering (cf. Hanrahan and Krueger [31]) by using (in its most general form) a multiresolution 8-tensor. 7.3 Implementation Choices For reasons given in Section 6.8, we have chosen to use Haar wavelets both for WRT and for Lucifer. In this section, we will discuss the other choices we've made for the Lucifer Tenderer in the context of Section 2.5. 116 7.3.1 Integration Scheme Considering the results for WRT found in Section 6.10, it is tempting to favour use of the 2D trapezoidal method in transport coefficient computation. However, doing so for Lucifer revealed a lack of robustness. Accuracy dropped considerably when the geometry departed from the simple parallel unit square ge-ometry used for the above results. Monte Carlo and quasirandom methods did not exhibit this problem. For this reason, we chose a quasirandom integration method in Lucifer. It remains to be seen, however, if other methods are possible which might be more suitable for an integrand with such discontinuities as we have. 7.3.2 Storing Transport Coefficients Chapter 6, which considered radiative transport between two polygons, did not require the storage of transport coefficients. In other situations, however, storing them may be desirable. In this section we will outline the considerations to be made in deciding whether or not they should be stored, and if so, how? Like form factors in radiosity, the transport coefficients are dependent only upon geometry. If the same geometry occurs in several different instances, precom-puting the transport coefficients means they only have to be computed once. Counterbalancing the argument for storing coefficients is the observation that, assuming wavelet compressibility, b? is sparse. This means that only those elements of Tjk for which we have nonzero 6| values need to be computed. If we decide to store transport coefficients, we face another problem: Tjk is large. If / m a x is the maximum level of resolution on the source surface, there are 16 ' m a x + 1 possible values of j (for Haar). Assuming the same resolution is required 117 for the destination surface, Tjk would require 256 m a x + 1 elements. If represented as a conventional array, / m a x = 3 would exhaust a 32-bit address space! Fortunately, Tjk is also sparse in this sense: for a given j , the number of destination basis functions Bk that it affects is much smaller than the number of possible values of k. If we do decide to put T j k in a lookup table, that table will have to be sparse, in which case the nodal compression strategies we discussed in Section 6.4.3 would also apply to stored transport coefficients. 7.3.3 Luminaire Model The most natural luminaires to implement in Lucifer are area light sources, the reason being that since we are collecting, interacting, and re-radiating W C T s on surfaces, it is simple to add an emissive W C T prior to re-radiation. Because they are represented by W C T s , these area light sources are much more flexible than the usual area light sources. In fact, they are equivalent in functionality to Levoy and Hanrahan's [44] "light fields": Any light distribution within the limit of resolution can be represented. This is truly "light though a window". It would not be difficult to add point and directional light sources to Lucifer. A point light source would be treated like an object. When a cell contained a single point light source, its light could be discretized and wavelet compressed on the cell walls. Thereafter, it would be indistinguishable from any other light source. Linear luminaires could be treated similarly, with segmentation done when they spanned multiple cells, as could directional luminaires coming from within the root cell (as from a collimated source). Directional light sources coming from infinity (approximating sunlight, for instance) could be treated as light incident on the root cell, with the light discretized 118 and wavelet compressed on the root cell's walls. The fact that a directional light source is of uniform radiance in a single direction would lead to a very compact representation. • 7.3.4 Object Model Our initial object models in Lucifer are polygonal. Other parametric forms are possible, but the difficulty of transport coefficient computation (see Equation 5.27) dictates that the geometry for this initial effort be as simple as possible. 7.3.5 Illumination Model Lucifer's illumination model is very flexible: any B R D F or B T D F that can be represented as in Equation 5.33. The easy way to do this is to sample any B R D F or B T D F in four-dimensional Nusselt coordinates at any desired resolution, perform a (four-dimensional) Haar transform on the samples, and (for efficiency) threshold the result before storing. This only has to be done once for each distinct material and may be done prior to rendering. 7.3.6 Shading Model It is possible to create an image from Lucifer using W R T directly. We can collect light on a surface W C T placed anywhere in the scene. This would have the advantage of being usable as a light field for additional rendering. The drawback, of course, is that it would not be a very good looking image, primarily because it would not be practical to render such an image at a high enough resolution to be "visually pleasing". Just as we noted in Section 2.5.7.4, the criteria for an acceptable image differ from those for physically-correct global illumination. 119 For this reason, as part of our "final gather" (see Section 2.5.7.4) we will make use of FIATs "light update buffer" . This is a conventional ray traced image buffer whose pixel values are determined by casting primary rays onto the emissive input and iucffer-generated W C T s attached to objects in the scene. Whenever a cell is updated, all pixels whose primary rays intersect objects within that cell have their radiances incremented. 7.3.7 Surface Model In this implementation of Lucifer, we assume that reflective properties of objects do not vary with position: i.e., that texturing is constant. This means that our BRDFs and B T D F s are functions only of incident and reflected or transmitted directions. At some point in the future, it would be a natural extension to allow two additional positional degrees of freedom, but our need to reduce storage requirements demands that we defer a full treatment of texturing. 7.4 Surface Interaction with Haar Wavelets In Section 6.8, we described the benefits of Haar wavelets for WRT. In this sec-tion, we will demonstrate an additional benefit that pertains to Lucifer: the fast computation of surface interaction. Recall from Equation 5.47 that we can compute surface interaction directly from wavelet representations of the incident radiance and the B R D F or B T D F if we can compute the multiresolution delta tensors AJ^". For most wavelets, these tensors need to be computed numerically, since they rely upon such things as the inner product of a smoothing function and its dual of differing level. As we will shortly see, orthogonality helps somewhat, but not completely because of the multidimensional 120 nature of the inner products. tensors for Haar wavelets, on the other hand, do have a closed form representation: A t T T 2 - W = 0 i (mi+m")(Qm^) (Qmi)(m*+m") 2^m"=0 > i + m " ) ( g m * ) = 3.-1 u(mJa+m")(Qm*) 0 0 ^ m » = 0 ^(Qmi)(m*+"i") 7 fr (a) <r 7 fr 7 fr (<?) ( T 7 fr I// ,1 / / (<r) ( T 00 A P > / f c 00 A V < lk 01 A P > / f c 01 A V < lk 10 A V > lk 10 A V < lk 11 (7.1) where Q = 2y~l I. Of the seven cases here, the fourth, fifth, and seventh would be true for any wavelet, not just Haar. Equations 5.47 and 7.1 tell us how to evaluate the inner products for a given j , k, and n, but that is not what is required algorithmically. We start with sparse sets of coefficients |/|} and {bk} and want to generate a new set {bn}. We need to organize our implementation of Equation 5.47 in particular so that it performs the mapping (j, k) —>• {n, bn} using 7.1 to eliminate as much as possible the generation of n values which do not contribute. The pseudocode for this is lengthy and therefore included as Appendix B. Future work should consider the application of these techniques to the general problem of multidimensional inner products: /(x)y) = (j(x,-)|%1-))1 (7.2) 121 cljransportO wct_tmnsport() Figure 7.2: Simplified Functional Decomposition for Lucifer Convolution is a special case of this, although we expect it to be more useful in other cases such as spatially-variant filtering. 7.5 Functional Decomposition of Lucifer Figure 7.2 shows a simplified functional decomposition for Lucifer.. In this section, we will describe these building blocks individually. As in Section 6.11, we proceed in a bottom-up order. 7.5.1 wct-transportConfined() This is the same as wct.transport(), with one enhancement: the ability for the user to confine all transport computations to the bounding box of a given cell. This allows Lucifer to deal with objects that span cell boundaries, the need for which was described in Section 7.1. 122 Figure 7.3: Unblocked, Partially Blocked, and Completely Blocked Propagation. The source wavelet support is red, the destination wavelet support is green, and the blocking rectangle is blue. 7.5.2 wct_transportRestricted() This is another modification of wctJransport() which allows the caller to restrict both source and destination W C T s to subtrees of specified wavelet indices. Since it is always used for wall-to-wall transport, the confinement to a cell's bounding box is implicit. 7.5.3 wi-blocking() This function determines the degree to which an object is capable of blocking the transfer of radiance from the support of one Haar wavelet coefficient (and its chil-dren) to another (and its children). Figure 7.3 illustrates how blocking is determined. Given a source coefficient (on the vertical wall), and a destination coefficient (on the horizontal wall), the function constructs the "blocking hull": the polyhedron which is the convex hull bounded by their (square) spatial supports, which we may refer to as the "end caps". Once the hull is constructed, wi-blocking() determines the nature of the in-123 tersection of the hull and the object. If the object lies entirely outside the blocking hull, every point in the source coefficient support can "see" every point in the des-tination support, so the function returns NONE. If the object "slices" the blocking hull into two regions, each of which wholly contains a single end cap, no point in the source support can "see" a point in the destination support, so wLblocking() returns FULL. Otherwise, at least one source point can see at least one destination point, so blocking Type () returns PARTIAL. This query does not take into account the directional components of the supports: that is left for the actual transport computation done by wctJ,ransport() and its variants. 7.5.4 cLcfcTransblock() Given source and destination cell walls, this function, the pseudocode for which is in Figure 7.5, is responsible for transporting that part of the source wall's W C T that is not obstructed by the objects contained within the cell to the destination wall. For each object in the cell, it calls wLblocking(). If any single object totally blocks transport, no transport takes place, so it returns. If no object blocks transport at all, it calls wct-transportRestricted() to perform the actual transport and returns. Otherwise, there is partial blocking, so cLcfcTransblockf) acts recursively. If either the source or destination wavelet indices are above the maximum wavelet level, it refines the coarser of them and invokes itself on the children. If the source index is the coarser one, it also needs to push the source W C T before recurring. If the source and destination wavelet coefficients are at the maximum wavelet level, it assumes total blocking and does nothing. This rationale deserves justifica-tion. If there is a single object inside the cell, the choice is arbitrary as we are at the 124 cLcfcTransblock(Cell cl, WaveletCoeffcientTree wctSWall, TransportGeometry tgSD, WaveletCoefficientTree wctDWall, Waveletlndex wiS, Waveletlndex wiD): blocking = NONE f o r each object obj in cl, s w i t c h (wi-blocking(wiS, obj, wiD)) case NONE case PARTIAL 10 i f blocking = NONE blocking = PARTIAL case TOTAL r e t u r n / / no transport s w i t c h (blocking) case NONE wctJransportRestricted(wiS, wctSWall, tgSD, wiD, wctDWall, cl—>bx) 20 case PARTIAL i f wiS.l and wiD.I are the maximum wavelet level j j assume total blocking e l s e i f wiS.I < wiD.I wn = wctS[wiS] delete 'wctS[wiS] wn-push(wn) f o r each child wiSChild of wiS, wctS[wiSChild][wiSChild.nu = 0] = wn[wiSChild] 30 cLcfcTransblock(cl, wctSWall, tgSD, wctDWall, wiSChild, wiD) e l s e f o r each child wiDChild of wiD, cLcfcTransblock(cl, wctSWall, tgSD, wctDWall, wiS, wiDChild) Figure 7.4: cLcfcTransblock() Pseudocode maximum resolution and expect sampling effects to show up at this level in any case. If there is more than one object in the cell, that means that the cell itself is also at maximum resolution. The most likely circumstance in this case is that the objects are joined as, for example, polyhedral facets. In this case, allowing transport to take place would permit "holes" in the polyhedron, which would be an undesirable 125 wct_transportWallToObj(WaveletCoefficientTable wctSWall, Cell cl) for each object obj contained within cl, if obj has a BRDF or BTDF on its upper side, wct-transportConfined(cl.wiTop, wctSWall, tgSD, obj.wiTop, obj.wctlncomingAbove, cl.bx) if obj has a BRDF or BTDF on its lower side, wctJransportConfined(cl.wiTop, wctSWall, tgSD, obj.wiTop, obj.w'ctlncomineBelow, cl.bx) Figure 7.5: wctJransportWallToObj() Pseudocode artifact. We therefore assume total blocking in this case. 7.5.5 wctAransportWallToObj() Given a cell wall, this function, the pseudocode for which is in Figure 7.5, calls wct-transportConfined() for each object contained within the cell, restricting trans-port to the bounding box of the cell in case an object extends outside the cell. Since Lucifer allows objects to have distinct B R D F s and B T D F s on both sides of their surfaces, we need to collect incoming radiance on both sides as well. Notice that this function does not allow for blocking of one object by another. Inter-object blocking is addressed in Lucifer by choosing a maximum cell division level such that the maximum number of objects in a cell is approximately one. 7.5.6 wct_addInteraction() This is an implementation of Haar wavelet surface interaction as discussed in Sec-tion 7.4. Pseudocode for this is in Appendix B . 126 wct-transportObjToWall(Cell cl, TransportGeometry tgSD, WaveletCoefficientTree wctDWall): for each object obj contained within cl, wct_transportConfined(obj.wiTop, obj.wctOutgoingAbove, tgSD, cl.wiTop, wctDWall, cl.bx) wct-transportConfined(obj.wiTop, obj.wctOutgoingBelow, tgSD, cl.wiTop, wctDWall, cl.bx) Figure 7.6: wctAransportObjToWall() Pseudocode cl-transport(Cell cl): for each channel chan, bS0[chan] = 0 wiRoot = { 0, 0, 0, 0, "0, 0 } for each incoming wall wS of cl, for each outgoing wall wD of cl, tgSD — tg-init(cl, wS, wD) wctJransport(wiRoot, bSO, cl.wetWall[wS], tgSD, wiRoot, cl.wctWall[wD]) Figure 7.7: cLtransportf) Pseudocode 7.5.7 wct_transportObjToWall() As with wctJransportWallToObjf) (see Section 7.5.5), this function, whose pseu-docode is shown in Figure 7.6, calls wctJransportConfined() for each object con-tained within the cell, restricting transport to the bounding box of the cell in case an object extends outside the cell. In the case of radiance that has been collected on the object surfaces by wct-transportWallToObj(), this restriction is redundant, but it is necessary in the case of emissive objects, whose emissive-WCTs may extend outside the cell walls. Also like wctJransportWallToObj(), this function allows from transport from both surfaces of the object. 127 cLtransblock(Cell cl): for each incoming wall wS, for each outgoing wall wD != wS, tg = tg_init(cl, wS, wD) cLcfcTransblock(cl, cl.wetWall[wS], tgSD, cl.wctWall[wD], wiRoot, wiRoot) for each incoming wall wS, wct-transportWallToObj(cl.wctWall[wS], cl) for each object obj in cl, obj.wctOutgoingAbove = wctJnteract(obj.wctIncidAbove, obj.above.BRDF) + wctJ,nteract(obj.wctIncidBelow, obj.below.BTDF) bbj.wctOutgoingBelow = wctJnteract(obj.wctIncidAbove, obj.above.BTDF) + wctJnteract(obj.wctIncidBelow, obj.below.BRDF) for each outgoing wall wD, wctJransportObjToWall(cl, cl.wetWall[wD]) Figure 7.8: cLtransblock() Pseudocode 7.5.8 cLtransport() Figure 7.7 shows the pseudocode for this function. It is quite straightforward: for every distinct pair of source and destination walls, invoke wet^transport(). Note that since the source and destination W C T s are by definition contained within c/, it is not necessary to check for confinement, so wctJransport() can be used directly instead of wct.transportConfined(). 7.5.9 cLtransblock() Wavelets provide us with the ability to selectively refine a light representation so that we can realize the transblock() function described in Chapter 3. Figure 7.8 is the pseudocode for cLtransblock(). One of two things happens to light that enters a cell containing an object: either it passes through the cell unimpeded from incoming cell wall to outgoing cell wall or it impinges on the surface of an object and is (possibly) re-radiated to an 128 cLbalance(Cell cl): if cl contains no objects, cl-transport (cl) else cLtransblock(cl) Figure 7.9: cLbalance() Pseudocode outgoing cell wall. To deal with the first case, this function calls cLcfcTransblock() for each distinct pair of cell walls. To deal with the second case, it first uses wct_transportWallToObj() to trans-fer the radiance of all walls to all objects in the cell. Then it invokes wctJnteract() to compute the outgoing object surface W C T s , allowing for two-sided objects if specified. Finally, it uses wctJransport0.bjToWall() to transfer the object surface W C T s to the outgoing walls. 7.5.10 cLbalance() Figure 7.9 shows the pseudocode for this function, which is quite simple: To balance a cell, call cLtransportQ if the cell is empty and cLtransblockf) if it is not. 7.6 Results To demonstrate Lucifer, we will use a simple example: a source emissive square illuminating a destination triangle with an optional blocking triangle between them, as shown in Figure 7.10. Directionally, the emissivity of the square is a cosine lobe similar to 6.12, but positionally, it is uniform over the square, giving the light distribution the character of a spotlight. 129 front view side view Figure 7.10: Configuration for Lucifer Example Several discretionary variables determine how long Lucifer takes to run, how much memory it requires, and the visual quality of its output. They are: • /"" x, the maximum subdivision level of the cell walls. • whether or not the blocker is in place. This allows us to test the efficiency of our blocking algorithm. • ^mlx' the maximum subdivision level of light source and surface BRDFs. We could, in principle, select a different / m a x for each surface and light source, but we will set them all to the same value in this example to keep the number of test cases manageable. In general, it makes sense to restrict / 3 u r f < / c e " . o o ' max — max Otherwise, the higher resolution is lost during transfer. • whether or not the oracle is used at the lowest level as an alternative to numerical quadrature (as described in Section 6.11.4). We have constructed several tests cases that vary these parameters, as shown 130 Uses Low-Level Case Blocker? /surf max Oracle? 1 yes 2 no 2 yes 3 no 3 yes 3 yes 4 yes 4 no 5 no 2 no 6 no 3 yes 7 no 3 no Table 7.1: Lucifer Test Cases in Table 7.1. As it is likely to have the greatest influence on the results, we will take Z""x as the principal independent variable for each test case. We run each case until the unbalanced energy remaining in the scene is approximately 2% of the original luminaire energy. Integration in each case was done using the Halton quasirandom integration method (see Section 6.10.2) with 16 samples. We used a four-dimensional scheme instead of the two-dimensional approach discussed in Section 6.10.1. Since the need to be able to restrict transport as required by wctJ,ransportConfined() requires additional clipping to guarantee that the destination or projected source support rectangles lie within the cell, this increased the cost of the two-dimensional approach. In four dimensions, confining the integration requires a simple bounding box test of each point. As we will discuss below, however, the choices for integration scheme turn out not to be critical performance factors. Lucifer allows the user to save as a sequence of frames the intermediate contents of the light update buffer while the algorithm runs. Figure 7.11 shows these for Case 4 (Z""x = 1™*X = 4, blocker in place, no low-level oracle). In this figure, we have omitted frames that show no visible change, such as cell splits or transport through empty cells. 131 132 Figure 7.12: Lucifer Example Figure 7.12 is an enlarged view of the final frame: the end result of Lucifer. There are obvious visual artifacts, but these are comparable to those produced by radiosity and other global illumination schemes in the absence of final gather (cf. Figure 4, Color Plate 7(d) in Christensen, et al. [11]). Figures 7.13 and 7.14 show the timing and peak memory usage for all test cases examined. These were collected on a Micron Paradigm XLI workstation with a 200MHz Intel Pentium Pro C P U and 256MB of memory. Each plot shows a comparison line which is proportional to 16'™"x, which would be the space complexity in the absence of wavelet compression. 133 100000 1 i 1 'day: 1 hour' Case 1 -e— 2 + 3 • 4 x 5 6 * 7 o comparison 2 3 Maximum Cell Subdivision Level 4 Figure 7.13: Lucifer C P U Time vs. Z""x. (see Table 7.1 for test case configurations) 7.7 Analysis In order to explain the excessive resource requirements of Lucifer, we instrumented the program and made several profiling runs with various sets of parameters. We made two important discoveries. First of all, Lucifer spends a considerable amount of time using the oracle, even when the low-level oracle is turned off. As / c e " and ls"!l increase, the mean ITiaX III tlx ' depth of a source coefficient increases, and since the oracle must be consulted on each level, the amount of oracular work required per source coefficient goes up. Within the oracle, most of the time is spent projecting points from one surface to 134 100 x + i — Case 1 -e— 2 + 3 • 4 x 5 -6 x 7 © comparison 2 3 4 Maximum Cell Subdivision Level Figure 7.14: Lucifer Peak Memory Usage vs. Z""x. (see Table 7.1 for test case configurations) the other. Recall that each invocation of the oracle requires the projection of 16 points (the corners of the destination support hypercube). After we made this observation, we tried to minimize the number of pro-jections required by caching destination-to-source projections1. Unfortunately, the destination coefficient space prior to compression is not sparse, and the cache, even for Z^e"x = 4 used an exorbitant amount of memory. The second discovery we made is that when dealing with the largest cell, fewer than 2% of the integrations produced a non-zero result. Recall that in order 'Recall that integration is being done in destination space, so most projections start from there. 135 for an integration to be made (always at the finest pure smoothing level, /^e"x + 1 or 1™TJX + 1, depending on the destination), the oracle on all coarser levels must indicate that the two coefficients can interact. This low level of accuracy could in part be due to the choice of integration scheme, but since we have observed the same phenomenon with several such schemes (particularly varying the number of samples) it is much more likely due to the nature of our oracle. Attempts to devise a more accurate one have so far been unrewarding. As to Lucifer's space requirements, as we noted above, the advantage of wavelet compression is only available after we have collected destination smoothing coefficients. We cannot compress a W C T while it is being accumulated. There is another aspect in which compression doesn't help. Consider a wavelet source node, which can contain up to 16 wavelet coefficients (perhaps includ-ing a pure smoothing component from the node's parent). Recall from Section 6.11.6 that to transport the contents of the node, wct-transport() pushes the node contents and then invokes tc.propagate() on the resulting pure smoothing coefficients. As long as the node contains one non-pure smoothing coefficient, the amount of time required to transport the contents of the node is the same. While it reduces space re-quirements, compression only helps time requirements when it removes entire nodes, not individual coefficients. 136 C h a p t e r 8 Conclusions and Future Work While we have shown light-driven global illumination with wavelets to work in prin-ciple, extrapolating the Lucifer performance results in Section 7.6 to what would be considered a typical rendered scene (Nobi ~ 1000 and /^"x > 6) makes it clear that the current implementation still has a long way to go to compete with other ren-dering techniques such as photon maps in time efficiency, space efficiency, or visual quality (although the last might not be so bad with an appropriate final gather such as Lalonde and Fournier [42]). Nevertheless, in implementing Lucifer we have learned much about the wavelet representation of radiance and have made several discoveries along the way that are worthy of note: • We have examined the physical plausibility of several illumination models and suggested some useful new ones. • We have shown how a wavelet representation of radiance in Nusselt coordinates leads to the simplified computation of radiative transport and other physical quantities. 137 • We have presented a scheme for computing transport coefficients that uses a combined geometric and numerical algorithm. • We have investigated compression strategies for radiance and transport coef-ficients. • We have developed a way to compute surface interaction that acts directly on the wavelet-transformed coefficients of an incident radiance and a B R D F or B T D F . The algorithm used for this shows promise of having applicability in more general numerical inner product computations. • We have shown that the representation of radiance on a wall leads directly to an efficient representation of light fields. We have spent a considerable amount of time trying to improve the speed of Lucifer. While we have gotten speedups of better than an order of magnitude (not including faster hardware) as a result, there is clearly still much room for improve-ment. Even using the oracle, the program spends far too much time generating null smoothing coefficients and (partly as a result of the oracle itself) makes too many source-to-destination (or vice-versa) projections per coefficient. To make Lucifer viable, we believe the most promising direction for further study would be to move from a numerical/geometrical integration scheme to one that was entirely geometrical. The goal of this work would be to find an algorithm to determine the four-dimensional volume of the intersection of the projection of a source hypercube into destination space with a destination hypercube. Such an algorithm would speed up Lucifer dramatically. The "spin-off" work mentioned above done that resulted from Lucifer explo-ration has its own set of future directions. 138 8.1 Fitting Plausible Shaders to Physical Data Now that we are able to classify shaders as to their plausibility, there is a need to fit them to some real reflectance data as Larson [78] did with his own illumination model. 8.2 Fitting Plausible Shaders to Physical Data Since Lucifer deals with four-dimensional radiance representations of light intrinsi-cally, it should be possible to incorporate Levoy and Hanrahan [44] light fields or Gortler, et al. [27] lumigraphs directly into it, both as inputs and outputs., 8.3 Curved Surfaces In this thesis, for reasons of simplicity we have considered only planar surfaces. With appropriate extensions of the parametric-to-object mappings, resulting complication of the transport computations1, and some way to treat self-shadowing, it should be possible to use W R T to transport light from one curved parametric surface to another. 8.4 Importance In much the same way that importance can be incorporated into other global illu-minations schemes (cf. Smits, et al. [72] and Aupperle and Hanrahan [3], it would fit well into the Lucifer approach. We would take the observer to be an emissive 'To do this, we would probably have to give up on geometric integration schemes and fall back on entirely numerical ones. 139 source of importance (probably a point source confined to the viewing frustrum) and propagate it through cells separately from, but in the same manner as, radiance. We would modify the algorithm of Chapter 3 so that the selection of cell to be balanced would take into account not just the cells' undistributed power, but also their (integrated) importances. 8.5 Unifying Illumination Modelling and Surface Mod-elling The parameterization of BRDFs and B T D F s we have used here was four-dimensional. It is possible to incorporate the spatial variation of reflectance properties (i.e., tex-turing) by adding two additional dimensions to these distribution functions. This would represent the unification of illumination and surface modelling. Adding two more positional dimensions to represent an incident-point-to-emergent-point offset could incorporate subsurface scattering. Extending our fast inner product algorithm to work in this case would be trivial. We could even use Lucifer itself to generate such positional BRDFs and B T D F s on detailed microstructural models for later incorporation in larger-scale image synthesis. 140 A p p e n d i x A One-Dimensional Wavelet Properties In this appendix, we will review some of the properties of wavelets that make them particularly suitable for the representation of radiance. The standard reference on wavelets is Daubechies [16], from which much of this appendix is derived. We will confine ourselves here to one-dimensional wavelets. Section 5.3 will extend these properties to multiple dimensions. A . l Scaling Functions and Wavelets Wavelets are built from scaling functions, which we define by enumerated dilations and translations of a base scaling function <f>(x) of the form: <l>im(x) = 2ll24>(2lx-m) ( A . l ) 141 Figure A . l : Haar Wavelet and Smoothing Functions. Haar wavelets are the simplest possible wavelet and are the most commonly used wavelets in graphics. They have only one vanishing moment. each level / corresponds to a function space Vj, which is part of a nested sequence of subspaces . . . C V _ i C V 0 C Vi C V 2 . . - 1 . Scaling functions have the property that /+ o o (f)(x)dx / 0 (A.2) - o o We define a wavelet function space W\ as composed of those functions that need to be added to a given space Vj to span the next finer space Vj+i: V/ + i = V/ © W\. The basis functions for W\ are then dilations and translations of a mother wavelet ip(x): iblm{x) = 2l/2ip{2lx - m) (A.3) Wavelets have the property /+ o o ip(x)dx = 0 (A.4) - o o Figures A.1-A.3 show some commonly-used wavelets and their corresponding smoothing functions. 'This notation is a departure from that of Daubechies [16] that is more suitable for discrete data: an increasing subscript of V corresponds to an increasing number of samples. 142 Figure A.2: Daubechies-4 Wavelet and Smoothing Functions. Daubechies wavelets were the first compact orthogonal wavelets discovered. This particular wavelet has two vanishing moments. Figure A.3: Linear Spline Wavelet And Smoothing Functions. These wavelets are biorthogonal and have two vanishing moments. Unfortunately, while the dual wavelets and smoothing functions have a simple closed form, this is not true for the primals. A.2 Multiresolution Refinement Equations Since 4>(x) £ Vb and V 0 C V I , we can write <f>(x) as a linear combination of the basis functions 4>(2x — j) for V i : . <f>{x) = y/2y£hj<l>{2x-j) (A.5) This also holds for tp: rf;(x) = V2^2gj<t>(2x-j) (A.6) These are the dilation or refinement equations. Wavelet bases differ princi-143 pally in their choices of {hj} (which determines {gj}). Using the enumerated bases: j (A.T) j A.3 Orthogonal Wavelets We define the inner product of two functions / and g with respect to a x: f + OO (f\9)X= • f(x)g{x)dx (A.8) J — o o (Using x as a subscript is non-traditional, but will become useful when we consider multidimensional wavelets.) Some wavelets are orthogonal: But these have undesirable features: except for Haar, they are not symmetric and they do not include useful functions like splines. A.4 Biorthogonal Wavelets We can construct biorthogonal bases by using four functions instead of two: wavelets (pim and 4>im and smoothing functions ipim and ipim. These are defined defined, respectively, by four sets of coefficients: {hj}, {hj}, {gj}, and {gj}. {hj} determines {gj} and {hj} determines {gj}. If done consistently, the primal and dual components are entirely interchangeable. (A.9) 144 For these, i>lm | i>l'm') = Sll'Smmi (A.10) In the rest of this section, we'll assume the more general biorthogonality, since we can always treat orthogonal wavelets as a special case of biorthogonal wavelets. A.5 Wavelet Projections and Approximation Let us discuss the ability of a wavelet representation to approximate an arbitrary function / . Let Vif be the projection of a function / £ L? into the subspace Vj: Vif(x)=^2(f\4>im)x<f>im(x) ( A . l l ) m It can be shown I / - Vif \\2 < C2- , J V» || / W | | 2 (A.12) ' y n where Nv is the number of vanishing moments of the wavelets, i. e. for n = 0,...,Nv-l J xn^(x)dx = (xn | ^oo)^ = 0 - (A.13) What Equation (A.12) means is that we can always decrease the I? error of a wavelet approximation by increasing the number of levels / or by choosing wavelets with a higher number of vanishing moments. A.6 The Fast Wavelet Transform The fast wavelet transform, first developed by Beylkin, et al. [5], starts with a set of data s / m , m = 0 . . .Af — 1, where N = 2l. We treat these as coefficients of 4>im and can compute a wavelet transform in O(N) operations. Figure A.4 shows 145 Figure A.4: The Fast Wavelet Transform this schematically: starting on the right, each step operates on half the data of the previous step and there are / steps. This is another advantage of wavelet methods over Fourier techniques, which typically require 0(N log N) operations. The transform leaves us with 2l - 1 wavelet coefficients (dark gray in Fig-ure A.4) and 1 smoothing coefficient (light gray). The hierarchical arrangement of these coefficients is sometimes referred to as the "wavelet pyramid". A.7 Wavelet Compression Even with the wavelet pyramid, we are still dealing with N coefficients. We have not saved any storage (yet). It can be shown that if / is of the form (A.14) l,m 146 and we approximate it by F(x) = Z~2 wlm^lm{x) (A.15) where the W. is a subset of the set of coefficients {wim}, that l l / - ^ l l 2 = E K » i 2 ( A- 1 6) This gives us a convenient error metric. It also tells us that the optimal compres-sion scheme discards the coefficients with smaller magnitudes first - a thresholding process. 147 A p p e n d i x B Pseudocode for Haar Surface Interaction Figures B.1-B.7 show the pseudocode for surface interaction with Haar wavelets. This implements Equations 5.47 and 7.1 with several optimizations, particularly "short circuiting": if one of the factors in Equation 5.47 is zero, the remaining factors are not computed. 148 # computing bR[N] = sum(J,K) f[J] bS[K] < < F[J] | B[K] > \ B[N] > # This attempt takes advantage of the sparsity of both /[] and bS[]. compute-bR: for all N bR[N] = 0 for all J such that f[J] != 0 for all K such that bS[K] != 0 Delta-ks(f[J] * bS[K], J, K) Figure B . l : Pseudocode for Surface Interaction I: top level DeltaJks(s, J, K): if J.m/<0> == 0 LL K.nu<2> == 0 if K.l > J.l Q = 2'{K.l - J.l) if Q * J.m[Q] <= K.m[2] < Q * [J.m[0] + 1) DeltaJs(s / sqrt(Q), J, K) else Q = 2'(J.l - K.l) if Q * K.m[2] <= J.m[0] < Q * (K.m[2] + 1) DeltaJs(s / sqrt(Q), J, K) else if J.nu<0> =- 0 K.nu<2> == 1 if J.l > K.l Q = 2'(J.l - K.l) Qh = Q I 2 if 0 <= J.m[0] - Q * K.m[2] < Qh DeltaJs(s / sqrt(Q), J, K) else if Qh <= J.m[0] - Q * K.m[2] < Q DeltaJs(-s / sqrt(Q), J, K) else if J.nu<0> == 1 LL K.nu<2> =- 0 if J.l < K.l Q = 2-(K.l - J.l) Qh = Q I 2 if 0 <= K.m[0] - Q * J.m[2] < Qh DeltaJs(s / sqrt(Q), J, K) else if Qh <= K.m[Q] - Q * J.m[2] < Q DeltaJs(-s / sqrt(Q), J, K) else if J.l == K.l LL K.m[2] == J.m[0] DeltaJs(s, J, K) Figure B.2: Pseudocode for Surface Interaction II: DeltaJts 149 DeltaJs(s, J, K): if J.nu<\> - - 0 LL K.nu<3> == 0 if K.l > J.l Q = 2'(K.l - J.l) if 0 <= K.m[3] - Q * J.m[l] < Q DeltaJir(s / sqrt(Q), J, K) else Q = 2-(J.l - K.l) if 0 <= J.m[l] - Q * K.m[3] < Q DeltaJkr(s / sqrt(Q), J, K) else if J.n«<l> == 0 LL K.nu<3> ==1 if J.l > K.l Q = 2'(J.l - K.l) Qh = Q I 2 if 0 <= J.m[l] - Q * K.m[3] < Qh DeltaJcr(s / sqrt(Q), J, K) else if Qh <= J.m[l] - Q * K.m[3] < Q DeltaJcr(-s / sqrt(Q), J, K) else if J.n«<l> == 1 LL K.nu<3> == 0 if J.l < K.l Q = 2"(K.l - J.l) Qh = Q I 2 if 0 <= J.m[l] - Q * K.m[3] < Qh DeltaJcr(s j sqrt(Q), J, A) else if Qh <= J.m[l] - Q * K.m[3] < Q DeltaJcr(-s / sqrt(Q), J, A') else if K.l == J.l LL K.m[3] == J.m[l] DeltaJcr(s, J, K) Figure B.3: Pseudocode for Surface Interaction III: DeltaJs 150 Delta Jzr(s, J, K): for N.nu = 0 to 15 if J.nu<2> == 0 LL N.nu<2> == 0 INmax = ( N.nu == 0 ? 0 : maximum possible level ) for N.l = 0 to INmax if TV./ > J . / Q = 2-(7VJ - J./) for delM = 0 to Q - 1 N.m[2] = Q * J.m[2] + delM DeltaJr(s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) else Q = 2'(J.l - N.l) N.m[2] = floor(J.m[2] / Q) DeltaJr(s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) else if J.nu<2> == 0 LL N.nu<2> == 1 INmax = ( N.nu == 0 ? 0 : J.l - 1 ) for N.l = 0 to INmax Q = 2-(J./ - AT./) Qh = Q / 2 N.m[2] = floor(J.m[2] / Q) if 0 <= J.m[2] - Q * N.m[2] < Qh DeltaJr(s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) else DeltaJr(-s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) else if J.nu<2> == 1 LL N.nu<2> == 0 INmax = ( N.nu == 0 ? 0 : maximum possible level ) for N.l = J.l + 1 to INmax Q = 2'(J.l - N.l) Qh = Q I 2 for delM = 0 to Qh - 1 TV.m[2] = Q * J.m[2] + delM DeltaJr(s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) for delM = Qh to Q - 1 yV.m[2] = g * J.m[2] + de/M DeltaJr(-s / sqrt(Q), J, K, N.nu, N.l, N.m[2]) else if N.nu != 0 11 J.l == 0 N.l = J.l TV. ra [2] = J.m[2] DeltaJr(s, J, K, N.nu, N.l, N.m[2\) Figure B.4: Pseudocode for Surface Interaction IV: Delta_kr 151 DeltaJr(s, J, K, N.nu, N.l, N.m[2]): if J.nu<3> == 0 && N.nu<3> == 0 if N.l > J.l Q = 2'(N.l - J.l) for delM = 0 to Q - 1 N.m[3] = Q * J.m[3] + delM DeltaJX(S / sqrt(Q), J, K, N.nu, N.l, N.m[2], N.m[3]) else Q = 2~(J.l - N.l) N.m[3] = floor(J.m[3] / Q) Delta_x(s J sqrt(Q), J, K, N.nu, N.l, Nm[2], N.m[3]) else if J.nu<3> == 0 && N.nu<3> == 1 if J.l > N.l Q = 2'(J.l - N.l) Qh = Q I 2 N.m[3] = floor(J.m[3] / Q) if 0 <= J.m[3] - Q * N.m[3] < Qh Delta-x(s / sqrt(Q), J, K, N.nu, N.l, N.m[2], N.m[3]) else Delta_x(-s / sqrt(Q), J, K, N.nu, N.l, N.m[2], N.m[3]) else if J . n « < 3 > == 1 A^.nw<3> == 0 if J.l < N.l Q = 2~(N.l - J.l) Qh = Q I 2 for delM = 0 to Qh - 1 N.m[3] = Q * J.m[3] + delM Delta-x(s / sqrt(Q), J, K, N.nu, N.l, N.m[2], N.m[3]) for delM = Qh to Q - 1 N.m[3] = Q * J.m[3] + delM DeltaJC(-S j sqrt(Q), J, K, N.nu, N.l, N.m[2], N.m[3]) else if N.l - - J.l N.m[3] = J.m[3] Delta_x(s, J, K, N.nu, N.l, N.m[2], N.m[3]) Figure B.5: Pseudocode for Surface Interaction V: DeltaJr 152 Delta.x(s, J, K, N.nu, N.l, N.m[7}): if K.nu<0> == 0 LL N.nu<0> == 0 if N.l > K.l Q = 2~(N.l - K.l) for delM = 0 to Q - 1 Nm[0] = Q * K.m[0] + delM Delta..y(s / sqrt(Q), J, K, N.nu, N.l, N.m[Q], N.m[2], N.m[3]) else Q = 2~(K.l - N.l) N.m[0] = floor(K.m[0] / Q) 10 Delta.y(s j sqrt(Q), J, K, N.nu, N.l, N.m[Q], N.m[2], N.m[3]) else if K.nu<0> == 0 LL N.nu<0> == 1 if K.l > N.l Q = 2~(K.l - N.l) Qh = Q I 2 N.m[0] = floor(K.m[tr] / Q) if 0 <= K.m[Q] - Q * N.m[0] <. Qh Deltajy(s / sqrt(Q), J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3\) else Delta-y(-s / sqrt(Q), J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3]) 20 else if K.nu<0> == 1 LL N.nu<0> == 0 if K.l < N.l Q = 2~(N.l - K.l) Qh = Q j 2 for delM = 0 to Qh - 1 N.m[0] = Q * K.m[0] + delM Delta^y(s / sqrt(Q), J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3]) for delM = Qh to Q - 1 N.m[0] = Q * K.m[0] + delM Delta^y(s / sqrt(Q), J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3]) 30 else if K.l == N.l N.m[Q] = K.m[0] Delta.y(s, J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3]) Figure B.6: Pseudocode for Surface Interaction VI: Delta_x 153 Delta.y(s, J, K, N.nu, N.l, N.m[0], N.m[2], N.m[3]): i f K.nu<\> == 0 LL N.nu<\> 0 i f N.l > K.l Q = 2~(N.l - K.l) for delM = 0 to Q - 1 N.m[l] = Q * K.m[l] + delM bR[N] = bR[N] + s / sqrt(Q) else Q = 2-(K.l - N.l) N.m[l] = floor(K.m[l] / Q) • 10 bR[N] = bR[N] + s / sqrt(Q) else i f 7\.nw<l> == 0 LL N.nu<l> == 1 i f K.l > N.l Q = 2~(K.l - N.l) Qh = Q I 2 N.m[l] = floor(K.m[l] / Q) i f .0 <= K.m[l] - Q * N.m[l] < Qh bR[N] = bR[N] + s j sqrt(Q) else bR[N] = bR[N] - s I sqrt(Q) 20 else i f K.nu<l> == 1 LL N.nu<l> == 0 i f K.l < N.l Q = 2~(Nl - K.l) Qh - Q I 2 for delM = 0 to Qh - 1 N.m[l] = Q * K.m[l] + delM bR[N] = bR[N] + s / sqrt(Q) for delM = Qh to Q - 1 N.m[l] = Q * K.m[l] + delM bR[N] = bR[N] - s / sqrt(Q) 30 else i f K.l == N.l N.m[l] = K.m[l] bR[N] = bR[N] + s Figure B.7: Pseudocode for Surface Interaction V I I : Delta_y 154 Bibliography [1] James Arvo and David Kirk. A survey of ray tracing acceleration techniques. In Andrew S. Glassner, editor, An Introduction to Ray Tracing, pages 201-262. Academic Press, 1989. [2] Ian Ashdown. Near-field photometry: The helios approach. In Graphics Inter-face '92 Workshop on Local Illumination, pages 1-14, May 1992. [3] Larry Aupperle and Pat Hanrahan. Importance and discrete three point trans-port. In Michael F . Cohen, Claude Puech, and Francois Sillion, editors, Fourth Eurographics Workshop on Rendering, pages 85-94. Eurographics, June 1993. held in Paris, France, 14-16 June 1993. [4] P. Beckmann and A. Spizzichino. The scattering of electromagnetic waves from rough surf aces. Artech House Inc, 1987. [5] G . Beylkin, Ronald R. Coifman, and V . Rokhlin. Fast wavelet transforms and numerical algorithms i. Communications on Pure and Applied Mathematics, 44:141-183,1991. [6] James F . Blinn. Models of light reflection for computer synthesized pictures. In James George, editor, Computer Graphics (SIGGRAPH '77 Proceedings), volume 11, pages 192-198, San Jose, California, July 1977. [7] James F . Blinn. Simulation of wrinkled surfaces. In Computer Graphics (SIG-GRAPH '78 Proceedings), volume 12, pages 286-292, August 1978. [8] OpenGL Architecture Review Board, Mason Woo, Jackie Neider, and Tom Davis. OpenGL Programming Guide. Addison Wesley Developers Press, 1997. [9] Edwin E . Catmull. A Subdivision Algorithm for Computer Display of Curved Surfaces. PhD thesis, Dept. of CS, U. of Utah, December 1974. [10] S. Chandrasekar. Radiative Transfer. Oxford Univ. Press, 1950. 155 [11] Per H. Christensen, Eric J . Stollnitz, David H. Salesin, and Tony D. DeRose. Wavelet radiance. In Fifth Eurographics Workshop on Rendering, pages 287-302, Darmstadt, Germany, June 1994. [12] Per Henrik Christensen, Eric J . Stollnitz, David H. Salesin, and Tony D. DeRose. Global illumination of glossy environments using wavelets and im-portance. ACM Transactions on Graphics, 15(1):37-71, January 1996. [13] F . J . J . Clarke and D. J . Parry. Helmholtz reciprocity: its validity and appli-cation to reflectometry. Lighting Research & Technology, 17(1):1—11, 1985. [14] Michael F . Cohen and John R. Wallace. Radiosity and Realistic Image Synthe-sis. Academic Press Professional, San Diego, C A , 1993. [15] R. L. Cook and K. E . Torrance. A reflectance model for computer graphics. ACM Transactions on Graphics, l(l):7-24, January 1982. [16] Ingrid Daubechies. Ten Lectures on Wavelets, volume 61 of CBMS-NSF Re-gional Conference Series in Applied Mathematics. SIAM, Philadelphia, PA, 1992. [17] Eugene L. Fiume. The Mathematical Structure of Raster Graphics. Academic Press, April 1989. [18] James D. Foley, Andries van Dam, Steven K. Feiner, and John F . Hughes. Computer Graphics, Principles and Practice, Second Edition. Addison-Wesley, Reading, Massachusetts, 1990. Overview of research to date. [19] A . Fournier. Separating reflection functions for linear radiosity. In Eurographics Rendering Workshop 1995. Eurographics, June 1995. [20] A . Fournier, E . Fiume, M . Ouellette, and C . K. Chee. FIAT L U X : Light driven global illumination. Technical Memo DGP89-1, Dynamic Graphics Project, Department of Computer Science, University of Toronto, 1989. [21] Alain Fournier. Normal distribution functions and multiple surfaces. In Graph-ics Interface '92 Workshop on Local Illumination, pages 45-52, May 1992. [22] Alain Fournier. From local to global illumination and back. In P. M . Hanrahan and W . Purgathofer, editors, Rendering Techniques '95 (Proceedings of the Sixth Eurographics Workshop on Rendering), pages 127-136, Dublin, Ireland, June 1995. Springer-Verlag. 156 [23] Akira Fujimoto, Takayuki Tanaka, and Kansei Iwata. Arts: Accelerated ray-tracing system. IEEE Computer Graphics and Applications, pages 16-26, April 1986. [24] Andrew S. Glassner. Principles of Digital Image Synthesis. Morgan Kaufmann, 1995. [25] Cindy M . Goral, Kenneth E . Torrance, Donald P. Greenberg, and Bennett Bat-taile. Modelling the interaction of light between diffuse surfaces. In Computer Graphics (SIGGRAPH '84 Proceedings), volume 18, pages 212-22, July 1984. [26] Steven J . Gortler, Michael F . Cohen, and Phillip Slusallek. Radiosity and relaxation methods: Progressive refinement is Southwell relaxation. Technical Report CS-TR-408-93, Department of Computer Science, Princeton University, February 1993. to appear in I E E E Computer Graphics & Applications. [27] Steven J . Gortler, Radek Grzeszczuk, Richard Szeliski, and Michael F . Cohen. The lumigraph. In Holly Rushmeier, editor, Proceedings of SIGGRAPH '96 (New Orleans, LA, August 4-9, 1996), Computer Graphics Proceedings, Annual Conference Series, pages 43-54. A C M S I G G R A P H , A C M Press, August 1996. [28] Steven J . Gortler, Peter Schroder, Michael F . Cohen, and Pat Hanrahan. Wavelet radiosity. In Computer Graphics Proceedings, Annual Conference Se-ries, 1993, pages 221-230, 1993. [29] Eric A . Haines. A proposal for standard graphics environments. IEEE Com-puter Graphics and Applications, 7(11):3-5, November 1987. [30] Roy Hall. Illumination and Color in Computer Generated Imagery. Springer-Verlag, New York, 1989. [31] Pat Hanrahan and Wolfgang Krueger. Reflection from layered surfaces due to subsurface scattering. In James T. Kajiya, editor, Computer Graphics (SIG-GRAPH '93 Proceedings), volume 27, pages 165-174, August 1993. [32] Pat Hanrahan, David Salzman, and Larry Aupperle. A rapid hierarchical ra-diosity algorithm. In Thomas W . Sederberg, editor, Computer Graphics (SIG-. GRAPH '91 Proceedings), volume 25, pages 197-206, July 1991. [33] Xiao D. He, Kenneth E . Torrance, Francois X. Sillion, and Donald P. Greenberg. A comprehensive physical model for light reflection. In Thomas W. Sederberg, editor, Computer Graphics (SIGGRAPH '91 Proceedings), volume 25, pages 175-186, July 1991. 157 [34] Berthold Klaus Paul Horn. Robot Vision. McGraw-Hill, 1986. [35] David S. Immel, Michael F . Cohen, and Donald P. Greenberg. A radiosity method for non-diffuse environments. In David C. Evans and Russell J . Athay, editors, Computer Graphics (SIGGRAPH '86 Proceedings), volume 20, pages 133-142, August 1986. [36] American National Standards Institute. ANSI standard nomenclature and def-initions for illuminating engineering,. ANSI/IES RP-16-1986, Illuminating En-gineering Society, 345 East 47th Street, New York, N Y 10017, June 1986. [37] Henrik Wann Jensen. Global illumination using photon maps. In Seventh Eurographics Workshop on Rendering, pages 22-30, Porto, Portugal, June 1996. [38] James T . Kajiya. The rendering equation. In David C . Evans and Russell J . Athay, editors, Computer Graphics (SIGGRAPH '86 Proceedings), volume 20, pages 143-150, August 1986. [39] Alexander Keller. Quasi-monte carlo radiosity. In Seventh Eurographics Work-shop on Rendering, Porto, Portugal, June 1996. [40] Renate Kempf and Chris Frazier, editors. OpenGL Reference Manual. Addison Wesley Developers Press, 1997. [41] Paul Lalonde and Alain Fournier. Filtered local shading in the wavelet domain. In J . Dorsey and Ph. Slusallek, editors, Rendering Techniques '97 (Proceedings of the 8th Eurographics Workshop on Rendering), pages 163-174. Springer-Verlag, June 1997. [42] Paul Lalonde and Alain Fournier. A wavelet representation of reflectance func-tions. IEEE Transactions on Visualization and Computer Graphics, 3(4):329-336, October-December 1997. [43] J . Lane, L. Carpenter, T . Whitted, and J . Blinn. Scan line methods for display-ing parametrically defined surfaces. Communications of the ACM, 23(l):23-34, 1980. [44] Marc Levoy and Pat Hanrahan. Light field rendering. In Holly Rushmeier, editor, Proceedings of SIGGRAPH '96 (New Orleans, LA, August 4-9, 1996), Computer Graphics Proceedings, Annual Conference Series, pages 31-42. A C M S I G G R A P H , A C M Press, August 1996. 158 [45] Robert Lewis. Making shaders more physically plausible. In Michael F . Cohen, Claude Puech, and Francois Sillion, editors, Fourth Eurographics Workshop on Rendering, pages 47-62. Eurographics, June 1993. held in Paris, France, 14-16 June 1993. [46] Robert R. Lewis. Solving the classic radiosity equation using multigrid tech-niques. In Proceedings of the 1992 Western Computer Graphics Symposium, pages 157-164, Banff, Alberta, Canada, April 1992. [47] Robert R. Lewis. Making shaders more physically plausible. Computer Graphics Forum, 13(2):109 - 120, June 1994. [48] William E . Lorensen and Harvey E . Cline. Marching cubes: A high resolution 3D surface construction algorithm. In Maureen C . Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 163-169, July 1987. [49] Martti J . Mantyla. Introduction to Solid Modeling. Computer Science Press, Rockville, M D , December 1986. [50] Steve McConnell. Code Complete. Microsoft Press, 1993. [51] Dimitri Mihalas. Stellar Atmospheres. W. H . Freeman, San Francisco, 1978. [52] M . Minnaert. The reciprocity principle in lunar photometry. Astrophysical Journal, 93:403-410, May 1941. [53] Curtis D. Mobley. Light in Water: Radiative Transfer in Natural Waters. Academic Press, San Diego, 1994. [54] Hans Moravec. 3D graphics and the wave theory. In Computer Graphics (SIG-GRAPH '81 Proceedings), volume 15, pages 289-96, August 1981. [55] L . Neumann and A. Neumann. Photosimulation: Interreflection with arbitrary reflectance models and illuminations. Computer Graphics Forum, 8(l):21-34, March 1989. [56] W. Nusselt. Graphische Bestimmung des Winkelverhaltnisses bei der Warme-Strahlung. Zeitschrift des Vereines Deutscher Ingenieure, 72(20):673, 1928. [57] Michael Oren and Shree K. Nayar. Generalization of Lambert's reflectance model. In Andrew Glassner, editor, Proceedings of SIGGRAPH '94 (Orlando, Florida, July 24~29, 1994), Computer Graphics Proceedings, Annual Confer-ence Series, pages 239-246. A C M S I G G R A P H , A C M Press, July 1994. ISBN 0-89791-667-0. 159 [58] Bui-T. Phong. Illumination for computer generated pictures. Communications of the ACM, 18(6):311-317, June 1975. [59] Pierre Poulin and John Amanatides. Shading and shadowing with linear light sources. In C. E . Vandoni and D. A. Duce, editors, Eurographics '90, pages 377-386. North-Holland, September 1990. [60] Pierre Poulin and Alain Fournier. A model for anisotropic reflection. In Forest Baskett, editor, Computer Graphics (SIGGRAPH '90 Proceedings), volume 24, pages 273-282, August 1990. [61] Mark C. Reichert. A two-pass radiosity method driven by lights and viewer position. Master's thesis, Program of Computer Graphics, Cornell University, January 1992. [62] Lawrence G . Roberts. Machine perception of three-dimensional solids. T R 315, Lincoln Lab, MIT, Lexington, M A , May 1963. [63] Holly Rushmeier. Rendering participating media: Problems and solutions from application areas. In Fifth Eurographics Workshop on Rendering, pages 35-56, Darmstadt, Germany, June 1994. [64] Holly E . Rushmeier and Kenneth E . Torrance. Extending the radiosity method to include specularly reflecting and translucent materials. ACM Transactions on Graphics, 9(l):l-27, January 1990. [65] Peter Schroeder, Steven Gortler, Michael Cohen, and Pat Hanrahan. Wavelet projections for radiosity. In Michael F . Cohen, Claude Puech, and Francois Sillion, editors, Fourth Eurographics Workshop on Rendering, pages 105-114. Eurographics, June 1993. held in Paris, France, 14-16 June 1993. [66] Peter Schroeder and Pat Hanrahan. Wavelet methods for radiance computa-tions. In Fifth Eurographics Workshop on Rendering, pages 303-311, Darm-stadt, Germany, June 1994. [67] Michael Shantz and Sheue-Ling Lien. Shading bicubic patches. In Maureen C. Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 189-196, July 1987. [68] Robert Siegel and John R. Howell. Thermal Radiation Heat Transfer. Hemi-sphere Publishing Corp., Washington, D C , 1981. 160 [69] Francois X. Sillion, James R. Arvo, Stephen H. Westin, and Donald P. Green-berg. A global illumination solution for general reflectance distributions. In Thomas W. Sederberg, editor, Computer Graphics (SIGGRAPH '91 Proceed-ings), volume 25, pages 187-196, July 1991. [70] Charles Simonyi. MetaTprogramming: A software production method. Techni-cal Report CSL-76-7, Xerox Palo Alto Research Center, 1976. [71] Brian Smits, James Arvo, and Donald Greenberg. A clustering algorithm for radiosity in complex environments. In Andrew Glassner, editor, Proceedings of SIGGRAPH '94 (Orlando, Florida, July 24~29, 1994), Computer Graph-ics Proceedings, Annual Conference Series, pages 435-442. A C M S I G G R A P H , A C M Press, July 1994. ISBN 0-89791-667-0. [72] Brian E . Smits, James R. Arvo, and David H. Salesin. An importance-driven radiosity algorithm. In Edwin E . Catmull, editor, Computer Graphics (SIG-GRAPH '92 Proceedings), volume 26, pages 273-282, July 1992. [73] Barton T . Stander and John C . Hard. Guaranteeing the topology of an implicit surface polygonalization for interactive modeling. In SIGGRAPH '97 Confer-ence Proceedings, pages 279-286, Los Angeles, California, August 3-8 1997. [74] Toshimitsu Tanaka and Tokiichiro Takahashi. Shading with area light sources. Eurographics '91, pages 235-46, 535-7, September 1991. [75] K. E . Torrance and E . M . Sparrow. Theory for off-specular reflection from roughened surfaces. Journal of Optical Society of America, 57(9), 1967. [76] T . S. Trowbridge and K. P. Reitz. Average irregularity representation of a roughened surfaces for ray reflection. Journal of Optical Society of America, 65(3), 1967. [77] C. P. Verbeck and D. P. Greenberg. A comprehensive light source description for computer graphics. IEEE Computer Graphics and Applications, 4(7):66-75, July 1984. [78] Gregory J . Ward. Measuring and modeling anisotropic reflection. In Edwin E . Catmull, editor, Computer Graphics (SIGGRAPH '92 Proceedings), volume 26, pages 265-272, July 1992. [79] J . Warnock. A hidden-surface algorithm for computer generated half-tone pic-tures. Technical Report T R 4-15, NTIS AD-733 671, University of Utah, Com-puter Science Department, 1969. 161 [80] Alan Watt. Fundamentals of Three-Dimensional Computer Graphics. Addison-. Wesley, Wokingham, England, 1989. [81] K. Weiler and K. Atherton. Hidden surface removal using polygon area sorting. Computer Graphics (SIGGRAPH '77 Proceedings), l l (2 ) :214-222, July 1977. [82] Stephen H. Westin, James R. Arvo, and Kenneth E. Torrance. Predicting re-flectance functions from complex surfaces. In Edwin E. Catmull, editor, Com-puter Graphics (SIGGRAPH '92 Proceedings), volume 26, pages 255-264, July 1992. [83] Turner Whitted. An improved illumination model for shaded display. Commu-nications of the ACM, 23(6):343-349, June 1980. [84] Lance Williams. Casting curved shadows on curved surfaces. In Computer Graphics (SIGGRAPH '78 Proceedings), volume 12, pages 270-274, August 1978. [85] R. J. Woodham and T. K. Lee. Photometric method for radiometric correction of multispectral scanner data. Canadian Journal of Remote Sensing, 11(2): 132-161, December 1985. [86] Brian Wyvill, Craig McPheeters, and Geoff Wyvill. Data structure for soft objects. The Visual Computer, 2(4):227-234, 1986. [87] Daniel Zwillinger. Handbook of Integration. Jones and Bartlett Publishers, Inc., Boston, MA, 1992. 162 Index ANSI/IES, 6, 10, 60 bidirectional reflectance distribution func-tion (BRDF), 5, 14, 16, 25-28, 31-33, 37-39, 47, 59, 63, 67-69, 82, 116, 119, 120, 126, 137 bidirectional transmittance distribution function (BTDF) , 5, 15, 16, 37, 47, 84, 116, 119, 120, 126, 137 B R D F , see bidirectional reflectance dis-tribution function B T D F , see bidirectional transmittance distribution function FIAT, 41, 54, 55, 76, 120 final gather, 38, 120. illumination global, 1, 8, 35, 38-41 local, 8, 13, 16, 36, 38, 40, 44 model, 17, 23-33, 36, 37, 40 illumination model, 14, 44, 48, 58-72, 119 index multiresolution, 77, 78 wavelet, see wavelet indices Lucifer (algorithm), 39, 41-57 Lucifer (implementation), 85-87,114-136 Nusselt coordinates, 74-75, 82, 119, 137 • shading model, 17, 33-34, 119-120 subsurface scattering, 32, 116 surface interaction, 5, 16, 35, 37, 47, 58, 75, 76, 81-84, 120-122, 126, 137, 148 model, 17, 32-33, 120 U D C , see unit directional circle 163 unit directional circle, 52, 74, 79, 82, 104 wavelet basis, 54 coefficient trees (WCTs), 89-92, 115, 118, 119, 123, 124, 127-129 compression, 50, 117, 118 indices, 86-89, 91, 116, 123, 124 luminaires, 19 radiance, 55-57 radiative transport (WRT), 73-86 radiosity, 38 refinement, 128 wavelets Haar, 116, 120-123, 126, 148 one-dimensional, 141-147 W C T , see wavelet coefficient trees W R T (algorithm), see wavelet radia-tive transport WRT (implementation), 85-110, 116, 117, 120 164
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Light-driven global illumination with a wavelet representation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Light-driven global illumination with a wavelet representation of light transport Lewis, Robert R. 1998
pdf
Page Metadata
Item Metadata
Title | Light-driven global illumination with a wavelet representation of light transport |
Creator |
Lewis, Robert R. |
Date Issued | 1998 |
Description | This thesis considers the problem of global illumination: the modelling of light as it travels through a scene interacting with the objects contained within the scene. Starting with a description of the problem and a discussion of previous work, we explore a new approach called light-driven global illumination that offers several advantages over its predecessors: a lower asymptotic complexity, a wider range of representable surface interaction phenomena, and an absence of the need for "meshing" - object surface subdivision needed primarily to represent shadows. Light-driven global illumination is intermediate between local and global illumination. Representing light with wavelet basis functions, we are able to treat both the interaction between two surfaces and the interaction of a surface with a radiation field in a source-to-destination model that applies to whole surfaces, not just small elements. We have found this "wavelet radiative transfer" to be a valid way to generate and store complex global light field data as four-dimensional textures for incorporation in local illumination solutions. Wavelets can considerably reduce the otherwise substantial storage and reconstruction problems associated with doing this. We include several examples of this. We also discuss plausible illumination models, which are required to make light-driven global illumination work theoretically. Like wavelet radiative transfer, these models have application in other areas of rendering besides global illumination. Finally, we develop the theory behind light-driven global illumination and apply it successfully to some simple examples. While we find the algorithm to be quite slow compared to other well-known rendering algorithms, we analyze what is needed to make it competitive. In conclusion, we find that representing light with wavelets has a set of advantages that are independent of the comparative inefficiency of the light-driven algorithm. |
Extent | 13645814 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051615 |
URI | http://hdl.handle.net/2429/9529 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1998-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1998-345793.pdf [ 13.01MB ]
- Metadata
- JSON: 831-1.0051615.json
- JSON-LD: 831-1.0051615-ld.json
- RDF/XML (Pretty): 831-1.0051615-rdf.xml
- RDF/JSON: 831-1.0051615-rdf.json
- Turtle: 831-1.0051615-turtle.txt
- N-Triples: 831-1.0051615-rdf-ntriples.txt
- Original Record: 831-1.0051615-source.json
- Full Text
- 831-1.0051615-fulltext.txt
- Citation
- 831-1.0051615.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051615/manifest