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

CHAPTER 4. THE FORM FACTOR

4.3. THREE FORMULATIONS OF THE FORM FACTOR

Figure 4.2: Differential solid angle.

4.3 Three Formulations of the Form Factor

The full form factor integral for constant bases, equation 4.4, involves performing an integration over the support of two basis functions (i.e., the two element areas). Three different formulations of this form factor integral between elements can be derived. Each is used in the algorithms that will be described in the following sections.

The first formulation results from inserting the differential form factor (equation 4.6) into the form factor integral, equation 4.4. Thus, the element-to-element form factor is the double area integral, with geometric kernel, Gij ,

Fij =	1	òAi òAj	cos θi cos θ j	Vij dAj dAi	(4.7)
	A		π r2
	i

where Vij is the visibility term between differential areas dAi and dAj .

The second formulation results from replacing the inner integration over area with an equivalent integration over the hemisphere around dAi . The differential

solid angle (see Figure 4.2) dwj from dAi to dAj is
dω		=	cos θ j	dA		(4.8)
	j				j
			r2

Thus, the inner integral of equation 4.7 can be rewritten over the hemisphere of directions, Ω, and over dAi , resulting in the area-hemisphere integral,

Fij =	1	òAi òΩ	cos θi	Vij dω j dAi	(4.9)
	Ai
			π

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

CHAPTER 4. THE FORM FACTOR

4.4. COMPUTING THE FORM FACTOR

where Vij is now 1 if element j is visible from dAi in direction d ω j . The new geometric kernel, Vij cos θi/π, will be denoted Giω. This alternate form of the integral will be used in some algorithms to compute the form factors from dAi to all elements at once.

If the elements are assumed to be fully visible to one another, a third variation of the form factor integral can be derived by applying Stokes’ theorem to convert the double area integral into a double contour integral [143, 222]. The result (not derived here) is

Fij	=		1	ln r dxi dx j + ln r dyi dyj + ln r dzi dz j	(4.10)
			π Ai
	2			òCi òCj

where Ci and Cj are the boundaries of elements i and j. This third form of the integral will also be used by some algorithms when visibility can be determined

4.4 Computing the Form Factor

Both closed form analytic and numerical methods have been applied to solving the form factor integral.4

4An excerpt from Schröder and Hanrahan [206]: The history of computing the amount of light impinging on a diffusely reflecting surface from some light source is very long. A closed form expression for the form factor between a differential surface element and a polygon had already been found by Lambert in 1760 [143]. Lambert proceeded to derive the form factor for a number of special configurations among them the form factor between two rectangles sharing a common edge with an angle of 90 degrees between them. He writes about the latter derivation:

Although this task appears very simple its solution is considerably more knotted than one would expect. For it would be very easy to write down the differential expression of fourth order which one would need to integrate four fold; but the highly laborious computation would fill even the most patient with disgust and drive them away from the task. The only simplification which I was able to achieve was to reduce the expression to a second order differential using [the formula for differential surface element to polygon form factor (equation 4.15)] with which I was able to perform the computation.

Lambert also formulates the reciprocity principle in his theorem 16 and uses form factor algebra to compute unknown factors from known ones. The first use of Stokes’ theorem [224] to solve for the form factor between two arbitrary surfaces can be found in a book by Herman in 1900 [126]. Through two applications of Stokes theorem he reduces the form factor between two arbitrary surfaces to a double contour integral. He uses this result to give the form factor for two parallel quadrilaterals in an exercise. A similar derivation can be found in an article by Fock in 1924 [83]. Fock proceeds by applying the formulation to elliptical disks for which he derives a closed form solution. In 1936 Moon [169] aware of Fock’s works, derives closed form solutions for a number

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

CHAPTER 4. THE FORM FACTOR

4.4. COMPUTING THE FORM FACTOR

Figure 4.3: A taxonomy of form factor algorithms.

