Bounds and error estimates for radiosity.
We present a method for determining a posteriori bounds and estimates for local and total errors in radiosity solutions. The ability to obtain bounds and estimates for the total error is crucial fro reliably judging the acceptability of a solution. Realistic estimates of the local error improve the efficiency of adaptive radiosity algorithms, such as hierarchical radiosity, by indicating where adaptive refinement is necessary.First, we describe a hierarchical radiosity algorithm that computes conservative lower and upper bounds on the exact radiosity function, as well as on the approximate solution. These bounds account for the propagation of errors due to interreflections, and provide a conservative upper bound on the error. We also describe a non-conservative version of the same algorithm that is capable of computing tighter bounds, from which more realistic error estimates can be obtained. Finally, we derive an expression for the effect of a particular interaction on the total error. This yields a new error-driven refinement strategy for hierarchical radiosity, which is shown to be superior to brightness-weighted refinement.
INTRODUCTION
During the past decade, radiosity has become the method of choice
for simulating global illumination in environments consisting entirely
of Lambertian (ideal diffuse) reflectors and emitters. Global
illumination in such environments is governed by a Fredholm integral
equation of the second kind:
Z
k(x, y)B(y) dy, (1)
E-mail: {danix | bes | dpg}@graphics.cornell.edu
where B(x) is the total radiosity at point x, E(x) is the emission, and
ae(x) is the reflectivity 1 . The kernel of the integral operator
xy
describes the point-to-point transfer from y to x: ' x and ' y are the
angles between the surface normals and the line connecting x and
y, r xy is the distance between the two points, and V(x, y) is 1 if y is
visible to x and 0 otherwise.
Radiosity algorithms compute an approximate solution b B(x) to
equation (1) by projecting the continuous integral equation into a
finite dimensional function space. Most radiosity algorithms use
spaces of piecewise constant functions, although recently several
investigators described radiosity algorithms that use higher order
basis functions [12, 21, 22]. The projection of the continuous equation
results in a discrete system of n linear equations, where n is the
dimension of the finite dimensional space:
Since radiosity algorithms produce only approximate results, reliable
error bounds and realistic error estimates are crucial for assessing
the acceptability of a particular solution, as well as for automatic
adaptive refinement.
The theory of integral equations [1, 10] provides an a priori
error analysis of computational methods for Fredholm equations of
the second kind. This general framework is valuable for analyzing
convergencerates [10], and for boundingvarious components of the
error in terms of the corresponding operator norms [3]. However,
a priori error bounds are typically pessimistic and often difficult
to evaluate, requiring information about the operators or the exact
solution, which may not be available. In this paper we discuss a
posteriori methods that bound the combined effect of the various
error sources for a particular instance of the problem. Because the
analysis takes place after or during the computation of the solution,
such methods can provide tighter bounds more easily.
During the 1980's the issues of a posteriori error estimation and
automatic adaptive refinement received considerable attention in the
finite elements literature [4, 19]. Most of the work has related to
partial differential equations. More recently, results have become
available also for boundary integral equations andboundaryelement
methods [5], which are essentially similar to radiosity methods.
In the radiosity literature, however, very little has been written
regarding a posteriori error analysis. In particular, current radiosity
algorithms can provide the user with neither guaranteed bounds on
the total error any function norm k \Delta k, nor reliable
1 All of these quantities depend on the wavelength, but for simplicity we
restrict the discussion to a single wavelength throughout this paper.
estimates of the total error. The goal of this paper is to address these
deficiencies.
We begin with a brief account of past attempts to estimate the
error in radiosity approximations.
In Section 3 we describe an algorithm that computes two piece-wise
constant functions that are guaranteed to bound the exact solution
B(x) from below and from above. This can be done simultaneously
with the computation of the approximate solution b B(x).
These bounds pertain to both B(x) and b B(x), so they can be used
to compute a conservative upper bound on the error
Our algorithm resembles the general bracketing technique suggested
by Brown [6], but we utilize special properties of equation (1) to
produce tighter bounds.
In many instances,a conservative error bound is not required,or is
too pessimistic, and a good error estimate is preferable. In Section 4
we describe a non-conservative bounding algorithm that computes
tighter "bounds", which are no longer guaranteed to contain the
exact solution B(x) everywhere. These non-conservative bounds
are much cheaper to compute, and result in more realistic error
estimates, i.e., estimates which are closer to the actual errors.
Our final contribution is the derivation of an expression that
describes the effect that a particular interaction has on the total error
bound. This expression uses boundson importance and on radiosity
that can be computed using our algorithms. This gives rise to a new
error-driven refinement strategy for hierarchical radiosity, which is
shown to be superior to brightness-weighted refinement.
Cohen's adaptive subdivision algorithm [8] can be viewed as the
first attempt at a posteriori error estimation. In this algorithm,
elements are subdivided if their radiosity differs from that of their
neighbors by more than a certain threshold. Thus, the difference
between radiosities of adjacent elements is treated as an indicator of
the local error. Note that this heuristic may fail to identify elements
that should be subdivided.
A more conservative approach is described by Campbell [7], who
uses numerical optimization to find the minimum and the maximum
radiosity over each element. The difference between these extrema
is used as a criterion for adaptive subdivision.
Among the various radiosity methods,hierarchical radiosity (HR)
comes closest to bounding the errors in the solution: a bound is
computed on the error in each interaction (link) between patches.
Hanrahan et al. [13] compute an approximate upper bound on the
form-factor corresponding to the link, and use it to obtain an approximate
upper bound on the transferred energy.
Smits et al. [17] obtain an estimate of the error in each interaction
by point sampling the kernel function over the areas of the two
interacting patches. The difference between the maximum and the
minimum values thus obtained is taken to be the error in the form-
factor. This error estimate is less conservative than Hanrahan's, but
it is generally more accurate.
Given error estimates for all the links in an environment, one
can approximate the total error over each patch by summing up the
errors associated with all its links. However, this does not take into
account propagation of errors, which parallels the propagation of
light.
Several other error estimation techniques are surveyed by Cohen
and Wallace [9]. These techniques canbe classified into two groups:
(i) comparing the approximation to a higher order approximation,
and (ii) computing the residual
Z
dy.
Both of the above are computationally expensive and do not provide
foolproof subdivision criteria.
3 A RADIOSITY-BOUNDING ALGORITHM
For simplicity we start with a radiosity-bounding algorithm for full
matrix radiosity. Later in this section, we extend this algorithm to
work within a more efficient radiosity algorithm, namely hierarchical
radiosity.
The goal of the bounding algorithm is to compute two piecewise
constant functions B(x) and B(x) that bound the exact radiosity B(x)
from below and from above, respectively. In other words, given a
discretization of the environment into a set of n elements
we , such that
Roughly, this is achieved by replacing the coefficients of the discrete
linear system (2) by infima and supremabounding the corresponding
continuous functions over the areas of the elements
3.1 Computing Lower Bounds
the infima on the reflectivity and on the emis-
sion, respectively, over S i . Also, let F ij denote the infimum on the
(point-to-area) form-factor from a point on S i to
Z
and let K
Then the lower bounds vector
can be obtained by solving a linear system of n equations
using a standard method such as Jacobi or Gauss-Seidel iteration.
Proof. It is easy to see that the matrix I \Gamma K is strictly diagonally
dominant, from which it follows [11] that the spectral radius of K is
strictly less than 1. Thus, the inverse matrix
can be expressed as a Neumann series
The solution vector B can then be written as
Alternatively, B can be expressed as the limit of the series of vectors
B (k) that are defined by the recurrence relation
Note that the above corresponds to Jacobi iteration in matrix nota-
tion. It is easy to see by induction that B (k) is a lower bound on
the radiosity after k light bounces. Therefore, the limit B is a lower
bound on the total radiosity function B(x). 2
3.2 Computing Upper Bounds
Unfortunately, we cannot use the same straightforward approach for
computing the upper bounds
. Let K be defined in
the same way as K, except that suprema are taken instead of infima.
Then the matrix I \Gamma K is not necessarily diagonally dominant: the
entries of K are upper bounds on the form-factors and thus the sum
of the absolute values of the entries in a row of I \Gamma K can exceed 1.
In such a case, iterative solution methods such as Jacobi or Gauss-Seidel
could diverge. Furthermore, it is possible for I \Gamma K to be
singular.
We now describe an iterative algorithm that transforms K into a
strictly diagonally dominant matrix K   , such that the solution to the
linear system of equations
is a vector of upper bounds. The solution B is computed simultaneously
with the transformation of K into K   .
The algorithm starts by forming the matrix K. We then perform
Jacobi iterations using a modified matrix: because the sum of the
form-factors from a given element cannot exceed 1, we can zero
out elements in each row of K until the form-factors corresponding
to the remaining non-zero entries sum to at most 1. To ensure
that the modified matrix still yields an upper bound on the current
light bounce,we zero out those entries correspondingto the dimmest
elements. Since the brightness ordering of the elements may change
from one iteration to another, different entries on the matrix may be
zeroed out in each iteration. However, as the iterates converge to
the solution B, the order of the elements becomes fixed, at which
point we obtain K   .
More precisely, at each iteration we sort the current solution vector
in order of decreasing brightness, and permute the columns
of the matrix accordingly. Each entry in the solution vector is
updated by:
l
where l is the largest index in row i such that
It can be shown by induction that, if the iterates are sorted as described
is an upper bound on the radiosity after k inter-
reflections. Thus, if the iterates converge, they converge to a total
upper bound on the radiosity function B(x).
Proof. To prove convergence, we must examine the behavior of
each entry in the solution vector as the algorithm progresses. Because
all the entries in the matrix K are non-negative, if B (0) is set to
E, each entry in the solution vector increases at each iteration, i.e.,
for all i, B (k)
. Because we do not allow the form-factors in
a row to exceed 1, each entry of B (k) is bounded from above:
where Emax is the maximum emission and ae max the maximum reflectivity
in the environment. From elementary analysis, we know
that a monotonically increasing series that is bounded from above
must converge. Since the vectors B (k) are of a finite dimension, it
follows that the vector sequence converges. 2
GatherLowerBounds(node, lower)
foreach link 2 node.links do
lower += node.rho   link.Fmin   link.source.lower
end for
if IsLeaf(node) then
else
foreach child 2 node.children do
GatherLowerBounds(child, lower)
end for
GatherUpperBounds(node, contribList)
foreach link 2 node.links do
add the pair (link.Fmax, link.source.upper) to contribList
end for
if IsLeaf(node) then
foreach pair (ff, upper) 2 contribList do
ffSum += ff
node.upper+= node.rho   ff   upper
else
break
end for
else
foreach child 2 node.children do
GatherUpperBounds(child, contribList)
end for

