BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR93-1346
ENTRY:: 1993-10-14
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Discontinuity Meshing for Radiosity Image Synthesis
AUTHOR:: Tampieri, Filippo
DATE:: May 1993
PAGES:: 185
COPYRIGHT:: Filippo Tampieri - All RIghts Reserved
ABSTRACT::
The simulation of global illumination is one of the most fundamental problems 
in computer graphics, with applications is a wide variety of areas. This 
problem studies the light energy transfer between reflective surfaces in an 
environment. Initially derived from the field of thermal engineering, 
radiosity has emerged over the past several years as one of the most promising 
solution methods.

Despite having produced some of the most realistic-looking computer generated 
images to date, radiosity methods have not yet met with widespread 
acceptance. The main obstacle has been their need for very careful and time 
consuming user intervention, without which, current techniques are prone to 
generating a wide range of annoying visual artifacts. These artifacts are 
generally due to poor surface meshing, resulting in insufficient sampling 
density and ineffective sample placement.

This thesis investigates the roots of this problem by taking a step back from 
the traditional finite element formulation of radiosity and examining the more
general integral equation formulation. An analysis of the radiance functions 
described by this equation shows how umbra and penumbra boundaries as well as 
other sharp changes in illumination actually correspond to discontinuities in 
the radiance function and its derivatives. The results of this analysis have 
led to the concept of discontinuity meshing, whereby accurate approximations 
to the radiance functions are computed by explicitly representing their 
discontinuities as boundaries in the mesh.

This concept has been applied to the design of a discontinuity meshing 
algorithm for polyhedral environments. The algorithm is embedded in a 
progressive refinement radiosity system and uses piecewise quadratic 
interpolation to reconstruct a smooth radiance function while preserving 
discontinuities where appropriate.

