BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR94-1443
ENTRY:: 1994-08-26
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Efficient Hierarchical Radiosity in Complex Environments
AUTHOR:: Smits, Brian Edward 
DATE:: August 1994
PAGES:: 120
ABSTRACT::
This thesis presents methods for speeding up the global illumination 
computations by using bounds on error to eliminate work that is not needed 
for a solution of a given accuracy.  This work makes the hieerarchical 
radiosity approach feasible for complex environments.

First, a new radiosity algorithm for efficiently computing global solutions 
with respect to a constrained set of views is presented.  Radiosities of 
directly visible surfaces are computed to high acccuracy, while those of 
surfaces having only an indirect effect are computed to an accuracy 
commensurate with their contribution.  The algorithm uses an adaptive 
subdivision scheme that is guided by the interplay between two closely 
related transport processes:  one propagating power from the light sources, 
and the other propagating importance from the visible surfaces.  By 
simultaneously refining approximate  solutions to the dual transport 
equations, computation is significantly reduced in areas that contribute 
little to the region of interest. This approach is very effective for 
complex environments in which only a small fraction is visible at any time.  
Our statistics show dramatic speedups over the fastest previous radiosity 
algorithms for diffuse environments with details at a wide range of scales.

A new approach for accelerating hierarchical  radiosity by clustering objects 
is also presented.  Previous approaches constructed effective hierarchies by 
subdividing surfaces, but could not exploit a hierarchical grouping on 
existing surfaces.  This limitation resulted in an excessive number of 
initial links in complex environments.  Initial linking is potentially the 
most expensive portion of hierarchical radiosity algorithms, and constrains 
the complexity of the environments that can be simulated.  The clustering 
algorithm presented here operates by estimating energy transfers between 
collections of objects which maintaining reliable error bounds on each 
transfer.  Two methods of bounding the transfers are employed with different 
tradeoffs between accuracy and time.  In contrast with the $O(s^2)$ time and 
space complexity of the initial linking in previous hierarchical radiosity 
algorithms, the new methods have complexities of $O(s$ log $s)$ and $O(s)$ for 
both time and space.  Using these methods we have obtained speedups of two 
orders of magnitude  for environments of moderate complexity while 
maintaining comparable accuracy.