Figure

1: Pseudocode for gathering lower and upper bounds
3.3 A Radiosity-Bounding Algorithm for HR
The algorithms in Sections 3.1 and 3.2 are presented in the context
of the full matrix radiosity algorithm. They involve O(n
to compute the bounds on the form-factors, and therefore are not
practical. However, they can be modified to work more efficiently
within the HR framework. This requires a few changes to the
standard HR algorithm. We need to compute and store upper and
lower bounds Fmin and Fmax on the form-factor associated with
each link. We can then gather, push, and pull radiosity bounds, in a
way similar to that of the patch radiosities themselves. See Figure 1.
Note that gathering, pushing, and pulling are all performed in a
single sweep of the hierarchy. Thus, the bounds are updated in place,
which makes this a Gauss-Seidel iteration. This iteration converges
faster than the Jacobi iteration used in previous HR algorithms. The
speedup applies not only to the computation of bounds, but to the
computation of the radiosities as well.
Gathering upper bounds is a bit more involved than gathering
lower bounds, as can be seen from the pseudocode in Figure 1.
GatherUpperBounds essentially pushes links down to the leaves,
where they are sorted according to the brightnesses of the corresponding
source nodes. The upper bound of the leaf node is then
updated using equation (3).
One undesirable feature of this algorithm is the sorting of the
links, which implies that O(k log work must be done to update the
upper bound at each leaf node, where k is the number of the contributing
links. However, complete sorting can be avoided in most
cases. Sorting is only necessary when the sum of the contributing
form-factors exceeds 1. Even then, only a few of the dimmest
sources must be removed from the list. This can be done efficiently
using a heap data structure with a DeleteMin operation. Tarjan [20]
provides a detailed description of such data structures.
3.4 Bounding Form-Factors
Thus far we have assumed that the bounds F ij and F ij are available
to us. How do we obtain them? We know of no analytical
way of computing such bounds. Thus, we must resort to numerical
optimization to find the minimum and the maximum values
of the point-to-polygon form-factor from points on S i to S j . Such
numerical methods are discussed in detail by Campbell [7]. These
methods generally require the computation of form-factor gradients
in the presence of occluders, which can be done either using finite
differences or analytically, as described by Arvo [2].
In our current implementation we estimate form-factor bounds
by evaluating the analytical point-to-polygon form-factor [18] at the
center and at the vertices of each receiving patch. In the partially
occluded case, the analytical formula is applied only to the visible
parts of the source [16]. The upper and lower bounds on the form-factor
are then set to the maximum and the minimum of these values.
While these bounds are not conservative, as elements decrease in
size the point-to-area form-factor function becomesmonotonic over
most elements, and the heuristic yields accurate bounds. This is
particularly sowhena discontinuity-driven subdivision strategy [15]
is employed.
The plots in Figures 2, 3, and 4 show the piecewise-constant bounding
functions computed by the hierarchical bounding algorithm for
three simple environments shown in Figure 5. The exact radiosity
function is also plotted for each of the three cases. The functions are
plotted along the dotted lines across the floor of each environment.
The bounds in Figure 2 were found in a single iteration, since
there are no interreflections in this environment, and errors do not
propagate. Because of this, the bounds are very tight: in fact, tighter
piecewise constant bounds are only possible if more elements are
used.

Figures

3 and 4 show two cases with interreflections, which results
in looser bounds due to error propagation. As illustrated in the
figures, the boundscomputedby the algorithm become tighter as the
accuracy of the solution increases. Note that the radiosity function
and the bounds are not symmetric, because in the corresponding
environments the reflectivity of the left wall is higher than that of
the right wall.
Having obtained bounds on the radiosity of each element, we can
obtain an upper bound on the local error there. Assuming that the0.50.60.7Radiosity
bounds
exact

Figure

2: Upper and lower bounds for the configuration shown in
Figure 5a. The smooth curve in the middle is the exact radiosity
coarse
fine
exact

Figure

3: Nested lower and upper bounds corresponding to two HR
solutions for the configuration in Figure 5b. One solution is coarse
and the other is fine. The smooth curve in the middle is the exact
radiosity function.0.20.611.4
Radiosity
coarse
fine
exact

Figure

4: Nested lower and upper bounds corresponding to two HR
solutions for the configuration in Figure 5c. One solution is coarse
and the other is fine. The smooth curve in the middle is the exact
radiosity function.
(a) (b) (c)

Figure

5: Three simple test cases: (a) direct illumination only; (b) inter-reflections without occlusion; (c) inter-reflections with occlusion. All
the surfaces are grey with reflectivities between 0.3 and 0.7.
approximate radiosity B lies halfway between B i and B , the errors
in the L1 , L 1 , and L 2 norms are:
R
R
However, these error bounds represent a worst case scenario, and
do not give an accurate error estimate even when the bounds on
radiosity are as tight as in Figure 2. More realistic estimates are
obtained if we assume that the values of B(x) over S i are uniformly
distributed between B i and B i . Under this assumption the L 1 and L 2
norm errors are estimated as:
We compute local error estimates at all the leaves of the hierarchy.
These estimates are then pulled up to the roots of the hierarchies,
and combined to obtain the estimate for the total error. The error
at a parent node is the sum of the children's errors in the L 1 norm,
and the maximum over the children's errors in L1 norm. For the
norm, the squares of the children's errors are summed and the
square root of the sum is taken.
4.1 Obtaining Tighter Non-Conservative Bounds
The conservative algorithm produces bounds which are quite loose
even for simple environments, resulting in pessimistic error esti-
mates. For most practical purposes good non-conservative error
estimates could prove more useful than conservative error bounds.
To quote Delves and Mohamed [10]:
Given two error estimates of equal realism and equal cost and
given that one is also a bound, we would clearly prefer the
bound. In practice bounds usually prove expensive and pes-
simistic. Most users will then usually prefer a cheaper and
more realistic error estimate, recognizing and accepting that
estimates can sometimes err on the side of optimism.
We now present a way for computing tighter, although not conser-
vative, bounds on radiosity. Error estimates computed from these
radiosity bounds are not guaranteed error bounds, but they are much
closer to the true errors. To compute these tighter bounds we approximate
This is equivalent to first setting both B and B to b B and then performing
a single gathering iteration over the environment using the same
GatherLowerBounds and GatherUpperBounds routines as before.
Because only one iteration is performed, error propagation is no
longer accounted for, but the computation is significantly simpli-
fied. More conservative estimates can be obtained by using two or
more iterations, yielding a trade-off between speed and conservativity

Two additional changes to the algorithm have been used to produce
even tighter bounds. First, instead of bounding each form-factor
separately and summing the bounds, we directly bound the
radiosity contributed to each node through all its links: the contribution
is evaluated at the center and at the vertices of the patch, and the
extrema values are used as bounds. Second, we set node.lower and
node.upper to the area-weighted averages of the children bounds,
rather than to the minimum and the maximum, respectively.
4.2 Results
To test the quality of our error estimates we computed an extremely
fine radiosity solution for each of the three environmentsin Figure 5.
The floor polygon in eachenvironmentwassampledon a 400 by 400
grid to serve as a reference. We then computed a series of radiosity
solutions of various degrees of accuracy for each environment. The
floor polygon was sampled as before, and the difference between
each solution and the reference solution was computed. This difference
is referred to as the measured error. During each solution, two
error estimates were computed: one from the conservative bounds,
and the other from the tighter non-conservative bounds. In both
cases equations (4), (7), and (8) were used to estimate the errors.

Figure

6 shows plots of the error estimates together with the
measured error for each of the three test cases. The left column
corresponds to L1 norm and the right column to L 1 norm. Error
plots for L 2 norm are not given, since they look almost exactly the
same as the L 1 plots for these test cases.
In the top pair of plots (corresponding to Figure 5a) the conservative
and the non-conservative bounding algorithms yield the
same bounds, and therefore the corresponding error estimates are
always the same. For the L1 norm these estimates follow tightly
the measured error, and for the L 1 norm they give realistic upper
bounds.
In the remaining two cases there is a significant difference between
the two types of error estimates. In the L1 norm, the non-conservative
estimate does not always yield a bound on the error
(see lower left plot). In the L 1 norm, both estimates bound the
measured error from above, but the conservative error estimates are
much more pessimistic than their non-conservative counterparts.
h (element size)
conservative
non-conservative
h (element size)
conservative
non-conservative
h (element size)
conservative
non-conservative
h (element size)
conservative
non-conservative
h (element size)
conservative
non-conservative
h (element size)
conservative
non-conservative
measured

Figure

Estimated and measured errors as a function of the element size for the three environments shown in Figure 5. The top row
corresponds to Figure 5a, the middle row to 5b, and the bottom row to 5c. L1 errors are shown on the left and L 1 errors on the right.
Both types of estimates become closer to the measured error as the
accuracy of the solutions increases. Both the estimated and the
measured errors vary linearly in the element size h, as expected for
piecewise-constant approximations.
The original hierarchical radiosity algorithm [13] adaptively improves
the accuracy of a solution by estimating the error associated
with each link. All links for which this error estimate exceeds a
given threshold ffl are refined and the solution is recomputed. More
specifically, Hanrahan et al. suggest to refine all those links for
which the estimated upper bound on the transferred energy exceeds
ffl. A link is refined by subdividing the node which has the greater
form-factor, as seen from the center of the other node. This strategy
is referred to as brightness-weighted refinement, as it gives priority
to links transferring energy from bright sources.
We make two observations with respect to brightness-weighted
refinement: (i) The links carrying the largest amount of energy are
not necessarily the ones with the greatest errors in the transferred
energy. (ii) Due to propagation of errors, refining the link between
nodes i and j reduces not only the error over patch i, but may also
reduce errors on other patches receiving light from it. Thus, to make
the most benefit of each link refinement, priority should be given to
links with the greatest effect on the total solution error.
The first observation suggests using the
estimate of the link form-factor error, rather than using an upper
bound on the form-factor. The error in the transferred energy is
then given by
To address the second observation, we derive an expression that
relates DK ij to the overall value of which is an upper
bound on the overall error in L 1 norm. This is done by using the
concept of importance.
5.1 Importance and Linear Functionals
Recall that importance, which we shall denote by Z, is the solution
to the adjoint of the transport equation [14, 17]:
where M T is the adjoint of the real matrix M. In the context of
radiosity, is the receiver vector, which is dual to the
emission vector E.
Any linear functional f (B) of the radiosities can be expressed as
an inner product R T B, where R i gives the contribution of B i to the
value of f (B). It can be shown [17] that
Thus, the i-th element of the importance vector gives the contribution
made by a unit of emission at patch i to the value of R T B.
Note that because all the entries in the vector are non-
negative, the value of its L 1 norm is given by a linear functional:
is the vector of patch areas. We can now
derive the main result of this section:
.
This expression relates to the link errors DK and to errors
in the emission. The term Z T DKB is a sum over all links. The values
summed are the link form-factor error DK ij weighted by the lower
bound on the importance of the receiver Z i and the upper bound on
the radiosity of the source B j . Intuitively, is an upper bound
on the incorrect radiosity that patch i receives from patch j, while
Z i gives the contribution of this incorrect radiosity on the total error
bound. The Z T E) term gives the effect non-constant emission
has on the total error bound.
Proof. Let K. The system for lower
boundson the radiosities and its adjoint equation can now be written
as
and the upper bounds system can be written as Noting
that
from which it follows that
in terms of DK as
follows:
Relative
brightness-weighted
error-driven

Figure

7: Relative L 1 error as a function of the number of links:
brightness-weighted vs. error-driven refinement.
Substituting equation (10) into the above expression we have
E).Equation (9) is an exact expression relating errors in the transferred
energy to the value of any linear functional R T In
particular this expression applies to the view-dependent functionals
used by Smits et al. [17]. Thus, the expression that Smits et
al. used to assess the view-dependent error can be viewed as an
approximation to our more conservative bound.
5.2 A New Refinement Strategy
Equation suggests a new error-driven refinement strategy for
HR. In this strategy we refine all links for which Z
since these links have the greatest effect on the total error bound.
In order to evaluate this expression for each link, the routine
GatherLowerBounds must be modified to solve for the lower bounds
on importance in addition to lower bounds on radiosity. These
changes are straightforward, and they do not significantly increase
the amount of work done in the routine.
When we refine a link between a receiver node i and a source
node j, we need to decide which node to subdivide to obtain the
greatest reduction in the error. We use the following heuristic: node
i is subdivided if
and node j is subdivided otherwise. Intuitively, this means that we
subdivide the receiver if the error due to the non-constant form-factor
is greater than the error due to non-constant radiosity on the
source.
5.3 Results
We implemented both brightness-weighted and error-driven refinement
strategies and tested them on the simple environment shown in
Figure 5c. The plot in Figure 7 shows how the measured relative
error on the floor decreasesas a function of the number of links. This
plot demonstrates clearly that error-driven refinement is much more
efficient: it is capable of achieving much smaller errors using only
a fraction of the links required by brightness-weighted refinement.
Similar results were observed on several other test cases.
6 CONCLUSIONS
We have described a hierarchical radiosity algorithm that computes
conservative bounds on the exact radiosity function simultaneously
with the computation of the approximate solution. In contrast to
previous attempts to estimate the errors in radiosity, our algorithm
properly accounts for the propagation of errors due to interreflec-
tions, and provides conservative upper boundson the error. We have
also presented a non-conservative version of this algorithm that provides
tighter bounds and more realistic error estimates. Finally,
we have derived a new error-driven adaptive refinement strategy
for hierarchical radiosity that significantly outperforms brightness-
weighted refinement.
We hope that the results described in this paper will prove useful
not only for reliable image synthesis, but also for other fields in
which similar integral equations arise. Examples include radiative
heat transfer, illumination engineering, and remote sensing. In all
these fields numerical accuracy of the results is even more crucial
than in image synthesis.
Much further work is needed on a posteriori error analysis for
radiosity. For instance, we need techniques for efficiently computing
reliable tight bounds on form-factors, perhaps using tools such
as irradiance Jacobians [2]. Also, we would like to extend the techniques
described in this paper to higher order basis functions. While
we are able to bound the exact radiosity function, our bounds will
not necessarily apply to higher order approximations, as they may
exhibit oscillations [22]. Even if the approximation is contained
within our bounds, equations (7) and (8) will probably produce pessimistic
error estimates, since they assume that the approximation
is piecewise-constant.
Prior to rendering an image, the radiosity solution is usually
projected onto a new set of basis functions, typically piecewise-linear
or piecewise-quadratic [15]. We need to analyze the effect of
such projections on the error in order to extend our error estimates
from solutions to images.
Finally, we would like to emphasizethe need for the development
of perceptual error metrics for image synthesis. Fast convergencein
a quantitative error metric, such as the L 1 norm, does not necessarily
imply fast convergence of the resulting images, as perceived by a
human observer. Thus, we expect that a perceptual error-driven
refinement strategy would be more appropriate for image synthesis.

