Cohen M.F., Wallace J.R. - Radiosity and realistic image synthesis (1995)(en)

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION

3.4 APPROXIMATING RADIOSITY ACROSS A SURFACE

at the associated node. The radiosity function, B(x), is thus approximated by

ˆ (x), where

B

		ˆ	n
	B(x) ≈		åBi Ni (x)	(3.2)
		B(x) =

i =1

Evaluating ˆ (x) at a particular point x involves summing only those basis

B j functions with their support over the point.

There are many possible basis functions. Low-order polynomials are the most common, including constant, linear, quadratic, and cubic functions. The constant (or “box”) basis, which is often used in radiosity, is defined by

Ni	ì0	if x is outside element
	(x) = í	if x is inside element	(3.3)
	î1

The number and placement of nodes within elements will depend on the order of the polynomial. In the case of constant elements, defined by constant basis functions, a single node located in the center of the element is commonly used.

The radiosity function across linear elements is defined by nodes and associated basis functions at the boundaries or corners (in two dimensions) of the element. A simple example of function approximation using linear bases is shown in Figure 3.4 for a function of a single variable. In this example, a linear basis function is associated with each node, and the nodes are located at element boundaries. The linear (or “hat”) basis function in one dimension is defined by

	ì	x – xi –1
	ï	xi		– xi –1
	ï xi				1 – x
Ni	( x ) = í			+
		x			– x
	ï		i	+1 i
	ï0
	î

for	xi –1 < x < xi
for	xi	< x < xi +1	(3.4)

otherwise

In Figure 3.4 the nodal values, Bi , are determined by simply evaluating the function to be approximated at the node locations. In general, the function to be approximated is not known, thus the coefficient values are computed to minimize an error estimate or provide an overall best fit according to some criterion.

The generalization of finite element approximations to functions defined over a two-dimensional domain, such as the radiosity function, is straightforward (see Figure 3.5). The elements are typically simple shapes, usually triangles and/or convex quadrilaterals. Nodes for constant, linear, and quadratic elements are shown in Figure 3.6. The approximation of a function for a single, two-dimensional element using constant and linear basis functions is shown in Figure 3.7. Note that the tensor product of the two one dimensional linear basis functions results in a bilinear basis in two dimensions and is thus curved along the diagonals.

Radiosity and Realistic Image Synthesis	49
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION 3.4 APPROXIMATING RADIOSITY ACROSS A SURFACE

Figure 3.4: Finite element approximation of a function using linear basis functions.

Radiosity and Realistic Image Synthesis	50
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION 3.4 APPROXIMATING RADIOSITY ACROSS A SURFACE

Figure 3.5: Finite element approximation of the radiosity function.

Figure 3.6: Dividing the surface into elements with nodes.

Radiosity and Realistic Image Synthesis	51
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION

3.4 APPROXIMATING RADIOSITY ACROSS A SURFACE

In two dimensions, the linear basis function, Ni (x), has the properties:

ì1 at node i

ï

ï0...1 within adjacent elements

Ni (x) = í (3.5)

ï0 at all other node points

ïî0 outside adjacent elements

This definition ensures that the basis function associated with a particular node is zero outside of the elements adjacent to the node, as well as at any other nodes within those elements.

The variation of the linear basis function from a value of 1 at the node to 0 across the element depends on the parametric mapping of a generic basis defined on a standard element to the triangular or quadrilateral elements. In the case of triangles, barycentric coordinates are often used to define the 2D linear basis. In the case of quadrilaterals a tensor product of two one-dimensional linear basis functions defines one quarter of a bilinear basis within each element. For example, the value of the hat function is determined by the interpolation, Ni (x) = (1 — u)(1 — v), where (u, v) is the (0, 0) to (1, 1) parametric coordinate of the point x within the element with node i at the parametric origin (see Figure 3.9). Parametric mappings are discussed briefly in section 3.8 and in more detail in most finite element texts [273].