A taxonomy of form factor computation methods is shown in Figure 4.3. The discussion of these methods will begin with an examination of closed form (analytic) solutions. Although there is no closed form solution to the general form factor integral, if the two elements are assumed (or can be determined) to be fully visible to each other, there are some useful analytic formulae.

of specialized configurations. In the same year Gershun [93] puts various photometric quantities on a vector calculus footing and gives an especially elegant derivation of the double contour integration using differential forms. Sparrow [221] in 1963 used the double contour form to derive closed form solutions for the case of parallel disks and parallel quadrilaterals. However none of these sources or any since that we are aware of has iven a closed form solution of the form factor between two general polygons.”

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

CHAPTER 4. THE FORM FACTOR

4.5 FORMULAE FOR SIMPLE SHAPES

II.Closed Form Solutions for the Form Factor

4.5Formulae for Simple Shapes

Analytic formulae for specialized geometries can be found in the appendices of most radiative heat transfer texts [131, 155, 216, 222]. The formulae given in Figure 4.4 for opposing and perpendicular rectangles are typical. Although the geometries are simple and visibility is not an issue, the analytic formulae ale far from straightforward.

Analytic formulae are often used in conjunction with form factor algebra, which allows the form factors for the union or difference of simple areas to be computed from the form factors to the individual areas. As shown in Figure 4.5, the form factor from element i to element j plus the form factor from element i to element k must be equal to the form factor from element i to a new element made up of the union of j and k. Similarly, the form factor from the union of j and k to element i is the area average of the two individual elements. Thus, if an element shape can be decomposed into simple shapes for which analytic form factors are known, the form factor algebra can be used to determine the full form factor.

Closed form formulae are also available for a differential area to various finite geometries. An example is shown in Figure 4.6 for a parallel, axially aligned disk. The differential area to finite area form factor arises naturally in point collocation, where the radiosity equation is evaluated at element nodes. It is also used in some numeric algorithms, since the outer integral in equations 4.7 and 4.9 can be performed numerically by evaluating the inner integral at one or more locations on region i and averaging the result.

4.6 Differential Area to Convex Polygon

The analytic formula for the form factor from a differential area to a polygon is particularly useful, since polygonal models are often encountered in image synthesis [19, 130, 143, 175]. The geometry for this formula is given in Figure 4.7. The form factor is computed as a sum around the perimeter of the polygon:

						n
FdAi → Aj			=	1		åβi		cos αi	(4. 15)
				2π
						i =1
or equivalently,
			n
FdAi → Aj =	1		åβi Ni			•	(Ri	× R(i +1)%n )	(4.16)
	2π
			i =1
Radiosity and Realistic Image Synthesis									72
Edited by Michael F. Cohen and John R. Wallace

CHAPTER 4. THE FORM FACTOR

4.6 DIFFERENTIAL AREA TO CONVEX POLYGON

Figure 4.4: Analytic form factors between rectangles.

where n is the number of sides on the polygon, bi is the angle between Ri and R(i+1)%n in radians, ai is the angle between the plane of differential area dAi and the triangle formed by dAi and the ith edge of the polygon, and Ni is the unit normal to dAi.

This formula does not take into account occlusion. However, in conjunction with the appropriate visibility algorithm, the form factor algebra can be used to compute an exact result. In an algorithm by Nishita and Nakamae [175], the form factor to the original polygon is first computed, ignoring occlusion.

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

CHAPTER 4. THE FORM FACTOR

4.7 GENERAL POLYGON TO POLYGON

Fi,( j + k ) = Fi, j + Fi,k	F( j + k ),i =	F j,i Aj	+ Fk,i	Ak	(4.13)
		Aj	+ Ak

Figure 4.5: Form factor algebra.

Fij =	r2	(4.14)
	h2 + r2

Figure 4.6: Analytic form factor from point to disk.