The radiosity solutions produced by the new algorithm are compared against a 
photograph of a physical environment, an analytical solution, and a 
conventional, yet state-of-the-art, radiosity system, and its performance on 
architectural models of medium complexity is measured. The results are 
remarkably accurate both numerically and visually. The new discontinuity 
meshing algorithm drastically reduces, and in many cases eliminates, many of 
the annoying artifacts typical of conventional radiosity meshes, producing 
images of previously unattained quality. Moreover, the meshing is completely 
automatic and produces solutions that are highly view-independent.
END:: CORNELLCS//TR93-1346
BODY::
Discontinuity Meshing for
Radiosity Image Synthesis
Filippo Tampieri
Ph.D Thesis
93-1346
May1993
Department of Computer Science
Cornell University
Ithaca NY 14853-7501
DISCONTINUITY MESHING
FOR
RADIOSITY IMAGE SYNTHESIS
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
Filippo Tampieri
May 1993
Qc Filippo Tampieri 1993
ALLRIGHTSRESERVED
DISCONTINUITY MESHING
FOR
RADIOSITY IMAGE SYNTHESIS
Filippo Tampieri, Ph.D.
Cornell University 1993
The simulation of global illumination is one of the most fundamental prob-
lems in computer graphics, with applications in a wide variety of areas. This
problem studies the light energy transfer between reflective surfaces in an envi-
ronment. Initially derived from the field of thermal engineering, radiosity has
emerged over the past several years as one of the most promising solution meth-
ods.
Despite having produced some of the most realistic-looking computer gener-
ated images to date, radiosity methods have not yet met with widespread accep-
tance. The main obstacle has been their need for very careful and time consum-
ing user intervention, without which, current techniques are prone to generating
a wide range of annoying visual artifacts. These artifacts are generally due to
poor surface meshing, resulting in insufficient sampling density and ineffective
sample placement.
This thesis investigates the roots of this problem by taking a step back from
the traditional finite element formulation of radiosity and examining the more
general integral equation formulation. An analysis of the radiance functions de-
scribed by this equation shows how umbra and penumbra boundaries as well
as other sharp changes in illumination actually correspond to discontinuities in
the radiance function and its derivatives. The results of this analysis have led
to the concept of discontinuity meshing, whereby accurate approximations to the
radiance functions are computed by explicitly representing their discontinuities
as boundaries in the mesh.
This concept has been applied to the design of a discontinuity meshing al-
gorithm for polyhedral environments. The algorithm is embedded in a progres-
sive refinement radiosity system and uses piecewise quadratic interpolation to
reconstruct a smooth radiance function while preserving discontinuities where
appropriate.
The radiosity solutions produced by the new algorithm are compared against
a photograph of a physical environment, an analytical solution, and a conven-
tional, yet state-of-the-art, radiosity system, and its performance on architectural
models of medium complexity is measured. The results are remarkably accurate
both numerically and visually. The new discontinuity meshing algorithm drasti-
cally reduces, and in many cases eliminates, many of the annoying artifacts typi-
cal of conventional radiosity meshes, producing images of previously unattained
quality Moreover, the meshing is completely automatic and produces solutions
that are highly view-independent.
Biographical Sketch
The author was born in Novara, Italy, on December 1st, 1962. He was first ex-
posed to Computer Graphics while studying Computer Science at the Urnversita'
degli Studi di Milano, where he spent his senior year implementing the obliga-
tory ray tracer, the distinctive mark of any respectable computer graphnic.
During that time1 he grew very fond of his creature and, worried that the
world of industry might require him to abandon it, he decided to embrace the
profession of eternal student. In 1985, he entered a competition for a Fulbright
scholarship which eventually lead him to enter the Program of Computer Graph-
ics at Cornell University in August 1986.
While at Cornell, the author did research on global illumination algorithms
for parallel computer architectures. During this investigation, he was exposed
to many new and wild ideas and, despite his resistance, by the time he received
his Master's, he was showing the first symptoms of the radiosity virus.
In an attempt to fight his illness, in September 1988, the author returned to
Italy, where he worked as a computer graphics consultant on geometric model-
ing and more ray tracing.
The radiosity virus, though, proved too resilient to overcome and, by Sum-
mer 1989, the author was back at Cornell as a PhD candidate in search for a cure.
iii
There, he spent his last years researching an antidote that could save future fel-
low graphnics from spending their life making pictures of the dreadful "Cornell
box."
iv
Dedicato a Molly,
ai miei genitori
e a mio fratello Andrea
v
Acknowledgements
First of all, I would like to thank my parents for all their love, encouragement,
and support; they never faltered despite the many years, the long distance, and
my many "surprise announcements."
I am deeply grateful to Professor Don Greenberg. His enthusiasm and energy
have always been contagious; he has been a caring friend and his continuous
trust and support have helped through many hard times. As an advisor, he has
provided many suggestions that have greatly improved the quality of this the-
sis. As the director of the Program of Computer Graphics, he has provided the
terrific environment in which my work, and that of many other fellow students,
was carried out.
My appreciation goes to Professor Roy Hall for his help reviewing this thesis
and for the many discussions, personal and professional, that we have had over
the past few years.
Thanks also go to Professor Charles Van Loan for serving on my committee.
I would also like to acknowledge and thank Professor Daniele Marini for in-
troducing me to computer graphics, many years ago. Had he not gone through
the trouble of getting me onto a graphics workstation, I think my past few years
would have been very different.
vi
Next, I would like to thank my research partner, officemate, and friend, Dani
Lischinski, without whom this work would not have been possible. We thought
we were going to take a month and extend radiosity to curved surfaces, and
ended up two years later with beautiful images of the flattest models in the lab.
Working together has been instructive, productive, and especially fun. Most of
the ideas I present in this thesis are actually our ideas, born and matured out of
our endless discussions.
My gratitude also goes to all the past and present students and staff at the lab
who contributed in many ways to make my time here a truly enriching experi-
ence.
Many thanks to Ellen French and Fran Brown for protecting me from bu-
reaucracy and feeding me on many occasions. To Jonathan Corson-Rikert for
keeping the lab running so well in Don's absence. To Hurf Sheldon for keeping
the system runnmg so smoothly and for his patience and prompt help the many
times I needed assistance. To Xiao Dong He, sailor Sheldon, my back-to-back
officemate, for his easy-going attitude. To Ben Trumbore, DeathBed expert, for
being always ready to help. To Brian Smits for shedding light on the mysteries
of color and spectral curves. To Chris Schoeneman and Rick Pasetto for writing
graphics libraries that people want to use; to Rick for the colored padding noise
and to Chris for Widget Lisp. To Dan Kartch for the many wonderful cakes that
brought me back to the lab so many Mondays. To Harold Zatz for volunteering
to be our T??( guru. To Patrick Heynen for helping me find my way on the NeXT
machines.
I would also like to thank Greg Spencer for modeling the simple scene used
vii
to demonstrate the visual artifacts due to uniform meshing shown in Figures 3.1
and 3.2; Ted Hirnian for helping Dani and I set up the simple physical scene
shown in Figure 6.1; Matt Hyatt for modeling the Victorian Room, the White
Room, the Green Room, and the Pink Room used in Chapter 6; and Charles
D7Autremont for modeling the Mackintosh Room used in the same chapter.
Cornell has been a wonderful place to learn. It has been great for my aca-
demic education, but has also been the exploratory playground for many other
interests. I would like to thank sifu Maurice Haltom and all the instructors and
students at the Agape School of Movement Studies for teaching me Chinese mar-
tial arts and trying to teach me African Dance. Thanks to David Pearson for in-
troducing me to the rock climbing experience; to Karl Bohringer for the many
hours spent clinging to "the wall," and to Dan Tillemans and the other instruc-
tors for taking me to climb real rock. Thanks to Jen Whiting, Iris, and Ken for
pushing me through the "birth canal" in J4. Thanks to Sara for teaching the "for-
eigners office" how to sail catamarans. And thanks to Maya for squeezing my
big nose.
I would also like to take the opportunity to thank the Italian friends who kept
in touch during these many years away from home. Many thanks and lots of love
to Andrea, my brother, for growing closer with distance and time, for writing
the funniest letters I have ever read, and for keeping at it despite my infrequent
replies. Thanks to Paolo Bucci and Roberto Mameli, the other two fish out of
water, for sharing our moods and experiences. Thanks to the "cucadores" for
sharing their Saturday night whenever I was in town. And thanks to Stefano and
Laura Piccardi, Luigi Pergolini, Manuela Bandini, Lucilla Consalez, and Enrico
viii
and Gabriella Molinari for never forgetting their friend across the Ocean.
Lastly, I would like to thank my friend, Molly Pedersen1 for her love, friend-
ship, understanding, and support. She has been a wonderful companion and an
inexhaustible source of energy, helping me grow, taking me away from the office
on the rare sunny days, and showing me how to tell a dog's smile.
This work was supported by the National Science Foundation grant "Inter-
active Input and Display Techniques" (CCR-8617880) and by the NSF/DARPA
Science and Technology Center for Computer Graphics and Scientific Visualiza-
tion (ASC-8920219). I would also like to acknowledge the generous equipment
grant from Hewlett Packard Corporation on whose workstations this research
was conducted.
ix
Table of Contents
1 Introduction
2
3
Radiance Functions
2.1 Definition of Terms and Symbols
2.2 The Rendering Equation . .
2.3 The Radiosity Equation
2.4 Characteristics of Radiance Functions
2.4.1 D0 Discontinuities .
2.4.2 D' and D2 Discontinuities .
2.4.3 Propagation of Discontinuities
2.4.4 Singularities .
Previous Work
3.1			Radiosity. .
3.1.1 Diffuse Radiosity
3.1.2 Full Matrix Solution
3.1.3 Progressive Refinement Radiosity
3.2			Meshing
3.2.1 Uniform Meshing
3.2.2 Adaptive Subdivision .
3.2.3 Hierarchical Meshing
3.2.4 Discontinuity Meshing
3.3 Direct Energy Transfer
3.3.1 Numerical Form Factors
3.3.2 Analytical Form Factors
3.3.3 Non-Constant Exitance .
3.4 Radiance Function Reconstruction
3.4.1			Constant Elements			. .
3.4.2 Gouraud Shading
3.4.3 Bilinear Interpolation . .
3.4.4 Higher Order Interpolation
5
5
8
10
12
15
16
21
23
25
25
25
26
27
28
29
34
35
37
44
45
49
51
52
53
53
54
55
x
4
5
6
7
Single Source Discontinuity Meshing
4.1 Single Source Discontinuity Meshing--HAlgorithm Overview
4.2			Locating the Discontinuities
4.2.1			Locating D0 Discontinuities .			. .			.			. .
4.2.2			Locating D1 and D2 Discontinuities
4.3			Constructing the Discontinuity Mesh			.			.
4.3.1			Implementation.
4.4			Sampling and Reconstruction on the			Receiver			.
4.5			Computing the Source Contribution to a			Point
4.5.1			Implementation . . .			. .			.
4.6			Adaptive Subdivision . . . . . .			.			. .			.
4.6.1			When to Subdivide
4.6.2			How to Subdivide . . . . .			. .			.
Discontinuity Meshing Radiosity
5.1 Discontinuity Meshing Radiosity--HAlgorithm Overview
5.2			Locating the Discontinuities . 			. 			.			.
5.2.1			Implementation . . 			.			.
5.3			Computing the Source Contribution to a Point
5.4			Combining Radiance Functions . . 			. . 			. .			. .
5.5			Efficiency Considerations. 			. .
5.5.1			Computation of Radiance Discontinuities
5.5.2 Computation of the Contribution of a Source to a Point
Results
6.1			Physical vs. Simulated World .
6.2			Single Source Discontinuity Meshing vs. Analytical . . .			.
6.3			Accumulating Radiance Functions by Resampling vs. Exact			Ad-
dition
6.4			Discontinuity Meshing Radiosity vs. Conventional Radiosity
6.5			Statistics. . . . . . . . . . . . . . . .			.
Summary and Conclusion
7.1			Summary and Main Contributions			.
7.2			Further Research
7.2.1			Robustness . . . .			. .			.
7.2.2			Data Structures
7.2.3			Performance and			Flexibility. .			.
7.2.4			Solution Accuracy and Perceptual			Issues
7.2.5			Extensions . . . .			. .			. .
Conclusion
7.3
Bibliography
xi
56
57
60
60
64
79
82
87
92
94
96
97
100
104
104
106
107
109
110
112
115
117
119
119
121
125
139
145
154
154
158
158
159
160
162
162
165
166
ListofTables
6.1
6.2
6.3
6.4
6.5
Discontinuity meshing vs. analytical solution126
Statistics for the discontinuity meshing solution of the Victorian
Room model144
Timings for the discontinuity meshing solution of the Victorian
Room model. . . . 144
Statistics for the discontinuity meshing solutions of the room modelsi52
Timings for the discontinuity meshing solutions of the room models152
xii
List of Figures
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
Geometry for BRDF's. .			.
Geometry for Equations (2.1 ) and (2.2)			.
Photograph of a simple scene			.
Discontinuity lines in a simple scene 			.
An example of a D0 discontinuity . .
A D1 discontinuity caused by an edge-vertex (EV) event
A D2 discontinuity caused by an edge-vertex (EV) event
An edge-edge-edge (EEE) event
Critical curves on receiver R for four simple scenes.
Example of function singularity . . . . .
Visual artifacts due to uniform meshing
More visual artifacts
Hierarchical meshing on a pair of perpendicular polygons
Computing D0 discontinuities prior to meshing .
Shadow volumes and umbra and penumbra volumes
The hemicube algorithm
Ray-traced form factors . .			. . .			. .
Analytical form factors
Pseudocode for the single source discontinuity meshing algorithm
Pseudocode for the computation of D0 discontinuities (Part 1)
Pseudocode for the computation of D0 discontinuities (Part 2)
A wedge corresponding to a VE event and its interaction with the
surfaces in the environment
A wedge corresponding to an EV event and its interaction with
the surfaces in the environment
Pseudocode for the generation and classification of visual events
Inconsequential visual events
Pseudocode for the location of the discontinuities generated by a
VE event (Part 1)
Pseudocode for the location of the discontinuities generated by a
VE event (Part 2) . .			. .			. .			. . .			.
6
9
14
14
17
18
19
20
22
23
30
31
36
39
41
46
48
50
58
61
63
67
69
70
72
o+ 74
75
4.9
xiii
4.10 Pseudocode for the location of the discontinuities generated by a
VE event (Part 3). . . 76
4.11 Organization of coplanar surfaces in a BSP-tree node 77
4.12			Incremental construction of a DM-tree. . . 81
4.13 Pseudocode for the DM-tree data structures 83
4.14 Pseudocode for the insertion of a critical curve into the disconti-
nuity mesh of a surface84
4.15 Pseudocode for the insertion of a discontinuity segments into the
leaves of a DM-tree86
4.16 Splitting a face to capture a discontinuity segment . . . . . . . 86
4.17 Triangulating a mesh element87
4.18 Radiance reconstruction by piecewise polynomial interpolation . 89
4.19 Quadratic interpolation over a triangular mesh element 90
4.20			Treatment of D0 vertices. . . 90
4.21 Pseudocode for the evaluation of the radiance function at a point x
on a surface. . . . . . . . . 92
4.22 Analytical form factors . . . . . . . . . 93
4.23 Pseudocode for the computation of form factors in the presence
of occlusions. . . 95
Rivara's 2-triangles mesh refinement algorithm .102
106
4.24
5.1
5.2
5.3
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
6.9
6.10
6.11
6.12
6.13
Pseudocode for the discontinuity meshing radiosity algorithm
Pseudocode for the generation and classification of the visual events
associated with the discontinuity curves on secondary light sources108
Combining the illumination from two light sources in the pres-
ence of a simple obstacle . .			. .			113
Photograph of a simple scene. . .			. .			. .			.			120
Computer simulation of a simple scene			120
Plot of the radiance values across a scanline . .			. . .			. .			. .			122
Discontinuity curves and triangular mesh for the simple			scene			126
Full view. Discontinuity meshing vs. analytical			solution			. .			.			127
Floor only. Discontinuity meshing vs. analytical			solution			128
Plot of the radiance function along room scanline			. . .			. .			.			129
Plot of the radiance function along floor scanline			130
Accumulating radiance functions by resampling vs. exact addition 132
Error introduced by the radiance accumulation algorithm . . . 			134
Combined effect of progressive refinement and radiance accumu-
lation algorithms			136
Discontinuity meshes and individual contributions for a discon-
tinuity meshing radiosity solution .			.			137
Growth in mesh size as a function of progressive refinement iter-
ations			. .			. .			.			138
xiv
6.14			A room with a view . .			. .			140
6.15			A closer look at the wall			. 			. . . . . .			. 141
6.16			A closer look at the chair . .. .			. 142
6.17			White Room. (a) A ray-traced view of a			discontinuity meshing
radiosity solution. (b) The corresponding discontinuity mesh. .			. 146
6.18			Green Room. A ray-traced view of a discontinuity			meshing ra-
diosity solution. .			. 147
6.19			Green Room. The final discontinuity mesh.			148
6.20			Pink Room. (a) A ray-traced view of a discontinuity			meshing ra-
diosity solution. (b) The corresponding discontinuity			mesh149
6.21			Mackintosh Room. (a) A ray-traced view of a			discontinuity mesh-
ing radiosity solution. (b) The corresponding			discontinuity mesh.			150
xv
Chapter 1
Introduction
Realistic image synthesis is one of the oldest and most fundamental problems
of computer graphics. The initial attempts at creating convincing renderings of
reality were based on somewhat empirical techniques [GOUR71,PHON75j, but
their impact was so profound that those algorithms are now an essential compo-
nent of any modern graphics workstation. As the field matured, however, many
researchers have recognized the need for a more rigorous approach, whereby
image synthesis is not just the art of making "pretty pictures," but is instead a
faithful simulation of physical laws.
Recently, much research in computer graphics has focused on global illumi-
nation, the complex phenomena describing the interaction of light energy be-
tween reflective surfaces in an environment. Over the past decade, two different
methods have emerged: Radiosity and Monte Carlo methods.
Initially derived from the field of thermal engineering [SPAR63j, radiosity re-
duces the global illumination problem to a system of linear equations employing
techniques similar to those used in finite element analysis [GORA84,COHE85,
NISH85a].
2
Monte Carlo methods, instead, use the integral equation formulation of the
global illumination problem, and try to solve it by point sampling the integrand
using ray tracing techniques [KAJ186,WARD88j.
Both methods have generated some of the most realistic images to date and
both are very computationally intensive. However, there is a fundamental differ-
ence between the two: radiosity approaches compute solutions in object space,
while Monte Carlo methods yield solutions in image space. As a consequence,
the same radiosity solution can be used to render images from many different
viewer's positions, while Monte Carlo methods must recompute the entire sim-
ulation every time the viewer's position is changed. We classify algorithms ac-
cording to these properties into view-independent and view-dependent methods re-
spectively.
The ability to generate interactive "walk-throughs" from view-independent,
physically accurate, global illumination simulations, makes radiosity methods
attractive for a wide variety of applications, particularly enclosed spaces such as
in architectural design and lighting analyses.
Despite all the areas that could benefit from radiosity, this method has not
yet met with widespread acceptance. This is due for the most part to the lack
of satisfactory meshing algorithms. The meshing techniques used in current ra-
diosity systems are prone to generating a wide range of annoying visual artifacts
that can only be eliminated at the cost of tedious and time consuming manual
adjustments to the mesh, often made surface by surface and repeated on a trial-
and-error basis until a proper mesh density is found.
Clearly, new meshing techniques are needed to generate high quality solu-
3
tions with the degree of reliability and automation necessary for practical appii-
cations. This thesis investigates a new approach to meshing, discontinuity mesh-
ing, that is completely automatic and drastically reduces the visual artifacts that
afflict conventional radiosity algorithms. This approach promises to play a sig-
nificant role in bringing radiosity methods closer to the demands of useful ap-
plications.
An analysis of the characteristics of radiance functions provides the insight
and rationale behind the concept of discontinuity meshing. The analysis shows
how umbra and penumbra boundaries, as well as other sharp changes in illu-
mination actually correspond to discontinuities in the radiance function and its
derivatives. The shape, location, and order of these discontinuities is determined
by the interaction of the light sources and obstacles in the environment. Dis-
continuity meshing explicitly represents radiance discontinuities as boundaries
in the mesh. This is in contrast to standard techniques which generate mesh
elements by geometric subdivision alone. When one of these discontinuities
crosses a mesh element, a large error may occur because the simple interpola-
tion schemes used to reconstruct the radiance function cannot reproduce sharp
illumination changes within an element.
This idea has been applied to the design of a discontinuity meshing algo-
rithm for polyhedral environments. The algorithm is embedded in a progres-
sive refinement radiosity system for ideal diffuse reflectance distributions and
uses piecewise quadratic interpolation to reconstruct a smooth radiance func-
tion while preserving discontinuities where appropriate.
The actual implementation has generated some remarkably accurate solu-
4
tions for environments of moderate complexity showing significant improve-
ments in both the numerical and the visual quality of the simulations. Although
the system is not yet suitable for very complex models, it does prove the validity
of the concepts and ideas, and our results suggest that discontinuity meshing
may become an important component of future radiosity systems.
A complete description of the algorithm and experimental results and com-
parisons is presented in the following sections.
Chapter 2 introduces the necessary terminology, discusses the global illumi-
nation problem with a formulation that is cast in terms of radiance functions,
and investigates the characteristics of the latter ones.
Chapter 3 reviews previous work, discussing favorable properties and limi-
tations of the meshing, direct energy transfer, and radiance reconstruction algo-
rithms used in current radiosity systems.
Chapter 4 discusses the design and implementation of a discontinuity mesh-
ing algorithm capable of capturing the illumination cast by a single area light
source of constant emission onto the surfaces of a polyhedral environment.
Chapter 5 presents a modified progressive refinement radiosity system that
uses the new discontinuity meshing algorithm, extending it to correctly compute
light transfers from secondary light sources.
Chapter 6 reports on a series of experiments designed to test the performance
of the new system and assess the validity of the ideas presented in the thesis.
Finally, Chapter 7 discusses advantages and limitations of the proposed meth-
ods and provides guidelines for future research.
Chapter 2
Radiance Functions
This chapter discusses the global illumination problem with a formulation that is
cast in terms of radiance functions. According to this formulation, each surface
in a scene has an associated radiance function that describes the illumination at
any poffit on the surface.
The characteristics of radiance functions are then investigated. This study
provides the tools for the analysis of past and current radiosity algorithms pre-
sented in the next chapter as well as the insight into the problem that led to the
new solution algorithms presented later in this thesis.
2.1 Definition of Terms and Symbols
This section introduces the symbols and definitions that will be used through-
out the rest of this thesis. Most of the units, symbols, and defining equation of
the fundamental radiometric quantities are taken from the lES Lighting Hand-
book [1ES91] as given by the Illuminating Engineering Society.
5
6
Th
4)
dA?
Y
Figure 2.1: Geometry for BRDF's
Bidirectional reflectance-distribution function (BRDF): the ratio of the radiance
reflected in direction (Or, ?r) to the radiant flux density incident from direc-
tion (O?, ?j) through a unit solid angle (see Figure 2.1). The BRDF describes
the reflecting properties of a surface element placed at x and will be de-
noted by:
fr(X', x, x").
The latter notafion will be used to indicate radiant energy arriving at point
x from the direction of x' and leaving x in the direction of x"
orby
fr(x;Oi, (ki;Or, ?r)
z
7
Irradiance (E): the density of radiant flux incident on a surface (d?/dA). It is
denoted by E and measured in watts per square meter (Wm-2).
Radiance (L): the radiant flux per unit solid angle per unit projected area. More
precisely1 the radiance leaving a differential surface area in direction (0, ?)
(relative to the surface normal) is equal to d2?/(dwdA cos 0). It is denoted
by L and is measured in watts per steradian per square meter (Wsr-1m-2).
Radiant energy (Q): energy traveling in the form of electromagnetic waves. It
is denoted by Q and is measured in joules U)
Radiant exitance (M): the density of radiant flux leaving a surface (d?/dA). It
is denoted by M and measured in watts per square meter (Wm-2).
Radiant flux (?): the radiant energy flowing through an area per unit time (dQ/dt).
It is denoted by ? and is measured in watt (W).
Radiant flux density (d?/dA): the radiant energy flowing through a unit area
per unit time (d?/dA). It is referred to as radiant exitance when the flow
is leaving the surface and as irradiance when the flow is incident on the
surface.
Radiant intensity (I): the radiant flux emanating from a point light source through
a unit solid angle around a given direction (d?/dw). It is denoted by I and
measured in watts per steradian (Wsr-1).
Radiosity (B): the density of radiant flux leaving a Lambertian reflector. This is
the name commonly used for radiant flux density in the computer graph-
8
ics literature. It is denoted by B and measured in watts per square meter
(Wm-2).
Reflectivity (p): the fraction of incident radiant energy that is reflected (as op-
posed to absorbed or transmitted.) It is denoted by p(x). This quantity is
commonly used to specify the reflecting properties of a Lambertian surface
and is related to a constant BRDF by:
p(x) =
cos Ordw
cos Ordw)fr(x) = 7rfr(X)
All of the above quanfities are actually wavelength dependent. Wherever needed,
this dependency will be made explicit by specifying the wavelength A.
2.2 The Rendering Equation
The transfer of radiant energy is governed by the laws of physics. The complex-
ities of the interaction of light and surfaces in an environment were abstracted
by Kajiya (KAJI86] in 1986 in a compact formulation of the global illumination
problem, generally known as the rendering equation. This is a Fredholm integral
equation of the second kind [HFCK9laJ and is presented here in a slightly mod-
ified form [SHIR9Ob]:
L(x, x", A) = Le(X, x", A) + ffr(M? x? x", A)L(x', x, A) cos$i2cos9i v(x, x')dA(x') (2.1)
where
L(x? x", A) is the radiance at point x in the direction of point x" at wavelength A.
9
6.
4X
Figure 2.2: Geometry for Equations (2.1) and (2.2).
Le(X, x", A) is the radiant exitance per unit solid angle emitted from x towards x"
at wavelength A.
S is the union of all the surfaces in the environment.
v(x, x') is a visibility term; it is 1 if x and x' are visible to each other, and 0 other-
wise;
dA(x') is a differential surface area centered at x'.
r is the distance between x and x';
Ot and O? are the angles between the surface normals at x and x' and the line
connecting the two points (see Figure 2.2).
The rendering equation as given above is well suited for view-dependent
global illumination algorithms like ray tracing and Monte Carlo techniques. In
fact, these algorithms typically discretize the integration by tracing rays from x
in a number of random directions, chosen perhaps according to the distribution
10
offr. Instead, view-independent algorithms, such as radiosity, typically work by
computing the energy transfer between pairs of surfaces (or surface elements)
in the environment. For this latter class of algorithms, it is therefore more con-
venient to break the domain of integration & into the surfaces s1 ..... . ,S? of the
scene, to get:
Lt'(X, x", A) =
L?j(X, X11) A) + ? $ f?i(M? x, x", A)Lj(X', ?, ? ?cos O? cos
S1ES `Es1			r2
"I			v(x.1 1aA(x')
(2.2)
According to this formulation, each surface Sin the environment has an as-
sociated radiancefrnction Lj. Solving for the set of Lj'5, i ? [1. .n], yields a view-
independent solution to the global illumination problem.
If all the surfaces are ideal diffuse (Lambertian) reflectors, the general render-
ing equation can be simplified to:
L,(x, A) = L?j(X, A) +fri(x, A) ?			A)C0S O? cos Oj
L1(x',			r2			v(x, x')dA(x')			(2.3)
S)ES x'Es1
This equation will be referred to as the diftse rendering equation and most of the
remainder of this thesis will be dedicated to its solution.
2.3 The Radiosity Equation
The computer graphics literature has traditionally formulated the diffuse global
illumination problem in terms of radiosity [GORA84]. In this section we will
show how the tradifional radiosity equation can be derived from the diffuse ren-
dering equation given above. In order to simplify the equations, the dependency
on the wavelength A will be omitted.
11
By definition, the radiosity Bi(X), is the radiant flux density leaving surface Sj
at point x. This quantity can be computed from the radiance, Li(X) by integrating
the latter over a hemisphere centered above x. For an ideal diffuse (Lambertian)
reflector, the radiance is equal in all direction, thus:
Bt(x) =
Li(X) cos Odw
cos Odw)Li(x) = ?Li(x).
The traditional radiosity formulation requires that the environment be dis-
cretized into patches of constant radiosity and constant reflectivity. For simplic-
ity we will assume that these patches correspond to our surfaces Sj'5. Then, the
area averaged radiosity Bi and reflectivity Pi for a finite area surface Sj of area Aj
are given by
and
Bj = 1			7rLi(x)dA(x)
Aj XESi
1 7rfri(x)dA(x).
Pi = AT xEst
We can now derive the traditional radiosity equation by integrating each side
of Equation (2.3) over the area of surface s?:
fESi 7rLj(X)dA(X)
--H fES ?LCj(X)dA(X) +
? Lj(X') cos O? cos
fES 7rfri(X) ?			r2			v(x, x')dA(x')dA(x)
S1ES x'ES
--H fESi ?LCj(x)dA(x) +
?iZf L1(x') cosO;osQ v(x? x')dA(x')dA(x)
S1ES xES1 x'ES1
BiAi =
12
BetAt + pt ? Bj f $ cosOt cos 01v(x, x')dA(x')dA(x)
S1ES XESi x'Es1 irr2
BetAt + pi ? BiAtFit
S1ES
(2.4)
where Bet is the emitted radiosity and Ftj is the form factor from surface ?t to
surface St. Dividing each side of the equation by At yields the traditional radiosity
equation:
Bj =Bet+pt?BtFtt			(2.5)
S1ES
2.4 Characteristics of Radiance Functions
In trying to evaluate existing global illumination algorithms or to design effec-
tive new ones it is useful to examine actual radiance functions. Determining
which characteristics of radiance functions are mostly responsible for the "real-
istic look" of actual objects and surfaces can help concentrate research efforts on
subproblems whose solution will yield the highest payoffs. Furthermore, under-
standing the mechanisms responsible for the most crucial properties of radiance
funcfions can provide important clues to the design of effective simulation tech-
niques.
This section is concerned with the characteristics of the radiance functions
related to the diffuse rendering equation (see Equation (2.3)). In the real world,
these functions exhibit an abundance of subtle and complex illumination details
that are hard to simulate, even for the simplest environments. Figure 2.3 shows
a photograph of a cardboard triangle suspended over a matte white floor illu-
minated by a square diffuse light source. The most obvious characteristic of the
radiance function across the floor is the contrast between the smoothness of the
13
unocciuded areas and the sharpness of the boundaries of the triangular shadow.
The shadow itself is divided into a darker region (umbra) that is completely oc-
cluded from the light source and a surrounding area (penumbra) that provides
a transition to the lighter regions outside the shadow. Umbras and penumbras
are marked by boundaries of different sharpness; for example, the variation in
radiance at the lower and outer boundary of the penumbra is so sharp that it
gives rise to a bright Mach band, while that towards the upper left of the picture
is much less noticeable.
An analysis of radiance functions will show that these illumination details
correspond to discontinuities in the radiance function and/or its derivatives (see
Figure 2.4). An analysis and classification of these discontinuities first appeared
in Heckbert's thesis [HECK91a]. Following Heckbert's definition:
A function L has a Dk discontinuity at x, if it is Ck-l but not Ck there.
The analysis that follows will be limited to planar polygonal non-interpene-
trating surfaces. This class of surfaces was chosen for a number of reasons:
1. The analysis is simpler than in the case of more general geometries.
2. The results of the analysis are immediately applicable to the design and
understanding of radiosity algorithms.
3. The vast majority of current radiosity systems are limited to this class of
surfaces and use polyhedral approximations to handle curved surfaces.
The diffuse rendering equation described by Equation (2.3) is a linear system
of equations. Each equation describes the radiance function over a surface s? as
the sum of simpler radiance functions, each one representing the contribution of
J4
Figure 2.3: Photograph of a simple scene. The light is a small square area source.
Figure 2.4: Discontinuity lines in a simple scene
15
a single surface in the environment. More formally:
where
Li(X) = L?j(X) + ? Lij(X)			(2.6)
SjES
Lq(x) =ft(X)!xtEsiLj(X') cosO$osOi v(x,x')dA(x')			(2.7)
denotes the contribution of surface &j to surface Sj. Due to the principle of su-
perposition, it is then sufficient to study the properties of a single Lij to give a
characterization of the radiance function Li.
The first part of the analysis will assume that surfaces Sj and sj do not touch
each other and that the radiance function Lj of the source Sj is smooth, i.e. C?.
These constraints will be lifted at end of the discussion.
Under these assumptions, the integrand in Equation (2.7) is a smooth func-
tion. Therefore, if the visibility of the source Sj is constant or varies smoothly,
the radiance function Ljj will be smooth as well. Abrupt changes in visibllity?
however, will introduce discontinuities of various orders into Lij. The rest of this
section will introduce the possible discontinuities and explain their causes.
2.4.1 D0 Discontinuities
By definition, D0 discontinuities are discontinuities in the function itself (value
discontinuities). They are introduced by edges or vertices of occluders (or the
light source itself) lying on the receiver Sj. Thus, they can occur either along line
segments, or in a pointwise fashion. An example is shown in Figure 2.5.
D0 discontinuities also occur along shadow boundaries cast by point light
sources. In this case no penumbra regions are present, and the transition from
the umbra into the unoccluded area is a discontinuous one. For conciseness, only
16
finite area light sources will be discussed. The analysis of point light sources is
only simpler.
In an environment with m edges, no interpenetrating polygons, and no point
light sources, there can be 0(m2) D0 discontinuities, as it is possible for each of
0(m) edges to lie on 0(m) faces.
2.4.2 D1 and D2 Discontinuities
Discontinuities in the first and second derivatives of the function are referred to
as D1 and D2 discontinuities, respectively. They are caused by the interaction
of a light source with the occluding objects intervening between the source and
the receiver, causing abrupt changes in the visibility of the former one. These
changes (also known as visual events) have been studied in the computational
geometry and computer vision literature [GIGU9O,GIGU91j. There, it is shown
that these visibility changes occur along surfaces defined by the edges and the
vertices of the objects in the scene (including light sources.) The visual events
in a polyhedral environment can be classified into two types: edge-vertex (EV)
events and edge-edge-edge (EEE) events [GIGU9O].
EV events occur on a subset of the plane defined by an edge and a vertex fully
or partially visible to each other. The set of points on that plane from which the
vertex can be seen coinciding with the edge is called the critical surface. Visibility
changes on a receiver occur along the curves of intersection between the receiv-
ing plane and the critical surfaces. For polygonal environments, these critical
curves are line segments in the case of an FV event. They correspond to either
or D2 discontinuities, as demonstrated in Figures 2.6 and 2.7. D' discontinu-
ities are perceived as Mach bands on the receivet Strong D2 discontinuities are
17
xo+
A
s
B
(a)
source is not visible
from x
1*1
x			y
(d)
Figure 2.5: An example of a D0 discontinuity. (a) Edge AB of the occluder 0 lies
on the receiver R. (b) The radiance function over R. (d) The radiance function
along the lIne through x and y.
Points on R immediately to the left of AB cannot see the source S and the radi-
ance there is zero. However, immediately to the right of AB, the entire source is
visible, and the radiance there is non-zero. Thus, the radiance function is discon-
tinuous along AB. The points A and B are points of singularity in the radiance
function.
18
Figure 2.6: A D1 discontinuity caused by an edge-vertex (EV) event. (a) Edge AB
of the light source S is coplanar to edge CD of the occiuding polygon 0. The
critical surface (shaded area) defined by A and CD intersects the receiving plane
along FE. (b) The radiance function over R. (c) The occluder 0 and the light
source S as seen from x. (d) The radiance function along the line through x and y.
From point y on R, none of the source is visible, hence the radiance there is zero.
As we move from y towards x, part of the source adjacent to AB becomes re-
vealed. The visible area grows linearly in the displacement from FE towards x.
Thus, along FE the radiance function has a D1 discontinuity
(c)
0
ci			ID
E			LF
R
(a)
visible portion of the source
A?B
x y
19
M
R
(a)
visible portion of the source
A
0
x y
Figure 2.7: A D2 discontinuity caused by an edge-vertex (EV) event. (a) The
critical surface (shaded area) defined by vertex A and edge BC intersects the re-
ceiving plane along DE. (b) The radiance function over R. (c) The occiuder 0
and the light source S as seen from x. (d) The radiance function along the line
through x and y.
From point y on R none of the source is visible, hence the radiance there is zero.
As we move from y towards x, part of the source adjacent to vertex A becomes
revealed. The visible area grows quadratically in the displacement from DE to-
wards x. Thus, along DE the radiance function has a D2 discontinrnty?
20
R
F
(a)
visible portion
of the source			F
c
E
2
-t
x y
(c)			(d)
Figure 2.8: An edge-edge-edge (EEE) event. (a) The quadric critical surface
(shaded area) defined by the three edges AB, CD, and EF intersects the receiver R,
resulting in a conic critical curve (shown dotted). (b) The radiance function
over R. (c) The occiuders Oi and O2 and the light source S as seen from x. (d) The
radiance function along the line through x and y.
As we move from y towards x, part of the source becomes revealed. A displace-
ment from the critical curve towards x results in quadratic growth in the visi-
ble source area1 hence the discontinuity along that curve is D2. This example
illustrates that when several occluding obstacles are involved, the boundaries
between umbra and penumbra regions on a receiver may be curved.
21
noticeable as well.
EFE events occur on a subset of a ruled quadric surface defined by the family
of lines that go through three skew edges, each visible to the others. The critical
surface is the set of all the points from which the three participating edges can
be seen intersecting in one point. The critical curves are conics and generally
correspond to D2 discontinuities as demonstrated in Figure 2.8.
In the examples shown in the figures all the visual events involve either a
vertex or an edge of the source. However, visual events can be defined by edges
and vertices belonging to any object in the environment. Thus, in an environ-
ment with m edges, there can be 0(m2) EV critical surfaces, and 0(m3) FEE criti-
cal surfaces. These can give rise to 0(m3) and 0(m4) critical curves, respectively
In general, the event will cause a discontinuity in Li1 if it is defined by edges
and vertices located between the receiver i and the source j and if it is visible
from both. Figure 2.9 shows the critical curves resulting on the receiver for the
examples depicted in Figures 2.5,2.6, 2.7, and 2.8.
2.4.3 Propagation of Discontinuities
The preceding discussion assumed that the radiance function Lj was smooth.
This is generally true for primary light sources, but, as the analysis above just
showed, secondary light sources can have radiance functions exhibiting discon-
tinuities of various orders. The effect of these discontinuities was summarized
by Heckbert [HECK91a] as the "Discontinuity Propagation Law", which states
that Dk discontinuities on the source can result in Dk+l and Dk+2 discontinuities
on the receivers. In general, discontinuities of higher orders are less noticeable
than low order ones.
22
(a)
(c)
(d)
Figure 2.9: Critical curves on receiver R for the simple scenes shown in: (a) Fig-
ure 2.5, (b) Figure 2.6, (c) Figure 2.7, and (d) Figure 2.8
23
P
Figure 2.10: Example of function singularity
2.4.4 Singularities
Radiance functions can exhibit singularities. These may rise at the endpoints of
D0 discontinuity lines where the radiance gradient becomes infinite.
Figure 2.10 shows an example of a singularity. The dashed lines identify the
umbra and penumbra boundaries of the shadow cast by the cube onto the floor.
A penumbra region lies between the two boundaries and provides a continuous
transition between the unoccluded and totally occluded regions of the floor. As
we move in towards the cube, though, the width of the penumbra gets smaller
and smaller and the gradient of the radiance function becomes larger and larger.
In the limit, under the corner of the cube, the gradient becomes infinite.
This corresponds to a natural phenomena we have all experienced in the en-
vironment around us: shadow boundaries are sharper close to the object casting
24
the shadow and get fuzzier moving away from it. Furthermore, the quality of
the shadow in the immediate vicinity of the cube helps us determine whether
or not the cube is actually in contact with the floor. Thus, correctly reproducing
these singularities increases the realism of our simulations.
Chapter 3
Previous Work
3.1 Radiosity
The radiosity approach is a physically-based method that attempts to solve the
global illumination problem according to the principles of transfer and conser-
vation of energy; it was derived from techniques used in thermal engineering
for the determination of radiative heat exchange in enclosures [SIEC8l] and was
first introduced to the computer graphics community by Coral et al. [CORA84j.
Since its first appearance in 1984, radiosity methods have been extended to
occiuded environments [COHE85], textured surfaces [COHE86], specular sur-
faces [IMME86,WALL87,SILL89j, light sources of arbitrary directional distri-
bution [WALL89,DORS91], arbitrary bidirectional reflectance distribution func-
tions [51LL91], dynamic environments [BAUM86,CEOR9O,CHEN9Ob], curved
surfaces [ZATZ92], and complex environments [COHF88b,HANR9lb,SMIT92].
3.1.1 Diffuse Radiosity
Most radiosity algorithms are restricted to ideal diffuse reflectors and attempt
to provide a view-independent solution to the global illumination problem. Thus,
25
26
JJwaW?though5fl of static diffuse environments can be often generated at inter-
active rates from a single radiosity solution.
Radiosity algorithms use a discrete formulation of the diffuse global illumi-
nation problem. The traditional radiosity equation and its derivation were given
in Section 2.3. When the environment is subdivided into a set of small surface ar-
eas, or elements, of constant radiosity and constant reflectivity, the radiant energy
transfer between surfaces is governed by a simple system of linear equations.
This can be expressed in matrix notation as:
where
(I--HM)b=e			(3.1)
M ? ????, M = [M4] = [piFij], is the form factor matrix;
b E ??, b = [Bi]1 is a vector of unknown element radiosities;
e E ??, b = [Be11, is a vector of emitted element radiosities.
The matrix I --H M is a so-called `,M-matrix," i.e. all the diagonal elements are
positive and all the off-diagonal elements are negative. Furthermore, for phys-
ical reflectivities (p? ? [0,1]), the matrix is also strongly diagonally dominant.
Many nice properties are known about these matrices [VARG62], including non-
singularity and the convergence of Jacobi and Gauss-Seidel iterations.
3.1.2 Full Matrix Solution
The early radiosity algorithms were organized in three major steps:
1. Compute the matrix M of form factors.
27
2. Solve forb = (I--HM)-1e.
3. Display the results.
The first radiosity implementation in computer graphics [CORA84] used Gaus-
sian elimination to invert the matrix I --H M. Later, Cohen et al. [COHE85] greatly
reduced the cost of Step 2 by using Causs-Seidel iteration to solve for the un-
known radiosities. Unfortunately, though, solving the system of linear equations
is not the most expensive step in the process. Computing the form factors be-
tween all pairs of elements in an environment, even when employing a fast hard-
ware z-buffer for the hemicube form factor calculation algorithm [COHE85], is
very computationally intensive. Furthermore, the cost of storing the form factor
matrix, even when employing compression techniques to exploit the sparsity of
the matrix, is prohibitive for anything but very simple scenes.
Precomputing the entire form factor matrix is not strictly necessary. The form
factors could be computed one row at a time during the Gauss-Seidel iteration.
This strategy, though, would require recomputing each row of form factors at
each iteration and thus greatly increase the computation time.
3.1.3 Progressive Refinement Radiosity
The shortcomings of the full matrix approach were overcome by the introduc-
tion of progressive refinement radiosity [COHE88b]. This algorithm is based on
a refinement operator that uses a single column of the matrix to update the ra-
diosities of all elements in the environment with the contribution from a chosen
source.
The difference between conventional radiosity and progressive radiosity is
28
best explained in terms of gathering and shooting. With Gauss-Seidel iteration,
each pass over row i of the form factor matrix results in computing a new esti-
mate of radiosity Bj based on the current estimates for all the other radiosities in
b. In physical terms, the estimate Bi is computed by "gathering" in the radiant en-
ergy contribution from the rest of the environment. With progressive refinement
radiosity, each pass over column j of the form factor matrix results in computing
a new estimate of all the radiosities in b based on the current contribution from
element j. In effect, element j "shoots" its energy out into the environment.
Since each refinement operation effects the radiosities of the entire environ-
ment, useful intermediate results can be displayed even at the very beginning of
the solution process. The initial rate of convergence can be improved greatly by
shooting at each iteration from the element with the highest remaining energy
contribution rather than visiting each column in a strictly sequential order.
This algorithm has the same asymptotic convergence rate as Jacobi and Gauss-
Seidel iterations. However, by intelligently observing that the global illumina-
tion in a scene is mostly determined by the dominant effect of primary and a
few secondary sources, progressive radiosity comes close to the final solution in
a fraction of the time that it would take for convergence. In practice, the iteration
process is seldom carried to convergence and only a small fraction of the form
factor matrix is ever computed.
3.2 Meshing
Most radiosity algorithms to date attempt to capture the variations in illumi-
nation across a surface in a scene by discretizing the surface into pieces small
29
enough that the radiance function across each piece can be closely approximated
by a simple polynomial [GORA&4]. This process is referred to as meshing.
Meshing is perhaps the most problematic module of any radiosity system
and is responsible for many of the visual artifacts typical of radiosity images.
The next section will use uniform meshing as a way of example to illustrate these
problems. The following sections will discuss different meshing algorithms and
their ability to eliminate or reduce such problems.
3.2.1 Uniform Meshing
The first radiosity algorithms introduced by Goral et al. [GORA84], Cohen and
Greenberg [COHE85], and Nishita and Nakamae [NISH85a] were based on uni-
form meshing. Each surface in the environment was meshed in parametric space
into a set of rectangular elements by splitting each surface along the lines of a
regular grid. Radiance samples were then taken at the center of each element.
This simple technique is easy to implement, but, as many researchers pointed
out [HAIN9?b,BAUM91,WALL92], it results in a large variety of visual artifacts:
Missing shadows. If the shadow cast by an obstacle in the scene onto a receiving
surface falls within sample points, it will be entirely missed (see the legs of
the table closest to the back wall in Figure 3.1.)
Blocky shadows. If a shadow area is undersampled, its shape will be distorted
to resemble the underlying mesh structure. Also, due to the interpolation
methods used for reconstruction, penumbra areas will usually expand to
the sample points immediately outside the true penumbra area (see, for
30
Figure 3.1: Visual artifacts due to uniform meshing. Missing shadows (the legs
of the table closest to the back wall.) Blocky shadows (the shadow of the table on
the floor and back wall.) Jagged shadow boundaries (the shadow of the leftrnost
leg of the table.) Shadow leaks (the shadow under the cube and behind the table.)
Floating objects (the legs of the table closest to the back wall.) Mach banding (the
region of the floor under the table.)
31
Figure 3.2: More visual artifacts. Shadow leaks (the shadow under the doorway.)
Seams (the region of the wall above the doorway.)
32
example, the shadow of the table on the floor and back wall in Figure 3.1.)
Reducing the element size will reduce, but not eliminate, these problems.
Jagged shadow boundaries. Unless the shadow boundaries are parallel to one
of the two axis of the mesh, shadow boundaries will appear jagged (see
the shadow of the leftmost leg of the table in Figure 3.1.) This aliasing can
be eliminated, for a given view, by ensuring that the projected area of each
mesh element be smaller that that of a screen pixel.
Shadow/light leaks. If a sharp variation in the illumination falls between sam-
ple points, the reconstruction process will smooth it and stretch it, dis-
tributing the variation uniformly along the area between the two samples.
This can result in dark shadow regions "leaking" into bright illuminated
regions and vice versa (see the shadow leaks under the cube and behind
the table in Figure 3.1 as well as that under the doorway in Figure 3.2.)
Floating objects. Missing shadow and shadow and light leaks can cause objects
to look as if they were floating above the floor (see the legs of the table
closest to the back wall in Figure 3.1.)
Mach banding. Most radiosity algorithms to date use linear interpolation to re-
construct the radiance function over a surface. This method introduces
unwanted discontinuities in the illumination gradient along the mesh el-
ement boundaries. These discontinuities are perceived by our visual sys-
tem as Mach bands (see the Mach bands on the region of the floor under the
table in Figure 3.1.) Increasing the mesh resolution usually helps reduce
these artifacts. However, better reconstruction methods could be used, e.g.
33
higher order interpolation, that would greatly reduce Mach banding with-
out needing to increase the sampling rate.
Seams. A complex surface can be modeled as a set of polygons. Fach polygon
is then independently meshed into elements. If the polygon meshes are
not aligned or have different granularities, radiance discontinuities may
result along the boundaries between neighboring polygons (see the region
of the wall above the doorway in Figure 3.2.) The problem of misaligned
meshes has been thoroughly examined by Haines [HAIN9?bj, and Baum
et al. [BAUM9?], The interested reader is invited to refer to their work for
an overview of the possible solutions.
Missing and incorrect shadows, jagged shadow boundaries, shadow/light
leaks, and floating objects are all due to undersampling of the radiance function.
These artifacts occur in areas where the illumination gradient is highest and thus
undersampling is most noticeable.
Increasing the sampling rate by reducing the mesh element size will reduce
the artifacts at the cost of increased storage and computing time. Since regular
sampling is used, increasing the sampling rate to avoid undersampling shadow
boundaries can also lead to oversampling in regions that are uniformly illumi-
nated and where the radiance varies slowly, thus resulting in unnecessary com-
putations.
Deciding what the best mesh resolution is for each surface in a scene is gener-
ally generally left to the user. Generating high quality radiosity solutions takes
many tedious and time consuming iterations and requires the user to have an
unnecessarily high understanding of the workings of meshing and radiosity al-
34
gorithms. Despite these manual adjustments, the resulting radiosity solutions
are generally accurate for only a small range of views and thus, are not truly
view-independent. For example, zooming in on a shadow is likely to reveal a
jagged boundary.
3.2.2 Adaptive Subdivision
Cohen et at. [COHE86] introduced an adaptive subdivision scheme for refining
the mesh resolution in regions where the radiance function displays high inten-
sity gradients. Their method meshes the environment into a two-level hierarchy:
the first level consists of a coarse subdivision of the environment into patches,
the second level further subdivides patches into smaller elements. First, all the
patch-to-patch form factors are computed using the hemi-cube method [COHE85]
and a system of linear equations is solved to yield the patch radiosities. Then
element-to-patch form factors are computed and each element radiosity is cal-
culated as a simple linear combination of the patch radiosities. If the radiosity
gradient over an element exceeds a user-specified threshold, the element is sub-
divided into a set of new smaller elements. The radiosities at these new elements
are then calculated exactly as done previously for the other elements. Notice that
since the subdivision of an element does not affect the geometry of the environ-
ment, the radiosities of all the other elements do not change. This subdivision
process can be applied repeatedly to refine the quality of the solution.
This substructuring scheme represents an important step forward over the
full-matrix solution method based on uniform meshing [CORA84,COHE85].
Substructuring makes it possible to discretize the environment into many more
elements than could ever be handled by the previous scheme, both because of
35
the cost of computing all the element-to-element form factors and because of the
enormous memory required to store the full element-to-element matrix.
Provided the size of the initial patches is not too coarse, this adaptive mesh-
ing scheme automatically refines the mesh around umbra and penumbra bound-
aries without unnecessarily refining uniformly illuminated areas. This approach
reduces both the need for user intervention and the time and memory costs re-
quired to produce an accurate radiosity solution. Most of the problems related
to uniform meshing are somewhat reduced, but none are truly eliminated.
The two-level hierarchical organization of the mesh and the adaptive subdi-
vision scheme have also been used within the context of progressive refinement
radiosity [COHE88b,WALL92].
3.2.3 Hierarchical Meshing
Hanrahan et al. [HANR9O,HANR9lbj generalized the substructuring approach
to a full hierarchical meshing scheme. The surface mesh is stored as a quadtree
and each level of the quadtree represents a different level of subdivision. The
distinction between source and receiver is blurred and energy transfers can take
place between nodes at any level in the respective hierarchical meshes of a source-
receiver pair (see example in Figure 3.3.)
If energy is transferred to an internal node, this energy is "pushed down" to
the leaves of the structure by recursively distributing the energy to the node's
children in proportions determined by their relative areas. The energy stored at
the leaves is then `6pulled up" and the energy of each interior node is computed
as an area weighted average of its children's energies.
Transferring energy between large patches (nodes close to the root of their re-
36
Figure 3.3: Hierarchical meshing on a pair of perpendicular polygons. The line
segments show the interactions between pairs of elements at different levels in
the quadtree hierarchies (after Hanrahan [HANR9lb].)
37
spective hierarchies) reduces the number of source-receiver interactions needed
to produce a radiosity solution; transferring energy between small patches re-
duces errors associated with the computation of form factors and visibility terms.
Hanrahan et al. compute an upper bound on the error incurred in transferring
energy between two nodes. Energy between a source-receiver pair is then trans-
ferred by subdividing source and receiver until all the interactions fall within
the error bound.
With this scheme, the number of interactions for an environment of m patches
and n elements (leaf nodes) is reduced from O(mn) to only 0(n). Furthermore,
this meshing scheme is completely automatic.
On the other hand, however, requiring that all energy transfers occur in quan-
tities smaller than a given threshold limits the efficiency of the algorithm. Also,
the errors incurred gathering energy from different patches can compound, but
no error bound for the final solution was shown. Furthermore, although the ra-
diosity solution is quantitatively accurate, because of the point sampling nature
of the visibility calculations and of the constant elements used to subdivide the
input surfaces, the images computed with this technique can still exhibit all the
visual artifacts produced by its predecessors.
3.2.4 Discontinuity Meshing
All of the meshing techniques discussed above mesh the input surfaces accord-
ing to local geometric considerations only. Therefore, shadow boundaries and
other sharp variations in the illumination will generally fall within the inte-
rior of elements; since reconstruction is usually done by interpolating the radi-
ance values at the element vertices, these illumination details will be distorted,
38
smoothed, or even entirely missed.
These problem can be eliminated by explicitly representing discontinuities in
the radiance function and its derivatives as boundaries in the mesh. Discontinu-
ity meshing methods implement this idea and can greatly improve the quality of
radiosity images1 even though these methods only capture a subset of the dis-
continuities described in the Section 2.4.
Capturingvaluediscontinuities
Baum et al. [BAUM9l] preprocess the environment to split surfaces along lines of
intersection or contact. Thus, the corresponding D0 discontinuities are marked
by edges in the mesh. The radiance samples stored with these edges are the
values corresponding to the illuminated side of the edge. The problem of storing
both the occluded and unoccluded radiance values is avoided by making all D0
edges lie at the perimeter of the mesh. If a surface A lies on another surface
B, the entire area covered by A is dropped from the mesh of B (see example in
Figure 3.4.)
This approach has several advantages over previous meshing techniques:
o+ The quality of radiosity images is significantly improved by eliminating
floating objects and the most obvious shadow and light leaks.
o+ The need for user intervention is reduced since value discontinuities are
captured automatically, independently of the initial size of the mesh ele-
ments.
o+ Time and storage costs are reduced since adaptive subdivision is no longer
needed to narrow in on value discontinuities.
39
(a)
(b)
Pigure 3.4: Computing D0 discontinuities prior to meshing. (a) The initial en-
vironment has object A resting on surface B. (b) The area covered by object A
is removed from surface B; the resulting D0 edges are shown as thicker lines
(object A is not shown.)
40
This meshing algorithm, however, is still liable to miss small illumination
details and generate shadows with imprecise coverage and poor definition.
Capturing umbra and penumbra boundaries
Campbell and Fussell [CAMP9O] suggested meshing the environment by pro-
jecting each object onto each surface in the environment from the light source
(see example in Figure 3.5a.) For point light sources, shadow volumes were used
to mesh the environment. Finite area light sources were discretized into smaller
parts and each part was approximated by a point light source. This method
works well for small light sources, but as the number of point sources needed
for an adequate approximation grows, subdivision within regions of penumbra
becomes too fine.
Recently, Campbell and Fussell [CAMP9la,CAMF9lb,CAMP9lc] have pro-
posed an alternative method for polygonal environments whereby umbra and
penumbra volumes are used to classify each surface into regions of umbra, penum-
bra, and total visibility (see example in Figure 3.5b.) Each region is further split
into elements, and the illumination at the vertices of each element is computed
analytically. Numerical optimization techniques are then used to determine the
correct element density inside unoccluded and penumbra regions.
The surface mesh is stored as a two-dimensional BSP tree [FUCH8O] with the
mesh elements as the leaves of the tree (see examples in Figure 3.5.) Since vertices
are not shared among elements, illumination calculations may be recomputed
several times at each location; on the other hand, because of this organization,
D0 edges do not require any special treatment.
Umbra and penumbra volumes were previously used by Nishita and Naka-
41
s
(a)			(b)
0
10
Figure 3.5: Shadow volumes (a) and umbra and penumbra volumes (b) con-
structed from source S and object 0 are used to classify surface R into regions of
umbra, penumbra, and total visibility. The resulting meshes (c-d) are stored as
two-dimensional BSP trees.
42
mae [NISH83,NISH85aj to classify regions in a scene according to light source
visibility. However, this information was not used to mesh the environment in
object-space, but rather to speed up the illumination calculations within a screen-
space, scanline-based renderer (see discussion in Section 3.3.2.) In fact, Camp-
bell's approach is, to the author's knowledge, the first object-space algorithm to
accurately compute shadows cast by area light sources. Many of the visual arti-
facts typical of previous radiosity algorithms are eliminated without any manual
intervention. All the shadows are captured, no matter how small, umbra and
penumbra regions have the correct shape and sharp straight boundaries, and
floating objects and shadow and light leaks are eliminated.
However, despite the impressive improvements in image quality, a number
of problems still remain:
o+ Umbra and penumbra volumes only capture a subset of the FV events re-
sponsible for the discontinuities on a receiving surface. More precisely,
only umbra and penumbra boundaries are computed; discontinuities falling
within a penumbra area are ignored.
o+ As Teller pointed out [TELL92], shadow and discontinuity meshing algo-
rithms based on umbra and penumbra volumes only take into considera-
tion visual events originating from the interaction of a light source and a
single polygonal occluder. Since multiple interactions are ignored, some
umbra areas may be incorrectly classified as penumbras.
o+ For the same reason, EFE events are ignored. The discontinuity curves
generated by these events are conics and can not be correctly represented
as edges of a polygonal mesh.
43
o+ Mach bands can still appear in the resulting images since linear interpola-
tion is used for reconstruction.
o+ The problem of radiance singularities is still unaddressed.
Discontinuity meshing in 2D
The idea of including discontinuities in radiance functions and their derivatives
as boundaries in the mesh was developed independently by Heckbert [HECK9la,
HECK92c] and by Lischinski and the author at Cornell Program of Computer
Graphics [LISC9I]. In both cases, the ideas were implemented and tested in 2D,
showing promising results.
Discontinuity meshing in 3D
Recently, both Heckbert [HECK92a] and Lischinski, Tampieri, and Greenberg
[L1SC92j extended their ideas to 3D polygonal environments.
Heckbert [HECK92aj constructs a discontinuity mesh that contains all the
discontinuity lines arising from EV events in which primary light sources di-
rectly participate. First, the critical lines are computed on all the surfaces in the
environment, and then Delaunay triangulation and mesh relaxation are used to
produce the final discontinuity mesh on each surface.
Lischinski et al. [LISC92] construct a different discontinuity mesh for each
light source. Each mesh contains the discontinuity lines arising from EV events
in which that light source directly participates. A two-dimensional BSP-tree
coupled with a topological data structure is used to construct and represent the
mesh. Adaptive subdivision is used to increase the element density where neces-
sary. The contributions due to different light sources are merged into a separate
44
data structure which accumulates the total illumination over each surface. A
detailed presentation of this method can be found in the next two chapters.
Both algorithms still ignore EEE events and only consider VE (or EV) events
involving a light source. As opposed to Campbell's algorithm [CAMP91c], these
discontinuity meshing methods correctly capture discontinuities falling within
penumbra regions.
Furthermore, in the Cornell work higher order interpolation schemes are
used to reconstruct the radiance function and correctly resolve the different or-
ders of discontinuity along mesh edges. This technique reproduces the correct
degree of sharpness of various illumination details, effectively resolves radiance
singularities, and eliminates the annoying Mach banding typical of linear inter-
polation schemes.
3.3 Dfrect Energy Transfer
The ability to compute the light energy contribution of one mesh element to an-
other is one of the most important building blocks of any radiosity algorithm.
As was shown in Equation (2.7), this computation requires evaluating a com-
plex integral. Analytical solutions have only been offered for very simple cases,
such as that of a source of constant emission, and in the absence of any occluding
obstacles.
In order to simplify the problem, traditional radiosity algorithms have dis-
cretized the diffuse rendering equation into a simple system of linear equations
(see Equation (2.5)). Then, the contribution of element Sj to the illumination of
45
element Sj becomes:
where
B1 = piBjF4			(3.2)
Fij= $iIAf)V(X?X') cos ($??50idA4Ai			(3.3)
This formulation relies on the assumption that the radiance must be constant
across each element in the scene. If this assumption is met, the problem of com-
puting the energy transfer between two elements is reduced to that of computing
a form factor, which, in the case of Lambertian reflectors, is a purely geometrical
term.
3.3.1 Numerical Form Factors
The computation of form factors is complicated by the presence of a visibility
term under the integral. Considerable research efforts have been devoted to the
investigation of accurate and efficient methods for computing form factors.
Goral et ai. [GORA84], who introduced radiosity to the computer graphics
community, transformed the area integrals to contour integrals and evaluated
them numerically by discretizing the contours into a set of short line segments.
This technique was computationally expensive and only simple unoccluded en-
vironments were demonstrated.
Cohen and Greenberg [COHE85] were the first to present an efficient algo-
rithm for computing form factors in the presence of occiusions. Their hemicube
algorithm was a numerical approximation of the geometric analog for the form
factor integral developed byNusselt [SIEC8l]. The inner integral of Equation (3.3)
was approximated by a weighted sum of the hemicube pixels covered by the
46
Th
Sj
Figure 3.6: The hemicube algorithm. The form factor from surface s? to the cen-
ter of surface Sj is approximated by the weighted sum of the hemicube pixels
covered by &j's projection.
47
projection of element &j relative to the center of element Sj (see Figure 3.6.) The
outer integral was simply dropped by assuming that the distance between the
two elements is large compared to their size and that the visibility of &j is inde-
pendent on the location on Sj (constant visibility.) The hemicube algorithm can
be implemented using a hardware z-buffer and is thus could be very efficient.
Unfortunately, the introduction of progressive refinement radiosity (COHE88b]
and its use on larger, more complex, environments, uncovered many problems
and inaccuracies connected to the use of the hemicube. As Baum et al. [BAUM89j
pointed out, in order to accurately compute form factors with the hemicube al-
gorithm, a number of assumptions must be met. The distance between elements
must be large compared to their area, the visibility of one element must be the
same as seen from any point on the other element, and the area of the projec-
tion of one element onto the hemicube must be correctly represented by an inte-
gral number of hemicube pixels. The full matrix solution projected patches onto
hemicubes centered over mesh elements. Since the patches were large compared
to elements and the environments being solved were generally very simple, the
hemicube computations were relatively accurate. In progressive refinement ra-
diosity, however, the role of patches and elements was reversed. The Hemicube
was placed on large patches, often violating the constant visibility assumption,
and small elements were projected onto the hemicube, causing significant alias-
ing problems and resulting in the patch-to-element form factors being severely
over/under-estimated. Furthermore, progressive radiosity made it possible to
solve more complex models than could be handled by the full matrix approach,
thus making the hemicube artifacts even more apparent.
48
dA ?
Figure 3.7: Ray-traced form factors. Surface s? is subdivided into smaller regions;
the contribution from each region AAj to the form factor dFAj,dAj is approximated
by a disc-to-point form factor.
When elements are close to each other or when one of them is a primary
light source, Baum et al. [BAUM89] proposed replacing the hemicube with an
analytical formula for the computation of the inner integral [HOTT67] coupled
with a numerical evaluation of the outer contour integral. This method works
well for unoccluded form factors. Partial visibility is resolved by recursively
subdividing the elements and adding up the new smaller form factors.
A different solution was offered by Wallace et al. [WALL89]. Rather than
using the traditional finite element formulation of radiosity [CORA84], the au-
49
thors proposed a point sampling approach. Thus, instead of computing the en-
ergy transfer between pairs of elements, the authors computed the radiance con-
tributed by element &j to a differential area on a receiving surface s?;
where
Bj(x)dAj(x) = PiBjAjdFAj,dAj
(3.4)
cos O? cos 0
dFAj?dAz =			v(x, x')			7rr2 ?dAj(x')			(3.5)
A1
The integral was approximated by the formula for a disc-to-point form factor and
adaptive subdivision of element sj was used to reduce the error (see Figure 3.7.)
Visibility was estimated using ray tracing.
The ray-traced form factor method is efficient and flexible. While the hemicube
method requires that a complete row or column of form factors be computed at
once, with the same resolution, the other method allows individual control on
the accuracy and computational cost of each form facton
3.3.2 Analytical Form Factors
All of the methods presented above, though to different degrees, are subject to
aliasing due to the point sampling nature of the visibility computation. Nishita
and Nakamae [NISH85a] proposed an analytical, exact method of computing
the contribution from a source of constant emission to a point. The algorithm
was used within the context of a scanline rendering system and was applied to
points spaced at intervals of a few pixels across each scanline. The algorithm
started by building umbra and penumbra volumes for each source-obstacle pair
so that each point on a receiving surface could be easily classified as being in the
umbra1 penumbra, or unoccluded region with respect to the light source. The
50
(a)
Oi
:0
(c)
(b)
Figure 3.8: Analytical form factors. (a) Contour integration along the path indi-
cated by the arrows is used to compute the contribution of source &j to point x.
(b) Only the part of the source lying in front of the receiver is considered in corn-
puting the energy transfer. (c) Obstacles intervening between the source and
receiver affect the integration path.
illumination at an unocciuded point was then computed analytically using the
contour integral method [HOTT67] (see Figure 3.8a.) If the light source polygon
was intersected by the supporting plane of the receiving point, the illumination
was computed by integrating along the contour of the part of the light source
lying in the front of the receiving plane (see Figure 3.8b.) If the receiving point
was in the penumbra area, the illumination was computed by integrating over
segments of the contours of the source and obstacles in a way that amounted
to integrating across the portions of the source that were visible from the re-
ceiving point (see Figure 3.8c.) Computing the illumination at points lying in
the penumbra was expensive but yielded the correct result. No approximations
were made.
3.3.3 Non-Constant Exitance
In the preceding algorithms, the expression relating the energy contribution of
one element to another element or differential area is based on the assumption
that the exitance of the source be constant. This assumption is often violated,
especially due to the subdivision of the input surfaces into patches and elements
first introduced by Cohen et al. [COHE86] and now adopted by most radiosity
systems. While elements are often small enough to satisfy the constant radiosity
assumptions, patches, being clusters of (possibly many) elements, often violate
this premise. Furthermore, since patches, rather than elements, are often used
as a source, a number of problems may result.
These problems were discussed by Tampieri and Lischinski [TAMP9l], who
also presented a simple solution, based on ray-traced form factors [WALL89]
(see Figure 3.7.) The source patch was recursively subdivided into smaller re-
52
gions and ray-traced form factors were used. Rather than assuming constant
emission across the entire patch, though, the underlying element mesh was in-
terrogated to compute separate radiance values for each region. This method
has all the efficiency and flexibility of ray-traced form factors, while at the same
time providing an accurate treatment of arbitrary radiance distributions.
The scheme of Hanrahan et al. [HANR9O,HANR9lb] is also immune to the
problems connected to violations of the constant radiosity assumptions. By us-
ing the hierarchical mesh and subdividing interacting pairs of elements until the
energy transferred between the two falls below a given threshold, a source of ar-
bitrary distribution can be subdivided until each piece can be safely treated as a
constant emitter.
3.4 Radiance Function Reconstruction
Radiosity algorithms typically compute radiance values only at selected sample
locations on the surfaces of the scene. Reconstruction techniques are then used
to provide radiance values at any desired location.
Radiance functions must be reconstructed for display and, in the case of pro-
gressive refinement radiosity, used for selecting the patch with the highest un-
shot radiant flux. The former task requires that the reconstructed radiance func-
tions reproduce as many of the characteristics discussed in Section 2.4 as possi-
ble. The latter task only requires a rough approximation since only the integral
of the radiance over a surface is needed and choosing the wrong patch would
not affect the accuracy of the solution anyway.
The traditional radiosity formulation presented in Equation (2.5) discretizes
53
the continuous domain of the diffuse global illumination problem into a finite
set of elements using constant shape functions. Higher order shape functions,
though, should be used at reconstruction time because they yield more realistic
pictures.
3.4.1 Constant Elements
When using constant shape functions, the mesh elements are displayed as flat
shaded polygons. From a subjective, or visual, standpoint, this piecewise con-
stant approximation is usually very poor because the value discontinuities in-
troduced by the reconstruction process are very easily perceived by our visual
system as sudden jumps11 in the illumination of a surface. These artifacts reveal
the artificial discretization of the surface into elements.
3.4.2 Gouraud Shading
In order to improve the visual appearance of radiosity solutions, Coral et al.
[GORA84] suggested using Couraud shading [COUR7l]. First, vertex radiosi-
ties are computed from the element radiosities, then linear interpolation in screen
space is used to compute the radiosities corresponding to the locations within
mesh elements that map onto pixel centers.
This technique yields significantly better results that simple flat shading, but
is not without problems. The reconstructed radiance functions are now con-
tinuous everywhere across a surface, but exhibit gradient discontinuities across
element boundaries, which are perceived as Mach bands [RATL72j. Further-
more, as Duff pointed out [DUFF79], when interpolating in screen space, the
shading across a surface may appear to change with the view due to both the
54
non4inear mapping of distances effected by the viewing projection and the axis-
dependence of the Gouraud shading algorithm for polygons other than trian-
gles.
Despite these drawbacks, this technique is widely adopted in current radios-
ity systems because Gouraud shading is now available in hardware on most
graphics workstations.
3.4.3 Bilinear Interpolation
Cohen and Greenberg [COHE85] reconstructed the radiance function using bi-
linear interpolation in object space. Unlike Couraud shading, this method yields
a unique value at every point on an element, independent of the view and orien-
tation of the element, Bilinear interpolation, though, like Gouraud shading, in-
troduces unwanted first derivative (slope) discontinuities across element bound-
aries.
These techniques require radiance values at the element vertices, but, in fact,
when they were first proposed [GORA&4,COHF85], radiosity systems computed
radiance values at the center of the elements. Cohen and Greenberg [COHE85]
computed vertex radiance values as a simple average of the radiances of the
surrounding elements. This technique works reasonably well for regular rectan-
gular meshes, but for triangular or arbitrary quadrilateral meshes its not at all
clear what choice of weights should be used in computing the averages.
Furthermore, reconstructing the radiance function from the radiance samples
now takes two interpolation stages, each one introducing further errors and ap-
proximations into the solution. A simple solution to this problem was offered.by
Wallace et al. [WALL89], who proposed computing the radiance values directly
55
at the element vertices, thus reducing the reconstruction process to a simple in-
terpolation step.
3.4.4 Higher Order Interpolation
Higher order interpolation schemes have been proposed as an alternative to
linear interpolation. Reichert [REIC92j described a C1 quadratic interpolation
scheme, based on work done by Cendes and Wong [CEND87], that imposes first
derivative continuity across the mesh edges and thus greatly reduces visual ar-
tifacts. Unfortunately, this technique also has the effect of washing out sharp
illumination details such as umbra and penumbra boundaries.
All of the interpolation methods discussed above yield a reconstructed ra-
diance function that is infinitely differentiable inside each individual mesh ele-
ment. Thus, no shadow boundary crossing an element can be correctly recon-
structed. This results in the jagged shadow boundaries, expanded or exagger-
ated penumbrae, floating objects, and shadow and light leaks, so typical of ra-
diosity images.
As first demonstrated by Campbell and Fussell [CAMF9O], these problems
can be greatly reduced and sharp umbra and penumbra boundaries can be con-
vincingly rendered if the surfaces in the scene are split into mesh elements along
these boundary lines.
Chapter 4
Single Source Discontinuity Meshing
The previous chapter showed that the conventional meshing algorithms used
by most radiosity systems to date are the cause of many inaccuracies in radios-
ity solutions. In particular, umbra and penumbra boundaries and other sharp
variations in the illumination across a surface are poorly reproduced if not en-
tirely missed.
The analysis of radiance functions presented in Section 2.4 showed that these
illumination changes correspond to discontinuities in the radiance function and
its derivatives and occur along critical curves determined by the geometric in-
teraction of the objects and light sources in the environment.
Conventional radiosity algorithms ignore these discontinuity curves. If a crit-
ical curve crosses a mesh element, the reconstructed radiance will be inaccurate.
In fact, the piecewise polynomial interpolation methods used to reconstruct the
radiance function over a surface mesh are unable to reproduce a discontinuity
that crosses a mesh element.
The ability to predict the location of critical curves and the order of the cor-
responding discontinuities can be used to design a new algorithm less prone
56
57
to the problems that afflict conventional radiosity. The new approach explic-
itly represents the discontinuities in the radiance function and its derivatives as
boundaries in the mesh and uses a piecewise polynomial interpolation scheme to
correctly reproduce these discontinuities in the reconstructed radiance function.
This concept, known as discontinuity meshing, was first introduced by the au-
thor and Lischinski at the Cornell Program of Computer Graphics [LISC91] and,
independently, by Heckbert [HECK91a,HECK92c].
In the succeeding two chapters a radiosity algorithm based on discontinuity
meshing and capable of producing accurate solutions to the diffuse global illu-
mination problem for three-dimensional polygonal environments is discussed.
The presentation is organized in two main parts: the rest of this chapter dis-
cusses a discontinuity meshing algorithm that computes the direct illumination
due to a single primary light source of finite area and constant emission. The
next chapter extends the basic algorithm to secondary light sources of arbitrary
emission and uses the result as a refinement operator within a progressive-style
radiosity algorithm to produce a diffuse radiosity solution.
4.1 Single Source Discontinuity Meshing--HAlgorithm Overview
The direct illumination problem that we want to solve is stated more rigorously
as follows:
Given a set S=fsi,s2,... ,s?1 of convex polygonal surfaces, their dif-
fuse BRDF'sfri,fr2,... ,frn, and a convex polygonal light source So of
constant emission L0, compute the set of radiance functions L10, ....... ,Ln0
58
foreach visual event Vk involving So do
trace Vk through S and insert resulting critical curves in the meshes
end for
fors? e Sdo
triangulate mesh for Sj
foreach element tk in mesh for Sj do
compute L10 at sample locations within tk
adaptively subdivide tk as needed
end for
end for
Figure 4.1: Pseudocode for the single source discontinuity meshing algorithm
such that
Lj0(x) =fnL0!,Eso cosOi cos00v(x, x')dA(x')
r2
(4.1)
for all i = 1,2,..., n.
Of course, solving Equation 4.1 for every point x on surface 5? is not feasi-
ble. Instead, we use a sampling approach whereby the radiance function LjO
is evaluated at a few selected locations and then reconstructed using a suitable
interpolant to yield an approximation Li0.
In computer graphics, this kind of approach has roots that go back to the ear-
liest shading algorithms such as Couraud shading [GOUR71j. The originality
of the algorithm described here, though, rests in its ability to predict the loca-
tion and order of the radiance discontinuities, and use these results to guide the
choice of sample locations and reconstruction interpolant.
Pseudocode for the single source discontinuity meshing algorithm is given
in Figure 4.1. The algorithm starts by tracing all the visual events involving
source So through the environment S. Whenever a visual event results in a criti-
59
cal curve on a surface Sj, one or more edges are created in the discontinuity mesh
for Sj. These edges store the location of the critical curve as well as the order of
discontinuity associated with it.
Once all the visual events involving source ?o are processed, each surface in
the environment has been discretized into a discontinuity mesh where each im-
portant discontinuity in the radiance function and its derivatives is explicitly
represented by one or more edges.
At this point, the mesh is triangulated and the direct illumination due to
source 50 is computed at a set of selected locations within each mesh element.
Once sample values are computed over an element, an error estimate is com-
puted for that element, and, if necessary, the element is adaptively subdivided
to yield a more accurate approximation to Lj0. Since no discontinuity is allowed
to cross an element, though, the radiance function inside each element is well
behaved and can in most cases be captured by a relatively small number of sam-
ples.
Finally, a piecewise quadratic interpolant, Lj0, is used to reconstruct an ap-
proximation to the actual radiance function Lj0. The interpolation scheme pro-
vides smooth behavior inside mesh elements while still preserving the appro-
priate order of discontinuity across each mesh edge.
The following sections describe each of the steps outlined above in greater
detail and present the data structures and techniques used for an efficient im-
plementation of the algorithm.
60
4.2 Locating the Discontinuities
In Section 2.4, we showed that the direct illumination due to a single area light
source can cause value (zero) and first and second derivative discontinuities in
the radiance functions of the receiving surfaces.
This section describes the techniques used to locate these discontinuities in a
polygonal environment with no interpenetrating surfaces.
4.2.1 Locating D0 Discontin?ities
Value, or D0, discontinuities are responsible for the highest illumination contrast.
Failure to capture these discontinuities by global illumination simulations has
resulted in the shadow/light leaks and the floating objects that are evident in so
many radiosity images (see Figures 3.1 and 3.2.)
D0 discontinuities are due to touching surfaces and point light sources. Point
light sources can be handled similarly to area light sources and their implemen-
tation will not be discussed here.
When surfaces touch, one surface may act as an occluder preventing light
from the source from reaching part of the other surface, or receiver. If the oc-
cluder lies on the receivet? i.e. the two surfaces are coplanat then the D0 discon-
tinuities will occur along the perimeter of the occluder; otherwise, the occluder
may have either an edge or a vertex lying on the receiver, resulting in a D0 dis-
continuity occurring along a line segment or at a single point respectively.
The location of these discontinuities is then independent of the position of
the light sources and can be computed in a separate preprocessing step.
61
class Surface has
Plane plane
Polygon polygon
list of Segment D0edges
list of Point D0 vertices
method ComputeD0 (in Segment)
end class
class Node has
130UndingBox box
list ofNode children
list of Surface surfaces
method ComputeD0 (in Segment)
end class
Figure 4.2: Pseudocode for the computation of D0 discontinuities (Part 1)
Implementation
Given the set 5 = fSl,S2,... ,Sn) of input surfaces, we could easily find all D0
discontinuities by checking all possible pairs of surfaces, (Si, Sj), i # I. Assuming
that the number of edges per polygonal surface is bounded, this strategy would
result in 0(n2) time complexity.
A more efficient algorithm can be designed by organizing the input surfaces
into a hierarchy of bounding volumes, a popular scheme used by many ray trac-
ing systems [WHIT8O].
Pseudocode for the algorithm is given in Figures 4.2 and 4.3. An object of
type Surface is used to represent each surface Sj in the environment; the object
provides access to the polygonal representation of the surface as well as its plane
equation; also, the object is used to store a list of the D0 discontinuities falling on
62
its associated surface.
The environment is organized into a hierarchy of bounding volumes. A pop-
ular technique for constructing such a hierarchy was described by Goldsmith
and Salmon [GOLD87]. Class Node is used to represent the nodes of the hierar-
chy. Each node has its bounding box, a list of its children, and, if it is a leaf node,
a list of the surfaces it contains. The hierarchy of bounding volumes is accessed
through the special node root.
The algorithm checks each edge in the environment against the hierarchy of
bounding volumes, creating D0 segments and points as it progresses.
Given an edge e, the algorithm traverses the hierarchy of bounding volumes
visiting only the nodes whose bounding boxes contain at least part of the edge.
If a leaf node is reached, then all the surfaces it contains are checked against e.
When an edge e is checked against a surface s, two cases may arise. If the
edge lies in the supporting plane of the surface, we compute the part, if any, of e
contained within the boundaries of s and insert it into the list of D0 discontinuity
segments falling on the surface. If, instead, the edge intersects the surface, we
insert the point of intersection into the list of D0 discontinuity points falling on
the surface.
The performance of this algorithm is strictly dependent on the quality of the
hierarchy of bounding boxes used. In the degenerate case, the hierarchy reduces
to a linear list, giving a worst time complexity of 0(n2). In general, though, for a
relatively balanced hierarchy, we can expect a O(n log n) performance.
63
foreach edge e in the environment do
root. ComputeD0 (e)
end for
method Node::ComputeD0 (in Segmente) is
if e intersects boxthen
foreach Surfaces in surfaces do
s. ComputeD0 (e)
end for
foreach Nodendin children do
nd. ComputeD0 (e)
end for
end if
end method
method Surface::ComputeD0 (in Segmente) is
if e lies on planethen
Segmente' H intersection of e with polygon
ife' #NlLthen
add e' to D0edges
end if
else
Polntp & intersection of e with plane
ifp lies inside or on the boundary of polygon then
add p to D0vertices
end if
end if
end method
Figure 4.3: Pseudocode for the computation of D0 discontinuities (Part 2)
64
4.2.2 Locating D' and D2 Discontinuities
Discontinuities in the first and second derivative of the radiance function, re-
ferred to as D1 and D2 discontinuities respectively, occur along the umbra and
penumbra boundaries cast by finite area light sources as well as along other crit-
ical curves inside penumbra regions.
Conventional radiosity algorithms ignore these discontinuities (see Chap-
ter 3) and yield images with jagged shadow boundaries and missing or distorted
shadows such as those presented in Figures 3.1 and 3.2. The ability to determine
the location of these discontinuities can significantly improve the accuracy of
radiosity simulations.
Choosing which discontinuities to represent
The D1 and D2 discontinuities are generated by the visual events discussed in
Section 2.4. There are two different types of visual events: EV and FEE events.
In a polygonal environment with m edges, these events can result respectively
in 0(m3) and 0(m4) critical curves. Trying to compute all the D1 and D2 discon-
tinuities would result in time and storage costs so high as to limit the system to
anything but the simplest environments.
Not all visual events, though, result in discontinuities of the same visual im-
portance. A careful selection may result in an algorithm of lower complexity
capable of accepting most simulated environments.
As a first simplification, we decided to ignore FEE events. These events can
result in discontinuities along curved segments. Capturing these discontinuities
would require a discontinuity mesh capable of representing curved edges and
65
non-polygonal elements. Also, a reconstruction method capable of accepting
curved boundary elements would be needed. Fortunately, the boundary be-
tween penumbra and unoccluded regions is marked by ?V events only; thus,
ignoring EEE events does not affect the main shape of the shadows.
As a further simplification, we considered only the EV events to which the
light source contributes either the participating edge or participating vertex. Other
EV events are only significant if the critical surfaces they define intersect the light
source. These events are relatively few, but considering them would add an or-
der of magnitude to the cost complexity of the algorithm.
The number of remaining EV events is only 0(m) and the corresponding crit-
ical curves are straight line segments for polygonal environments. In particular,
these segments include all the boundaries separating totally unoccluded regions
from partially occluded ones, as well as additional important discontinuities in-
side the penumbrae.
The set of EV events that we considered is not as limited as it may appear.
In fact, it is a superset of the visual events considered by other shadowing and
discontinuity meshing algorithms [CAMP9ic,CHIN92,NISH83].
Classifying EV visual events
The order of the discontinuity associated with a critical curve is determined by
the rate of change in the amount of light source area that is visible from the re-
ceiver when moving across the critical curve (see Figures 2.6 and 2.7.) If the rate
of change is lineat the discontinuity is D1; otherwise, it is D2. If we consider the
visual event and critical surface that generated the critical curve, we will notice
that the above rate of change is the same everywhere across the critical surface;
66
this means that all critical curves corresponding to a given visual event share the
same order of discontinuity.
Determining the order of discontinuity associated with an EV visual event
is simple. Let V and E be respectively the vertex and edge defining the visual
event. If v is the endpoint of an edge E' that is coplanar with E, the event should
be classified as D'; otherwise, it should be classified as D2
It is convenient to distinguish between the EV events involving a source ver-
tex and the ones involving a source edge. In what follows, we will refer to the
former kind as VE events and to the latter as EV events.
Computing critical curves for VE visual events
A VE visual event captures the interaction of a vertex V of the light source and an
edge E of an occluding surface. The critical surface corresponding to the event
is a planar wedge W of infinite extent with its apex at vertex V and bounded by
rays through the endpoints of the edge E (see Figure 4.4(a).)
As shown in Section 2.4, discontinuities in the derivatives of the radiance
function happen along critical curves corresponding to the intersection of a crit-
ical surface with the surfaces in the environment. The discontinuities are caused
by abrupt changes in the visibility of the light source. Not all intersections, how-
ever, are critical curves. Only the intersections that are visible from the apex V
of W, i.e. the critical surface, are actual critical curves. If fact, if an intersection
segment is not visible from V, then no change in the visibility of the light source,
and therefore no discontinuity, may occur along the segment.
Computing the set of critical curves corresponding to a given VE event is then
similar to a visible surface determination problem in 2D (see Figure 4.4(b).) The
67
source
vertex
(a)			(b)
Figure 4.4: (a) A wedge corresponding to a VE event. (b) Interaction of the wedge
with the surfaces in the environment. The wedge is clipped against surfaces A,
B, and C, but only surfaces B and C gain a new discontinuity boundary. Discon-
tinuities are not added to surfaces intersecting the wedge in the region enclosed
between V and E, since the visual event is not visible from such surfaces. The
dark areas show the clipped portion of the wedge.
68
surfaces in the environment are sorted in order of increasing distance from V.
If a surface is intersected by W along a visible segment and it faces V, then the
segment of intersection is inserted in the surface mesh as a discontinuity edge. If,
however, the surface is closer to V than the occluding edge E is, then the surface
is not effected by the visual event and no critical curves should be generated.
Computing critical curves for EV visual events
An EV visual event captures the interaction of an edge E of the light source and
a vertex V of an occluding surface. The critical surface corresponding to this
event is slightly more complex than that of a VE event; the surface is bounded
by edge E and by semi-infinite rays originating at the endpoints of the edge and
passing through vertex V (see Figure 4.5(a).)
Computing the set of critical curves corresponding to a given EV event is
done similarly to the tracing of VE events except that the surfaces in the envi-
ronment are sorted in order of increasing distance away from E rather than V
(see Figure 4.5(b).) If a surface is intersected by the event along a segment that is
visible from E and it faces V, the the segment is inserted in the surface mesh as
a discontinuity edge. If, however, the intersected surface lies between E and V,
then no discontinuity may result.
Implementation
The pseudocode given in Figure 4.6 shows how the selected visual events dis-
cussed above are generated and classified with the appropriate order of discon-
tinuity.
Three sets of visual events are generated. First, VE events are generated from
69
source edge
occiuding
vertex			V
(a)
v
Figure 4.5: (a) A wedge corresponding to an EV event. (b) Interaction of the
wedge with the surfaces in the environment. The wedge is clipped against sur-
faces A and B, but only surface B gains a new discontinuity boundary. Disconti-
nuity boundaries are not added to surfaces intersecting the wedge in the region
enclosed between E and V, since the visual event is not visible from such surfaces.
The dark areas show the clipped portion of the wedge.
70
class VEevent has
Point V
Segment E
DiscOrder order
method Set(in Point, in Segment, in DiscOrder)
end class
class EVevent has
same as above
end class
Point C H centroid of source So
foreach edge E of source So do
VEevent ve.Set(C, E, D2)
modeL ThaceVE( ve)
end for
foreach vertex V of source So do
foreach edge E in set of potentially occiuding edges do
if 50 has an edge E' incident on V and coplanar with E then
VEevent ve.Set(V, E, D')
else
VEevent ve.Set(V, E, D2)
end if
modeL TracevE( ye)
end for
end for
foreach edge E of source 5o do
foreach vertex V in set of potentially occluding vertices do
if V is the endpoint of an edge E' that is coplanar with E then
EVevent ev.Set(E, V, D1)
else
EVevent ev.Set(E, V, D2)
end if
modeL TraceEv(ev)
end for
end for
Figure 4?6: Pseudocode for the generation and classification of visual events
71
the centroid of the light source through each of its edges. These events are used
to locate the discontinuity in the second derivative of the radiance function on
those surfaces in the environment that are intersected by the supporting plane of
the light source. In fact, the contribution of the light source is maximum in the
direction normal to the source, decreases away from the normal, reaches zero
for points coplanar to the source, and is of course zero for any point behind the
source.
Second, we trace VE events from each of the vertices of the light source. For
each vertex V, there are as many possible visual events as there are edges in
the environment. Only a subset of these edges, however, has an impact on our
simulation. The set of potentially occluding edges contains only edges that have
at least one endpoint in front of the light source. Edges that are shared between
two or more polygons are counted only once. If the endpoints of an edge lie on
opposite sides of the light source, the edge is clipped and only the part lying in
front of the source is considered.
If the environment is made of closed non-interpenetrating polyhedra, this set
can be made smaller by including only edges incident to at least one polygon fac-
ing the apex V. In fact, a visual event whose participating edge E' did not belong
to any polygon facing V would only see polygons that lie between V and E', and
would therefore be unable to generate any critical curves (see example in Fig-
ure 4.7.)
If the environment contained only convex polyhedra, the set of potentially
occluding edges could be made even smaller by including only those edges that
appear as silhouette edges when seen from V. The detection of silhouette edges,
72
Fl
1E
Figure 4.7: Visual event defined by an edge E' whose incident faces are all facing
away from V. In a polyhedral environment, this event is inconsequential since
it is completely arrested by the faces, F? and F2, intervening between V and E'.
however, was not implemented.
Determining the order of discontinuity carried by these VE events is simple.
We take the two edges of the light source, E' and E11, that are incident on v; if
either edge lies in the plane defined by V and E, then the visual event will create
D' discontinuities; otherwise, it will create D2 discontinuities.
Finally, we trace EV events from each of the edges of the light source. Simi-
larly to VE events, we compute a set of potentially occluding vertices by includ-
ing only those vertices in the environment that lie in front of the light source and
belong to a polygon facing the light source. Vertices sharing the same location
are inserted only once.
In order to determine the order of discontinuity of an EV event, each occlud-
ing vertex V stores a list of its incident edges. If any of these edges is coplanar
with the edge of the light source generating the EV event, then the event is clas-
73
sified as D1; otherwise, the event is D2.
As mentioned before, computing the critical curves associated to a VE visual
event is essentially a visible surface determination algorithm. We order to im-
plement this process efficiently, we organize the surfaces in the environment in a
BSP-tree [FUCH80]. This data structure, in fact, provides a simple and fast way
of sorting its surfaces by distance with respect to any point in space.
The pseudocode given in Figures 4.8 to 4.10 shows how to compute the crit-
ical curves (or segments) generated by a VE event ve.
The environment is stored as a BSP-tree. Each node stores a plane equation
plane and a list of the surfaces lying on that plane, as well as subtrees to the
part of the environment in front of the plane (Childout) and that behind the plane
(Childin).
Given a VE event ve, the BSP-tree is traversed in front-to-back order along
the wedge defined by the participating vertex and edge of ye, starting from its
apex ve.V (see Figure 4.9.) If the wedge faces the plane of a BSP-tree node, then
the wedge is clipped against the surfaces associated with that node (see Fig-
ure 4.4.) Only the unclipped parts of the wedge are traced further. The traversal
can stop as soon as all portions of the wedge have been clipped. In order to keep
track of which parts of a VE event remain to be traced, the data structure for the
event stores a list of active intervals. At first, this list contains a single interval,
representing the entire wedge. When tracing the event through the environ-
ment, however, portions of the wedge may be clipped and old intervals may be
dropped from the list of active intervals or replaced by new smaller ones. Should
the list of active intervals become empty, the recursive traversal of the BSP-tree
74
class VEevent has
list of Interval intervals
end class
class Surface has
Plane plane
Polygon polygon
method CipVE(inout VEevent)
end class
class BSPnooe has
Plane plane
list of Surface surfaces
BSPnode &hild??ut, Childin
method ThaceVE(inout VEevent)
method ClipVE(inout VEevent)
end class
Figure 4.8: Pseudocode for the location of the discontinuities generated by a VE
event (Part 1)
75
method BSPnode::ThaceVE(inout VEevent ye) is
if ye,interyals = NIL then
return
end if
if ve. V is in front of planethen
Childout. TraceVE( ye)
if ve faces planethen
Op'VE( ye)
childm. ThaceVE( ye)
end if
else if ye. V is behind planethen
childm. TraceVE( ye)
if ve faces planethen
CIpVE( ye)
Childout. ThaceVE( ye)
end if
else
if veE lies on planeoreither endpoint of veE is in front of planethen
Childout. TraceVE( ye)
end if
if either endpoint of ye,E is behind planethen
child1?. TracevF( ye)
end if
end if
end method
method BSPnode::ChpVE(inout VEevent ye) is
foreach Surfaces in surfaces do
s. Op'VE( ye)
if ye.intervals = NIL then
break
end if
end for
end method
Figure 4.9: Pseudocode for the location of the discontinuities generated by a VE
event (Part 2)
76
method Surace::CiipVE(inout VEevent ve) is
iistoflntervalintervals' & NIL
foreach Intervall in ye.intervals do
clip I against polygon to get Ii, 12,13
ifi1 #NiLthen
intervals' & intervals'+f11)
end if
ifplanefaces veVandi2 ? NIL and 12 is not between ve.V and veEthen
add 12 to surface's discontinuity mesh
end if
ifi3 #NiLthen
intervals' & intervals' + ?I3?
end if
end for
ye.intervals & intervals'
end method
Figure 4.10: Pseudocode for the location of the discontinuities generated by a VE
event (Part 3)
may be safely interrupted.
The pseudocode given in Figure 4.10 shows in more detail how a VE event
ve can be clipped against a surface in the environment. Each active interval I
of ve is examined in turn. First, I is intersected with the surface polygon. Since
the polygons in a BSP-tree are guaranteed to be non-concave, I can be split at
the most into three subintervals. I1 and 13, if non-empty, fall outside the polygon
and are therefore the new active intervals; 12, if it exists, lies inside the polygon
and is a candidate for insertion in the discontinuity mesh of the surface. As each
interval I is examined, a new list of active segments, intervals', is created which
will replace the old list in ye.intervals at the end of the process.
A number of speedups were used to obtain an efficient implementation. For
77
A
Rj?
:.Rout
Figure 4.11: Coplanar surfaces in a BSP-tree node are grouped in two sets, ac-
cording to their orientation (light surfaces face outwards, dark surfaces face in-
wards.) Each set is further organized in a two-level hierarchy of bounding rect-
angles.
clarity, however? these details were omitted from the pseudocode. One of these
efficiency schemes is nonetheless worth mentioning.
The coplanar surfaces stored with a BSP-tree node are grouped into two sets:
one containing all the surfaces that face towards the outside of the supporting
plane of the node and one containing all the surfaces that face towards its inside
(think of the plane as a halfspace.) Each set of surfaces is further organized in a
two-level hierarchy of bounding rectangles (see Figure 4.11.) At the lower level,
each surface is bounded by its own rectangle; at the higher level, a single rect-
78
angle bounds all the surfaces. All of these rectangles are aligned to the axis of a
two-dimensional coordinate system local to the plane of the BSP-tree node.
Before any of the active intervals of a VE event are examined, the entire wedge
representing the event is intersected with the plane of the BSP-tree node. If the
intersection segment falls completely outside a higher level bounding rectan-
gle, then its contents are skipped. Otherwise, the intersection segment is further
checked against the lower level bounding rectangles before the active intervals
of the visual event can be clipped against the surfaces.
This simple scheme is very effective due to the following reasons:
1. Large groups of coplanar surfaces may be distant from the path of the vi-
sual event and can be culled efficiently using the higher level bounding
rectangles.
2. Polygons in the BSP-tree may have many vertices and be therefore costly
to intersect. Using the lower level bounding rectangles can avoid needless
computations.
3. Visual events can be fragmented into many intervals. Checking the event
first in its entirety may quickly cull out many needless intersection tests.
The implementation of EV events is not too different from that of VE events.
An EV event evcan be looked at as two VE wedges sharing the same apex V (see
Figure 4.5.) First, the BSP-tree is traversed in front-to-back order from V in the
direction of E; then, once E is reached, a new front-to-back traversal is started
from V in the direction opposite to E. The list of active intervals of ev is modi-
fied by the surfaces intersected during the first traversal, but no critical curves
79
are generated. The same list is then used as the initial list of active intervals
for the second traversal. This time1 any visible intersections with polygons fac-
ing V will result in new discontinuity segments being added to the appropriate
discontinuity meshes.
4.3 Constructing the Discontinuity Mesh
The critical curves computed in the previous section are used to construct a dis-
continuity mesh capable of capturing all the important discontinuities in the ra-
diance function and its derivatives. Each surface maintains its own mesh. This
discontinuity mesh is represented by a data structure called a Discontinuity Mesh
Tree, or DM-tree.
A DM-tree consists of two parts: a two-dimensional BSP-tree and a data
structure of topologically interconnected faces (elements), edges, and vertices.
Each inner node of the tree contains the line equation of a discontinuity segment.
Each leaf contains a pointer to a face bounded by the intersection of the halfs-
paces encountered along the path down from the root of the tree. This data struc-
ture is constructed incrementally: each discontinuity segment is filtered down
the tree, eventually splitting only the intersected leaves. The faces attached to
these leaves are also split, and the newly added edges are labeled with the ap-
propriate order of discontinuity.
The topological mesh is arranged similarly to a winged-edge data structure
so that no duplication of edges and vertices is necessary and each element can
easily determine its neighbors. In this way the contribution of the source can be
computed once on each vertex and shared between the incident elements. Since
80
edges are shared between neighboring elements, no T-vertices are introduced
when elements are split.
Two-dimensional BSP-trees were also used to represent radiosity meshes by
Campbell and Fussell [CAMP90]. Note however, that the elements stored at the
leaves of their element BSP-tree were not topologically connected. This results in
redundant illumination computations, as well as in a need for a separate pass to
eliminate T-vertices [CAMP9lcj.
Figure 4.?2 shows an example of how a DM-tree is constructed incrementally
from a sequence of visual events. The simple configuration given in (a) shows
the discontinuity lines on a receiving surface R, generated by the interaction of
vertex V of source S with obstacle 0. Initially, (b) the DM-tree of the receiver
consists of a single node pointing to a single face F. When the first edge of the
obstacle is projected on the receivet, (c) the line equation a of the resulting dis-
continuity line is used to split face F into faces F? and F21 and the DM-tree is
updated so that its root now stores line equation a and points to the two new
faces on either side of a. Notice that two collinear edges have been created along
line a: one corresponds to the actual discontinuity segment, while the other is
due to our splitting scheme and is not marked as a discontinuity. Next, the sec-
ond VE event is processed (d). A new discontinuity segment is generated and
the DM-tree is searched from the root to determine whether any face is crossed
by this segment. Face F2 is found and is therefore split into faces F21 and F22.
Finally, the third VE event is processed to yield the DM-tree depicted in (e).
81
v
1
21 22
(d)
(a)
`
` I
121
1			2
1			b
21
2l)222
Figure 4.12: Incremental construction of a DM-tree
82
4.3.1 Implementation
The pseudocode given in Figure 4.13 shows the data structures used to imple-
ment a DM-tree.
The hierarchical portion of the DM-tree is made of DMnode nodes. If the
node is a leaf, then face points to the face in the topological mesh corresponding
to that node. Otherwise, the node partitions the surface into two halfspaces;
variable line stores the equation of the dividing line, and Childout and Childm
point to the subtrees of mesh elements that are respectively to the outside and
inside of line.
The topological part of the mesh is represented by a winged-edge data struc-
ture. Each face points to a doubly linked ring of edges. Each edge has pointers
to its two endpoints, neighboring edges, and incident faces (possibly NIL.) The
boundary of face[i] can be traversed either clockwise or counter-clockwise by
following the chain of prev[i] or next[ij pointers respectively. Each vertex has its
coordinates and a pointer back to one to its incident edges. The order of discon-
tinuity of an edge or a vertex is stored in variable order.
To start with, each surface has a discontinuity mesh consisting of a single
node and a single face. This mesh is then refined incrementally by inserting the
critical curves and their associated order of discontinuity. The D0 discontinuity
segments computed in Section 4.2.1 are inserted first. Then, D1 and D2 disconti-
nuity segments are computed as shown in Section 4.2.2 and inserted in the mesh
on the fly.
The pseudocode given in Figure 4.14 shows how a discontinuity segment is
inserted in a discontinuity mesh by filtering it down starting at the root of the
83
class Surface has
DMnode mesh
method InsertSegment(in Segment, in DiscOrder)
end class
class DMnode has
pointer to Face face
Line line
DMnode Childout, Childin
method InsertSegment(in Segment, in DiscOrder)
end class
class Face has
pointer to Edge boundary
method AdjustRoundary(in Segment, in DiscOrder)
method Sp]it(in Segment, in DiscOrder, out Face, out Face)
end class
class Edge has
DiscOrder order
pointer to Vertex endpoint[2]
pointer to Edge next[2], prev[2]
pointer to Face facd2]
end class
class Vertex has
Point coord
DiscOrder order
pointer to Edge edge
end class
Figure 4.13: Pseudocode for the DM-tree data structures
84
method Surface::Insert5egment(in Segmentseg, in DiscOrderord) is
mesh.InsertSegment(seg, ord)
end method
method DMnode::InsertSegment(in Segmentseg, in DiscOrderord) is
if this is an interior node then
if seg crosses line then
split seg across line to yield seg0?? and seg??
Childout.ThsertSegment(seg0??, ord)
child??.ThsertSegment(seg1.?, ord)
else if seglies on the outside of linethen
childout.InsertSegment(seg, ord)
else
childm.Insert5egment(seg, ord)
end if
else
if segfalls on the boundary of facethen
face.AdjustBoundary(seg, ord)
else
line H supporting line of seg
Childout create new leaf node
ChIldin H create new leaf node
face.Split(seg, ord, Childout.face, Childin face)
dispose of face and mark this node as non-leaf
end if
end if
end method
Figure 4.14: Pseudocode for the insertion of a critical curve into the discontinuity
mesh of a surface
85
DM-tree. Notice that the input segment is generated by the algorithms given in
the previous section and is guaranteed to lie completely inside (possibly on the
boundary) the surface being meshed.
If the node being visited is an interior node, we compare the discontinuity
segment against the line equation stored at that node. If the segment crosses
the line, we split it across the line and attempt to insert each of the resulting
pieces independently by recursively searching the Childout and Childout subtrees
respectively. If the segment falls completely on one side of the dividing line,
then we continue the search on that side only. If the segment lies on the dividing
line, then we continue the search down the hierarchy by arbitrarily choosing to
follow the child? subtree.
When a discontinuity segment reaches a leaf node, it must be inserted into the
topological part of the mesh. The discontinuity segment must either fall on the
boundary of the face associated with the leaf node, or be completely contained
within the face. The pseudocode given in Figure 4.15 describes these two cases
in detail.
In the former case1 the endpoints of the segment are added to the boundary
of the face possibly splitting one or two edges. The edges covered by the seg
ment are then reclassified according to the order of discontinuity of the segment.
Notice that if two discontinuities overlap, the one with the lower ordet i.e. the
sharper one, will prevail.
In the latter case, the face is split along the supporting line of the discontinuity
segment resulting in a new edge E (see parts a and b of Figure 4.16.) The edge is
initially classified as D?. Then, E is split in up to three parts at the endpoints of
86
method Face::AdjustBoundaiy(inSegmentseg,in DiscOrderord)is
foreach endpoint V of segdo
if v lies on an edge E of boundajythen
split E and create a new vertex corresponding to V
end if
end for
foreach EdgeEcovered by segdo
if ord< E.orderthen
E.order			ord
end if
end for
end method
method Face::Split(inSegmentseg,in DiscOrderord,out FaceF0??,out FaceF1?)is
split face along segto create new faces Fout and Fjn joined by edge E
E.order
clip E against segto get E1, E2, E3 (E1 and E3 may be NIL)
E2.order & ord
end method
Figure 4.15: Pseudocode for the insertion of a discontinuity segments into the
leaves of a DM-tree
E3
(a)			(b)
Figure 4.16: Splitting a face to capture a discontinuity segment. (a) The initial
face and discontinuity segment S. (b) The face is split along the supporting line
of S and a new edge E is created. (c) E is split into three parts, with edge E2
corresponding to the discontinuity segment S.
87
(a)			(b)			(c)
Figure 4.17: Triangulating a mesh element. (a) The initial element. (b) The new
vertex is added halfway between the two most distant vertices. (c) New edges
are created between each vertex on the element boundary and the new vertex.
the segment and the part that is covered by the segment is reclassified with the
order of discontinuity of the segment (Figure 4.16c.)
4.4 Sampling and Reconstruction on the Receiver
Once the discontinuity mesh is built, sample locations must be chosen at which
to compute the radiance contribution of the current light source. The contributed
radiance function is then reconstructed from the samples.
The discontinuity mesh built using the techniques described in the previous
sections partitions a surface into regions within which the radiance function is
continuous. Because of the binary space partitioning scheme used, the resulting
mesh elements are convex polygons, but may have many vertices. Reconstruct-
ing a smooth radiance function across such general elements while meeting the
boundary conditions is difficult. Therefore, each region is triangulated first. The
triangulation algorithm avoids creating new vertices on the boundaries of the
mesh elements so that each element can be triangulated independently. The al-
88
gorithm is illustrated in Figure 4.?7 by way of a simple example. We start with
a convex element with possibly collinear vertices (a). First, we choose the pair
of vertices Vj and Vj that are the farthest apart and create a new vertex halfway
between Vj and Vj (b). And then, we triangulate the element by creating new
edges between each vertex on the element boundary and the new vertex (c).
Triangular elements of different orders can now be used to represent the ra-
diance function. Constant and linear elements are only able to resolve D0 and D1
discontinuities respectively and generate annoying visual artifacts such as block-
iness and Mach bands. Higher order elements can correctly resolve D01 D11
and D2 boundary discontinuities, and can more accurately approximate the ra-
diance function inside each element.
We use quadratic elements with each mesh triangle defined by six radiance
values: one at each vertex and one at the midpoint of each edge (see Figure 4.?9.)
Rather than actually computing six values for each triangular element, radiance
values are shared between neighboring elements across non-D0 edges. Correct
treatment of D0 edges and vertices requires some attention. The radiance func-
tion is not uniquely defined along these edges and vertices and the radiance val-
ues in their vicinity depend on the side from which the edge is approached. The
correct values for a given element are the limit values of the radiance function at
the edge midpoints and vertices when moving towards them from inside the el-
ement. In practice, a good approximation to these limit values can be computed
by slightly displacing the sample locations corresponding to the vertices and the
midpoints of the edges inward, so that they lie strictly inside the element.
A special problem is posed by elements which are incident on a D0 vertex V,
89
(a)
(c)
Fb)
Figure 4.18: Radiance function reconstruction by piecewise polynomial interpo-
lation. (a) The underlying discontinuity mesh. (b) Constant elements. (c) Linear
elements. (d) Quadratic elements.
90
Figure 4.19: Quadratic interpolation over a triangular mesh element
o+--H D0
____ non- D0
Figure 4.20: A simple discontinuity mesh with D0 edges and D0 vertices marked
by thicker lines. Face F, among others, is incident on a D0 vertex, V, but has no
D0 edges. The radiance function over F has a singularity at V.
91
but do not contain D0 edges (e.g. elementF incident on V in Figure 4.20.) There is
a need to correctly represent the singularity that the radiance function may have
at V. Assigning a single radiance value to V may cause D0 discontinuities in the
interpolated function along the edges incident on it. To solve this problem we fit
a degenerate bi-quadratic patch over the triangular domain of the element. The
degenerate patch has three control points corresponding to V and thus a range
of values can be represented simultaneously there, while maintaining continuity
along the incident edges.
Our implementation correctly resolves D0 and D1 discontinuities. The actual
radiance function is C1 continuous across D2 boundaries. Currently, we do not
impose C1 continuity across these boundaries. It is possible to do so using the
piecewise cubic interpolation scheme described by Salesin et a!. [SALE92].
Certain operations, e.g. ray tracing, require the ability to quickly inquire for
the radiance value at a given location x on a surface s?. The pseudocode given
in Figure 4.21 shows how this radiance value can be efficiently evaluated using
the DM-tree data structure representing the discontinuity mesh for Si.
First, the DM-tree is descended recursively, starting at the root, until the leaf
node, or element, containing point x is located. Then, point x is transformed
to the (u, v) parametric coordinates of this element. Finally, the quadratic patch
representing the radiance function over the element is evaluated at (u, v) and the
desired radiance value Lj(x) is stored in variable L.
Assuming that the DM-tree is relatively balanced, this operation takes O(log n),
where n is the number of elements in the discontinuity mesh of a single surface,
St.
92
method DMnode::EvaluateRadiance(in Pointx, RadianceL) is
if this is a leaf then
Realu, v & parametric coordinates of x in this face
L & radiance value of quadratic patch at (u, v)
else
if x is on the out side of linethen
Childout.EvaluateRadiance(x, L)
else
Childm'.EvaluateRadiance(x, L)
end if
end if
end method
Figure 4.21: Pseudocode for the evaluation of the radiance function at a point x
on a surface
4.5 Computing the Source Contribution to a Point
The radiance contributed by a source ?o of constant emission L0 to a point x is
given by Equation (4.1).
If no occluding surfaces intervene between x and So, the visibility term v is
everywhere 1, and the area integral can be transformed into a countour integral
using Stoke's theorem. This, in turn, can be evaluated analytically using the
formula provided by Hottel and Sarofin [HOTT67j (see Figure 4.22a):
1 n--H1 cos--H?(R?Rj?i)
RjeiII			`			(4.2)
FdA(X)?AO=TirZjo			Rj? 1(R?xRj?i)?Ni
where
FdAj(x)?A0 is the form factor notation for the integral in Equation (4.1);
represents addition modulo n;
n is the number of vertices of the polygonal source so;
93
ffl\.Rj?i
RiD).
5
-Q
(a)
(c)
(b)
Figure 4.22: Analytical form factors. (a) Hottel and Sarofin's method for com-
puting unocciuded form factors. (b) Only the part of the source lying in front of
the receiver is considered in computing the energy transfer. (c) Only the visible
parts of the source contribute to the form factor.
94
Rj is the unit vector from point x on surface Sj to the ith vertex of source So;and
Nj is the unit normal to surface St.at point x.
As Nishita and Nakamae first pointed out [NISH83j, this formula requires
that the light source So be completely in front of the receiving surface 51 If, in-
stead, the supporting plane of surface Sj intersects the source So,then the source
must be split and the integration must be carried out on that part of Sothat is in
front of Sj only (see Figure 4.22b.)
The contribution of source So to a point x can be computed analytically in
the presence of occlusions, too. In fact, it is sufficient to determine which parts
of So are visible from x, and integrate the contribution of each part separately.
To compute these visible regions, we project all the potential obstacles from the
point of view ofx onto the supporting plane of the source 5o and discard any parts
of the source that are covered by any of the projected obstacles (see Figure 4.22c.)
Point sampling algorithms for testing the visibility of the light source [COHE85,
HANR9lb,WALL89] can also be used, but they are not nearly as accurate, and
are prone to aliasing.
4.5.1 Implementation
Figure 4.23 gives the pseudocode for the computation of the form factor from a
point x on a surface 5 to a polygonal light source So.
Each occluding surface in the environment is projected onto the supporting
plane of the source So and its projection is intersected with the source to deter-
mine which parts of Soare visible from the receiving point x.
Considering all the surfaces in the environment as potential occluders would
95
compute shadow shaft for x and So
foreach Surface s? in the environment do
ifs? is not culled by the shadow shaft then
project s? onto supporting plane of So
intersect projection with BSP-tree for So
end if
end for
Real ff			0
foreach visible leaf polygon p in BSP tree for So do
ff & ff+ form factor FdAj(x)?Ap
end for
Figure 4.23: Pseudocode for the computation of form factors in the presence of
occlusions
be terribly inefficient. Instead, we use a shaft culling technique similar to that
proposed by Haines and Wallace [HAIN9?c] to drastically reduce the number of
potential occluders that must be projected onto the source So
Our shadow shaft is a pyramid with the source So as its base and the receiving
point x as its apex. More formally, the shaft is defined as the intersection of
halfspaces: one halfspace for each of the edges of the source So, defined as the
plane passing through the receiving point x and the edge itself, and one halfspace
for the supporting plane of the source itself.
Checking a surface sj against our shaft is quick and easy: if all the vertices ofs?
lie outside any of the halfspace defining the shaft, then s? can be safely culled out
of consideration as its projection onto the supporting plane of So could never
intersect the source.
In order to keep track of which parts of the source are visible from x we use
a two-dimensional BSP-tree. A boolean variable is associated with each of the
96
leaves of the BSP-tree to indicate whether or not the region of the source corre-
sponding to that leaf is visible or occiuded. Each time the projection of a new
surface &j is intersected with So, the BSP-tree is updated to reflect the change in
the visible portions of the source. The intersection algorithm uses the techniques
presented by Naylor et al. [NAYL9O], who show how polyhedral set operations
can be implemented by merging BSP-trees.
Once the visible parts of the source are determined, the form factor from x
to So, FdAj(x)?+A0? is computed as the sum of the unocciuded form factors from x
and each of the visible parts of the source. Each of these form factors, in turn, is
computed using the analytical formula reported in Equation (4.2).
4.6 Adaptive Subdivision
Discontinuity meshing guarantees that no mesh elements will be crossed by
strong discontinuities, thus making simple polynomial functions suitable to ap-
proximate the radiance distribution across a receiving surface. However, there
are still instances in which the radiance across a mesh element varies too quickly
for the simple quadratic interpolation scheme presented in Section 4.4 to provide
a sufficient approximation. Large or badly shaped mesh elements and the close
proximity of a strong light source are often responsible for this problem. In these
cases, adaptive mesh subdivision can be used to effectively reduce the errors.
In order to design a mesh refinement algorithm, we must decide when and
how to subdivide a mesh element.
97
4.6.1 when to Subdivide
Intuitively, a mesh element should be subdivided whenever its approximation
to the correct radiance function falls below the desired accuracy according to a
predetermined error metric.
Defining the best error metric to use is difficult. Each metric quantifies a
particular kind of error. In the following, we will present three different error
metrics and discuss their effect on the accuracy of the solution as well as their
ease of implementation and computational cost.
2-norm error metric
The 2-norm error over a surface element e is defined as:
where
IlL --H Li12 =			--H L(x))2dA(x))?
L(x) is the actual radiance at point x on element e;
(4.3)
L(x) is the estimated radiance at point x on element e as computed from the
quadratic interpolant defined over e;
dA(x) is a differential surface area centered at x.
This metric measures the overall error over the entire surface e. Using the
2-norm error metric to control adaptive subdivision has the effect of improving
the quantitative accuracy of the solution. However the distribution of the error
across the surface is not taken into account; for example, the same error could be
localized to a small subregion or be equally disftibuted across the whole surface
98
element Visually, the former case would cause much more noticeable artifacts
than the latter; thus, the 2-norm error metric is not very well suited for qualitative
analysis. Furthermore, the error measured by this metric is proportional to the
area of the surface element and can therefore lead to unnecessary subdivision of
large elements.
Infinity-norm error metric
The infinity-norm error over a surface element e is defined as the maximum ab-
solute error over the element:
IL--HII? =max, IL(x)--HL(x)I
(4.4)
Using the oc-norm error metric to trigger the refinement process has the ef-
fect of improving the accuracy of the solution qualitatively, i.e. visually, as well
as quantitatively. Also, unlike the 2-norm error metric, this metric is indepen-
dent of the element size and thus allows the simulation to capitalize on large
elements whenever possible. This independence, though, coupled with numer-
ical approximations, could lead to infinite subdivision. This problem can be eas-
ily solved by preventing subdivision whenever the element area falls below a
prespecified minimum.
Computing the exact error according to either the 2-norm or the oo-norm er-
ror metric is not feasible since it would require knowledge of the actual radiance
at infinitely many points across a surface. An approximation, though, can be
computed using a finite number of samples. The more the samples, the better
the approximation; however, the cost of estimating the error should, of course,
be less than that of computing the radiance samples necessary to define the in-
99
terpolant used to approximate the radiance function over an element.
On the average, excluding boundary conditions, samples taken on a mesh
vertex will typically be shared among about six elements and those taken on a
mesh edge will be shared between two elements; thus, using a piecewise quadratic
interpolant to reconstruct the radiance function will take about two samples per
element (3 x 61 + 3 x 21). This suggests that a single radiance sample be used to
estimate the error and decide whether or not to trigger adaptive subdivision. We
take this sample at the centroid of the element. Although the error estimate thus
computed is not always reliable and can cause the adaptive subdivision to stop
prematurely; in practice, this strategy has yielded good results.
Smooth-boundary error metric
The piecewise quadratic interpolation scheme presented in Section 4.4 does not
impose C' continuity across elements boundaries. This may lead to unwanted
Mach bands along mesh edges classified as D2 or higher. Adaptive subdivision
can help eliminate these artifacts.
In fact, we can measure the angle of the cusp formed by the reconstructed
radiance function along the edge between two elements, and subdivide if this
angle is too far from flat. More formally, let angle(f, e) be the function that takes
as input a triangular face f and an edge e incident on f, and whose behavior is
described by the following pseudocode:
x H midpoint(e)
g & neighbor(f,e)
if order(e) > D2 and g # NIL then
angle			arccos(tangent(f, e, x) tangent(g, e, x))
else
angle
loo
endif
return angle
where
inidpoint(e) is the midpoint of edge e;
neighbor(f, e) is the neighbor of face f across from edge e;
tangent(f, e, x) is the tangent vector to the reconstructed radiance function for
element f, at point XI and in the direction perpendicular to edge e.
Then, we define the smooth-boundary error of a triangular face f as:
max angle(f, e?) --H
(4.5)
where E = (e?) is the set of edges incident on
This error metric has the effect of forcing subdivision in areas where the piece-
wise quadratic interpolant fails to reconstruct a smooth approximation to the
radiance function.
4.6.2 HowtoSubdivide
The mesh subdivision strategy we have adopted uses the 2-triangles mesh re-
finement algorithm presented by Rivara (RIVA87].
The meshes produced by her method have a number of desirable properties
that make them particularly suitable to the task of capturing radiance functions:
1. No T-vertices are introduced by the refinement process. This property is
important since T-vertices could introduce value discontinuities in the re-
constructed radiance function.
?o1
2. The transition between large and small triangles is smooth; i.e. it happens
gradually over intermediate-sized triangles rather than abruptly. The re-
suiting distribution of sample points is less prone to generating visual ar-
tifacts at reconstruction time.
3. The refined mesh tends to be made of well shaped triangles. More pre-
dsely, the algorithm yields triangles whose longest edge is shorter than
twice the length of either of the other edges [RIVA87]. This property is
important because the binary space partitioning scheme used to construct
the discontinuity mesh is liable to produce elements with bad aspect ratios.
Well shaped elements1 in fact, can reduce errors in the radiance approxima-
tion [BAUM9lj.
2-triangles mesh refinement algorithm
Rivara's 2-triangles mesh refinement algorithm is based on simple bisection op-
erations and is illustrated in Figure 4.24. The algorithm starts by bisecting an
elements by the midpoint of its longest edge e (a). This operation creates a new
vertex v. If element s has a neighbor t across edge e, then v is a T-vertex and the
triangulation is no longer considered valid. In order to restore validity, element
t is also bisected along the midpoint of its longest edge e, (b). If e and e, are the
same edge, then we are done; otherwise, v and v', the midpoint of edge e,, are
connected to form new triangles (c). If v' is, in turn, a T-vertex, the process is
iterated until all mesh elements are triangles once again.
?O2
(a)
Figure 4.24: Rivara's 2-triangles mesh refinement algorithm
Implementation details
In our implementation, whenever a triangular element is bisected, the corre-
sponding dm-tree node is refined as well. In this way, point location within a
mesh can be carried out efficiently using the same binary search described in
Section 4.4. Also, maintaining a one-to-one correspondence between the mesh
elements and the leaf nodes of the dm-tree, enables us to implement adaptive
subdivision without any change to the existing data structures.
The new edges created by the refinement process are marked as D?. This
information can be used at radiance reconstruction time to maintain the highest
possible order of continuity across these edges. In the current implementation,
howevet the reconstructed radiance function is only guaranteed to be ? across
these edges.
Every time a new triangular element is created during the refinement process,
103
new radiance samples must be taken at its vertices and edge midpoints. Most
of the sample locations, however, will fall on the vertices and edge midpoints of
the elements of the original mesh. Recomputation of the same values is avoided
by means of a simple flagging scheme.
This adaptive refinement algorithm has proved to work very well. In partic-
ular, even a small amount of subdivision can improve the quality of the discon-
tinuity mesh considerably.
Chapter 5
Discontinuity Meshing Radiosity
In the previous chapter, we introduced the concept of discontinuity meshing. We
applied this concept to the design of an algorithm capable of accurately captur-
ing the characteristics of the radiance functions over a set of polygonal surfaces
resulting from the direct illumination of a finite area light source of constant
emission.
Despite its accuracy, this algorithm, as it is, is not very useful. Real environ-
ments typically contain multiple light sources and often a significant portion of
their illumination is due to interreflections between surfaces. In this chapter, we
show how to develop a radiosity algorithm that uses the single source discon-
tinuity meshing algorithm as a building block to compute a global illumination
solution for diffuse polygonal environments.
5.1 Discontinuity Meshing Radiosity--HAlgorithm Overview
The global illumination problem was discussed in Section 2.2 and its formulation
for diffuse environments was given in Equation (2.3). Following the conventions
established before, we restate our problem as follows:
104
105
Given a set S = fsi, ...... ,Sn J of convex polygonal surfaces, their dif-
fuse BRDF'sfri,fr2, . `fm, and constant emissions U1, L2e,. . ,L?n,c0m?
pute the set of radiance functions L1, ..... . ,Ln such that
Lj(x) = L?t +fit			cos Oi cos
? !x'Es Lj(x')			r2			v(x, x')dA(x')			(5.1)
sjES
foralli=I,2,...,n.
As with the single source discontinuity meshing algorithm presented in the
previous chapter, we will use a meshing approach to sample the radiance func-
tions at a finite set of locations and reconstruct an approximation L1 to the actual
solution.
The new discontinuity meshing radiosity algorithm is based on progressive
refinement radiosity [COHE88b]. This approach essentially reduces the prob-
lem of solving the diffuse rendering equation into that of computing the direct
illumination from a single light source. Since we have already developed an
algorithm that solves this last problem, the progressive refinement approach is
only a natural choice.
Pseudocode for the discontinuity meshing radiosity algorithm is given in Fig-
ure 5.1. The cumulative effect of multiple light sources is maintained in L1. Ini-
tially, the radiance Li of a surface Si is simply equal to the emission Uj (possibly
zero) of the surface. Then, the approximation Lt is refined by adding, at each
iteration, the direct illumination Lij contributed by a new light source sj, whether
primary or secondary.
In order to implement this approach, two major steps are necessary. First,
the single source discontinuity meshing algorithm must be extended to handle
secondary sources whose emission is typically not a constant. Second1 a method
106
foreach Surt??ces? in S do
Li & L?j
?Lj & Uj
end for
repeat until convergence
find surface &` with maximum unshot power
compute Lij for all Sj in S
foreach Sur(aces? in S do
L1 & Li+Lij
ALj & AL + Lij
end for
ALj & 0
end repeat
Figure 5.1: Pseudocode for the discontinuity meshing radiosity algorithm
of combining radiance functions implicitly defined by sample locations on dif-
ferent discontinuity meshes must be devised.
In the following sections, we describe these two steps in detail and conclude
with a discussion on the algorithmic efficiency.
5.2 Locating the Discontinuities
As shown in Section 2.4, the direct illumination due to an area light source of
constant emission results in radiance functions that may exhibit D0, D1, and D2
discontinuities. The ability to predict the order and location of these discontinu-
ities was used in the previous chapter to design an algorithm capable of accu-
rately capturing the direct illumination due to such a source on a set of polygonal
surfaces.
When the indirect illumination is also considered, additional discontinuities
are created. According to the Discontinuity Propagation Law presented in Sec-
107
tion 2.4.3, in fact, discontinuities of order k on a secondary light source are prop-
agated into the environment, resulting in discontinuities of order k +1 and k +2.
Thus, the total number of discontinuity curves in the environment grows rapidly
with the number of interreflections considered. In the limit, a full global illumi-
nation solution would exhibit an infinite number of discontinuities.
Trying to locate all these discontinuities would be impractical even for a finite
number of interreflections. Fortunately, though, as the order of a discontinuity
increases, its effect on both the quantitative and qualitative accuracy of our sim-
ulations diminishes.? Following this simple observation, it is possible to drasti-
cally cut down on the number of discontinuities considered without noticeably
affecting the quality of the simulations.
5.2.1 Implementation
In our implementation, we consider D0, D1, and D2 discontinuities only. Captur-
ing higher order discontinuities would add an order of magnitude to the time
and storage costs of the simulation and would require modifying our reconstruc-
tion scheme. The piecewise quadratic interpolation technique presented in Sec-
tion 4.4, in fact, is capable of resolving D0, D1, and D2 discontinuities only.
In order to propagate discontinuities on a secondary source that might result
in D1 or D2 discontinuities on a receivet we consider visual events in which the
participating edge is either a D0 or D1 discontinuity segment on the source.
The pseudocode given in Figure 5.2 shows how to generate and classify the
1The order of a discontinuity is actually only one of several factors whose interplay deter-
mines the visual importance of a discontinuity. Developing a metric of visual importance in-
volves the study of complicated perceptual issues and was not attempted in this research. A
brief discussion of the problem, howevet is offered in the concluding chapter.
108
foreach discontinuity segment E on source &` do
if order of E = D0 then
foreach vertex V in set of potentially occiuding vertices do
if V is the endpoint of an edge E' that is coplanar with E then
EVevent ev.Set(E, V, D1)
else
EVevent ev.Set(E, V, D2)
end if
model. ThaceEV(ev)
end for
else if order of E = D1 then
foreach vertex V in set of potentially occiuding vertices do
if V is the endpoint of an edge E' that is coplanar with E then
EVevent ev.Set(E, V, D2)
modeL ThaceEV(ev)
end if
end for
end if
end for
Figure 5.2: Pseudocode for the generation and classification of the visual events
associated with the discontinuity curves on secondary light sources
visual events associated with the discontinuity curves of a secondary light source
Additional visual events are generated by the edges and vertices of the light
source as it was done for primary light sources (see pseudocode in Figure 4.6.)
We look at each discontinuity segment on source &j, in turn. D0 discontinu-
ities interact with the set of potentially occluding vertices (see discussion in Sec-
tion 4.2.2) to generate D1 and D2 EV visual events which are then traced through
the environment. D1 discontinuities1 instead, can potentially generate both D2
and D3 EV visual events, but only D2 events are actually created and traced.
D2 discontinuity segments on the source, finally, can only propagate to D3 and
?O9
higher order discontinuities and are therefore ignored.
5.3 Computing the Source Contribution to a Point
In a global illumination simulation, any surface in the environment may act as
a secondary light source by reflecting part of the incident light back into the
environment. These sources, as opposed to the primary light sources discussed
in the previous chapter, typically have non-uniform emittance.
The radiance contributed by a source to a point x on a receiving surface 5j
is given by:
cos O? cos			(5.2)
Lij(x) = fri x'Es1 Lj(x')			r2			v(x, x')dA(x')
As was shown in Section 4.5, the visibility term v(x, x') can be taken out from
under the integral by determining which parts of sj are visible from point x.
However, the presence of the arbitrary Lj(x') term still precludes the possibility
for an analytical solution.
Our algorithm uses numerical integration techniques to compute the contri-
bution from each of the visible parts of the light source. These regions are convex
polygons of arbitrarily many sides. Rather than trying to integrate L1 over these
complicated domains, we first triangulate each region and compute the radiance
contribution as:
where
O?cos0jdA(x')
Lij(x) = fri ? tt !x'Et Lj(x')dA(x') ixIEt cos r2
tET
T is the set of visible triangles for &j; and
At is the area of triangle t.
(5.3)
110
Notice that the radiance contribution of each individual triangle is computed
as the product of the area average radiance of the triangle and its form factor.
Splitting the contribution integral in this fashion is only accurate if the radiance
across a triangle is relatively constant. In order to meet this condition, our al-
gorithm adaptively subdivides the triangle into four similar triangles until the
relative error between two successive estimates falls below a specified threshold.
The fx'Et Lj(x')dA(x') term in the equation above is the radiant flux of source S?j
over the area of triangle t. Assuming that the behavior of Lj over t can be well ap-
proximated by a quadratic polynomial, it is possible to evaluate the integral an-
alytically using the quadrature formula for quadratic triangular elements given
by Brebbia and Dominguez [BREB89], page 269:
Lj(x')dA(x') = ?31(Lj(x1) + Lj(x2) + Lj(x3)) (5.4)
where Xl, X2, and x3 are the midpoints of the edges of t.
The form factor integral is also evaluated analytically using Equation (4.2) as
was shown in Section 4.5.
5.4 Combining Radiance Functions
Conventional radiosity systems use a single mesh throughout the entire simula-
tion. The mesh can be adaptively subdivided in an attempt to capture fine illumi-
nation details, but basically the same set of sample locations is used to compute
the contribution from all the light sources, both primary and secondary, as well
as to store the accumulated radiance values.
Such a scheme would be impractical for discontinuity meshing, since every
surface in the environment may become a source, participating in a different
?i1
set of visual events. Accounting for all these visual events would result in a
very large number of critical curves and very complex meshes. Computing the
radiance contribution from every source to such a mesh would be very costly.
Our solution, instead, consists of building a separate discontinuity mesh for
each source &j using only the critical curves resulting from visual events involv-
ing this source. In this way, the number of points at which the contribution Lij
from a source S?j to a receiving surface Sj is computed is substantially reduced.
The cumulative radiance function Lj is also stored as a discontinuity mesh
since, of course, any discontinuity in a Lij will also be a discontinuity in Li. Our
progressive radiosity system refines Li incrementally by adding in at each itera-
tion the latest contribution L4:
Ln?ew & Lo1ld + L4			(5.5)
In this way, only the storage used to represent the current Li needs to be main-
tained, while that required by the L4 is released at the end of each iteration.
Since for each source &j the discontinuity mesh of the interpolant Lij is typi-
cally substantially different from the meshes used for previous sources, the prob-
lem of merging contributions is not trivial.
In order to compute Ln?ew we start by creating a new discontinuity mesh con-
taining all the significant discontinuity boundaries from LOjid and L4. These bound-
aries can be identified by inspecting the corresponding discontinuity meshes and
can be inserted into the new discontinuity mesh using the algorithm outlined in
Figure 4.14. The new mesh is then triangulated and quadratic elements are com-
puted as described in subsections 4.4 and 4.6. Note however, that the value at
any point p is given by LOjld(p) + L4(p). Thus, rather than integrating over the
?12
source at each sample location of the new mesh, we only need to evaluate LOjld(p)
and 4(p). Each value is computed by locating the element containing p in the
appropriate mesh and evaluating the quadratic interpolant. The key steps of the
algorithm are depicted by way of a simple example in Figure 5.3.
Since there is no need to process visual events nor to compute the contribu-
tion of a source at each location, constructing the new interpolant Ln?ew is typically
much faster than computing the contribution Lij. As will be shown in Chapter 6,
in fact, the accumulation process takes up only a small percentage of the total
computation time and thus does not compromise the advantages gained by us-
ing a separate discontinuity mesh for each source.
Our scheme for combining radiance functions is essentially a resampling
method, and therefore a certain loss of accuracy can be expected. Notice, though,
that since the important critical curves are reproduced in the new discontinuity
mesh, the resampling occurs inside regions where the radiance function can be
considered smooth, and thus the additional error is small. Furthermore, since
the time spent combining radiance functions constitutes a very small part of the
total simulation time, we can afford to increase the level of adaptive subdivision
used by the combination algorithm so as to provide the desired level of accuracy.
5.5 Efficiency Considerations
The discontinuity meshing radiosity algorithm presented in this chapter is very
accurate, but also very computationally intensive.
Despite the restrictions on the kind of visual events that we chose to trace
through the environment (see Sections 4.2.2 and 5.2,) the process of locating the
113
(a)			(b)			(c)			(d)
Figure 5.3: Combining the illumination from two light sources in the presence
of a simple obstacle (the cube floating above the floor). Columns (b), (c), and (d)
show a plan view of the floor. Row 1. (a) Surface A is selected as the current
source. (b) Discontinuity boundaries due to source A are computed. D? and D2
discontinuities are drawn respectively in red and yellow. (c) A triangular mesh
which resolves the discontinuity boundaries is produced. (d) The source con-
tribution Lfloor?A(X) is computed at sample locations within the mesh to yield the
interpolant L?oor,A. Row 2. The process is repeated with surface B selected as
the source. Row 3. (b) The discontinuity boundaries due to A and B are com-
bined. (c) A new triangular mesh with the combined discontinuity boundaries
is produced. (d) Interpolants L?or,A and Looor,B are evaluated at sample locations
within the new mesh and the sample values are added to yield the combined
contribution of A and B.
114
radiance function discontinuities takes a large portion of the execution time. The
visual events generated by secondary light sources are mainly responsible for
this cost. The number of visual events associated with a secondary source, in
fact, may be very large due to both the possible large size of the surface compared
to that of typical primary sources and to the possible large number of D0 and D1
discontinuity segments falling on it (for example, the floor of a furnished room
which contains many shadows.)
Computing the radiance contribution of the light source to a point is also a
computationally intensive process. Determining the visible portions of a source
for each sample point we want to take, in fact, can be much more expensive
than using a more conventional point sampling approach such as the hemicube
[COHE85j or ray-traced form factors [WALL89].
Another possible source of inefficiency can be found in our solution to the
problem of combining radiance functions. Any time that two radiance functions
are merged, all their discontinuities are inserted in the resulting discontinuity
mesh. Since the set of discontinuities generated on a receiver is typically different
for every source, an attempt to compute an accurate global illumination solution
for even a medium size environment might result in storage problems.
Many of the above inefficiencies can be alleviated by noticing that many com-
putations are carried out to a much higher degree of precision than is required
to achieve a final radiosity solution of a given accuracy. The same observation
led Pat Hanrahan et al. [HANR91bj to the design of a rapid hierarchical radiosity
algorithm, which is an order of complexity faster than its predecessors.
In the following sections, we use this observation to improve the efficiency
115
of the discontinuity meshing radiosity algorithm.
5.5.1 Computation of Radiance Discontinuities
In Chapter 4 we showed how incorporating lines of discontinuity in the radiance
function (including its first and second derivatives) as edges in a surface mesh
may result in an illumination algorithm capable of eliminating shadow leaks,
jagged shadow boundaries, and other artifacts that afflict conventional radiosity
systems.
Some of these discontinuity edges, howevet may sometimes be omitted from
a surface mesh without resulting in any visible artifacts. For example, if a light
source s? contributes but a very small percentage to the overall illumination of
a receiving surface Sj, any discontinuity lines Sj may cause on 5 are likely to be
unnoticeable.
Limiting the range of visual events
The radiance contribution of a source &j to a surface Si decreases quadratically
with the distance between Sj and Sj. When the two surfaces are sufficiently dis-
tant, then, computing discontinuity lines becomes unnecessary. In our disconti-
nuity meshing radiosity algorithm, discontinuity lines are computed by tracing
visual events from a source through the environment, visiting the receiving sur-
faces in sorted order, away from the source. This computation may be shortened
by terminating the traversal as soon as we reach a surface Sj such that Lij falls be-
low a specified threshold. The precise value of Lij is not critical for termination
criteria as long as we can compute an upper bound for it. Such an upper bound
can be computed much faster than Lij itself by simply assuming total visibility
?16
between the source and receiver.
Dropping weak discontinuities when combining radiance functions
In Section 5.4, we described an algorithm for combining radiance functions that
takes two piecewise quadratic interpolants, LOjid and Lij, and yields a new inter-
polant, Ln?ew, representing their combined effect. Each of the initial interpolants
is defined over its own discontinuity mesh. The discontinuity mesh for n?ew is
constructed starting with a single DM-tree node and adding in all the disconti-
nuity segments from both L191d and Lij.
As we observed in the introduction to this section, this algorithm can result in
an excessive growth of the DM-tree data structures as the progressive refinement
radiosity progresses. This growth can be controlled, often without significant
losses in accuracy, by noticing that not all discontinuities are equally noticeable.
Some of the discontinuities of L4, in fact, may become invisible once they are
combined with LOjid, and vice versa. Developing a metric capable of assessing
the visual impact of a discontinuity in the illumination across a surface, how-
ever, involves complicated perceptual issues that were beyond the scope of the
research described in this thesis. A brief discussion of the problem is attempted
in the concluding chapter.
Although we recognized the need for perceptual metrics, only a simple, purely
quantitative, approach was actually implemented. Rather than trying to mea-
sure the relative strength of a discontinuity segment c, we measure its relative
contribution to the radiance of the receiving surface; assuming c is coming from
?17
a radiance function Lij, we compute:
max			L4(x)			(5.6)
xEc Lo1ld(x)+Lq(x)
and insert the segment in the discontinuity mesh for Ln?ew only if its strength ex-
ceeds a specified threshold eL. Only a few samples taken along the discontinuity
segment c are needed to obtain a good estimate of its relative contribution. By
controlling the magnitude of CL we can control the growth of our data structures
and are therefore able to apply our radiosity algorithm to more complex envi-
ronments.
5.5.2 Computation of the Contribution of a Source to a Point
In Section 5.3, we described an algorithm for computing the contribution from a
source sj of arbitrary emission to a point x on a receiving surface Sin the presence
of occiusions.
Despite the efficient implementation using shaft culling (see Section 4.5), a
significant part of the computation time is still spent projecting obstacles onto
the supporting plane of the source and intersecting them with s? to compute the
parts of the light source that are visible from each sample point x.
When the contribution Lq(x) is very small, its calculation by means of the
above algorithm is carried out to unnecessary accuracy. Thus, whenever the
contribution is estimated to fall below a specified threshold, we can evaluate
Ljj(x) by means of a faster ray-traced form factor algorithm [WALL89] without a
significant loss in the overall accuracy of our solution Li.
Once again, a precise estimate of the contribution L4(x) is not necessary as
long as we can compute an upper bound for it. This upper bound can be com-
118
puted very quickly by numerically integrating Lij(x) over the entire source S?j US
ing the algorithm presented in Section 5.3, but omitting the visibility calcula-
tions.
When the contribution from &j to a receiving surface s? is computed using
different algorithms for nearby sample locations, annoying visual artifacts may
result. These artifacts are due to a sudden change in the accuracy of the contri-
bution estimate when switching from the analytical projection-based method to
the ray-traced form factor algorithm or vice versa. This problem can be avoided
by requiring that the same method be used for all the sample points Xk on the sur-
face Si necessary to capture the contribution Lij of the source &j. In other words,
a single method is used to compute the energy transfer within a given source-
receiver pair.
If each input polygon is treated as an independent surface, radiance disconti-
nuities may occur at the boundary between two contiguous surfaces, resulting in
the visual artifacts mentioned above. Our implementation avoids this problem
by preprocessing the input scene so that each surface is actually a connected set
of coplanar polygons of homogeneous surface reflectance characteristics. Fur-
thermore, polygons within the same surface are topologically connected into a
single surface discontinuity mesh. This way, mesh conformity is maintained
across neighboring polygons even during adaptive subdivision of the mesh.
Chapter 6
Results
6.1 Physical vs. Simulated World
This section presents a qualitative comparison between a photograph of a simple
physical environment and a computer generated rendering of the same scene
based on discontinuity meshing radiosity.
The test environment consists of three surfaces only: a receiver, a light source,
and an occiuding surface. The receiver is a square piece of wood covered with
several coats of matte white paint. The source is a lamp centered above the re-
ceiver and casting its light through a square-shaped diffusing panel. The obsta-
cle is a triangle cut out of black cardboard and suspended over the receiver using
fish line.
The receiver was surrounded on three sides (one side was missing to allow
for photographs) by matte black panels so as to reduce interreflections and allow
us to compare the effect of direct illumination only.
Figure 6.1 shows the reference photograph of this simple scene. Figure 6.2
shows a rendering of the same view, computed from a view-independent so-
119
120
Figure 6.1: Photograph of a simple scene
Figure 6.2: Computer simulation of a simple scene
121
lution generated using the single source discontinuity meshing algorithm de-
scribed in Chapter 4. The two images are remarkably similar. Even subtle illu-
mination details are captured. Note, for example, the dark bands extending from
the corners of the umbra and into the penumbra region as well as the bright hor-
izontal Mach band along the outer edge of the penumbra near the bottom of the
picture.
The behavior of the illumination across the receiving surface can be perhaps
better understood by looking at a cross section of the radiance function. Fig-
ure 6.3 shows a cross section taken along a horizontal scanline passing through
the middle of the triangular umbra of Figures 6.1 and 6.2. Only an overall quali-
tative comparison is meant here since the views in the original images are aligned
only very approximately. The physical values reflect noise introduced by the
photographic and scanning processes. The scanline shown in the figure tra-
verses both the umbra and penumbra regions cast on the receiver. The plot of the
simulated radiance clearly shows how shadow boundaries correspond to sharp
changes in an otherwise smooth radiance function.
These sharp changes in the illumination occur along discontinuity bound-
aries. Figure 6.4a shows their location in a top view of the receiving surface.
Figure 6.4b, instead, shows how these lines can be incorporated into a disconti-
nuity mesh to yield the solution displayed in Figure 6.2.
6.2 Single Source Discontinuity Meshing vs. Analytical
This section presents a comparison between a discontinuity meshing solution
and an analytical solution of a simple environment.
122
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
physical
simulated
0
200			400			600
pixel
800			1000			1200
Figure 6.3: Plot of the radiance values across a scanline of the images shown in
Figures 6.1 and 6.2. The scanline crosses the middle of the triangular umbra in
the original images. Only an overall qualitative comparison is meant here since
the views in the original images are aligned only very approximately. The phys-
ical values reflect noise introduced by the photographic and scanning processes.
?23
Unfortunately no closed-form solution to the radiosity equation is known
for arbitrary environments. It is possible, however, to analytically compute the
contribution of one or more polygonal light sources of constant emission to any
point in an arbitrary polygonal environment (see Section 4.5.)
Therefore, we chose to compute our reference solution in a view-dependent
fashion, by sampling the illumination at a finite set of sample points arranged
in a regular grid on a projection plane, and compare it to a rendered view of a
discontinuity meshing solution computed using direct illumination only.
Two views were selected for the comparisons presented in this section. The
first, shown in Figure 6.5, is a view of a simple environment illuminated by a
single square light source placed at the center of the ceiling. The second, shown
in Figure 6.6, is a top view of the floor of the same simple environment.
Because of the view-dependent method, it was necessary to compute a dif-
ferent reference solution for each of the two figures shown in the upper right
corners (b). The discontinuity meshing solution, instead, was computed only
once and then rendered twice, to obtain the two different views shown in the
upper left corners (a).
The resulting mesh for the discontinuity meshing solution is shown in im-
age (c) of the two figures. Image (d), instead, shows a pseudocolored image of
the absolute error between discontinuity meshing and reference radiance val-
ues, computed on a pixel-by-pixel basis. Finally, image (e) shows a histogram of
the pseudocolored error image. The histogram is colored to show the mapping
between error values and colors.
The pseudocolored images are very helpful for a qualitative analysis of the
124
error. A look at the mesh and error images shown side-by-side in Figure 6.5, re-
veals a correlation between the magnitude of the error and the underlying mesh
structure. In general, in fact, the error is small along mesh edges and is larger in
the middle of mesh elements. Furthermore, inspection of the mesh and error im-
ages in Figure 6.6 shows that the largest errors, corresponding to the red regions,
occur near singular points. In fact, the largest error found occurs around a sin-
gularity where the lower left corner (as seen from above) of the tall box touches
the floor (see Figure 6.6d.) The magnitude of this error is 0.099 Wsr-1 m?2, while
the brightest radiance value in the reference images was 4.72 Wsr-1m-2.
Figure 6.8a shows a plot of the radiance function along the scanline passing
through such largest error. The match between the discontinuity meshing and
analytical solutions is so close that it was necessary to zoom in in order to reveal
any discrepancies (Figure 6.8b.)
In order to give a quantitative measure of the error, we used 2-norms [COLU89j.
Given an image A of size n-by-m, we define its 2-norm, IA112, as:
All2 = n'm A?1 s?1(A?? (6.1)
We can then define the absolute error E and the relative error e between two
images A and B respectively as:
and
E=IIA--HBII2= ?;mZ?i$i(A4?B4?
(6.2)
e= llAllBJIB2ll2			?z)?f;?ffl(A)1$)i?4)2			(6.3)
This definition of absolute and relative errors is useful because it gives a mea-
sure of the average error over the entire images while at the same time giving
125
particular weight to large localized errors. In contrast, computing the overall
error as a simple unweighted average of the individual pixel errors would have
had the effect of washing out localized, but possibly large, errors. Note also that
the Aij'5 and B4,s are actual radiance values, i.e. values which have not yet been
mapped to colors for display on a CRT device.
Table 6.1 reports the absolute error E and relative error e for both the full
view shown in Figure 6.5 and the top view of the floor shown in Figure 6.6.
The 2-norm of the reference image, 11B112, is also given for reference. As it can
be seen, the discontinuity meshing algorithm produced a solution accurate to
within about half of one percent erron
6.3
Accumulating Radiance Functions by Resampling vs. Exact
Addition
This section examines the error introduced by the resampling scheme used to
combine radiance function due to different sources.
The accumulation algorithm was described in detail in Section 5.4. In order to
add two radiance functions their underlying discontinuity meshes are inspected
and a new discontinuity mesh is created. The new mesh contains as edges the
discontinuity segments of each of the two input meshes, but can otherwise be
very different from them. Radiance values for each of the new mesh elements are
then determined by sampling the two input radiance functions at corresponding
locations. This mechanism essentially introduces an additional resampling step
in our radiosity algorithm and is likely to introduce additional errors in our ra-
diance approximations. Because the discontinuity segments are reproduced in
the combined mesh and because the radiance function generally varies slowly
126
(a)
(b)
Figure 6.4: Top view of the floor of the simple scene. (a) The critical curves result-
ing from the interaction of the light source and the triangular obstacle. (b) The
resulting discontinuity mesh.
11B112 (Wsr--H'm			 (Wsr--H'm--H2)			e
full view			0.54676			0.00326			0.00597			I
single surface			0.47820			0.00215			0.00450
Table 6.1: Discontinuity meshing vs. analytical solution. Error measurements
for the comparisons shown in Figures 6.5 and 6.6.
127
0
(c)
O?OO625
0.0125
(e)
0.025
0.05			0.1
sr?1 ?
W			m
Figure 6.5: Full view. Discontinuity meshing vs. analytical solution. (a) The dis-
continuity meshing solution. (b) The pixel-by-pixel analytical solution. (c) The
mesh underlying the discontinuity meshing solution. (d) Pseudocolored image
of the absolute error between discontinuity meshing and analytical solutions.
(e) Histogram of the absolute error image. For better resolution, error values are
mapped to colors on a logarithmic scale. The maximum and average brightness
of the reference image are respectively 4.7219 and 0.4015 Wsr-1 m?2
128
0
(c)			(d)
0.00625
0.0125
(e)
0.025
0.05
W
0.1
m?2
Figure 6.6: Floor surface only. Discontinuity meshing vs. analytical solution.
(a) The discontinuity meshing solution. (b) The pixel-by-pixel analytical solu-
tion. (c) The mesh underlying the discontinuity meshing solution. (d) Pseu-
docolored image of the absolute error between discontinuity meshing and an-
alytical solutions. (e) Histogram of the absolute error image. For better res-
olution, error values are mapped to colors on a logarithmic scale. The maxi-
mum and average brightness of the reference image are respectively 0.8809 and
0.3620 Wsr?1m-2
129
0.8
0.6
0.4
0.2
disc. mesh.
analytical-
0
0.5
0A
0.3
0.2
0.1
0
0 50 100 150 200 250 300 350 400
pixel
(a)
disc. mesh.
analytical-
40			70
45			50			55			60			65
pixel
(b)
Figure 6.7: Plot of the radiance function corresponding to the discontinuity
meshing and analytical solutions along the scanline containing the largest er-
ror for the images in Figure 6.5. (a) The entire scanline. (b) A blow-up of the
boxed region, showing the largest error.
C\I
130
0.9
0.8			disc. mesh.
0.7			analytical
0.6
0.5
0.4
0.3
0.2 :
0.1
0 50 100 150 200 250 300 350 400
pixel
(a)
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
disc. mesh.
0.1 analytical
0
140			150			160			170			180			190 200			210
pixel
(b)
Figure 6.8: Plot of the radiance function corresponding to the discontinuity
meshing and analytical solutions along the scanline containing the largest er-
ror for the images in Figure 6.6. (a) The entire scanline. (b) A blow-up of the
boxed region, showing a radiance singularity at about pixel 200.
13?
within each of our mesh elements, we would expect our accumulation algorithm
to only introduce minor errors.
In order to test the accuracy of the accumulation algorithm, we compare it
against a simple, but generally impractical, reference method that combines ra-
diance functions in screen space to completely eliminate any resampling errors.
The reference method stores the contribution of each shooting source in a differ-
ent discontinuity mesh and keeps them all in memory, never trying to combine
them into a single mesh. To evaluate the radiance at a given point x on a surface,
all these meshes are interrogated and the returned radiance values are summed
to yield the accumulated result. To render a radiosity solution for a given view,
we repeat this process for each pixel in our image.
Clearly, this method is impractical. The amount of memory required to store
all these meshes is prohibitive and the time needed to evaluate the radiance func-
tion increases linearly with the number of shooting sources processed. However,
this method does not introduce any resampling errors and can therefore be used
to produce reference solutions to test the radiance accumulation algorithm de-
scribed in Section 5.4.
Figure 6.9 shows a comparison of the two methods for a simple environment.
Images (a) and (b) show renderings computed respectively from the discontinu-
ity meshing radiosity and reference solutions after twenty progressive refine-
ment iterations.
The resulting mesh for the discontinuity meshing solution is shown in im-
age (c). Image (d), instead, shows a pseudocolored image of the absolute differ-
ence between the discontinuity meshing and reference solutions, computed on a
132
0
(c)
0.0125
0.025
0.05
0.1
(e)			Wsr? n?2
0.2
Figure 6.9: Accumulating radiance functions by resampling vs. exact addition.
(a) Solution computed by the discontinuity meshing radiosity algorithm. (b) Ref-
erence solution computed by adding in each contribution at each pixel. (c) The
final mesh for the discontinuity meshing radiosity solution. (d) Pseudocolored
image of the absolute difference between discontinuity meshing radiosity and
reference solutions. (e) Histogram of the absolute difference image. For better
resolution, difference values are mapped to colors on a logarithmic scale. The
maximum and average brightness of the reference image are respectively 50.420
and 0.658 Wsr-1rn-2
133
pixel-by-pixel basis. Finally, image (e) shows a histogram of the pseudocolored
difference image. The histogram is colored to show the mapping between dif-
ference values and colors. The largest difference found was 0.2295 Wsr-1m-21
while the brightest radiance value in the reference image was 50.42 Wsr-1m-2.
The overall absolute and relative differences for the view shown were com-
puted using Equations (6.2) and (6.3) to yield values of 0.00831 Wsr-1m-2 and
0.00454 respectively. For reference, the 2-norm of the reference image was
1.83035 Wsr-1m-2. Thus, after twenty iterations the relative difference is less
then half of one percent.
Figure 6.10 gives graphs showing the absolute error E and relative error e
introduced by the combination method as a function of progressive refinement
iterations. These errors were computed by comparing the radiosity solutions
generated by the two combination methods for the same number of iterations.
As it could be expected, both the absolute and relative errors grow with the num-
ber of iterations; however, both their magnitude and the rate of growth are very
small.
The error introduced by the combination method, of course, is just one com-
ponent of the overall error resulting from a progressive refinement radiosity ap-
proach. This latter error still decreases rapidly with the number of progressive
refinement iterations. Figure 6.11 gives graphs showing the absolute and relative
errors resulting from the combined effect of the progressive refinement approach
and radiance accumulation algorithm. These errors were computed by compar-
ing the intermediate radiosity solutions generated by the discontinuity meshing
radiosity algorithm after each iteration with the final reference solution shown
134
c\I
0
0
0.009
0.008
0.007
0.006
0.005
0.004
0.003
0.002
2			4			6			8			10 12 14 16 18 20
number of iterations
(a)
0.005
0.0045
0.004
0.0035
0.003
0.0025
0.002
0.0015
0.001
2			4			6			8			10 12 14 16 18 20
number of iterations
Figure 6.10: Error introduced by the radiance accumulation algorithm as a func-
tion of progressive refinement iterations. The absolute error E (a) and the relative
error e (b) were computed from a sequence of discontinuity meshing radiosity
and reference images (the final pair is shown in Figure 6.9) using Equations (6.2)
and (6.3).
135
in Figure 6.9b.
Figure 6.12 shows the discontinuity meshes for the first six iterations of our
simulation. Each discontinuity mesh corresponds to a different source. The ra-
diance contribution corresponding to each source is also shown in the figure.
In order to make the contributions visible, a different mapping of radiance
values to RCB displayed colors was used than that shown in Figure 6.9. How-
ever, the same mapping was used consistently for all of the six images shown.
The sources can be easily inferred from the pictures. Going down the rows,
left to right: the primary source in the middle of the ceiling (not shown), the back
wall, the floor, the left wall, the right wall, and the top of the right box. Each
mesh is significantly different from the others. Each one contains a different set
of discontinuity segments. Also, adaptive mesh subdivision was triggered in
different areas for each of the discontinuity meshes; in each of the pictures the
mesh is denser in those areas were the rate of change of the radiance contribution
is higher.
Could combining such different meshes result in an excessive growth of the
accumulated mesh? Figure 6.13 shows that this is not generally the case. The
graph plots the number of elements in the accumulated mesh as a function of
progressive refinement iterations. For reference, the size of the discontinuity
mesh capturing the contribution of the current shooting source is also shown for
each iteration. As can be observed from the graph, the size of the accumulated
mesh grows rapidly during the first iterations, but then it tends to level out and
reach an approximate steady state.
136
0.25
0.2
0.15
0.1
0.05
c\I
0
0
w
0.12
0.1
0.08
0.06
0.04
0.02
0			2			4			6			8 10 12 14 16 18 20
number of iterations
(a)
u
0			2			4			6			8 10 12 14 16 18 20
number of iterations
(b)
Figure 6.11: Combined effect of the progressive refinement approach and the
radiance accumulation algorithm as a function of progressive refinement itera-
tions. The absolute error E (a) and the relative error e (b) were computed from a
sequence of discontinuity meshing radiosity images (the last image is shown in
Figure 6.9a) and a reference solution (shown in Figure 6.9b) using Equations (6.2)
and (6.3).
137
Figure 6.12: Discontinuity meshes and individual contributions for the first six
iterations of the discontinuity meshing radiosity solution shown in Figure 6.9.
Going down the rows, left to right, the sources are: the primary source in the
middle of the ceiling (not shown), the back wall, the floor, the left wall, the right
wall, and the top of the right box.
138
1400
w
w
w
1200
1000
800
600
400
200
accumulated contrib.
single contrib.
(1
0			2			4			6			8 10 12 14 16 18 20
number of iterations
Figure 6.13: Growth in mesh size as a function of progressive refinement iter-
ations for the discontinuity meshing solution shown in Figure 6.9. One curve
shows the size of the mesh representing a single source contribution; the other
shows the size of the mesh representing the accumulated effect of multiple con-
tributions.
139
6.4 Discontinuity Meshing Radiosity vs. Conventional Radiosity
This section presents a qualitative comparison between discontinuity meshing
radiosity (DMR) and conventional progressive refinement radiosity (PRR).
A rendering system based on this latter algorithm was implemented during
the past several years at the Program of Computer Graphics [TRUM91]. It in-
cludes ray-traced form factors [WALL89j and adaptive mesh subdivision. Users
have control over parameters such as initial patch size, minimum element size
and the tolerance used for adaptive subdivision. Meshing can be controlled both
globally and on a surface-by-surface basis.
The model used for the comparison is a simple reproduction of a Victorian-
style room consisting of 1382 polygons. The illumination originates from two
light sources: the sun, approximated by a polygon placed at a considerable dis-
tance, and an area light source at the center of the room's ceiling. Both systems
were allowed to perform two iterations, thus the only illumination present in
the solutions is the direct illumination from the two light sources. Allowing sec-
ondary sources to shoot would have complicated the comparison significantly
and would have required software changes to the PRR system. In fact, secondary
sources correspond to input polygons in DMR and to smaller mesh patches in
PRR; thus, if a secondary source A were processed in DMR, we would have to
make sure that all the patches corresponding to surface A in PRR were processed
as well. The solution produced by DMR was computed without any need for
user intervention. Only an error threshold for adaptive subdivision needed to
be specified. The solution resulted in 35,540 elements (more statistics are given
in Tables 6.2 and 6.3.) Then, a solution with a similar number of elements (37,846)
140
(a)
Figure 6.14: Victorian Room. A ray-traced view of a discontinuity meshing ra-
diosity solution (a) and the corresponding conventional progressive refinement
radiosity solution (b).
141
Figure 6.15: A closer look at the wall. A ray-traced view of a discontinuity mesh-
ing radiosity solution (a) and the corresponding conventional progressive refine-
ment radiosity solution (b).
142
Thb)
Figure 6.16: A closer look at the chair. A ray-traced view of a discontinuity mesh-
ing radiosity solution (a) and the corresponding conventional progressive refine-
ment radiosity solution (b).
143
was produced by the PRR system. This solution was obtained through a process
of trial and error: we experimented with various initial meshing densities and
minimum element sizes to obtain reasonable quality. To make a fair comparison
no surface-by-surface tweaking was performed; only global meshing parame-
ters were used. For accuracy of the illumination computations, 64 rays per form
factor were used by the PRR system.
Three images, each with a different view of the model, were rendered from
each of the two solutions. The rendering algorithm used was ray-tracing (eye
rays only) with adaptive supersampling for anti-aliasing.
A general view of the room is shown in the images in Figure 6.14. Note
how DMR captures the thin shadows from the window on the left wall, many of
which are entirely missed by PRR. The PRR solution also exhibits some shadow
leaks from under the window frames.
The next view (Figure 6.15) is a closer look at the left wall of the room. Note
the high quality of the shadows in the image produced from the DMR solution
and the realistic softening of the shadows on the wall as we move away from the
window. On the other hand, the artifacts in the PRR solution are now even more
noticeable. The captured shadows are severely aliased: some appear jagged;
others are too wide and blurry.
Finally, Figure 6.16 features a closer view of the chair. Note the high quality of
the shadow cast by the chair on the floor in image (a). The shadow is sharp and
crisp near the legs of the chair, and it becomes softer farther away. In image (b)
this shadow exhibits artifacts that were not visible in the general view of the
room. Also note that the shadow cast by the armrest on the seat is captured by
144
Victorian Room
Input polygons			1,382
Patches			1,382
BSP tree nodes			3,412
BSP tree polygons			5,406
D edges			2,702
VEevents			12,636
EVevents			10,532
D discontinuities			1,349
D discontinuities			10,278
Final elements			35,540
Table 6.2: Statistics for the discontinuity meshing solution of the Victorian Room
model
F			Victorian Room
Task			time			T
Build discont. mesh			11 : 52			36.81
Illumination comp.			17; 10			53.26
Combining rad. funcs.			02: 08			6.62
Miscellaneous			01: 04			3.31
Total			32:14			100.00
Table 6.3: Timings for the discontinuity meshing solution of the Victorian Room
model. All timings are given in minutes and seconds on an HP720 workstation.
145
DMR, but is entirely missing in the PRR solution. The fine self shadowing on
the chair and the window frame is correctly captured by the new method, yet
almost completely absent from the PRR solution. Another important difference
is visible in the penumbra cast by the table top on the floor. The intensity changes
smoothly in (a), while in (b) it appears to change in several steps. This is due
to the use of ray-traced form factors in the PRR system, which computes the
visibility of the source from each point of the floor by point sampling, in contrast
to the accurate hidden surface removal algorithm employed by the new method.
The sets offigures clearly demonstrate the superiority of the new algorithm in
computing solutions of high visual quality, independent of the viewer position.
6.5 Statistics
This section examines the performance of the discontinuity meshing radiosity
algorithm presented in Chapter 5 when applied to environments of "reasonable"
complexity.
Four different architectural interiors were modeled and fed to the discontinu-
ity meshing radiosity system to generate diffuse global illumination solutions.
Ray-traced renderings of these solutions are shown in Figures 6.17, 6.18, 6.20,
and 6.21 along with their corresponding discontinuity meshes.
Table 6.4 gives some statistics generated during the radiosity simulations for
the four models. Reported are the size of the input, the size of the BSP tree used
to sort the surfaces in the environment, the number of visual events generated,
the number of resulting discontinuity segments, and the size of final mesh.
Two important observations can be made from the results shown in the ta-
146
(`7)
Figure 6.17: White Room. (a) A ray-traced view of a discontinuity meshing ra-
diosity solution. (b) The corresponding discontinuity mesh.
147
Figure 6.18: Green Room. A ray-traced view of a discontinuity meshing radios-
ity solution.
148
Figure 6.19: Green Room. The final discontinuity mesh.
149
Figure 6.20: Pink Room. (a) A ray-traced view of a discontinuity meshing ra-
diosity solution. (b) The corresponding discontinuity mesh.
15)
Figure 6.21: Mackintosh Room. (a) A ray-traced view of a discontinuity meshing
radiosity solution. (b) The corresponding discontinuity mesh.
ble. First, the ratio of BSP tree nodes and polygons to input polygons is very
high. In fact, the algorithm used to build the BSP tree is very simplistic. Any
improvements capable of reducing these ratios would immediately improve the
performance of the entire discontinuity meshing radiosity system since it would
reduce the number of BSP polygons that must be intersected with visual events
when the latter are traced through the environment.
Second, although in principle there could be 0(n) discontinuity segments for
each visual event, in practice there seem to be many, many fewer. This is due
to the fact that in any realistic environment many patches are not visible to each
other.
Table 6.5 gives timings for the four test environments. The discontinuity
meshing radiosity algorithm is broken down into its constituent modules and
computation times are reported for each module.
Most of the time is taken up by the illumination computations, i.e. by comput-
ing the light energy contribution from a source to a point on a receiving surface.
This is due to the compute-intensive visibility calculations used at that step. The
cost of the illumination computations could be reduced using the techniques
discussed in Section 5.5.2 at the cost of some accuracy. The results reported in
Table 6.5, however, were computed using the more expensive approach.
Building the discontinuity mesh for each of the shooting sources amounts
to approximately one quarter to one third of the total computation time. The
higher percentage taken in the case of the Mackintosh Room environment is due
to the particularly large ratio of BSP nodes and BSP polygons to input polygons
(see Table 6.4.) A better, more balanced, BSP tree could substantially reduce the
152
WhiteRoomIGreenModels
Room Pink Room Mackintosh
Inputpolygons			1,311			1,688			1,688			1,154
Patches			1,311			1,688			1,688			1,154
BSP tree nodes			1,760			3,885			5,004			8,197
BSP tree polygons			3,103			5,522			7,303			9,500
D edges			2,561			3,092			3,092			1,879
VEevents			16,834			15,238			15,840			9,992
EVevents			13,952			13,456			13,604			8,816
D discontinuities			6,707			1,482			0			517
discontinuities			7,211			14,528			8,952			22,166
Final elements			69,139			55,348			59,947			48,187
Table 6.4: Statistics for the discontinuity meshing solution of the White Room,
Green Room, Pink Room, and Mackintosh Room models
Models
White			Room			Green			Room
Task			time			time
Build discont. mesh			07: 20			24.33			13 : 50			32.46
Illumination comp.			16: 09			53.60			24 : 30			57.49
Combining rad. funcs.			OS : 20			17.70			02 : 50			6.65
Miscellaneous			01: 19			4.37			01: 27			3.40
Total			30: 08			100.00			42 : 37			100.00
Models
PinkRoomMackintoshI
TasktimeI?timeI?I
Build discont. mesh			10 : 43			29.55			14 : 48			46.74
Illumination comp.			22: 26			61.86			12 : 50			40.53
Combining rad. funcs.			01: 49			5.01			02 : 54			9.16
Miscellaneous			01:18			3.58			01: 08			3.57
Total			36:16			100.00			31: 40			100.00
Table 6.5: Timings for the discontinuity meshing solution of the White Room,
Green Room, Pink Room, and Mackintosh Room models. All timings are given
in minutes and seconds on an HP720 workstation.
153
computation time for this module.
Finally, the algorithm used to combine the individual radiance contribution
into a single, accumulated, radiance function takes a smaller, but yet non-negligible
portion of the total computation time.
Chapter 7
Summary and Conclusion
7.1 Summary and Main Contributions
In the ten years since their introduction to the computer graphics community
[GORA84], radiosity methods have produced some of the most realistic-looking
computer generated images to date. Furthermore, radiosity is not merely an
image synthesis technique: it is a physically-based simulation of light energy
transfer that can be used to predict the value of actual physically measurable
quantities.
The ability to simulate global illumination, as well as the possibility of gener-
ating interactive "walk-throughs" from the view-independent solutions, makes
radiosity methods attractive for a wide variety of applications, ranging from ar-
chitectural design and lighting analysis to remote sensing and flight simulators.
Despite all the areas that could benefit from radiosity, this method has not
yet met with widespread acceptance. This is not surprising since with current
radiosity systems users have to go through a very tedious and time consum-
ing iterative process in order to generate a surface mesh capable of producing a
154
155
global illumination solution to the desired level of detail. Without very careful
user intervention, these meshes are prone to generating a wide range of annoy-
ing visual artifacts. In fact, as many researchers have pointed out, poor meshes
are actually responsible for most of the visual artifacts appearing in radiosity
images.
Clearly, new automatic meshing techniques are needed to generate high qual-
ity radiosity solutions with the degree of reliability necessary for practical ap-
plications. This thesis has investigated the roots of the problems with current
meshing techniques and proposed a new approach to meshing that is completely
automatic and that, used within a progressive refinement radiosity system, has
generated radiosity solutions of previously unattained quality.
In order to gain new insight into the problem of generating highly accurate
radiosity solutions and to befter understand the flaws and limitations of the
meshing techniques used in conventional radiosity systems, we have taken a
step back from the traditional finite element formulation of radiosity [GORA84]
(see Equation (2.5).) Lifting the assumption that radiosity is constant across each
element results in the more general diffuse rendering equation (Equation (2.3)),
a formulation of radiosity as a system of integral equations in the style of the
rendering equation first proposed by Kajiya [KAJI86j. Each equation describes
the radiance function over a surface Sj as the sum of simpler radiance functions,
each one representing the contribution of a single surface in the environment
(see Equation (2.7).)
An analysis of the radiance functions described by this last equation has
shown how umbra and penumbra boundaries as well as other sharp changes
156
in illumination within the penumbra regions actually correspond to discontinu-
ities in the radiance function and its derivatives. The shape, location, and order
of these discontinuities is determined by the interaction of the light sources and
obstacles in the environment. The results of this analysis have led to the concept
of discontinuity meshing, whereby shadow boundaries and other sharp varia-
tions in the illumination across a surface are captured by explicitly representing
the corresponding discontinuities as boundaries in the mesh.
Based on this concept, we have designed and implemented a discontinu-
ity meshing algorithm capable of producing accurate solutions to the radios-
ity problem for three-dimensional polyhedral environments illuminated by area
light sources. The algorithm stores a radiance function as a set of sample points
organized in a discontinuity meshing tree. This data structure allows discontinu-
ity segments to be inserted into the mesh incrementally while maintaining mesh
conformity and provides fast access to the radiance values. Piecewise quadratic
interpolation is used in conjunction with discontinuity meshing to reconstruct a
smooth radiance function while preserving discontinuities where appropriate.
Finally, in Chapter 6, we examined the accuracy of the new algorithm by
comparing it against a photograph of a physical environment, an analytical so-
lution, and a conventional, yet state-of-the-art, radiosity system, and measured
its performance on architectural models of medium complexity.
The results have been extremely promising. The new discontinuity mesh-
ing algorithm improves upon previous techniques in many respects. The major
advantages are described below.
157
Mesh generation
As opposed to conventional radiosity meshing algorithms that require specify-
ing an initial sampling density the new discontinuity meshing algorithm is corn-
pletely automatic. The resulting mesh is guaranteed to capture all the important
shadow boundaries as well as other sharp illumination details resulting from
discontinuities in the penumbra areas.
Identifying the location of these trouble areas allows the discontinuity mesh-
ing algorithm to distribute sample points very effectively In fact, a small set
of sample points suffices to capture fine illumination details that would have
required a much larger set with more conventional radiosity methods. Further-
more, while the visual quality of shadow boundaries in conventional radiosity
solutions is tightly dependent on the projected size of the mesh elements onto
the image plane, discontinuity meshing solutions exhibit a very high degree of
view-independence.
While conventional radiosity systems maintain a single mesh (possibly re-
fining it by adaptive subdivision) throughout the entire simulation, out new al-
gorithm creates a fresh discontinuity mesh at each progressive refinement iter-
ation. Within the progressive refinement radiosity framework, this strategy has
shown to require less radiance samples for a given accuracy than the more tra-
ditional approach. The possible loss of accuracy introduced by the resampling
technique used to accumulate radiance functions has proven to be negligible. As
we reported in Section 6.2, in fact, the errors introduced by this scheme are very
small (less than half of one percent in our test.)
158
Radiance reconstruction
The simple piecewise quadratic interpolation scheme used to reconstruct the ra-
diance function has proved to be superior to the more popular bilinear interpola-
tion and Couraud shading techniques. For a given visual accuracy, the quadratic
formulation requires less elements than the bilinear approach and in most cases
eliminates unwanted Mach banding without incurring the computational costs
and complexity of more sophisticated methods [SALE92j.
General radiance distributions
Removing the assumption that elements must have constant radiosity from the
radiosity formulation allows for the correct treatment of light sources of non-
uniform radiance distribution (generally the case for secondary sources.) This
eliminates any restriction on patch size resulting in more accurate and reliable
physical simulations, better efficiency, and reduced need for manual interven-
tion.
7.2 Further Research
Hopefully, the algorithm presented in this thesis and the results generated should
convince the reader of the power of the discontinuity meshing approach to ra-
diosity image synthesis. While the concept and ideas are fundamentally correct,
however, some implementation issues deserve more attention.
7.2.1 Robustness
Tracing visual events through the environment and building the discontinuity
mesh requires performing many intersection and clipping operations. Since all
159
the operations are carried out in finite precision, numerical accuracy becomes an
issue. Different visual events, for example, may sometimes result in nearly coin-
cident discontinuity segments; whether these segments are classified as parallel,
coincident, or intersecting, may make the difference between an accurate simu-
lation and one that contains streaks of light and other annoying visual artifacts.
Small errors in the arithmetic computations not only result in numerical er-
rors, but can propagate to cause inconsistencies in our data structures. The topo-
logical data structures used to represent the discontinuity meshes are particu-
larly vulnerable to this kind of errors. Special care has to be taken to maintain
consistency at all times if frustrating program crashes are to be avoided.
In our implementation, we alleviated these problems significantly by veri-
fying and enforcing consistency after numerically sensitive calculations and by
using "representatives" and Epsilon-Geometry [SALE9lj.
However, despite these precautions, the current implementation of the dis-
continuity meshing algorithm is not completely reliable. Occasionally, small in-
consistencies occur that result in visual artifacts in the displayed radiosity solu-
tion. Clearly, a more robust solution is needed before the algoritlim can be used
in practical applications.
7.2.2 Data Structures
The discontinuity meshing radiosity system we have implemented makes exten-
sive use of BSP trees. A 3D BSP tree is used to sort the input polygons so that
visual events can be traced through the environment very efficiently. A 2D BSP
tree is used as part of the DM tree data structure to represent the discontinu-
ity mesh on a surface, support incremental updates of the mesh, and allow for
160
efficient point location.
These were just two of several possible implementation choices and their se-
lection is not an inherent characteristic of discontinuity meshing algorithms. In
retrospect, we would recommend the choice of a 3D BSP tree for tracing visual
events, but would suggest investigating data structures other than the 2D BSF
tree for the surface mesh. Use of a 2D BSP tree in this context can often lead
to elements with very bad aspect ratios. Our adaptive subdivision scheme can
improve these aspect ratios when applied to these elements. However, the re-
finement step is driven by energy considerations rather than aspect ratios, and
offers no guarantee that all badly shaped elements will be corrected.
The quality of the discontinuity meshes is also affected by the order in which
discontinuity segments are inserted in the 2D BSP tree. Rather than inserting
these segments on the fly as described in Section 4.3, it might be better to collect
them all and then choose an insertion order based on some heuristics, e.g. trying
to split a node so as yields almost equal areas [LISC93].
The discontinuity mesh, however, could be represented by an altogether dif-
ferent data structure. Heckbert [HECK92a] uses constrained Delauney triangu-
lation in his 3D implementation of discontinuity meshing. This triangulation
technique is also used by Lischinski, Tampieri, and Greenberg [LISC93] to yield
remarkably well-shaped discontinuity meshes.
7.2.3 Performance and Flexibility
The new discontinuity meshing algorithm is remarkably accurate when comput-
ing the direct illumination from a small number of primary light sources. Com-
puting the illumination due to secondary light sources, however, has proven
161
to be particularly expensive. The secondary sources are typically larger than
primary light sources and their radiance functions may have several disconti-
nuities; as a consequence, many more visual events are associated with a sec-
ondary light source than are used for a primary light source. This can result
in finer meshes and can require higher numbers of radiance samples than are
needed for primary light sources. Clearly, these computations are carried out
to an unnecessary accuracy since secondary sources are typically much weaker
than primary light sources. Furthermore, a significant effort is spent propagating
the discontinuities that cross secondary light sources; since this process increases
the order of discontinuity, however, their effect on the receiving surfaces is often
hardly noticeable.
The techniques discussed in Section 5.5 were developed to speed up the com-
putation of secondary light source contributions at the cost of some accuracy. Al-
though these techniques have proved useful in computing radiosity solutions of
moderately complex environments, the benefits accrued by trading off accuracy
for speed are limited. Many energy transfers are still computed with unnecessar-
ily high accuracy, thus preventing use of the system with more complex models.
The flexibility of the discontinuity meshing algorithm is further limited by
the progressive refinement radiosity framework in which it was implemented.
Combining discontinuity meshing with hierarchical radiosity [LISC93] has re-
cently shown how discontinuity meshing techniques can be successfully applied
to the solution of radiosity problems of reasonable complexity, yielding results
of unprecedented quality while at the same time improving the performance of
the original hierarchical radiosity approach [HANR9ibJ.
162
7.2.4 Solution Accuracy and Perceptual Issues
As mentioned in the preceding section, some of the computations associated
with the discontinuity meshing algorithm are carried out to unnecessarily high
accuracy. Ideally, if the goal of the radiosity simulation is to produce photoreal-
istic results, no computing time and memory should be spent refining a solution
unless doing so resulted in visible improvements in the generated images.
Predicting whether a certain operation, e.g. propagating a discontinuity seg-
ment or subdividing a mesh element, will make a visual difference is very diffi-
cult since it would necessarily involve perceptual issues that are hard to quantify.
Our ability to perceive a radiance discontinuity as a Mach band, for example, de-
pends not only on the order of the discontinuity, but also on the brightness of the
radiance function in the neighborhood of the discontinuity, its rate of change, its
color, and the position of the observer.
Developing a perceptual metric capable of assessing the visual impact of a
given illumination detail, e.g. a shadow or a radiance discontinuity, however,
could potentially result in very significant computational savings. In fact, any
rendering algorithm that allows for trade-offs between computational resources
and solution accuracy would greatly benefit from such a metric.
7.2.5 Extensions
The discontinuity meshing radiosity algorithm presented in this thesis is limited
to ideal diffuse reflectance distributions and polyhedral environments. It is only
natural, then, to think of ways to extending it to general reflectance distributions
and curved geometries.
163
General reflectance distributions
The conventional progressive refinement radiosity algorithm was initially ex-
tended to integrate specular and diffuse reflection [WALL87,SILL89] and later
modified to accept general reflectance distributions [SILL9i]. Augmenting these
algorithms with discontinuity meshing can be done using the same techniques
presented in this thesis since the location and order of the radiance disconti-
nuities are generally not affected by the shape of the reflectance distributions.
Specular spikes in the bidirectional reflectance distribution, though, may result
in discontinuities on a receiving surface that would have to be accounted for
using different techniques.
Caustics, i.e. specular to diffuse energy transfers, also need special attention.
In order to simulate the transfer of light from a light source to a receiving surface
through an intermediate specularly reflective surface, radiosity algorithms use
the concept of "virtual source." A virtual source is created by reflecting the origi-
nal source about the mirror surface, so that the indirect transfer from the original
source to the receiving surface can be computed as the direct transfer from the
original source, through the mirror, to the receiving surface. The resulting radi-
ance contribution may exhibit discontinuities similar to the more conventional
diffuse to diffuse transfers. These discontinuities are computed by tracing visual
events involving the virtual source and treating the mirror surface as a hole, or
portal [TELL92j, into the environment.
If the reflecting surface is not a perfect mirror, however, the boundaries of
the causfic on the receiving surface will be blurry. In this case, computing the
location of the resulting discontinuities using analytical methods may not be
164
possible; in fact1 it is not even clear whether the caustic causes discontinuities
and, if it does, what their order should be.
Curved geometries
Extending discontinuity meshing algorithms to curved geometries appears to
be a difficult problem. Although the concept of discontinuity meshing is not
limited to planar geometries, the algorithms we used to trace visual events and
locate critical curves, and the data structures we used to represent discontinuity
meshes are only applicable to polyhedral environments.
Designing a discontinuity meshing algorithm for curved surfaces requires so-
lutions to non-trivial problems. Curved surfaces interact in very complex ways
to generate visual events. Computing critical curves from these visual events
requires intersecting curved surfaces. Representing such general critical curves
as edges in a discontinuity mesh requires describing mesh elements with curved
boundaries, and even if parametric representations were used, these boundaries
would be curved in the parametric space of the surface. Furthermore, recon-
structing a radiance function from such discontinuity meshes requires satisfying
continuity constraints across curved element boundaries.
Many of these problems could be addressed by discretizing the curves into
line segments. How to decide on the appropriate level of subdivision, though,
is not clear. One of the achievements of the discontinuity meshing algorithm, in
fact, has been to increase the degree of view-independence of radiosity solutions.
Approximating a curve with straight line segments would certainly reduce the
view-independence of the solutions. Future useful research could investigate
what families of curved surfaces would allow for analytical solutions to some of
?65
these problems. Some work in this direction has already been done in the field
of computer vision to determine intervisibility between curved objects [KRIE89,
KRIE9O,SRIP89J.
7.3 Conclusion
This thesis has proposed a new concept, discontinuity meshing, as a way of over-
coming the limitations of current meshing techniques for radiosity. The ratio-
nale, and mathematical justification, for this idea was provided by an analysis
of the characteristics of the radiance functions.
A discontinuity meshing radiosity algorithm was presented and results were
shown indicating that discontinuity meshing can generate remarkably accurate
solutions while at the same time drastically reducing the need for user interven-
tion.
Although much work remains to be done in order to design a robust and ef-
ficient discontinuity meshing algorithm, we believe that discontinuity meshing
represents an important step towards making radiosity practical for a wide vari-
ety of applications; the concept and ideas are sound and the benefits they provide
suggest that they should become an integral part of future radiosity systems.
Bibliography
[AIRE9O]
[AMAN84]
Airey, John M., John H. Rohif, and Frederick P. Brooks Jr. "Towards
Image Realism with Interactive Update Rates in Complex Virtual
Building Environments," Proceedings of the 1990 Symposium on
Interactive 3D Graphics (Snowbird, Utah, March 25-28, 1990), in
Computer Graphics, 24(2), March 1990, pages 41--H50.
Amantides, J. "Ray Tracing with Cones," Proceedings of SIG-
GRAPH'84 (Minneapolis, Minnesota, July 2?27, 1984), in Corn-
puter Graphics, 18(3), July 1984, pages 129--H135.
[ARVO86] Arvo, James. "Backward Ray Tracing," SIGGRAPH'86 Develop-
ments in Ray Tracing Course Notes, August 1986.
[ARVO87]
[ARVO9O]
Arvo, James and David Kirk. "Fast Ray Tracing by Ray Classifi-
cation," Proceedings of SIGGRAPH'87 (Anaheim, California, July
27--H31, 1987), in Computer Graphics, 21(4), July 1987, pages 55--H64.
Arvo, James and David Kirk. "Particle Transport and Image Syn-
thesis," Proceedings of SIGGRAPH'90 (Dallas, Texas, August ?10,
1990), in Computer Graphics, 24(4), August 1990, pages 6?6.
[ARVO9lj Arvo,James, editor. Graphics Gems II, Academic Press, Inc., Boston,
Massachusetts, 1991.
[BAUM86]
Baum, Daniel R., John R. Wallace, Michael F. Cohen, and Donald P.
Greenberg. "The Back-Buffer Algorithm; An Extension of the Ra-
diosity Method to Dynamic Environments," The Visual Computer,
2(5), September 1986, pages 29?308.
Baum, Daniel R., Holly E. Rushmeier, and James M. Winget. "Im-
proving Radiosity Solutions Through the Use of Analytically De-
termined Form-Factors," Proceedings of SIGGRAPH'89 (Boston,
Massachusetts, July 31--HAugust 4, 1989), in Computer Graphics,
23(3), July 1989, pages 325--H334.
[BAUM89]
166
167
[BAUM91]
[BERG86j
[BLIN77]
[BLIN88]
[BoUV9O]
[BREB89]
[BUCK89]
[CAMP9O]
Baum, Daniel R., Stephen Mann, Kevin P. Smith, and James M.
Winget. 6,Making Radiosity Usable: Automatic Preprocessing
and Meshing Techniques for the Generation of Accurate Radios-
ity Solutions," Proceedings of SIGGRAPH'91 (Las Vegas, Nevada,
July 2?August 2, 1991), in Computer Graphics, 25(4), July 1991,
pages 51?0.
Bergman, Larry, Henry Fuchs, Eric Grant, and Susan Spach. "Im-
age Rendering by Adaptive Refinement," Proceedings of SIG-
GRAPH'86 (Dallas, Texas, August 1?22, 1986), in Computer Graph-
ics, 20(4), August 1986, pages 29--H38.
Blinn, Jim F. "Models of Light Reflection for Computer Syn-
thesized Pictures," Proceedings of SIGGRAPH'77 (San Jose, Ca-
lifornia, July 2()?22, 1977), in Computer Graphics, 11(3), July 1977.
pages 192--H198.
Blinn, James F. "Jim Blinn's Corner: Me and My (Fake)
Shadow," IEEE Computer Graphics and Applications, 8(1), January
1988, pages 82--H86.
Bouville, Christian, Kadi Bouatouch, Pierre Tellier, and Xavier
Pueyo. "Theoretical Analysis of Global Illumination Models," in
Proceedings of Eurographics Workshop on Photosimulation, Real-
ism and Physics in Computer Graphics (Rennes, France, June 11--H
12, 1990), June 1990, pages 53A)6.
Brebbia, C. A. and J. Dominguez. Boundary Elements: An Introduc-
tory Course, Computational Mechanics Publications, Southampton
Boston, Massachusetts, 1989.
Buckalew, Chris and Donald Fussell. "Illumination Networks:
Fast Realistic Rendering with General Reflectance Functions,"
Proceedings of SIGGRAPH'89 (Boston, Massachusetts, July 31--H
August4,1989), in Computer Graphics, 23(3), July 1989, pages 89--H98.
Campbell, A. T., III and Donald 5. Fussell. "Adaptive Mesh
Generation for Global Diffuse Illumination," Proceedings of SIG-
GRAPH'90 (Dallas, Texas, August 6--H10, 1990), in Computer Graph-
ics, 24(4), August 1990, pages 155--H164.
Campbell, A. T., III and Donald 5. Fussell. An Analytic Approach
to Illumination with Area Light Sources, Technical Report TR-91-25,
Department of Computer Sciences, University of Texas, Austin,
Texas, August 1991.
[CAMP91a]
168
[CAMP91b]
[CAMP91c]
[CEND87]
[CHAT87j
[CHEN90aj
Campbell, A. T., III and Donald 5. Fussell. Analytic Illumination
with Polygonal Light Sources, Technical Report TR-91-15, Depart-
ment of Computer Sciences, University of Texas, Austin, Texas,
April 1991.
Campbell, III, A. T. Modeling Global Difuse Illumination for Image
Synthesis, PhD dissertation, University of Texas at Austin, Austin,
Texas, December 1991.
Cendes, Zoltan J. and Steven H. Wong. "C' Quadratic Interpola-
tion Over Arbitrary Point Sets," IEEE Computer Graphics and Appli-
cations, 7(11), November 1987, pages ?16.
Chattopadhyay, Subdeb and Akira Fujimoto. "Bi-directional Ray
Tracing," in Proceedings of CC International'87, 1987, pages 335
343.
Chen, Hong and En-Hua Wu. "An Adapted Solution of Progres-
sive Radiosity and Ray-Tracing Methods for Non-diffuse Environ-
ments," in Proceedings of CC International'90: Computer Graph-
ics Around the World1 1990, pages 477?90.
[CHEN90b] Chen, Shenchang Eric. "Incremental Radiosity: An Extension of
Progressive Radiosity to an Interactive Image Synthesis System,"
Proceedings of SIGGRAPH'90 (Dallas, Texas, August 6--H10, 1990),
in Computer Graphics, 24(4), August 1990, pages 135-144.
[CHEN91]
[CHIN89]
[CHIN92j
Chen, Shenchang Eric, Holly E. Rushmeier, Gavin Miller, and Dou-
glass Turner. "A Progressive Multi-Pass Method for Global Il-
lumination," Proceedings of SIGGRAPH'91 (Las Vegas, Nevada,
July 2?August 2, 1991), in Computer Graphics, 25(4), July 1991,
pages 165174.
Chin, Norman and Steven Feiner. "Near Real-Time Shadow Gen-
eration Using BSP Trees," Proceedings of SIGGRAPH'89 (Boston,
Massachusetts, July 31--HAugust 4, 1989), in Computer Graphics,
23(3), July 1989, pages 99--H106.
Chin, Norman and Steven Feiner. "Fast Object-Precision Shadow
Generation for Area Light Sources Using BSP Trees," in Proceed-
ings of 1992 Symposium on Interactive 3D Graphics (Cambridge,
Massachusetts, March 29--HApril 1, 1992), March 1992.
[COHE85] Cohen, Michael F. and Donald P. Greenberg. "The Hemi-Cube:
A Radiosity Solution for Complex Environments," Proceedings
169
[COHE86]
[COHE88a]
[COHE88b]
[CCOK81]
of SIGGRAPH'85 (San Francisco, California, July 22--H26, 1985), in
Computer Graphics, 19(3), July 1985, pages 31?O.
Cohen, Michael F, Donald P. Greenberg, and David 5. Immel. "An
Efficient Radiosity Approach for Realistic Image Synthesis," IEEE
Computer Graphics and Applications, 6(2), March 1986, pages 2?35.
Cohen, Michael F. "A Consumer's and Developer's Guide to Ra-
diosity," SIGGRAPH'88 A Consumer's and Developer's Guide to
Image Synthesis Course Notes, August 1988.
Cohen, Michael F., Shenchang Eric Chen, John R. Wallace, and
Donald P. Greenberg. "A Progressive Refinement Approach to
Fast Radiosity Image Generation," Proceedings of SIGGRAPH'88
(Atlanta, Georgia, August 1--H5, 1988), in Computer Graphics, 22(4),
August 1988, pages 75--H84.
Cook, Robert L. and Kenneth E. Torrance. "A Reflectance Model
for Computer Graphics," Proceedings of SIGGRAPH'81 (Dallas,
Texas, August 3--H7, 1981), in Computer Graphics, 15(3), August 1981,
pages 307--H316.
[C00K82] Cook, Robert L. "A Reflectance Model for Computer Graphics,"
ACM Transactions on Graphics, 1(1), January 1982, pages 7--H24.
[COOK84a]
[COOK84bj
[COOK86a]
Cook, Rob L. "Shade Trees," Proceedings of SIGGRAPH'84 (Min-
neapolis, Minnesota, July 23--H27, 1984), in Computer Graphics, 18(3),
July 1984, pages 223--H231.
Cook, Robert L. "Distributed Ray Tracing," Proceedings of SIG-
GRAPH'84 (Minneapolis, Minnesota, July 2?27, 1984), in Com-
puter Graphics, 18(3), July 1984, pages 137--H146.
Cook, Robert L. "Practical Aspects of Distributed Ray Tracing,"
SIGGRAPH'86 Developments in Ray Tracing Course Notes, Au-
gust 1986.
[COOK86b] Cook, Robert L. "Stochastic Sampling in Computer Graphics,"
ACM Transactions on Graphics, 5(1), January 1986, pages 51--H72.
[CORT82] Cortey, Noel E. "Interpolation of Arbitrary Spaced Points by
Closed Surfaces," computer & graphics, 6(1), 1982, pages 351--H364.
Crow, Franklin C. "Shadow Algorithms for Computer Graphics,"
Proceedings of SIGGRAPH'77 (San Jose, California, July 2022,
1977), in Computer Graphics, 11(2), July 1977, pages 242--H248.
[CROW77]
170
[CROW8l]
[CROW84]
[DIPP85]
[DL79]
[DL8lj
[DL82]
[DORS9l]
[DUFF79]
[FARI9O]
[FOLE9O]
Crow, Franklin C. `6A Comparison of Antialiasing Techniques,"
IEEE Computer Graphics and Applications, 1(1), January 1981,
pages 4048.
Crow, Franklin C. "Summed-Area Tables for Texture Mapping,"
Proceedings of SIGGRAPH'84 (Minneapolis, Minnesota, July 23
27, 1984), in Computer Graphics, 18(3), July 1984, pages 207--H212.
Dippe, Mark A. Z. and Erling Henry Wold. "Antialiasing Through
Stochastic Sampling," Proceedings of SIGGRAPH'85 (San Fran-
cisco, California, July 22--H26, 1985), in Computer Graphics, 19(3), July
1985, pages 69--H78.
Di Laura, David L. "On a New Technique for Interreflected Com-
ponent Calculations," Journal of the illumination Engineering Society,
1979, pages 5359.
Di Laura, David L. "On Radiative Transfer Calculations in Un-
empty Rooms," Journal of the illumination Engineering Society, 1981,
pages 12O?126.
Di Laura, David L. "On the Simplification of Radiative Transfer
Calculations," Journal of the Illumination Engineering Society, 1982,
pages 12--H16.
Dorsey, Julie O'B., Fran?is X. Sillion, and Donald P. Greenberg.
"Design and Simulation of Opera Lighting and Projection Ef-
fects," Proceedings of SIGGRAPH'91 (Las Vegas, Nevada, July 2?
August 2, 1991), in Computer Graphics, 25(4), July 1991, pages 41--H50.
Duff, Tom. "Smoothly Shaded Renderings of Polygonal Objects
on Raster Displays," Proceedings of SIGGRAPH'79 (Chicago, Illi-
nois, August ?10, 1979), in Computer Graphics, 13(2), August 1979,
pages 270--H275.
Farin, Gerald. Curves and Surfaces for Computer Aided Geometric
Design, Computer Science and Scientific Computing, Academic Press,
Inc., San Diego, California, second edition, 1990.
Foley, James D., Andries van Dam, Steven K. Feiner, and John F.
Hughes. Computer Graphics, Principles and Practice, Addison-
Wesley Publishing Company, Reading, Massachusetts, second edi-
tion, 1990.
171
[FOUR88j
Fournier, Alain and Eugene Fiume. "Constant-Time Filtering with
Space-Variant Kernels," Proceedings of SIGGRAPH'88 (Atlanta,
Georgia, August 1--H5, 1988), in Computer Graphics, 22(4), August
1988, pages 229--H238.
[FRAN82] Franke, Richard. "Scattered Data Interpolation: Tests of Some
Methods," Mathematics of Cornputation, 38(157), January 1982.
[FUCH80]
[FUJ186]
[GEOR90]
[GIGU88J
[GIGU90]
[GIGU91]
[GLAS84]
[GLAS86j
Fuchs, Henry, Zvi M. Kedem, and Bruce Naylor. "On Visible Sur-
face Generation by a Priori Tree Structures," Proceedings of SIG-
GRAPH'80 (Seattle, Washington, July 1?18, 1980), in Computer
Graphics, 14(3), June 1980, pages 17?181.
Fujimoto, Akira, Tanaka Takayuki, and Iwata Kansei. "ARTS: Ac-
celerated Ray-Tracing System," IEEE Computer Graphics and Appli-
cations, 6(4), April 1986, pages 16--H26.
George, David W., Francois X. Sillion, and Donald P. Greenberg.
"Radiosity Redistribution for Dynamic Environments," IEEE Com-
puter Graphics and Applications, 10(4), July 1990, pages 2?35.
Gigus, Ziv, John Canny, and Raimund Seidel. "Efficiently Com-
puting and Representing Aspect Graphs of Polyhedral Objects," in
Proceedings of the Second International Conference on Computer
Vision, December 1988, pages 30--H39.
Gigus, Ziv and Jitendra Malik. "Computing the Aspect Graph for
Line Drawings of Polyhedral Objects," IEEE Transactions on Pattern
Analysis and Machine Intelligence, 12(2), February 1990, pages 113?
122.
Gigus, Ziv, John Canny, and Raimund Seidel. "Efficiently Com-
puting and Representing Aspect Graphs of Polyhedral Objects,"
IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(6),
June 1991, pages 542--H551.
Glassner, Andrew 5. "Space Subdivision for Fast Ray Trac-
ing," IEEE Computer Graphics and Applications, 4(10), October 1984,
pages 15--H22.
Glassner, Andrew. "Adaptive Precision in Texture Mapping," Pro-
ceedings of SIGGRAPH'86 (Dallas, Texas, August 1?22, 1986), in
Computer Graphics, 20(4), August 1986, pages 297--H306.
[GLAS89] Glassner, Andrew 5., editor. An Intrnduction to Ray Tracing, Aca-
demic Press, Inc., San Diego, California, 1989.
172
[GLAS90] Glassner, Andrew 5., editor. Graphics Gems, Academic Press, Inc.,
Boston, Massachusetts, 1990.
[GOLD87]
[GOLU89]
[GORA84]
[GORD91]
Goldsmith, Jeffrey and John Salmon. 1'Automatic Creation of Ob-
ject Hierarchies for Ray Tracing," IEEE Computer Graphics and Ap-
plications, 7(5), May 1987, pages 1?20.
Golub, Gene H. and Charles F. Van Loan. Matrix Computations,
The Johns Hopkins University Press, Baltimore, Maryland, second
edition, 1989.
Goral, Cindy M., Kenneth E. Torrance, Donald P. Greenberg, and
Bennett Battaile. "Modeling the Interaction of Light Between Dif-
fuse Surfaces," Proceedings of SIGGRAPH'84 (Minneapolis, Min-
nesota, July 23--H27, 1984), in Computer Graphics, 18(3), July 1984,
pages 213222.
Gordon, Dan and Shuhong Chen. "Front-to-Back Display of BSP
Trees," IEEE Computer Graphics and Applications, 11(5), September
1991, pages 79--H85.
[GOUR71] Gouraud, H. "Continuous Shading of Curved Surfaces," IEEE
Transactions on Computers, C-20(6), June 1971, pages 623628.
[GREE86]
[GREE89]
[GREE90]
[GREE91]
[HAIN86]
Greenberg, Donald P., Michael F. Cohen, and Kenneth E. Torrance.
"Radiosity: A Method for Computing Global Illumination," The
Visual Computer, 2(5), September 1986, pages 291--H297.
Greenberg, Donald R, Michael Cohen, Roy Hall, Holly Rushmeier,
and John Wallace. "SIGGRAPH'89 Radiosity Course Notes," Au-
gust 1989.
Greenberg, Donald P., Michael Cohen, David George, Holly Rush-
meier, John Wallace, and Greg Ward. "SIGGRAPH'90 Radiosity
Course Notes," August 1990.
Greenberg, Donald P., Michael Cohen, Holly Rushmeier, Fran?is
Sillion, and John Wallace. "SIGGRAPH'91 Radiosity Course
Notes," August 1991.
Haines, Eric A. and Donald P. Greenberg. "The Light buffer: a
Shadow Testing Accelerator," IEEE Computer Graphics and Applica-
tions, 6(9), September 1986, pages ?16.
[HAIN91a] Haines, Eric A. "Beams 0' Light: Confessions of a Hacket" SIG-
GRAPH'91 Frontiers in Rendering Course Notes, July 1991.
173
[HAIN91b] Haines, Eric A. "Ronchamp: A Case Study for Radiosity,', SIG-
GRAPH'91 Frontiers in Rendering Course Notes, July 1991.
[HAIN91c]
[HALL83]
Haines, Eric A. and John R. Wallace. "Shaft Culling for Efficient
Ray-Traced Radiosity," in Proceedings of the Second Eurographics
Workshop on Rendering (Barcelona, Spain, May 1315,1991), May
1991.
Hall, Roy A. and Donald P. Greenberg. "A Testbed for Realistic
Image Synthesis," IEEE Computer Graphics and Applications, 3(8),
November 1983, pages 1(Y?20.
[HALL89] Hall, Roy. Illumination and Color in Computer Generated Image??,
Springer-Verlag, New York, New York, 1989.
[HANR90]
Hanrahan, Pat and David Salzman. "A Rapid Hierarchical Radios-
ity Algorithm for Unoccluded Environments," in Proceedings of
Eurographics Workshop on Photosimulation, Realism and Physics
in Computer Graphics (Rennes, France, June 11--H12, 1990), June
1990, pages 151--H172.
[HANR9la] Hanrahan, Pat. "Rapid Hierarchical Global Illumination Algo-
rithms," SIGGRAPH'91 Frontiers in Rendering Course Notes, July
1991.
[HANR91b] Hanrahan, Pat, David Salzman, and Larry Aupperle. "A Rapid
Hierarchical Radiosity Algorithm," Proceedings of SIGGRAPH'91
(Las Vegas, Nevada, July 2?August 2,1991), in Computer Graphics,
25(4), July 1991, pages 197--H206.
[HE91]
[HECK84]
[HECK86a]
He, Xiao D., Kenneth E. Torrance, Fran?ois X Sillion, and Donald P.
Greenberg. "A Comprehensive Physical Model for Light Reflec-
tion," Proceedings of SIGGRAPH'91 (Las Vegas, Nevada, July 2?
August 2, 1991), in Computer Graphics, 25(4), July 1991, pages 17?
186.
Heckbert, Paul 5. and Pat Hanrahan. "Beam Tracing Polygonal
Objects," Proceedings of SIGGRAPH'84 (Minneapolis, Minnesota,
July 2327, 1984), in Computer Graphics, 18(3), July 1984, pages 119--H
128.
Heckbert, Paul 5. "Filtering by Repeated Integration," Proceed-
ings of SIGGRAPH'86 (Dallas, Texas, August 1?22, 1986), in Com-
puter Graphics, 20(4), August 1986, pages 315--H322.
174
[HECK86b] Heckbert, Paul 5. "Survey of Texture Mapping," IEEE Computer
Graphics and Applications, 6(11), November 1986, pages 5?7.
[HECK9O]
[HECK91a]
[HECK91b]
[HECK92a]
Heckbert, Paul 5. "Adaptive Radiosity Textures for Bidirec-
tional Ray Tracing," Proceedings of SIGGRAPH'90 (Dallas, Texas,
August ?10, 1990), in Computer Graphics, 24(4), August 1990,
pages 145--H154.
Heckbert, Paul 5. Simulating Global Illumination Using Adaptive
Meshing, PhD dissertation, Department of EECS, UC Berkeley, Ca-
lifornia, June 1991.
Heckbert, Paul 5. and James M. Winget. Finite Element Methods for
Global illumination, Technical Report UCB/CSD 91/643, CS Divi-
sion, University of California at Berkeley, Berkeley, California, July
1991.
Heckbert, Paul 5. "Discontinuity Meshing for Radiosity," in Pro-
ceedings of the Third Eurographics Workshop on Rendering (Bris-
tol, UK, May 1820, 1992), May 1992, pages 203--H216.
[HECK92b] Heckbert, Paul 5. "Finite Element Methods for Radiosity," SIG-
GRAPH'92 Global Illumination Course Notes, July 1992.
[HECK92c] Heckbert, Paul 5. "Radiosity in Elatland," in Proceedings of Euro-
graphics'92, September 1992.
[H0RN75] Hornbeck, Robert W. Numerical Methods, Quantum Publishers,
New York, New York, 1975.
[HOTT67j Hottel, Hoyt C. and Adel F. Sarofim. Radiative Transfrr, McGraw-
Hill, New York, New York, 1967.
[HOU78]
Hou, Hsieh 5. and Harry C. Andrews. "Cubic Splines for Image
Interpolation and Digital Filtering," IEEE Transactions on Acous-
tic, Speech, and Signal Processing, ASSP-26(6), December 1978,
pages 508517.
[HOWE82j Howell, John R. A Catalog of Radiation Configuration Factors,
McGraw-Hill, New York, New York, 1982.
[IES9l] lES. lES Lighting Handbook, Illuminating Engineering Society of
North America, New York, New York, 1991.
[IMME86] Immel, David 5., Michael F. Cohen, and Donald P. Greenberg. "A
Radiosity Method for Non-Diffuse Environments," Proceedings of
175
[KAJ186]
[KAY86]
[KEYS8l]
[KIRK91]
[KOK9O]
[KOK91j
[KRIE89]
[KRIE90]
SIGGRAPH'86 (Dallas, Texas, August 1?22, 1986), in Computer
Graphics, 20(4), August 1986, pages 133142.
Kajiya, James T. "The Rendering Equation," Proceedings of SIG-
GRAPH'86 (Dallas, Texas, August 1?22, 1986), in Computer Graph-
ics, 20(4), August 1986, pages 143-150.
Kay, Timothy L. and James t Kajiya. "Ray Tracing Complex
Scenes," Proceedings of SIGGRAPH'86 (Dallas, Texas, August 1?
22, 1986), in Computer Graphics, 20(4), August 1986, pages 269--H278.
Keys, Robert G. "Cubic Convolution Interpolation for Digital Im-
age Processing," IEEE Transactions on Acoustic, Speech, and Signal
Processing, ASSP-29(6), December 1981, pages 1153-1160.
Kirk, David and James Arvo. "Unbiased Sampling Techniques
for Image Synthesis," Proceedings of SIGGRAPH'91 (Las Vegas,
Nevada, July 2?August 2, 1991), in Computer Graphics, 25(4), July
1991, pages 153156.
Kok, Arjan J. E, Celal Yilmaz, and Laurens H. J. Bierens. "A Two-
Pass Radiosity Method for Bezier Patches," in Proceedings of Eu-
rographics Workshop on Photosimulation, Realism and Physics in
Computer Graphics (Rennes, France, June 11--H12, 1990), June 1990,
pages 117--H126.
Kok, Arjan J. F. and Frederik Jansen. "Source Selection for the Di-
rect Lighting Computation in Global Illumination," in Proceedings
of the Second Furographics Workshop on Rendering (Barcelona,
Spain, May 13-15, 1991), May 1991.
Kriegman, D. J. and J. Ponce. "Computing Exact Aspect Graphs of
Curved Objects: Solids of Revolution," in Proceedings of the IEEE
Workshop on the Interpretation of 3D Scenes, IEEE Society Press,
New York, New York, November 1989, pages 116--H121.
Kriegman, D. J. and J. Ponce. "Computing Exact Aspect
Graphs of Curved Objects: Parametric Surfaces," in Proceed-
ings of the Eighth National Conference on Artificial Intelligence,
AAAI Press/MIT Press, Cambridge, Massachusetts, June 1990,
pages 107?1079.
Lee, Mark E., Richard A. Redner, and Samuel P. Uselton. "Statisti-
cally Optimized Sampling for Distributed Ray Tracing," Proceed-
ings of SIGGRAPH'85 (San Francisco, California, July 22--H26,1985),
in Computer Graphics, 19(3), July 1985, pages 61--H68.
[LEE85]
176
[LEE90]
[L1SC90]
[LISC91]
[LISC92]
[LISC93]
[LS90j
[MAX90]
[MAXW86]
[MIST87]
[MITC87]
Lee, Mark E. and Richard A. Redner. "A Note on the Use of Non-
linear Filtering in Computer Graphics," IEEE Computer Graphics
and Applications, 10(3), May 1990, pages 23--H29.
Lischinski, Daniel and Jakob Gonczarowski. "Improved Tech-
niques for Ray Tracing Parametric Surfaces," The Visual Computer,
6(3), 1990, pages 13?152.
Lischinski, Dani, Filippo Tampieri, and Donald P. Greenberg. Im-
proving Sampling and Reconstruction Techniques for Radiosity, Techni-
cal Report 91-1202, Department of Computer Science, Cornell Uni-
versity, Ithaca, New York, August 1991
Lischinski, Dani, Filippo Tampieri, and Donald P. Greenberg.
"Discontinuity Meshing for Accurate Radiosity," IEEE Computer
Graphics and Applications, 12(6), November 1992, pages 2539.
Lischinski, Dani, Filippo Tampieri, and Donald P. Greenberg.
"Combining Hierarchical Radiosity and Discontinuity Meshing,"
Proceedings of SIGGRAPH'93 (Anaheim, California, August 16,
1993), in Computer Graphics, 27(4), August 1993.
Le Saec, Bertrand and Christophe Schlick. "A Progressive Ray
Tracing Based Radiosity with General Reflectance Functions," in
Proceedings of Eurographics Workshop on Photosimulation, Real-
ism and Physics in Computer Graphics (Rennes, France, June 11--H
12, 1990), June 1990, pages 103--H116.
Max, Nelson. "Cone-Spheres," Proceedings of SIGGRAPH'90
(Dallas, Texas, August 6--H10, 1990), in Computer Graphics, 24(4), Au-
gust 1990, pages 59--H62.
Maxwell, Gregory M., Michael J. Bailey, and Victor W. Gold-
schmidt. "Calculations of the Radiation Configuration Factor Us-
ing Ray Casting," Computer-Aided Design, 18(7), September 1986,
pages 371--H379.
Mistrick, R. G. and D. L. Di Laura. "A New Finite Orthogonal
Transform Applied to Radiative Transfer Calculations," !ournal of
the Illumination Engineering Society, 1987, pages 115128.
Mitchell, Don P. "Generating Antialiased Images at Low Sampling
Densities," Proceedings of SIGGRAPH'87 (Anaheim, California,
July 27--H31, 1987), in Computer Graphics, 21(4), July 1987, pages 65
72.
177
[MITC88]
[MITC91]
[NAIM87]
[NAKA90]
[NAYL90]
[NEUM90j
[NEUM91]
Mitchell, Don P. and Arun N. Netravali. "Reconstruction Filters
in Computer Graphics," Proceedings of SIGGRAPH'88 (Atlanta,
Georgia, August 1--H5, 1988), in Computer Graphics, 22(4), August
1988, pages 221--H228.
Mitchell, Don P. "Spectrally Optimal Sampling for Distributed
Ray Tracing," Proceedings of SIGGRAPH'91 (Las Vegas, Nevada,
July 2?August 2, 1991), in Computer Graphics, 25(4), July 1991,
pages 157--H164.
Naiman, Avi and Alain Fournier. "Rectangular Convolution for
Fast Filtering of Characters," Proceedings of SIGGRAPH'87 (Ana-
heim, California, July 27--H31, 1987), in Computer Graphics, 21(4),July
1987, pages 233--H242.
Nakamae, Eihachiro, Kazufumi Kaneda, Takashi Okamoto, and
Tomoyuki Nishita. "A Lighting Model Aiming at Drive Simula-
tors," Proceedings of SIGGRAPH'90 (Dallas, Texas, August 6--H10,
1990), in Computer Graphics, 24(4), August 1990, pages 39?04.
Naylor, Bruce, John Amanatides, and William Thibault. "Merging
BSP Trees Yields Polyhedral Set Operations," Proceedings of SIG-
GRAPH'90 (Dallas, Texas, August 6--H10, 1990), in Computer Graph-
ics, 24(4), August 1990, pages 11?124.
Neumann, Laszlo and Attila Neumann. "Efficient Radiosity Meth-
ods for Non-Separable Reflectance Models," in Proceedings of Eu-
rographics Workshop on Photosimulation, Realism and Physics in
Computer Graphics (Rennes, France, June 11--H12, 1990), June 1990,
pages 83--H102.
Neumann, Laszlo and Csaba Kelemen. "Solution of Interreflec-
tion Problem for Very Complex Environments by Transillumina-
tion Method," in Proceedings of the Second Eurographics Work-
shop on Rendering (Barcelona, Spain, May 13--H15,1991), May 1991.
[NEWM79] Newman, William M. and Robert F. Sproull. Principles ofInteractive
Computer Graphics, McGraw-Hill, Tokyo, Japan, 1979.
Nishita, Tomoyuki and Eihachiro Nakamae. "Half-Tone Repre-
sentation of 3-D Objects Illuminated by Area Sources or Polyhe-
dron Sources," in Proceedings of the IEEE Computer Society's
International Computer Software and Applications Conference
(COMPSAC83) (Chicago, Illinois, November 1983), November
1983, pages 237--H241.
[NISH83j
178
[NISH85a]
[NISH85b]
[NISH86]
[NISH87]
[PAIN89j
[PARK83]
[PEAC85j
[PERL85]
Nishita, Tomoyuki and Eihachiro Nakamae. "Continuous Tone
Representation of Three-Dimensional Objects Taking Account of
Shadows and Interreflections," Proceedings of SIGGRAPH'85 (San
Francisco, California, July 22--H26, 1985), in Computer Graphics, 19(3),
July 1985, pages 2330.
Nishita, Tomoyuki, Isao Okamura, and Eihachiro Nakamae.
"Shading Models for Point and Linear Light Sources," ACM Trans-
actions on Graphics, 4(2), April 1985, pages 12?146.
Nishita, Tomoyuki and Eihachiro Nakamae. "Continuos Tone
Representation of Three-Dimensional Objects Illuminated by Sky
Light," Proceedings of SIGCRAPH'86 (Dallas, Texas, August 1?
22, 1986), in Computer Graphics, 20(4), August 1986, pages 12?132.
Nishita, Tomoyuki, Yasuhiro Miyawaki, and Eihachiro Nakamae.
"A Shading Model for Atmospheric Scattering Considering Lumi-
nous Intensity Distribution of Light Sources," Proceedings of SIG-
GRAPH'87 (Anaheim, California, July 27--H31, 1987), in Computer
Graphics, 21(4), July 1987, pages 303-310.
Painter, James and Kenneth Sloan. "Antialiased Ray Tracing by
Adaptive Progressive Refinement," Proceedings of SIGGRAPH'89
(Boston, Massachusetts, July 31--HAugust 4, 1989), in Computer
Graphics, 23(3), July 1989, pages 281--H288.
Park, Stephen K. and Robert A. Schowengerdt. "Image Re-
construction by Parametric Cubic Convolution," Computer Vision,
Graphics, and Image Processing, 23(3), September 1983, pages 25?
272.
Peachey, Darwyn R. "Solid Texturing of Complex Surfaces," Pro-
ceedings of SIGGRAPH'85 (San Francisco, California, July 22--H26,
1985), in Computer Graphics, 19(3), July 1985, pages 279--H288.
Perlin, Ken. "An Image Synthesizer," Proceedings of SIC-
GRAPH'85 (San Francisco, California, July 22--H26, 1985), in Com-
puter Graphics, 19(3), July 1985, pages 287--H296.
[PHON75j Phong, Bui-Tuong. "Illumination for Computer-Generated Pic-
tures," Communications of the ACM, 18(6), June 1975, pages 311--H317.
Press, WilliamH., Brian P. Flannery, Saul A. Teukolsky, and
William T. Vetterling. Numerical Recipes in C, Cambridge Univer-
sity Press, Cambridge, Massachusetts, 1990.
[PRES90j
179
[RATL72j Ratliff, Floyd. "Contour and Contrast," Scientific American, June
1972, pages 90--H101.
[REIC92]
[RIVA86]
[RIVA87]
[RUSH87]
[RUSH90a]
[RUSH90b]
[SABL87]
Reichert, Mark C. A Two-Pass Radiosity Method Driven by Lights and
Viewer Position, Master's thesis, Program of Computer Graphics,
Cornell University, Ithaca, New York, January 1992.
Rivara, Maria-Cecilia. "Adaptive Finite Element Refinement and
Fully Irregular and Conforming Triangulations," in Babuska, I., J.
Gago, E. R. de Oliveira, and 0. C. Zienkiewicz, editors, Accuracy
Estimates and Adaptive R?nements in Finite Element Computations,
John Wiley and Sons, 1986, pages 359--H370.
Rivara, Maria-Cecilia. "Numerical Generation of Nested Series of
General Triangular Grids," in Chui, C. K., L. L. Schumaket and
F. I. Utreras, editors, Topics in Multivariate Approximation, Academic
Press, Inc., Orlando, Florida, 1987, pages 193--H206.
Rushmeier, Holly E. and Kenneth E. Torrance. "The Zonal Method
for Calculating Light Intensities in the Presence of a Participating
Medium," Proceedings of SIGGRAPH'87 (Anaheim, California,
July 27--H31,1987), in Computer Graphics, 21(4),July 1987, pages 29?
302.
Rushmeier, Holly E., Daniel R. Baum, and David E. Hall. "Accel-
erating the Hemi-Cube Algorithm for Calculating Radiation Form
Factors," in 5th AIAA/ASME Thermophysics and Heat Transfer
Conference (Seattle, Washington, June 1990), 1990.
Rushmeier, Holly E. and Kenneth E. Torrance. "Extending the
Radiosity Method to Include Specularly Reflecting and Translu-
cent Materials," ACM Transactions on Graphics, 9(1), January 1990,
pages 1--H27.
Sablonniere, Paul. "Composite Finite Elements of Class C2," in
Chui, C. K., L. L. Schumaker, and F. I. Utreras, editors, Topics in
MultivariateApproximation, Academic Press, Inc., Orlando, Florida,
1987, pages 207--H218.
Salesin, David and Jorge Stolfi. "The ZZ-Buffer: A Simple and Effi-
cient Rendering Algorithm with Reliable Antialising," in Proceed-
ings of the PIXIM'89 Conference (Hermes Editions, Paris), Septem-
ber 1989, pages 451M66.
[SALE89]
180
[SALE91]
[SALE92j
Salesin, David. Epsilon-Geometry: Building Robust Algorithms from
Imprecise Computations, PhD dissertation, Department of Com-
puter Science, Stanford University, Palo Alto, California, 1991.
Salesin, David, Dani Lischinski, and Tony DeRose. "Reconstruct-
ing Illumination Functions with Selected Discontinuities," in Pro-
ceedings of the Third Eurographics Workshop on Rendering (Bris-
tol, UK, May 1?20, 1992), May 1992, pages 99--H112.
[SAME90a] Samet, Hanan. Applications of Spatial Data Structures, Addison-
Wesley Publishing Company, Inc., Reading, Massachusetts, 1990.
[SAME90b] Samet, Hanan. The Design and Analysis of Spatial Data Struc-
tures, Addison-Wesley Publishing Company, Inc., Reading, Mas-
sachusetts, 1990.
[SCHM86]
[SCHU87j
[SEDE84j
[SEDE90]
[SEGA90]
[SEQU89]
Schmitt, Francis J. M., Brian A. Barsky, and Wen-Hui Du. "An
Adaptive Subdivision Method for Surface-Fitting from Sampled
Data," Proceedings of SIGCRAPH'86 (Dallas, Texas, August 1?
22, 1986), in Computer Graphics, 20(4), August 1986, pages 179--H188.
Schumaker, Larry L. "Triangulation Methods," in Chui, C. K., L. L.
Schumaker, and F. I. Utreras, editors, Topics in Multivariate Approx-
imation, Academic Press, Inc., Orlando, Florida, 1987, pages 219--H
232.
Sederberg, Thomas and David C. Anderson. "Ray Tracing of
Steiner Patches," Proceedings of SIGGRAPH'84 (Minneapolis,
Minnesota, July 23--H27,1984), in Computer Graphics, 18(3), July1984,
pages 159--H164.
Sedeberg, Thomas W. "Techniques for Cubic Algebraic Surfaces,"
IEEE Computer Graphics and Applications, 10(4), July 1990, pages 1?
25.
Segal, Mark. "Using Tolerances to Guarantee Valid Polyhedral
Modeling Results," Proceedings of SIGGRAPH'90 (Dallas, Texas,
August 6--H10, 1990), in Computer Graphics, 24(4), August 1990,
pages 105--H114.
Sequin, Carlo H. and Eliot K. Smyrl. "Parametrized Ray-Tracing,"
Proceedings of SICCRAPH'89 (Boston, Massachusetts, July 31--H
August 4, 1989), in Computer Graphics, 23(3), July 1989, pages 307--H
314.
181
[SHAO88]
[SHIN87]
[SHIR9Oa]
[SHIR9Ob]
[SHIR9Oc]
[SHIR91j
[SIEG81
[SILL89]
[SILL91]
Shao, Min-Zhi, Qun-Sheng Peng, and You-Dong Liang. "A New
Radiosity Approach by Procedural Refinements for Realistic Im-
age Synthesis," Proceedings of SIGGRAPH'88 (Atlanta, Geor-
gia, August 1--H5, 1988), in Computer Graphics, 22(4), August 1988,
pages 93102.
Shinya, Mikio, Tokiichiro Takahashi, and Seiichiro Naito. "Prin-
ciples and Applications of Pencil Tracing," Proceedings of SIG-
GRAPH'87 (Anaheim, California, July 27--H31, 1987), in Computer
Graphics, 21(4), July 1987, pages 4?54.
Shirley, Peter. Physically Based Lighting Calculations for Computer
Graphics, PhD dissertation, Department of Computer Science, Uni-
versity of Illinois, Urbana-Champaign, Illinois, 1990.
Shirley, Peter. "Physically Based Lighting Calculations for Com-
puter Graphics: a Modern Perspective," in Proceedings of ?uro-
graphics Workshop on Photosimulation, Realism and Physics in
Computer Graphics (Rennes, France, June 11--H12, 1990), June 1990,
pages 67--H82.
Shirley, Peter. "A Ray Tracing Method for Illumination Calcula-
tion in Diffuse-Specular Scenes," in Proceedings of Graphics Inter-
face'90 (Halifax, Nova Scotia, May 1990), May 1990, pages 205--H212.
Shirley, Peter. "Time Complexity of Monte Carlo Radiosity," in
Proceedings of Eurographics'91 (Vienna, Austria, September 16,
1991), September 1991, pages 459-t65.
Siegel, Robert and John R. Howell. Thermal Radiation Heat Transfer,
Hemisphere Publishing Corp., Washington D.C., 1981.
Sillion, Fran?is and Claude Puech. "A General Two-Pass Method
Integrating Specular and Diffuse Reflection," Proceedings of SIG-
GRAPH'89 (Boston, Massachusetts, July 31--HAugust 4, 1989), in
Computer Graphics, 23(3), July 1989, pages 335--H344.
Sillion, Fran?ois X, James R. Arvo, Stephen H. Westin, and Don-
ald P. Greenberg. "A Global Illumination Solution for General Re-
flectance Distributions," Proceedings of SIGGRAPH'91 (Las Vegas,
Nevada, July 2?August 2, 1991), in Computer Graphics, 25(4), July
1991, pages 187--H196.
Smith, Alvy Ray. "Planar 2-Pass Texture Mapping and Warping,"
Proceedings of SIGGRAPH'87 (Anaheim, California, July 27--H31,
1987), in Computer Graphics, 21(4), July 1987, pages 263272.
[SMIT87]
182
[SMIT92]
[SPAR63j
Smits, Brian E., James R. Arvo, and David H. Salesin. "An
Importance-Driven Radiosity Algorithm," Proceedings of SIG-
GRAPH'92 (Chicago, Illinois, July 2?31, 1992), in Computer Graph-
ics, 26(4), July 1992, pages 273282.
Sparrow, Ephraim M. "On the Calculation of Radiant Interchange
between Surfaces," in Ibele, Warren E., editor, Modern Developments
in Heat Transfer, Academic Press, New York, New York, 1963.
[SPAR78] Sparrow, E. M. and R. D. Cess. Radiation Heat Transfer, Hemisphere
Publishing Corp., Washington D.C., 1978.
[SRIP89]
[STEW88j
[STRO86j
[SUTH74]
[TAKA9O]
[TAMM85]
Sripradisvarakul, Thawach and Ramesh Jain. "Generating Aspect
Graphs for Curved Objects," in Proceedings of the IFEE Work-
shop on the Interpretation of 3D Scenes (Austin, Texas, November
1989), IEEE Society Press, New York, New York, November 1989,
pages 109--H115.
Stewman, John and Kevin Bowyer. "Creating the Perspective Pro-
jection Aspect Graph of Polyhedral Objects," in Proceedings of the
Second International Conference on Computer Vision, December
1988, pages 49?500.
Stroustrup, Bjarne. The C++ Prngramming Language, Addison-
Wesley Series in Computer Science, Addison-Wesley Publishing
Company, Reading, Massachusetts, 1986.
Sutherland, Ivan E., Robert F. Sproull, and R. A. Schumacker. "A
Characterization of Ten Hidden-Surface Algorithms," ACM Com-
puting Surveys, 6(1), March 1974.
Takagi, Atsushi, Hitoshi Takaoka, Tetsuya Oshima, and Yoshi-
nori Ogata. "Accurate Rendering Technique Based on Colori-
metric Conception," Proceedings of SIGGRAPH'90 (Dallas, Texas,
August ?10, 1990), in Computer Graphics, 24(4), August 1990,
pages 263272.
Tamminen, Markku and F. W. Jansen. "An Integrity Filter for
Recursive Subdivision Meshes," computer & graphics, 9(4), 1985,
pages 351--H364.
Tampieri, Filippo and Daniel Lischinski. "The Constant Radiosity
Assumption Syndrome," in Proceedings of the Second Eurograph-
ics Workshop on Rendering (Barcelona, Spain, May 13-15, 1991),
May 1991.
[TAMP9l]
183
[TANA91
[TARJ83i
[TELL92]
[THIB87]
[TRUM91]
Tanaka, Toshimitsu and Tokiichiro Takahashi. "Shading with Area
Light Sources," in Proceedings of Eurographics'91 (Vienna, Aus-
tria, September 16, 1991), September 1991, pages 235246.
Tarjan, Robert Endre. Data Structures and Network Algorithms, So-
ciety for Industrial and Applied Mathematics, Philadelphia, Penn-
sylvania, 1983.
Teller, Seth J. "Computing the Antipenumbra of an Area Light
Source," Proceedings of SIGGRAPH'92 (Chicago, Illinois, July 2?
31, 1992), in Computer Graphics, 26(4), July 1992, pages 139--H148.
Thibault, William C. and Bruce F. Naylor. "Set Operations on Poly-
hedra Using Binary Space Partitioning Trees," Proceedings of SIG-
CRAPH'87 (Anaheim, California, July 27--H31, 1987), in Computer
Graphics, 21(4), July 1987, pages 153--H162.
Trumbore, Ben, Wayne Lytle, and Donald P. Greenberg. "A
Testbed for Image Synthesis," in Proceedings of Eurograph-
ics'91 (Vienna, Austria, September 1--H6, 1991), September 1991,
pages 467?8O.
[VARG62] Varga, Richard 5. Matrix Iterative Analysis, Prentice-Hall, Engle-
wood Cliffs, New Jersey, 1962.
[vH87]
[vH89]
[vL91]
von Herzen, Brian and Alan H. Barr. "Accurate Triangulation of
Deformed, Intersecting Surfaces," Proceedings of SIGGRAPH'87
(Anaheim, California, July 27--H31, 1987), in Computer Graphics,
21(4), July 1987, pages 103--H110.
von Herzen, Brian. Applications of Surface Networks to Sampling
Problems in Computer Graphics, PhD dissertation, California Insti-
tute of Technology, Pasadena, California, 1989.
van Liere, R. "Divide and Conquer Radiosity," in Proceedings
of the Second Eurographics Workshop on Rendering (Barcelona,
Spain, May 1?15, 1991), May 1991.
Wallace, John R., Michael F. Cohen, and Donald P. Greenberg. "A
Two-Pass Solution to the Rendering Equation: A Synthesis of Ray
Tracing and Radiosity Methods," Proceedings of SIGGRAPH'87
(Anaheim, California, July 27--H31, 1987), in Computer Graphics,
21(4), July1987, pages 311--H320.
[WALL87]
184
[WALL89]
[WALL9Oa]
[WALL9Ob]
Wallace, John R., Kells A. Elmquist, and Eric A. Haines. "A Ray
Tracing Algorithm for Progressive Radiosity," Proceedings of SIG-
GRAPH'89 (Boston, Massachusetts, July 31--HAugust 4, 1989), in
Computer Graphics, 23(3), July 1989, pages 315--H324.
Wallace, John. "Trends in Radiosity for Image Synthesis," in Pro-
ceedings of Eurographics Workshop on Photosimulation, Realism
and Physics in Computer Graphics (Rennes, France, June 11--H12,
1990), June 1990, pages 1--H14.
Wallace, John R. "Radiosity and Ray Tracing: A comparison of
Shading Strategies," SIGGRAPH'90 Advanced Topics in Ray Trac-
ing Course Notes, August 1990.
[WALL92] Wallace, John R. "Introduction to Meshing Algorithms for Radios-
ity," SIGGRAPH'92 Radiosity Course Notes, July 1992.
[WANG91a]
[WANG9Ib]
[WARD88]
[WATT9O]
[WEGH84]
[WHIT8O]
Wang, Changyaw. "Direct Lighting Calculation by Monte Carlo
Integration," in Proceedings of the Second Eurographics Work-
shop on Rendering (Barcelona, Spain, May 13--H15,1991), May 1991.
Wanger, Leonard R. Perceiving Spatial Relationships in Computer
Generated Images, Master's thesis, Program of Computer Graphics,
Cornell University, Ithaca, New York, May 1991.
Ward, Gregory J., Francis M. Rubinstein, and Robert D. Clear. "A
Ray Tracing Solution for Diffuse Interreflection," Proceedings of
SIGGRAPH'88 (Atlanta, Georgia, August 1--H5, 1988), in Computer
Graphics, 22(4), August 1988, pages 85--H92.
Watt, Mark. "Light-Water Interaction using Backward Beam Trac-
ing," Proceedings of SIGGRAPH'90 (Dallas, Texas, August 61Q
1990), in Computer Graphics, 24(4), August 1990, pages 377--H386.
Weghorst, H., G. Hooper, and Donald P. Greenberg. "Improved
Computational Methods for Ray Tracing," ACM Transactions on
Graphics, 3(1), January 1984, pages 52--H69.
Whitted, Turner. "An Improved Illumination Model for Shaded
Display," Communications of the ACM, 23(6), June 1980, pages 34?
349.
Wolberg, George and Terrance E. Boult. "Separable Image Warp-
ing with Spatial Lookup Tables," Proceedings of SIGGRAPH'89
(Boston, Massachusetts, July 31--HAugust 4, 1989), in Computer
Graphics, 23(3), July 1989, pages 369--H378.
[WOLB89]
185
[WOLB9O] Wolberg, George. Digital Image Warping, lEFE Computer Society
Press, Los Alamitos, California, 1990.
[XU89]
[YEN56]
Xu, Hau, Qun-Shang Peng, and You-Dong Liang. "Accelerated
Radiosity Method for Complex Environments," in Proceedings of
Eurographics'89 (September 1989), September 1989, pages 51?I.
Yen, J. L. "On Nonuniform Sampling of Bandwidth-Limited
Signals," IRE Transactions on Circuit Theory, 3, December 1956,
pages 251--H257.
Zatz, Harold Robert Feldman. Galerkin Radiosity: A Higher Order
Solution Method for Global Illumination, Master's thesis, Program of
Computer Graphics, Cornell University, Ithaca, New York, August
1992.
[ZATZ92]