Elements constructed using other types of basis functions are also possible. For example, cubic Hermite elements provide continuity of derivatives across element boundaries by using basis functions that interpolate function derivatives at the nodes as well as values. Given the radiosity values and their parametric derivatives at the four corners of a quadrilateral, one can derive a bicubic variation of radiosity across the element. Hybrid sets of basis functions with constant, linear, and higher-order functions can also be employed. Other polynomial bases (e.g., Jacobi polynomial sets [270]) and non-polynomial functions are other possible choices. For a given approximation accuracy, there is typically a tradeoff between the number of elements required and the complexity of the basis functions.

To date, radiosity implementations have used constant basis functions almost exclusively. Linear functions have been explored in one dimension by Heckbert [122]. Max and Allison [163] have implemented radiosity using linear elements. Lischinski et al. [154] have used quadratic elements in a radiosity solution. For radiosity, as will be discussed at length in Chapter 9, elements of linear or higher order have more commonly been used during the rendering stage [203].

Radiosity and Realistic Image Synthesis	52
Edited by Michael F. Cohen and John R. Wallace

ˆ

B (x) to agree

53

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION

3.5 ERROR METRICS

Figure 3.7: Basis functions.

3.5 Error Metrics

A set of discrete elements, with a total of n nodal values at n nodal points, xi , and n basis functions, Ni , defines a finite dimension function space (see box on page 45). This finite function space has n degrees of freedom, the n nodal

values, B ,thus ˆ (x) is everywhere defined by the finite set of the n coefficients,

i B

Bi.

Ideally, one would like the approximate radiosity solution

Radiosity and Realistic Image Synthesis

Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION

3.5 ERROR METRICS

with the actual solution B(x) everywhere. This is not possible in general since

ˆ (x) is restricted to the finite subspace of functions representable as the linear

B

sum of the selected basis functions. Given the choice of basis functions, the goal is then to find the nodal values that generate a solution minimizing the error according to some measure or metric.

The actual error, ε(x), is the difference between the approximate and exact solutions, that is,

ˆ	(3.6)
ε(x) = B(x) – B(x)

Unfortunately, it is impossible to determine ε(x) directly since the exact solution is not known.

An alternative characterization of the error is provided by the residual. If

ˆ (x) is substituted for both occurrences of B in the original radiosity equation

B

ˆ	ˆ	(3.7)
B(x) = E(x) + ρ(x) òs	B(x´) G(x, x´) dA´

then the residual function is given by the difference of the left and right hand ides

ˆ	ˆ	(3.8)
r(x) = B(x) – E(x) –ρ(x) òs	B(x´) G(x, x´) dA´

An exact solution will make r(x) zero everywhere, but this is unobtainable since,

as we have seen, ˆ is restricted to the subspace of functions realizable using the

B

basis functions. Although ˆ is restricted to a finite function space, the residual,

B

r, is not. However, a similar approach to that taken in approximating B can be taken to project the residual into a finite function space. In this way, rather than seek a solution technique that makes r(x) zero everywhere, the goal will instead be to minimize r(x) according to a finite dimensional error metric.

The general approach is to choose n independent weighting functions, Wi(x), each having local support, much like the earlier basis functions. The norm or size of the residual is then approximated as a finite sum:3

n
r(x) = å < r(x), Wi (x) >	(3.9)

i =1

This residual norm can be driven to zero by finding radiosity values such that each of the n terms < r(x), Wi(x) > is equal to zero. As will be derived below,

3The notation < f(x), g(x) > indicates the inner product of the two functions. The inner product of two functions is analogous to the dot product of two vectors and is defined to be

< f(x), g(x) > = òs f(x) g(x) dx

Radiosity and Realistic Image Synthesis	54
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION

3.5 ERROR METRICS

If < f(x), g(x) > = 0, then f(x) and g(x) are orthogonal.