Other polygons in the scene are clipped against this polygon (as viewed from the differential area) and against each other. The form factor to each of the clipped polygons is computed and subtracted from the unoccluded form factor, giving the form factor to the visible portion of the original polygon.

4.7 General Polygon to Polygon

Schröder and Hanrahan give a closed form solution for general polygon-to- polygon form factors, ignoring occlusion [206]. The formulation is non-elemen- tary, as it is based on the dilogarithm [151] arising in the integration of the contour integral form of the form factor, equation 4.10. The specifics of the closed form solution are not repeated here as they involve a long series of complex terms.

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

CHAPTER 4. THE FORM FACTOR

4.7. GENERAL POLYGON TO POLYGON

Figure 4.7: Differential area to polygon form factor.

III. Numerical Solutions for the Form Factor

Closed form analytic formulae do not lend themselves to the direct evaluation of form factors between complex shapes or where occlusion is a factor. For more general situations, numerical approaches are required to approximate the form factor.

Numerical methods of evaluating integrals (known as quadrature methods) generally involve sampling the kernel at various points and approximating the integral as a weighted sum of the kernel evaluated at the sample points. In general, the more sample points selected, the more accurate the approximation. However, the cost of approximating the integral is directly related to the number of kernel evaluations required (each of which generally requires a visibility calculation). Thus, the goal in developing a numerical algorithm (or quadrature rule) to solve the form factor integral is to get the most accuracy with the fewest (and/or cheapest) kernel evaluations.

There are a number of choices to make in designing a quadrature rule for form factor evaluation. First, in the case of constant basis functions, one can

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

CHAPTER 4. THE FORM FACTOR

4.8. NUMERICAL INTEGRATION IN GENERAL

choose to evaluate any of the different formulations of the form factor integral given in section 4.3. The area-area form (equation 4.7) and the area-hemisphere form (equation 4.9) are equivalent, and the contour form (equation 4.10) is also a suitable choice if one can determine a priori that the elements are fully visible to one another.

In addition, there is not just a single form factor to be evaluated, but rather a matrix of form factors. Each entry of a row (or column) of the form factor matrix shares a common element as the domain of integration, thus one may take advantage of this coherence by simultaneously solving for a complete row and/or column of form factors.

Finally, one is free to a great extent to choose where to sample the kernel and how to weight the samples so long as it can be shown that as the number of samples increases, the weighted sum converges to the true solution.

After a brief discussion of the general topic of numerical integration, the chapter will proceed with a description of a variety of numerical algorithms that have been developed for form factor evaluation. Other surveys of form factor algorithms such as [78, 187] provide additional insights. The approaches are broadly classified according to which form of the form factor integral they use, how the sample points are selected and weighted, and in what order the form factors are evaluated (e.g., one at a time or a row at a time).

4.8 Numerical Integration in General

Generically, quadrature rules approximate some integral H with kernel h(x) over the domain X as a weighted sum :

	n
ˆ	åwk h(xk )	(4.17)
H = òX h(x) dx ≈ H =

k =1

One is free to use any information available about h(x) in choosing the quadrature points, xk .

Normally, one would like to make n as small as possible to limit the cost incurred in the evaluation of the h(xk ), without compromising accuracy. The simplest methods, like the trapezoidal rule or Simpson’s rule, sample the domain at evenly spaced intervals, evaluate the kernel at these points, and sum the results weighted by the size of the interval. Clearly, as the number of samples increases, and the size of the intervals decreases the approximation approaches the true integral.

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

CHAPTER 4. THE FORM FACTOR

4.8. NUMERICAL INTEGRATION IN GENERAL

4.8.1 Gaussian Quadrature

Simple methods like the trapezoidal rule ignore any knowledge of the integrand. A more efficient selection of quadrature points and weights can be made given the available knowledge about the nature of the integrand h(x). For example, if one knows h(x) is constant across the limits of the integral, then one quadrature point anywhere is sufficient and the weight is simply the difference of the upper and lower limits of integration (in our case, the area of the elements). In general, the smoother the integrand, the fewer quadrature points required to approximate the integral to within some given error.

