Wavelet Radiosity in Computer Graphics by P h i l i p p Ziegler Dipl.Inform., Universitat Kaiserslautern, Germany, 1996 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF Master of Science in THE FACULTY OF GRADUATE (Department STUDIES of C o m p u t e r Science) we a c c e p t this thesis as conforming to the r e q u i r e d s t a n d a r d The University of B r i t i s h 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. s i 2?, tm e( Date Computer Science The University of'British Co(um6ia 2366 Main maCC "Vancouver, "EC Canada yeT 1Z4 Abstract T h e thesis presents a n o v e r v i e w o f the recent d e v e l o p m e n t o f r a d i o s i t y m e t h o d s i n c o m p u t e r g r a p h i c s i n a u n i f o r m m a t h e m a t i c a l f r a m e w o r k . T h e focus is o n h i e r a r c h i c a l m e t h o d s u s i n g wavelets. T h e thesis e x p e r i m e n t a l l y a n a l y z e s the b e h a v i o r o f h i g h e r - o r d e r wavelet bases i n hiera r c h i c a l m e t h o d s . T h e f u n c t i o n s a p p l i e d are m u l t i w a v e l e t s a n d a f a m i l y o f wavelets p r o p o s e d by C o h e n , Daubechies and V i a l . T h e l a t t e r wavelets have o v e r l a p p i n g s u p p o r t . Generally, we find t h a t h i g h e r - o r d e r wavelet bases save m e m o r y c o m p a r e d t o H a a r wavelets, w h i l e t h e y r e q u i r e m o r e t i m e for the c o m p u t a t i o n . F u r t h e r m o r e , we i n v e s t i g a t e h o w K r y l o v s u b s p a c e m e t h o d s c a n b e e m p l o y e d to solve t h e d i s c r e t e s y s t e m o f e q u a t i o n s a r i s i n g i n h i e r a r c h i c a l m e t h o d s . W e show t h a t the G e n e r a l i z e d M i n i m a l R e s i d u a l m e t h o d ( G M R E S ) is a d v a n t a g e o u s c o m p a r e d to the u s u a l l y e m p l o y e d P i c a r d iteration. 11 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 Radiosity M o d e l 2.2 D i s c r e t i z i n g the R a d i o s i t y E q u a t i o n 8 2.2.1 Fundamental Problems 8 2.2.2 Basic Methods 9 2.2.3 Hierarchical Methods 13 2.2.4 Discontinuity Meshing 23 2.2.5 Hybrid Methods 24 2.3 . . 4 Solution Methods 25 2.3.1 25 Properties iii 2.4 2.5 3 Gauss-Seidel-like Methods 2.3.3 Other Methods 2.3.4 Comparison 26 • • 33 34 Computing the Discrete System 35 2.4.1 Form Factor Approximation 36 2.4.2 Form Factor Estimation 39 2.4.3 Comparison 39 Summary 40 Exploring Wavelet Radiosity 3.1 4 2.3.2 42 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 3.3.4 Memory Representation of the Sparse Matrix 55 3.3.5 Solution Method 56 • • Numerical Experiments 4.1 54 57 Comparison of Wavelet Bases 57 4.1.1 59 Test Scenes iv 4.2 5 4.1.2 Form Factor Computation . . . 4.1.3 Measurements 60 4.1.4 Interpretation 61 C o m p a r i s o n of S o l u t i o n M e t h o d s 67 4.2.1 T e s t Scenes 67 4.2.2 Experiments 4.2.3 Interpretation • • • 5.2 69 72 Conclusions 5.1 59 73 Results :-. 73 5.1.1 Wavelet Bases 73 5.1.2 Solution Methods 74 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 Transform 82 A.2 N o n - S t a n d a r d Representation 83 v : Acknowledgements I w o u l d like to t h a n k P r o f e s s o r U r i A s c h e r for his w o r k as m y research s u p e r v i s o r , i n p a r t i c u l a r for his fairness a n d for p r o v i d i n g a v e r y f r i e n d l y research e n v i r o n m e n t . F o u r n i e r gave c o n s t r u c t i v e c r i t i c i s m a n d p o i n t e d o u t v a l u a b l e resources. Professor Alain O n e of these b e i n g Ian A s h d o w n ' s g l o b a l i l l u m i n a t i o n b i b l i o g r a p h y w i t h o u t w h i c h w r i t i n g t h i s thesis w o u l d have t a k e n m u c h longer. L a s t b u t n o t least I w o u l d like to m e n t i o n L o r i n e C h a n g , w h o e n d u r e d all m y c o m p l a i n i n g , D a v e G r a v e s , w h o p a t i e n t l y answered a l l m y q u e s t i o n s r e g a r d i n g E n g l i s h style a n d g r a m m a r , as well as M i k e H o r s c h a n d A s h l e y W i j e y e r a t n a m , w h o c o n v i n c e d me t h a t a thesis has to be c o m p l e t e r a t h e r t h a n perfect. PHILIPP The University of British Columbia September 1998 vi ZIEGLER to / \0U5Y V 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 discretizations (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 equation 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. v e r y g e n e r a l f r a m e w o r k . W e use t h e i m p l e m e n t a t i o n to a n a l y z e the b e h a v i o r , o f h i g h e r - o r d e r wavelet bases c o m p a r e d to the c o m m o n l y u s e d H a a r basis, w h i c h consists of piecewise c o n s t a n t functions. T h e h i g h e r - o r d e r wavelets we use are :[alpe93] a n d the wavelets p r e s e n t e d i n [cohe93a]. T h e l a t t e r wavelets have a h i g h e r degree o f s m o o t h n e s s t h a n m u l t i w a v e l e t s . 2. T h e a p p l i c a b i l i t y o f K r y l o v m e t h o d s multiwavelets • •< '' as solvers for the sparse s y s t e m r e s u l t i n g f r o m t h e wavelet r a d i o s i t y m e t h o d a n d t h e i r convergence p r o p e r t i e s 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 t e r a t i o n . W e find t h a t h i g h e r - o r d e r wavelet basis f u n c t i o n s t e n d to have a p e r c e p t u a l a d v a n t a g e over the, H a a r basis a n d r e d u c e t h e r e q u i r e d a m o u n t o f m e m o r y . O n the o t h e r h a n d the s a m e error in.the. L\ n o r m c a n be r e a c h e d faster u s i n g a H a a r basis. , O u r . e x p e r i m e n t s s h o w t h a t the. K r y l o v m e t h o d G M R E S (see e.g. [saad95]) is a n easy- t o - r e a l i z e a l t e r n a t i v e to P i c a r d i t e r a t i o n w h i c h i n o u r e x p e r i m e n t s always p e r f o r m e d b e t t e r t h a n .the P i c a r d i t e r a t i o n . R o u g h l y s p e a k i n g we c a n r e a c h the s a m e a c c u r a c y as the P i c a r d i t e r a t i o n w i t h at m o s t h a l f t h e n u m b e r of i t e r a t i o n s . F o l l o w i n g this. i n t r o d u c t i o n t h e thesis presents a g e n e r a l f o r m a l f r a m e w o r k a n d a n overview o f e x i s t i n g r a d i o s i t y m e t h o d s i n C h a p t e r 2. In C h a p t e r 3 we present different t y p e s of wavelet bases a n d s o l u t i o n m e t h o d s . In C h a p t e r 4 we d e s c r i b e a n d j u s t i f y d e t a i l s o f o u r i m p l e m e n t a t i o n of wavelet r a d i o s i t y i n o r d e r to p r o p e r l y specify the n u m e r i c a l e x p e r i m e n t s conducted. C h a p t e r 5 c o n c l u d e s w i t h a s u m m a r y a n d e v a l u a t i o n o f the o b t a i n e d results. 3 Chapter 2 The Radiosity P r o b l e m and Its Numerical Solution 2.1 The Radiosity M o d e l T h e p r o b l e m o f g l o b a l i l l u m i n a t i o n c o n s i d e r e d here is h o w to o b t a i n t h e (i.e. angle-dependent radiance L(p,u>) e m i t t e d p o w e r of light p e r u n i t a r e a for a fixed w a v e l e n g t h ) leaving a p o i n t p o f t h e surface S o f a scene a t a s o l i d angle LO G £1 (tt is t h e set of s o l i d angles). T h e scene is d e s c r i b e d b y g e o m e t r i c independent. a n d physical properties. we a s s u m e t h a t computing c o l o r values c o r r e s p o n d s t o t h e c o m p u t a t i o n o f r a d i a n c e values for a sufficient number of wavelength of reflection B y c o m p u t i n g r a d i a n c e for a fixed w a v e l e n g t h T h i s f o r m u l a t i o n is v i e w - samples. is c a l l e d Lambertian emitters. terms of radiosity is defined as u(p) W e assume that emitted Lambertian reflection. r a d i a n c e is i n d e p e n d e n t type L i g h t - s o u r c e s w i t h t h i s p r o p e r t y are c a l l e d T h i s a s s u m p t i o n allows c h a r a c t e r i z i n g the i l l u m i n a t i o n o f t h e scene i n (i.e. t h e e m i t t e d p o w e r o f l i g h t p e r u n i t a r e a ) . — o f LO. T h i s J L(p,Lo)dLo. n T h e radiosity W e restrict o u r considerations 4 u:S — > R to interaction between surface p o i n t s a n d e x c l u d e effects t a k i n g place i n the m e d i a i n w h i c h the l i g h t is t r a n s f e r r e d . F u r t h e r m o r e , we w i l l n o t i n c l u d e t r a n s l u c e n t o b j e c t s 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 d e s c r i p t i o n of this p h y s i c a l m o d e l is g i v e n b y the radiosity equation, a n i n t e g r a l e q u a t i o n of the s e c o n d k i n d : u(p) = e(p) + JJ r(p, q)u(q)dS(q) (2.1) W e w i l l also w r i t e the r a d i o s i t y e q u a t i o n i n the s h o r t - h a n d n o t a t i o n u = e + Tu (2.2) r where T r is the c o r r e s p o n d i n g i n t e g r a l o p e r a t o r . T h e q u a n t i t i e s S, T h e set S C R 3 surfaces S^. e a n d r are defined as follows: is the surface of the scene. S = T h e surface is c o m p o s e d of m u t u a l l y d i s j o i n t |JSj. T h e f u n c t i o n e : S —> R d e s c r i b e s the l i g h t - s o u r c e s of t h e scene i n t e r m s o f e m i t t e d r a d i o s i t y . T h e s u p p o r t of this f u n c t i o n is u s u a l l y s m a l l c o m p a r e d to S since o n l y few areas of the surface o f a scene are l i g h t - s o u r c e s . The function r : S x S —> R is the • '- non-negative radiosity kernel given by r(p,q) = P(p)?~i(p,<7), w h e r e ri(2,q) = v (p,q)r (p,q) s cos Z ( n , q x 2 — p) cos l(n , p — q) y T2{p,q) |2 7T where v$ is the so c a l l e d o f t h e scene. visibility function. It satisfies v${p, p-q \ 2 T h e function vs d e s c r i b e s the geometry q) = 0 i f the l i n e s e g m e n t f r o m p to q intersects w i t h ( w h i c h m e a n s t h a t t h e r e are o b s t r u c t i o n s b e t w e e n p a n d q see F i g u r e 2.1), vs(p, otherwise. 5 q) — 1 F i g u r e 2.1: V i s i b i l i t y f u n c t i o n . T h e g r a y line shows p o i n t s T h e function p : S —> [0,1)] is the reflectivity. y vs{x,y) such that It d e s c r i b e s = 1. the p h y s i c a l p r o p e r t i e s of the scene, specifically, it describes the f r a c t i o n o f r a d i o s i t y reflected i n a p o i n t o f the surface. T h e vector n x 6 R 3 is the n o r m a l i n p p o i n t i n g to t h e t h e r e b y defined i n t e r i o r o f S. S i l l i o n a n d P u e c h [sill94] a r g u e b a s e d o n p h y s i c a l g r o u n d s t h a t \\T T \\ L < 1, w h i c h g u a r a n t e e s a u n i q u e s o l u t i o n for the r a d i o s i t y e q u a t i o n . A m o r e r i g o r o u s a n a l y s i s is g i v e n i n [gwin87]. F o r sufficiently s m o o t h surfaces S the s i n g u l a r i t y i n r for p = q c a n be r e m o v e d b y setting r(p,p) = 0. T h e k e r n e l is g e n e r a l l y d i s c o n t i n u o u s d u e t o the v i s i b i l i t y f u n c t i o n . S i n c e e is t y p i c a l l y d i s c o n t i n u o u s the s o l u t i o n u c a n be d i s c o n t i n u o u s as well. of the kernel r is expensive: n u m b e r o f surfaces. T h e computation N a i v e l y t h e cost o f c o m p u t i n g one v a l u e o f Vs is l i n e a r i n the T h i s w i l l t u r n o u t to be the m a j o r c h a l l e n g e i n s o l v i n g t h e e q u a t i o n . T h e s o l u t i o n o f e q u a t i o n (2.1) c a n be expressed i n t e r m s o f the N e u m a n n series: oo (2.3) i=0 6 T h e components Te l r o f t h e s u m have a n i n t u i t i v e p h y s i c a l i n t e r p r e t a t i o n : t h e y r e p - resent t h e c o n t r i b u t i o n o f the l i g h t - s o u r c e s after i reflections o n t h e surface. I m p o r t a n t effects t h a t c a n n o t be m o d e l e d w i t h the r a d i o s i t y e q u a t i o n are s p e c u l a r h i g h l i g h t s a n d i d e a l m i r r o r s (due to the a s s u m p t i o n o f diffuse light-sources reflection a n d L a m b e r t i a n [cohe93b]). T h e r a d i o s i t y e q u a t i o n has b e c o m e a s t a n d a r d m o d e l i n C o m p u t e r G r a p h i c s . E x h a u s tive i n t r o d u c t i o n s c a n be f o u n d i n [cohe93b] a n d [si!194]. N u m e r o u s m e t h o d s have b e e n s u g g e s t e d to solve the r a d i o s i t y e q u a t i o n . T h o s e m e t h o d s t y p i c a l l y have f o u r objectives: b e i n g fast, a c c u r a t e , r e q u i r i n g l i t t l e m e m o r y a n d b e i n g r o b u s t (i.e., r e q u i r e l i t t l e user i n t e r a c t i o n ) . H o w e v e r , s p e e d a n d a c c u r a c y c a n have different i n t e r p r e t a t i o n s i n c o m p u t e r g r a p h i c s t h a n i n the c o m m o n u n d e r s t a n d i n g o f n u m e r i c a l a n a l ysis. A m e t h o d w h i c h o b t a i n s a final result m o r e s l o w l y b u t p r o d u c e s useful results m a y be preferable t o a n o v e r a l l faster m e t h o d . intermediate S u c h m e t h o d s are called' progressive: A less a c c u r a t e result w h i c h r e p r o d u c e s c e r t a i n v i s u a l effects, like s m a l l shadows, well m i g h t be m o r e useful t h a n a m e t h o d w i t h a n o v e r a l l s m a l l e r error i n a s i m p l e error n o r m . In the f o l l o w i n g sections we w i l l d e s c r i b e p r i n c i p l e s o f i m p o r t a n t d e t e r m i n i s t i c m e t h o d s for t h e s o l u t i o n o f the r a d i o s i t y e q u a t i o n . F o r a n o v e r v i e w o f s t o c h a s t i c m e t h o d s we refer to [veac97]. kernel. S e c t i o n 2.2 is c o n c e r n e d w i t h h o w m e t h o d s d i s c r e t i z e t h e r a d i o s i t y f u n c t i o n a n d S e c t i o n 2.3 presents m e t h o d s for the s o l u t i o n o f the d i s c r e t i z e d s y s t e m . In s e c t i o n 2.4 we discuss how t h e d i s c r e t i z e d s y s t e m itself, the f o r m factors, c a n be c o m p u t e d . chose to present the t h r e e c o m p o n e n t s We of discretization, solution m e t h o d a n d form factor c o m p u t a t i o n i n d e p e n d e n t l y , since it allows us to p o i n t o u t m o r e c l e a r l y v a r i o u s a n d d i s a d v a n t a g e s o f single c o m p o n e n t s of a r a d i o s i t y m e t h o d . 7 advantages 2.2 Discretizing the Radiosity Equation B e f o r e p u b l i c a t i o n s of [gwin87] a n d e s p e c i a l l y [heck91], w h i c h i n t r o d u c e d G a l e r k i n m e t h o d s ; for t h e s o l u t i o n o f the r a d i o s i t y e q u a t i o n , r a d i o s i t y m e t h o d s were b a s e d p u r e l y o n p h y s i c a l c o n s i d e r a t i o n s o f transfer o f energy b e t w e e n surface p a t c h e s g o v e r n e d b y t h e p r i n c i p l e of energy balance. T h e p r e s e n t a t i o n here gives a u n i f o r m o v e r v i e w i n t h e f r a m e w o r k of t h e . n o t a t i o n i n t r o d u c e d i n s u b s e c t i o n 2.2.1. In the f o l l o w i n g the d o m a i n S w i l l be c o m p o s e d of p l a n a r surfaces Si (see [scha97] for a n o v e r v i e w o f m e t h o d s for c u r v e d surfaces). : Fundamental Problems 2.2.1 T h e u n d e r l y i n g p r i n c i p l e o f a l l the p r o c e d u r e s p r e s e n t e d here is the G a l e r k i n m e t h o d . ; T h i s method approximates d i m e n s i o n a l space V functions. Si. the e x a c t s o l u t i o n u of e q u a t i o n = (2.1) span{</> • • • 4>n-i}., where fa . : S 0 T h e s u p p o r t ^ o f each fa, also c a l l e d a patch, T h e reflectivity p is a s s u m e d t o have c o n s t a n t by a function u in a —>• R are l i n e a r l y finite independent is e n t i r e l y c o n t a i n e d i n a surface v a l u e pi o n p a t c h Ai (for the efficient t r e a t m e n t of n o n - c o n s t a n t reflectivities see [gers94]). W e define the p r o j e c t i o n P onto V v E for L (S): 2 n-l (2.4) where <v,w> = f u n c t i o n v = Pv J v(p)w(p)dS(p) s is the s c a l a r - p r o d u c t i n m i n i m i z e s || v — v \\ L2 in L (S). 2 It c a n be s h o w n t h a t the 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 a p p r o x i m a t e s u b y u € V s u c h t h a t the r e s i d u a l o f the i n t e g r a l e q u a t i o n p r o j e c t e d i n t o V approximation: 8 is m i n i m i z e d . T h i s results i n the following D e f i n i t i o n 2.2.1 The equation Pu = Pe + is called B u b n o v - G a l e r k i n With u = approximation Y^=o i t i^ P x ( ) PT Pu r of the integral equation u = e + T u. r YTiZl h^i, e — a d e x p a n s i o n of the p r o j e c t i o n s , the B u b n o v - n G a l e r k i n a p p r o x i m a t i o n results i n a discrete s y s t e m of equations: 71-1 • Xi = h + ^2fijXj T h e coefficients fa = .9ij = ^ , ,<T (j) ^ > <o o , > T 3 forrri factors are c a l l e d = l i i: i = 0 , 1 , . . . ,n - 1 1 «pi,.0i> where if j JJ JJ Ai {fij)nxn, k= (b ,... 0 W e call ,6„_i) . T F (gij) xn n is defined as follows: r (p,q)(j) {p)(f) (q)dA (q)dA (p) - 1 Aj W e w i l l w r i t e the s y s t e m i n m a t r i x n o t a t i o n as = G= (2.5) i j J t Mx — b w i t h M = (mij) nxn the m a t r i x r e p r e s e n t a t i o n of • (2.6) = / — F, PT P. r T h e f o r m factors p h y s i c a l l y relate to the transfer of energy f r o m p a t c h Ai to p a t c h T h e i r c o m p u t a t i o n is w h a t m a k e s the r a d i o s i t y p r o b l e m so t i m e - c o n s u m i n g . m a t i o n of gij involves the c o m p u t a t i o n of s a m p l e s of the v i s i b i l i t y f u n c t i o n c o s t l y (as w i l l be d i s c u s s e d i n S e c t i o n 2.4). s i m p l e scenes (e.g. [cohe86]) to 1 0 6 2.2.2 T h e approxi- v {p_,q) s w h i c h is N u m b e r s of patches range f r o m a b o u t 1000 for for c o m p l e x scenes (e.g. the m e t h o d s p r e s e n t e d is to c o m p u t e a n a p p r o x i m a t i o n M the coefficients g^ Af. [tell94]). A m a j o r o b j e c t i v e of of M w h i c h requires as l i t t l e of as p o s s i b l e . Basic Methods Classical Radiosity C l a s s i c a l r a d i o s i t y , also c a l l e d i n t o d i s j o i n t surface p a t c h e s full matrix radiosity, A. k T h e operator 9 T r starts w i t h a s u b d i v i s i o n of the surface is a p p r o x i m a t e d by PT P, r S hereby directly r e a l i z i n g the G a l e r k i n a p p r o a c h d e s c r i b e d i n s u b s e c t i o n 2.2.1. T h e m e s h has to be a d a p t e d t o t h e s o l u t i o n s u c h t h a t s h a d o w b o u n d a r i e s a n d h i g h r a d i o s i t y g r a d i e n t s c a n be a d e q u a t e l y r e p r e s e n t e d . T h e m e t h o d uses basis f u n c t i o n s w h i c h are c o n s t a n t o n the patches: { : p e A 1 k (2-7) 0 : otherwise W i t h these basis f u n c t i o n s we o b t a i n for the f o r m factors: In case of piecewise c o n s t a n t basis f u n c t i o n s the f o r m f a c t o r s represent the f r a c t i o n o f e n e r g y t r a n s f e r e d f r o m p a t c h k to p a t c h / . C l a s s i c a l r a d i o s i t y r e q u i r e s user i n t e r a c t i o n for the g e n e r a t i o n of the m e s h . T h e n u m ber of f o r m f a c t o r s t h a t have to be c o m p u t e d is q u a d r a t i c i n the n u m b e r o f p a t c h e s . c o m p u t e d s o l u t i o n has d i s c o n t i n u i t i e s d u e to the t y p e o f basis f u n c t i o n s t h a t we use. The These d i s c o n t i n u i t i e s are s m o o t h e d o u t i n a p o s t - p r o c e s s i n g step t h a t l i n e a r l y i n t e r p o l a t e s b e t w e e n t h e centers of the c o m p u t e d patches. T h e p o s t - p r o c e s s i n g step is n o t i n t e n d e d to decrease the n u m e r i c a l error b u t r a t h e r a i m s at m a k i n g t h e o b t a i n e d result v i s u a l l y m o r e p l e a s i n g . P a t c h e s c a n have a r b i t r a r y s h a p e a n d so t h e y c a n be f l e x i b l y a d j u s t e d to the p r o b l e m . T h i s m e t h o d was first p r e s e n t e d for c o n v e x scenes (scenes w i t h o u t o b s t r u c t i o n ) b y [gora84]. It was t h e n g e n e r a l i z e d to n o n - c o n v e x scenes b y [nish85] a n d [cohe85]. Substructuring Sub structuring a i m s at r e d u c i n g the n u m b e r o f f o r m f a c t o r s to c o m p u t e . F o r m f a c t o r s are i n i t i a l l y c o m p u t e d o n a coarse m e s h . A n a p p r o x i m a t e s o l u t i o n is f o u n d o n t h i s coarse m e s h . T h e n the s o l u t i o n is t r a n s f e r e d to a n i n i t i a l fine m e s h , w h i c h is s u b s e q u e n t l y 10 iteratively refined b a s e d o n the g r a d i e n t (e.g. c o m p u t e d b y a c e n t e r e d difference) o f the s o l u t i o n f o u n d so far. :;• T h e m e t h o d uses piecewise c o n s t a n t basis f u n c t i o n s 0 { level: O n the finer level there are k < no, 0 < q < q 1 : Pe A 0 : otherwise O A . : 0 < k < n 0 o n the coarse 0tk ri\ = X ^ l o * Qk piecewise c o n s t a n t basis f u n c t i o n s 4>i,(k,q) such that: k 9/fc-l <t>o,k = ^2 01,(/c, ) 9 T h e refinement process is l o c a l a n d the i n i t i a l m e s h is sufficiently fine t h a t q = k for m o s t k. T h e projections P a n d Pi 0 f u n c t i o n s <^n,. a n d </> . T h e o p e r a t o r lr h,(k,q),i = are defined as i n e q u a t i o n (2.4) P\T P r 0 b a s e d o n the 1 basis is d e s c r i b e d b y t h e f o r m factors: ^l\M',{k,q)\'ll II r(p,q)dAo^q)dA ( (p) h T h e f o r m factors o f PQT PQ c a n b e c o m p u t e d b a s e d o n A;i,(i ). T l9 k>q) J" Oi-l k j = l/\A \ 0ii The a p p r o x i m a t e s o l u t i o n ui o f e q u a t i o n s UQ U\ 0>l = Pof + P T P uo 0 = Pif + PiT P u . r 0 r 0 ^ ^i,(»,g)jl^l,(z,g)l is t h e n o b t a i n e d b y a n first s o l v i n g the discrete d t h e n t r a n s f e r r i n g t h i s s o l u t i o n o n t o the T h i s l a t t e r s t e p is r e p e a t e d w h i l e r e f i n i n g 0 g r a d i e n t decreases b e l o w a t h r e s h o l d v a l u e e. P i c a r d i t e r a t i o n for PIT PQ T Pi system fine mesh u n t i l the r a d i o s i t y In t o t a l t h e s e c o n d s t e p is j u s t one step o f a n d as s u c h i m p r o v e s t h e s o l u t i o n u\ i n c o m p a r i s o n to UQ. T h e p h y s i c a l i n t u i t i o n o f t h i s a l g o r i t h m is t h a t UQ is a sufficiently g o o d a p p r o x i m a t i o n o f the average r a d i o s i t y o n the coarse m e s h . T h i s coarse m e s h however is n o t able to r e p r o d u c e s h a r p s h a d o w edges. T o increase t h e s t r u c t u r e o f t h e r a d i o s i t y s o l u t i o n w i t h i n a coarse 11 0 < patch A 0ti Ai^ y we transfer r a d i o s i t y f r o m t h e o t h e r coarse p a t c h e s to the fine p a t c h e s iq S o m e fine p a t c h e s m i g h t n o t receive r a d i o s i t y since t h e y are o b s t r u c t e d b y a n o b j e c t . g o a l o f t h i s o p e r a t i o n is to i m p r o v e the s h a d o w The representation. S u b s t r u c t u r i n g was i n t r o d u c e d i n [cohe86]. T h e n u m b e r o f f o r m factors c o m p u t e d is s t i l l q u a d r a t i c , 0(nl -f-nnni), however, no is m u c h s m a l l e r t h a n i n f u l l - m a t r i x r a d i o s i t y . T h e p r o c e s s s t i l l requires user i n t e r a c t i o n to d e t e r m i n e the coarse m e s h a n d the i n i t i a l fine m e s h . T h o s e meshes have to be fine e n o u g h a n d a d a p t e d to the s o l u t i o n . m e s h is not p r o p e r l y a d a p t e d to the s o l u t i o n t h e n refinement radiosity gradient. If the the i n i t i a l fine takes p l a c e i n areas o f s m a l l If the coarse m e s h is not fine e n o u g h v i s u a l a r t i f a c t s m i g h t a p p e a r like the reflection o f r a d i o s i t y f r o m s h a d o w s . T h e final step w h i c h transfers the s o l u t i o n o n t o a fine g r i d c a n o n l y r e d i s t r i b u t e r a d i o s i t y w i t h i n a coarse p a t c h b u t n o t a m o n g coarse T h i s m e t h o d is a n a t t e m p t to a p p r o x i m a t e the d i s c r e t i z e d k e r n e l b y a n n\ n 0 x n\ PiT P\ r patches. (represented m a t r i x ) b y a sparser k e r n e l r e q u i r i n g fewer f o r m f a c t o r s (represented x no a n d a n ri\ x n 0 matrix). by an T h e c o m p u t a t i o n o f this a p p r o x i m a t i o n is d r i v e n b y the approximate solution. Higher-Order Bases G a l e r k i n m e t h o d s w h i c h use h i g h e r - o r d e r basis f u n c t i o n s have the p o t e n t i a l o f decreasing t h e n u m b e r o f necessary f o r m factors. T h e y p r o v i d e a b e t t e r a p p r o x i m a t i o n o f the r a d i o s i t y s o l u t i o n i n areas w h e r e it is s m o o t h . A s m a l l n u m b e r of h i g h e r - o r d e r basis f u n c t i o n s a p p r o x i m a t e s u c h a f u n c t i o n well, w h e r e a s for piecewise c o n s t a n t finement w o u l d be necessary. basis f u n c t i o n s t r o n g re- F u r t h e r m o r e , w i t h h i g h e r - o r d e r basis f u n c t i o n s it is to use a p p r o x i m a t i o n spaces w i t h c e r t a i n s m o o t h n e s s p r o p e r t i e s . can possible Such an approximation m a k e s the p o s t - p r o c e s s i n g step r e d u n d a n t . W h i l e h i g h e r - o r d e r m e t h o d s w i t h basis f u n c t i o n s w i t h w i d e s u p p o r t are useful i n areas 12 of s m o o t h v a r i a t i o n p a t c h e s have to be refined c a r e f u l l y a r o u n d d i s c o n t i n u i t i e s . : H i g h e r - o r d e r m e t h o d s have b e e n a p p l i e d to t h e r a d i o s i t y p r o b l e m b y [heck91], [trou93] a n d [zatz93]. T h e m e t h o d p r e s e n t e d b y [zatz93] uses L e g e n d r e p o l y n o m i a l s as basis f u n c t i o n s a n d is c a l l e d 2.2.3 Galerkin radiosity. Hierarchical Methods T h i s section gives a b r i e f m o t i v a t i o n for h i e r a r c h i c a l m e t h o d s . m e t h o d s is to find T h e general idea of a h i e r a r c h i c a l r e p r e s e n t a t i o n o f the a p p r o x i m a t e o p e r a t o r PT P. r the The h i e r a r c h i c a l r e p r e s e n t a t i o n c o n t a i n s f o r m factors b e t w e e n h i e r a r c h i e s o f p a t c h e s . W h i l e there w i l l g e n e r a l l y be significant l i g h t t r a n s f e r , a n d c o n s e q u e n t l y b i g f o r m f a c t o r s , b e t w e e n b i g p a t c h e s w h i c h are close to each o t h e r , the m a g n i t u d e o f transfer d e c a y s w h e n p a t c h e s are s m a l l e r a n d are f u r t h e r a p a r t . F u r t h e r m o r e we w a n t to find a basis s u c h t h a t f o r m factors b e t w e e n p a t c h e s w i t h s m a l l v a r i a t i o n o f the k e r n e l are s m a l l themselves. W e c a n t h e n set a l l f o r m factors b e l o w a c e r t a i n t h r e s h o l d e to zero, t h e r e b y o b t a i n i n g a sparse r e p r e s e n t a t i o n M. £ A p a r t i c u l a r c h a l l e n g e is to o b t a i n t h e sparse r e p r e s e n t a t i o n w i t h o u t h a v i n g to c o m p u t e a l l the f o r m f a c t o r s b e f o r e h a n d , i.e. to d e c i d e w h i c h f o r m factors w o u l d be b e l o w the t h r e s h o l d without computing them first. F o r the sake o f c l a r i t y we t e m p o r a r i l y s w i t c h to the a b s t r a c t case o f a o n e - d i m e n s i o n a l i n t e g r a l e q u a t i o n w i t h k e r n e l s : R x R —» R s u c h t h a t the f o l l o w i n g s m o o t h n e s s r e q u i r e m e n t s are fulfilled: \s(x,y)\ < 7 - f 1 V >y>\ - \x-y\ T h e s e r e q u i r e m e n t s a l l o w for one d i s c o n t i n u i t y a l o n g the d i a g o n a l x = y. i n g d i s c u s s i o n follows [beyl91]. 13 ^ T h e follow- H i e r a r c h i c a l basis functions W e define h i e r a r c h i c a l basis f u n c t i o n s w i t h c e r t a i n p r o p e r t i e s r e q u i r e d for the p r o c e d u r e d e s c r i b e d below. F i r s t we define a h i e r a r c h y of finite d i m e n s i o n a l spaces Vo C V\ C ... C V . Further- m m o r e we choose spaces Wj s u c h t h a t ( © is the d i r e c t s u m ) : w ®v J T h e space = v j J+l Vj is defined b y basis f u n c t i o n s 4>j^'- Vj = span E q u i v a l e n t l y w e choose basis f u n c t i o n s "0j,fc ' 0 < k < rij for the spaces : 0 < dilation fij}. T h i s h i e r a r c h y is equation: $j,k = YALV ^j,k = J2ZV~ 1 k,l&j+l,l ^ h L 9k,i(f>j+i,i W e represent f u n c t i o n s Vj € VJ; VJ = X ^ L o * 3,k<t>3,k '• 0 < j c coefficients c- = < Wj. T h e , h i e r a r c h y of spaces i m p l i e s a h i e r a r c h y o f basis f u n c t i o n s . expressed b y t h e k (cj^oKkKnj the v e c t o r o f coefficients and functions = (dj^)o</c<n If we s t a r t u p w i t h a f u n c t i o n v e WJ; WJ = YTk=o dj, 4>j,k Wj k '• b y t h e . v e c t o r of 0 < j < m , by r represented by c m < m, m i n the finest space V t h e n we c a n m successively d e c o m p o s e it i n the f o r m v = u > _ i + u> -2 + • • • + v~>o + VQ. T h e c o m p u t a t i o n m oi Cj 0 < : j If the f u n c t i o n s Cj '• 0< reconstruction. j < m < a n d rf- : m 0 < j < from m <f>j k are o r t h o g o n a l for a fixed and d, t 3 : 0 < j < m m is p o s s i b l e . If we r e q u i r e t h a t for fixed k j c m is c a l l e d a wavelet then the c o m p u t a t i o n of c decomposition. m from given T h i s reverse o p e r a t i o n is c a l l e d a o n l y a c o n s t a n t n u m b e r o f t h e coefficients wavelet hk i t a n d gkj be n o n - z e r o t h e n b o t h c o m p u t a t i o n s c a n be c o n d u c t e d i n a s t r a i g h t f o r w a r d f a s h i o n in 0 ( n ) m o p e r a t i o n s v i a t h e 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 D N>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 C to ^ do D N,mQ.n (2.12) D N.m tm-1 \ dm-l \ J The matrices Ds^Dg^, and D^ m dm-l J are never formed explicitly. They merely serve as a convenient notation for wavelet decomposition and reconstruction. Ds, m transfers the representation c rn 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 decomposition. Dyv im Note that is called a standard 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 k^l, then P j+1 0^ for fixed j, i.e. <<j>j,k, ^>j,k> = Pj + Qj. It is possible to define sets of functions ct>j^ and ipj,k such that: J ip (x)x^dx j>k = OVfj, : 0 < fi< 15 v = ^-,< f j,k, <fij,i> — ( ) tpj^ T h e function is t h e n s a i d to have v vanishing moments. W e w i l l give e x a m p l e s of s u c h f u n c t i o n s i n S e c t i o n 3.1. The qt>j k functions t and ipj k, t w h i c h we have defined here, fit i n t o the f r a m e w o r k o f o r t h o g o n a l wavelets d i s c u s s e d , e.g., i n [daub88] a n d [jawe94]. I n t h i s f r a m e w o r k the f u n c t i o n s tpj k t are c a l l e d wavelets, a n d (pj^ the f u n c t i o n s scaling functions. are c a l l e d multiresolution analysis wavelets a n d s c a l i n g f u n c t i o n s is c a l l e d T h e pair of (MRA). W e i g n o r e the m o r e g e n e r a l f r a m e w o r k of b i o r t h o g o n a l wavelets (e.g. [jawe94]) since we d o n o t a p p l y it i n this thesis. Wavelet Compression of Integral Operators W i t h the a b o v e p r o p e r t i e s we c a n express P as: m 771—1 Pm Po + ^2 Qj — j=0 A s a c o n s e q u e n c e we c a n r e w r i t e P TP : m 771—1 P TP m m = P TP 0 0 + P T(J2 j m 771— 1 771—1 Qj) + ($2 ^ j 0 Q 0 0 T P + C j O A c c o r d i n g to e q u a t i o n (2.13) a n o t h e r way to w r i t e d = m w h e r e the m a t r i x Fs QJ) (J2 Qi) T 0 P TP v m 771—1 m (- ) 2 13 j=0 i n m a t r i x n o t a t i o n is s)n sDs,mC D F m is defined as F,s and F( °' °} p p w i t h entries w i t h entries o f the f o r m <(f>o,k,Tipjj>, and the standard representation o f <T(p , </> , >, F^ °' ^ P 0>l F^'^ 0 fc w i t h entries P TP . m m 16 Q w i t h entries <Tipj,i,tpi,k>- <Tip <po, >, F^ ' ^ Q P jih T h e matrix k Fs is c a l l e d i=m-7 QTQ | PTQ QTP QTQ PTQ QTP Figure 2.2: Non-standard representation F of N Another way to rewrite P TP m m m = P TP 0 + 0 m m (see Appendix A.2) is: m m—1 P TP P TP . i=o m—1 QjTQj- + m—1 Qi i + 52 i ®i TP j=o • m p (<T^- /,Vj fc>)o<fc<no,0<Z<no. i d ) <3 m N Ntin m m Qi F >' > = (<Tlp (Q P TP : = (<T<f> i, <l>o,k>)o<k<n ,o<i<no, F&'W = Po N n 14 D^ F D c The matrix F is composed of blocks F^ °' ^ a 2 j=o With equation (2.14) we get the non-standard representation of dm = (- ) P T ) jth 0 Ipj,k>)o<k<n ,0<l<n 0 0 • The structure of Fjv is shown in Figure 2.2. Analyzing the sparsity of the matrices F, Fs and F we obtain the following numbers N of non-zero entries for the typical case of a dyadic hierarchy (n i = representation matrix non-zero entries full matrix F standard F n non-standard F n J+ n 2 s 2 m 2 N •m 11 17 - n l = n l ( l - 2 - ^ ) 2rij): In t e r m s of s p a r s i t y t h e wavelet r e p r e s e n t a t i o n has n o t b r o u g h t a n y a d v a n t a g e i n c o m p a r i s o n to t h e full m a t r i x r e p r e s e n t a t i o n . <Ti/jj k, <p> m o m e n t s we see t h a t entries i n a n a r e a where t and t h e kernel is sufficiently H o w e v e r , i f we use wavelets w i t h v a n i s h i n g <T(p, il>j,k> smooth. will be small if ipj,k h a s its s u p p o r t T h e f o l l o w i n g result is p r e s e n t e d i n [beyl91]: Theorem 2.2.2 Assume an integral operator T with kernel s is given which satisfies equation (2.9). Consider e > 0 and wavelets with local support. Then: • only O ( n m log(n )) m of the entries of F$ are larger than e • only 0 ( n ) of the entries of are larger than e m T h e i n t u i t i o n b e h i n d t h i s result is t h a t o n l y those coefficients w h i c h are b a s e d o n wavelets t h a t are close to or o v e r l a p t h e d i s c o n t i n u i t y a l o n g t h e d i a g o n a l x = y w i l l have significant v a l u e . T h e e n u m e r a t i o n o f s u c h wavelets for t h e case o f t h e H a a r M R A is s h o w n i n F i g u r e 2.3 a n d F i g u r e 2.4. F o r t h e s t a n d a r d r e p r e s e n t a t i o n we c a n c o m p u t e a m a t r i x Fs c b y s e t t i n g elements o f Fs w h i c h are s m a l l e r t h a n e to zero. T h i s p r o c e d u r e is c a l l e d k e r n e l c o m p r e s s i o n . W e d e n o t e the operator represented by Dg Fs Ds m e to c o m p u t e a n a p p r o x i m a t i o n o f to 0(n ) m Pu m >m P TP v m m with Ts . c For v £ V t h i s p r o c e d u r e allows us m with only 0 ( n l o g ( n ) ) m m entries i n c o m p a r i s o n entries o f t h e full m a t r i x a p p r o a c h . If u is t h e s o l u t i o n o f t h e i n t e g r a l e q u a t i o n = P e + P TP u m m and s h o w n t h a t || u — u m e u is s o l u t i o n o f Pu m e = P e.+ m Ts u t e t h e n it c a n be easily | | = c „ e , for sufficiently s m a l l e (e.g. [atki76j). T h e s a m e h o l d s for t h e 2 m n o n - s t a n d a r d r e p r e s e n t a t i o n , however, t h e o p e r a t o r here has o n l y 0 ( n ) m T o s u m m a r i z e we have f o u n d a p p r o x i m a t i o n s to t h e m a t r i x entries. F, Dg^Fs^Ds^ ^ j v . m - ^ V t - D j v . m ) , w h i c h allow c o m p u t i n g a m a t r i x v e c t o r p r o d u c t i n 0 ( n l o g ( n ) ) m 18 m (and (and m— w=n 1 jt=n 12 3 a; B. teQ, £=4.5.6,7 (c=2.3 Jc=0 1.2,3 f l tf ¥ F i g u r e 2.3: i,j,k,l ¥ V S t a n d a r d representation, s u c h t h a t the s u p p o r t o f f="—1 j = 2 1=0 V j = 2 1=1 H a a r M R A : supports of h(x,y) = ipi k{x)tpj,i(y} t for o v e r l a p s the d i s c o n t i n u i t y . 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 r e p r e s e n t a t i o n , H a a r M R A : s u p p o r t s j for i,j,k,l s u c h t h a t the s u p p o r t of h overlaps 19 oih(x,y) = xl)i^{x)il)j^{y),i = the discontinuity. 0(n )) operations. For this method to work we need hierarchical basis functions with the m following properties: • Wavelets with vanishing moments. This property is the main requirement necessary 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 h i and g^. This property is necessary to be able to ky compute the decomposition Ds,m and in 0(n ) operations. rn 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 <T(pj i, ipj k> : t and <T(j>j i, (f>j k> for k, t : I : 0< k,l < rij 20 <Tipj i, t ipj k>, t <Tipj i,(f)j k>, t ! can be computed recursively using the dilation equation: rij+i — l <Ti/jj,i,ipj, > = k 52 t O <Tip (p > jh = hk 52 1=0 <T4>j,i,ipj,k> = 52 1=0 nj i-l + 52 9i,Wk,k<T(f) ij,^j i'k> j+ + k=0 52 dijK^T^j+ijAj+i;^ fc=0 52 \i9k,k< <t>j+\j.Aj+i,k> T k=0 T i j + i - l n j + i - i ' <T<t>j,i,<l>j,k> = 52 i,i k,k<T(t> ij,(t)j i;k> h r=o fc=o h J+ + (2-15) The algorithm outlined above in the presented form is not useful for radiosity problems. It reduces the cost for solving the discrete system to 0 ( n ) , but the cost of computing m Fs or F^ still requires 0(n ) e c m 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-presented 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. Clustering Hierarchical methods require computation of the form factors of the operator Pr,T Po, r 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 ) . Smits et. al. [smit94] present a clustering algorithm which reduces 2 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 ) form factors. Them 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(n log(n )) m m form factors. However, it allows computations without the normaliza- tion and denormalization steps (D^ , l m -Djv ) since Ds, >m 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 nonstandard 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(l s ) 2 2 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 functions 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 Aj k t of a basis function fij^ into sets k and A? k . Then (frj^ is replaced by two functions </>]. and ' <p\ which are obtained by restricting the support of 4>j,k to A A ) l k k and A . This newly established basis is not orthogonal anymore resulting in more complicated 2 k 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 d o not c o n s i d e r d i r e c t m e t h o d s i n t h i s section since i n the c o n t e x t of the r a d i o s i t y e q u a t i o n t h e y are c l e a r l y i n f e r i o r to i t e r a t i v e m e t h o d s . T h e a d v a n t a g e o f i t e r a t i v e m e t h o d s here s t e m s less f r o m the desire o f e x p l o i t i n g s p a r s i t y t h a n f r o m the s t r o n g d i a g o n a l d o m i n a n c e of t h e d i s c r e t i z e d s y s t e m , w h i c h i m p l i e s fast convergence o f s i m p l e i t e r a t i v e m e t h o d s . Moreover, the e l e m e n t s of the m a t r i x relate t o a d i s c r e t i z a t i o n o f s m o o t h q u a n t i t i e s , r a i s i n g the p r o s p e c t of efficient m u l t i l e v e l t r e a t m e n t . 2.3.1 Properties C o n s i d e r the l i n e a r s y s t e m of e q u a t i o n s (2.5) w r i t t e n as Mx R = d i a g ( p , • • • , Pn-i) a 0 n = b, where M — I — RG d G as i n t r o d u c e d i n e q u a t i o n (2.6). with A s m e n t i o n e d i n 2.2.1 the size o f M c a n b e c o m e t o o large t o fit i n t o a c o m p u t e r ' s p r i m a r y m e m o r y . M is b l o c k w i s e sparse if t h e scene has o b s t r u c t i o n s s u c h t h a t entire p a t c h e s are not v i s i b l e f r o m each o t h e r . S p a r s i t y is significant i f the scene is c o m p l e x , i.e. i f it is c o m p o s e d of different i s o l a t e d r o o m s . M is g e n e r a l l y not s y m m e t r i c . H o w e v e r , it c a n be easily m a d e s y m m e t r i c by m u l t i p l i c a t i o n w i t h D A R ~ l , where D A = diag(<0n, 0o>) • • • , <4 n-ii > A l s o , b y d e f i n i t i o n two p o i n t s i n a p l a n e are not v i s i b l e f r o m each other, so g k k hence rrikk = l,k := 0 < = 4>n-i>)0 and k < n. N o w we restrict o u r c o n s i d e r a t i o n s to the case o f piecewise c o n s t a n t basis f u n c t i o n s . H e r e the i d e n t i t v du(x) = ^-^^-MdA(y) \\x-y\\ (ui is t h e s o l i d angle) yields Y\,g i k = 1 (e.g. 2 [cohe93b]). S i n c e 0 < pk < 1, M is s t r i c t l y row d i a g o n a l l y 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 t h a t D R~ M l A is n o t o n l y s y m m e t r i c b u t also p o s i t i v e definite p r o p e r t y c a n be u s e d t o solve t h e S P D s y s t e m D R~ Mx l A = D R~ b l A (SPD). This instead. However, this m e t h o d does not allow for v e r y s m a l l or zero pi, w h i c h is a significant d r a w b a c k since s u c h 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 following properties: 1. moderate sparsity; 2. it can be made symmetric by multiplication with the diagonal matrix D^R~ l if p > k 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\ > \ ij\ m 2. rriij < 0,i / j 3. D R- M l A 2.3.2 is SPD (if p > 0,k : 0 < k < n) k 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)" 6 + {I - d i a g ( M ) " M ) x 1 1 26 (fc) In case of the r a d i o s i t y e q u a t i o n d i a g ( M ) = /. Recalling F = RG we o b t a i n the simpler form: x (k i) + = T h i s f o r m is also k n o w n as one step o f b + .(k) (2.16) F? Picard iteration (e.g. [hack95]). P i c a r d i t e r a t i o n s i m p l y a p p r o x i m a t e s the N e u m a n n series f r o m e q u a t i o n (2.3). S o each sweep c o r r e s p o n d s to one a d d i t i o n a l level o f reflection. l i g h t - s o u r c e after k reflections. If we choose xj°) = b then F b k A l s o note f r o m e q u a t i o n (2.16) t h a t one step o f the i t e r a t i o n is e q u i v a l e n t to the a d d i t i o n o f the c u r r e n t r e s i d u a l ^(fc+l) {k) = x is the c o n t r i b u t i o n of the = b — Mx^ Picard resulting in (k)_ + r Gauss-Seidel Iteration L e t / c o u n t single u p d a t e s r a t h e r t h a n sweeps. T h e n Gauss-Seidel iteration updates the u n k n o w n s X{ i n the f o l l o w i n g way: (2.17) The i n d e x i(l) (O'ew = specifies w h i c h v a r i a b l e is u p d a t e d (0,1, 2 , . . . ,n — 1, 0 , 1 , 2 , . . . ,n — 1,...). i n the l-th step. T h e u s u a l o r d e r is: T o i n c o r p o r a t e the r e s i d u a l e q u a t i o n (2.17) c a n be easily r e w r i t t e n y i e l d i n g the f o l l o w i n g a l g o r i t h m : A l g o r i t h m 2.3.1 Gauss-Seidel U p d a t e ( M , 6, i:= next i in turn Vj : j / i do: 27 i,x,r) explicitly ri :.= 0 return (i, x, r) If we choose xj ^ = 0 as the initial guess we can compute the initial residual: r^ = b. 0 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^) solution at any give stage . The quantity x — where x ^ ' is the current approximate 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 final step corresponds to exactly one complete sweep of the Jacobi iteration based on For this sweep no new entries of M have to be computed. 28 This xj \ k W e r e w r i t e the a l g o r i t h m m a i n t a i n i n g the i n v a r i a n t x^ = x^ + (x^ i t e r a t i o n v a r i a b l e i n P R ) . T h i s f o r m of the P R a l g o r i t h m gives it a n i m m e d i a t e interpretation. T h e a p p r o x i m a t i o n x^ d i s p l a y p u r p o s e s after each Algorithm 2.3.2 is m o r e a c c u r a t e the physical W e c a n use x^ k for update. Progressive Refinement choose i such that ri is maximized Vy / t h a n x^ \ is Update(M, b, x, r) • / do: Tj . Xi J . Tj Xj r e t u r n (i,x, TTljiTi TThJ7T'1 r) N o t e t h a t rj^ is not the r e s i d u a l r e g a r d i n g x^ the S o u t h w e l l i t e r a t i o n . S o as i n i t i a l guess x^ b u t r a t h e r the r e s i d u a l of xj ^ k we h a v e to choose x^ := xj ^ 0 from — b. + A} 0 S o i n c o n t r a s t to the G a u s s - S e i d e l m e t h o d the r e s i d u a l of the a p p r o x i m a t e s o l u t i o n is not available in this m e t h o d . T h e p h y s i c a l i n t e r p r e t a t i o n is i n s o m e sense e x a c t l y o p p o s e d to the one of the G a U s s Seidel m e t h o d . (shooting). In P R u n s h o t r a d i o s i t y f r o m one p a t c h i is d i s t r i b u t e d to a l l o t h e r p a t c h e s j C h o o s i n g the m a x i m u m r e s i d u a l (sorting), w h i c h is i n i t i a l l y c o n c e n t r a t e d i n the light-sources, g u a r a n t e e s t h a t significant light t r a n s p o r t s are c o n d u c t e d e a r l y i n the i t e r a t i o n . PR was i n t r o d u c e d b a s e d o n p h y s i c a l a r g u m e n t s i n [coheSS] . 1 [gort94] s h o w e d a b o v e r e l a t i o n to G a u s s - S e i d e l i t e r a t i o n . Hhe original version of progressive refinement chooses an i such that |>li||r;| is maximized 29 the method PR Shao Gortler Feda reference [cohe88] [shao93] [gort94] [feda92] [xu94] Xu [xu94] selection of i M DAL I DAL I DAL £ j i = 1 +F = I « I + F F (I - F)D L A (I - +F *• F)D r_ A I + F + = ZM** Pavg 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. Overrelaxation Methods 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 x := x + t { A Vy' d o : Tj := Tj — TUjiA return (i, x, L) Special cases of this algorithm are the Gauss-Seidel method and progressive refinement, 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 choose 2 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~ by the leading components of the Neumann series. l ei is the i-th vector of unity 2 30 The quantities in Table 2.1 are defined as follows: p is a matrix with all entries equal to 1) and fikl = Fi ^ = Pkfki • k I and 0 : otherwise = YliPil i n avg F = T\\R1DA (1 (fiki)nxn' (k = i or I = i) \ 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 corresponds 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 overrelaxation 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 approximate 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 r iy 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 ambient . 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 S P D system MM T be computed as x = M y. T and solves for y. The solution of the original system can then 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 S P D 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 S P D system: D R~ Mx l A = D R~ b l A 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 (p avg = 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 p avg > 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 (fc) ) by refinement of F - ~^ (x} ~^) e( k k = solution of (T- F <*>(xW))£( > = 6 e fc+1 : 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 ) entries of the matrix before starting the solution process. Such a procedure 2 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 samp n < 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 i n t e g r a n d to [zatz93] a n d [schr93]. 2.4.1 Form Factor Approximation Piecewise Constant Basis Functions We discuss the m o s t c o m m o n l y u s e d t e c h n i q u e s here. O t h e r suggested m e t h o d s are a t r a n s - f o r m a t i o n of the surface i n t e g r a l i n a c o n t o u r i n t e g r a l (contour integral method), analytic c o m p u t a t i o n for u n o c c l u d e d scenes b y a p p r o x i m a t i o n w i t h s t a n d a r d g e o m e t r i e s , a n d M o n t e C a r l o integration. F o r piecewise c o n s t a n t basis f u n c t i o n s we have to c o m p u t e the i n t e g r a l : 1 9ki = ff J-T-7 / / \A-k\ 1 JjA ff // N cos L(n ( Pv, ^q- — p ) p) c ocos s Z L(n ( n q,,p p -— 4 )q) ( C O s Z v (p,q) 2 JJA, k ff " 7T|| p - dA^qjdA^p) '2 q II ff cos L(n ,q - p) " vs(p,q) ^ ~ du(q)dA (p) x X PJ k k 9ki k 9 ~ s If we choose p to be the center of p a t c h A If p a t c h e s A n ff = JJ is: k cos s(P'Q) v n t h e n a s i m p l e a p p r o x i m a t i o n for g i L(n , q — p) * " -dcu(q) x (2.18) 1 a n d Ai are s m a l l , far a p a r t a n d there are no o c c l u d e r s close t o A k a p p r o x i m a t i o n is j u s t i f i e d . then this S t i l l , one d o u b l e i n t e g r a l i n v o l v i n g the v i s i b i l i t y f u n c t i o n vs has to be c o m p u t e d . T h i s c a n be d o n e b y s i m p l e n u m e r i c a l i n t e g r a t i o n or b y one o f t h e m e t h o d s described below. Hemicube T h i s technique subdivides il that the part of p a t c h A s u b t e n d e d b y Qi is c o m p l e t e l y v i s i b l e or i n v i s i b l e f r o m p t h e n t i n t o a fixed n u m b e r n herni o f subsets fij 0 < i n t e g r a l of e q u a t i o n (2.18) is i n d e p e n d e n t o f t h e l o c a t i o n of A t 36 i < n . hemi Assuming the itself. If the fij are sufficiently s m a l l t h e n t h i s a s s u m p t i o n a p p r o x i m a t e l y h o l d s . T h e m e t h o d c o m p u t e s a c o m p l e t e row o f f o r m factors f 0 kh < I < n. p a t c h a n d a d i s t a n c e di. W i t h each Qi we store a coefficient 9hemi,i, a n index k x of a I n i t i a l l y di is set t o oo. T h e a l g o r i t h m is s h o w n 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 9hemi,i 9hemi,i :« // q - p) dco(q) is computed based on just one sample oj the kernel. f o r all I: 9ki: =0 f o r all 0 < i < nh i : em 9kki • ~ 9kki ~t~ 9hemi,i T h i s m e t h o d has obvious disadvantages. F a r a w a y or s m a l l p a t c h e s m i g h t n o t be p r o p e r l y r e p r e s e n t e d b y the d i s c r e t i z a t i o n o f Q a n d r e s u l t i n a l i a s i n g i n the o b t a i n e d s o l u tion. O n the o t h e r h a n d the m e t h o d c a n use e x i s t i n g Z - b u f f e r h a r d w a r e . i n t r o d u c e d b y [cohe85]. T h e name hemicube is b a s e d o n t h e fact t h a t T h e m e t h o d was Qi are t h a t t h e y c o r r e s p o n d to r e c t a n g l e s o n the s u r f a c e o f a c u b e w i t h center i n p. 37 chosen such Higher-Order basis functions and Hierarchical Methods T h e r e are two m a j o r differences b e t w e e n the q u a d r a t u r e for piecewise c o n s t a n t basis f u n c t i o n s a n d h i g h e r - o r d e r basis f u n c t i o n s a n d h i e r a r c h i c a l m e t h o d s . F i r s t , the patches f o u n d in higher- o r d e r m e t h o d s a n d o n coarse levels o f h i e r a r c h i c a l m e t h o d s are p o t e n t i a l l y bigger t h a n those f o u n d w i t h piecewise c o n s t a n t basis f u n c t i o n s . Secondly, simple approximations assuming c o n s t a n t r a d i o s i t y over a p a t c h w o u l d c l e a r l y defeat the p u r p o s e of h i g h e r - o r d e r m e t h o d s . visibility-in-quadrature W i l l m o t t a n d H e c k b e r t [will97b] use and fractional visibility to c o m p u t e the f o r m factors. V i s i b i l i t y - i n - q u a d r a t u r e a p p r o x i m a t e s t h e f o r m f a c t o r i n t e g r a l b y n u m e r i c a l q u a d r a t u r e . V a r i o u s q u a d r a t u r e rules are d i s c u s s e d i n [gers95]. T h e s e m e t h o d s are e x p e n s i v e since each s a m p l e requires e v a l u a t i o n o f the v i s i b i l i t y f u n c t i o n vsv i s i b i l i t y is the m e t h o d i n t r o d u c e d i n [hanr91]. B o t h patches A k c o m p u t a t i o n of the f o r m f a c t o r are s u b d i v i d e d i n t o a 4 x 4 g r i d . and A t involved in k to a g r i d p o i n t o n Ai. s a m p l e s we o b t a i n a n e s t i m a t e a o f the f r a c t i o n o f A k f a c t o r is t h e n c o m p u t e d w i t h o u t v i s i b i l i t y f u n c t i o n multiplied by a. the N e x t , for each g r i d p o i n t of p a t c h Ak the v i s i b i l i t y f u n c t i o n w i t h o n l y one g r i d p o i n t of Ai is c o m p u t e d . a fixed s c h e m e to a l l o c a t e g r i d p o i n t s o n A Fractional [hanr91] use F r o m these v i s i b i l i t y t h a t is v i s i b l e f r o m A\. (unoccluded formfactor) T h e form and F r a c t i o n a l v i s i b i l i t y reduces the n u m b e r of s a m p l e s r e q u i r e d f r o m n afterwards A samp to n 2 samp' In h i e r a r c h i c a l m e t h o d s we c a n c o m p u t e f o r m f a c t o r s o n coarse levels a c c u r a t e l y b y u s i n g the d i l a t i o n e q u a t i o n a n d base the c o m p u t a t i o n o n t h e f o r m factors c o m p u t e d for the n e x t finer level as g i v e n b y e q u a t i o n (2.15). 38 2.4.2 Form Factor Estimation T o be a b l e to realize the t o p - d o w n a p p r o a c h t a k e n b y h i e r a r c h i c a l m e t h o d s we need a f u n c t i o n t h a t c a n e s t i m a t e the m a g n i t u d e or a n u p p e r b o u n d , o f a f o r m f a c t o r i n v o l v i n g wavelets w i t h o u t h a v i n g to e x p e n s i v e l y c o m p u t e it. T h e i r c o m p u t a t i o n o n coarse levels is p r o h i b i t i v e l y e x p e n s i v e since t h e y w o u l d r e q u i r e a large n u m b e r of k e r n e l s a m p l e s to be able to represent the kernel accurately. F o r piecewise c o n s t a n t basis f u n c t i o n s [hanr91] give a n u p p e r b o u n d for the u n o c c l u d e d f o r m f a c t o r (this c o r r e s p o n d s to the f o r m factor i n v o l v i n g o n l y s c a l i n g f u n c t i o n s ) . V i s i b i l i t y is incorporated using fractional visibility yielding an estimate the v a l u e of f o r m factors i n v o l v i n g the H a a r wavelets, f i k if A , ^ is o n l y p a r t i a l l y v i s i b l e f r o m ' ' //c/ for <T(pj i, (j)j k>t t T o estimate is m u l t i p l i e d b y a p a r a m e t e r £ > 1 Ajj. In t h e g e n e r a l case o f h i g h e r - o r d e r basis f u n c t i o n s there is no m e t h o d k n o w n to b o u n d the u n o c c l u d e d f o r m f a c t o r . T h e r e f o r e we have to t a k e a different a p p r o a c h . L e t the n u m b e r of v a n i s h i n g m o m e n t s o f the c o n s i d e r e d wavelets be v. T h e n we k n o w t h a t a f o r m f a c t o r is s m a l l if the kernel b e h a v e s s i m i l a r to a p o l y n o m i a l o f degree v — 1. [gort93b] suggest to e s t i m a t e the f o r m f a c t o r b y e s t i m a t i n g the d e v i a t i o n o f the k e r n e l f r o m a p o l y n o m i a l o f degree v — 1. T h e k e r n e l is s a m p l e d a n d a n i n t e r p o l a t i n g p o l y n o m i a l is c o n s t r u c t e d o n these s a m p l e s . T h e d e v i a t i o n is c o m p u t e d b y c o m p a r i n g t h e p o l y n o m i a l to based samples c o m p u t e d i n b e t w e e n the p r e v i o u s l y c o m p u t e d s a m p l e s . T h i s gives a q u a l i t a t i v e m e a s u r e for t h e m a g n i t u d e o f the f o r m f a c t o r . 2.4.3 Comparison F o r m f a c t o r c o m p u t a t i o n a n d e s t i m a t i o n for h i g h e r - o r d e r basis f u n c t i o n s are m o r e e x p e n sive t h a n for piecewise c o n s t a n t basis f u n c t i o n s . T h e c h e a p e r m e t h o d o f f r a c t i o n a l v i s i b i l i t y 39 p r o d u c e s i n a c c u r a t e results i n t h i s case. Q u a d r a t u r e f o r m u l a e need h i g h e r a c c u r a c y , conse- quently a higher n u m b e r of visibility samples. F u r t h e r m o r e , hardware-supported^ methods like the h e m i c u b e m e t h o d are n o t a v a i l a b l e here. 2.5 Summary progressive radiosity. T h e m e t h o d m o s t c o m m o n l y u s e d i n c o m p u t e r g r a p h i c s is This method c o m b i n e s a d i s c r e t i z a t i o n w i t h piecewise c o n s t a n t basis f u n c t i o n s , s u b s t r u c t u r i n g a n d p r o gressive refinement as the s o l u t i o n m e t h o d . F o r m factors are t y p i c a l l y c o m p u t e d using, the h e m i c u b e a p p r o a c h . T h e h e m i c u b e m e t h o d c a n c o m p u t e c o m p l e t e rows or c o l u m n s o f f o r m factors as n e e d e d b y progressive refinement. N u m e r o u s ways of p a r a l l e l i m p l e m e n t a t i o n s o n different p a r a l l e l a r c h i t e c t u r e s have b e e n p r e s e n t e d [guit91], [garc97]). execution time. 2 [chen89], [bric95], [baum90], T h e s e m e t h o d s e x p l o i t t h a t f o r m factors c a n be c o m p u t e d i n d e p e n d e n t l y . P r o g r e s s i v e r a d i o s i t y requires o n l y O(n) 0(n ) ([reck90], m e m o r y (n is the n u m b e r o f p a t c h e s ) , b u t requires T h e h i e r a r c h i c a l m e t h o d w i t h piecewise c o n s t a n t basis f u n c t i o n s is d i s c u s s e d f r e q u e n t l y i n the l i t e r a t u r e . It requires O(n) m e m o r y a n d O(n) execution time. H o w e v e r , the c o n s t a n t factors i n these t i m e a n d m e m o r y e s t i m a t e s are m u c h bigger t h a n t h o s e for progressive r a d i o s i t y . A l l c o n s i d e r e d m e t h o d s t r y to r e d u c e t h e s y s t e m o f e q u a t i o n s (2.5) be c o m p u t e d a n d s o l v e d efficiently. s u c h t h a t it c a n O n t h e one h a n d this is d o n e b y e x p l o i t i n g s m o o t h n e s s p r o p e r t i e s of the s o l u t i o n . D i s c o n t i n u i t y m e s h i n g a n d h i g h e r - o r d e r G a l e r k i n m e t h o d s choose a basis w h i c h a p p r o x i m a t e s t h e s o l u t i o n well. O n the other h a n d hierarchical methods a n d s u b s t r u c t u r i n g use s m o o t h n e s s p r o p e r t i e s o f the k e r n e l . T h e u l t i m a t e g o a l is to find a sparse r e p r e s e n t a t i o n f u l l y a u t o m a t i c a l l y w i t h o u t user i n t e r a c t i o n b y a n a d a p t i v e process. T h i s process c a n be b a s e d p u r e l y o n t h e k e r n e l itself, or i n a n i t e r a t i v e process it c a n be 40 b a s e d o n the a p p r o x i m a t e s o l u t i o n f o u n d so far. H i e r a r c h i c a l m e t h o d s r e d u c e t h e n u m b e r of p a r a m e t e r s to choose to the c o m p r e s s i o n p a r a m e t e r e a n d the n u m b e r o f levels. E f f o r t s have b e e n m a d e to i n t e g r a t e c o n c e p t s f r o m d i s c o n t i n u i t y m e s h i n g w i t h h i g h e r order hierarchical methods. H o w e v e r , n o c o m p l e t e a n a l y s i s o f s u c h a s y s t e m has b e e n p r e - sented. C l u s t e r i n g has b e e n a p p l i e d to h i e r a r c h i c a l m e t h o d s functions. W e do not know of reports o n the w i t h piecewise c o n s t a n t basis a p p l i c a t i o n to h i e r a r c h i c a l m e t h o d s with h i g h e r - o r d e r basis f u n c t i o n s . F o r piecewise c o n s t a n t basis f u n c t i o n s n u m e r o u s p r o v a b l y convergent s o l u t i o n m e t h o d s exist. T h e s e s o l u t i o n m e t h o d s are also a p p l i e d to s y s t e m s b a s e d o n h i g h e r - o r d e r basis f u n c t i o n s a n d the c o m p r e s s e d s y s t e m s a r i s i n g i n h i e r a r c h i c a l m e t h o d s . k n o w of a n y c a r e f u l m a t h e m a t i c a l a n a l y s i s o f t h e i r convergence. 41 H o w e v e r , we d o n o t 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 ( M R A ) 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 : <ipi k, "4>j,i> — t ®,<'4>i,k,' r i,k> = !)• We do not consider flatlets here (flatlets [gort93b] are based on the Haar l , 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 = /ii = l/y/2: (3.1) (=0,1 The Haar wavelet is defined based on orthogonality considerations with g = l / \ / 2 , g i = 0 (3.2) /=o,i 43 Haar MW 0 F i g u r e 3.1: 1 2 3 4 5 6 7 D e p e n d e n c i e s b e t w e e n different levels o f t h e M R A . T h e n o d e s i n d i c a t e s c a l i n g f u n c t i o n s a n d wavelets. W e define "00,0 = Tp(x)- the s c a l i n g f u n c t i o n a n d wavelet o n the coarsest level $ ,o 0 (x) = 4>{x), T h e h i e r a r c h y of basis f u n c t i o n s is t h e n o b t a i n e d b y the d i l a t i o n e q u a t i o n , e q u a t i o n (2.10), w h e r e be define h kfik •= h , h k+\ •= h 0 kfi u g ^k •= 9o a n d g , k+i •= 9i k k 2 ( c o r r e s p o n d i n g t o F i g u r e 3.1). F u n c t i o n s c o n s t a n t o n the entire d o m a i n c a n be r e p r e s e n t e d e x a c t l y i n t h e H a a r basis. T h e H a a r wavelets have one v a n i s h i n g m o m e n t . W e show t h e g r a p h of the H a a r s c a l i n g f u n c t i o n a n d wavelet i n F i g u r e 3.2. 3.1.2 Multiwavelets T h e M u l t i w a v e l e t M R A was p r e s e n t e d b y [alpe93]. It is b a s e d o n L e g e n d r e p o l y n o m i a l s . H e r e we are c o n c e r n e d w i t h m u l t i w a v e l e t s b a s e d o n l i n e a r f u n c t i o n s . expressions for these f u n c t i o n s ( F i g u r e 3.3). W e c a n give e x p l i c i t T h e h i e r a r c h y is c o n s t r u c t e d a n a l o g o u s l y 44 to Scaling Function Wavelet phi psi 0 O -1 F i g u r e 3.2: H a a r s c a l i n g f u n c t i o n a n d wavelet. the case o f the H a a r wavelets u s i n g the d i l a t i o n e q u a t i o n . T h e wavelets have two v a n i s h i n g moments. They are a b l e to r e p r o d u c e l i n e a r f u n c t i o n s exactly. T h e space spanned by multiwavelet scaling functions contains discontinuous functions. T h e coefficients for the g e n e r a t i o n o f the h i e r a r c h y are l i s t e d i n the f o l l o w i n g table: MW k Scaling Functions MW Wavelets I 'hki 0 0 0 1/V2 1 l/S/2 3 l/y/2 0 l/y/8 2 0 -\/3/y/8 -l/y/2 1 1/V8 1 V3/VS 2 3/78 2 -1/V8 3 l/y/8 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 wavelet. 4 The C D V scaling functions and wavelets in the interior of the interval, i.e. the functions with indices k : 2 < k < n m On the boundaries, — 2, are identical to the D scaling function and wavelet. 4 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 -2,V'n -im m 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: h = (1 + v 3)/(4 2), / / 0 h = (3 + v 3)/(4v 2), h = (3 / / x 92 / 2 v v 3 ) / ( 4 v 2 ) , h = (1 - v 3)/(4v 2), So = -h , / / / 3 3 = -hi, g = h . 3 Q For the functions at the boundaries we have the coefficients: C D V Scaling Functions k k - n - I- n - 1 h I hk, 0 6.033325e-01 -0 8.705088e - 01 1 6.908955e-01 -1 4.348970e - 01 -2 2.303890e - 01 m 2 -3.983130e - 01 1 m k -1 0 3.751746e-02 -o. — 1.942334e - 01 1 4.573277e-01 -1' 1.901514e - 01 2 8.500881e - 01 -2. 3.749553e-01 3 2.238204e - 01 . -3 7.675567e - 01 4 -1.292227e- 01 -4 4.431490e - 01 47 g = h, x 2 C D V Wavelets k I k- n - 1 i - n 9k,i m - l g ktl -0 -2.575129195e - 01 1 5.463927140e - 01 -1 8.014229620e-01 2 -2.587922483e - Ol -2 -5.398225007e - 01 0 -7.965435169e-01 -0 0 1.003722456e-02 -0 3.717189665e - 01 1 1.223510431e-01 -1 -3.639069596e - 01 2. 2.274281117e - 01 -2 -7.175799994d- 01 3 -3 4.010695194e - 01 -4 2.315575950d- 01 -8.366029212e - 01 4 4.830129218e - 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 approximation 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 = 48 (3.3) 49 f M W G D \ / Figure 3.5: Projection of the function / of (3.3) into different spaces. A l l bases are w i t h 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. In Figure 3.6 we show results for an equivalent experiment w i t h a function: 3(x + y) < 4 fxy f(x,y) (3.4) = 0 otherwise Here is the table for the compression rate of the two-dimensional case. We show the proportion of coefficients of wavelets that were below a given threshold e: Haar MW CDV -l 100% 100% 99% 10" -2 96% 96% 88% 10"-3 73% 78% 71% 10" -4 38% 68% 52% 10"-5 24% 56% 33% -6 16% 42% 25% 12 16% 13% 9% e 10" 10" 10" 1 50 4 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 L error obtained for the function / defined by (3.4) to the x relative L error obtained for the function x g(x, y) = sjxy: Haar MW CDV 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 discontinuity. This results in larger coefficients of such basis functions. In Chapter 4 we will experimentally analyze how these three types of wavelet bases 51 p e r f o r m i n h i e r a r c h i c a l m e t h o d s for the s o l u t i o n o f the r a d i o s i t y e q u a t i o n , w h e r e the process involves s o l v i n g a n i n t e g r a l e q u a t i o n , n o t j u s t f u n c t i o n a p p r o x i m a t i o n . 3.2 Solution Methods for Wavelet Radiosity A s d i s c u s s e d i n S e c t i o n 2.3.3 K r y l o v subspace m e t h o d s t e c h n i c a l l y present a n a l t e r n a t i v e t o the P i c a r d i t e r a t i o n for W a v e l e t R a d i o s i t y . In m a n y a p p l i c a t i o n s t h e y p e r f o r m b e t t e r than Picard's. T h e matrix system i n h i e r a r c h i c a l m e t h o d s w i l l g e n e r a l l y n o t be S P D . In p a r t i c u l a r , for h i g h e r - o r d e r basis f u n c t i o n s we d o n o t have a t h e o r e m w h i c h w o u l d a l l o w us to m a k e t h e s y s t e m our analysis. radiosity. S P D . T h e r e f o r e we choose the m e t h o d s T by and G M R E S for W e w i l l a n a l y z e the p e r f o r m a n c e of these m e t h o d s for the case o f wavelet In o u r i m p l e m e n t a t i o n o f t h e m e t h o d s we follow [saad95], m a t r i x vector m u l t i p l i c a t i o n Mx CGNR Mj^ DN x, m e Mx by M D x m Nc N>m the a n d the m a t r i x v e c t o r m u l t i p l i c a t i o n as d e s c r i b e d i n S e c t i o n 2.2.3. tm T h e i n v e s t i g a t i o n o f fast s o l u t i o n m e t h o d s wavelet r a d i o s i t y is of interest. entire sparse s y s t e m where we r e p l a c e Mry e for the l i n e a r a l g e b r a s y s t e m i n case of O n t h e one h a n d , wavelet r a d i o s i t y requires k e e p i n g i n m e m o r y u n l i k e progressive a n d c o m p u t a t i o n o f f o r m factors are i n t e r m i x e d . r a d i o s i t y where s o l u t i o n the process S o the s o l u t i o n process is a n a d d i t i o n a l c o m p u t a t i o n a l step, the cost o f w h i c h seems t o be w o r t h w h i l e to m i n i m i z e . O n the other h a n d , i n c o n n e c t i o n w i t h t h e m u l t i l e v e l - l i k e m e t h o d s f r o m S e c t i o n 2.3.3 a r e d u c t i o n o f t h e n u m b e r o f i t e r a t i o n s c o u l d also r e d u c e the n u m b e r o f c o s t l y refinement steps n e e d e d . 3.3 Implementation Aspects H i e r a r c h i c a l m e t h o d s for r a d i o s i t y i n t r o d u c e the f o l l o w i n g t y p e s o f a p p r o x i m a t i o n e r r o r i n t o the c o m p u t a t i o n : 52 • d i s c r e t i z a t i o n e r r o r c a u s e d b y the G a l e r k i n d i s c r e t i z a t i o n . • discretization error caused by compressing the G a l e r k i n discretization. • error caused by the solution m e t h o d . • e r r o r c a u s e d b y the f o r m factor a p p r o x i m a t i o n . • error c a u s e d b y a f a u l t y e s t i m a t i o n of f o r m factors. In o u r a n a l y s i s we w a n t to focus o n the first three t y p e s o f error. W e therefore d e s i g n e d o u r a l g o r i t h m a n d test e n v i r o n m e n t s to m i n i m i z e the last two t y p e s o f error. W e d o n o t use a n y o f the m u l t i l e v e l - l i k e m e t h o d s since this w o u l d m a k e t h e i n d e p e n d e n t a n a l y s i s o f s o l u t i o n m e t h o d s a n d d i s c r e t i z a t i o n difficult. In the f o l l o w i n g we d e s c r i b e h o w we r e a l i z e d t h e v a r i o u s c o m p o n e n t s o f t h e a l g o r i t h m for wavelet r a d i o s i t y i n the n o n - s t a n d a r d r e p r e s e n t a t i o n . 3.3.1 Galerkin Discretization In o u r e x p e r i m e n t s separate we use H a a r wavelets, M R A for each surface Si. multiwavelets a n d CDV wavelets. W e use H a a r wavelets have one v a n i s h i n g m o m e n t a b i l i t y to a p p r o x i m a t e f u n c t i o n s c o n s t a n t o n t h e surface exactly. and a the Multiwavelets and C D V wavelets have b o t h two v a n i s h i n g m o m e n t s a n d c a n a p p r o x i m a t e a l l f u n c t i o n s l i n e a r o n the surface exactly. A l l the wavelets a n d s c a l i n g f u n c t i o n s were m a d e o r t h o n o r m a l i n 1/2(5) b y multiplication with y^rjq ( s is the s u p p o r t o f X^/clo* ^o,k, s = 1 for H a a r a n d m u l t i w a v e l e t s , s = 4 for C D V wavelets). 3.3.2 Form Factor Estimate T h e f o r m f a c t o r is e s t i m a t e d as d e s c r i b e d i n S e c t i o n 2.4.2 for h i g h e r - o r d e r b a s i s f u n c t i o n s . F o r each o f the f o u r d i m e n s i o n s we c o m p u t e 53 the average o f the d e v i a t i o n over t h e three remaining dimensions. T h e e s t i m a t e is t h e m a x i m u m o f these four averages. To compute d e v i a t i o n f r o m a c o n s t a n t f u n c t i o n (for a basis w i t h one v a n i s h i n g m o m e n t ) we c o m p u t e difference o f the kernel values at the b o u n d a r i e s o f the i n t e r v a l . the F o r the d e v i a t i o n f r o m a l i n e a r f u n c t i o n (for bases w i t h two v a n i s h i n g m o m e n t s ) we i n t e r p o l a t e b y a l i n e a r f u n c t i o n b a s e d o n t h e kernel values at the b o u n d a r i e s o f t h e i n t e r v a l a n d c o m p u t e the difference of t h e v a l u e at the center of t h e i n t e r v a l a n d the kernel v a l u e at t h e center o f the i n t e r v a l . 3.3.3 We Form Factor Computation d o n o t d e a l w i t h the p r o b l e m o f s i n g u l a r i t i e s . W e design our experiments to exclude singular f o r m factor integrals. Unlike [gort93b] a n d [will97b] i n m o s t of t h e i r e x p e r i m e n t s , q u a d r a t u r e for t h e c o m p u t a t i o n o f f o r m factors. we use visibility-in- W e t h i n k t h a t i n the c o n t e x t o f h i g h e r o r d e r basis f u n c t i o n s it is p a r t i c u l a r l y i m p o r t a n t to have a g o o d a p p r o x i m a t i o n to the f o r m factors a l o n g d i s c o n t i n u i t i e s since we have larger p a t c h e s t h a n w i t h piecewise c o n s t a n t function. basis F r a c t i o n a l v i s i b i l i t y does n o t allow to e x p l o i t t h e b e t t e r a p p r o x i m a t i o n o f h i g h e r o r d e r basis f u n c t i o n s along discontinuities since it e s s e n t i a l l y assumes c o n s t a n t visibility w i t h i n ;a p a t c h . W e w i l l s u p p o r t t h i s c l a i m e x p e r i m e n t a l l y i n c h a p t e r 4. We i l l u s t r a t e the d i m e n s i o n a l case J p r o c e d u r e o f h o w we a p p r o x i m a t e f(x)w(x)dx, where / t h e f o r m factors c o r r e s p o n d s to t h e k e r n e l , w for t h e one- to a basis f u n c t i o n . D u e to t h e i r r e g u l a r i t y of t h e C D V basis f u n c t i o n s a s i m p l e r u l e o f s a m p l i n g t h e i n t e g r a n d f(x)g(x) is not a p p l i c a b l e i n o u r case. W e therefore approximate / f u n c t i o n , a n d c o m p u t e t h e i n t e g r a l s i n v o l v i n g w separately. b y a piecewise l i n e a r If we a p p r o x i m a t e / di + biX i f x : Xi < x < Xi \ t h e n t h e i n t e g r a l to c o m p u t e is: + fw(x) = 52 a i / w(x)dx + bi 54 xw(x)dx b y f(x) = 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 n — 3. The Haar and multiwavelet coefficients can be computed samp 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].).. 51! fxo Haar • MW D 4 x ) ^ x Ixl x^Wte 0 1 2 I I 2 8 fro 1 2 1 2 1 8 HI xfr(x)dx 3 8 3 8 1^3 01 CDV 0 ( 2 4 ^ 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 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 Representation of the Sparse 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 i n the m a t r i x was d e t e r m i n e d to be n o n - z e r o . N o t e t h a t the h i e r a r c h y o f the C D V M R A is n o t a tree s t r u c t u r e . T h e r e f o r e we d e s i g n e d the s t r u c t u r e s u c h t h a t it is p o s s i b l e to check w h e t h e r a l i n k h a d a l r e a d y b e e n c o m p u t e d to prevent r e d u n d a n t c o m p u t a t i o n s d u r i n g the t o p - d o w n form factor c o m p u t a t i o n . 3.3.5 Solution M e t h o d Unless otherwise n o t e d we use P i c a r d i t e r a t i o n w i t h 80 i t e r a t i o n s for t h e s o l u t i o n o f the system of equations. and P i c a r d i t e r a t i o n is t h e s t a n d a r d m e t h o d used for h i e r a r c h i c a l m e t h o d s it c o n v e r g e d s t a b l y i n a l l o u r e x p e r i m e n t s . 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. p = 0, and constant emission e\ = 100. Note the proximity of the light-source to the surface S , which causes a steep radiosity gradient. x 0 „ . - - - 0.4 1 Figure 4.2: "Occluded". Surface S 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 S has zero reflectivity and has no emission. 0 2 58 4.1.1 Test Scenes T h e scene " U n o c c l u d e d " is the s a m e scene t h a t was used b y [gort93b] for t h e i r m e a s u r e m e n t s . It allows us to c o m p u t e a n a n a l y t i c s o l u t i o n . T h e p r o x i m i t y o f e m i t t e r a n d receiver in a h i g h b u t s m o o t h r a d i o s i t y g r a d i e n t . w i t h the b a s i s - f u n c t i o n s . T h e l i g h t - s o u r c e it r o t a t e d to p r e v e n t results alignment T h e r a d i o s i t y o n face <Si is e q u a l to t h e e m i s s i o n e, therefore, we get t h e f o l l o w i n g a n a l y t i c e x p r e s s i o n for the r a d i o s i t y : u(p) = p(p) r (p,q)e(q)dSi(q) 2 JSi F o r o u r reference 128 x 128 s o l u t i o n we a p p r o x i m a t e d t h i s i n t e g r a l b y a m i d p o i n t r u l e w i t h samples. T h e scene " O c c l u d e d " has a b l o c k e r t h a t causes a s h a d o w o n SQ. L i k e i n the u n o c c l u d e d scene we r o t a t e d the b l o c k e r to p r e v e n t a l i g n m e n t o f t h e s h a d o w the supports of the basis-functions. boundaries with W e d o n o t have a n a n a l y t i c s o l u t i o n for this scene so for q u a n t i t a t i v e results we c o m p a r e w i t h a m a s t e r s o l u t i o n o b t a i n e d w i t h a d i s c r e t i z a t i o n of h i g h g r a n u l a r i t y . T h e m a s t e r s o l u t i o n was c o m p u t e d u s i n g m u l t i w a v e l e t s w i t h e = I O with m= 5 levels ( r e s u l t i n g i n n$ = 4096 b a s i s - f u n c t i o n s ) and n samp - 4 , = 3. In b o t h scenes a l l the i n t e r e s t i n g events take p l a c e o n t h e surface SQ. T h e r e f o r e , i n the subsequent d i s c u s s i o n we w i l l o n l y d i s p l a y the c o m p u t e d i m a g e s for t h i s surface. Errors are also o n l y c o m p u t e d for SQ. 4.1.2 Form Factor Computation T o s u p p o r t the c l a i m m a d e i n s u b s e c t i o n 3.3.3 t h a t the m e t h o d o f f r a c t i o n a l v i s i b i l i t y is n o t s u i t a b l e for h i g h e r o r d e r b a s i s - f u n c t i o n s b o t h techniques. we c o m p u t e d a s o l u t i o n w i t h m u l t i w a v e l e t s using W e show the results i n F i g u r e 4.3. T h e y s h o w t h a t f r a c t i o n a l v i s i b i l i t y fails to a d e q u a t e l y represent s h a d o w b o u n d a r i e s . A s i m i l a r o b s e r v a t i o n was m a d e b y [will97b]. 59 Visibility-ln-Quadrature Fractional Visibility Figure 4.3: Result for the scene "Occluded" obtained using a basis of multiwavelets with n = 1024 functions. The left image was obtained using visibility-in-quadrature, the right image using fractional visibility. 4 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 Unoccluded Figure 4.4: "Unoccluded". Relative L solution. x error. Obtained by comparison to an analytical our computations: e Haar IO MWi lCT CDV 4 - 5 10~ n 777. 4 4 3 5 n 5 = 1024 3 4 n 4 = 1024 3 3 ™3 = 1024 Following [gort93b] we calculate the relative error in the L norm. These errors are x 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 basis61 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 CDV 20000 MW 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 1024 = 1048 576 form factors. 2 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 unoccluded 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 approx63 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 i n the upper left corner. master MW 100 CDV 100 Figure 4.9: "Occluded". Images computed w i t h the basis and number of form factors as indicated. The master solution is shown in the upper left corner. imation. T h e error for m u l t i w a v e l e t s is also h i g h i n the s m o o t h areas. A t first g l a n c e this seems s u r p r i s i n g since m u l t i w a v e l e t s s h o u l d be a b l e to a p p r o x i m a t e the s m o o t h function- a d equately. H o w e v e r , it is a k n o w n p h e n o m e n o n i n G a l e r k i n d i s c r e t i z a t i o n s w h e n t h e s o l u t i o n c a n n o t be a p p r o x i m a t e d well. In o u r e x p e r i m e n t C D V wavelets are able to c o p e w i t h t h a t p r o b l e m better. Perceptually multiwavelets a n d C D V wavelets p e r f o r m v e r y s i m i l a r l y . W e leave the j u d g e m e n t to the reader. T h e C D V wavelets i n t r o d u c e a n u n e x p e c t e d g r i d p a t t e r n for s m a l l n u m b e r s of f o r m f a c t o r s (see F i g u r e 4.8 d i s c r e t i z a t i o n itself. CDV2000). It c a n n o t be seen i n the s i m p l e p r o j e c t i o n s h o w n i n F i g u r e 3.6. also t r i e d a h i g h e r n u m b e r o f f o r m f a c t o r s a m p l e s change. T h i s g r i d p a t t e r n is d u e to t h e G a l e r k i n n , samp We however, this d i d n o t result i n a n y B o t h h i g h e r o r d e r bases p r o d u c e c l e a r l y b e t t e r results t h a n t h e H a a r basis w i t h o u t post-smoothing. C o n s i d e r i n g t h e w i d e s u p p o r t o f the C D V wavelets t h e i r p e r f o r m a n c e i n o u r e x p e r i m e n t s for the o c c l u d e d scene is s u r p r i s i n g l y g o o d . F o r the case o f m u l t i w a v e l e t s i n the o c c l u d e d case we note the difference i n p e r c e p t u a l e r r o r a n d m e a s u r e d error. Perceptually multiwavelets wavelets, n o t so however i n t e r m s o f L\ error. G e n e r a l l y we f o u n d it h a r d i n o u r e x p e r i m e n t s e for t h e scene w i t h o c c l u s i o n . a p p e a r to p e r f o r m b e t t e r t h a n H a a r to choose the a p p r o p r i a t e p a r a m e t e r W e noticed a pretty a b r u p t change between a n inaccurate s o l u t i o n a n d a s o l u t i o n w h i c h r e q u i r e s a large n u m b e r o f f o r m factors, a n d c o n s e q u e n t l y o v e r l y long computations. CDV T h e s e l e c t i o n o f the p a r a m e t e r was m o r e difficult for m u l t i w a v e l e t s wavelets t h a n for H a a r . W e also t r i e d the m u l t i l e v e l like refinement method s u b s e c t i o n 2.3.3. H e r e we also f o u n d it h a r d to devise a g o o d s c h e m e for r e d u c i n g e. 66 and from 4.2 Comparison of Solution Methods We apply different solution methods to solve the system of equations resulting from a compressed kernel and evaluate their convergence. Our main interest is in the behavior of certain Krylov subspace methods ( C G N R , 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 p avg — 0.52 and with low average reflectivity p surface high- 1 low = 0.35: avg Sio Sii 0.4(a)/0.8(6) 0.4 0 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 p; = So Si s 0.8 0.8 0.8 2 s 3 0.8 s 4 0.4 s 5 0.4 s 6 0.4 s 7 0.4 s 8 0.4 s 9 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 y / / A 1 5 "0 V 2 10 F i g u r e 4.10: B o x scene w i t h o u t o c c l u s i o n . A n u m b e r i n d i c a t e s t h e i n d e x i of a surface Si. A l l surfaces have c o n s t a n t r e f l e c t i v i t y pi as specified i n t h e t e x t . T h e p o s i t i o n o f the l i g h t - s o u r c e is as specified i n F i g u r e 4.2. F i g u r e 4.11: B o x scene w i t h o c c l u s i o n . A l l s p e c i f i c a t i o n s are as for F i g u r e 4.10. c u b e is l o c a t e d at p o s i t i o n ( 1 , 1 , 1 ) i n s i d e t h e b i g c u b e . 68 T h e small 4.2.2 Experiments In o r d e r t o get a clear d i s t i n c t i o n b e t w e e n e r r o r c a u s e d b y t h e d i s c r e t i z a t i o n a n d error c a u s e d b y the s o l u t i o n m e t h o d we need to c o m p a r e the s o l u t i o n xj ^ o b t a i n e d i n the k-th k t h e i t e r a t i o n t o the e x a c t s o l u t i o n x o f t h e s y s t e m . step o f S i n c e s u c h a s o l u t i o n is n o t a v a i l a b l e we c o m p a r e w i t h a s o l u t i o n x o b t a i n e d after 80 steps o f the P i c a r d i t e r a t i o n . T h e c h a n g e i n the s o l u t i o n after 80 i t e r a t i o n s was i n the o r d e r o f 1 0 the l a t t e r was i n t h e o r d e r o f 1 0 - 3 - 1 1 , w h i c h is way b e l o w the a n a l y z e d error: or l a r g e r . W e c o m p u t e t h e error w i t h o u t w e i g h t i n g b y a r e a since here we are o n l y i n the c o n v e r g e n c e o f the solver, i n d e p e n d e n t of the r a d i o s i t y e q u a t i o n itself. interested W e use the following formula: error Efc ( ^ - x y/n k A l l c o m p u t a t i o n s were c o n d u c t e d w i t h a H a a r basis. A s i n i t i a l guess for the i t e r a - t i o n we chose t h e discrete r e p r e s e n t a t i o n b o f the l i g h t - s o u r c e s . T h e results p r e s e n t e d were c o m p u t e d for e = for e = 10~ IO - 4 , however we also m a d e m e a s u r e m e n t s a n d o b t a i n e d s i m i l a r results. computed with n samp = 3. T h e m a x i m u m n u m b e r o f levels is m = 5. 0.1 and e = 3 F o r m factors were T h e t o t a l n u m b e r of l i n k s c o m p u t e d i n these c o m p u t a t i o n s is g i v e n i n the f o l l o w i n g table: Scene # of potential links reflectivity # of c o m p u t e d links Unoccluded Occluded 2 516 580 10 066 320 High Low High Low 484 452 283 272 898016 606 398 W e present g r a p h s o f the c o m p u t e d e r r o r i n F i g u r e 4.12, F i g u r e 4.13, F i g u r e 4.14 a n d F i g u r e 4.15. 69 Unoccluded High Picard o GMRES CGNR rQ: t 0.0001 1 e - 0 6 b- # iterations F i g u r e 4.12: C o n v e r g e n c e o f i t e r a t i v e m e t h o d s for a test scene w i t h o u t o c c l u s i o n a n d h i g h average reflectivity. Occluded t High 0 . 0 0 1 r- # iterations F i g u r e 4.13: C o n v e r g e n c e o f i t e r a t i v e m e t h o d s for a test scene w i t h o c c l u s i o n a n d h i g h average reflectivity. 70 Unoccluded Low 0.1 i i icard ^IRES ;GNR'' i 3 ar 0.01 3~rr~rrr-^_^^. 0.001 0.0001 1 e-OS 1 e-06 1 e-07 1e-08 1 e-09 1 e-1 O ' I 1 # iterations F i g u r e 4.14: C o n v e r g e n c e o f i t e r a t i v e m e t h o d s for a test scene w i t h o u t o c c l u s i o n a n d l o w average reflectivity. Occluded Low # iterations F i g u r e 4.15: C o n v e r g e n c e o f i t e r a t i v e m e t h o d s for a test scene w i t h o c c l u s i o n a n d low average reflectivity. 71 4.2.3 The Interpretation g r a p h s c l e a r l y show t h a t C G N R t i c u l a r , for the o c c l u d e d scene C G N R has n o a d v a n t a g e over the P i c a r d i t e r a t i o n . has a s i g n i f i c a n t l y larger error. In p a r - O n the o t h e r h a n d , G M R E S has a clear a d v a n t a g e over P i c a r d i t e r a t i o n . A f t e r three i t e r a t i o n s the error is o n l y ^ ^ of the error of P i c a r d i n the u n o c c l u d e d cases, | a n d w | for the o c c l u d e d scenes. In a l l scenes three i t e r a t i o n s o f G M R E S are sufficient for a p e r c e p t u a l c o n v e r g e n c e . The m a j o r cost i n each step o f P i c a r d i t e r a t i o n a n d G M R E S vector multiplication. In wavelet is one sparse m a t r i x - r a d i o s i t y t h i s step i n c l u d e s wavelet decomposition and wavelet r e c o n s t r u c t i o n . G M R E S a d d i t i o n a l l y has a n u m b e r of v e c t o r - v e c t o r m u l t i p l i c a t i o n s the cost o f w h i c h however is n e g l i g i b l e c o m p a r e d to the o p e r a t i o n s m e n t i o n e d before. So i n t e r m s o f cost G M R E S has a clear a d v a n t a g e over P i c a r d . C G N R o n the o t h e r h a n d p e r f o r m s c l e a r l y worse t h a n P i c a r d since it is a b o u t twice m o r e e x p e n s i v e p e r i t e r a t i o n t h a n P i c a r d . 72 Chapter 5 Conclusions W e have p r e s e n t e d a f r a m e w o r k of r a d i o s i t y m e t h o d s w h i c h we h o p e w i l l h e l p researchers i n the areas o f n u m e r i c a l a n a l y s i s a n d c o m p u t e r g r a p h i c s to g a i n b e t t e r access t o the p r o b l e m s a r i s i n g i n the s o l u t i o n o f the r a d i o s i t y e q u a t i o n . W e c o m p a r e d s o l u t i o n m e t h o d s for wavelet r a d i o s i t y as well as the b e h a v i o r o f different b a s i s - f u n c t i o n s . W e believe t h a t t h i s analysis w i l l c o n t r i b u t e to a b e t t e r u n d e r s t a n d i n g of wavelet r a d i o s i t y . W e w i l l first give a b r i e f s u m m a r y of results o b t a i n e d i n o u r e x p e r i m e n t s . Subsequently we list a n u m b e r of q u e s t i o n s w h i c h arose i n o u r research, w h i c h we c o u l d n o t find a n s w e r e d i n the l i t e r a t u r e a n d w h i c h r e m a i n t o be p u r s u e d . 5.1 Results 5.1.1 Wavelet Bases In o u r c o m p a r i s o n b o t h e x a m i n e d h i g h e r o r d e r bases give p e r c e p t u a l l y b e t t e r results than t h e H a a r basis. than C D V wavelets t e n d to p r o d u c e p e r c e p t u a l l y m o r e a c c u r a t e results multiwavelets. 73 We observed that C D V wavelets allow computation of a more accurate solutions without 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 inaccurate 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 r e m a i n a n u m b e r o f u n a n s w e r e d questions a n d s t r a i g h t f o r w a r d e x t e n s i o n s to the w o r k presented i n the f o r e g o i n g c h a p t e r s . In the f o l l o w i n g we present a list of those w h i c h we c o n s i d e r to be m o s t i m p o r t a n t . • W e d i d not analyze G M R E S subsection ! • 2.3.3. i n the f r a m e w o r k o f the m u l t i l e v e l - l i k e m e t h o d s Employing G M R E S from i n t h i s f r a m e w o r k m i g h t l e a d to a significant s p e e d - u p since it c a n p o t e n t i a l l y result i n a c o m p u t a t i o n o f less f o r m factors, however, with possibly more iterations. A n a t u r a l q u e s t i o n t h e n w o u l d also be t h e a n a l y s i s o f p r e c o n d i t i o n e r s for G M R E S . A s i m p l e p r e c o n d i t i o n e r w o u l d be the m a t r i x R m e n t i o n e d i n s u b s e c t i o n 2.3.1. • M o d e l i n g the s o l u t i o n w i t h a m i x of different bases a p p e a r s to be a n i n t e r e s t i n g a p proach. A r e a s w i t h s m a l l g r a d i e n t i n the s o l u t i o n c a n be well a p p r o x i m a t e d w i t h a small n u m b e r of higher order basis-functions w h i l e a c c o r d i n g to o u r e x p e r i m e n t s g r a d i e n t s i n the s o l u t i o n are best a p p r o x i m a t e d w i t h the H a a r basis. Different bases a p p e a r to be easy to i n t e g r a t e as l o n g as each surface has o n l y one t y p e of basis. multiwavelets a n d C D V wavelets o f h i g h e r t h a n first-order, big c a n be i n t e g r a t e d Also, without m u c h effort. v • A n i m p o r t a n t e l e m e n t m i s s i n g i n the a n a l y s i s is h o w well t h e o r a c l e a c t u a l l y e s t i m a t e s the f o r m f a c t o r , e s p e c i a l l y i n the case o f t h e scene h a v i n g o c c l u s i o n s . 75 Bibliography [alpe93] B. K . Alpert. " A class of bases in L for the sparse representation of integral operators". SIAM J. Math. Anal, Vol. 24, No. 1, pp. 246-262, January 1993. [appe85] A . A . Appel. " A n 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. SIAM, 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. 2 [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 Workstation 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 Distributed Processing Environment". M.Sc. thesis, Program of Computer Graphics, Cornell University, Ithaca, N Y , January 1989. [chia97] G . Chiavassa and P. Liandrat. "On the effective Construction of Compactly Supported 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 Proceedings), Vol. 19, No. 3, pp. 31-40, July 1985. [cohe86] M . F. Cohen, D. Greenberg, D. Immel, and P. Brock. " A n Efficient Radiosity Approach for Realistic Image Synthesis". IEEE Computer Graphics and Applications, Vol. 6, No. 3, pp. 26-35, March 1986. [cohe88] M . F. Cohen, S. Chen, J. Wallace, and D. Greenberg. " A Progressive Refinement 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. SIAM, 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 Architecture" . Proc. First Eurographics Workshop on Parallel Graphics and Visualisation, September 1997. 77 [gers94] . ; R. Gershbein, P. Schroder, and P. Hanrahan. "Textures and Radiosity: Controlling 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 SIGGRAPH '84 Proceedings), Vol. 18, No. 3, pp. 212-222, July 1984. j [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, N J , 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 integral 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 Illumination". Technical report, Computer Science Division; University of California, Berkeley, July 1991. 78 [heck92] P. S. Heckbert. "Discontinuity Meshing for Radiosity". Third Eurographics Workshop 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. " A n 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 '86 Proceedings), Vol. 20, No. 4, pp. 143-149, 1986. SIGGRAPH [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 Radiosity 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 Definite Radiosity Matrix". Photorealistic Rendering Techniques (Fifth Eurographics Workshop on Rendering), pp. 227-243, 1995. [niev97] Y . Nievergelt. "Making Any Radiosity Matrix Symmetric Positive Definite". Journal of the Illuminating Engineering Society, Vol. 26, No. 1, pp. 165-172, 1997. [nish85] T. Nishita and E . Nakamae. "Continuous Tone Representation of ThreeDimensional 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 Include 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. P W S 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. 187192, 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, Princeton, N J , 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, P A 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] ford, [stoe80] Relaxation Methods in Theoretical Physics. R. Southwell. 1946. Introduction to Numerical Analysis. J . Stoer and R . Bulirsch. New Y o r k Berlin Heidelberg, [stol96] C l a r e n d o n Press, O x - Springer Verlag, 1980. E . Stollnitz, T . Derose, a n d D . Salesin. Wavelets for Computer Graphics. Morgan K a u f m a n n P u b l i s h e r s , S a n F r a n c i s c o , C a l i f o r n i a , 1996. [tell94] S. T e l l e r , C . F o w l e r , T . F u n k h o u s e r , a n d P . H a n r a h a n . Computer Graphics (ACM SIGGRAPH ing Large Radiosity Computations". '., Proceedings), [trou93] pp. 443-450, '94 1994. R. Troutman and N . L . Max. Element Methods". "Partitioning and Order- "Radiosity A l g o r i t h m s Using Higher O r d e r Finite Computer Graphics (ACM SIGGRAPH '93 Proceedings), pp. .209-212,1993. [veac97] E . Veach. Robust Monte Carlo Methods for Light Transport Simulation. thesis, S t a n f o r d U n i v e r s i t y , D e c e m b e r [whit80]. T . Whitted. [will97a] 1997. " A n i m p r o v e d i l l u m i n a t a t i o n m o d e l for s h a d e d d i s p l a y " . cations of the ACM, V o l . 23, N o . 6, p p . 3 4 3 - 3 4 9 , A . W i l l m o t t a n d P. S. H e c k b e r t . Ph.D. Communi- 1980. " 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" . Technical Report C M U - C S - 9 7 , Carnegie M e l l o n University School C o m p u t e r Science, [will97b] of 1997. . A . W i l l m o t t a n d P . S. H e c k b e r t . " 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 Algo- Rendering Techniques 97 (Proceedings of the Eurographics Workshop on Rendering in St. Etienne, p p . 1 7 5 - 1 8 6 , J u n e 1997. rithms". [xu94] W . X u a n d D . S. F u s s e l l . " C o n s t r u c t i n g S o l v e r s for R a d i o s i t y E q u a t i o n S y s t e m s " . Photorealistic Rendering Techniques (Proceedings of the Fifth Eurographics Workshop on Rendering), p p . 2 0 7 - 2 1 7 , J u n e 1994. [zatz93] H . R. Zatz. Illumination". 220, A u g u s t "Galerkin Radiosity: A Higher-Order Solution Computer Graphics (ACM SIGGRAPH 1993. 81 M e t h o d for G l o b a l '93 Proceedings), pp. 213- Appendix A A.l Two-Dimensional Wavelet Transform We generalize the one-dimensional M R A presented i n section 3.1 by forming products of scaling functions and wavelets. T h i s results i n a new two-dimensional scaling function and three new two-dimensional wavelets: ^toM^iV) = ^koM&y) = 0 kl MxtykM ^k (x)tpkAy) = ^loM^'y) ^k (x)'</> (y) 0 W i t h these definitions we get new two-dimensional dilation equations for scaling functions and wavelets: 7lj <t>j,k M 0 — 52 52 l - l + ; =o 0 Tlj + l ^JMM = —1 = -1 1 /i=o Tlj+ hko,i hki,h<l>i+\,i ,h 0 l —1 52 ; =o 52 9ko,l hk h<t>j+\,l ,h /i=o 52 52 ; =o /i=o 0 ^jMM + 0 0 82 u 0 h ,io9k h<!>j+i,i M ko u 0 0 rij + i-1 n i — 1 J + ; i =o 0 Bekaert et. al. 1 = o [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 P m m—1 P TP m = P TP m 0 m—1 m—1 + J2 QJTQJ 0 w i t h T can be written as + ] T Q TP 3 3=0 + 3 PJTQJ 3=0 • ( - ) A j=0 — Pj and the identities T h i s equality can be easily obtained using the identity Q = P i 3 3+ m— 1 P TP m m = P TP 0 0 + J 2 P J ^ T P J + i - p j T P J 3=0 and P TP -P TP 3+1 3+l 3 3 = (P -P )T(P -P ) = QJTQJ 3+l 3 3+l \ + P TP J 3 PJTQJ 83 rQjTPj. 3+1 + P TP -2P TP 3+l 3 3 3 1
- 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 |
File Format | 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 |
Graduation Date | 1998-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051668/manifest