Advanced Computer Graphics March 31, 2003 G. Farin

Radiosity

1 Form factors

Assume a scene is composed of n patches P_{i}. A subset of these patches encloses the scene. Each patch receives light from other patches and reflects light back into the scene. All reflections are assumed to be diffuse only. The interaction between patch i and patch j is given by a form factor f_{ij}. It determines how much light P_{i} receives from patch P_{j}:

f_{ij} = \frac{light from P_{j} to P_{i}}{all light from P_{j}} .

We immediately have:

∑_{i=1}^{n} f_{ij} = 1; j = 1, …, n.

The f_{ij} are given by

f_{ij} = \frac{1}{∥P_{i}∥} ∫_{P_{i}} ∫_{P_{j}} \frac{cos Φ_{i} cos Φ_{j}dP_{j}}{πr^{2}}  (1)

where the meaning of the involved terms is from Figure 1.

In order to motivate (??), assume that P_{i} has very small area. Then, again using Figure 1, we would have

f_{ij} = \frac{cos Φ_{i} cos Φ_{j}∥P_{j}∥}{r^{2}}.

From this we arrive at (1) by integrating over the areas of P_{i} and P_{j}. Note that we have

f_{ij}∥p_{i}∥ = f_{ji}∥P_{j}∥.  (2)

An easier way to determine the form factors is by the use of a hemicube, a simplified model of a hemisphere. Subdivide each of the five faces of the hemicube into squares (pixels), resulting in a total of M pixels. Center the upper half of a cube at a subpatch

Figure 1: The form factor geometry.

P_{i,ρ} of P_{i} such the top is parallel to P_{i}. We assume that are R such subpatches.

Now project P_{j} onto the hemicube. Count all pixels visible from P_{i,ρ}. Let the number of these pixels be n_{i,j,ρ}. Keep in mind that we only count pixels which are not occluded by other objects. We repeat this for all subpatches P_{i,ρ}, which means a recentering of the hemicube. Then

f_{ij} = \frac{1}{∥P_{i}∥} ∑_{ρ=1}^{R} \frac{∥P_{i,ρ}∥n_{ijρ}}{M}

2 Setting up the linear system

Let b_{i} the brightness of patch P_{i}. Let e_{i} be the amount of light emitted by P_{i}, let ρ_{i} be the reflectivity of P_{i}, and let a_{i} be its area. Now the total intensity of light leaving P_{i} is given by

a_{i}b_{i} = a_{i}e_{i} + ρ_{i} ∑_{j=0}^{n} f_{ij}a_{j}b_{j}; i = 1, …, n.  (3)

The first term gives how much light is emitted by P_{i}. The second term collects all contributions from the other patches; the factor ρ_{i} shows how much of that is being reflected by P_{i}.

Invoking (2), we may rewrite (3) as

a_{i}b_{i} = a_{i}e_{i} + ρ_{i} ∑_{j=0}^{n} f_{ji}a_{i}b_{j}; i = 1, …, n

and simplify to

b_{i} = e_{i} + ρ_{i} ∑_{j=0}^{n} f_{ji}b_{j}; i = 1, …, n.  (4)

In matrix form:

b = e + RF b  (5)

or

e = [I − RF]b  (6)

where R is a diagonal matrix holding the ρ_{i} and F is the matrix of form factors. The matrix I − RF has 1’s on the diagonal since we can assume f_{ii} = 0.

3 Solving the system

This linear system for b may be huge, but it will always be sparse. For such systems, an iterative solution is most effective. In an iterative scenario, one produces a first guess b^{0} for b, e.g., by setting all b_{i} = 0.5. Then a next guess b^{1} is found from

b^{1} = e + RF b^{0}

and this process is continued until convergence happens, which is guaranteed for this scenario. In practice a vector b^{k} is called the solution if ∥b^{k} − b^{k−1}∥ < ϵ for some tolerance epsilon. This iterative process is known as Gauss − Jacobi iteration.

In that method, all elements of b^{k} is updated from b^{k−1} in one step. A different strategy is usually more effective. Once the first element of b^{k} has been computed, the second element of b^{k} may be computed by using this new b_{1}^{k} instead of the old b_{1}^{k−1}. Similarly, the b_{3}^{k} may be computed using b_{1}^{k} and b_{2}^{k}. This process of immediate updating is known as Gauss-Seidel iteration.

The matrix I − RF is strictly diagonally dominant, with all eigenvalues less than unity in absolute value. For such matrices, a power expansion exists for the inverse:

[I − RF]^{−1} = ∑_{i=0}^{∞}[RF]^{i}.

Thus

b = ∑_{i=0}^{∞}[RF]^{i}e = e + [RF]e + [RF]^{2}e + …

Each term corresponds to one level of tracking reflections.