Finally, the thesis describes a method for reconstructing the radiance 
functions across the visible surfaces given a global solution to the energy 
balance equations.  This approach greatly reduces artifacts resulting from 
the choice of constant basis functions used for the global solution.
END:: CORNELLCS//TR94-1443
BODY::
Efficient Hierarchical Radiosity In
Complex Environments
Brian Edward Smits
Ph.D Thesis
TR 94-1443
August 1994
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
EFFICIENT HIERARCHICAL RADIOSITY IN COMPLEX ENVIRONMENTS
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Brian Edward Smits
August 1994
Qc Brian Edward Smits 1994
ALL RIGHTS RESERVED
EFFICIENT HIERARCHICAL RADIOSITY IN COMPLEX ENVIRONMENTS
Brian Edward Smits, Ph.D.
Cornell University 1994
This thesis presents methods for speeding up the global illumination computa-
tions by using bounds on error to eliminate work that is not needed for a solution
of a given accuracy. This work makes the hierarchical radiosity approach feasi-
ble for complex environments.
First, a new radiosity algorithin for efficiently computing global solutions
with respect to a constrained set of views is presented. Radiosities of directly vis-
ible surfaces are computed to high accuracy, while those of surfaces having only
an indirect effect are computed to an accuracy commensurate with their contri-
bution. The algorithm uses an adaptive subdivision scheme that is guided by
the interplay between two closely related transport processes: one propagating
power from the light sources, and the other propagating importance from the visi-
ble surfaces. By simultaneously refining approximate solutions to the dual trans-
port equations, computation is significantly reduced in areas that contribute lit-
tle to the region of interest. This approach is very effective for complex environ-
ments in which only a small fraction is visible at any time. Our statistics show
dramatic speedups over the fastest previous radiosity algorithms for diffuse en-
vironments with details at a wide range of scales.
A new approach for accelerating hierarchical radiosity by clustering objects
is also presented. Previous approaches constructed effective hierarchies by sub-
dividing surfaces, but could not exploit a hierarchical grouping on existing sur-
faces. This limitation resulted in an excessive number of initial links in complex
environments. Initial linking is potentially the most expensive portion of hierar-
chical radiosity algorithms, and constrains the complexity of the environments
that can be simulated. The clustering algorithm presented here operates by esti-
mating energy transfersbetween collections of objects while maintaining reliable
error bounds on each transfer. Two methods of bounding the transfers are em-
ployed with different tradeoffs between accuracy and time. In contrast with the
0(s2) time and space complexity of the initial linking in previous hierarchical ra-
diosity algorithms, the new methods have complexities of O(s log s) and 0(s) for
both time and space. Using these methods we have obtained speedups of two or-
ders of magnitude for environments of moderate complexity while maintaining
comparable accuracy.
Finally, the thesis describes a method for reconstructing the radiance func-
tions across the visibles surfaces given a global solution to the energy balance
equations. This approach greatly reduces artifacts resulting from the choice of
constant basis functions used for the global solution.
Biographical Sketch
The author was born, lived, and will probably die.
Along the way he spent some time at the University of Oregon and collected
bachelor's degrees in Math and Computer Science. He then ended up at Cornell
doing computer graphics. Beyond this, little is known about his life.
iii
To my parents and sister
iv
Acknowle dgements
I would like to thank my adviser Don Greenberg. He has been very helpful and
supportive throughout my stay in the lab. As the director of the program he has
created a wonderful research environment both in terms of the people he has
brought here and the equipment he has obtained for us.
I would like to thank my minor committee members Nick Trefethan and Os-
car Rothaus for their assistance and flexibility in scheduling my exams.
I would especially like to thank Gary Meyer for introducing me to computer
graphics and spending a lot of time teaching me color theory and how to do re-
search.
Jim Arvo has been the most influential person to my research at Cornell. He
taught me the value and power of mathematical abstraction as well as being an
incredible co-author and friend.
Much thanks to the people that keep the lab running, Fran, Ellen, Jonathan,
Jim, Hurf, and Ben. Thanks Hurf for keeping my machine running, and always
being willing to change things to make it easier to get my work done. Many
thanks to Ben Trumbore for always being willing to give assistance in all things.
Greg Spencer has been a great office mate. PR made looking at my images
a much easier task, and the discussions on perception and high dynamic range
v
images were a lot of fun. Thanks to Jed Lengyel for always asking me how to
render a billion polygons, as well as for providing me a bottle of wine on March
2, 2017. Peter Pruyn was a great office mate, his views on life kept things inter-
esting, as well as giving me a great overview of Ithaca and New York. Thanks
to Chris Schoeneman for writing WM and for making the Painting with Light
project an incredible amount of fun. Thanks to Filippo Tampieri for encourag-
ing me to finish my thesis and sparing me the agony of writing a resume. Dani
Lischinski has been a great person to talk to and bounce ideas off of. Dan Kartch
was a great person to overlap with. His seminars were exceptional, especially
Chocolate 701.
Thanks to Jim, Erin, Filippo, Molly, Greg, Julia, Jed, Ben, Kevin, Mike, and
Katrina for being supportive over the last year. It was much appreciated.
Last, but not least, I would like to thank my parents for their encouragement
and confidence.
vi
Table of Contents
2
3
Introduction
1.1 Photorealism and Computer Graphics
1.2			Contributions. .			.
1.3			Organization
Global Illumination and Hierarchical Methods
2.1 Radiative Transport
2.1.1 Discretization . . .
2.1.2 Solving the Discrete System.
2.2			Norms and Error. .
2.3			Hierarchical Methods
2.3.1			N-body Problem
2.3.2 Radiosity and the N-body Problem
2.4			Hierarchical Radiosity . . .			.			.
2.4.1			Surface Hierarchies
2.4.2			Bounding Error . .			.			. .
2.4.3			Hierarchical			Algorithm
2.4.4			Complexity Analysis			.
2.4.5			Deficiencies in HR
2.5			Global vs Local Passes .			.			.
Importance
3.1 Introduction and Motivation
3.2			Importance .			. .			. .			.
3.3			Mathematics			. . .			.
3.3.1 Duality of Radiosity and Importance
3.3.2 Importance-Driven Refinement
3.4 Algorithm. .
3.4.1			Overview . . . .			. . .			.
3.4.2 Solving for Radiosity and Importance
3.4.3 Refining the Interactions .
3A.4 Estimating the Form Factors
vii
2
4
5
6
10
13
16
16
17
17
18
19
21
24
30
32
33
35
35
37
40
40
43
46
46
48
50
51
3.5
3.6
3.4.5 Assigning Initial Importance .
Results
Other Uses for Importance			. .			.
3.6.1 Initial Importance
3.6.2 Bounding Error Conservatively with Importance
Clustering
4.1 Introduction/Motivation
4.1.1 Overview
4.2 Mathematics
4.2.1 &links .
4.2.2 P-Links
4.2.3 Strategies for linking
4.3 Other Linking Criteria
4.4 Implementation
4.5 Results .
4.5.1 Link Complexity.
4.5.2 Refinement Complexity
4.5.3 Importance
4.6 Better Hierarchies
4
5
6
Local Pass
5.1			Introduction
5.2			Mathematical Background.
5.3 Reconstruction Using Hierarchical Radiosity Solution
5.3.1			Unbiased Reconstruction . . . .			.
5.4			Reconstruction with Texture Mapping
5.5			Results			. .			.
Conclusions and Future Work
6.1 Summary
6.2 Limitations and Future Work
6.2.1 Arbitrary Reflectance Functions
6.2.2 Arbitrary Geometries
6.2.3 Input and Output . .
Bibliography
viii
52
52
57
57
58
62
62
65
65
66
71
72
73
74
80
81
81
84
86
87
87
88
90
91
94
95
99
99
101
101
102
103
104
List of Figures
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
Quantities involved in the integral equation .			. . . .			.
Representation of? and K at a single point
Transfer between differential areas is a fttnction of the distance 4 and
angles Oi and 02 as well as the visibility between the two points
A patch and corresponding elements			. .			. . . .			. .
A hierarchical patch structure. . .			. .
Interactions between nearby surfaces occur at higher level ofdetail than
distant interactions. . .			. .
Pseudocodefor refining a link.			. . .			. . .			.
Pseudocodefor linking the children. . . .
Pseudocodefor solving the system described by the links. . . .
Pseudocodefor redistributing energy over all representations.
Pseudocodefor hierarchical algorithm			. .			. . .
Pseudocodefor creating the initial links. .
Two levels in a hierarchy for a iD patch. All valid interactions are
shown for shaded patches. . .			. .
Tree showing resulting interactionsfor patches at different levels in the
hierarchy.			. . .
7
9
10
19
20
24
25
26
27
28
29
29
30
31
Left, a radiosity solutionfor a labyrinth. Right, an importance solution
for the same model. . 			. . .			. 39
The radiosity and importance solutions together. . 			40
The duality ofradiosity and importance. . .						42
Transporting radiance and importance in a hierarchy. 			. 			45
Pseudocodefor an importance driven hierarchical algorithm			47
Solvingfor both importance and radiosity. .			. .			48
Gathering and shooting ofradiosity and importance.			49
Sweeping importance in a hierarchy. .			. .			49
Pseudocodefor refining an interaction using importance			.			51
Left, a view-dependent radiosity solution. Right, the same solution after
reconstruction53
The radiosity solution (left) and importance solution (right) of Fig-
ure 3.10, seen fromft&rther back. . .			. .			. .			54
ix
3.11
3.12
3.13
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
5.1
5.2
5.3
5.4
5.5
5.6
5.7
The radiosity solution (left) and importance solution (right)from even
flirther back
A B-only solution (left) and a comparable BI solution (right).
54
56
Two collections ofobjects can interact at different levels: a) conventional
HR links, b) an o-link requiring linear time and space, and c) a P-link
requiring constant time and space. . .			. . . . 64
Straightforward approach to bounding error. Each link represents a
maximum valuefor k			. 			. 			67
More efficient a-link between two clusters. 						68
fl-link between clusters			. 						72
Pseudocodefor creating an a-link. . . . . . . 			. 			. 						75
Pseudocodefor computing error bounds on an a-link			76
Pseudocodefor gathering energy across an a-link. 			. 						. 						77
Pseudocodefor sweeping energy throughout a hierarchy						78
Pseudocodefor creating a fl-link. . 			. 												78
Pseudocodefor bounding error on a fl-link						79
Pseudocodefor gathering energy across a fl-link. . 			. 			. 									79
Pseudocodefor refining two clusters			80
Link costfor environments ofdifferent sizes compared with theflinction
slog(s)82
Timefor solutions of increasing accuracy, both with and without clus-
tering. .			83
Solutions at different accuracies. With clustering (a-e), without clus-
tering (frg)
Image resultingfrom local pass (above) and solution used for the local
pass (below).			. . .			. .
Four different approaches to reconstruction.
The local pass uses the radiances L computed during the global solution
to determine the radiance at a pixel
patches contributing to point x in a hierarchy. .
Local pass and corresponding global pass. . . 			. .
Local passes with varying number ofsamples			. .			.
Unbiased (left) and biased (right) reconstruction
Local passes with texture mapping and transparent			surfaces.			.
x
84
85
88
89
90
95
96
97
98
Chapter 1
Introduction
1.1 Photorealism and Computer Graphics
The field of computer graphics generally involves the display of some sort of im-
age on a CRT. One area of computer graphics is devoted to producing photo-
realistic images. Ideally, these images would be indistinguishable from a pho-
tograph taken of the same scene. Photorealism is useful for many applications.
Designers, especially large-scale designers such as architects, often want to see
what the design will look like once it is built. By creating photorealistic images
of the object being designed, problems can be solved before it is built. Photoreal-
istic rendering is also very useful for simulation. Flight or driving simulators are
more believable if the image being looked at is impossible to distinguish from re-
ality. The techniques are also important for the entertainment industry, both for
virtual reality and for computer-generated effects for the flim industry.
In addition to achieving realism, computer graphics should also be concerned
with accuracy. It is often useful to know the illumination levels on various sur-
faces in an environment. This is especially true for illumination engineering,
2
which is concerned with determining how to adequately light buildings. Guar-
anteed accuracy is one of the main goals in photo-realistic rendering.
In order to produce photo-realistic images1 it is necessary to determine the
intensities of light on various surfaces. Solutions can be achieved by simulat-
ing how light is scattered throughout an environment. This type of simulation
is known as global illumination. Global illumination is difficult because light scat-
tered from a point on any surface in an environment can potentially reach a point
on any other surface. To exactly simulate all of these transfers of light through-
out an environment is computationally intractable. Because of this inherent dif-
ficulty in the problem, it is necessary to estimate the actual transfers by grouping
many transfers together. These approximations introduce error into the simula-
tion, which needs to be accounted for.
Much of the research in global illumination has been to reduce the computa-
tional complexity of the simulation. The philosophy of this thesis is that we don't
need exact solutions, we only need solutions that are accurate to within some
given error tolerance. We want to take as many short cuts in the simulation as is
possible, as long as we are within the error tolerance for the simulation. In order
to do this, we need to know the effects of the short cuts on the simulation. This
will be done by bounding the errors introduced by various approximations and
only using approximations that fall within the bound for the simulation.
1.2 Contributions
The main contribution of this thesis is efficiency. Global illumination is a diffi-
cult problem that becomes even more difficult as the size and complexity of the
environment increases. Model environments often consist of several hundred
3
thousand surfaces. Currently, the best physically accurate approach to comput-
ing global illumination solutions is hierarchical radiosity This approach is not
feasible for environments of more than several thousand surfaces because the
complexity of the problem grows quadratically with the number of initial sur-
faces. This limitation on the size of the environments that are able to be solved
makes radiosity impractical as a design tool for real world problems.
This thesis addresses two ways of increasing the efficiency of global illumi-
nation simulations. The first is a method of determining the amount of accuracy
needed in various parts of the environment based on the user's focus of inter-
est and the desired accuracy. By finding areas where less accuracy is needed, the
computation can be reduced significantly. Areas where less accuracy is needed
is determined by a closely related transport process, the transport of importance.
This approach is very effective in complex environments where only a small sub-
set of the environment is of interest to the user.
The second approach introduces a method for bounding the errors on very
coarse transfers of energy, thereby enabling coarse transfers of energy to be used
in areas where these errors will have little effect. Energy will be transferred be-
tween collections of objects, thereby eliminating the main deficiency of hierar-
chical radiosity; it's inability to allow coarser transfers other than those between
the initial surfaces. This approach changes the complexity of the initial phase of
the hierarchical radiosity algorithm from quadratic to O(s log s) in the number of
initial surfaces. Together, these approaches can speed up a global illumination
simulation by several orders of magnitude, and make hierarchical radiosity fea-
sible for environments more than an order of magnitude larger than previously
possible.
4
1.3 Organization
The next chapter of the thesis will describe the global illumination problem more
fully. It will give a mathematical formulation of the problem and the restrictions
that will be used for this thesis. In addition, it will describe the more relevant
previous work in global illumination.
The third and fourth chapters will describe new algorithms for speeding up
global illumination simulations by reducing the amount of work needed to corn-
pute a solution to a given accuracy. The third chapter describes a way of speed-
ing up radiosity solutions when only part of the solution is desired. This is useful
in large environments where only a few images are required. The fourth chapter
describes a way of reducing the costly initialization stage in hierarchical radios-
ity by clustering objects together.
The fifth chapter of the thesis describes a method for creating images based
on global illumination. The solution of a global illumination simulation is not
suitable for display. Further processing must be done to create an image without
noticeable artifacts.
Chapter 2
Global Illumination and
Hierarchical Methods
The simulation of global illumination involves modeling the interactions of light
in an environment. In the most general case, the global illumination model must
account for the emission, reflection, transmission, and absorption by the objects
in the model, and the emission and scattering by the media through which the
light is traveling. This computation must be done at all wavelengths within the
visible spectrum. In this thesis we have made several simplifying assumptions to
reduce the computational complexity. We assume that the media does not inter-
act with the light in any way. Also, we assume that there is no energy transferbe-
tween wavelength bands (no fluorescence) which allows each wavelength band
to be solved for independently. We can describe these processes by equations of
transfer, following the development presented in [ATS94j. Once a mathematical
model of the process has been defined, it can then be used to derive algoritluns
for solving approximations to the process.
5
6
Some of these algorithms will be described later in this chapter. In particular,
emphasis will be placed on hierarchical approaches as they are the most promis-
ing in terms of accuracy and efficiency. Most of the work in this thesis will draw
upon hierarchical approaches to global illumination.
2.1 Radiative Transport
The goal in global illumination is to solve for the radianceftinction of a given envi-
ronment. Radiance has units of watts per meter squared per steradian so it can
be thought of as the energy density at all points and all directions. Let ? be a
collection of surfaces and S2 be the unit sphere. The space ? x S2 consists of all
directions at every point on a surface of the environment. The radiance functions
X, for a given environment are real valued functions defined on ? x S2. Given an
emission function g E X which describes the emission of energy at every pomt
and in every direction on the surfaces of the the environment, we wish to deter-
mine the radiance function that satisfies the physical model we are simulating.
This radiance function f satisfies
f(x,w) = g(x,w) ? ?? p(x;w' ? w)f(x',w') cos(O')dw', (2.1)
where H?? is the hemisphere of incoming directions, 0' is the angle of the inci-
dent beam contained in the solid angle dw', and x' is a point on a distant surface
determined by x and w'. The function, p(x; w' ? w), is the bidirectional reflectance
distributionflinction or BRDF which relates the incident irradiance from direction
w' to the exitant radiance at point x in direction w. Irradiance has units of watts
per meter squared, and is the result of weighting the incident radiance by the
solid angle over which the incident radiance occurs. This equation is similar in
7
0)
Figure 2.1: Quantities involved in the integral equation
spirit to the rendering equation introduced by Kajiya [Kaj86] but uses radiomet-
ric quantities.
In this thesis, we will only be concerned with diffuse environments. That is,
environments where p(x; w' ? c with c ? ). Radiance functions are now
only a function of position, as the radiance is the same in all outgoing directions.
This gives a simplified integral equation which can be written as
f(x) g(x) + fHin p(x)f(x') cos(O')dw'			(2.2)
Much of the mathematics behind radiative transport becomes easier when
described using linear operators. The integral in equation (2.2) (or in equation
(2.1)) is often represented as an operator, giving an operator equation as follows
f(x) g(x) + Mf(x)
(2.3)
8
or
(I --H M)f(x,w) = g(x,w).
(2A)
The M operator encapsulates several different operations. In addition to the in-
tegration over the hemisphere, the x' in the integral requires determining the
closest point in ? along the ray determined by x and w'. This operation is equiva-
lent to projecting the radiance function over the environment onto a hemisphere
at x such that the radiance on the hemisphere in direction w' is the radiance at
point x'. As in [ATS94], we can define a linear operator ? which takes a radiance
function over an environment and returns a function expressing the incidentfield
radiance at each point, by projecting the environment onto the hemisphere at each
point.
We can now express the linear operator M in terms of ? and a new linear op-
eratorK. Again following [ATS94] we shall call K the local reflection operator. Let
h be a function describing the incident field radiance at each point in the envi-
ronment. The surface radiance function, ICh, can be expressed as follows
= f, p(x;w' ? w)h(x,w')cos(O')dw'			(2.5)
Figure 2.2 shows how ? and K operate on a single point.
Letting AT = 1 --H K? gives us the following integral equation for global illu-
mination
ATf =g.			(2.6)
There are several advantages to this formulation. The g operator encapsulates all
of the global interactions and the determination of visibility between two points,
and knows nothing about the local reflection model at a particular point. The K
9
G			K
Figure 2.2: Representation ofU and IC at a single point
operator encapsulates the local reflection model, but needs no information about
other points in the environment. This decomposition of M makes error analysis
easier [ATS94j, as well as making some of the developments in this thesis more
complete. Also, any analysis that is done with respect to these operators can be
more easily generalized to the similar operators for non-diffuse environments.
In general, it is easier to design algorithms which integrate over the surfaces,
?, rather than over all incident directions at x. This requires a change of variables
from a differential solid angle to a differential area. This can be expressed as
dw vis(x,x') cos(O)2dA
Ix --H
(2.7)
where x is the point at which the integration is taking place, x' is the point at
the center of dA, and the xCO%(811)2 is the change of variables from surface to solid
angle. The function vis(x, x') is 1 if no other surface in ? touches the line segment
between x and x' and 0 otherwise. This allows only the closest surface point in
a given direction to contribute to the radiance at point x (See Figure 2.3).
10
4
Figure 2.3: Transfer between differential areas is afltnction of the distance d and angles
Oi and 02 as well as the visibility between the two points
2.1.1 Discretization
The integral equation that needs to be solved is an infinite dimensional prob-
lem. The most common way to bring the problem into a finite dimensional one
is through boundary element methods. Boundary element methods restrict the set
of possible solutions to a finite-dimensional subspace Xn C X with n degrees of
freedom. For example, the environment may be discretized into n pieces over
which the radiance function is constant. However the discretization is done,
any f? c Xn can be expressed as a linear combination of the basis functions,
11
tui, u2, ..., uni spanning Xn. This linear combination can be written
--H--H Zajuj.
5='
(2.8)
We now need pick an f? E Xn that is close to the true solution, which corre-
sponds to picking the right values for each a5. This can be done by minimizing
the residual error, given by
ATfn --H g.			(2.9)
The residual error is a function and cannot be minimized directly. Instead, we
impose a finite number of conditions on the residual error. These conditions take
the form of linearfiinctionals. A linear functional is a scalar-valued linear func-
tion. In this case, the linear functional will map radiance functions to scalars.
We choose n linear functionals ? and set
--H--H0.			(2.10)
for i 1,2,...n. Equations (2.8) and (2.10) define a system of n linear equations
for the n unknown coefficients aj, which can be expressed as
or in matrix form as
?%fl,??5?g
--H--H0.			(2.11)
?jATu1 ... ?iATun
?ATu, ... ?ATun			a?
This equation can be expressed more compactly in matrix notation as
(2.12)
NL E.			(2.13)
12
A Concrete Example of Discretization
The first approach to radiosity in computer graphics was done by Goral [GTGB84].
This approach subdivided the environment into n planar elements. Over each el-
ement, Sj, the radiance was assumed to be constant, so
u5(x)			Aj ` ? Element5			(2.14)
u5(x)			O,x ? Element5			(2.15)
Each linear functional ?(f) computes the average value over element i. If Aj is
the area of Pi, the jth element, this can be expressed as
f f(?)ui(x)			A1j A f(x)dx			(2.16)
This linear functional and choice ofbasis functions is commonly known as Galerkin
methods. If the basis functions were expanded to be a polynomial basis over
each element, we have the Galerkin radiosity approach described by Heckbert
[Hec91j and Zatz [Zat93j. Each entry in the matrix then becomes (<ki, Afu5?
which can be expanded to
(?, Afu5? A1j LtU5(X)d? A1i At is3 p(x,y)cosOicosO2vis(x,y) dti.dx (2.17)
Ir--Hy2
Since the support of any two elements does not overlap, the first integral is zero
except when i j, in which case it is one. Since p(r, y) is ideal diffuse, and as-
sumed to be constant over Pi, it can be pulled out of the integral. The remaining
double integral encapsulates the geometric terms and is often called the form fac-
tor or configuration factor [SH81] and wriften
Fi5			A1i At ?j coscosO2vis(x,?? d?dx.			(2.18)
jx --H y12
13
Therefore, letting Sij be the Dirac delta function, we have the following form for
zith entry of the matrix
N?'j=Sij pjF??.			(2.19)
The resulting discrete system of equations has the form
1 --H piF11 ... --Hpi F1n			O[i
pnFni ... 1 --H pnFnn
Each entry of the matrix determines the radiance on R? due to S5. When the linear
functionals are applied to the source function g the result is the average radiance
emitted from each element. If we let Ej = ?(g) then Ej is the average emittance
of element?. Now we can write the system of linear equations as follows,
Pi Fin			?i			E1
= :			(2.21)
1 --H pnFnn			?n			En
1--HpiFii
?pnFn1
or in matrix notation as
NL = E.			(2.22)
Once the system of equations has been constructed, it can be solved in a va-
riety of ways.
2.1.2 Solving the Discrete System
The resulting system of equations can now be solved directly or by iterative tech-
niques. Direct approaches involve inverting the matrix and have a time complex-
ity of 0(n3). In general inverting the matrix is prohibitively expensive as well as
14
not being necessary? and iterative techniques can be used. These methods all re-
quire time 0(n2) per iteration. Usually a relatively few number of iterations is
needed before the solution is satisfactory.
In general1 either Gauss-Siedel or Jacobi iteration is used to solve the system
of equations. Jacobi iteration has a very intuitive physical analog. Energy starts
out at the emifters. Each iteration of Jacobi transfers energy one more bounce
through the environment. So Jacobi iteration corresponds to energy flowing out
from the lights and into the environment, being attenuated by the BRDF at each
surface it hits. Jacobi iteration computes the first k terms of the Neumann series
expansion of the inverse of the matrix N = I --H M! in other words, it computes
(I+MQM2tM3+...tMk)E.			(2.23)
This corresponds to the first k bounces of energy m the environment.
This approach to global illumination is too inefficient. We would like to be
able to solve environments with hundreds of thousands of patchs. A solution
for an environment this large would require a matrix with over 1010 entries, and
even iterative methods in problems this size are impossibly slow, even if the ma-
trix could be stored in memory.
In an aftempt to produce solutions quickly, Cohen et ai.[CCWG88j proposed
a solution technique they called progressive radiosity Their approach was based
on the observation that in general, many of the energy transfers in an environ-
ment are insignificant, and a comparatively few transfers, such as those from
light sources and some bright secondary surfaces, account for most of the illumi-
nation, in an environment. They proposed computing the energy transfers in or-
der of significance. Their approach was to transfer energy from the brightest sur-
face in the environment to the rest of the environment. This is considered one shot
15
in progressive radiosity. By repeatedly choosing the surface with the most un-
shot energy and shooting it to the rest of the environment, energy is propagated
throughout the environment. This method has been shown to converge and to
be equivalent to a matrix technique known as Southwell relaxation [GCS93].
In practice, progressive radiosity never reaches a converged solution. Each
shot takes linear time to compute, and every surface must shoot at least once rn
order to converge, so the method is still 0(n2) for a converged solution. In most
cases, liffle more than the direct illumination is computed before the progressive
solution is stopped.
There is a further problem inherent in the progressive radiosity approach.
Patches with the most unshot energy are shot first, so a patch which has a strong
local effect, but is too small to have a significant global effect is usually not shot.
This results in visual artifacts. In general, progressive radiosity solutions are too
dark because most secondary interactions are not computed.
Furthermore, computing the entries of the matrix requires a lot of work as the
double integral needs to be evaluated for each entry. Usually this is done using
various quadrature rules because until recently there was no analytic solution for
the form factor between two constant polygonal elements and there is still no an-
alytic solution for arbitrary surfaces or non-constant distributions. The analytic
solution of [GCS93] is expensive, and does not hold in the presence of occlud-
ing objects. The visibility term inside the integral is also an expensive function
to compute. Because of this, computing the matrix entries requires significantly
more effort than computing a solution to the system of equations.
16
2.2 Norms and Error
Norms provide a measure of the "size" of a function. We can define a norm over
the space of all radiance functions and get a normed linear space. This will allow
us to determine the `1size" of a radiance function. One common set of norms for
function spaces are the L?-norms defined by
If			[f J?? f(x,w) Pcos(O)d(w)dj'??			(2.24)
Three of these have particular interest. These are the cases where p is 1, 2, or
00. The 2-norm corresponds to the standard Euclidean distance for vectors. The
other two have physical interpretations. The i-norm of a radiance function is the
total power in the environment. The oo-norm is the maximum radiance over all
surfaces and in all directions.
Given two different radiance functions f, and h, we will often want to deter-
mine how far apart they are. We can compare the two functions by looking at the
norm of the difference I If --H h ?. For the Oc-norm this gives the maximum differ-
ence in radiance between f and h. For the i-norm we get the total power of the
difference function. For diffuse radiance functions, the i-norm can be written
IIf--Hg ? --H
ffs2 jf(x,w) --Hg(x,w) cos(O)d(w)dx
7rf f(x)--Hg(x)jdx.
2.3 Hierarchical Methods
(2.25)
(2.26)
Hierarchical radiosity provides an efficient way of computing, storing, and re-
fining the radiosity matrix. These features will be necessary for much of the re-
maining chapters of this thesis. This approach was developed from rapid hierar-
17
chical algorithms for the n-body problem and is based on the idea of computing
a solution to within a given error tolerance.
2.3.1 N-body Problem
The N-body problem is the computation of forces between n particles. Each par-
ticle exerts a force on all of the other particles, resulting in n(n --H 1) interactions.
The force of the interactions decreases as the inverse square of the distance, so
distant interactions have little effect. Rapid hierarchical algorithms for the N-
body problem are built on the idea that numerical computations are subject to
error and that solutions need only be computed to within a given accuracy. At
best, the solution accuracy is limited by the accuracy of the machine the com-
putation is running on. Because solutions need not (and cannot) be exact, and
because the effects of particles fall off rapidly with distance, the force due to a
cluster of particles can be computed as if there were only a single particle present.
The first of these hierarchical algorithins was developed by Appel [App85]. This
algorithm ran in O(n log n), but it was shown by Esselink [Ess89j that this com-
plexity was due to the preprocessing time required to build the hierarchy. The
first 0(n) algorithm was developed by Greengard and Rokhlin (Gre88] and in-
volves a multipole expansion for each cluster and the ability to quickly combine
the multipole expansions of nearby clusters to get the expansion of the combined
cluster.
2.3.2 Radiosity and the N-body Problem
Radiosity is similar to the N-body problem in several ways. The number of po-
tential interactions is the same, as any two surfaces can potentially interact, re-
18
suiting is the 0(n2) behavior of early radiosity algorithms. Like the force be-
tween two clusters, the energy transferred between two surfaces falls off as the
inverse square of the distance between them. This means that distant transfers
of energy between patches may have little effect on the receiver.
Although the two problems have similar features, there are several major dif-
ferences between them. Interactions between two particles can be computed
solely with the information about the two particles. In radiosity visibility be-
tween the two patches must be computed, and this requires determining if any
other object is positioned between the two patches. This makes computing the
strengths of the interactions much more difficult. Radiosity however, often re-
quires much less accuracy than the N-body problem, particularly when the re-
sults of the radiosity solution are intended for use in making an image. These
solutions need be no more accurate than the visual system is capable of detect-
ing. Thus, the ability to solve for the radiance in an environment to within some
error tolerance is even more useful for radiosity.
2.4 Hierarchical Radiosity
Hierarchical algorithms for radiosity require three properties.
o+ The environment has a hierarchy of possible representations.
o+ Errors in the interactions between two parts of the environment at some
level in the hierarchy can be quickly bounded.
o+ Interactions can be estimated.
We now show how to create a hierarchy of representations for a surface and
how to bound errors in the interactions between representations of two surfaces.
19
Figure 2.4: A patch and corresponding elements
2.4.1 Surface Hierarchies
Hierarchical representations for surfaces began with the patch-element structur-
ing of Cohen et al. [CGIB86j. The initial discretizafion of the environment creates
a set of patches with constant radiance. Each of these patches is further subdi-
vided into elements, again with constant radiance. (See Figure 2.4). This subdi-
vision can be automated during the solution process in areas where the radiance
function changes rapidly. This technique is known as adaptive meshing. This two
level hierarchy is used in a very constrarned way. Energy transfers occur from
patches to elements. Elements capture the details of the radiance functions and
there are many fewer source patches in the environment. The radiosities of the
elements are used to determine the radiosities of the patches. Each patch has a
radiosity equal to the area weighted average of its children's radiosities.
The patch-element structure is the beginning of a hierarchy, but is not flexible
enough for hierarchical radiosity? The two level hierarchy was extended by Han-
20
1'
Figure 2.5: A hierarchical patch structure
rahan et a!. [HSA91j to be a tree of arbitrary (or large) depth. Each initial surface
has a hierarchy of representations starting with a single patch at the top level,
with each lower level being a subdivision of the patches above it. (See Figure
2.5).
In both Hanrahan's description and the rest of this work, the tree structure
used is a quadtree. This is not a requirement of the hierarchical approach, but
simply for ease of implementation. Other implementations [LTG93] use other
tree structures.
One of the consequences of having multiple representations for a surface is
that we need a way of converting the radiance function over one representation
to the radiance function over another. It is necessary to determine the radiance
on a child I?? given the radiance on the parent R1. Also, given the radiance on all
21
the children, we must be able to determine the radiance on the parent. This can
be done using the linear functionals that are part of the projection. Let ? be the
linear functional associated with Ri. If oi is the coefficient of the basis function
for R1 then the value of ?? is
(2.27)
A similar approach works for transforming the values from the children to obtain
the value for the parent. The only difference is that we must sum over all of the
children of R1. This operation can be expressed as
i=1
(2.28)
The case we will use in this thesis will be for constant basis functions and the
linear functionals used for Galerkin radiosity (equation (2.16)). Therefore, the
radiance on the parent is pushed uniformly `6down" to all of the children, so each
child's radiance is the same as it's parent. In the other direction, the radiance
on the parent, drawn "up" from its children, is the area weighted average of the
radiance on the children.
2.4.2 Bounding Error
We now need a way of bounding the error in the interactions between two
patches. This can be done by looking at the expression for the transfer from a
surface to a point, and then determining how to place bounds on the transfer.
From equation (2.2) and equation (2.7) the value of the radiance function L on a
surface at a point x due to another surface Si can be expressed as
L(x)			k(x, y)L(y)dy
(2.29)
22
where the kernel k(x, y) expresses the radiance at x due to a differential area at
the point y. The kernel can be written
COS
k(x,y)--Hp(x) ???_y112vis(x,y)			(2.30)
where Oi is the angle formed by the normal of the differential area at x and the di-
rection given by y --H x and O2 is the angle formed by the normal of the differential
area at y and the the direction given by x --H y. We can use these equations to com-
pute the interaction between two surfaces. To compute the exact radiance func-
tion across the receiver, L(x) must be evaluated at every point on the receiver. In
general, some set of basis functions are used to represent the radiance functions
and quadrature rules are used to compute the coefficients for the basis functions.
This reduces the problem to evaluating L(x) at some fixed number of points on
the receiver.
We can now bound the oo-norm of the error by bounding the difference be-
tween the maximum and minimum transfers. The maximum transfer can be
computed by finding the maximum value of k(w, y) over all points on the two
patches and then integrating the product of this with the radiance of the source.
We will denote the maximum value of a function f over some range A x B as
[flAB --H WE%f(X?Y)			(2.31)
?EB
Using this notation and moving the maximum value of the function outside the
integral, the bound on the energy transfer may be written
ma?xL(x) ? [kims, ft L(y)dy.			(2.32)
This expression allows for arbitrary distributions of radiance across the source.
The minimum can be computed similarly, with the minimum set to zero if any
two points are occluded.
23
Computing the maximum and minimum values of k exactly is a difficult
problem. In general, all we really need is a good estimate of these values. This
can be done by taking a set of sample points on both the source and the receiver
and computing the maximum and minimum values for k over all pairs of points.
For the maximum values we should assume that all pairs are unoccluded. This
handles the case where there is only a small unoccluded path between the source
and the receiver. In general, since the unoccluded kernel values change smoothly
over the two surfaces, this method gives a good estimate of the error. Difficulties
can occur when two surfaces are very close to each other relative to their sizes.
In these cases, the pairs of points may not estimate the maximum value of the
kernel very well. In general, though, when two patches are very close together,
there is enough error that the surfaces will not interact at this level anyway, so
the problem disappears. Another potential problem is that when two patches
are very close together or touching, the kernel value becomes very large. This
can cause subdivision to occur indefinitely. This problem in hierarchical radios-
ity is dealt with by determining a minimum size for subdivided patches. This
prevents excessive subdivision in the corners and other places where surfaces
touch.
The i-norm of the error can be bounded very similarly. As the i-norm of the
error is the amount of incorrect energy on the receiver due to the source, we can
bound this by taking the bound on the radiance and weighting it by area of the
receiver. This gives a conservative bound on the amount of incorrect energy on
the receiver.
24
Figure 2.6: Interactions between nearby surfaces occur at higher level ofdetail than dis-
tant interactions
2.4.3 Hierarchical Algorithm
We can now describe the hierarchical radiosity algorithm. First the bounds on
an interacfion between two coarse representations is checked. If this bound is
less than some epsilon, then the interaction is accurate enough, and it is left as
is. If the bound is greater than the epsilon, then more detailed representations of
the surfaces are used and the bounds are checked on all of the new interactions.
In general, distant interactions are handled at coarse levels in the hierarchy and
close interactions are handled at more accurate levels (See Figure 2.6).
Now we describe how the level of interaction is determined for two patches.
Assume we have an interaction L?i between a receiving patch Ri and a source
25
Refine(I5?1?, I,?)
s ? FindError(i, j)
if S> ?
I?I\ tIj?1
SubdivideAndRefine(I???, I, c)
Figure 2.7: Pseudocodefor refining a link.
patch Sj. First we find the error (5 in the interaction. If this error is less than the
given error tolerance ? then the interaction is accurate enough and is kept. Other-
wise we need to allow the interactions to occur at a more accurate level in the hi-
erarchy. This is done by subdividing the larger of S5 and Ri and creating new in-
teractions between the children of the subdivided patch and the non-subdivided
patch. These new interactions are then recursively refined in the same way the
original interaction was refined. This is expressed as pseudocode in Figure 2.7.
The original algorithm presented by Hanrahan et al.(HSA9l] used a slightly
different approach. Rather than using a bound on the error as determined by
the FindError function as the refinement criterion, they use the magnitude of the
form factor between the two surfaces. As the error decreases faster than the mag-
nitude of the form-factor, they are able to control the error in the solufion.
The SubdivideAndRefine routine creates a link between the subdivided patch
and the the original patch if the two patches interact. Patches interact if some part
of the source is visible from some part of the receiver. CreateLink computes the
form factor between the source and the receiver so that energy can be transferred
across the link. Figure 2.8 contains the pseudocode for this procedure.
The collection of all interactions 1 represents the transport matrix. Since links
26
SubdivideAndRefine(I???, I, ?)
if A? > Aj and Aj > Arnin
Subdivide(R?)
foreach Child i' of Ri
if Interad(i',5)
I ? Iu Crea?Link(i', j)
Refine(I5??', I, ?)
else if A5>? Aj and Aj > Amin
Subdivide(S??)
foreach Child j' of Sj
if In?ract(i, j')
I ? Iu CreateLink(i, j')
Refine(Lt??, I, ?)
else I ? Iu CreateLink(i, j)
Figure 2.8: Pseudocodefor linking the children.
27
SolveSystem(I, L,E)
repeat
foreach I??? ini
ti?L'i+PiThh
L?L'
foreach toplevel patch Rj
Sweep(R?, 0)
until convergence
Figure 2.9: Pseudocodefor solving the system described by the links.
at higher levels in the hierarchy connect patches with children together, these
links account for more than a single entry of the original matrix. It is a block
structured matrix. We can solve this system by using Jacobi iteration. Energy
is transferred by each link from the source patch to the receiving patch in the
same way as Jacobi iteration with a more standard representation of the matrix.
As each surface has multiple representations, and each link transfers energy to
a single one of these representations, the energy gathered at a some level of the
hierarchy needs to be distributed to other levels in the hierarchy. This is done
by "sweeping" the energy gathered at each level in the hierarchy to all the other
levels. Solving the system of equations described by the interactions 1 and the
emission function E is described with pseudocode in Figure 2.9.
The sweep function operates on the coarsest representation of each patch, and
accumulates the radiance at each level of the hierarchy to get the correct radiosi-
ties at the finest level of subdivision which are called leaves. Once this is done,
the radiance at any higher level patch is simply the area weighted average of
28
Sweep(R?,Ldown)
if Rj is a leaf
Lup? Li t Ej t Ldown
else
Lup?O
foreach Child i' of l??
x ? Sweep(I??i, Lj + Ldown)
Lup ? Lup t X * Ai?/Ai
Lj ? Lup
return Lup
Figure 2.10: Pseudocodefor redistributing energy over all representations.
the radiances of its children. The recursive function Sweep is expressed as pseu-
docode in Figure 2.10
Hierarchical radiosity is an iterative algorithm. Starting with a coarse repre-
sentation of the transport matrix, it refines the interactions and then solves the
new system of equations. By decreasing the error tolerance ? for each iteration,
solutions of higher and higher accuracy are achieved. This is expressed in pseu-
docode in Figure 2.11.
The CreatelnitialLinks routine creates the initial coarse representation of the
transport matrix. This is done by checking all pairs ofpatches in the environment
to determine if they interact, and if they do, then creating a link between them.
This set of initial links is the coarsest representation of the transport operator in
HR. Figure 2.12 shows the pseudocode for this operation.
This approach towards radiosity is highly efficient because only a linear num-
ber of interactions are needed to represent the transport matrix, as opposed to the
29
HRad(?, E)
L?E
CwateIrntia1Links(? I)
for ? decreasing from oc to e0
SolveSystem(I, L, E)
foreach I5?i of I
RefineInteracfion(I5??, I, ?)
Figure 2.11: Pseudocodefor hierarchical algorithm.
CreateInitiaWinks(?, I)
foreach Ri in g
foreach S5in?
if Interact(R?, Sj)
I ? Iu CwateLink(R?, Sj)
Figure 2.12: Pseudocodefor creating the initial links.
30
Levelk-1
Level k
Figure 2.13: Two levels in a hierarchyfor a iD patch. All valid interactions are shown
for shaded patches.
quadratic number required for full matrix techniques. This size reduction allows
the matrix to be stored. More importantly, the coarse approximations can be re-
fined into higher accuracy solutions without losing the work done to create the
original coarse approximation. We now show that this approach will only create
a linear number of interactions.
2.4.4 Complexity Analysis
In order to show that only a linear number of interactions are created we start
with a simpler case described in Hanrahan et al.[HSA9lj. Consider the iD prob-
lem of n patches represented as adjacent line segments. These n patches can be
grouped into a hierarchy by merging neighbors together to form a binary tree.
Now allow the patches to interact with each other. The error in this case is a func-
tion of the size of the patches and the distance between them. For a fixed error
tolerance and size, there is a minimum distance for which the two patches can
interact without exceeding the error threshold. For our simple case, we assume
that this means that a patch cannot interact with it's neighbors, but can interact
with patches that are further away than this (See Figure 2.13). Therefore, the par-
31
Figure 2.14: Tree showing resulting interactions for patches at different levels in the
hierarchy.
ents of a node take care of all interactions except for it's siblings, and the children
of it's parent's neighbors (See Figure 2.14). Two of these five nodes are neighbors,
but the other three can have links created to them. Therefore, each internal node
will be linked to three other nodes. At the leaves, a patch must interact with it's
neighbors, so leaf nodes interact with five other nodes. As this is a constant num-
ber of interactions per node, and there are a linear number of nodes, there are a
constant number of interactions in the environment.
In 3D, the proof is the similar. Each patch in the hierarchy can interact with all
patches such that the error in the interaction is less than some constant. Although
the orientation of the source and the receiver have an effect on the error in the
interaction, it is most strongly determined by the distance. The error decreases
with the square of the distance. Because of this, there is a distance beyond which
32
all patches can interact with the receiver. There is also a greater distance beyond
which all patches will interact with the receiver's parent. The only patches the
receiver will need to interact with are those that are further away than the mini-
mum distance and closer than the minimum distance for the parent. For each en-
vironment there is a constant such that the number of surfaces within this range
for all receivers is less than that constant. Therefore, each patch interacts with
only a constant number of patches in the environment. In most radiosity envi-
ronments there is occlusion. This is not a problem as occlusion will reduce the
number of interactions in the environment. Therefore there are still a linear num-
ber of interactions. Hanrahan et al.[HSA9lj show that for perpendicular patches,
the number of interactions is proportional to the square root of n, so it may be
plausible that the orientation of the surfaces reduces the number of potential in-
teractions.
2.4.5 Deficiencies in HR
Hierarchical radiosity employs multiple representations for the surfaces, and
bounds on the transfer between these surfaces. It is limited in several ways. First,
the coarsest representations are the initial surface descriptions. This means that
the algorithm is at least 0(s2) in the number of input surfaces. This is due to the
inability of hierarchical radiosity to create a coarser representation of the envi-
ronment than that originally given by the top level surfaces. This limitation is a
serious problem when the number of initial surfaces is large, as it is in complex
environments. This thesis will present a method for clustering surfaces together
to obtain a coarser representation for the environment. This will require bound-
ing the error in the transfers between two clusters.
33
A second limitation of hierarchical radiosity is that the size of the environ-
ment may be too large to get the accuracy needed in some region. It is not possi-
ble with the hierarchical radiosity algorithm to vary the accuracy of the solution
throughout the environment. In particular, if a solution is needed for a subset of
an environment such as a particular view or a room in a large environment, there
is no way of determining the accuracy needed in other parts of the environment
to guarantee the accuracy over the subset the user is interested in. This thesis de-
scribes a way to control the amount of accuracy wanted in various parts of the
environment.
2.5 Global vs Local Passes
Recent trends in realistic image synthesis have been towards a separation of the
rendering process into two or more stages[Rei92,CRMT9l,LTG93j. One of these
stages solves for the global energy equilibrium throughout the environment. Hi-
erarchical radiosity does this fairly well, and can be made even more effective.
Once the global solution has been computed, we often want to make an im-
age of the result. Simply displaying the basis functions used for the solution is
generally not good enough, especially if the basis functions are piecewise con-
stant. The eye is very sensitive to discontinuities in the radiance functions. Early
approaches to radiosity used Gouraud [Gou7l] shading across the patches to
smooth the radiance functions to eliminate the discontinuities. Unfortunately,
the eye is also very sensitive to first derivative discontinuities as well. Better ap-
proaches have been proposed by Reichert [Rei92j and Lischinski et al.[LTG93j.
The approach presented by Reichert is an image space approach that recon-
structs the radiance function at each pixel. A coarse global solution is first com-
34
puted. Then, at each pixel in the image the radiance is determined by finding the
surface point used for that pixel and evaluating
f(x) g(?) + Ain p(x)L(x') cos(O')dw'			(2.33)
where L is radiance function computed by the coarse global solution. This in-
tegral can be computed by determining the contribution to point x from all sur-
faces in the environment. Clearly this is an expensive operation to be performed
at every pixel. Its advantage is that it recomputes the radiance function across
the image plane with no direct meshing artifacts.
The approach presented by Lishinski et al.is an object space approach. First
a global solution is computed using a mesh that incorporates the discontinuities
in the radiance functions. Once a global solution is computed, a new mesh is cre-
ated containing the discontinuities and enough additional patches so that the ra-
diance functions are well behaved over each patch. Then the integral from equa-
tion (2.33) is evaluated at the corners of all the mesh elements, and the resulting
solution is displayed using either Gouraud shading or quadratic interpolation.
This approach is view-independent, making it well suited for walkthroughs and
other interactive uses. However, it currently only works for diffuse, polyhedral
environments. Also, for efficiency reasons, only discontinuities from primary
light sources are computed, so shadowing effects from secondary sources may
be missed.
Chapter 3
Importance
3.1 Introduction and Motivation
View-independent global radiosity algorithms have two major drawbacks: they
oversolve globally and undersolve locally. The algorithms oversolve globally in
that they afteffipt to compute radiosities to a uniform precision throughout the
environment--Heven on surfaces hidden from all useful points of view. They un-
dersolve locally in that a single global radiosity solution is seldom adequate un-
der close inspection--Hlocal effects such as shadowing and color bleeding among
small objects are often lost when the radiosity of the entire environment has to be
computed to a uniform precision. Thus, the utility of a purely view-independent
solution diminishes as the environment complexity or the required accuracy in-
creases.
To address this problem, several adaptive meshing schemes have been de-
vised to increase the level of approximation where significant intensity gradients
occur [BMSW91,CF9O,CGIB86j. These algorithms can achieve good results with
fewer surface elements. For complex environments, however, the cost of creat-
35
36
ing a fine mesh to captuw every illumination detail, whether visible or not, may
still be far too high.
In practice, radiosity implementations often rely on some form of additional
intervention by the user in order to handle complex scenes. For example, in or-
der to keep the total number of surface elements small, the user may need to
supply meshing hints based on the anticipated set of views. This approach has
the advantage of saving work in areas that are unimportant to the final image.
However, it also tends to sacrifice global accuracy, as there is no obvious way for
the user to predict the level of meshing required for distant objects to have their
proper effects on the visible parts of the scene.
The converse of this problem arises during the solution process: which of the
vast number of interactions in a complex environment are significant enough
to evaluate? Consider, for example, computing a radiosity solution for a large
building. In principle, light leaving a small surface on one floor could reach any
other floor by some circuitous path of stairwells and corridors. A radiosity algo-
rithm that accurately computed all such interactions would be highly impracti-
cal.
The LIR algorithm focuses effort on the significant energy transfers, quickly
approximating the insignificant interactions. However, because the algorithm is
still view-independent, it does more work than necessary in complex environ-
ments when only a single view or set of views is required. Indeed, it is not dif-
ficult to construct models for which even this rapid algorithm is impractically
slow.
One way to make radiosity practical for complex environments is to incor-
porate a notion of view-dependence. For instance, we might first use some
37
form of visibility preprocessing to determine the surfaces that are directly visi-
ble [ARBJ9O,TS91]. To compute the radiosities of these surfaces, all surfaces con-
tributing energy must be taken into account. But these contributing surfaces, in
turn, receive energy from surfaces that are still further away. However, in gen-
eral the effect of distant interactions will be less important to the visible scene.
In short, to make effective use of view-dependence it is necessary to consider all
potential interactions among surfaces, but to compute each interaction only to
an appropriate level of accuracy.
In this chapter, we describe such an algorithm: an extension to hierarchical
radiosity that refines the interactions contributing the most error to the view-
dependent solution. The algorithm makes use of importance fttnctions, which
have been studied extensively in neutron transport theory [KW86,Lew65]. In our
context, importance functions describe how the radiosity originating at a given
patch influences the visible surfaces. The algorithm we propose combines es-
timates of importance and radiosity to drive the global solution, allowing it to
exploit view-dependent information as part of an adaptive refinement scheme.
3.2 Importance
Illumination algorithms can be divided into two categories: those that simulate
the propagation of all light throughout an environment, typified by radiosity al-
goritlims; and those that simulate only the light reaching the eye, typified by ray
tracmg. These two strategies are in some sense dual: the former simulate the
process of photons emanating from sources of light, while the latter trace rays
that emanate from the eye, but behave very much like photons in every other re-
spect. The two strategies have advantages for modeling different modes of light
38
scattering. Indeed, the many bidirectional ray tracing and multi-pass radiosity
methods suggested in the last few years [Arv86,CRMT91,Hec90,SP89,SAWG91,
WCG87j exploit the complementary nature of these two processes.
An analogous duality appears in neutron transport theory, where equations
similar to those of radiative transfer are used to predict neutron flux [Dav57]. If
only the flux impinging on a small receiver is required, it is typically the adjoint
of the original transport equation that is solved, in effect reversing the direction
of neutron migration back toward the source. This strategy is closely related to
the various ray tracing approaches for global illumination [Kaj86,WRC88j.
The efficiency of these "backward" methods depends on the ability to quickly
find paths leading back to the source. But determining these paths is equiva-
lent to solving the transport equation in the forward direction. Thus, solving the
transport equation in one direction would seem to require its prior solution in
the other direction. However, through variational methods, nuclear engineers
have been able to use this interdependence to advantage. These methods allow
approximate solutions to the two transport equations to be combined, yielding
an overall solution with higher accuracy than either component alone [Dav57,
Lew65].
The effectiveness of these techniques suggest that in the realm of global illu-
mination it may also be possible to exploit the dual nature of light transport more
effectively than has previously been recognized. Rather than using these dual
processes in multiple passes to simulate different light scattering modes, we can
instead use the two processes together--Hsolving them simultaneously--Hin order
to produce an accurate result more quickly.
As an illustration of our approach, consider the labyrinth model depicted in
39
Figure 3.1: Left, a radiosity solutionfor a labyrinth. Right, an importance solutionfor
the same model.
Figure 3.1. The model is illuminated by several light sources, shown in white.
The radiosity solution due to these lights is shown in red. In an analogous fash-
ion, the camera can also "illuminate" the scene with a new quantity, "impor-
tance." Figure 3.1 also depicts the importance solution due to the camera in
green.
Our algorithm uses radiosity and importance together to accelerate the ra-
diosity solution for the visible scene. It does this by refining estimates of the
transport equations most where the interaction of radiosity and importance is
highest. Thus, in Figure 3.2 the algorithm uses the finest mesh for the parts of
the scene that are yellow, somewhat less meshing for the parts that are orange,
and even less for the parts that are red, green, or black--Hi.e., in regions of little
importance, little brightness, or little of eithet
40
Figure 3.2: The radiosity and importance solutions together.
The remainder of this section shows how the duality of radiosity and impor-
tance can be established more formally, and develops some mathematical tools
for using these quantities together to drive a hierarchical algorithm.
3.3 Mathematics
3.3.1 Duality of Radiosity and Importance
Importance is related to the adjoint of the global illumination operator AT, which
will be written AT* The adjoint of AT can be easily expressed using the K and ?
operators from Chapter 2. Arvo et al. [ATS94j show that both K and ? are self-
41
adjoint so K = )c? and g = ?*. Therefore, the adjoint of AT is
= (I--H			= I--H			= I --H
(3.1)
Therefore, the adjoint integral equation for global illumination can be expressed
as
(3.2)
Once we have projected the continuous problem into a discrete domain, we
are working with matrices. The adjoint of a matrix is simply its transpose, which
makes working with it relatively easy. This gives two sets of related equations
NL = E			(3.3)
NTz = v			(3.4)
where V is a receiver, dual to the source and the quantity Z is the importance vec-
tor dual to L. To understand the significance of R and Z, suppose that we wish
to compute a scalar function of the radiosity solution v(L), rather than the entire
solution L. For example, the function v might compute the average radiosity vis-
ible through a certain aperture, such as a pixel or the entire image plane. Since
any linear function of L can be expressed as an inner product, we can write v(L)
as the product vTL, where V? gives the contribution of Lj to the scalar function v.
We can then use equations (3.3) and (3.4) together to derive an alternate formu-
lation of v(L) with respect to Z:
v(L) = vTL = (NTZ)TL = zTNL = zTs.
(3.5)
Thus, each element Z5 of the importance vector Z gives the contribution made
by a unit of emittance at patch j to the scalar function v(L). In other words, if v is
42
Figure 3.3: The duality ofradiosity and importance.
chosen for a particular view, then Zj gives thefraction ofradiosity emitted at patch j
that ultimately reaches the eye.
Note that since the vector V plays a role dual to S in equations (3.3) and (3.4),
we can think of V as describing the inifial "emittance of importance." In this
sense, we can think of every patch j as having associated with it both a steady-
state radiosity Lj coming from all patches in the environment but originating at
the lights 5, and also a steady-state importance Zj coming from all patches but
originating at V, determined by the eye (Figure 3.3).
While radiosity and importance are similar in many respects, the two quan-
tities are not exactly the same. Radiosity is a flux density? measured in watts per
meter-squared, whereas importance is defined as a fraction, and is therefore a di-
mensionless quantity. This distinction has ramifications when we distribute the
two quantities in a hierarchical system. Because of the difference in quantities,
43
the importance of a parent is the sum1 not the average, of its children's impor-
tances, and the importance of the parent is distributed to it children according to
the relative area of each child compared to the parent [SAS92].
3.3.2 Importance-Driven Refinement
For most transport equations, we must use a discrete approximation to the ex-
act transport operator. The global illumination problem is no exception; the ma-
trix approximation of N generally contains imprecise form-factor estimates and
other assumptions that introduce error. Howevei, the approximate matrix can
be expressed as a perturbation of the exact operator, and refined using estimates
of radiosity and importance, as follows.
Let N be an approximation to N, with N = N + AN, and let L be the solution
to the corresponding transport equation, which is an approximation to the ex-
act radiosity vector L. Assuming exact emittances S, the approximate radiosity
transport equation can be written as
which is equivalent to
NL			(3.6)
NL = S --H ANL.			(3.7)
Thus, we can rewrite the approximate transport equation with exact emittances 5
in terms of an exact transport operator with perturbed emiftances 5 --H ANL.
Since the importance of a patch describes the fraction of emitted radiosity that
reaches the eye, we can now derive an expression for how the error in the emit-
tances affects the view-dependent function v. Writing the visible error as v(L --H
44
and using equations (3.5) and (3.7), we have:
v(L--HL) --H RTL?RTL
zTs - zTNL
--H--H zTs?zT(s?ANL)
zTANL.
(3.8)
Thus, the quantity zTANL is the error that the approximations to N and L con-
tribute to the view-dependent function.
Unfortunately, in practice it is not possible to compute the importance vec-
tor Z exactly, since computing the global importance solution is as difficult as
computing the global radiosity solution. However, we can compute an approxi-
mation Z to the exact importance vector Z, and refine it as we refine N and L. In
this case, we can use the quantity 2TANL as a close approximation to the actual
error zTANL.
To incorporate this insight into a hierarchical radiosity algorithm we begin by
considering the interaction between two patches i and j. In general, computation
of the energy transfer between any two patches will be approximate, due to the
approximations inherent in the discrete operator N. Let ?? denote the product
ZiANjiLd . The quantity % approximates the error contributed by the interac-
tion of patch i and j to the view-dependent function v. By reducing Si5 over all
pairs (i, j), we can make the magnitude of the overall error zTANL arbitrarily
small.
We reduce the error in an interaction by subdividing one of the two patches
involved and computing new interactions for the refined system. The net effect
of this refinement is that transfers of radiosity from bright patches and transfers
45
Important
Important
Bright
Bright
Figure 3.4: Transporting radiance and importance in a hierarchy
of importance from important patches will generally be treated with greater ac-
curacy at a lower level in the hierarchy.
Figure 3.4 illustrates this idea for two patches i and j with high importance
and high radiosity? respectively. The diagram on the top shows a single link from
patch? to patch j. The link carries radiosity from? to j, in the direction of the ar-
row, and importance from j to i, against the arrow. Since the effect of these inter-
actions on the error is not great, they take place at a high level in the hierarchy.
46
On the other hand, the diagram on the bottom shows the transfer of radiosity
from j to i and of importance from? to j. These transfers, which have greater
potential effect on the error, take place at a more refined level of the hierarchy.
In this way the algorithm can put the most work into refining the parts of the
transport process where the impact of error is greatest.
Because of the duality of radiosity and importance, the same error criterion
ZjANijLj used to refine radiosity can be used to refine the importance solution
as well. Moreover, since estimates of importance are used to drive the radiosity
solution and viceversa,we need to refine both solutions to the same level of accu-
racy. Therefore, the hierarchical system used for radiosity is also appropriately
refined for importance, so a single hierarchical system suffices for both.
3.4 Algorithm
The view-dependent radiosity algorithm we describe here is very srmilar to the
brightness-weighted hierarchical algorithm proposed by Hanrahan etal., with
the crucial difference that in the view-dependent algorithm, importance as well
as radiosity plays a role in refining interactions.
The algorithm iteratively computes an accurate radiosity solution for visible
patches by refining the interactions between any two patches? and jwhose es-
timated error ZiANijLj exceeds a given tolerance c.
3.4.1 Overview
The algorithm takes as input the initial emiftance of radiosity S and impor-
tance R for each patch, along with an initial set of interactions I among the
patches. The set I contains a link between every pair of patches that are not com-
47
ImportanceDrivenRadiosity(S,1?, I)
for ? decreasing from oo to 1%
SolveDualsystems(L, 2, s, R,I)
foreach 1j?? of I
RefineInterachon(I5??, I, ?)
Figure 3.5: Pseudocodefor an importance driven hierarchical algorithm
pletely occluded from one another. Each link carries radiosity in one direction
and importance in the other. The set I can be computed once for each model
and stored along with it.
The algorithm refines the initial set of interactions I by transforming this
single-level network into a hierarchy of interactions. At each iteration of the out-
ermost loop, the algorithm solves for both radiosity and importance using the
current set I. These solutions, in turn, are used to guide a refinement step, which
improves the accuracy of the radiosities and importances with respect to a partic-
ular view, by adding more links to I. The iteration continues until a preset error
tolerance i? is met. The basic algorithm is described with pseudocode in Figure
3.5.
The following sections describe the steps of this refinement process in more
detail.
48
SolveDualSystems(L, 2, s, 1?, I)
repeat
GatherAndShoot(L, 2,1)
foreach top-level patch p
SweepRadiosity(L, S, p, 0)
SweepImP0rtance(2?R1 p, 0)
until convergence of L and Z
Figure 3.6: Solvingfor both importance and radiosity
3.4.2 Solving for Radiosity and Importance
For a given hierarchical system, L and 2 are found by iteratively solving two sys-
tems of linear equations. This is an extension of what is done for solving the sys-
tem in standard hierarchical radiosity (Figure 2.9). Each iteration involves gath-
ering radiosity and shooting importance between every linked pair of patches.
Radiosity and importance are then distributed up and down the hierarchy so that
every patch receives the appropriate contributions from its ancestors and descen-
dants. This process is repeated until convergenc?that is! until the difference
between the radiosities and importances from one iteration to the next becomes
smaller than some tolerance. (See Figure 3.6.)
Gathering radiosity and shooting importance is implemented in Figure 3.7.
For each interaction between two patches j ? i, radiosity is "gathered" from
j to i, and importance is "shot" from i to j. Because the importance matrix is the
transpose of the radiosity matrix, the same matrix entry p?F?? appears in both
these operations. Here Fij denotes an estimate of the actual form factor F?j?.
49
GatherAndShoot(L, Z, I)
<? (0,0)
foreach I5?j0fI
? L?+pii?ij%
Zj ? Z?t pj%Zi
?
Figure 3.7: Gathering and shooting ofradiosity and importance.
Sweeplmportance(2,I?, i,Zdown)
ifi is a leaf then
Zup ? Ri + Zj + Zdown
else
Zup ?O
for each child ii of i do
x ? (2i+Zdown)*Ai/Ai
Zup ?Zup + Sweeplmportance(2, R, i', x)
Zj <--H Zup
return Zup
Figure 3.8: Sweeping importance in a hierarchy
Once radiosity and importance havebeen transferred between linked patches,
the two quantities are distributed up and down the hierarchy so that every patch
receives the appropriate contributions from its parents and children. This was
described for radiosity in Figure 2.10, and is described for importance in Figure
3.8.
50
The two quantities are distributed in slightly different ways. The radiosity
Li of a child patch i is the sum of its emitted radiosity S?, the radiosity it receives
directly, and the radiosity of its parent Ldown. The importance of a child patch i is
the sum of its emitted importance l?j, the importance it receives directly, and an
area-weighted fraction of the importance of its parent Zdown.
When pulling radiosity and importance back up the hierarchy the situation
is reversed. The radiosity of a parent patch is the area-weighted average of the
radiosities of its children, whereas its importance is the sum of the importances
of its children.
3.4.3 Refining the Interactions
We refine any interactions whose estimated error Z?PiE?rrLj exceeds the toler-
ance ?. For F??, we use an upper bound on the error in the form factor, computed
by taking the difference between upper and lower bounds Fi5 and Eij on the ac-
tual form factor Fij, as described in the next section. If the error in the interaction
exceeds the tolerance ?, the interaction is refined by subdividing the patch p with
greater area and creating new links directly from the children of p to the other
patch. This is much like the refinement code for hierarchical radiosity in Figure
2.7, and is shown in Figure 3.9.
The procedure SubdivideAndRefine is essentially the same as the refinement
routme described by Hanrahan et al.[HSA9l]. With importance-driven refine-
ment, however, the minimum patch-size criterion that guarantees termination is
rarely necessary, since the importance of a patch always decreases with its area.
51
Refinelnteraction(5 ? i,I,c)
F?rr ? Fij -
if Zipif?rrTh' > ? then
SubdivideAndReiine(p, I ? i, ?)
I ? I u fnew links created through subdivision1
end if
Figure 3.9: Pseudocodefor refining an interaction using importance.
3.4.4 Estimating the Form Factors
We estimate the form factor F?4 from patch to patch j by taking a number of sam-
ples across patch i and averaging the point-to-disk form factors [WEH89] over all
samples.
Care is needed in estimating upper and lower bounds on Fij, since the as-
sumptions required by point-to-disk form factors may not always hold. We
choose samples on both patch i and j and compute double-differential form
factors. We set the upper and lower bounds Fj5 and Eij to Aj maxtF?i,?51 and
Aj min tFdi,djl, respectively. These "bounds1, are not necessarily strict if the sam-
ples are poorly chosen; however, they seem to work well in practice. Note that if
any of the double-differential form factors encounters an occluding object, then
the lower bound Eij is 0, and the estimated form-factor error is set to the upper
bound Fj5.
52
Maze on table			Maze on table
Medium accuracy			High accuracy
B-only			RI			Gain			B-only			RI			Gain
Total patches			46778			2983			16x
Total links			585308			27055			22x
Time (minutes)			585			13			45x
163053+			8779			 19x
2803112+			79192			 35x
2337+			76			 30x
Table 3.1: B-only versus RI refinement. Statistics for two test cases.
3.4.5 Assigning Initial Importance
The vector R determines the initial emitters of importance. Different choices of
1? allow us to use importance in different ways.
Typically, for a single, static view, we define an infinitesimal patch for the
"camera" and make it the sole emitter of importance. In this case, we also en-
sure that an initial link between the camera and the rest of the environment is
set up only if it falls within the camera's viewing frustum. This choice of 1? has
the same effect as assigning an initial importance to each patch according to its
visible projected area.
3.5 Results
Figure 3.10 shows a small maze sitting on a table. The images were computed
using importance-driven refinement. The left image is flat-shaded to show the
meshing produced by the algorithm, while the right image is displayed using a
reconstruction algorithm described in Smits et al.[SAS92]. Note the color bleed-
53
Figure 3.10: Left, a view-dependent radiosity solution. Right, the same solution after
reconstruction.
ing from the red wall of the maze to the wall behind it. This effect is most pro-
nounced on the top of the wall, as the bottom half of the wall is illuminated
mostly by light reflected from the maze floor.
The next pair of images, Figure 3.11, shows the same computed solution, but
displayed from further back. The left image shows the radiosity of each patch,
while the right image depicts the importance of each patch, divided by its area.
In the lafter image, the visible parts of the scene show up very clearly as white
because all wavelengths of light emitted from these surfaces are equally impor-
tant to the camera. The other surfaces in the room also have importance, but to
a lesser degree. The color of each surface gives the fraction of light emitted at a
particular wavelength that finds its way to the camera. Due to the red and blue
walls, red light is more important on the left side of the floor and back wall, while
blue light is more important on the right. Note how the algorithm has done less
54
Figure 3.11: The radiosity solution (left) and importance solution (right) ofFigure 3.10,
seen fromfttrther back.
Figure 3.12: The radiosity solution (left) and importance solution (right)from evenflir-
ther back.
meshing in areas with little importance.
Figure 3.12 shows the radiosity and importance of the same solution--Hdisplayed
from even further back, showing the whole environment. From this vantage
point, the meshing done on the table is extremely fine. A view-independent al-
55
gorithm would have great difficulty solving the entire environment to such high
precision. Note that no meshing at all has taken place in the room at the lower
left, which contains a complex block sculpture, since the importance of the room
and everything inside it is negligible. Note also that the wall illuminated by the
bright light in the center does have some importance, and is therefore refined to
a certain extent. The importance of indirect interactions such as these would be
very difficult for a user to anticipate in giving meshing hints to a conventional
algorithm.
Some quantitative measurements of our tests are summarized in Table 3.1.
All tests were run on an Ell) 720, a 55 MIPS machine with 64 megabytes of phys-
ical memory.
Figure 3.13 compares brightness-weighted refinement (B-only) to brightness-
and-importance-weighted refinement (BI). In running the test, we refined the
B-only solution as much as possible before running out of physical memory We
then ran a BI algorithm on the same environment until the solution appeared
to be at a similar or slightly higher level of refinement. For this environment of
moderate complexity (1002 initial polygons), the B-oniy solution was 45 times
slower, requiring 16 times the total number of patches, and 22 times the total
number of links. We expect that for environments of greater complexity, the
speedups would be even more dramatic.
We have also tried to compare timings for the more refined solution of Figure
3.10 with a B-only solution of comparable accuracy. However, the extraordinary
memory requirements for a high-accuracy global solution makes this compari-
son difficult. We let the B-only solution run for about 40 hours, at which point the
virtual memory size was well over 100 megabytes. With 19 times the number of
56
Figure 3.13: A B-only solution (left) and a comparable BI solution (right).
patches, 35 times the number of links, and 30 times the number of CPU-minutes
as the BI solution, the B-only algorithm still had not achieved a solution close to
the same level of accuracy.
As an additional test, the camera was moved inside the small maze on the
table, in the same position with respect to the small maze as it was with respect
to the large maze in Figure 3.13. The BI solution was more accurate than B-only
and took only 0.5% of the time to compute.
Finally, in order to get a sense of how well the BI algorithm might perform
for a moving camera, we started with the solution for Figure 3.13 and then turned
the camera 10 degrees. Solving to the same level of accuracy took an additional
3?21 minutes, about a quarter the time for the original view.
The initial-linking time, which is identical for BI and B-only refinement, is
not mcluded in the statistics in Table 3.1. For the 1002-polygon environment used
in these tests, initial linking took 75 minutes and produced 18314 links. Note,
57
however, that since initial linking is independent of the level of refinement, ma-
terial characteristics, or view, it can performed once and stored along with the
model. Our current implementation of initial linking uses a brute-force algo-
rithm that checks every pair of patches to determine visibility; however, this step
could be easily optimized using a more sophisticated visibility testing scheme.
3.6 OtherUsesforImportance
So far we have mainly discussed using importance to accelerate a single view
of an environment by estimating the errors contributed by interactions in differ-
ent parts of the environment. There are other ways to use importance that make
it more flexible as well as approaches that give guaranteed bounds on the error
caused by each interaction.
3.6.1 Initial Importance
There is nothing, however, that limits us to a single, static view. For example, im-
portance could be used to speed the radiosity solution for an animation or walk-
through of an environment by making every patch that is visible at any time an
emifter of importance. Similarly, importance could be used when the scene is
visible from more than a single viewpoint, such as in stereoscopic views or in
stage design applications [DSG9i].
Finally, for some applications it may not be necessary for everything visible
to be important. For example, perhaps we would like to study a particular object
from all angles, such as a sculpture in a museum. If we are unconcerned about
the environment itself except in how it contributes to the illumination of this ob-
ject, we could make the object itself the sole emitter of importance.
58
Note that in all cases, the overall magnitude of I? is arbitrary. However, the
error tolerance F? must be scaled appropriately.
3.6.2 Bounding Error Conservatively with Importance
The previous derivation relied on estimates of radiance and importance and was
not exact. If we are willing to do more work, we can get conservative bounds on
the error each link contributes to a linear functional. This can be done by com-
puting upper and lower bounds T and L on the exact radiance functions across
the patch. Using these upper and lower bounds, we can determine bounds on
the error in the environment as described by Lischinski et al.(LSG94].
Upper and lower bounds on the radiance on each patch in the environment
give information as to which patches potentially contain the most error. This dif-
ference can be used to estimate where more work needs to be done. If the differ-
ence is caused by incorrect transfers of radiance to the patch, this estimate will
be effective. However, the difference could be a result of large differences on all
patches contributing to the receiver. In addition to not knowing the cause of the
error, the effect of that error on the rest of the environment is unknown as well.
These bounds give no way of determining how the errors propagate, and they do
not give precise information as to how much error in the environment is caused
by the error in a given transfer.
We want a way to determine the global effect on the solution error introduced
by a link error. The effect of link error is magnified due to the error being propa-
gated throughout the environment. We must find a way of determining how far
the error is propagated in order to determine the true effect of the link error on
the solution error. We can do this efficiently when we use i-norm for bounding
59
the error.
The i-norm of T --H L can be expressed as
T?LiIl?7rSITi?LiI*Ai			(3.9)
i=1
where Aj is the area of patch i. The absolute value is not needed since T > L.
Without the absolute value, the i-norm becomes a linear functional. We wish to
minimize the value of this linear functional over the environment. The impor-
tance driven radiosity approach presented by Smits et al. [SAS92] relates linear
functionals to the solution process of hierarchical radiosity. It was shown that
importance could be used to approximate the effect of a single interaction on a
linear functional of the environment. Here we will use importance to find an ex-
act relationship between the i-norm and the error on a link.
In order to apply the concept of importance to this problem we must intro-
duce some notation. The one-norm linear functional can be written
vTB			(3.10)
where Vi is the area of patch i.
Importance is the solution to a related system of transport equations. V be-
comes the emission vector for importance in the system. The solution of the sys-
tem gives the importance of each patch. The importance gives the total effect of
radiosity emitted from a patch to the linear functional. In this case, the product of
the importance and the emission of a patch gives the total effect of that emission
to the i-norm.
Let N = 1 --H M and N = 1 --H M. Now the set of transport equations for
radiosity and importance for the lower bounds can be written
NL = E			(3.11)
60
NTz = v			(3.12)
The difference AM = N --H N = M --H M is the error in links. It is this quantity
which we want to relate to the i-norm of the error bound. AM can be used to
relate the two matrices of bounds as follows: first
so
NL = (?N-AMff=E
(3.13)
NT=?E+ AML.			(3.14)
Now we derive an expression for the i-norm of the error which is in terms of
the errors in interactions. This derivation will make use of a relationship between
importance and radiosity as shown in [SAS92]
vTT=zTE.			(3.15)
Given these relationships, the i-norm of error bounds vT(T?L) can be expressed
in terms of link error as follows
Now using equation 3.14 we get
=			vTT?uL			(3.16)
=			vTT --H zTE			(3.17)
=			(?NTz)TT --H zTF			(3.18)
=			zTNT - zTE			(3.19)
vT(T?L) = zTNT?zTE
= zT(E t AML) - ZUE
= zTAML+?zT(E?E)
(3.20)
(3.21)
(3.22)
61
This expression relates the i-norm of the bounded error to the link errors
and the emission errors. The term ZTAML is a sum over all links. The values
summed are the link error weighted by the importance of the receiver and the
radiosity of the source. This gives exact information on the effect of each link to
the bounds on global error, as opposed to the approximations given earlier in the
chaptet
We can now use the conservative bounds to determine a bound on the i-norm
of the error, and we can know exactly where in the environment we will get the
most benefit from doing work, as we know which interactions contribute the
most to the error.
This result also shows that importance can be used with arbitrary functionals
to achieve conservative bounds. The use of importance as described earlier in the
chapter did not guarantee bounds on the result.
This approach is useful when guaranteed bound on the global error is nec-
essary. Computing the upper and lower bounds on the exact solution is com-
putationally costly, however, and the approximations used earlier work well in
practice. The discussion of non-conservative bounds and other changes can be
found in Lischinski et al.[LSG94].
Chapter 4
Clustering
4.1 IntroductionlMotivation
Currently, the best radiosity algorithms are analogous to linear-time algorithms
for charged particle simulation [Gre88j. These algorithms work by clustering
particles together so that the mutual effect of well-separated collections of parti-
cles can be computed through a single interaction. Hanrahan et al. [HSA9lj used
a similar strategy to reduce the number of interactions needed for a radiosity so-
lution by imposing a hierarchical structure on each surface in an environment
and allowing each interaction to occur between the appropriate levels on each.
This approach works well when the number of initial surfaces is small, as hier-
archical radiosity (HR) algorithms can "subdivide" large surfaces into smaller
ones, but cannot "group" smaller elements into larger ones. The initial linking
phase of HR must check all pairs of initial surfaces to see if they interact, there-
fore this inability to group surfaces together means that the algoritlim is 0(s2)
in the number of initial surfaces. Furthermore, as complexity in environments is
often a result of replacing large surfaces with many small surfaces, this inability
62
63
to group objects together is a major obstacle in conventional HR algorithms.
Several methods have been developed to increase the efficiency of FIR algo-
rithms for complex environments. Global visibility [TH93j is effective for com-
puting the initial links for environments in which only a small number of "cells"
see each other. Many environments, however, also contain large collections of
objects that are mutually visible. In such cases, global visibility algorithms do
not suffice.
Importance [SAS92j is another method that reduces computation in areas that
have little or no noticeable effect on the surfaces of interest. This approach is ef-
fective once the initial interactions have been computed, but does nothing to re-
duce the time of the initial linking phase. Importance driven hierarchical radios-
ity stands to gain even more from clustering than standard hierarchical radiosity.
Because it attempts to compute the radiance very coarsely in some regions of the
environment, some means of clustering is needed if it is to do less work than that
required for the initial linking.
An approach to clustering was developed by Rushmeier et al. [RPV93] in
which the effect of complex groups of surfaces is approximated by simpler rep-
resentations resembling BRDF's obtained through Monte Carlo sampling. One
disadvantage of this approach is that it is not automatic; clusters must be spec-
ified and simplified in advance. Also, no hierarchy is maintained, so the coars-
est representation used by the algorithm is that of the selected clusters. Another
clustering approach was proposed by Kok[Kok93], which extended progressive
radiosity so that patches could transfer energy to groups of surfaces at once. Nei-
ther of these approaches analyze the error of the approximations used for the
transfers.
64
Patch links
a-link
`3-link
Figure 4.1: Two collections ofobjects can interact at different levels: a) conventional HR
links, b) an ?-link requiring linear time and space, and c) a P-link requiring constant
time and space.
To illustrate why clustering is important, consider a simple configuration
consisting of two chairs, each containing 100 surfaces. HR would check the
10,000 potential interactions between the two chairs and would create in the
neighborhood of 2,500 patch links (Figure 4.1). (It is assumed that half the poly-
gons are not facing each other.) However, if the chairs are well separated, then
it is unlikely that the energy transfer between them will have a significant effect
on the illumination of either. If we can guarantee this throughout the entire so-
lution, we can gain efficiency by coarsely approximating the transfer of energy
between them. We shall present a new strategy for linking that would handle
65
this configuration by creating a single link with a cost comparable to either 1 or
200 conventional links, depending on the separation of the objects. By reducing
the total cost of the links created between two collections of objects, we reduce
the algorithmic complexity of the 0(32) initial linking step in HR.
4.1.1 Overview
Hierarchical algorithms for radiosity have three components: 1) a hierarchical
description of the environment, 2) a criterion for determining the level in the hi-
erarchy atwhich two objects can interact, and 3) a means of estimating the energy
transfer between the objects.
As our criterion for interaction, we use bounds on the potential error in the
transfer of energy between two objects. We first describe hierarchical radios-
ity using this approach, and then present two efficient techniques for bounding
the potential transfers between clusters. The first technique is a fairly accurate
approach which we call a-linking. The second is a faster but less accurate ap-
proach which we call P-linking. Corresponding to these methods for determin-
ing bounds are two ways of estimating the transfer of energy between clusters.
We also describe the method used to create the hierarchy of clusters, which is the
starting point for the cluster-based linking strategies. Finally we give theoretical
bounds on the complexity and results of our implementation demonstrating dra-
matic speedups over conventional HR in complex environments.
4.2 Mathematics
In this section we describe two strategies for computing bounds on the energy
transfer between clusters; one requiring linear time and space, and one requir-
66
ing constant time and space. The strategies are derived by systematically intro-
ducing approximations into the exact expression for energy transfer. The first
level of approximation results in a type of link we call an oL-link. By introducing
further approximations we produce ?-links. Although these bounds are coarse,
they often suffice for a very large fraction of the interactions.
4.2.1 ?-1inks
The transfer of energy from a cluster of patches 8 to a point can be computed by
summing over the n patches in the source cluster
Alsi k(x,y)L(y)dy.
(4.1)
Since k(x, y) is a function of the position and orientation of the receiver as well
as the source and of the visibility between them, computing the transfer of radi-
ance between two clusters requires evaluating L(x) at least once on each of the
m receivers, resulting in O(mn) work.
Applying the same approach used for HR to two clusters is relatively expen-
sive. To bound the transfer between the clusters with this approach, we can use
the maximum value of k(x? y) over pairs of patches to bound the transfer. The
maximum radiance over all receivers due to the source cluster can be bounded
by
maxL(w) ? m3a.xffi [k]m,& ?? L(y)dy. (4.2)
i=1
As the maximum value of the kernel function must be computed for each pair of
patches, the time complexity is O(rnn) (Figure 4.2).
To improve the time complexity, we split the kernel function into two simpler
67
Receiver
k
Source
Figure 4.2: Straightforward approach to bounding error. Each link represents a maxi-
mum valuefor k.
functions kr(x, y) and ks(x, y) given by
Then
kr(x,y)			p(x)cosOi
ks(x,y)			cos O2
x--Hy2.
k(x,y) kr(x,y)ks(x,y)vis(x,y).
(4.3)
(4.4)
(4.5)
Each of the new functions depends on two p0int5, but each requires informa-
tion about the orientation of only one surface. k? accounts for the projected area
and the reflectivity of a differential area around x and ks corresponds to the dif-
ferential solid angle subtended by a differential area around y. We now show that
the above separation of the kernel can be used to conservatively approximate the
transfer of energy from the source cluster to the receiver cluster. This is done by
bounding the energy that reaches the receiving cluster from the source, and then
determining how much of this energy potentially reaches each receiving patch.
68
Figure 4.3: More efficient ?-1ink between two clusters.
Z&?i fi kr(x, y)ks(x, y)vis(x, y)L(y)dy
We can apply this idea to obtain upper bounds on the transfer between the
two clusters. The quantity ks(x, y) can be maximized over all points in the receiv-
ing volume B(1Z) making ks(x, y) independent of any particular receiver. The
L(y)dy.			(4.9)
[k?]xB(s) ? [kS?x,s
t=1
Finally, assuming visibility is always 1, we arrive at the conservative upper
bound
(4.8)
(4.7)
x,B(S) Zj=1 `s ks(x, y)vis(x, y)L(y)dy
[krlx,B(s)li [kS]x,st f? vis(x,y)L(y)dy.
(See Figure 4.3).
We now use these pieces of the kernel to bound the energy transferred to a
point x. First we replace kr(X, y) by the maximum value it attains over all yin the
volume B(S) containing all the source patches. Then, for each of the sources, we
replace ks(x, y) by the maximum value it attains over all y on the source patch.
More precisely, these steps produce the following bound:
L(x) =
(4.6)
Receiver
kr
MM
Source
69
factor k? (x, y) can now be maximized over each receiving patch separately. Thus,
the energy transferred can be bounded as follows
max L(x) < m?a.x [kr]R,,B(s)? [kS]B(R),st			L(y)dy.
St
This expression can be separated into two pieces by defining
(4.10)
Lsmax			(4.11)
Then
))?1[kS1B(?),Si Si L(y)dy.
maxL(x) ? max[k?1Ri,B(s)Lsmax			(4.12)
Maximizing over all of the receiving patches in ? requires 0(m) work. There-
fore this bound can be computed in O(rn + n) instead of O(rnn).
This upper bound is one component needed to apply the mechanisms of the
HR algorithm to clusters, as HR requires bounds on the error in the interaction
between two objects. These bounds can be computed from upper and lower
bounds on the transfer. We have given upper bounds; a trivial lower bound of
zero can always be used. This conservative lower bound is feasible because we
do not know the visibility between all points in the two clusters. It is possible that
each patch is partly occluded and the visibility calculations are expensive. Now
we have a bound on the error in the interaction between two clusters, which al-
lows us to determine whether or not a given interaction is accurate enough. This
bound on clusters determines if an Q-link is an acceptable approximation for the
interaction.
In addition to the error bound, we require an estimate for the energy trans-
ferred to each receiving patch across the a-link. We do this by computing an es-
timate as well as the bound for the kr(X, y) associated with each receiving patch
70
and the k8(x, y) associated with each source patch. Let
(f?A,B ?V?xy?%Af(X? Y)
Now the radiance on receiving patch j can be estimated by
(4.13)
Lj (kr?Ri,B(s) Ls?g			(4.14)
where
Ls?g			(kS?B(?),si L(y)dy.			(4.15)
(4.16)
Each average is bounded above by the maximum transfer and therefore falls
within the error bounds hence, it is accurate enough for the algorithm; in fact
it is quite likely far more accurate than the conservative bounds indicate.
We now show that the a-link clustering approach has time and space corn-
plexity of O(s log s). First assume we have s initial patches stored in a hierarchy
of clusters with each cluster containing two smaller clusters or, at the leaves, a
single patch. This structuring results in a binary tree with a depth of log s. The
total number of clusters at level d is 2d and each cluster at level d has s/2d patches.
Now, following [HSA91], we assume each cluster is linked to a constant c1 num-
ber of other clusters resulting in ci2d a-links on level d. We can also assume that
each cluster is linked to other clusters at approximately the same level in the tree.
An a-link between two clusters requires space and time proportional to the size
of the clusters. Summing the space/time costs for each of the log s levels of the
hierarchy gives
logs			logs
? numAinksd * link?costd			? c12dc2??			Cs logs.
0			d=O
(4.17)
71
Therefore, rather than the 0(s2) work required for a direct approach, cluster-
ing using ?-links gives a time and space complexity of O(s tog s).
4.2.2 4-Links
For many transfers, the previous technique is still too expensive. In many envi-
ronments large numbers of interactions are insignificant and can be treated very
coarsely or ignored altogether [LTG93]. In this section we introduce a more effi-
cient but cruder method for bounding the interaction between two clusters. We
do this using a very simple bound on the kernel k(x, y) which we denote kd(x, y),
k(x,y)?kd(x,y)			1			(4.18)
IIx--Hy112
This bound requires no knowledge about the orientation of the surfaces. We can
bound the transfer of energy between two clusters by replacing k by the maxi-
mum value of kd over all points x in the receiving cluster and all points y in the
source cluster. Thus
max L(x) ? [kdlB(Th,B(s)f L(y)dy.
(4.19)
This is a worst-casebound on the transferbetween two clusters; it is only achieved
when all surfaces are as close as possible to the other cluster, and each source di-
rectly faces all the receivers. The virtue of this bound is that it requires no knowl-
edge of the surfaces. If the integral of the radiance over all the source patches
has been computed in advance, which can be done during a sweep operation
described in Section 4A, then the bound can be computed in constant time. As
with the a-links, we can use the average value of kd over the two clusters as an
estimate the energy transferred between them. Again, since the average value
72
Receiver
kd
Source
Figure 4.4: P-iink between clusters.
lies within the range given by the maximum transfer and a minimum transfer of
zero, it meets the error tolerance.
The time and space complexity of a system with ?-links is equivalent to
the complexity of the standard hierarchical radiosity method. As with linking
patches, the error term decreases rapidly with distance, so each cluster will be
linked to a constant number of other clusters. Since the number of clusters is lin-
ear m the number of initial surfaces, and the cost of a P-link is independent of
the size of the clusters, the complexity is 0(s) where 5 is the number of initial
surfaces.
4.2.3 Strategies for linking
We have described two different approaches for linking clusters. The more accu-
rate clustering technique tends to be too expensive for many clusters in a large
environment. The faster clustering technique is too coarse to eliminate many
of the negligible interactions. Our approach is to exploit the strengths of both
strategies. If the coarse approach does not produce an acceptable bound then
the more accurate technique is invoked. If neither are accurate enough, then the
73
clusters are recursively subdivided. Only when the level of individual patches
is reached does the algorithm resort to conventional HR.
4.3 Other Linking Criteria
Norms provide a measure of the "size" of a function. We shall use two different
norms to quantify error in the energy transfer between two patches or clusters.
In the previous section we computed a bound on the maximum possible differ-
ence between the approximated radiance and the actual radiance over all points
on the receiving object. This corresponds to computing a bound on the oo-norm
of the radiance function on the receiver due to the source. More formally, if L
is the exact radiance function over the receiver, and L is the computed radiance
function then
IIL--HL oo=mx??7z4L(X)?L(X))?M			(4.20)
where M is one of the previous bounds on the transfer.
The oo-norm gives a bound on the variation between the computed radiance
and the actual radiance. Another useful bound is the energy of the difference
between the computed and actual radiance functions due to a link. This bound
corresponds to the i-norm. More formally,
IIL--HL i?7rj?ILY)--HL(x)dx			(4.21)
The i-norm can be bounded by weighting each term of the receiver by the area of
the receiver and summing instead of finding the maximum value. We can write
this bound for Q-links as follows
IIL--HL ?` Aj [kr]R3,B(s) Zj? ([kslB(R),s f L(Y)du)
(4.22)
74
where Aj is the area of receiver j. Coarse bounds and patch-to-patch bounds are
computed similarly. In a hierarchical solution both norms are easily accommo-
dated, but the norm used affects the subdivision strategy After performing a lo-
cal pass to display a solution, the error at each pixel is important, which implies
that the oo-norm is most appropriate. However, during the global propagation
of energy throughout an environment, it may be more useful to minimize the
i-norm of the error. Intuitively, this is because large surfaces will often have a
much more significant effect on the local pass than the smaller ones.
Another useful norm is obtained from the importance weighting found in
Smits et al.[SAS92]. This norm is very similar to the i-norm; rather than weight-
ing each receiver by the area, the receiver is weighted by its importance.
4.4 Implementation
Most ray tracing acceleration techniques collect nearby objects together into
groups. Some of these, such as octrees and hierarchical bounding volumes, form
a natural hierarchy The same hierarchy can be used to identify clusters. For sim-
plicity, a useful constraint is that each object appear in only one cluster. Other-
wise, some method of eliminating copies of objects is needed when computing
links to prevent an object from contributing to some other cluster twice. For these
reasons we chose to use bounding volume hierarchies (GS87] to generate our hi-
erarchy of clusters.
The HR algorithm can be combined with a bounding volume hierarchy to
perform clustering fairly easily The bounds on the error given earlier fit natu-
rally into the brightness-weighted refinement scheme of Hanrahan et al. [HSA91].
This approach first links the surfaces together, then goes through several iter-
75
Create-?-link(Src,Rec)
o?-link T
vis = EstimateVisibility(Src,Rec)
foreach patch? in Src
T.SrcMax[i] ? [kS]B(Rec),s
T.SrcAvg[i] ? (kS\)B(Rec),s
foreach patch j in Rec
T.RecMax(j] ?
T.RecAvg[j] ? Vis *
if a-Bound(T) ? c then return TRUE
else return FALSE
Figure 4.5: Pseudocodefor creating an &-iink.
ations of energy propagation and refinement of the system. Each iteration is
done with a smaller acceptable error tolerance, until the final error tolerance is
reached. This approach is a variation on the `,multigridding', method for itera-
tively solving systems of linear equations.
Very few modifications need to be made to the bounding volume data struc-
tures to implement clustering. For the more accurate clustering approach, we
must loop over all of the patches in the bounding volume. Each link will con-
tain four arrays that hold both the maximum and average values of k8 for all the
source patches as well as the maximum and average values of k? for all of the re-
ceiving patches. In our implementation the maximum values are approximated
by taking a fixed number of jittered samples on the patch and bounding volume,
then finding the maximum value of the kernel at these points. This method does
not guarantee finding a bound, but in practice it produces a good approximation.
76
a-Bound(T)
SrcError <- 0
MaxError #. 0
foreach patch in T.Src
SrcError <? SrcError + T.SrcMax[i] * Aj *
foreach patch I in T.Rec
MaxError <? Max(MaxError, T.RecMax[j] * SrcError)
return MaxError
Figure 4.6: Pseudocodefor computing error bounds on an &-link.
We weight the average value of k? for each patch by an estimate of the visibility
between the two clusters. The pseudocode in Figure 4.5 shows the steps needed
to compute an a-link between two clusters.
Procedure a-Bound(T) computes the bound on the maximum amount of en-
ergy transferred across a-link T. We have estimates of the maximum values of kS
and k? stored in the link. The integral of the radiance function Lt over patch i for
constant patches is the product of the radiance Lj and the area Aj of the patch.
The procedure for computing the error incurred by an a-link between clusters is
shown in Figure 4.6.
In addition to computing the error on the links, energy must be transferred
between the two clusters. During the gather step of HR, energy is transferred
across each link, including all cluster links. The gather across an a-link T is pre-
sented in Figure 4.7.
For the coarse approach to clustering using P-links, we store the sum of the
areas of all the patches in the bounding volume. We also store the radiance of the
bounding volume for use as a source, and the irradiance striking the bounding
77
o?-Gather(T)
SrcRad ? 0
foreach patch in T.Src
SrcRad ? SrcRad + T.SrcAvg[ij * Aj *
foreach patch j in T.Rec
Lj ? Lj + T.RecAvg[ij * SrcRad
Figure 4.7: Pseudocodefor gathering energy across an o-link.
volume for use as a receiver. If the radiance on each patch of the receiving clus-
ter were updated each time a link was used in a gather, then gathering through a
single link would require linear time, resulting in worse complexity bounds for
coarse links. As in FIR, after transferring energy across the links the radiances
must be swept to each object's parents and children. For clusters, this sweep is
a slightly different from the sweep described in HR since each cluster holds the
irradiance 1 incident upon it, and the radiance leaving it. Irradiance is pushed
down the bounding volume hierarchy to the patches. There it is sent through
the BRDF of the patch and the resulting radiance is pushed down the patches
as in FIR. Pulling radiance up from the children to the parents is the same as in
FIR. Each parent receives the area-weighted average of the radiances of its chil-
dren. The pseudocode in Figure 4.8 is used when solving the linear system to
redistribute the energy gathered across the links. HSweep is the conventional
FIR sweep procedure.
The P-links between clusters are computed in a similar fashion to the ?-links.
They are easier to compute, however, as no looping over the contents of the
bounding volume is required. Also, the maximum value for kd is the inverse
78
Sweep(Cluster,L?0??)
1down ? 1down + 1Cluster
if Cluster is a leaf then
foreach patch Rj in Cluster
Lj ? Lj + Idown * Pi
Lup <7Lup+ HSweep(1??) * Aj/Ac1uster
else foreach Child in Cluster
Lup ? Lup + Sweep(Chlld,I?0wn) * AChild/ACluster
LCluster ? Lup
return Lup
Figure 4.8: Pseudocodefor sweeping energy throughout a hierarchy
Create-?-link(Src,Rec)
P-link T
vis ? Estimatevisibility(Src,Rec)
T.Max ? [kdlB(Rec),B(src)
T.Avg ? vis * (kd?B(Rec),B(src)
if fl-Bound(T) < ? then return TRUE
else return FALSE
Figure 4.9: Pseudocodefor creating a p-link.
of the square of the minimum distance between the clusters, which can be com-
puted with little work. Each fl-link stores only the maximum value and the best
computed value of kd for the two clusters.
The error in a fl-link between clusters is also straightforward, following di-
rectly from the bounds on the transfer and the previous discussion for higher
accuracy links. Computing the error on a fl-link T can be implemented as shown
79
P-Error(T)
return T.Max * LT,Src * AT.Src
Figure 4.10: Pseudocodefor bounding error on a (3-link.
(3-Gather(T)
JT.Rec = ii.Rec + T.Avg * LT.Src * AT.Src
Figure 4.11: Pseudocodefor gathering energy across a (3-link.
in Figure 4.10.
During the gather phase of the algorithm, each (3-link transfers energy to the
receiving cluster. It is stored there as irradiance until the sweep phase when it is
scaled by the BRDF. This is shown in Figure 4.11.
The above procedures allow the clusters and patches of the environment to
be linked togethet When two bounding volumes are to be linked it is possible
to first check for complete occlusion. If the two bounding volumes are not com-
pletely occluded, then the (3-bound between the clusters is checked. When that
is not sufficiently accurate the a-bound is checked. If this is still not accurate
enough, then the children of the larger cluster are recursively refined against the
smaller cluster. If neither of the clusters has children, then the patches in both
are refined against each other using the patch refinement from HR. This strategy
is summarized by the pseudocode in Figure 4.12.
The clustering algorithm begins by using the top-level bounding box both as
source and as receiver, which triggers the recursive refinement. Complete occlu-
sion is checked by checking the convex hull of the shaft between the two clusters
80
RefineCluster(Src,Rec)
if not Visible(Src,Rec) return
else if Create-(3-link(Src,Rec) return
else if Creat&o-link(Src,Rec) return
else if Children(Src) or Children(Rec)
Refinechildren(Src,Rec)
else foreach St E Src
foreach I?? ? Rec
HFefine(S?, Rj)
Figure 4.12: Pseudocodefor refining two clusters
agamst convex polygons in the environment. This is a conservative check, it can
consider two completely occiuded bounding volumes visible, but it will never
decide that two cluster are occluded if some point on one is visible from some
point on another.
4.5 Results
FIR has a complexity bound of O(s2 + p) where s is the number of initial surfaces
and p is the number of resulting patches. The 0(s2) term comes from computmg
the initial links between all initial surfaces. Using clustering we do not need to
check all pairs of surfaces even if they are mutually visible. Clustering replaces
the expensive initial linking step with an algorithm that is O(s log s) in the num-
ber of initial surfaces. This results in an algorithm with an overall complexity of
O(s log s + p). The sections that follow report results for several different envi-
ronments and refinement strategies.
81
4.5.1 Link Complexity
We first show how the cost of the links grows with the size of the environment.
As a first example, the inside of a sphere was tessellated into triangles. Each tri-
angle was given the same emissive power and reflectance. This allowed us to
vary s over a large range of values without really changing the geometry Also,
because every surface can see every other surface, it is easy to determine exactly
how many initial links would be needed without clustering. It is also a challeng-
ing environment for clustering because every cluster overlaps with several other
clusters. We assign a cost of 1 for patch links and ?-links. The cost for ?-links is
m + n where ni is the number of patches in the receiving cluster and n is the num-
ber of patches in the source cluster. The total link cost is the sum of the link costs
for each link computed for the environment. The model was refined using the
oc-norm criterion on the links. The graph in Figure 4.13 shows the link cost for
five tessellations ranging from 128 patches to 32768 initial patches, as well as the
function y s log s. The costs closely match the function. For the largest tes-
sellation a conventional FIR algorithm would need to create and store over one
billion initial links before it could compute a solution. With clustering the links
cost only about 0.2% of this, a reduction in both time and space of approximately
500.
4.5.2 Refinement Complexity
We now show how a solution refines over time both with and without cluster-
ing. We start with an environment of 4170 initial surfaces and refine the solution
from the initial linking stage through eight iterations of successively smaller er-
ror tolerances. In this example we used a 1-norm bound on the error for each in-
82
0
-J
le+07
le+06
100000
10000 ?
slogs
Trials
1000
100			100000
1000			10000
Initial Surfaces
Figure 4.13: Link cost for environments of different sizes compared with thefinction
slog(s).
teraction. A graph of the times for the global solution as the accuracy increases
appears in Figure 4.14. These solutions were computed on an HP 755 with suffi-
cient memory to prevent swapping. Images of the flat shaded patches are shown
for several different levels of accuracy in Figure 4.15. The long thin polygons in
the model do not create ideal clusters because of the large variation in position.
Neither algorithm grows linearly with the number of iterations, the algorithms
are simply linear with the number of resulting patches.
With clustering, FIR becomes more of a progressive algorithin since initial re-
sults can be seen relatively quickly. Without clustering it took 583 minutes for
the first solution to appear. Clustering reduced that time to 5.35 minutes, a speed
83
900
800			With Clustering
700			Without Clustering --*?
600			e			f
500
400
300
200
`00
0a
3			4
Iterations
0			1			2			5			6			7			8
Figure 4.14: Timefor solutions of increasing accuracy, both with and without cluster-
ing.
up of over one hundred times. Both solutions refine at about the same rate once
they are linked. After eight iterations, the clustering solution resulted in 114554
patches. Without clustering there were 114410 patches. These numbers are very
close because subdivision can only occur by refining a patch link, and this is done
the same way in both algorithms.
If an environment were refined to a very high level of accuracy, then all in-
teractions would be at the patch level, eliminating all the cluster links. It is not
feasible to refine a solution to this level of accuracy except for very simple en-
vironments (N100 initial surfaces) due to the time and memory requirements.
This does, however, give some insight into what clustering is doing. Clustering
84
Figure 4.15: Solutions at different accuracies. With clustering (a-e), without clustering
(frg).
amortizes the cost of initial linking over many refinements and except for simple
environments, the amortization continues well beyond the required levels of ac-
curacy. This amortization has some cost associated with it, but it is compensated
for by reduced link maintenance; that is, there is no need to maintain the 0(s2)
initial links through all the early iterations where high accuracy is unnecessary
4.5.3 Importance
We also tested our algorithm using importance weighted refinement. We used
the same environment as before, with a view of the group of chairs near the stairs.
The upper image shown in Figure 4.16 was created with the view-dependent lo-
cal pass that will be described in Chapter 5. The global solution used for the lo-
cal pass is shown below the image. Although our implementafion does not split
85
Figure 4.16; Image resultingfrom local pass (above) and solution usedfor the local pass
(below).
surfaces when they interpenetrate, hierarchical radiosity can be extended to han-
dle these problems [LTG93]; these techniques are independent of the clustering
algorithm. The global pass with clustering took about 3 minutes for the initial
86
solution and 53 minutes for the entire global solution, creating 27318 patches.
Without clustering, the initial linking took 728 minutes and the entire global pass
took 790 minutes. (In this trial, estimates of form factor and visibility error were
computed more carefully in order to make the local pass more effective, causing
the initial linking to take longer than it did in the previous example.)
4.6 Better Hierarchies
The hierarchy of clusters presented earlier were fairly straightforward modifica-
tions to traditional bounding volume hierarchies. It is beneficial to build SD hier-
archies in the spirit of [AK87], to bound direction of the patch normals as well as
position. This can be used to obtain tighter bounds on the transfers since patch
orientations are constrained. In the implementation described here, two clusters
with even one pair of patches facing each other may be linked with a relatively
expensive ?-link. A SD hierarchy seems to reduce this problem.
Chapter 5
Local Pass
5.1 Introduction
Once a global solution has been computed, we generally want to display the re-
suits. This usually involves smoothing the constant radiances computed for each
patch so that the artifacts of the patch boundaries are not visible. This can be
done by linear interpolation or Gouraud-shading the surfaces; however this is
very prone to artifacts in hierarchical systems [SAS92]. A better approach is to
use the techniques of Lischinski et al. [LTG93j, creating a separate mesh con-
taining discontinuities for display purposes, and then using the information ob-
tained by the global solution to recompute some parts of the illumination on each
surface. We used an approach based on this method as well as an image based
method proposed by Reichert [Rei92]. These approaches are illustrated in one
dimension in Figure 5.1.
87
88
Piecewise Constant
Linear with Discontinuities
Piecewise Linear
Image Based Reconstruction
Figure 5.1: Four different approaches to reconstruction.
5.2 Mathematical Background
Once the global solution has been computed we have an approximate radiance
function L defined over the environment. We now can use this to create an im-
age. We do this by determining the radiance coming through each pixel of an im-
age. Using ray tracing techniques, it is easy to determine a point x on a surface
that contributes to a particular pixel. (See Figure 5.2.) We then need to recon-
struct the radiance function L(x) at that point using the global solution L. The
simplest way to do this is to let L(x) L(x). The problem with this approach is
that the global solution is almost always computed with basis functions that do
not accurately capture all the illumination details across a patch. The best exam-
ple of this is the widely used constant basis function over a patch.
A better solution is to evaluate the integral equation for global illumination
at a single point. This can be expressed as
f(x) = g(x) + fHin p(x)L(x') cos(O')dw'			(5.1)
89
ImagePlane
Global Solution
Figure 5.2: The local pass uses the radiances L computed during the global solution to
determine the radiance at a pixel.
or
f(x) = g(x) + fo(x)L(?')cos(O')vis(x,x') cos(O)2dA (5.2)
Ix--Hx
Evaluating this integral directly can be expensive, because it potentially involves
determining the contribution to the radiance at x from each surface in the envi-
ronment. More importantly, the integration requires many expensive visibility
checks. This is the approach used by Reichert [Rei92].
90
Figure 5.3: patches contributing to point x in a hierarchy.
5.3 Reconstruction Using Hierarchical Radiosity
Solution
Hierarchical radiosity makes this view-dependent reconstruction much more ef-
ficient. Given x corresponding to a particular pixel we can quickly determine
which patches in the hierarchy x is contained within. The links connected to each
of these patches give the set of surfaces that contribute to this point. This reduces
the number of surfaces that need to be checked to a much smaller set. Also, each
patch linked to x is linked at the appropriate level of detail, keeping the number
of links at each point comparatively small (Figure 5.3).
The reconstruction process at point x is now fairly straightforward. Both the
91
form factor and the visibility need to be computed between each source and x.
The form factor between a polygon and a point has a simple analytic solution.
The analytic form factor must be used to eliminate the artifacts that result from
approximations to the form factor [Rei92]. Visibility can be estimated by shoot-
ing rays from x to the source and looking at the fraction that hit. Since visibility
error caused subdivision in the global solution1 partially occluded sources will
have been subdivided more than unocciuded sources, so it is easier to estimate
visibility accurately.
This reconstruction technique is still impractically slow. In order to recon-
struct shadows accurately, many rays mustbe cast towards the sources that cause
the shadows. If the same number of rays are east towards each source, it is pos-
sible that several thousand shadow rays will need to be shot for each pixel in the
image.
A similar problem arises with form factor computation. Many of the inter-
actions are very small, and have minimal error across the receiving patch. The
result is that far too much work is done for interactions that contribute little to
the resulting image. Fortunately, each of the links connected to x has estimates
of the visibility, form factor and error. This information was expensive to com-
pute and should be used to make the reconstruction more efficient. What follows
is a way of using Monte Carlo techniques to greatly reduce the amount of work
needed to reconstruct the radiance at x without introducing bias.
5.3.1 Unbiased Reconstruction
The first step in the process is to determine which interactions need to have the
most work done in the reconstruction phase. This can be done by looking at the
92
error St On each link. The absolute error will not work however, as the amount
of acceptable error is dependent upon the intensity of the radiance at point x.
For perceptual reasons, we want to do the same amount of work in dark areas of
the image as m bright areas because of the eye's ability to distinguish brightness
variations better at low level intensities. We can normalize the error of each link
by dividing by the sum of the errors on all the links contributing to x to get the
relative error 5 E [0,11 for each link.
si =
zisj
(5.3)
This relative error can be used to determine the amount of work to be done for
each interaction. For some interactions, we will compute either the analytic form
factor, the visibility, or both. For the others we will use the values stored in the in-
teraction. We will make this choice probablistically. In order to eliminate the bias
that would be introduced by using the incorrect value on the links, we compen-
sate for the error between the incorrect value and the newly computed value by
adjusting the newly computed value as in Arvo et ai.[AK9O]. We shall describe
this approach briefly.
Assume we have an estimate z of some quantity z. With probability p we will
compute the "exact" value z. Exact means that the expected value, E(z), of z is
Z. If we use z with probability p and z with probability 1 --H p the expected value
E can be written
E = E(pz+(1--Hp)z?			(5.4)
= pE(z)+(1--Hp)z			(5.5)
= pZt(1?p)zffl			(5.6)
In general, this approach will give the wrong expected value. This is known as a
93
biased approach. This bias can be corrected for by using a corrected term instead
of z whenever the value is recomputed. This corrected term can be written
z--H(l--Hp)z			(5.7)
p
Now the approach is not biased, since the expected value is
E
P)Z)t(1?p)zw			(5.8)
p
= E(z--H(l--Hp)z+(?--Hp)zW			(5.9)
= E(z).			(5.10)
This correction can be used to make the local pass much more efficient while
keeping the expected value at each pixel the same as if every interaction was re-
computed completely.
The most expensive computation for each link is the visibility. Shooting a
fixed number of rays per link results in more rays being used for points that have
more links. The more links connected to a point, the less each link contributes on
average. A better alternative is to shoot a fixed number of rays r per point. We
can distribute these rays based on the relative error of each link, using S?r visibil-
ity rays for link i. When S?r < 1, we shoot a single ray with probability p =
otherwise the value for the visibility stored in the interaction is used. The unbi-
ased approach described above is used to adjust the computed visibility.
The same approach can be used for recomputing the unoccluded form-factor
F . If a maximum of f form-factors on average are desired, then compute the an-
alytic form factor to x whenever fSj> 1. Again, if fSj < 1, compute the analytic
form factor to x with probability f??, otherwise use the estimate on the link. As
with the visibility, bias must be corrected for cases when fSj < 1.
94
In complex environments there can be several thousand links contributing
to the radiance at x. Many of these links carry very little energy and many of
them have very low error. Rather than computing thousands of form factors and
shooting thousands of rays per pixel we can compute tens or at worst hundreds
of form factors and visibility checks per pixel enabling vast computational sav-
ings without changing the expected solution.
5.4 Reconstruction with Texture Mapping
There are several interesting extensions that can be added using this form of re-
construction. The first of these is that transparency and specularity can be sim-
ulated. Specular and transparent surfaces often add much to the realism of an
environment, although there are still no good ways of capturing caustics in the
global pass. (An approach by Chen et al.[CRMT9l] does make an attempt.) Since
the ray tracing method handles specular transmission and reflection very well,
these procedures can be incorporated into the reconstruction. By propagating
these additional rays throughout the environment, and doing the pixel recon-
struction for any surface in the ray tree with a diffuse component, these effects
can be modeled. This approach is commonly used in ray tracing-based recon-
struction [WCG87].
Another advantage that is unique to the approach presented in this chapter as
well as in Reichert [Rei92] is the natural inclusion of many types of texture map-
ping. The radiance at the point can be found by using a BRDF from a texture
map, allowing images or surface properties to be texture mapped onto the sur-
face. When the irradiance is being computed, the normal at x can be displaced,
allowing bump mapping of the surfaces. These variations in the surface will
95
Figure 5.4: Local pass and corresponding global pass.
be shaded correctly because the shading at the point is being recomputed using
these new parameters.
Because these perturbations of the surface are not present in the global pass,
the global solution may have unknown errors. Many of these, such as minor
bump mapping, will generally be insignificant. Even image textures will not
cause many problems assuming the original color of the surface is close to the
average color of the texture map.
5.5 Results
The reconstruction approach presented here can produce images with very few
noticeable artifacts. This section will show global and local solutions to show the
effectiveness of this technique.
The first example shows the ability of the local pass to reconstruct accurate
96
Figure 5.5: Local passes with va?ing number ofsamples.
shadow boundaries. In Figure 5.4 we can see an image and the global solution
that was used to compute it. The shadows on the floor and on the table are accu-
rately captured, both the soft shadows from the ceiling light and the sharp shad-
ows from the sunlight.
The second example shows the degradation of image quality as the number
97
Figure 5.6: Unbiased (left) and biased (right) reconstruction.
of rays and form factors used for each pixel decreases. As can be seen in Figure
5.5, the noise in each image increases as the number of samples decreases. This
gives an acceptable way to continuously vary the quality of the final solution.
For example for this image, when fewer than ten visibility rays are cast for each
pixel, meshing artifacts are apparent. The number of visibility rays and form fac-
tors is responsible for reconstructing the incident irradiance over the entire hemi-
sphere, and ten samples over the hemisphere are not adequate. The expected
value is still correct because of the unbiased sampling.
In order to see why we need to correct for the bias1 two images were rendered
with and without bias correction. In Figure 5.6, the left image is computed with-
out correcting for bias, while the right image is computed with the bias correc-
tion from Equation (5.7). The mesh lines in the left image are very apparent to
the eye. This is because the errors on patches tend not to cancel each other out.
98
Figure 5.7: Local passes with texture mapping and transparent surfaces.
It is very common for most of the interactions to err in the same direction, either
brighter or darker. This causes the error in the biased approach to be noticeable.
In Figure 5.7 the surface properties were modified by texture maps and spec-
ular and transparent surfaces were added. The table top is now glass. The sur-
face of the glass has been altered with a bump map as can be seen in the re-
flected and refracted images, and the reflectance of the walls and floor have been
modified by procedural textures [Per85j. In addition1 there are several glass sur-
faces which contain reflections. This technique is very general and could handle
curved surfaces in addition to arbitrary BRDF's and texture mapping. However,
it is an expensive approach because of this generality.
Chapter 6
Conclusions and Future Work
6.1 Summary
The error bounding framework enables efficient global illumination algorithms
to be developed. Through the use of error bounds it is possible to determine
the effect of various approximations on the solution. Once we are able to deter-
mine the effects, it then becomes possible to use the coarsest of these approxima-
tions that remain within the desired accuracy. These coarse approximations are
much cheaper to compute, and therefore much unnecessary computation can be
avoided. Avoiding unnecessary computation reduces the solution time dramat-
ically.
The clustering approach speeds up the computation by enabling interactions
to occur between collections of surfaces rather than between the individual sur-
faces. This can be done much more cheaply than computing interactions be-
tween surfaces directly. The use of clustering reduces the quadratic complex-
ity of the initial phase of the hierarchical radiosity algorithm to a complexity of
O(s log s). The result is that hierarchical radiosity with clustering can now be
99
loo
used for environments that are more than an order of magnitude larger than the
environments that were feasible with hierarchical radiosity alone. Solution times
for sample scenes have been reduced by over two orders of magnitude.
The use of importance makes solutions of complex environments faster as
well. The advantage of solving for only a subset of the environment allows rapid
and accurate solutions of particular views in complex environments. For very
complex environments, solving for a subset of the environment allows solutions
that would not be possible to compute for the entire environment. The statistics
in this thesis show dramatic speedups over the fastest previous radiosity algo-
rithms.
Once a global solution has been computed, it is necessary to create an im-
age of the result. This thesis presents an efficient approach that produces im-
ages of very high visual accuracy. This reconstruction approach can also auto-
matically handle arbitrary reflectance functions, curved surfaces, and texturing.
The efficiency is obtained because most of the work is spent recomputing inter-
actions that contained the highest errot The estimates of these errors as well as
the global solution for the environment were computationally feasible because
importance restricted the amount of work needed in many parts of the environ-
ment. Clustering eliminated the initial overhead required for a hierarchical ra-
diosity solution as well as allowing very little work to be done in areas that could
not affect the view. These three advances in global illumination allow accurate,
photo-realistic images of environments that are more than an order of magnitude
larger than previous approaches could handle.
101
6.2 Limitations and Future Work
6.2.1 Arbitrary Reflectance Functions
Almost no materials are ideal diffuse. Although the assumption of diffuse sur-
faces is appropriate for many environments, in order to get realistic images for
many scenes it is necessary to account for the transfer of energy between non-
diffuse surfaces. There have been a several attempts to incorporate non-diffuse
surfaces into radiosity [1CG86,SAWG91,AH93]. The approaches by Sillion and
Auperle are both feasible, although both have some limitations. Either approach
will have difficulty with specular surfaces. Specular and near specular surfaces
are difficult because most of the energy coming in from a particular direction
goes out along a single reflected direction. This makes representing the distribu-
tion of energy over the hemisphere very costly. In addition, bounding the energy
transfer between specular surfaces will not result in tight bounds.
Importance becomes even more useful in non-diffuse environments. Impor-
tance can be used to determine the effect that light leaving in a particular direc-
tion will have on the initially important surfaces. This means that light leaving in
other directions can be ignored. As mentioned by Aupperle [AH93j, importance
is necessary in a non-diffuse environment because the number of interactions is
much higher.
Clustering with arbitrary reflectance functions should not be overly difficult,
assuming an approach used by Sillion [SAWG91] where reflectance distributions
are stored on the surfaces. Clustering requires storing both an incident distribu-
tion and an exitant distribution on a cluster. For non-diffuse environments the
same techniques that store distributions over patches can be used to store dis-
102
tributions for each cluster. The difficulty comes in bounding the errors in rep-
resenting arbitrary distributions, as well as in bounding the transfers between
non-diffuse surfaces. For surfaces that are close to diffuse, the problem is lim-
ited. For surfaces that are close to specular, the problem will be difficult.
6.2.2 Arbitrary Geometries
With very few exceptions [Zat93], radiosity approaches currently work only with
polygonal environments. This should not really be a problem when using error
bounding approaches to radiosity Transfers between curved surfaces can still
be estimated and bounded, soit should be relatively straightforward to compute
global solutions for curved surfaces. The main difficulties seem to be obtaining
good bounds on the transfers, and trying to accurately estimate the transfers so
as to reduce the number of artifacts that are present.
A more interesting problem with geometries comes when environments are
represented procedurally, or when they are so large that they can not be com-
pletely stored in memory. The clustering approaches presented in this thesis start
with the initial surfaces and build the hierarchy up from them. This is not a vi-
able approach with procedural models as it would require evaluating all of the
geometry at the beginning of the rendering process. It will be necessary to de-
velop new clustering approaches that do not require complete information about
the surfaces contained within the clusters before procedural models can be used
in radiosity Importance may prove a necessity in situations like this. It may be
impossible to solve for the radiances everywhere, simply because it is impossible
to create all of the geometry. In these cases, solving for some subset of the geom-
etry is necessary to minimize the amount of procedural geometry that needs to
103
be evaluated.
6.2.3 Input and Output
An often overlooked part of global illumination is the data used in the sim-
ulation. An accurate rendering is not useful unless the geometry surface re-
flectances, and light sources have been described accurately. It will be necessary
to have highly accurate descriptions of the reflectance characteristics of materi-
als. For lights it will be necessary to know the spectral and spatial distributions of
the energy they emit. Without this information, global illumination algorithms
are effectively solving the wrong problem [ATS94].
Another overlooked part of the global illumination process is the display of
the images we create. These images are stored with physical intensities at each
pixel. These intensities need to be mapped to the intensities producible on a dis-
play screen or on film or paper. The ranges of intensities that are displayable
by these output devices are much smaller than the range of intensities that ex-
ist in real world scenes. Because of this, some sort of compression of the dy-
namic range of the image must be done in order to display it. This compression
should preserve the perceptual characteristics of the image as much as possible.
Some work in this area has been done tTR91,CHS+93], however much more is
still needed.
Bibliography
[AH93]
Larry Aupperle and Pat Hanrahan. A hierarchical illumination al-
gorithm for surfaces with glossy reflection. In Computer Graphics
Proceedings, Annual Conference Series, ACM SIGGRAPH, pages
155--H162, 1993.
[AK87] James Arvo and David Kirk. Fast ray tracing by ray classification.
Computer Graphics, 21(4):5564, July 1987.
[AK9Oj James Arvo and David Kirk. Particle transport and image synthesis.
Computer Graphics, 24(4):6366, August 1990.
[App85] A. A. Appel. An efficient program for many-body simulation. SlAM
I. Sci. Stat. Computing, 6(1):8?103, 1985.
[ARBJ9Oj
John M. Airey, John H. Rohlf, and Frederick P. Brooks Jr. To-
wards image realism with interactive update rates in complex vir-
tual building environments. In Proceedings of the 1990 Symposium
on Interactive 3D Graphics (Snowbird, Utah, March 25-28, 1990), vol-
ume 24, pages 41--H50. ACM, March 1990.
[Arv86j James Arvo. Backward ray tracing. In Developments in Ray Tracing,
SIGGRAPH `86 Course Notes, volume 12, August 1986.
[ATS94]
James Arvo, Kenneth Torrance, and Brian Smits. A framework for
the analysis of error in global illumination algorithms. Tech Report,
January 1994.
Daniel R. Baum, Stephen Mann, Kevin P. Smith, and James M.
Winget. Making radiosity usable: Automatic preprocessing and
meshing techniques for the generation of accurate radiosity solu-
tions. In Proceedings of SIGGRAPH'91 (Las Vegas, Nevada, July 28--H
August 2, 1991), volume 25, pages 5160. ACM, July 1991.
[BMSW91]
104
105
[CCWG88j
[CF90j
[CGIB86]
[CHS+93]
[CRMT91j
Michael F. Cohen, Shenchang Eric Chen, John R. Wallace, and Don-
ald P. Greenberg. A progressive refinement approach to fast radios-
ity image generation. Computer Graphics, 22(4):75--H84, August 1988.
A. T. Campbell, III and Donald 5. Fussell. Adaptive mesh genera-
tion for global diffuse illumination. In Proceedings of SIGGRAPH'90
(Dallas, Texas, August 6--H10, 1990), volume 24, pages 155--H164. ACM,
August 1990.
Michael F. Cohen, Donald P. Greenberg, David 5. Immel, and
Philip J. Brock. An efficient radiosity approach for realistic im-
age synthesis. IFEE Computer Graphics and Applications, 6(2):2635,
March 1986.
K. Chiu, M. Herf, Peter Shirley, 5. Swamy, C. Wang, and K. Zimmer-
man. Spatially nonuniform scaling functions for high contrast `in-
ages. In Proceedings of Graphics Interface`92 (May 1993), pages 240--H
250, May 1993.
Shenchang Eric Chen, Holly E. Rushmeier, Gavin Miller, and Dou-
glass Turner. A progressive multi-pass method for global illumi-
nation. In Proceedings of SIGGRAPH'91 (Las Vegas, Nevada, July 28--H
August 2, 1991), volume 25, pages 165--H174. ACM, July 1991.
[Dav57] B. Davison. Neutron Transport Theory. Oxford University Press, New
York, 1957.
[DSG91]
[Ess89j
[GCS93j
[Gou71]
[Gre88j
[GS87]
Julie O'B. Dorsey, Fran?is X. Sillion, and Donald P. Greenberg. De-
sign and simulation of opera lighting and projection effects. Com-
puter Graphics, 25(4):41--H50, August 1991.
E. Esselink. About the order of appel's algorithm. Technical report,
University of Groningen, Department of Computer Science, 1989.
Steven Gortler, Michael Cohen, and Philipp Slusallek. Radiosity and
relaxation methods. Technical report, Princeton University, Depart-
ment of Computer Science, 1993.
H. Gouraud. Continuous shading of curved surfaces. IEEE Transac-
tions on Computers, C-20(6):623628, June 1971.
Leslie Greengard. The Rapid Evaluation of Potential Fields in Particle
Systems. MIT Press, Cambridge, Massachusetts, 1988.
Jeffrey Goldsmith and John Salmon. Automatic creation of object
hierarchies for ray tracing. lEFE Computer Graphics and Applications,
7(5):1?20, May 1987.
106
[GTGB84]
Cindy M. Goral, Kenneth E. Torrance, Donald P. Greenberg, and
Benneft Baftaile. Modeling the interaclion of light between diffuse
surfaces. Computer Graphics, 18(3)213--H222, July 1984.
[Hec9Oj Paul 5. Heckbert. Adaptive radiosity textures for bidirectional ray
tracing. Computer Graphics, 24(4):145--H154, July 1990.
[Hec91]
[HSA91]
[ICG86j
[Kaj86j
[Kok93]
Paul 5. Heckbert. Simulating Global Illumination UsingAdaptiveMesh-
ing. Ph.D. dissertation, Department of EECS, UC Berkeley, Califor-
nia, June 1991.
Pat Hanrahan, David Salzman, and Larry Aupperle. A rapid hier-
archical radiosity algorithm. Computer Graphics, 25(4):197--H206, July
1991.
David 5. Immel, Michael F. Cohen, and Donald P. Greenberg. A ra-
diosity method for non-diffuse environments. Computer Graphics,
20(4):133--H142, August 1986.
James T. Kajiya. The rendering equation. Computer Graphics,
20(4):143150, August 1986.
Aijan J.F. Kok. Grouping of patches in progressive radiosity. In
Proceedings of the Fourth Furographics Workshop on Rendering (Paris,
France, June 1416, 1993), pages 221--H231, June 1993.
[KW86j M. H. Kalos and Paula A. Whitlock. Monte Carlo Methods, Volume I:
Basics. John Wiley & Sons, New York, 1986.
[Lew65j
[LSG94]
[LTG93]
Jeffery Lewins. Importance, The Adjoint Function: The Physical Basis of
Variational and Perturbation Theory in Transport and Diffitsion Problems.
Pergamon Press, New York, 1965.
Dani Lischinski, Brian Smits, and Donald Greenberg. Bounds and
error estimates for radiosity. Computer Graphics, July 1994.
Dani Lischinski, Filippo Tampieri, and Donald P. Greenberg. Com-
bining hierarchical radiosity and discontinuity meshing. In Com-
puter Graphics Proceedings, Annual Conference Series, ACM SIG-
GRAPH, pages 199--H208, 1993.
Ken Perlin. An image synthesizet In Proceedings of SIGGRAPH'85
(San Francisco, California, July 22--H26, 1985), volume 19, pages 287--H
296. ACM, July 1985.
[Per85]
107
(Rei92]
[RPV93]
[SAS92]
[SAWG91]
[SH81j
tSP89]
[TH93]
[TR91]
[TS91j
[WCG87]
Mark C. Reichert. A two-pass radiosity method driven by lights and
viewer position. Masters dissertation, Program of Computer Graph-
ics, Cornell Umversity? Ithaca, New York, January 1992.
Holly E. Rushmeier, Charles Patterson, and Aravin-
dan Veerasamy. Geometric simplification for indirect illumination
calculations. Graphics Interface `93, pages 227--H236, May 1993.
Brian Smits, James Arvo, and David Salesin. An importance-driven
radiosity algorithm. Computer Graphics, 26(4):273--H282, July 1992.
Francois Sillion, James Arvo, Stephen Westin, and Donald P. Green-
berg. A global illumination solution for general reflectance distribu-
tions. Computer Graphics, 25(4):187--H196, July 1991.
Robert Siegel and John R. Howell. Thermal Radiation Heat Transfer.
Hemisphere Publishing Corp., New York, second edition, 1981.
Francois Sillion and Claude Puech. A general two-pass method in-
tegrating specular and diffuse reflection. Computer Graphics, 23(4),
August 1989.
Seth Teller and Pat Hanrahan. Global visibility for illumination
computations. In Computer Graphics Proceedings, Annual Confer-
ence Series, ACM SIGGRAPH, pages 239--H246, 1993.
Jack Tumblin and Holly E. Rushmeier. Tone reproduction for real-
istic computer generated images. Technical Report GIT-GVU-91-13,
Georgia Institute of Technology, Atlanta, Georgia, 1991.
Seth J. Teller and Carlo H. Sequin. Visibility preprocessing for in-
teractive walkthroughs. In Proceedings of SIGGRAPH'91 (Las Vegas,
Nevada, July 28--HAugust 2, 1991), volume 25, pages 61--H70. ACM, July
1991.
John R. Wallace, Michael F. Cohen, and Donald P. Greenberg. A two-
pass solution to the rendering equation: A synthesis of ray tracmg
and radiosity methods. In Proceedings of SIGGRAPH'87 (Anaheim,
Caltfi?rnia, July 27--H31, 1987), volume 21, pages 311--H320. ACM, July
1987.
John Wallace, Kells Elmquist, and Eric Haines. A ray tracing algo-
rithm for progressive radiosity. Computer Graphics, 23(3):315--H324,
July 1989.
[WEH89]
108
[WRC88]
Gregory J. Ward, Francis M. Rubinstein, and Robert D. Clear. A
ray tracing solution for diffuse interreflection. In Proceedings ofSIG-
GRAPH'88 (Atlanta, Georgia, August 1--H5, 1988), volume 22, pages 85--H
92. ACM, August 1988.
Harold Zatz. Galerkin radiosity: A higher order solution method
for global illumination. In Computer Graphics Proceedings, Annual
Conference Series, ACM SIGGRAPH, pages 213--H220, 1993.
[Zat93]