ACKNOWLEDGEMENTS

Special thanks go to Jim Arvo: without his encouragement and advice this
paper would not have been written. We would also like to thank Richard
Lobb and Charlie Van Loan for helpful discussions and for reviewing the
manuscript. The comments provided by the anonymous reviewers have
helped greatly to improve this paper. This work was supported by the
NSF grant, "Interactive Computer Graphics Input and Display Techniques"
(CCR-8617880), by the NSF/DARPA Science and Technology Center for
Computer Graphics and Scientific Visualization (ASC-8920219), and by
generous donations of equipment from Hewlett-Packard.



--R

Convergence and Error Bounds for Approximate Solutions of Integral and Operator Equations.
The Irradiance Jacobian for Partially Occluded Polyhedral Sources.
A Frame-work for the Analysis of Error in Global Illumination Algorithms

Adaptive Finite and Boundary Element Methods
Upper and Lower Bounds for Solutions of Integral Equations.

BROCK An Efficient Radiosity Approach for

Computational Methods for Integral Equations
LOAN. Matrix Computa- tions
A Rapid Hierarchical Radiosity Algorithm.
The Adjoint Function: The Physical Basis of Variational and Perturbation Theory in Transport and Diffusion Problems
GREENBERG. Combining Hierarchical Radiosity and Discontinuity Meshing.