setting each of these terms to zero results in a linear equation. The set of n such linear equations (one for each weighting function) of the n nodal values, Bi , can be solved simultaneously as discussed in Chapter 5.

This approach defines a general class of methods called weighted residual methods. Two different choices for the weighting functions Wi(x) lead to two contrasting formulations, point collocation and the Galerkin method, which will be discussed in the following sections.

3.5.1 Point Collocation

The simplest set of weighting functions are delta functions:

Wi(x) = δ(x –xi )

(3.10)

which are zero unless x is coincident with the node at xi .

Using these weighting functions, the norm of the residual defined by equations 3.8 and 3.9 is minimized (i.e., zero) when r(x) is exactly zero at all the node points, xi (i.e., r(xi ) = 0, "i). This clearly differs from requiring r(x) to be zero everywhere, since r(x) is free to vary away from zero in between the nodes. However, as the number of nodes increases this difference diminishes. In general, the points at which the residual function is minimized are selected to be located at the nodes used in approximating the radiosity function. This technique is known as point collocation.

If there are n nodes, there is now a finite number of conditions that must be met. These are captured by the n simultaneous linear equations, one for each node i located at location xi :

ˆ	ˆ	"i	(3.11)
B (xi) –E(xi) –ρ(xi)	òs B (x´) G(xi, x´) dA´ = 0,

Note that the only change from equation 3.8 to equation 3.11 is that equation 3.11 is defined only at the n node locations, xi , as opposed to the entire surface

domain, S. Expanding ˆ using equation 3.2 gives

B

n	Bj Nj (xi ) –E(xi ) –ρ(xi) òs	n
å		å Bj Nj (x) G(xi, x´) dA´ = 0 (3.12)
j =1		j =1

Grouping terms and factoring Bj, which is independent of x´, out of the integral leaves

é	n	é	´ùù	0	(3.13)
ê	å Bj	êN j (xi ) – ρ(xi ) òS N j (x¢ ) G(xi	, x¢ ) dA'úú – E( xi ) =
ê j =1		ê	úú
ë		ë	ûû
Equation 3.13 can be restated simply as the set of n linear equations
Radiosity and Realistic Image Synthesis					55
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION 3.5 ERROR METRICS

n
å Bj Kij	–E(xi ) = 0	(3.14)
j =1
or in matrix/vector form as
K B = E		(3.15)
The coefficients Kij are given by
Kij = Nj(xi ) –ρ(xi )	òs Nj(x´) G(xi, x´) dA´	(3.16)

These coefficients are independent of the nodal radiosities and are defined solely by the geometry and material properties of the environment. The evaluation of these coefficients will be the topic of the following chapter. However, it should be noted here that for any particular term, Kij, the integral only needs to be evaluated over the nonzero support of the jth basis function. Also, the term Nj(xi ) is simply the Kronecker delta δij , (i.e., one if i = j and zero otherwise), for constant and linear basis functions with local support limited to adjacent elements. Equation 3.14 completes the objective of deriving a linear system of equations that can be solved for the unknown nodal radiosities, Bi .

3.5.2 Galerkin Form of Weighted Residuals

The point collocation formulation defines the error metric to be zero when the approximate residual is zero at the node points only. An alternative approach to deriving a linear system of radiosity equations defines the approximate residual to be zero when a set of weighted integrals of the residual are zero.

The weighting functions in this case are selected to have local support much like the basis functions used to approximate the radiosity function. In this case the weighted residual method seeks a solution for which there is a weighted “average” zero residual over each small region of the domain. If n weighting functions, Wi(x), are defined, a solution is found if

< Wi(x), r(x) > = òs Wi(x) r(x) dA = 0, ι

or expanding r(x) above with equation 3.8

òs Wi(x) ˆ (x) dA –

0 = B

òsWi(x) E(x) dA – òs Wi(x) ρ(x) òs ˆ (x´) G(x, x´) dA´ dA

B

Radiosity and Realistic Image Synthesis

Edited by Michael F. Cohen and John R. Wallace

(3.17)

(3.18)

56

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION 3.6 CONSTANT ELEMENT RADIOSITIES

The Galerkin formulation selects the same basis functions used to approximate the radiosity function as the weighting functions, (i.e., Wi(x) = Ni(x)), thus, equation 3.18 becomes,

			ˆ
0 = òs Ni(x) B (x) dA –
				ˆ	(3.19)
	òs Ni(x) E(x) dA – òs Ni(x) ρ(x) òs B (x´) G(x, x´) dA’ dA
Finally, expanding and grouping terms as before results in:
é	n	é	òS Ni (x) N j (x) dA –	òS Ni (x) ρ(x) òS N j (x¢ ) G(x , x	ù
ê	å Bj	ê			¢ ) dA¢ dAú
ê j =1		ê			ú
ë		ë			û
			–òs E(x) Ni(x) dA = 0		(3.20)

The unknowns, Bj , have been isolated in this expression. There are n such expressions, one for each node i and, as with point collocation, the linear system of equations can be expressed as the matrix equation

K B =	E	(3.21)
The entries of K are given by
Kij = òs Ni(x) Nj(x) dA – òs Ni(x) ρ(x)	òs Nj(x) G(x, x´) dA´dA	(3.22)
and the entries in E by
Ei = òs E(x) Ni(x) dA		(3.23)

The Galerkin method has been explored for one-dimensional “flatland” radiosity by Heckbert [122, 123].

3.6 Constant Element Radiosities

The previous two sections have taken two different approaches to achieving the same basic objective, a linear system of equations for the discrete nodal values. However, the relationship of these equations to the physical problem of global illumination may still be somewhat obscure. The complexity in the

Radiosity and Realistic Image Synthesis	57
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 3. DISCRETIZING THE RADIOSITY EQUATION 3.6 CONSTANT ELEMENT RADIOSITIES

above formulations is reduced if constant basis functions are used, with the reflectivity and emission of each element assumed to be constant. This set of assumptions has been made in most radiosity implementations to date and provides the clearest physical intuition.

For constant basis functions, Ni , the Galerkin formulation can be rewritten by performing the integration for each term in equation 3.20. First, using the fact that the box basis functions have values of only 11 and 0, and do not overlap with one another,

òs	Ni(x) Nj(x) dA = íìAi	if i = j ýü	= δij Ai	(3.24)
	î 0	otherwiseþ

where δij is the Kronecker delta, and Ai is the area of element i. Note that the integration only has to be performed across elements i and j since the corresponding basis functions are zero elsewhere. Likewise,

òs E(x) Ni(x) dA = Ei Ai

(3.25)

where Ei is the area average emission value for element i.

Finally, since the basis functions have unit value within the element, the basis functions themselves can be dropped and it is possible to integrate explicitly over element areas. Thus, assuming the reflectivity ρ(x) is a constant ρi over Ai ,

òs Ni(x) r(x) òs Nj(x´) G(x, x´) dA´ dA = ρi òAi òAj	G(x, x´) dAj dAi
	(3.26)

Making these substitutions into equation 3.20 results in, for all i:

é n	é	ρi	òAi	òAj		ùù
êå Bj	êδij Ai –				G (x, x¢ ) dAj	dAi úú – Ei Ai = 0	(3.27)
ê j =1	ê					úú
ë	ë					ûû

Dividing through by Ai and moving the emission term to the right side gives

n

é

1

òAi

òAj G (x,

ù

å Bj êδij –

ρi

x¢ ) dAj

dAi ú = Ei

(3.28)

A

j =1

ê

i

ú

ë

û

or

é

n

[

]

ù

(3.29)

ê

å

Bj

δij

– ρi Fij

ú = Ei

ê j =1

ú

ë

û

Radiosity and Realistic Image Synthesis	58
Edited by Michael F. Cohen and John R. Wallace