Wavelet Radiosity in Computer Graphics by P h i l i p p Ziegler D i p l . I n f o r m . , Univers i ta t Kaisers lautern , G e r m a n y , 1996 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Depar tment of C o m p u t e r Science) we accept this thesis as conforming to the required s tandard The University of Br i t i sh Columbia September 1998 © P h i l i p p Ziegler, 1998 In presenting this thesis/essay in partialfuCfiCCment of the requirements for an advanced degree at the University of'British CoCumbia, I agree that the Library shaCC make it freeCy avaiCabCe for reference and study. I further agree that permission for extensive copying for this thesis for sc ho CarCy purposes may be granted by the 3-Cead of my department or by his or her representatives. It is understood that copying or pubdcation of this thesis for financiaCgain shaCCnot be aCCowed without my written permission. se(i 2?, tm Date Computer Science The University of'British Co(um6ia 2366 Main maCC "Vancouver, "EC Canada yeT 1Z4 Abstract T h e thesis presents an overview of the recent development of radiosity methods in computer graphics in a un i form m a t h e m a t i c a l framework. T h e focus is on hierarchical methods using wavelets. T h e thesis experimental ly analyzes the behavior of higher-order wavelet bases in hier-archical methods . T h e functions appl ied are multiwavelets and a fami ly of wavelets proposed by C o h e n , Daubechies and V i a l . T h e latter wavelets have over lapping support . General ly , we find that higher-order wavelet bases save m e m o r y c o m p a r e d to H a a r wavelets, while they require more t ime for the c o m p u t a t i o n . Furthermore , we investigate how K r y l o v subspace methods can be employed to solve the discrete system of equations aris ing in hierarchical methods . W e show that the G e n e r a l -ized M i n i m a l Res idua l m e t h o d ( G M R E S ) is advantageous c o m p a r e d to the usual ly employed P i c a r d i teration. 1 1 Contents Abstract ii Contents iii Acknowledgements vi Dedication vii 1 Introduction 1 2 The Radiosity Problem and Its Numerical Solution 4 2.1 T h e Rad ios i ty M o d e l . . 4 2.2 Di scre t i z ing the Radios i ty E q u a t i o n 8 2.2.1 F u n d a m e n t a l P r o b l e m s 8 2.2.2 Bas ic M e t h o d s 9 2.2.3 Hierarch ica l M e t h o d s 13 2.2.4 Discont inu i ty M e s h i n g 23 2.2.5 H y b r i d M e t h o d s 24 2.3 So lut ion M e t h o d s 25 2.3.1 Propert ies 25 i i i 2.3.2 Gauss-Seidel-like Methods 26 2.3.3 Other Methods • • 33 2.3.4 Comparison 34 2.4 Computing the Discrete System 35 2.4.1 Form Factor Approximation 36 2.4.2 Form Factor Estimation 39 2.4.3 Comparison 39 2.5 Summary 40 3 Exploring Wavelet Radiosity 42 3.1 Wavelet Bases for Radiosity 42 3.1.1 Haar Wavelets 43 3.1.2 Multiwavelets 44 3.1.3 C D V Wavelets 46 3.1.4 Comparison 48 3.2 Solution Methods for Wavelet Radiosity 52 3.3 Implementation Aspects 52 3.3.1 Galerkin Discretization 53 3.3.2 Form Factor Estimate 53 3.3.3 Form Factor Computation • • 54 3.3.4 Memory Representation of the Sparse Matrix 55 3.3.5 Solution Method 56 4 Numerical Experiments 57 4.1 Comparison of Wavelet Bases 57 4.1.1 Test Scenes 59 iv 4.1.2 F o r m F a c t o r C o m p u t a t i o n . . . 5 9 : 4.1.3 Measurements 60 4.1.4 Interpretat ion 61 4.2 C o m p a r i s o n of So lut ion M e t h o d s 67 4.2.1 Tes t Scenes 67 4.2.2 E x p e r i m e n t s • • • 69 4.2.3 Interpretat ion 72 5 Conclusions 73 5.1 Results :-. 73 5.1.1 Wavele t Bases 73 5.1.2 So lut ion M e t h o d s 74 5.2 F o r F u r t h e r Investigation 75 Bibliography 76 Appendix A 82 A . l T w o - D i m e n s i o n a l Wavelet T r a n s f o r m 82 A . 2 N o n - S t a n d a r d Representat ion 83 v Acknowledgements I would like to thank Professor U r i Ascher for his work as m y research supervisor, in part icu lar for his fairness and for prov id ing a very fr iendly research environment . Professor A l a i n Fourn ier gave construct ive cr i t i c i sm and pointed out valuable resources. O n e of these being Ian Ashdown's g lobal i l l u m i n a t i o n b ib l iography wi thout which wr i t ing this thesis would have taken much longer. L a s t but not least I would like to ment ion L o r i n e C h a n g , who endured all m y compla in ing , Dave Graves , who pat ient ly answered all m y questions regarding E n g l i s h style and g r a m m a r , as well as M i k e Horsch a n d A s h l e y W i j e y e r a t n a m , who convinced me that a thesis has to be complete rather t h a n perfect. P H I L I P P Z I E G L E R The University of British Columbia September 1998 vi to /V\0U5Y vii Chapter 1 Introduction The problem of global illumination is of major importance in computer graphics. Global illumination is concerned with computing radiance (i.e. brightness values) in a scene from its geometrical and physical description. This description contains objects (such as walls, furniture etc.) as well as light-sources. Global illumination problems can widely vary in size, ranging from illumination of single rooms with box-like furniture to highly complex scenes of entire buildings with the environment modeled into detail (e.g. [tell94]). A well known method for the solution of a global illumination problem is ray-tracing, which is concerned with the computation of radiance of the area of the surface of the scene seen from one given viewpoint ([whit80]). In contrast, methods based on Galerkin discretiza-tions (e.g. radiosity methods) compute radiance for each point of the surface. This allows for rapid change between views without expensive re-computations. Solutions of this kind are useful in architectural design and animations. A hybrid method which computes a solution for a restricted set of views is importance-driven radiosity presented in [smit92]. In this thesis we consider the radiosity problem. Surface properties are restricted to Lambertian reflection and light-sources are assumed to be Lambertian emitters." Lambertian" 1 means incoming light is equally likely to be scattered in any direction independent of the incident direction. So the amount of light reflected is independent of the outgoing direction of the light. The radiosity problem is modeled by the radiosity equation, a Fredholm integral equa-tion of the 2nd kind (e.g. [kaji86]). The kernel of the integral operator of the radiosity problem is positive, very smooth in a wide range of the domain on which it is defined; but has a certain number of discontinuities. < The accurate solution of the radiosity equation is a hard problem mainly because the evaluation of the kernel is very expensive. Currently known techniques are slow, generally inaccurate or may require large amounts of memory. The solution approach this thesis is concerned with is called wavelet radiosity [gort93b], [schr94], [chri94], [gort95], [schr96]. The idea of wavelet radiosity is essentially to employ a Galerkin approach where basis functions are chosen to be wavelets with vanishing moments. The discrete system generated by the Galerkin method is generally a blockwise sparse sys-: tern of linear equations, Mx = b. Using wavelet basis functions with the vanishing moment property results in numerous elements of the matrix M being very small. Wavelet radiosity follows the general method for the solution of integral equations presented by [beyl91] and exploits this property in two ways: • Small kernel coefficients are set to zero, with the remaining system being:sparse. This allows for a fast solution method. • A hierarchical approach allows not having to compute the small elements of M at all. This is essential since the computation of each element of M is costly. This thesis investigates two components of the wavelet radiosity method: 1. The choice of the wavelet basis. A wavelet radiosity algorithm is implemented in 2 a. very general framework. W e use the implementa t ion to analyze the behavior, of higher-order wavelet bases c o m p a r e d to the c o m m o n l y used H a a r basis, which consists of piecewise constant functions. T h e higher-order wavelets we use are multiwavelets :[alpe93] and the wavelets presented in [cohe93a]. T h e latter wavelets have a higher degree of smoothness t h a n multiwavelets. • •< '' 2. T h e appl icabi l i ty of K r y l o v methods as solvers for the sparse system result ing from the wavelet radiosity m e t h o d and their convergence properties as c o m p a r e d to the c o m m o n l y used P i c a r d i terat ion. W e find that higher-order wavelet basis functions tend to have a perceptual advantage over the, H a a r basis and reduce the required amount of memory. O n the other h a n d the same error in.the. L\ n o r m can be reached faster using a H a a r basis. , Our . experiments show that the. K r y l o v m e t h o d G M R E S (see e.g. [saad95]) is an easy-to-realize alternative to P i c a r d i terat ion which in our experiments always performed better t h a n .the P i c a r d i terat ion. R o u g h l y speaking we can reach the same accuracy as the P i c a r d i terat ion wi th at most hal f the n u m b e r of iterations. Fo l lowing this. in troduct ion the thesis presents a general formal framework and an overview of exist ing radiosity methods in C h a p t e r 2. In C h a p t e r 3 we present different types of wavelet bases and solut ion methods . In C h a p t e r 4 we describe and justify details of our implementa t ion of wavelet radiosity in order to proper ly specify the numer ica l experiments conducted . C h a p t e r 5 concludes w i th a s u m m a r y a n d evaluat ion of the obta ined results. 3 Chapter 2 The Radiosity Problem and Its Numerical Solution 2.1 The Radiosity Model T h e prob lem of g lobal i l l u m i n a t i o n considered here is how to obta in the radiance L(p,u>) (i.e. angle-dependent emit ted power of light per unit area for a fixed wavelength) leaving a point p of the surface S of a scene at a sol id angle LO G £1 (tt is the set of sol id angles). T h e scene is descr ibed by geometric a n d phys ica l properties . T h i s formulat ion is view-independent . B y c o m p u t i n g radiance for a fixed wavelength we assume that c o m p u t i n g color values corresponds to the c o m p u t a t i o n of radiance values for a sufficient n u m b e r of wavelength samples. W e assume that emi t ted radiance is independent of LO. T h i s type of reflection is cal led Lambertian reflection. Light-sources w i th this property are cal led Lambertian emitters. T h i s assumpt ion allows character iz ing the i l l u m i n a t i o n of the scene in terms of radiosity (i.e. the emit ted power of l ight per unit area). T h e radios i ty u : S —> R is defined as u(p) — JnL(p,Lo)dLo. W e restrict our considerations to interact ion between 4 surface points and exclude effects tak ing place in the m e d i a in which the light is transferred. Fur thermore , we wi l l not include translucent objects into our model . F o r a discussion of these effects, see [lang95] a n d [rush90]. T h e m a t h e m a t i c a l descr ipt ion of this physical m o d e l is given by the radiosity equation, an integral equat ion of the second k ind: u(p) = e(p) + JJ r(p, q)u(q)dS(q) (2.1) W e wil l also write the radiosity equat ion in the short -hand notat ion u = e + Tru (2.2) where T r is the corresponding integral operator. T h e quantit ies S, e a n d r are defined as follows: T h e set S C R 3 is the surface of the scene. T h e surface is composed of m u t u a l l y disjoint surfaces S^. S = |JSj . T h e funct ion e : S —> R describes the light-sources of the scene in terms of emit ted radiosity. T h e support of this funct ion is usual ly smal l c o m p a r e d to S since only few areas of the surface of a scene are light-sources. • '-T h e function r : S x S —> R is the non-negative radios i ty kernel given by r(p,q) = P(p)?~i(p,<7), where ri(2,q) = vs(p,q)r2(p,q) cos Z ( n x , q — p) cos l(ny, p — q) T2{p,q) |2 7T p-q \2 where v$ is the so cal led visibility function. T h e funct ion vs describes the geometry of the scene. It satisfies v${p, q) = 0 if the line segment from p to q intersects w i th 5 (which means that there are obstruct ions between p a n d q see F i g u r e 2.1), vs(p, q) — 1 otherwise. F i g u r e 2.1: V i s i b i l i t y funct ion. T h e gray line shows points y such that vs{x,y) = 1. T h e funct ion p : S —> [0,1)] is the reflectivity. It describes the physical properties of the scene, specifically, it describes the fract ion of radiosity reflected in a point of the surface. T h e vector nx 6 R 3 is the n o r m a l in p po in t ing to the thereby defined interior of S. Si l l ion and Puech [sill94] argue based on phys ica l grounds that \\TT \\L < 1, which guarantees a unique solut ion for the radiosity equat ion. A more rigorous analysis is given in [gwin87]. F o r sufficiently s m o o t h surfaces S the s ingular i ty in r for p = q can be removed by sett ing r(p,p) = 0. T h e kernel is generally discontinuous due to the vis ibi l i ty funct ion. Since e is typica l ly discontinuous the so lut ion u can be discont inuous as well. T h e c o m p u t a t i o n of the kernel r is expensive: Na ive ly the cost of c o m p u t i n g one value of Vs is l inear in the n u m b e r of surfaces. T h i s wi l l t u r n out to be the m a j o r challenge in solving the equat ion. T h e solut ion of equat ion (2.1) can be expressed in terms of the N e u m a n n series: o o (2.3) i=0 6 T h e components Tlre of the s u m have an intuit ive phys ica l interpretat ion: they rep-resent the contr ibut ion of the light-sources after i reflections on the surface. Important effects that cannot be mode led wi th the radiosity equat ion are specular highl ights and ideal mirrors (due to the assumpt ion of diffuse reflection and L a m b e r t i a n light-sources [cohe93b]). T h e radiosity equat ion has become a s tandard m o d e l in C o m p u t e r G r a p h i c s . E x h a u s -tive introduct ions can be found in [cohe93b] and [si!194]. N u m e r o u s methods have been suggested to solve the radios i ty equat ion. T h o s e m e t h -ods typica l ly have four objectives: be ing fast, accurate, requir ing little m e m o r y and be ing robust (i.e., require l itt le user interact ion) . However, speed a n d accuracy can have different interpretat ions in c o m p u t e r graphics t h a n in the c o m m o n unders tanding of numerica l anal -ysis. A m e t h o d which obtains a final result more slowly but produces useful intermediate results m a y be preferable to an overall faster m e t h o d . Such methods are called' progressive: A less accurate result which reproduces certain v isual effects, like smal l shadows, well might be more useful than a m e t h o d wi th an overall smaller error in a s imple error n o r m . In the fol lowing sections we wi l l describe principles of i m p o r t a n t determinist ic methods for the solut ion of the radios i ty equat ion. F o r an overview of stochastic methods we refer to [veac97]. Sect ion 2.2 is concerned wi th how methods discretize the radiosity funct ion a n d kernel. Sect ion 2.3 presents methods for the solut ion of the discret ized system. In section 2.4 we discuss how the discretized system itself, the form factors, can be computed . W e chose to present the three components of discret izat ion, solut ion m e t h o d and form factor c o m p u t a t i o n independently, since it allows us to point out more clearly various advantages a n d disadvantages of single components of a radios i ty m e t h o d . 7 2.2 Discretizing the Radiosity Equation Before publ icat ions of [gwin87] and especially [heck91], which in troduced G a l e r k i n methods ; for the so lut ion of the radiosity equat ion, radiosity methods were based purely on physical considerations of transfer of energy between surface patches governed by the pr inc iple of energy balance. T h e presentation here gives a un i form overview in the framework of the . no ta t ion in troduced in subsection 2.2.1. In the fol lowing the d o m a i n S wi l l be composed of p lanar surfaces Si (see [scha97] for an overview of methods for curved surfaces). : 2.2.1 Fundamental Problems T h e under ly ing pr inc iple of al l the procedures presented here is the G a l e r k i n method.; T h i s m e t h o d approximates the exact solut ion u of equat ion (2.1) by a funct ion u in a finite d imens iona l space V = span{</>0 • • • 4>n-i}., where fa . : S —>• R are l inearly independent functions. T h e s u p p o r t ^ of each fa, also cal led a patch, is entirely contained in a surface Si. T h e reflectivity p is assumed to have constant value pi on patch Ai (for the efficient treatment of non-constant reflectivities see [gers94]). W e define the project ion P onto V for v E L 2 ( S ) : funct ion v = Pv minimizes || v — v \\L2 in V. T h e B u b n o v - G a l e r k i n a p p r o x i m a t i o n approximates u by u € V such that the resid-ual of the integral equat ion projected into V is m i n i m i z e d . T h i s results in the fol lowing approx imat ion : n - l (2.4) where <v,w> = Jsv(p)w(p)dS(p) is the sca lar-product in L2(S). It can be shown that the 8 Defini t ion 2.2.1 The equation Pu = Pe + PTrPu is called B u b n o v - G a l e r k i n a p p r o x i m a t i o n of the integral equation u = e + Tru. W i t h u = Y^=o xi(t)i^ Pe — YTiZl h^i, a n d expansion of the projections, the B u b n o v -G a l e r k i n a p p r o x i m a t i o n results in a discrete system of equations: 71-1 • Xi = h + ^2fijXj i = 0 , 1 , . . . ,n - 1 (2.5) T h e coefficients fa = are cal led forrri factors where G = (gij)nxn is defined as follows: .9ij = ^ , l,<TT(j)3^i> = 1 if j r1(p,q)(j)i{p)(f)j(q)dAJ(q)dAt(p) (2.6) <oi: o ,> «pi,.0i> JJAiJJAj - - - - - - • W e wi l l write the system in m a t r i x notat ion as Mx — b w i th M = (mij)nxn = / — F , = {fij)nxn, k= (b0,... , 6 „ _ i ) T . W e cal l F the m a t r i x representation of PTrP. T h e form factors physical ly relate to the transfer of energy from patch Ai to patch Af. T h e i r c o m p u t a t i o n is what makes the radios i ty p r o b l e m so t ime-consuming . T h e approx i -m a t i o n of gij involves the c o m p u t a t i o n of samples of the v is ibi l i ty funct ion vs{p_,q) which is costly (as wi l l be discussed in Section 2.4). N u m b e r s of patches range f rom about 1000 for s imple scenes (e.g. [cohe86]) to 10 6 for complex scenes (e.g. [tell94]). A m a j o r objective of the methods presented is to compute an a p p r o x i m a t i o n M of M which requires as l itt le of the coefficients g^ as possible. 2.2.2 Basic Methods Classical Rad ios i ty Class ica l radiosity, also cal led full matrix radiosity, starts wi th a subdiv i s ion of the surface S into disjoint surface patches Ak. T h e operator Tr is approx imated by PTrP, hereby direct ly 9 real iz ing the G a l e r k i n approach described in subsection 2.2.1. T h e mesh has to be adapted to the solut ion such that shadow boundaries and h igh radiosity gradients can be adequately represented. T h e m e t h o d uses basis functions which are constant on the patches: {1 : p e A k (2-7) 0 : otherwise W i t h these basis functions we obta in for the form factors: In case of piecewise constant basis functions the form factors represent the fract ion of energy transfered f rom patch k to patch / . Class ica l radios i ty requires user interact ion for the generation of the mesh. T h e n u m -ber of form factors that have to be c o m p u t e d is quadrat i c in the number of patches. T h e c o m p u t e d solut ion has discontinuit ies due to the type of basis functions that we use. These discontinuities are smoothed out in a post-processing step that l inearly interpolates between the centers of the c o m p u t e d patches. T h e post-processing step is not intended to decrease the numerica l error but rather aims at m a k i n g the obta ined result v isual ly more pleasing. Patches can have arb i t rary shape and so they can be flexibly adjusted to the prob lem. T h i s m e t h o d was first presented for convex scenes (scenes wi thout obstruct ion) by [gora84]. It was then generalized to non-convex scenes by [nish85] and [cohe85]. Substructuring Sub structuring a ims at reduc ing the number of form factors to compute . F o r m factors are in i t ia l ly c o m p u t e d on a coarse mesh. A n approx imate solut ion is found on this coarse mesh. T h e n the solut ion is transfered to an in i t ia l fine mesh, which is subsequently i teratively 10 refined based on the gradient (e.g. c o m p u t e d by a centered difference) of the solut ion found so far. :;• T h e m e t h o d uses piecewise constant basis functions 0 O A . : 0 < k < n0 on the coarse level: {1 : P e A0tk 0 : otherwise O n the finer level there are ri\ = X ^ l o * Qk piecewise constant basis functions 4>i,(k,q) 0 < k < no, 0 < q < qk such that: 9/fc-l <t>o,k = ^2 01,(/c, 9) T h e refinement process is local a n d the in i t ia l mesh is sufficiently fine that qk = 1 for most k. T h e project ions P0 a n d Pi are defined as in equat ion (2.4) based on the basis functions <^ n,. and </>lr. T h e operator P\TrP0 is described by the form factors: h,(k,q),i = ^l\M',{k,q)\'ll II r(p,q)dAo^q)dAh(k>q)(p) J " T h e form factors of PQTTPQ can be c o m p u t e d based on A;i,(i l9). Oi-l k0iij = l/\A0>l\ ^ ^ i , ( » , g)jl^ l,(z , g)l T h e approx imate so lut ion ui is then obta ined by first so lv ing the discrete system of equations UQ = Pof + P0TrP0uo a n d then transferring this so lut ion onto the fine mesh U\ = Pif + PiTrP0u0. T h i s latter step is repeated while refining Pi unt i l the radiosity gradient decreases below a threshold value e. In total the second step is just one step of P i c a r d i terat ion for PITTPQ and as such improves the so lut ion u\ in compar i son to UQ. T h e physical in tu i t ion of this a l g o r i t h m is that UQ is a sufficiently good a p p r o x i m a t i o n of the average radiosity on the coarse mesh. T h i s coarse mesh however is not able to repro-duce sharp shadow edges. T o increase the structure of the radios i ty so lut ion wi th in a coarse 11 patch A0ti we transfer radiosity from the other coarse patches to the fine patches Ai^iqy Some fine patches might not receive radios i ty since they are obstructed by an object. T h e goal of this operat ion is to improve the shadow representation. S u b s t r u c t u r i n g was in troduced in [cohe86]. T h e number of form factors c o m p u t e d is st i l l quadrat ic , 0(nl -f-nnni), however, no is m u c h smaller t h a n in fu l l -matr ix radiosity. T h e process st i l l requires user interact ion to determine the coarse mesh a n d the in i t ia l fine mesh. T h o s e meshes have to be fine enough a n d adapted to the solut ion. If the the in i t ia l fine mesh is not proper ly adapted to the so lut ion then refinement takes place in areas of smal l radios i ty gradient. If the coarse mesh is not fine enough v isual artifacts might appear like the reflection of radiosity from shadows. T h e final step which transfers the solut ion onto a fine gr id can only redistr ibute radiosity w i th in a coarse patch but not a m o n g coarse patches. T h i s m e t h o d is an a t tempt to approx imate the discretized kernel PiTrP\ (represented by an n\ x n\ matr ix ) by a sparser kernel requir ing fewer form factors (represented by an n 0 x no a n d an ri\ x n 0 m a t r i x ) . T h e c o m p u t a t i o n of this a p p r o x i m a t i o n is dr iven by the approx imate solut ion. Higher-Order Bases G a l e r k i n methods which use higher-order basis functions have the potent ia l of decreasing the number of necessary form factors. T h e y provide a better a p p r o x i m a t i o n of the radiosity solut ion in areas where it is smooth . A smal l number of higher-order basis functions can approx imate such a funct ion well, whereas for piecewise constant basis funct ion strong re-finement would be necessary. F u r t h e r m o r e , w i th higher-order basis functions it is possible to use a p p r o x i m a t i o n spaces wi th certain smoothness properties. Such an a p p r o x i m a t i o n makes the post-processing step redundant . W h i l e higher-order methods wi th basis functions wi th wide support are useful in areas 12 of s m o o t h var ia t ion patches have to be refined carefully a r o u n d discontinuit ies . : Higher-order methods have been appl ied to the radiosity p r o b l e m by [heck91], [trou93] and [zatz93]. T h e m e t h o d presented by [zatz93] uses Legendre po lynomia l s as basis functions and is cal led Galerkin radiosity. 2.2.3 Hierarchical Methods T h i s section gives a brief mot iva t ion for hierarchical methods . T h e general idea of the methods is to find a hierarchical representation of the approx imate operator PTrP. T h e hierarchical representation contains form factors between hierarchies of patches. W h i l e there wi l l generally be significant light transfer, and consequently b ig form factors, between b ig patches which are close to each other, the magni tude of transfer decays when patches are smaller a n d are further apart . F u r t h e r m o r e we want to find a basis such that form factors between patches w i t h smal l var iat ion of the kernel are smal l themselves. W e can then set al l form factors below a certain threshold e to zero, thereby obta in ing a sparse representation M£. A part icu lar challenge is to obta in the sparse representation wi thout hav ing to compute al l the form factors beforehand, i.e. to decide which form factors would be below the threshold wi thout c o m p u t i n g t h e m first. F o r the sake of c lari ty we t e m p o r a r i l y switch to the abstract case of a one-dimens ional integral equat ion wi th kernel s : R x R —» R such that the fol lowing smoothness requirements are fulfilled: \s(x,y)\ < 7 - f 1 V >y>\ - \x-y\ ^ These requirements allow for one discont inui ty a long the d iagonal x = y. T h e follow-ing discussion follows [beyl91]. 13 H i e r a r c h i c a l basis funct ions W e define hierarchical basis functions w i th certain properties required for the procedure described below. F i r s t we define a hierarchy of finite d imens iona l spaces Vo C V\ C ... C Vm. F u r t h e r -more we choose spaces Wj such that ( © is the direct sum): wJ®vj = vJ+l T h e space Vj is defined by basis functions 4>j^'- Vj = span : 0 < k < fij}. E q u i v a l e n t l y w e choose basis functions "0j,fc ' 0 < k < rij for the spaces Wj. T h e , hierarchy of spaces implies a hierarchy of basis functions. T h i s h ierarchy is expressed by the dilation equation: $j,k = YALV 1 hk,l&j+l,l ^ ^j,k = J2ZV~L 9k,i(f>j+i,i W e represent functions Vj € VJ; VJ = X^Lo* c3,k<t>3,k '• 0 < j < m, by the .vector of coefficients c- = (cj^oKkKnj and functions Wj e WJ; WJ = YTk=o dj,k4>j,k '• 0 < j < m , by the vector of coefficients = ( d j ^ ) o < / c < n r If we start up wi th a funct ion vm represented by c m in the finest space Vm then we can successively decompose it in the form v = u> m _i + u>m-2 + • • • + v~>o + VQ. T h e c o m p u t a t i o n oi Cj : 0 < j < m a n d rf- : 0 < j < m f rom c m is cal led a wavelet decomposition. If the functions <f>jtk are orthogonal for a fixed j then the c o m p u t a t i o n of c m f rom given Cj '• 0 < j < m a n d d3, : 0 < j < m is possible. T h i s reverse operat ion is cal led a wavelet reconstruction. If we require that for fixed k on ly a constant number of the coefficients hkti a n d gkj be non-zero then b o t h computat ions can be conducted in a s traightforward fashion in 0 ( n m ) operat ions v i a the fast wavelet transform. 14 We can represent the wavelet decomposition by a non-singular square matrix Ds<r This allows expressing the decomposition as: \ do (2.11) \ 4 . - 1 / Analogously, a decomposition which includes Cj on all levels is denoted by D^m- The matrix DN>m is non-square. A reconstruction which adds up the functions represented by c-and d- and returns a representation c m is denoted by Djf m: to do D N,mQ.n D N.m ( C ^ t o tm-1 (2.12) \ dm-l J \ dm-l J The matrices Ds^Dg^, and D^m are never formed explicitly. They merely serve as a convenient notation for wavelet decomposition and reconstruction. Note that Ds,m transfers the representation crn of v G-V into a different basis while the representation obtained by TJjv.m is not a basis representation at all. For this reason Ds,m is called a standard decomposition. Dyv i m is called a non-standard decomposition. ' Analogously to equation (2.4) we define the projections Fj onto Vj and Qj onto Wj. If we 'assume orthonormality of the functions 0 ^ for fixed j, i.e. <<j>j,k, ^>j,k> = ^-,<(f)j,k, <fij,i> — 0 for k^l, then Pj+1 = Pj + Qj. It is possible to define sets of functions ct>j^ and ipj,k such that: J ipj>k(x)x^dx = OVfj, : 0 < fi< v 15 T h e funct ion tpj^ is then said to have v vanishing moments. W e wi l l give examples of such functions in Sect ion 3.1. T h e functions qt>jtk a n d ipjtk, which we have defined here, fit into the framework of or thogonal wavelets discussed, e.g., in [daub88] and [jawe94]. In this framework the functions tpjtk are cal led wavelets, a n d the functions (pj^ are cal led scaling functions. T h e pair of wavelets a n d scal ing functions is cal led multiresolution analysis ( M R A ) . W e ignore the more general framework of b iorthogonal wavelets (e.g. [jawe94]) since we do not app ly it in this thesis. Wavelet Compression of Integral Operators W i t h the above properties we can express Pm as: 771—1 Pm — Po + ^2 Qj j=0 A s a consequence we can rewrite PmTPm: 771—1 771— 1 771—1 771—1 PmTPm = P0TP0 + P0T(J2 Qj) + ($2 Q ^ T P O + C QJ)T(J2 Qi) (2-13) j 0 j 0 j 0 j=0 A c c o r d i n g to equat ion (2.13) another way to write PmTPmv in m a t r i x notat ion is dm = Ds)nFsDs,mCm where the m a t r i x Fs is defined as F, s and F(p°'p°} w i th entries of the form <T(p0>l, </>0,fc>, F^P°'Q^ w i th entries <Tipjih <po,k>, F^Q'P^ with entries <(f>o,k,Tipjj>, a n d F^'^ w i th entries <Tipj,i,tpi,k>- T h e m a t r i x Fs is cal led the standard representation of PmTPm. 16 i=m-7 QTQ | PTQ QTP QTQ PTQ QTP Figure 2.2: Non-standard representation FN of PmTPm. Another way to rewrite PmTPm (see Appendix A.2) is: m—1 m—1 m—1 P m T P m = P0TP0 + QjTQj- + QiTPi + 52 PiT®i (2-14) i=o j=o j=o With equation (2.14) we get the non-standard representation of PmTPm: dm = D^mFNDNtincm • The matrix FN is composed of blocks F^p°'Po^ = (<T<f>Qii, <l>o,k>)o<k<n0,o<i<no, F&'W = (<T^-i/,Vj)fc>)o<fc<no,0<Z<no. a n d F(Q>'<3>) = (<Tlpjth Ipj,k>)o<k<n0,0<l<n0 • The structure of Fjv is shown in Figure 2.2. Analyzing the sparsity of the matrices F, Fs and FN we obtain the following numbers of non-zero entries for the typical case of a dyadic hierarchy (nJ+i = 2rij): representation matrix non-zero entries full matrix F n2 standard Fs n 2 m - n l = n l ( l - 2 - ^ ) non-standard FN n2 11•m 17 In terms of sparsity the wavelet representation has not brought any advantage in compar i son to the full m a t r i x representation. However, if we use wavelets w i th vanish ing moments we see that entries <Ti/jjtk, <p> a n d <T(p, il>j,k> wi l l be smal l if ipj,k has its support in an area where the kernel is sufficiently smooth . T h e fol lowing result is presented in [beyl91]: Theorem 2.2.2 Assume an integral operator T with kernel s is given which satisfies equa-tion (2.9). Consider e > 0 and wavelets with local support. Then: • only O ( n m l o g ( n m ) ) of the entries of F$ are larger than e • only 0(n m) of the entries of are larger than e T h e in tu i t ion beh ind this result is that only those coefficients which are based on wavelets that are close to or overlap the discont inui ty a long the d iagonal x = y wi l l have significant value. T h e enumerat ion of such wavelets for the case of the H a a r M R A is shown in F i g u r e 2.3 a n d F i g u r e 2.4. F o r the s t a n d a r d representation we can compute a m a t r i x Fsc by sett ing elements of Fs which are smal ler than e to zero. T h i s procedure is cal led kernel compress ion. W e denote the operator represented by DgmFseDs>m w i th Tsc. F o r v £ Vm this procedure allows us to compute an a p p r o x i m a t i o n of PmTPmv w i th only 0 ( n m l o g ( n m ) ) entries in compar i son to 0(nm) entries of the full m a t r i x approach . If u is the solut ion of the integral equat ion Pmu = Pme + PmTPmu and u is so lut ion of Pmue = Pme.+ Tstue then it can be easily shown that || u — ue | | 2 = c „ m e , for sufficiently smal l e (e.g. [atki76j). T h e same holds for the non-s tandard representation, however, the operator here has only 0 ( n m ) entries. T o summar ize we have found approx imat ions to the m a t r i x F, Dg^Fs^Ds^ (and ^ j v .m-^Vt - D j v .m) , which allow c o m p u t i n g a m a t r i x vector produc t in 0 ( n m l o g ( n m ) ) (and 18 m— w=n 1 jt=n 12 3 B . a; teQ, 1.2,3 f l Jc=0 (c=2.3 £=4.5.6,7 tf j = 2 1=0 ¥ ¥ V V j = 2 1=1 Figure 2.3: S t a n d a r d representation, H a a r M R A : supports of h(x,y) = ipitk{x)tpj,i(y} for i,j,k,l such that the support of overlaps the discontinuity. f="—1 i -\ ^ \ ! \ =^0^ 3.4.5.6.7 F i g u r e 2.4: N o n - s t a n d a r d representation, H a a r M R A : supports oih(x,y) = xl)i^{x)il)j^{y),i = j for i,j,k,l such that the support of h overlaps the discontinuity. 19 0(nm)) operations. For this method to work we need hierarchical basis functions with the following properties: • Wavelets with vanishing moments. This property is the main requirement neces-sary for compression. It guarantees that coefficients in areas where the kernel is smooth are small. • Local support. Local support is necessary to guarantee a large number of wavelets with support in areas where the kernel is smooth, i.e. away from the discontinuity. • Dilation equation. The basis functions must fulfill equation (2.10) with a finite number of non-zero coefficients hkyi and g^. This property is necessary to be able to compute the decomposition Ds,m and in 0(nrn) operations. Technically the application of the method presented here to the radiosity equation is straightforward with two-dimensional scaling functions and wavelets (see Appendix A . l ) . Schroder [schr94] shows that the same complexity bounds hold by extending Thm. 2.2.2 to the integral operator that appears in the radiosity equation. Computation of the Compressed Operator In the framework of hierarchical basis functions the coefficients <Tipjti, ipjtk>, <Tipjti,(f)j!k>, <T(pj:i, ipjtk> and <T(j>jti, (f>j:k> for k, I : 0 < k,l < rij can be computed recursively using the 20 dilation equation: rij+i — l n j + i - l <Ti/jj,i,ipj,k> = 52 52 9i,Wk,k<T(f)j+ij,^j+i'k> t O k=0 <Tipjh(phk> = 52 52 dijK^T^j+ijAj+i;^ 1=0 fc=0 <T4>j,i,ipj,k> = 52 52 \i9k,k<T<t>j+\j.Aj+i,k> 1=0 k=0 T i j + i - l n j + i - i ' <T<t>j,i,<l>j,k> = 52 hi,ihk,k<T(t>J+ij,(t)j+i;k> (2-15) r=o fc=o The algorithm outlined above in the presented form is not useful for radiosity prob-lems. It reduces the cost for solving the discrete system to 0 ( n m ) , but the cost of computing Fse or F^c still requires 0(nm) operations. Furthermore, in radiosity problems the time spent on the solution of the system is usually negligible compared to the time it takes to compute the form factors. Therefore we have to find a method which allows us to compute all form factors above the threshold e without having to compute those which are smaller. The ideas proposed in [hanr91] (for piecewise constant basis functions) and [gort93b] (for higher-order basis functions) are based on a method known from n-body problems-pre-sented in [appe85]. The method computes form factors starting at the coarsest level. If a form factor between patches p and q is above the threshold e then we compute the form factors between all children of p and all children of q. This process is repeated recursively for the children. If a form factor is below the threshold e then we stop refining. This procedure assumes that if a form factor between patches p and q is small then this also will hold for form factors between subpatches of p and q. In addition it is too expensive to compute the form factors before comparing them. On coarse levels this would require a high number of kernel samples for the quadrature used for their approximation. Computing all these samples is exactly what the hierarchical 21 method tries to prevent. To avoid the full computation of form factors a procedure called' an oracle is applied which estimates form factors based on a small number of samples.or geometrical properties. These estimates are discussed in Section 2.4.2. Refinement, in the recursive procedure is then based on the estimate rather than a computation of the form factor. C lus te r ing Hierarchical methods require computation of the form factors of the operator Pr,TrPo, i^ e. the form factors between basis functions on the coarsest level. Supports of basis functions on the coarsest level are located within the surfaces that the scene is composed of. If s is the number of surfaces in a scene the complexity for the computation of the form factors on the coarsest level is 0 ( s 2 ) . Smits et. al. [smit94] present a clustering algorithm which reduces the cost to O(s). They use an octree to hierarchically cluster surfaces. The energy transfer between these clusters is estimated. Pairs of clusters with energy transfer above a certain threshold are hierarchically refined until the level of a single surface is reached. On that level the algorithm continues with the hierarchical method outlined in 2.2.3. The clustering technique can be applied in a framework with piecewise constant scaling functions since here form factors relate to the energy transfer between patches. Summary Hierarchical methods with hierarchical basis functions as required in subsection 2.2.3 find an approximation to the operator P n m T P n m of arbitrary accuracy with 0 ( n m ) form factors. The-accuracy is proportional to the threshold e. In an implementation these rigorous bounds do no longer hold due to the estimates used for a top-down computation of the form factors and the termination of recursion at a certain level. The meshing process is fully automatic up to 22 the selection of the threshold parameter e and the number of levels. These two parameters provide the user with an easy way to trade speed for accuracy. The realization with standard representation has a theoretical disadvantage because it requires O(nmlog(nm)) form factors. However, it allows computations without the normaliza-tion and denormalization steps (D^lm, -Djv>m) since Ds,m is applied to a basis representation of u and returns such a representation. The transformation into this representation is applied once in the beginning of the algorithm, the back-transformation applied once at the end. j Due to the hierarchical structure of the procedure, form factors cannot be computed independently. This implies that all form factors have to be stored until the final radiosity solution has been found. This typically results in high memory requirements. . A hierarchical method with piecewise constant scaling functions was first presented in [hahr91] under the name hierarchical radiosity. The general framework of wavelet methods for integral ;equations was presented by [beyl91]. Gortler et. al. [gort93b] applied that method to the radiosity problem and called it wavelet radiosity. While [gort93b] use the non-standard representation [chri94] and [stol96] describe a wavelet method using the standard representation. This method needs only one wavelet decomposition step in the beginning and one final wavelet reconstruction step. The entire solution procedure takes, place with a representation of the unknown in the wavelet basis, not in a basis of scaling functions as with the non-standard decomposition. 2.2.4 Discontinuity Meshing Discontinuity meshing constructs an appropriate basis based on the geometry of the scene. One computes potential locations of shadow boundaries by intersecting planes which touch edges of surfaces and corner vertices of light-sources with other surfaces in the scene. The lines generated by the intersection are potential locations of discontinuities. This process is 23 iterated using the newly generated lines in place of the light-sources. The areas generated by these intersections are used as supports for piecewise constant basis functions in the Galerkin discretization. Heckbert [heck92] and Lischinski et. al. [Iisc92] presented the first algorithms for discontinuity meshing. Discontinuity meshing can accurately resolve shadow boundaries. Following [lisc93] the number of basis functions computed in discontinuity meshing is 0(l2s2) if / is the number of light-sources and s is the number of surfaces in the scene. In [lisc93] discontinuity meshing is regarded as an expensive process. 2.2.5 Hybrid Methods Lischinski et. al. [Iisc93] combine hierarchical methods with piecewise constant basis func-tions and discontinuity meshing. They build the hierarchy of patches by subdividing along potential discontinuities. This process helps to terminate refinement on a coarse level since the'number of basis functions overlapping discontinuities is efficiently reduced. Bouatouch et. al. [boua95] experiment with a basis derived from multiwavelets which allows them to incorporate discontinuity information. Assume a discontinuity dividing the support Ajtk of a basis function fij^ into sets k and A? k . Then (frj^ is replaced by two functions </>])A. and ' <p\k which are obtained by restricting the support of 4>j,k to A l - k and A 2 k. This newly established basis is not orthogonal anymore resulting in more complicated expressions for wavelet decomposition and reconstruction. Bekaert [beka96] presents higher-order hierarchical basis functions with non-rectangular support which appear to be useful for combining discontinuity meshing and hierarchical methods. 24 2.3 Solution Methods W e do not consider direct methods in this section since in the context of the radiosity equat ion they are clearly inferior to iterative methods . T h e advantage of iterative methods here stems less f rom the desire of explo i t ing sparsity than from the strong diagonal dominance of the discretized system, which implies fast convergence of s imple iterative methods. Moreover , the elements of the m a t r i x relate to a discret izat ion of s m o o t h quantit ies, rais ing the prospect of efficient mult i leve l treatment . 2.3.1 Properties Cons ider the linear system of equations (2.5) wri t ten as Mx = b, where M — I — RG w i th R = d iag (p 0 , • • • , Pn-i) a n d G as in troduced in equat ion (2.6). A s ment ioned in 2.2.1 the size of M can become too large to fit into a computer's p r i m a r y memory. M is blockwise sparse if the scene has obstruct ions such that entire patches are not visible f rom each other. Spars i ty is significant if the scene is complex , i.e. if it is composed of different isolated rooms. M is generally not symmetr ic . However, it can be easily made symmetr i c by m u l t i p l i c a t i o n wi th D A R ~ l , where D A = diag(<0n, 0o>) • • • , <4>n-ii 4>n-i>)-A l s o , by definit ion two points in a plane are not visible f rom each other, so gkk = 0 and hence rrikk = l,k := 0 < k < n. Now we restrict our considerations to the case of piecewise constant basis functions. Here the identitv du(x) = ^-^^-MdA(y) (ui is the sol id angle) yields Y\,gki = 1 (e.g. \\x-y\\2 [cohe93b]). Since 0 < pk < 1, M is s tr ict ly row diagonal ly d o m i n a n t . F u r t h e r m o r e [neum95] a n d [niev97] show that DAR~lM is not only symmetr i c but also posit ive definite (SPD). T h i s property can be used to solve the S P D system DAR~lMx = DAR~lb instead. However, this m e t h o d does not allow for very smal l or zero pi, which is a significant drawback since such 25 situations can easily arise in computer graphics when objects are dark or black. Summary The system M arising from a Galerkin discretization of the radiosity equation has the fol-lowing properties: 1. moderate sparsity; 2. it can be made symmetric by multiplication with the diagonal matrix D^R~l if pk > 0, k : 0 < k < n; 3- rrikk = : 0 < k < n. Furthermore, for a discretization with piecewise constant basis functions the following holds: . 1. 'row diagonal dominance: Vz : \rrijj\ > \mij\ 2. rriij < 0,i / j 3. DAR-lM is SPD (if pk > 0,k : 0 < k < n) 2.3.2 Gauss-Seidel-like Methods In the following discussion we will use the terms update and sweep. The term update refers to the change of a single variable %i while a sweep means an update of all components of x. Jacobi Iteration One sweep of Jacobi iteration can be written as: x{k+l) = diag(M)" 16 + {I - d i ag (M)" 1 M)x ( f c ) 26 In case of the radiosity equat ion d i a g ( M ) = / . Reca l l ing F = RG we obta in the s impler form: T h i s form is also known as one step of Picard iteration (e.g. [hack95]). P i c a r d i terat ion s i m p l y approximates the N e u m a n n series from equation (2.3). So each sweep corresponds to one add i t iona l level of reflection. If we choose xj°) = b then Fkb is the contr ibut ion of the l ight-source after k reflections. A l s o note f rom equat ion (2.16) that one step of the P i c a r d i terat ion is equivalent to the add i t i on of the current res idual = b — Mx^ result ing in ^(fc+l) = x{k) + r(k)_ G a u s s - S e i d e l I t e r a t i o n L e t / count single updates rather t h a n sweeps. T h e n Gauss-Se ide l i terat ion updates the unknowns X{ in the following way: T h e index i(l) specifies which variable is u p d a t e d in the l-th step. T h e usual order is: (O'ew = (0,1, 2 , . . . , n — 1, 0 ,1 , 2 , . . . , n — 1 , ...). T o incorporate the res idual expl ic i t ly equat ion (2.17) can be easily rewrit ten y ie ld ing the fol lowing a lgor i thm: A l g o r i t h m 2.3 .1 G a u s s - S e i d e l U p d a t e ( M , 6, i,x,r) i:= next i in turn x (k+i) = b + F?.(k) (2.16) (2.17) V j : j / i d o : 27 ri :.= 0 return (i, x, r) If we choose xj0^ = 0 as the initial guess we can compute the initial residual: r^ = b. This gives us all the initial parameters we need for the algorithm. For one update of the approximate solution and for one update of the residual only one column of M is necessary. In the radiosity context equation (2.17) has a simple physical interpretation: in one Gauss-Seidel update light from all patches j is gathered and accumulates in patch i. The residual r^ can be written as r^ = M(x — x^) where x ^ ' is the current approximate solution at any give stage . The quantity x — is the energy surplus in the scene (we show later that this quantity is always positive for piecewise constant basis functions). For that reason the residual is also imprecisely called unshot energy in the graphics literature. . Gauss-Seidel in the graphics literature is also often referred to as gathering. Southwell Iteration The Southwell iteration ([sout46]) is a modification of Alg. 2.3.1. Here the index i of the next variable to update is not just chosen "in turn" but it is chosen such that jr-,| is maximized. Convergence of this method for the radiosity problem is shown in [gort93a]. Progressive Refinement Iteration Progressive refinement (PR) is a modification of the Southwell iteration. Here a certain number of Southwell updates is performed. After this the residual is added to the solution x^ obtained so far and the sum is returned as the approximate result x = xS^ + r^h\ This final step corresponds to exactly one complete sweep of the Jacobi iteration based on xjk\ For this sweep no new entries of M have to be computed. 28 W e rewrite the a lgor i thm m a i n t a i n i n g the invariant x^ = x^ + (x^ is the i terat ion variable in P R ) . T h i s form of the P R a l g o r i t h m gives it an immedia te physical interpretat ion. T h e approx imat ion x^ is more accurate t h a n x^k\ W e can use x^ for d isplay purposes after each update . Algorithm 2.3.2 P r o g r e s s i v e R e f i n e m e n t U p d a t e ( M , b, x, r) choose i such that ri is maximized • Vy / / d o : Tj . Tj TTljiTi Xi J . Xj TThJ7T'1 r e t u r n (i,x, r) Note that rj^ is not the residual regarding x^ but rather the residual of xjk^ from the Southwel l i terat ion. So as in i t ia l guess x^ we have to choose x^ := xj0^ + A0} — b. So in contrast to the Gauss-Se ide l m e t h o d the res idual of the approx imate so lut ion is not available in this m e t h o d . T h e phys ica l interpretat ion is in some sense exactly opposed to the one of the GaUss-Seidel m e t h o d . In P R unshot radiosity f rom one patch i is d i s tr ibuted to all other patches j (shooting). C h o o s i n g the m a x i m u m residual (sorting), which is in i t ia l ly concentrated in the light-sources, guarantees that significant light transports are conducted early in the i terat ion. P R was in troduced based on physical arguments in [coheSS] 1 . [gort94] showed the above re lat ion to Gauss -Se ide l i teration. Hhe original version of progressive refinement chooses an i such that |>li||r;| is maximized 29 method reference selection of i M P R [cohe88] DAL I = 1 Shao [shao93] DAL I + F = I + F Gortler [gort94] DAL £ j Fi « I + F Feda [feda92] [xu94] (I - F)DAL *• Pavg X u [xu94] (I - F)DAr_ I + F + ZM** = I + F + -£M9_p A Pavg Table 2.1: Overrelaxation methods. The index i is selected such that the absolute value of the i-th component of the given matrix is maximized. O v e r r e l a x a t i o n M e t h o d s A more general form of iterative methods is given by Alg. 2.3.3. A l g o r i t h m 2 .3 .3 U p d a t e ( M , b, i, x, r) Choose next i Choose A xt := x{ + A Vy' d o : Tj := Tj — TUjiA r e t u r n (i, x, L) Special cases of this algorithm are the Gauss-Seidel method and progressive refine-ment, with A = ri. Choosing A = uiri, 0 < to < 2 yields the SOR method (i.e. [stoe80]). The SOR method can be shown to converge for diagonally dominant systems. Ideally we would like to choose2 A = eiM~LL- Since M _ 1 is not known it has to be approximated. Several such approximations M have been suggested (see Table 2.1). They are based on approximation of M~l by the leading components of the Neumann series. 2ei is the i-th vector of unity 30 The quantities in Table 2.1 are defined as follows: pavg = YliPilni F = T\\R1DA (1 is a matrix with all entries equal to 1) and Fi = (fiki)nxn' Pkfki • k ^ I and (k = i or I = i) fikl = \ 0 : otherwise Note that with all the above methods A can be computed in O(n) operations. For one update one row and one column of the matrix F are necessary. Physically the computation of A using more components of the Neumann series cor-responds to not only accounting for unshot radiosity on a patch A i but also for the unshot radiosity reflected back to patch A i from other patches of the environment. [xu94] developed the mathematical framework for the methods presented above, which were originally derived based on physical arguments. For systems with piecewise constant basis functions they give convergence proofs for a modified version of the Feda method and the X u method. Comparison Let n be the number of unknown variables. Then all algorithms presented need O(n) time to compute one update. Southwell iteration, progressive refinement and the presented overre-laxation methods are more expensive per update since we have to find the maximum among all residual values. However, in the context of the radiosity equation the cost of computing the relaxation updates is in fact negligible compared to the cost of computing the entries of M. In the methods presented above the entries of M can be computed and stored while they are needed for an update. One hopes that such a procedure allows computation of rough ap-proximate solutions without having to compute all entries of the matrix M. Furthermore, we will never have to store more than individual columns and rows of M reducing the memory 31 complexity to O(n). This memory saving procedure, however, comes at the cost of having to recompute certain columns and rows. We assume the case of piecewise constant basis functions with properties outlined in subsection 2.3.1. For the light-sources we have rf^ = 6, > 0 ,Vz. Since < 0 for i / j we see that r\• > 0 in Gauss-Seidel and progressive refinement updates (the same was shown for the Gortler and the Shao overrelaxation methods). Thus, the iterates approach the solution monotonically, i.e. just by positive steps. Therefore choosing the i with the greatest riy like in progressive refinement, yields the fastest error reduction during the initial updates, as long as no variable is updated a second time. Based on physical considerations [cohe88] compute an estimate for the total amount of light that still has to be distributed. This amount, called ambient term, is distributed over all patches weighted by area, resulting in a contribution x a m b i e n t . The ambient term is not part of the iteration. It is added to the current approximate solution before it is output. This results in a progressive development of the output solution as required in Section 2.1. The ambient term was motivation for the development of the Feda overrelaxation method. The figures in [gort94] show that Southwell iteration initially has an advantage over Gauss-Seidel. The error obtained after n updates of the Southwell iteration is less than a third of the error obtained after n updates of Gauss-Seidel iteration; In the experiments presented by [gort94] a visual convergence is obtained after less than 6n Gauss-Seidel updates. Cohen et. al. [cohe88] obtain visual convergence with progressive refinement and an ambient term after | n updates. X u et. al. [xu94] compare progressive refinement and overrelaxation methods. They obtain a visually good result with the X u method after \n updates. Convergence of high accuracy of the X u method is obtained after about 2n updates, which is \ of the number of updates needed by the P R method to obtain the same accuracy. The other overrelaxation 32 methods perform inbetween these two extremes. 2.3.3 Other Methods Krylov Subspace Methods Krylov subspace methods often converge faster than the Gauss-Seidel method. A popular member of that family is the conjugate gradient (CG) method ([hest52], [saad95], [shew94]). It converges if the system M is SPD. A simple modification, C G N R , applies the C G method to the SPD system MMT and solves for y. The solution of the original system can then be computed as x = MTy. C G N R converges slower than the C G method but it does not require the matrix M to be SPD. Another method which does not require an SPD matrix is G M R E S ([saad95]). Generally, the convergence of any Krylov method cannot be better than the convergence of C G if the system is SPD. It requires to store the entire sequence of computed iterates. However, since iterative methods for radiosity problems converge usually in less than 10 iterations this is no significant drawback here. While the methods based on Gauss-Seidel iteration require information of individual entries of the matrix the only information required for C G and G M R E S is how to compute a matrix-vector multiplication. C G N R additionally requires the ability to compute a multiplication with the transposed matrix. [bara95], [niev97] apply the C G method to a discretization of the radiosity problem with piecewise constant basis functions. They solve the SPD system: DAR~lMx = DAR~lb Nievergelt [niev97] reports convergence of the conjugate gradient method in-10 times fewer iterations than Gauss-Seidel. [bara95] report slower convergence for scenes with low average reflectivity (pavg = 0.24). [will97a] apply the C G method to the non-symmetric system M. 33 They observe that the C G method still converges. They compare CPU-time and find that the C G method works as fast as Gauss-Seidel for scenes with pavg > 0.5, however slower for scenes with smaller reflectivity. Multilevel-like Methods A method which works with a sequence of approximations F^^x) to F is introduced in [hanr91] and further developed by [stol96]: Algorithm 2.3.4 while not converged compute F e ( / c ) (x ( f c ) ) by refinement of Fe(-k~^ (x}k~^) := solution of (T- F e<*>(xW))£( f c + 1> = 6 k:=k+l We require F^(x)x to converge to Fx with decreasing e > 0. The computation of a matrix F^ (x) which satisfies this requirement can be less expensive than the computation of an approximation to F itself. This method might need more iterations to converge, however, these additional costs are negligible compared to the potential savings of having to compute a smaller number of coefficients in total. 2.3.4 Comparison In the following paragraph we will compare the solution methods in the context of the discretization methods from section 2.2. Gauss-Seidel-like methods are suitable for standard Galerkin discretizations. To be efficient they require access to single rows and columns of the system to discretize. Krylov methods on the other hand are not well suited for this type of discretization since each 34 iteration involves a full matrix vector multiplication. This would require to compute and store all 0 ( n 2 ) entries of the matrix before starting the solution process. Such a procedure is only feasible for small systems and contradicts the requirement of being able to produce intermediate results during computation of the solution. Additionally reports on the behavior of the C G method do not indicate an advantage over Gauss-Seidel-like methods. We do not know about results obtained with other Krylov methods. Looking at the discretization obtained by the hierarchical methods we come to exactly the opposite conclusion. Here we do not have access to single elements of the matrix. The only Gauss-Seidel-like method applicable seems to be Picard iteration. As well all the Krylov methods are applicable. Storing the complete system with O(n) entries is inevitable in hierarchical methods, so this is no particular drawback of Krylov methods. The multilevel-like method is applicable to the hierarchical discretization since the computation of is a very natural process here, which can be realized by successive refinement of the hierarchy computed so far. 2.4 Computing the Discrete System The discussion so far has dealt with ways to discretize the radiosity equation and solving the discretized system. As shown in section 2.2 computing the system itself is generally the most costly part of the entire solution process. The entries of the linear system of equation are based on the quadruple integrals of the form factors. Due to speed considerations numerical integration will always be restricted to a small number of samples per dimension, typically nsamp < 4. Considerable effort has been spent in accelerating the computation of form factors. We give a brief overview of methods used in computer graphics. A more in depth discussion can be found in [cohe93b]. We leave the discussion of the handling of singularities 35 in the integrand to [zatz93] and [schr93]. 2.4.1 Form Factor Approximation Piecewise Constant Basis Functions W e discuss the most c o m m o n l y used techniques here. O t h e r suggested methods are a trans-format ion of the surface integral in a contour integral (contour integral method), analyt ic c o m p u t a t i o n for unocc luded scenes by a p p r o x i m a t i o n wi th s tandard geometries, a n d M o n t e -C a r l o integration. F o r piecewise constant basis functions we have to compute the integral: cos L(nv, q — p) cos L(nq,p — q) '2 1 ff ff ( N C O s Z ( n P ^ - p ) c o s Z ( n 9 , - 4 ) 9ki = J-T-7 / / / / vs(p,q) ~2 dA^qjdA^p) \A-k\ JjAk JJA, 7T|| p - q II 1 ff ff cos L(nx,q - p) " " vs(p,q) ^ X~ PJdu(q)dAk(p) If we choose p to be the center of patch Ak then a s imple a p p r o x i m a t i o n for gki is: ff cos L(nx, q — p) 9ki = JJnvs(P'Q) 1 * " -dcu(q) (2.18) If patches Ak a n d Ai are smal l , far apart a n d there are no occluders close to Ak then this a p p r o x i m a t i o n is justif ied. St i l l , one double integral invo lv ing the vis ibi l i ty funct ion vs has to be c o m p u t e d . T h i s can be done by s imple numer ica l integrat ion or by one of the methods described below. Hemicube T h i s technique subdiv ides il into a fixed number nherni of subsets fij 0 < i < nhemi. A s s u m i n g that the part of pa tch At subtended by Qi is complete ly visible or invisible f r o m p then the integral of equat ion (2.18) is independent of the locat ion of At itself. If the fij are sufficiently 36 smal l then this as sumpt ion approx imate ly holds. T h e m e t h o d computes a complete row of form factors fkh 0 < I < n. W i t h each Qi we store a coefficient 9hemi,i, a n index kx of a patch and a distance di. Ini t ia l ly di is set to oo. T h e a lgor i thm is shown below: A l g o r i t h m 2.4.1 ' f o r all patches A\: f o r all Qi which intersect the solid angle generated by A / / i f distance(p, A{) < di di-.=distance(p, Ai) ki'.=l cos Z ( n x , q - p) 9hemi,i : « / / dco(q) 9hemi,i is computed based on just one sample oj the kernel. f o r all I: 9ki: =0 f o r all 0 < i < nhemi: 9kki • ~ 9kki ~t~ 9hemi,i T h i s m e t h o d has obvious disadvantages. F a r away or smal l patches might not be proper ly represented by the discret izat ion of Q and result in al iasing in the obta ined solu-t ion. O n the other h a n d the m e t h o d can use exist ing Z-buffer hardware. T h e m e t h o d was in troduced by [cohe85]. T h e name hemicube is based on the fact that Qi are chosen such that they correspond to rectangles on the surface of a cube wi th center in p. 37 Higher-Order basis functions and Hierarchical Methods T h e r e are two major differences between the quadrature for piecewise constant basis functions a n d higher-order basis functions and hierarchical methods . F i r s t , the patches found in higher-order methods a n d on coarse levels of hierarchical methods are potent ia l ly bigger t h a n those found wi th piecewise constant basis functions. Secondly, s imple approx imat ions assuming constant radiosity over a patch would clearly defeat the purpose of higher-order methods . W i l l m o t t and Heckbert [will97b] use visibility-in-quadrature a n d fractional visibility to compute the form factors. V i s i b i l i t y - i n - q u a d r a t u r e approximates the form factor integral by numerica l quadrature . Var ious quadrature rules are discussed in [gers95]. These methods are expensive since each sample requires evaluat ion of the v is ibi l i ty funct ion vs- F r a c t i o n a l v is ibi l i ty is the m e t h o d in troduced in [hanr91]. B o t h patches Ak a n d At involved in the c o m p u t a t i o n of the form factor are subd iv ided into a 4 x 4 gr id . Next , for each gr id point of patch Ak the v is ibi l i ty funct ion wi th only one gr id point of Ai is computed . [hanr91] use a fixed scheme to al locate gr id points on Ak to a gr id point on Ai. F r o m these v is ibi l i ty samples we obta in an estimate a of the fract ion of Ak that is visible from A\. T h e form factor is then c o m p u t e d wi thout v is ibi l i ty funct ion (unoccluded formfactor) and afterwards m u l t i p l i e d by a. F r a c t i o n a l v is ib i l i ty reduces the n u m b e r of samples required from nAsamp to n2 samp' In hierarchical methods we can compute form factors on coarse levels accurate ly by us ing the d i lat ion equat ion a n d base the c o m p u t a t i o n on the form factors c o m p u t e d for the next finer level as given by equat ion (2.15). 38 2.4.2 Form Factor Estimation T o be able to realize the top-down approach taken by hierarchical methods we need a funct ion that can estimate the magni tude or an upper bound , of a form factor involv ing wavelets wi thout hav ing to expensively compute it. T h e i r c o m p u t a t i o n on coarse levels is prohib i t ive ly expensive since they would require a large number of kernel samples to be able to represent the kernel accurately. F o r piecewise constant basis functions [hanr91] give an upper b o u n d for the unocc luded form factor (this corresponds to the form factor involv ing only scal ing functions) . V i s i b i l i t y is incorporated using fract ional v is ibi l i ty y ie ld ing an estimate //c/ for <T(pjti, (j)jtk>- T o est imate the value of form factors invo lv ing the H a a r wavelets, fki is m u l t i p l i e d by a parameter £ > 1 if A , ^ is only part ia l ly visible from ' Ajj. ' In the general case of higher-order basis functions there is no m e t h o d known to b o u n d the unocc luded form factor. Therefore we have to take a different approach . L e t the n u m b e r of vanish ing moments of the considered wavelets be v. T h e n we know that a form factor is smal l if the kernel behaves s imi lar to a p o l y n o m i a l of degree v — 1. [gort93b] suggest to estimate the form factor by es t imat ing the deviat ion of the kernel f rom a p o l y n o m i a l of degree v — 1. T h e kernel is sampled and an interpolat ing p o l y n o m i a l is constructed based on these samples. T h e dev ia t ion is c o m p u t e d by c o m p a r i n g the p o l y n o m i a l to samples c o m p u t e d inbetween the previously c o m p u t e d samples. T h i s gives a qual i tat ive measure for the magni tude of the form factor. 2.4.3 Comparison F o r m factor c o m p u t a t i o n and es t imat ion for higher-order basis functions are more expen-sive t h a n for piecewise constant basis functions. T h e cheaper m e t h o d of fract ional v is ib i l i ty 39 produces inaccurate results in this case. Q u a d r a t u r e formulae need higher accuracy, conse-quent ly a higher number of v is ibi l i ty samples. Fur thermore , hardware-supported^ methods like the hemicube m e t h o d are not available here. 2.5 Summary T h e m e t h o d most c o m m o n l y used in computer graphics is progressive radiosity. T h i s m e t h o d combines a discret izat ion wi th piecewise constant basis functions, subs tructur ing and pro-gressive refinement as the solut ion m e t h o d . F o r m factors are typ ica l ly c o m p u t e d using, the hemicube approach . T h e hemicube m e t h o d can compute complete rows or co lumns of form factors as needed by progressive refinement. N u m e r o u s ways of paral le l implementat ions on different paral le l architectures have been presented ([reck90], [chen89], [bric95], [baum90], [guit91], [garc97]). These methods exploit that form factors can be c o m p u t e d independently. Progressive radiosity requires only O(n) m e m o r y (n is the number of patches), but requires 0(n2) execution t ime. T h e hierarchical m e t h o d wi th piecewise constant basis functions is discussed frequently in the l i terature. It requires O(n) m e m o r y a n d O(n) execution t ime. However, the constant factors in these t ime a n d m e m o r y estimates are m u c h bigger t h a n those for progressive radiosity. A l l considered methods t ry to reduce the system of equations (2.5) such that it can be c o m p u t e d a n d solved efficiently. O n the one h a n d this is done by explo i t ing smoothness properties of the solut ion. D i scont inu i ty mesh ing and higher-order G a l e r k i n methods choose a basis which approximates the so lut ion well. O n the other h a n d hierarchical methods a n d subs tructur ing use smoothness properties of the kernel. T h e u l t imate goal is to find a sparse representation fully automat ica l ly wi thout user interact ion by an adapt ive process. T h i s process can be based purely on the kernel itself, or in an iterative process it can be 40 based on the approx imate solut ion found so far. Hierarch ica l methods reduce the number of parameters to choose to the compress ion parameter e and the number of levels. Efforts have been made to integrate concepts from discontinuity meshing wi th higher-order hierarchical methods . However, no complete analysis of such a system has been pre-sented. C l u s t e r i n g has been appl ied to hierarchical methods wi th piecewise constant basis functions. W e do not know of reports on the appl ica t ion to hierarchical methods wi th higher-order basis functions. F o r piecewise constant basis functions numerous provably convergent solut ion meth-ods exist. These solut ion methods are also appl ied to systems based on higher-order basis functions and the compressed systems aris ing in hierarchical methods . However, we do not know of any careful m a t h e m a t i c a l analysis of their convergence. 41 Chapter 3 Exploring Wavelet Radiosity In this chapter we give background and motivation for the numerical experiments presented in Chapter 4. We will be concerned with the behavior of different wavelet bases in wavelet radiosity and the analysis of iterative methods for solution of the system of equations arising in that method. Section 3.1 presents and compares different wavelet bases, Section 3.2 motivates the application of G M R E S to the solution of our system. Finally Section 3.3 describes and justifies the algorithms used in our experiments. 3.1 Wavelet Bases for Radiosity We introduced the notion of multiresolution analysis (MRA) with scaling functions and wavelets in Section 2.2.3. M R A with various properties exist. We list some criteria that typically distinguish them: • Bounded or infinite support of scaling functions and wavelets. • Number of vanishing moments of the wavelets. 42 • Properties of the space spanned by the scaling functions. Normally M R A are defined on the real line. With wavelets on the interval we can span spaces defined on a bounded interval (e.g. [cohe93a], [alpe93]). A special type of wavelets on the interval are wavelets which satisfy boundary conditions on the boundaries of the interval ([chia97]). • Smoothness of the spaces spanned by the scaling functions. We will be concerned with three types of M R A : the Haar M R A , the multiwavelets, and an M R A introduced in [cohe93a] which we will denote as C D V M R A . A l l three types of M R A are orthonormal (VA; ^ I : <(j)j,k, <Pj,i> = 0, <<f>j,k, 4>j,k> = 1, Vi , j , k ^ I : <ipitk, "4>j,i> — ®,<'4>i,k,'lr,i,k> = !)• We do not consider flatlets here (flatlets [gort93b] are based on the Haar wavelet and have a higher number of vanishing moments). ' In Figure 3.1 we show the hierarchical structure of the three types of M R A . This allows us to use a simplified notation for indexing in the following exposition. Furthermore here we stay within a one-dimensional M R A . We show how to extend this framework to two-dimensional M R A in Appendix A . l . 3.1.1 Haar Wavelets The Haar scaling function is defined as the solution 0 of the following equation with ho = The Haar wavelet is defined based on orthogonality considerations with g0 = l / \ /2 ,g i = /ii = l/y/2: (3.1) (=0,1 (3.2) /=o,i 43 Haar M W 0 1 2 3 4 5 6 7 F i g u r e 3.1: Dependencies between different levels of the M R A . T h e nodes indicate scal ing functions and wavelets. W e define the scal ing funct ion and wavelet on the coarsest level $ 0,o (x) = 4>{x), "00,0 = Tp(x)- T h e hierarchy of basis functions is then obta ined by the d i la t ion equat ion, equat ion (2.10), where be define hkfik •= h0, hkfik+\ •= hu gk^k •= 9o a n d gk,2k+i •= 9i (corresponding to F i g u r e 3.1). Funct ions constant on the entire d o m a i n can be represented exactly in the H a a r basis. T h e H a a r wavelets have one vanish ing moment . W e show the graph of the H a a r scal ing funct ion and wavelet in F i g u r e 3.2. 3.1.2 Multiwavelets T h e Mult iwavele t M R A was presented by [alpe93]. It is based on Legendre po lynomia l s . Here we are concerned wi th multiwavelets based on l inear functions. W e can give explicit expressions for these functions (F igure 3.3). T h e hierarchy is constructed analogously to 44 Scal ing Funct ion phi Wavele t psi O -1 0 F i g u r e 3.2: H a a r scal ing funct ion and wavelet. the case of the H a a r wavelets using the d i la t ion equation. T h e wavelets have two vanishing moments . T h e y are able to reproduce l inear functions exactly. T h e space spanned by multiwavelet scal ing functions contains discontinuous functions. T h e coefficients for the generation of the hierarchy are l isted in the following table: M W Sca l ing Funct ions k I 'hki 0 0 1/V2 2 l/S/2 M W Wavelets 0 1 -l/y/2 3 l/y/2 0 -\/3/y/8 1 1/V8 2 3/78 3 l/y/8 0 l/y/8 1 V3/VS 2 -1/V8 3 \/3/s/S 45 Scaling Functions 3.1.3 C D V Wavelets The C D V M R A is presented in [cohe93a]. This M R A is derived from the Daubechies M R A [daub88] in order to make them span spaces with functions defined on a finite interval (the Daubechies wavelets are defined on 1R). Analogous to the Daubechies wavelets, C D V wavelets are available with different degrees of smoothness and different widths of support. Here we are solely concerned with the M R A derived from the Daubechies D 4 wavelet. The C D V scaling functions and wavelets in the interior of the interval, i.e. the func-tions with indices k : 2 < k < nm — 2, are identical to the D 4 scaling function and wavelet. On the boundaries, wavelets and scaling functions are modified, resulting in special edge' scaling functions <fo,0i A m - 2 A m - i and wavelets V ' o , V ' i > V ? n m - 2 , V ' n m - i -As in the case of Haar wavelets and multiwavelets, these coefficients also serve to generate the hierarchy using the dilation equation. The wavelets presented here have two vanishing moments. They can represent linear functions exactly. Functions in the spaces spanned by the wavelets are continuous but not differentiable. We show graphs of the C D V 46 scaling functions and wavelets in Figure 3.4. The D\ scaling function and wavelet are defined by a recurrence relation analogous to equation (3.1) and equation (3.2) with the following coefficients: h0 = (1 + v / 3)/(4 v / 2), hx = (3 + v /3)/(4v /2), h2 = (3 - v / 3)/(4v / 2), h3 = (1 - v /3)/(4v /2), So = -h3, gx = h2, 92 = -hi, g3 = hQ. For the functions at the boundaries we have the coefficients: k I hk, C D V Scaling Functions k - n m - 1 I - nm - 1 hk 0 6.033325e-01 1 6.908955e-01 2 -3.983130e - 01 0 3.751746e-02 1 4.573277e-01 2 8.500881e - 01 3 2.238204e - 01 . 4 -1.292227e- 01 - 1 -0 -1 -2 -o. -1' -2. -3 -4 8.705088e - 01 4.348970e - 01 2.303890e - 01 — 1.942334e - 01 1.901514e - 01 3.749553e-01 7.675567e - 01 4.431490e - 01 47 k I 9k,i C D V Wavelets k - n - 1 i - n m - l gktl 0 -7.965435169e-01 1 5.463927140e - 01 2 -2.587922483e - Ol 0 1.003722456e-02 1 1.223510431e-01 2. 2.274281117e - 01 3 -8.366029212e - 01 4 4.830129218e - 01 - 0 -0 -1 -2 -0 -1 -2 -3 -4 -2.575129195e - 01 8.014229620e-01 -5.398225007e - 01 3.717189665e - 01 -3.639069596e - 01 -7.175799994d- 01 4.010695194e - 01 2.315575950d- 01 3.1.4 Comparison Multiwavelets and C D V wavelets have higher order and more vanishing moments than Haar wavelets. They therefore potentially lead to a better compression and better approxima-tion than Haar wavelets. Multiwavelets can result in discontinuous approximations. C D V wavelets have the advantage of a higher degree of smoothness, which might result in better visual results, since the eye is particularly sensitive to discontinuities. For a qualitative comparison we show in Figure 3.5 how a function / : [0,1] —> R is approximated in the spaces spanned by the scaling functions. The function / is partially smooth but also has a discontinuity (we compare Haar, MW, and CDV bases): / ( * ) = atan(-Kx) 3x < 2 0 otherwise (3.3) 48 49 f M W G D \ / Figure 3.5: Projection of the function / of (3.3) into different spaces. A l l bases are with n = 32 functions. Haar basis on level j — 5, multiwavelet basis on level j = 4 and C D V basis on level j = 3. In the smooth part of / the approximation by higher-order bases is indistinguishable from the graph of / itself. f(x,y) = (3.4) In Figure 3.6 we show results for an equivalent experiment with a function: fxy 3(x + y) < 4 0 otherwise Here is the table for the compression rate of the two-dimensional case. We show the propor-tion of coefficients of wavelets that were below a given threshold e: e Haar MW1 CDV4 10" - l 100% 100% 99% 10" -2 96% 96% 88% 10" -3 73% 78% 71% 10" -4 38% 68% 52% 10" -5 24% 56% 33% 10" -6 16% 42% 25% 10" 12 16% 13% 9% 50 Figure 3.6: Projection of the function / of (3.4) into different spaces. A l l are bases with n — 1024 functions. Haar basis on level j = 5, multiwavelet basis on level j = 4, and C D V basis on level j = 3. We also compare the relative Lx error obtained for the function / defined by (3.4) to the relative Lx error obtained for the function g(x, y) = sjxy: Haar M W C D V f (discontinuous) 4.42e - 2 3.04e - 2 3.00e - 2 g (smooth) 1.72e - 2 1.87e - 3 1.31e - 3 For the discontinuous function the error in the approximation is very similar for all bases. The higher order bases have a small advantage since they approximate the smooth part of the function better. Both higher order bases produce wiggles along the discontinuity. The wiggles produced by CDV are less severe than those produced by MW. For the smooth function both higher-order bases show a clear advantage compared to the Haar basis. The results for the compression rate are also not unexpected. Due to the additional vanishing moment multiwavelets compress the smooth part of the function better. C D V wavelets, on the other hand suffer from the fact that more basis functions overlap the dis-continuity. This results in larger coefficients of such basis functions. In Chapter 4 we will experimentally analyze how these three types of wavelet bases 51 perform in hierarchical methods for the solut ion of the radiosity equation, where the process involves so lv ing an integral equation, not just funct ion a p p r o x i m a t i o n . 3.2 Solution Methods for Wavelet Radiosity A s discussed in Sect ion 2.3.3 K r y l o v subspace methods technical ly present an alternative to the P i c a r d i terat ion for Wavelet Radios i ty . In m a n y appl icat ions they per form better t h a n Picard's . T h e m a t r i x system in hierarchical methods wi l l generally not be S P D . In part icu lar , for higher-order basis functions we do not have a theorem which would allow us to make the system S P D . Therefore we choose the methods C G N R a n d G M R E S for our analysis. W e wi l l analyze the performance of these methods for the case of wavelet radiosity. In our implementa t ion of the methods we follow [saad95], where we replace the m a t r i x vector m u l t i p l i c a t i o n Mx by mMNcDN>mx a n d the m a t r i x vector m u l t i p l i c a t i o n MTx by mMj^eDNtmx, as described in Sect ion 2.2.3. T h e invest igation of fast so lut ion methods for the l inear a lgebra system in case of wavelet radios i ty is of interest. O n the one h a n d , wavelet radiosity requires keeping the entire sparse system Mrye in m e m o r y unlike progressive radiosity where so lut ion process and c o m p u t a t i o n of form factors are in termixed . So the solut ion process is an add i t i ona l c o m p u t a t i o n a l step, the cost of which seems to be worthwhile to m i n i m i z e . O n the other hand , in connect ion wi th the multi level- l ike methods from Section 2.3.3 a reduct ion of the number of i terations cou ld also reduce the n u m b e r of costly refinement steps needed. 3.3 Implementation Aspects Hierarch ica l methods for radiosity introduce the fol lowing types of a p p r o x i m a t i o n error into the c o m p u t a t i o n : 52 • d iscret izat ion error caused by the G a l e r k i n discret izat ion. • d iscret izat ion error caused by compress ing the G a l e r k i n discret izat ion. • error caused by the solut ion m e t h o d . • error caused by the form factor a p p r o x i m a t i o n . • error caused by a faulty es t imat ion of form factors. In our analysis we want to focus on the first three types of error. W e therefore designed our a lgor i thm a n d test environments to m i n i m i z e the last two types of error. W e do not use any of the mult i level- l ike methods since this would make the independent analysis of solut ion methods and discret izat ion difficult. In the fol lowing we describe how we realized the various components of the a lgor i thm for wavelet radiosity in the non-s tandard representation. 3.3.1 Galerkin Discretization In our experiments we use H a a r wavelets, mult iwavelets and CDV wavelets. W e use a separate M R A for each surface Si. H a a r wavelets have one vanishing m o m e n t and the abi l i ty to approx imate functions constant on the surface exactly. Mult iwavelets and C D V wavelets have b o t h two vanish ing moments and can approx imate all functions l inear on the surface exactly. A l l the wavelets and scal ing functions were made o r t h o n o r m a l in 1/2(5) by mul t ip l i ca t ion wi th y^ rjq ( s is the support of X^/clo* ^ o,k, s = 1 for H a a r a n d multiwavelets , s = 4 for C D V wavelets). 3.3.2 Form Factor Estimate T h e form factor is es t imated as described in Sect ion 2.4.2 for higher-order basis functions. F o r each of the four dimensions we compute the average of the dev iat ion over the three 53 r e m a i n i n g dimensions. T h e estimate is the m a x i m u m of these four averages. T o compute dev ia t ion from a constant funct ion (for a basis wi th one vanish ing moment ) we compute the difference of the kernel values at the boundaries of the interval . F o r the dev iat ion from a l inear funct ion (for bases wi th two vanishing moments) we interpolate by a l inear funct ion based on the kernel values at the boundaries of the interval a n d compute the difference of the value at the center of the interval and the kernel value at the center of the interval . 3.3.3 F o r m F a c t o r C o m p u t a t i o n W e do not deal w i th the p r o b l e m of singularities. W e design our experiments to exclude s ingular form factor integrals. Un l ike [gort93b] a n d [will97b] in most of their experiments , we use v i s ib i l i ty - in-quadrature for the c o m p u t a t i o n of form factors. W e t h i n k that in the context of higher order basis functions it is par t i cu lar ly impor tant to have a g o o d a p p r o x i m a t i o n to the form factors a long discontinuities since we have larger patches t h a n wi th piecewise constant basis funct ion . F r a c t i o n a l v is ib i l i ty does not allow to exploit the better a p p r o x i m a t i o n of higher order basis functions a long discontinuit ies since it essentially assumes constant v i s ib i l i ty w i t h i n ;a patch. W e wi l l support this c l a i m exper imental ly in chapter 4. W e il lustrate the procedure of how we approx imate the form factors for the one-d imens iona l case J f(x)w(x)dx, where / corresponds to the kernel, w to a basis funct ion . D u e to the irregularity of the C D V basis functions a s imple rule of s a m p l i n g the integrand f(x)g(x) is not appl icable in our case. W e therefore approx imate / by a piecewise l inear funct ion, and compute the integrals invo lv ing w separately. If we approx imate / by f(x) = di + biX if x : Xi < x < Xi+\ then the integral to compute is: fw(x) = 52 a i / w(x)dx + bi xw(x)dx 54 The following table gives values for the integrals on the right-hand-side for the Haar, multiwavelet and C D V scaling functions. We computed with equidistant samples over the support of 0 and nsamp — 3. The Haar and multiwavelet coefficients can be computed analytically. The C D V coefficients were computed numerically with a quadrature of high granularity based on an approximation to the C D V scaling functions computed via the cascade algorithm ([daub92].).. fxo 51! 0 ( x ) ^ x Ixl x^Wte HI xfr(x)dx Haar 0 1 2 I 2 I 8 3 8 • M W fro 1 2 1 2 1 8 3 8 01 1^3 2 4 ^ C D V 00 3.574485e - 01 4.603835e-03 -1.020265e - 01 6.665019e-03 01 1.078854e + 00 -7.740917e - 02 7.821774e - 01 -1.472866e-01 02 1.090248e + 00 -4.041994e - 04 8.495431e - 01 -1.586078e-01 03 1.384212e - 01 1.157058e + 00 1.806193e - 01 2.736973e + 00 D4 0. 1.077350e + 00 -7.734999c- 02 7.812263e-01 -1.472518e - 01 This table gives the integrals if w is a scaling function. If w is a wavelet then the integrals are computed via the dilation equation. In the two-dimensional case we have one scaling function and three wavelets. This results in 15 form factors to compute on levels j > 0 and 16 form factors on level j = 0 for each pair of patches. Following a computer graphics convention we call such a structure of 15 or 16 form factors a link. 3.3.4 M e m o r y R e p r e s e n t a t i o n o f t he Spa rse M a t r i x We associate a hierarchical graph corresponding to the graphs shown in Figure 3.1 with each surface. Each vertex of the graph is associated with a scaling function and three wavelets. Vertices of different hierarchical graphs can be connected by a link if the corresponding entry 55 in the m a t r i x was determined to be non-zero. Note that the hierarchy of the C D V M R A is not a tree structure. Therefore we designed the structure such that it is possible to check whether a l ink h a d already been c o m p u t e d to prevent redundant computat ions d u r i n g the top-down form factor c o m p u t a t i o n . 3.3.5 S o l u t i o n M e t h o d Unless otherwise noted we use P i c a r d i terat ion wi th 80 iterations for the solut ion of the system of equations. P i c a r d i terat ion is the s tandard m e t h o d used for h ierarchical methods a n d it converged stably in all our experiments. 56 Chapter 4 Numerical Experiments In contrast to the approach of deriving theoretical error bounds we try to analyze the behavior of the algorithms in the framework of characteristic scenes. Such an analysis can never be exhaustive due to the heterogeneous nature of possible input in computer graphics and the associated effects. Therefore we try to choose scenes with characteristic features encountered in practice which stress potential weaknesses and strengths of the analyzed algorithms. A particular behavior could then be a motivation for a more detailed analysis. 4.1 Comparison of Wavelet Bases We computed radiosity solutions for two simple test scenes. The scene "Unoccluded" is specified in Figure 4.1, the scene "Occluded" in Figure 4.2. We will measure error in the Li norm and additionally present images of the computed solution to allow for a visual comparison. 57 Figure 4.1: "Unoccluded". Surface So has constant reflectivity po = 0.4, no emission. Surface Si has zero reflectivity, i.e. px = 0, and constant emission e\ = 100. Note the proximity of the light-source to the surface S 0 , which causes a steep radiosity gradient. „ . - - - 0.4 1 Figure 4.2: "Occluded". Surface S 0 has constant reflectivity po = 0.4, no emission. Surface Si has constant reflectivity p\ = 0.4 and constant emission e\ = 100 in the indicated area. Surface S2 has zero reflectivity and has no emission. 58 4.1.1 Test Scenes T h e scene "Unocc luded" is the same scene that was used by [gort93b] for their measurements. It allows us to compute an analyt ic solut ion. T h e prox imi ty of emitter a n d receiver results in a h igh but smooth radiosity gradient. T h e light-source it rotated to prevent a l ignment w i th the basis-functions. T h e radiosity on face <Si is equal to the emission e, therefore, we get the fol lowing analyt ic expression for the radiosity: u(p) = p(p) r2(p,q)e(q)dSi(q) JSi F o r our reference so lut ion we approx imated this integral by a m i d p o i n t rule w i th 128 x 128 samples. T h e scene "Occluded" has a blocker that causes a shadow on SQ. L ike in the unoc-c luded scene we rotated the blocker to prevent al ignment of the shadow boundaries w i th the supports of the basis-functions. W e do not have an ana ly t i c so lut ion for this scene so for quant i tat ive results we compare wi th a master solut ion obta ined wi th a discret izat ion of high granularity . T h e master so lut ion was c o m p u t e d us ing multiwavelets wi th e = I O - 4 , w i th m = 5 levels (resulting in n$ = 4096 basis-functions) and nsamp = 3. In b o t h scenes all the interesting events take place on the surface SQ. Therefore , in the subsequent discussion we wil l only display the c o m p u t e d images for this surface. E r r o r s are also only c o m p u t e d for SQ. 4.1.2 Form Factor Computation T o support the c l a i m made in subsection 3.3.3 that the m e t h o d of fract ional v is ibi l i ty is not suitable for higher order basis-functions we c o m p u t e d a so lut ion wi th multiwavelets us ing b o t h techniques. W e show the results in F i g u r e 4.3. T h e y show that fract ional v is ibi l i ty fails to adequately represent shadow boundaries . A s imi lar observat ion was made by [will97b]. 59 V i s i b i l i t y - l n - Q u a d r a t u r e F r a c t i o n a l V i s i b i l i t y Figure 4.3: Result for the scene "Occluded" obtained using a basis of multiwavelets with n 4 = 1024 functions. The left image was obtained using visibility-in-quadrature, the right image using fractional visibility. 4.1.3 Measurements We measure the cost of the methods in terms of the number of form factors computed for the solution. The number of form factors reflects memory requirements as well as computational time. The cost of computing a form factor is determined by the number of kernel samples we need. For Haar wavelets we need at least two samples in each dimension to be able to estimate if the solution is constant. For first order basis functions we need at least three samples in each dimension to estimate linearity. Over four dimensions this means that the computation of a form factor for M W and C D V is more than 5 times more expensive than the computation of a form factor for the Haar basis. The solutions with a given number of form factors were obtained by first computing a sparse representation with a high number of form factors. Subsequently, the smallest coefficients were discarded until the desired number of form factors was reached. To be able to compare the quality of approximation we chose the number of levels such that we have the same number n of basis-functions. Furthermore we chose the following parameters for 60 U n o c c l u d e d Figure 4.4: "Unoccluded". Relative Lx error. Obtained by comparison to an analytical solution. our computations: e 777. n Haar I O - 5 3 5 n 5 = 1024 M W i l C T 4 3 4 n 4 = 1024 C D V 4 10~ 4 3 3 ™3 = 1024 Following [gort93b] we calculate the relative error in the Lx norm. These errors are shown for our two scenes in Figure 4.4 and Figure 4.5. We show images of the radiosity solution computed for surface So in Figure 4.8 and Figure 4.9. In Figure 4.6 and Figure 4.7 we show the distribution of the error for a computation with 20000 links. 4.1.4 Interpretation We observe that in case of the unoccluded scene acceptable solutions could be computed based on 2000 links (< 32000 form factors). That is far less than the number of form factors required for a standard Galerkin approximation with the same number of basis-61 Figure 4.5: "Occluded". Relative L\ error. Obtained by comparison to a master solution. Haar 20000 MW 20000 CDV 20000 Figure 4.6: "Unoccluded". Pointwise error for a computation with 20000 links. Obtained by comparison to an analytical solution. Black indicates small error, white indicates significant error. 62 Haar 20000 MW 20000 CDV 20000 Figure 4.7: "Occluded". Pointwise error for a computation with 20000 links. Obtained by comparison to a master solution. Black indicates small error, white indicates significant error. functions. A full matrix method in our case would require 10242 = 1048 576 form factors. This corresponds to a reduction to 3%. The compression looks worse in the occluded case. This is due to the fact that we chose the light-source for this scene not to cover the entire surface. Consequently numerous form factors were computed which do not contribute to the solution since the value of the solution is zero. Let us first focus on the error in the L i norm. Regarding the L\ error for the unoc-cluded scene multiwavelets perform best. They are followed by C D V wavelets. As expected the Haar basis shows the most significant error because the piecewise constant functions cannot approximate the smooth solution well. Multiwavelets approximate the solution well in areas with low gradient. In the area with high gradient the error is significantly higher. The same holds for C D V wavelets, however, due to the wider support the area with high error is larger. The behavior is different in the occluded scene. Here generally multiwavelets perform worst. C D V wavelets perform best, however for a higher number of links Haar is almost as good. The error for all three bases is high along the discontinuity. For Haar there is also error in the smooth areas simply because Haar is only a piecewise constant approx-63 exact MW 100 CDV 100 Haar 20000 MW 20000 CDV 20000 64 Figure 4.8: "Unoccluded". Images computed with the basis and number of form factors as indicated. The analytic solution is shown in the upper left corner. master MW 100 CDV 100 Figure 4.9: "Occluded". Images computed with the basis and number of form factors as indicated. The master solution is shown in the upper left corner. i m a t i o n . T h e error for mult iwavelets is also h igh in the s m o o t h areas. A t first glance this seems surpris ing since mult iwavelets should be able to approx imate the smooth function- a d -equately. However, it is a known phenomenon in G a l e r k i n discretizations when the solut ion cannot be approx imated well. In our experiment C D V wavelets are able to cope wi th that p r o b l e m better. Perceptual ly mult iwavelets and C D V wavelets per form very s imilarly. W e leave the judgement to the reader. T h e C D V wavelets introduce an unexpected gr id pat tern for smal l numbers of form factors (see F i g u r e 4.8 CDV2000). T h i s gr id pat tern is due to the G a l e r k i n discret izat ion itself. It cannot be seen in the s imple project ion shown in F i g u r e 3.6. W e also tr ied a higher n u m b e r of form factor samples nsamp, however, this d i d not result in any change. B o t h higher order bases produce clearly better results than the H a a r basis wi thout post - smoothing . C o n s i d e r i n g the wide support of the C D V wavelets their performance in our experiments for the occ luded scene is surpris ingly good . F o r the case of mult iwavelets in the occ luded case we note the difference in perceptual error and measured error. Perceptua l ly multiwavelets appear to per form better t h a n H a a r wavelets, not so however in terms of L\ error. Genera l ly we found it h a r d in our experiments to choose the appropr iate parameter e for the scene wi th occ lus ion. W e noticed a pretty abrupt change between an inaccurate so lut ion and a so lut ion which requires a large number of form factors, and consequently overly long computat ions . T h e selection of the parameter was more difficult for multiwavelets and C D V wavelets t h a n for H a a r . W e also tr ied the mult i leve l like refinement m e t h o d f rom subsection 2.3.3. Here we also found it h a r d to devise a good scheme for reducing e. 66 4.2 Comparison of Solution Methods We apply different solution methods to solve the system of equations resulting from a com-pressed kernel and evaluate their convergence. Our main interest is in the behavior of certain Krylov subspace methods (CGNR, G M R E S and CG) in comparison to the commonly applied Picard iteration. 4.2.1 Test Scenes As seen in Section 2.3 the Picard iteration generally converges rather rapidly for radiosity problems. One additional step of iteration corresponds to one additional reflection of light on a surface. So, scenes with significant degree of reflection potentially need a higher number of iterations. These scenes are where we anticipate other iterative methods to improve compared to the Picard iteration. We experiment with two geometries in two different physical configurations. The geometry of Figure 4.10 has no occlusion, the geometry of Figure 4.11 has. We compute the error for each geometry for surfaces with high average reflectivity pavg — 0.52 and with low average reflectivity pavg = 0.35: surface So Si s 2 s 3 s 4 s 5 s 6 s 7 s 8 s 9 Sio S i i h igh- 1 p; = 0.8 0.8 0.8 0.8 0.4 0.4 0.4 0.4 0.4 0.4(a)/0.8(6) 0.4 0 low Pi — 0.4 0.4 0.4 0.4 0.2 0.2 0.4 0.4 0.4 0.4(a)/0.8r» 0.4 0 Convergence of the Picard iteration in the high reflectivity setting is slower than in the low reflectivity setting due to the increased spectral radius of the system. We anticipate slower convergence of Picard for the scene with occlusion since there are surfaces which receive only indirect illumination and consequently need an additional step of reflection to receive light. 67 / A y / V 1 "0 5 2 10 F i g u r e 4.10: B o x scene without occlusion. A n u m b e r indicates the index i of a surface Si. A l l surfaces have constant reflectivity pi as specified in the text. T h e pos i t ion of the l ight-source is as specified in F i g u r e 4.2. F i g u r e 4.11: B o x scene wi th occlusion. A l l specifications are as for F i g u r e 4.10. T h e smal l cube is located at pos i t ion (1 ,1 ,1) inside the b ig cube. 68 4.2.2 Experiments In order to get a clear d is t inct ion between error caused by the discret izat ion and error caused by the solut ion m e t h o d we need to compare the solut ion xjk^ obta ined in the k-th step of the i terat ion to the exact solut ion x of the system. Since such a solut ion is not available we compare wi th a solut ion x obta ined after 80 steps of the P i c a r d i terat ion. T h e change in the solut ion after 80 iterations was in the order of 1 0 - 1 1 , which is way below the analyzed error: the latter was in the order of 1 0 - 3 or larger. W e compute the error without weighting by area since here we are only interested in the convergence of the solver, independent of the radiosity equat ion itself. W e use the fol lowing formula: E f c ( ^ - xky/n A l l computa t ions were conducted wi th a H a a r basis. A s in i t ia l guess for the itera-t ion we chose the discrete representation b of the light-sources. T h e results presented were c o m p u t e d for e = I O - 4 , however we also made measurements for e = 0.1 and e = 10~ 3 and obta ined s imi lar results. T h e m a x i m u m n u m b e r of levels is m = 5. F o r m factors were c o m p u t e d wi th nsamp = 3. T h e tota l n u m b e r of l inks c o m p u t e d in these computat ions is given in the fol lowing table: Scene U n o c c l u d e d O c c l u d e d # of potent ia l l inks 2 516 580 10 066 320 reflectivity H i g h L o w H i g h L o w # of c o m p u t e d l inks 484 452 283 272 898016 606 398 W e present graphs of the c o m p u t e d error in F i g u r e 4.12, F i g u r e 4.13, F i g u r e 4.14 and F i g u r e 4.15. error 69 t 0.0001 1 e - 0 6 b-U n o c c l u d e d H i g h # i t e r a t i o n s P i c a r d : o -G M R E S C G N R r Q -F i g u r e 4.12: Convergence of iterative methods for a test scene wi thout occlus ion a n d high average reflectivity. O c c l u d e d H i g h t 0 . 0 0 1 r-# i t e r a t i o n s F i g u r e 4.13: Convergence of iterative methods for a test scene wi th occlus ion a n d h igh average reflectivity. 70 0.1 0.01 0.001 0.0001 1 e-OS 1 e-06 1 e-07 1e-08 1 e-09 1 e-1 O Unoccluded Low i ar i 3icard ^IRES ;GNR' ' i 3~rr~rrr-^_^^. I ' 1 # iterations F i g u r e 4.14: Convergence of iterative methods for a test scene wi thout occlusion a n d low average reflectivity. Occluded Low # iterations F i g u r e 4.15: Convergence of iterative methods for a test scene w i t h occlus ion and low average reflectivity. 71 4.2.3 Interpretation T h e graphs clearly show that C G N R has no advantage over the P i c a r d i terat ion. In par-t icular , for the occ luded scene C G N R has a significantly larger error. O n the other h a n d , G M R E S has a clear advantage over P i c a r d i terat ion. Af ter three iterations the error is only ^ ^ of the error of P i c a r d in the unocc luded cases, | and w | for the occ luded scenes. In all scenes three iterations of G M R E S are sufficient for a perceptual convergence. T h e m a j o r cost in each step of P i c a r d i terat ion and G M R E S is one sparse matr ix -vector mul t ip l i ca t ion . In wavelet radios i ty this step includes wavelet decompos i t ion a n d wavelet reconstruct ion. G M R E S add i t iona l ly has a number of vector-vector mul t ip l icat ions the cost of which however is negligible c o m p a r e d to the operations ment ioned before. So in terms of cost G M R E S has a clear advantage over P i c a r d . C G N R on the other h a n d performs clearly worse t h a n P i c a r d since it is about twice more expensive per i terat ion t h a n P i c a r d . 72 Chapter 5 Conclusions W e have presented a framework of radios i ty methods which we hope wi l l help researchers in the areas of numerica l analysis and computer graphics to gain better access to the problems aris ing in the so lut ion of the radiosity equat ion. W e c o m p a r e d solut ion methods for wavelet radiosity as well as the behavior of different basis-functions. W e believe that this analysis wi l l contr ibute to a better unders tanding of wavelet radiosity. W e wi l l first give a brief s u m m a r y of results obta ined in our experiments. Subsequently we list a n u m b e r of questions which arose in our research, which we could not find answered in the l i terature a n d which remain to be pursued. 5.1 Results 5.1.1 Wavelet Bases In our compar i son b o t h examined higher order bases give perceptual ly better results than the H a a r basis. C D V wavelets tend to produce perceptual ly more accurate results t h a n multiwavelets . 73 We observed that C D V wavelets allow computation of a more accurate solutions with-out necessitating more links than needed for Haar wavelets. The number of links is reduced to up to half of the number of links required for Haar wavelets. In order to reduce memory consumption C D V wavelets appear to be a safe alternative to Haar wavelets. This advantage is important for hierarchical methods since we have to store the entire discretized system at a time. In terms of speed Haar wavelets outperform both higher-order wavelets even if we do not consider the faster methods available for the computation of Haar form factors. The computation of a higher-order link is intrinsically at least 5 times slower than the computation of a Haar link. This results in a clear disadvantage of higher-order wavelets in all observed cases. We confirm the observation made in [will97a] that fractional visibility produces inac-curate results with higher-order basis-functions. 5.1.2 Solution Methods We showed that G M R E S is an easy-to-realize alternative to the Picard iteration. For settings with high average reflectivity G M R E S converged in only 3 iterations to a solution with an error for which Picard needed more than 8 iterations. In low reflectivity settings G M R E S obtained a solution after 3 iterations that was reached by Picard after 5 — 6 iterations. In the context of wavelet radiosity one step of G M R E S and one step of Picard iteration have about the same cost. C G N R performed always slower than Picard iteration. 74 5.2 For Further Investigation T h e r e remain a number of unanswered questions and straightforward extensions to the work presented in the foregoing chapters. In the fol lowing we present a list of those which we consider to be most i m p o r t a n t . • W e d id not analyze G M R E S in the framework of the multi level- l ike methods f rom subsection 2.3.3. E m p l o y i n g G M R E S in this framework might lead to a significant ! speed-up since it can potent ia l ly result in a c o m p u t a t i o n of less form factors, however, • wi th possibly more iterations. A natura l quest ion then would also be the analysis of precondit ioners for G M R E S . A s imple precondit ioner would be the m a t r i x R ment ioned in subsection 2.3.1. • M o d e l i n g the solut ion wi th a m i x of different bases appears to be an interesting ap-proach. Areas w i th smal l gradient in the so lut ion can be well approx imated wi th a smal l number of higher order basis-functions while according to our experiments b ig gradients in the so lut ion are best approx imated w i t h the H a a r basis. Different bases appear to be easy to integrate as long as each surface has only one type of basis. A l s o , multiwavelets and C D V wavelets of higher t h a n first-order, can be integrated wi thout v m u c h effort. • A n impor tant element miss ing in the analysis is how well the oracle actual ly estimates the form factor, especial ly in the case of the scene h a v i n g occlusions. 75 Bibliography [alpe93] B. K . Alpert. "A class of bases in L2 for the sparse representation of integral operators". SIAM J. Math. Anal, Vol. 24, No. 1, pp. 246-262, January 1993. [appe85] A. A . Appel. "An efficient program for many-body simulation". SIAM J. Sci. Stat. Computing, Vol. 6, No. 1, pp. 85-103, 1985. [atki76] K . E. Atkinson. A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind. S IAM, Philadelphia, 1976. [bara95] G. V . Baranoski, R. Bramley, and P. Shirley. "Iterative Methods for Fast Radiosity Solutions". Technical report, Indiana University, Bloomington IN, http://www.cs.indiana.edu/ftp/techreports/index.html, 1995. [baum90] D. R. Baum and J. M . Winget. "Real Time Radiosity Through Parallel Processing and Hardware Acceleration". Computer Graphics (1990 Symposium on Interactive 3D Graphics), Vol. 24, pp. 67-75, March 1990. [beka96] P. Bekaert and Y . D. Willems. "Hirad: A Hierarchical Higher Order Radiosity Implementation". Proceedings of the Twelfth Spring Conference on Computer Graphics (SCCG '96), June 1996. [beyl91] G. Beylkin, R. Coifman, and V . Rohklin. "Fast Wavelet Transforms and Numerical Algorithms 1". Comm. Pure Appl. Math., Vol. 64, pp. 141-183, 1991. [boua95] K. Bouatouch and S. N . Pattanaik. "Discontinuity Meshing and Hierarchical Multiwavelet Radiosity". Proceedings of Graphics Interface '95, pp. 109-115, May 1995. [bric95] F. Bricout and E. Lepretre. "Distributed Progressive Radiosity on a Worksta-tion Network". Proceedings of Parallel and Distributed Processing Techniques and Applications (PDPTA '95), November 1995. 76 [chen89] Sh. E. Chen. "A Progressive Radiosity Method and its Implementation in a Dis-tributed Processing Environment". M.Sc. thesis, Program of Computer Graphics, Cornell University, Ithaca, NY, January 1989. [chia97] G. Chiavassa and P. Liandrat. "On the effective Construction of Compactly Sup-ported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval". Applied and Computational Harmonical Analysis, Vol. 4, pp. 62-73, 1997. [chri94] P. H. Christensen, E. J. Stollnitz, D. H. Salesin, and T. D. DeRose. "Wavelet Radiance". Fifth Eurographics Workshop on Rendering, pp. 287-302, June 1994. [cohe85] M . F. Cohen and D. P. Greenberg. "The Hemi-Cube. A Radiosity Solution for Complex Environments". Computer Graphics (ACM SIGGRAPH '85 Proceed-ings), Vol. 19, No. 3, pp. 31-40, July 1985. [cohe86] M . F. Cohen, D. Greenberg, D. Immel, and P. Brock. "An Efficient Radiosity Approach for Realistic Image Synthesis". IEEE Computer Graphics and Applica-tions, Vol. 6, No. 3, pp. 26-35, March 1986. [cohe88] M . F. Cohen, S. Chen, J. Wallace, and D. Greenberg. "A Progressive Refine-ment Approach to Fast Radiosity Image Generation". Computer Graphics (ACM SIGGRAPH '88 Proceedings), Vol. 22, No. 4, pp. 75-84, August 1988. [cohe93a] A. Cohen, I. Daubechies, and P. Vial . "Wavelets on the Interval and Fast Wavelet Transforms". Applied and Computational Harmonic Analysis, Vol. 1, pp. 54-81, 1993. [cdhe93b] M . F. Cohen and J. R. Wallace. Radiosity and Realistic Image Synthesis. Academic Press, Cambridge, 1993. [daub88] I. Daubechies. "Orthonormal Basis of Compactly Supported Wavelets". Comm. Pure. Appl. Math., Vol. 41, 1988. [daub92] I. Daubechies. Ten Lectures on Wavelets. S IAM, 1992. [feda92] M . Feda and W. Purgathofer. "Accelerating Radiosity by Overshooting". Third Eurographics Workshop on Rendering, pp. 21-32, May 1992. [garc97] B. Garcia and X . Pueyo. "Progressive Radiosity Solutions on SIMD Architec-ture" . Proc. First Eurographics Workshop on Parallel Graphics and Visualisation, September 1997. 77 [gers94]; R. Gershbein, P. Schroder, and P. Hanrahan. "Textures and Radiosity: Con-. trolling Emission and Reflection with Texture Maps". Computer Graphics (ACM SIGGRAPH '94 Proceedings), pp. 51-58, 1994. [gers95] R. Gershbein. "Integration Methods for Galerkin Radiosity Couplings". Rendering Techniques '95, pp. 264-273, June 1995. [gora84] C. M . Goral, K . E. Torrance, D. P. Greenberg, and B. Battaile. "Modelling the Interaction of Light between Diffuse Surfaces". Computer Graphics (ACM j SIGGRAPH '84 Proceedings), Vol. 18, No. 3, pp. 212-222, July 1984. [gort93a]-S. J. Gortler and M . F. Cohen. "Radiosity and Relaxation Methods". Technical report, Princeton University, 1993. [gort93b] S. J. Gortler, P. Schroeder, M . F. Cohen, and P. Hanrahan. "Wavelet Radiosity". Computer Graphics (ACM SIGGRAPH '93 Proceedings), August 1993. [gort94] S. J. Gortler, M . F. Cohen, and Philipp Slusallek. "Radiosity and Relaxation Methods". IEEE Computer Graphics and Applications, Vol. 14, No. 6, November 1994. [gort95] S. J. Gortler. Wavelet Methods for Computer Graphics. Ph.D. thesis, Technical Report, Department of Computer Science, Princeton University, Princeton, NJ , January 1995. [guit91] P. Guitton, J. Roman, and C. Schlick. "Two parallel approaches for a progressive radiosity". Second Eurographics Workshop on Rendering, May 1991. [gwin87] J. Gwinner. "On the Galerkin approximation of nonsmooth boundary inte-gral equations arising in radiative heat transfer". Boundary elements IX, Vol. 3 (Stuttgart, 1987), pp. 257-264. Comput. Mech., Ashurst Lodge, Ashurst, Southampton SO40 7AA, U K , 1987. [hack95] W. Hackbusch. Integral Equations. Birkhauser Verlag, Basel Boston Berlin, 1995. [hanr91] P. Hanrahan, D. Salzman, and L. Aupperle. "A Rapid Hierarchical Radiosity Algorithm". Computer Graphics (ACM SIGGRAPH '91 Proceedings), Vol. 25, No. 4, pp. 197-206, July 1991. [heck91] P. S. Heckbert and J. M . Winget. "Finite Element Methods for Global Illumi-nation". Technical report, Computer Science Division; University of California, Berkeley, July 1991. 78 [heck92] P. S. Heckbert. "Discontinuity Meshing for Radiosity". Third Eurographics Work-shop on Rendering, pp. 203-226, May 1992. [hest52] M . R. Hestenes and E. Stiefel. "Methods of conjugate gradients for solving linear systems". J. Res. National Bureau of Standards, Vol. 49, pp. 409-436, 1952.. [jawe94] B. Jawerth and W. Sweldens. "An Overview of Wavelet Based Multiresolution". SIAM Review, Vol. 36, No. 3, pp. 377-412, September 1994. [kaji86] J. T. Kajiya. "The Rendering Equation". Computer Graphics (ACM SIGGRAPH '86 Proceedings), Vol. 20, No. 4, pp. 143-149, 1986. [Iang95] E. Languenou, K. Bouatouch, and M . Chelle. "Global Illumination in Presence of Participating Media with General Properties". Photorealistic Rendering Tech-• •" niques (Proceedings of the Fifth Eurographics Workshop on Rendering), pp. 71-86, 1995. [Iisc92] D. Lischinski, F. Tampieri, and D. P. Greenberg. "Discontinuity Meshing for Accurate Radiosity". IEEE Computer Graphics and Applications, Vol. 12, No. 6, pp. 25-39, November 1992. [Iisc93] D. Lischinski, F. Tampieri, and D. P. Greenberg. "Combining Hierarchical Ra-diosity and Discontinuity Meshing". Computer Graphics (ACM SIGGRAPH '93 Proceedings), pp. 199-208, 1993. [neum95] L. Neumann and R. F. Tobler. "New Efficient Algorithms with Positive Defi-nite Radiosity Matrix". Photorealistic Rendering Techniques (Fifth Eurographics Workshop on Rendering), pp. 227-243, 1995. [niev97] Y . Nievergelt. "Making Any Radiosity Matrix Symmetric Positive Definite". Jour-nal of the Illuminating Engineering Society, Vol. 26, No. 1, pp. 165-172, 1997. [nish85] T. Nishita and E. Nakamae. "Continuous Tone Representation of Three-Dimensional Objects Taking Account of Shadows and Interreflection". Computer Graphics (ACM SIGGRAPH '85 Proceedings), Vol. 19, No. 3, pp. 23-30, July • 1985. [reck90] R. J. Recker, D. W. George, and D. P. Greenberg. "Acceleration Techniques for Progressive Refinement Radiosity". Computer Graphics (1990 Symposium on Interactive 3D Graphics), Vol. 24, pp. 59-66, March 1990. 79 [rush90] H. E. Rushmeier and K. E. Torrance. "Extending the Radiosity Method to In-clude Specularly Reflecting and Translucent Materials". ACM Transactions on Graphics, Vol. 9, No. 1, pp. 1-27, January 1990. [saad95] Y . Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 20 Park Plaza, Boston, M A , 1995. [scha97] S. Schaefer. "Hierarchical Radiosity on Curved Surfaces". Rendering Techniques '97 (Proceedings of the Eighth Eurographics Workshop on Rendering), pp. 187-192, 1997. [schr93] P. Schroder. "Numerical Integration for Radiosity in the Presence of Singularities". Fourth Eurographics Workshop on Rendering, Series E G 93 RW, pp. 177-184, June 1993. [schr94] P. Schroder. Wavelet Algorithms for Illumination Computations. Ph.D. thesis, Technical Report, Department of Computer Science, Princeton University, Prince-ton, NJ , November 1994. [schr96] P. Schroder. "Wavelet Radiosity: Wavelet Methods for Integral Equations". ACM SIGGRAPH '96 Course Notes - Wavelets in Computer Graphics, pp. 143-165. S I G G R A P H , 1996. [shao93] M . Shao and N . I. Badler. "Analysis and Acceleration of Progressive Refinement Radiosity Method". Fourth Eurographics Workshop on Rendering, E G 93 RW, pp. 247-258, June 1993. [shew94] J. R. Shewchuk. "An introduction to the conjugate gradient method without the agonizing pain". Technical report, Carnegie Mellon University, Pittsburgh, PA http://www.cs.cmu.edu/ quake/papers:html, 1994. [sill94] F. X . Sillion and C. Puech. Radiosity and Global Illumination. Morgan Kaufmann Publishers, San Francisco, 1994. [smit92] B. E. Smits, J. R. Arvo, and D. H. Salesin. "An Importance-Driven Radiosity Algorithm". Computer Graphics (ACM SIGGRAPH '92 Proceedings), Vol. 26, pp. 273-282, July 1992. [smit94] B. Smits, J. Arvo, and D. P. Greenberg. "A Clustering Algorithm for Radiosity in Complex Environments". Computer Graphics Proceedings, Annual Conference Series, 1994 (ACM SIGGRAPH '94 Proceedings), pp. 435-442, 1994. 80 [sout46] R . Southwel l . Relaxation Methods in Theoretical Physics. C l a r e n d o n Press, O x -ford, 1946. [stoe80] J . Stoer a n d R . B u l i r s c h . Introduction to Numerical Analysis. Spr inger Ver lag , N e w Y o r k B e r l i n Heide lberg , 1980. [stol96] E . Sto l ln i tz , T . Derose, and D . Salesin. Wavelets for Computer Graphics. M o r g a n K a u f m a n n Publ i shers , San Francisco , C a l i f o r n i a , 1996. [tell94] S. Tel ler , C . Fowler, T . Funkhouser , and P. H a n r a h a n . "Part i t ion ing and O r d e r -ing L a r g e R a d i o s i t y C o m p u t a t i o n s " . Computer Graphics (ACM SIGGRAPH '94 '., Proceedings), pp . 443-450, 1994. [trou93] R . T r o u t m a n a n d N . L . M a x . "Radios i ty A l g o r i t h m s U s i n g Higher O r d e r F i n i t e E l e m e n t M e t h o d s " . Computer Graphics (ACM SIGGRAPH '93 Proceedings), pp . . 2 0 9 - 2 1 2 , 1 9 9 3 . [veac97] E . Veach . Robust Monte Carlo Methods for Light Transport Simulation. P h . D . thesis, S tanford Universi ty , December 1997. [whit80]. T . W h i t t e d . " A n improved i l l u m i n a t a t i o n m o d e l for shaded display". Communi-cations of the ACM, V o l . 23, N o . 6, pp . 343-349, 1980. [will97a] A . W i l l m o t t a n d P. S. Heckbert . " A n E m p i r i c a l C o m p a r i s o n of R a d i o s i t y A l -gorithms" . Techn ica l R e p o r t C M U - C S - 9 7 , Carneg ie M e l l o n Univers i ty School of C o m p u t e r Science, 1997. [will97b] . A . W i l l m o t t and P . S. Heckbert . " A n E m p i r i c a l C o m p a r i s o n of Rad ios i ty A l g o -r i thms". Rendering Techniques 97 (Proceedings of the Eurographics Workshop on Rendering in St. Etienne, pp . 175-186, J u n e 1997. [xu94] W . X u and D . S. Fussell . "Cons truc t ing Solvers for Radios i ty E q u a t i o n Systems". Photorealistic Rendering Techniques (Proceedings of the Fifth Eurographics Work-shop on Rendering), pp . 207-217, June 1994. [zatz93] H . R . Z a t z . " G a l e r k i n Radios i ty : A H i g h e r - O r d e r Solut ion M e t h o d for G l o b a l I l luminat ion" . Computer Graphics (ACM SIGGRAPH '93 Proceedings), pp . 213-220, A u g u s t 1993. 81 Appendix A A . l Two-Dimensional Wavelet Transform We generalize the one-dimensional M R A presented in section 3.1 by forming products of scaling functions and wavelets. This results in a new two-dimensional scaling function and three new two-dimensional wavelets: ^toM^iV) = ^k0(x)'</>kl(y) ^koM&y) = MxtykM ^loM^'y) = ^k0(x)tpkAy) W i t h these definitions we get new two-dimensional dilation equations for scaling func-tions and wavelets: 7lj + l - l + 1 - 1 <t>j,k0M — 52 52 hko,i0hki,h<l>i+\,i0,h ;0=o /i=o T l j + l —1 T l j + l — 1 ^JMM = 52 52 9ko,l0hkuh<t>j+\,l0,h ;0=o /i=o ^jMM = 52 52 hko,io9kuh<!>j+i,i0M ;0=o /i=o 82 rij + i-1 n J + i — 1 i0=o ; 1 = o Bekaert et. al . [beka96] present a way to define two-dimensional M R A on non-rectangular grids. A.2 Non-Standard Representation For a given M R A the interaction of the projector Pm wi th T can be written as m—1 m—1 m—1 PmTPm = P0TP0 + J2 QJTQJ + ] T Q3TP3 + PJTQJ • ( A- 1) 3=0 3=0 j=0 This equality can be easily obtained using the identity Q3 = P3+i — Pj and the identities m— 1 PmTPm = P0TP0 + J 2 P J ^ T P J + i - p j T P J 3=0 and P3+1TP3+l-P3TP3 = (P3+l-P3)T(P3+l-PJ) + P3TP3+1 + P3+lTP3-2P3TP3 = QJTQJ \ PJTQJ rQjTPj. 83
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Wavelet radiosity in computer graphics
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Wavelet radiosity in computer graphics Ziegler, Philipp 1998
pdf
Page Metadata
Item Metadata
Title | Wavelet radiosity in computer graphics |
Creator |
Ziegler, Philipp |
Date Issued | 1998 |
Description | The thesis presents an overview of the recent development of radiosity methods in computer graphics in a uniform mathematical framework. The focus is on hierarchical methods using wavelets. The thesis experimentally analyzes the behavior of higher-order wavelet bases in hierarchical methods. The functions applied are multiwavelets and a family of wavelets proposed by Cohen, Daubechies and Vial. The latter wavelets have overlapping support. Generally, we find that higher-order wavelet bases save memory compared to Haar wavelets, while they require more time for the computation. Furthermore, we investigate how Krylov subspace methods can be employed to solve the discrete system of equations arising in hierarchical methods. We show that the Generalized Minimal Residual method (GMRES) is advantageous compared to the usually employed Picard iteration. |
Extent | 6627104 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-05-28 |
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.0051668 |
URI | http://hdl.handle.net/2429/8421 |
Degree |
Master of Science - MSc |
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-0684.pdf [ 6.32MB ]
- Metadata
- JSON: 831-1.0051668.json
- JSON-LD: 831-1.0051668-ld.json
- RDF/XML (Pretty): 831-1.0051668-rdf.xml
- RDF/JSON: 831-1.0051668-rdf.json
- Turtle: 831-1.0051668-turtle.txt
- N-Triples: 831-1.0051668-rdf-ntriples.txt
- Original Record: 831-1.0051668-source.json
- Full Text
- 831-1.0051668-fulltext.txt
- Citation
- 831-1.0051668.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-0051668/manifest