SALESIN. An Importance-Driven Radiosity Algorithm
A New and Simpler Formulation for Radiative Angle Factors.


MAX. Radiosity Algorithms Using Higher Order Finite Elements.

--TR
Data structures and network algorithms
Computational methods for integral equations
A rapid hierarchical radiosity algorithm
An importance-driven radiosity algorithm
Modeling global diffuse illumination for image synthesis
Radiosity and realistic image synthesis
Combining hierarchical radiosity and discontinuity meshing
Radiosity algorithms using higher order finite element methods
Galerkin radiosity
Wavelet radiosity
A framework for the analysis of error in global illumination algorithms
The irradiance Jacobian for partially occluded polyhedral sources

--CTR
Robert Garmann , Heinrich Mller, On the computational complexity of hierarchical radiosity, Computer science in perspective, Springer-Verlag New York, Inc., New York, NY,
Fraois Sillion , George Drettakis, Feature-based control of visibility error: a multi-resolution clustering algorithm for global illumination, Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, p.145-152, September 1995
M. Pellegrini, Monte Carlo approximation of form factors with error bounded a priori, Proceedings of the eleventh annual symposium on Computational geometry, p.287-296, June 05-07, 1995, Vancouver, British Columbia, Canada
James Arvo, The irradiance Jacobian for partially occluded polyhedral sources, Proceedings of the 21st annual conference on Computer graphics and interactive techniques, p.343-350, July 1994
Hualin Zhou , Min Chen , Mike F. Webster, Comparative evaluation of visualization and experimental results using image comparison metrics, Proceedings of the conference on Visualization '02, October 27-November 01, 2002, Boston, Massachusetts
Wu Enhua , Wang Wenping, Interleaving radiosity, Journal of Computer Science and Technology, v.17 n.1, p.1-8, January 2002
Wolfgang Heidrich , Philipp Slusallek , Hans-Peter Seidel, Sampling procedural shaders using affine arithmetic, ACM Transactions on Graphics (TOG), v.17 n.3, p.158-176, July 1998
Frdo Durand , George Drettakis , Claude Puech, The visibility skeleton: a powerful and efficient multi-purpose global visibility tool, Proceedings of the 24th annual conference on Computer graphics and interactive techniques, p.89-100, August 1997
Franois X. Sillion, A Unified Hierarchical Algorithm for Global Illumination with Scattering Volumes and Object Clusters, IEEE Transactions on Visualization and Computer Graphics, v.1 n.3, p.240-254, September 1995
Volevich , Karol Myszkowski , Andrei Khodulev , Edward A. Kopylov, Using the visual differences predictor to improve performance of progressive global illumination computation, ACM Transactions on Graphics (TOG), v.19 n.2, p.122-161, April 2000
Cyril Soler , F. X. Sillion, Texture-based visibility for efficient lighting simulation, ACM Transactions on Graphics (TOG), v.19 n.4, p.302-342, Oct. 2000
Doug L. James , Dinesh K. Pai, ArtDefo: accurate real time deformable objects, Proceedings of the 26th annual conference on Computer graphics and interactive techniques, p.65-72, July 1999
George Drettakis , Eugene L. Fiume, Structured Penumbral Irradiance Computation, IEEE Transactions on Visualization and Computer Graphics, v.2 n.4, p.299-312, December 1996
Frdo Durand , George Drettakis , Claude Puech, Fast and accurate hierarchical radiosity using global visibility, ACM Transactions on Graphics (TOG), v.18 n.2, p.128-170, April 1999
Kavita Bala , Julie Dorsey , Seth Teller, Radiance interpolants for accelerated bounded-error ray tracing, ACM Transactions on Graphics (TOG), v.18 n.3, p.213-256, July 1999
adaptive refinement;a posteriori error bounds and estimates;importance;hierarchical radiosity;global illumination