In one dimension, Gaussian quadrature methods [185] can be used to evaluate exactly integrals of polynomials up to order 2n + 1 with n carefully selected points and proper weights. The theory behind this observation is quite elegant [226]. The specific quadrature points and associated weights are tabulated and can be found in most numerical methods books or can be computed through recurrence relationships. Extensions to multiple dimensions as in the form factor problem are possible but can be difficult due to the exponential growth in the number of quadrature points required with respect to dimension.

4.8.2 Quadrature Points and the Form Factor Integral

In the double area integral form of the form factor (equation 4.7), the quadrature point is now defined in the four-dimensional space, R2 3 R2, of the combined elements.5 In other words, a quadrature point represents the selection of a pair of 2D points, x and x′, in elements i and j at which to evaluate the integrand. Similarly, in the area-hemisphere form (equation 4.9), a quadrature point is in the space R2 3 S2, that is, a point (x, ω ) is in the combined space of an element area and the hemisphere of directions above that point.

4.8.3 Monte Carlo Methods

Monte Carlo methods are a family of quadrature methods of very general applicability, often used for difficult equations where other methods are impractical or unavailable. Monte Carlo techniques use random numbers to select sample locations in the domain of integration at which to evaluate the integrand.6 The integral is then taken to be a weighted average of the kernel evaluation at sample points. The weights associated with each sample evaluation depend on how

5R is the space of the real number line and S is the space of directions on a circle. Thus R2 corresponds to a plane and S2 corresponds to the surface of a sphere.

6Quasi-random distributions such as the Poisson disk may also be used. The samples in a quasi-random distribution are not completely independent but have statistical properties that allow them to be used in place of random samples. See for example discussions by Shirley in [212].

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

CHAPTER 4. THE FORM FACTOR

4.8. NUMERICAL INTEGRATION IN GENERAL

the samples are selected. If, for example, samples are evenly distributed in the domain, then the weights are simply for n samples with A the size of the domain of integration.

One would like, however, to make as few evaluations of the integrand as possible to reach an answer with some level of confidence of being within some error bound. This process can be enhanced by taking advantage of whatever knowledge one has of the integrand. For example, if the integral H to be evaluated has a kernel h(x) that can be factored, h(x) = f(x) • g(x), where g(x) is a simple known positive function (g(x) > 0 for x X), then the integral can be written

H = òX f(x) g(x) dx

(4. 18)

In this case, one would like to find a distribution of samples xk in X that mimics the known function, g(x). This method of selecting samples is called importance sampling, since more samples are taken where g(x) is large (important) and fewer are taken where g(x) is small. More formally,

H	=	òX f(x) g(x) dx =		òX f(x) Gp(x) dx			(4.19)
where		òX g(x)dx
G	=		and	p(x) =	g( x )		(4.20)
					G
p(x) is essentially a	normalized g(x)		(i.e., ò–∞∞ p( x ) = 1) and is				thus a prob-

ability density function. The cumulative density function P(x) can be defined
as	x
	P(x) = ò— ∞ p(x) dx	(4.21)

Loosely speaking, if p(x) represents the odds of picking x, then P(x) represents the odds of picking some number less than or equal to x. If the inverse of the function P(x) is P—1 (x), then the Monte Carlo evaluation for the approximation

ˆ	≈ H is simply
H
	ˆ
	H = 0;
	for ( k = 1 to n) {				for n samples
	choose ξ ;				randomly in the interval from 0 to 1
	x = P—1	(ξ) ;			x will be chosen with probability p(x)
	ˆ	ˆ			sum the sample values of f(x)
	H =	H + f(x) ;
	}
	ˆ	ˆ	G
	H =	H *		;	normalize by G and divide by n samples
			n
Radiosity and Realistic Image Synthesis						78
Edited by Michael F. Cohen and John R. Wallace