Chapter 11

Global Illumination

Global illumination in a 3D scene results from the propagation of light from a light source in the scene through multiple inter-reflections. This chapter presents a few algorithms to compute color due to the global illumination in the scene. Most of these algorithms simulate this propagation in one form or another. Simulation should normally proceed with the distribution of light from the light source to all other surface patches in the environment. The simulation continues by distributing reflected light from those surface patches into the environment. This light distribution process is continued till the equilibrium is reached. Some of the simulation algorithms are designed to compute equilibrium light at only those points that are visible to the camera. Most such algorithms are based on ray tracing. We provided the idea of ray tracing in the first chapter. Here, we come back to it in order to provide more details and insights. Other algorithms compute equilibrium light for the whole scene, and thus compute equilibrium light at every surface of the scene, independent of whether the surface is visible or not. Two algorithms of this type that we describe here are the Photon Mapping and the Radiosity algorithms.

We would like to underline that this chapter is presented in quite a different way from the previous ones, since no practical implementation of the described algorithms is provided as in the usual Upgrade Your Client sections. This is because our main aim is to provide here only the basic ideas behind global illumination algorithms and not to explain how to implement them efficiently. In fact, this is a complex and highly specialized topic. For the readers interested in going deeper into this topic, we refer them to other textbooks such as the book on photon mapping by Henrik Wann Jensen [15], which includes an implementation of photon mapping in C language, and the global illumination books by Dutre et al. [9] and by Pharr and Humphreys [34], the latter of which contains many details about the art of ray tracing.

11.1 Ray Tracing

A ray is defined by an origin o, the point at which it originates and a direction d, the direction along which the ray extends. The actual extent of the ray can in principle be infinity. Any point p on the ray is defined by a scalar parameter t such as $\mathbf{\text{p}}=\mathbf{\text{o}}+\mathit{\text{t}}\mathit{\text{d}}$. Every $\mathit{\text{t}}\ge 0$ generates all valid points along the ray. If the direction is normalized then t represents the Euclidean distance of the point from the origin of the ray.

Ray tracing plays a very important role in global illumination computation:

• (a) For simulating the propagation of light rays originating at the light source through a scene (light ray tracing or photon tracing).
• (b) For computing the amount of light reaching the camera through a pixel by tracing a ray from the camera, following the ray through the scene, and collecting all the lighting information (classical ray tracing, Monte Carlo ray tracing, etc.).

Tracing a ray means extending the ray from its origin, and collecting some information along the ray. The exact information collected depends on the application and the type of scene. We restrict our discussion to scenes containing solid opaque objects. Independent of the application, the major computational effort in tracing the ray is the computation of the ray-scene intersection. As the scene is assumed to be composed of one or more solid objects, intersecting the ray with the scene involves computation of ray–object intersection. As parameter t represents a point along the ray, computing ray–object intersection can be resolved by computing the ray parameter t. Plenty of research has been devoted to finding efficient ray–object intersection algorithms for a large number of object types. We will restrict our discussion on ray–object intersection to only two classes of objects: algebraic and parametric surfaces.

11.1.1 Ray–Algebraic Surface Intersection

Each point p of an algebraic surface satisfies an algebraic equation of the type f(x, y, z) = 0 where f is a polynomial expression involving the coordinates of the point. For a ray to intersect such a surface there must be at least one point pi common between the surface and the ray. Or, in other words:

$\begin{array}{ccc}\mathit{\text{f}}\left({\mathit{\text{p}}}_{\mathit{\text{i}},\mathit{\text{x}}},{\mathit{\text{p}}}_{\mathit{\text{i}},\mathit{\text{y}}},{\mathit{\text{p}}}_{\mathit{\text{i}},\mathit{\text{z}}}\right)=0\text{and}{\mathbf{\text{p}}}_{\mathit{\text{i}}}=\mathbf{\text{o}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}\mathit{\text{d}}& & \left(11.1\right)\end{array}$

By substitution we get an algebraic equation of ti as

$\begin{array}{ccc}\mathit{\text{f}}\left({\mathit{\text{o}}}_{\mathit{\text{x}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{x}}},{\mathit{\text{o}}}_{\mathit{\text{y}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{y}}},{\mathit{\text{o}}}_{\mathit{\text{z}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{z}}}\right)=0& & \left(11.2\right)\end{array}$

Any polynomial root-finding method may be used to compute the value of ti. We detail here a couple of examples: planes and spheres.

11.1.1.1 Ray-Plane Intersection

The plane has the simplest of the algebraic equations:

$\begin{array}{ccc}\mathit{\text{a}}\mathit{\text{x}}+\mathit{\text{b}}\mathit{\text{y}}+\mathit{\text{c}}\mathit{\text{z}}+\mathit{\text{d}}=0& & \left(11.3\right)\end{array}$

Substituting the ray equation into the plane equation we get:

$\begin{array}{ccc}\mathit{\text{a}}\left({\mathit{\text{o}}}_{\mathit{\text{x}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{x}}}\right)+\mathit{\text{b}}\left({\mathit{\text{o}}}_{\mathit{\text{y}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{y}}}\right)+\mathit{\text{c}}\left({\mathit{\text{o}}}_{\mathit{\text{z}}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}{\mathit{\text{d}}}_{\mathit{\text{z}}}\right)+\mathit{\text{d}}=0& & \left(11.4\right)\end{array}$

or

$\begin{array}{ccc}{\mathit{\text{t}}}_{\mathit{\text{i}}}=-\frac{\left(\mathit{\text{a}}{\mathit{\text{o}}}_{\mathit{\text{x}}}+\mathit{\text{b}}{\mathit{\text{o}}}_{\mathit{\text{y}}}+\mathit{\text{c}}{\mathit{\text{o}}}_{\mathit{\text{z}}}+\mathit{\text{d}}\right)}{\left(\mathit{\text{a}}{\mathit{\text{d}}}_{\mathit{\text{x}}}+\mathit{\text{b}}{\mathit{\text{d}}}_{\mathit{\text{y}}}+\mathit{\text{c}}{\mathit{\text{d}}}_{\mathit{\text{z}}}\right)}& & \left(11.5\right)\end{array}$

Thus the ray–plane intersection computation is the simplest of all ray–object intersections.

11.1.1.2 Ray–Sphere Intersection

The sphere has the simplest of the algebraic equations of degree two. The algebraic equation for a sphere whose radius is r and whose center is located at c is

$\begin{array}{ccc}{\left(\mathit{\text{x}}-{\mathit{\text{c}}}_{\mathit{\text{x}}}\right)}^2+{\left(\mathit{\text{y}}-{\mathit{\text{c}}}_{\mathit{\text{y}}}\right)}^2+{\left(\mathit{\text{z}}-{\mathit{\text{c}}}_{\mathit{\text{z}}}\right)}^2-{\mathit{\text{r}}}^2=0& & \left(11.6\right)\end{array}$

which can be easily written as

$\begin{array}{ccc}\left(\mathbf{\text{p}}-\mathbf{\text{c}}\right)\cdot \left(\mathbf{\text{p}}-\mathbf{\text{c}}\right)-{\mathit{\text{r}}}^2=0& & \left(11.7\right)\end{array}$

For a ray to intersect a sphere, the point of intersection pi must satisfy the following equation:

$\begin{array}{ccc}\left({\mathbf{\text{p}}}_{\mathit{\text{i}}}-\mathbf{\text{c}}\right)\cdot \left({\mathbf{\text{p}}}_{\mathit{\text{i}}}-\mathbf{\text{c}}\right)-{\mathit{\text{r}}}^2=0\text{and}{\mathbf{\text{p}}}_{\mathit{\text{i}}}=\mathbf{\text{o}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}\mathit{\text{d}}& & \left(11.8\right)\end{array}$

Substituting of the ray equation into the sphere equation, we get:

$\begin{array}{ccc}\left(\mathbf{\text{o}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}\mathit{\text{d}}-\mathbf{\text{c}}\right)\cdot \left(\mathbf{\text{o}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}\mathit{\text{d}}-\mathit{\text{c}}\right)-{\mathit{\text{r}}}^2=0& & \left(11.9\right)\end{array}$

which on expansion gives us

$\begin{array}{ccc}{\mathit{\text{t}}}_{\mathit{\text{i}}}^2\left(\mathbf{\text{d}}\cdot \mathbf{\text{d}}\right)+2{\mathit{\text{t}}}_{\mathit{\text{i}}}\left(\mathit{\text{d}}\cdot \left(\mathbf{\text{o}}-\mathbf{\text{c}}\right)\right)+\left(\left(\mathbf{\text{o}}-\mathbf{\text{c}}\right)\cdot \left(\mathbf{\text{o}}-\mathbf{\text{c}}\right)\right)-{\mathit{\text{r}}}^2=0& & \left(11.10\right)\end{array}$

which is a quadratic equation in t. Quadratic equations have two roots, and hence ray–sphere intersection can result in two intersection points. The smallest of the positive real roots is chosen as the parameter of intersection.

11.1.2 Ray–Parametric Surface Intersection

We remind the reader that parametric objects’ surfaces are often expressed as a parametric equations of a form p = s(u, v), that is, every coordinate of the point on a parametric surface is a function of two parameters u and v. Hence, for a ray to intersect a parametric surface, the equality

$\begin{array}{ccc}\mathbf{\text{o}}+{\mathit{\text{t}}}_{\mathit{\text{i}}}\mathit{\text{d}}={\mathbf{\text{p}}}_{\mathit{\text{i}}}=\mathit{\text{g}}\left({\mathit{\text{u}}}_{\mathit{\text{i}}},{\mathit{\text{v}}}_{\mathit{\text{i}}}\right)& & \left(11.11\right)\end{array}$

must hold. So we get a system of three equations with three unknown parameters,

$\begin{array}{ccc}\mathit{\text{g}}\left({\mathit{\text{u}}}_{\mathit{\text{i}}},{\mathit{\text{v}}}_{\mathit{\text{i}}},{\mathit{\text{t}}}_{\mathit{\text{i}}}\right)=0& & \left(11.12\right)\end{array}$

Depending on the surface the equations can be non-linear. In such a case multivariate Newton iteration may be used to solve the system of equations.

A triangle T = {p1, p2, p3} may be considered the simplest of the parametric surfaces whose points satisfy the equation p = a + bu + cv, where u, v are scalar parameters in the range [0, 1], u + v is also in the range [0, 1] and a, b and c are related to the three vertices of the triangle as: a = p1, b = p2−p1 and c= p3−p1. For ray–triangle intersection, we get the following system of linear equations:

$\begin{array}{ccc}\left[\begin{array}{ccc}-\mathit{\text{d}}& \mathit{\text{b}}& \mathit{\text{c}}\end{array}\right]\times \left[\begin{array}{c}{\mathit{\text{t}}}_{\mathit{\text{i}}}\\ {\mathit{\text{u}}}_{\mathit{\text{i}}}\\ {\mathit{\text{v}}}_{\mathit{\text{i}}}\end{array}\right]=\left[\mathbf{\text{o}}-\mathbf{\text{a}}\right]& & \left(11.13\right)\end{array}$

An often-used method of computing ray–triangle intersection is solving this linear system. If the computed ti value is greater than zero, and ui, vi and their sum ui + vi are in the range [0, 1], then the ray intersects the triangle, and the point of intersection pi is equal to o + tid.

11.1.3 Ray–Scene Intersection

The simplest method of computing ray–scene intersection is intersecting the ray with each and every element of the scene, and choosing the point that corresponds to the smallest positive ray parameter ti. In this approach, the cost of per-ray intersection is linear in terms of the number of objects in the scene, and is acceptable for scenes composed of only a small number of objects. However, most realistic scenes are composed of hundreds to a hundred thousand or more objects. In such cases the linear approach to computing ray-scene intersection is unacceptably slow. So over the years a number of approaches have been proposed to accelerate this operation. Some of the widely used acceleration approaches (discussed in detail in [34]) are: uniform spatial subdivision (USS) and hierarchical structures such as kD tree and bounding volume hierarchy (BVH). All of these approaches rely on pre-computing acceleration structures that reduce the number of ray–object intersection tests, by intersecting the ray with only a few objects that are likely to lie along the ray. Independent of the acceleration approach used, computing an axis aligned bounding box (AABB) for the scene and checking its intersection with the ray is the starting point for every approach. So, we first describe what an AABB is and how the ray–AABB intersection test is calculated. Then, we continue discussing the USS and BVH methods for accelerating ray–scene intersection.

11.1.3.1 Ray–AABB Intersection

As previously stated, the bounding volume of a scene is a simple geometric primitive that entirely covers all the vertices of the scene. An axis-aligned bounding box (AABB) is the simplest bounding volume to compute. As the name suggests, this type of bounding volume is a rectangular box with its bounding planes aligned to the axis. In other words, the six bounding planes of an AABB are parallel to the three main axial planes. An oriented bounding box (OBB) is another type of bounding box that is more compact in enclosing the volume of the object, but is not constrained to have its planes parallel to the axial planes. Figure 11.1 shows an example of an AABB and an OBB. In the following we always refer to AABB because this box type is the one used by the algorithm we will describe.

The AABB for a scene is computed by finding the minimum and maximum of the coordinates of the objects in the scene. The minimum and maximum coordinates define the two extreme corner points referred to as cmin = AABB.min and cmax = AABB.max in the following.

Every pair of the faces of an AABB is parallel to a Cartesian plane. Let us consider the two planes parallel to the XY-plane. The Z-coordinates of every point on these planes are equal to cmin,z and cmax,z. So the algebraic equation of the two planes of the AABB parallel to the XY-plane are

$\begin{array}{ccc}\mathit{\text{z}}-{\mathit{\text{c}}}_{\mathrm{\text{min}},\mathit{\text{z}}}=0,\mathit{\text{z}}-{\mathit{\text{c}}}_{\mathrm{\text{max}},\mathit{\text{z}}}=0& & \left(11.14\right)\end{array}$

The parameters of intersection of a ray with these planes are simply:

$\begin{array}{ccccc}\frac{\left({\mathit{\text{c}}}_{\mathrm{\text{min}},\mathit{\text{z}}}-{\mathit{\text{o}}}_{\mathit{\text{z}}}\right)}{{\mathit{\text{d}}}_{\mathit{\text{z}}}}& \text{and}& \frac{\left({\mathit{\text{c}}}_{\mathrm{\text{max}},\mathit{\text{z}}}-{\mathit{\text{o}}}_{\mathit{\text{z}}}\right)}{{\mathit{\text{d}}}_{\mathit{\text{z}}}}& & \left(11.15\right)\end{array}$

The smallest of these two parameters, tmin,z, corresponds to the nearest point of intersection of the ray with the AABB parallel to the XY-plane and tmax,z corresponds to the farthest point of intersection. Similarly we can compute the ray parameters tmin,x, tmax,x and tmin,y, tmax,y corresponding to the nearest and farthest point of ray intersection with the planes parallel to the YZ-plane and ZX-plane, respectively. The intersection of a ray with a rectangular box can result in a pair of intersection points. These two points can be computed from the three pairs of ray–box–plane intersection points. If the intersection exists, then the nearest of the farthest intersection points will be the farthest of the ray–AABB intersection points, and the farthest of the nearest intersection points will be the nearest ray–AABB intersection point. Thus the ray parameter for the nearest point of ray–AABB intersection tmin is max(tminx,x, tmin,y, tmin,z), and tmax, that of the farthest point of intersection is min(tmax,x, tmax,y, tmax,z). If tmin < tmax then the ray does indeed intersect the AABB and the nearest point of ray–AABB intersection is given by the ray parameter tmin. The pseudo-code for ray–AABB intersection computation is given in Listing 11.1.

function ray‒AABB()

  INPUT ray, AABB

 // ray parameter for the point of intersection

 OUTPUT t_min,t_max

{

 tmin,x = (AABB.min.x‒ray.o.x)/ray.d.x

 tmax,x = (AABB.max.x‒ray.o.x)/ray.d.x

 if (tmin,x > tmax,x)

   swap(tmin,x,tmax,x)

 endif

 tmin,y = (AABB.min.y‒ray.o.y)/ray.d.y

 tmax,y = (AABB.max.y‒ray.o.y)/ray.d.y

 if (tmin,y > tmax,y)

   swap(tmin,y, tmax,y)

 endif

 tmin,z = (AABB.min.z‒ray.o.z)/ray.d.z

 tmax,z = (AABB.max.z‒ray.o.z)/ray.d.z

 if (tmin,z > tmax,z)

   swap(tmin,z, tmax,z)

 endif

 tmin = max(tmin,x, tmin,y, tmin,z)

 tmax = min(tmax,x, tmax,y, tmax,z)

}

LISTING 11.1: Ray-AABB intersection finding algorithm.

11.1.3.2 USS-Based Acceleration Scheme

As suggested by the name uniform spatial subdivision, in this scheme the space occupied by the scene is subdivided uniformly, and the objects in the scene are assigned to the created subspaces based on their space occupancy (see Figure 11.2). The actual space subdivided is the AABB of the scene. So the subdivision creates a uniform 3D grid structure. A list of objects is associated with every voxel of the 3D grid. If an object is fully or partially inside a voxel then the object is assigned to the list of objects of that cell. The creation of these object lists for the USS grid voxels is a first and important step in USS-based acceleration schemes. The pseudocode of a simple and often used algorithm for this objects list computation for a triangle-only scene is given in Listing 11.2.

function USSpreProcess()

  INPUT scene: AABB, triangleList

  OUTPUT USS:{AABB, N, objectList[N][N][N]}

  // assumes N3as the grid resolution

{

  for every T in triangleList

  pmin = min coordinates of the three vertices of T

  pmax = max coordinates of the three vertices of T

  index_min=ivec3((pmin − AABB.min)/(AABB.max − AABB.min))

  index_max=ivec3((pmax − AABB.min)/(AABB.max − AABB.min))

  for i = index_min.x to index_max.x

  for j = index_min.y to index_max.y

    for k = index_min.z to index_max.z

    append T to USS.objectList[i][j][k]

    endfor

  endfor

 endfor

 endfor

}

LISTING 11.2: USS preprocessing algorithm.

For USS-based ray tracing, the ray is first intersected with the AABB associated with the USS. If there is a valid intersection then the ray is marched through the 3D grid of the USS, one voxel at a time. For each voxel it traverses through, the ray is intersected with the list of triangles associated with the voxel objectList. If a valid intersection point is found then that point rep resents the nearest point of intersection of the ray with the scene, and the ray marching is terminated. Otherwise, the ray is marched through the next voxel, and the process is repeated until the intersection is found or the ray exits AABB. Two crucial steps of the ray marching operation are: find the entry voxel index, and then march from the current voxel to the next voxel along the ray. These operations must be computed accurately and efficiently. Simple use of ray-AABB intersection tests to find the voxel of entry and then finding the next voxel along the path will be highly inefficient. However, an adaptation of that approach has been shown to be very efficient. We describe it below.

11.1.3.3 USS Grid Traversal

Like in AABB, USS is composed of three sets of planes, each of which is parallel to a different Cartesian planes. As we have seen earlier, computing ray intersection with a plane parallel to a Cartesian plane is simple. The ray parameter of intersection with the sets parallel to the YZ, ZX and XY planes are, respectively:

$\begin{array}{lll}{\mathit{\text{t}}}_{\mathit{\text{x}},\mathit{\text{i}}}\hfill & =\hfill & \left(\text{USS}\text{.AABB}\text{.min}\text{.x}+\mathit{\text{i}}\mathrm{\text{Δ}}\mathit{\text{x}}-\text{ray}\text{.o}\text{.x}\right)\text{/ray}\text{.d}\text{.x}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}\mathit{\text{x}}=\left(\text{USS}\text{.AABB}.\mathrm{\text{max}}.\text{x}-\text{USS}\text{.AABB}\text{.}\mathrm{\text{min}}.\text{x}\right)/\mathit{\text{N}},\hfill \\ {\mathit{\text{t}}}_{\mathit{\text{y}},\mathit{\text{j}}}\hfill & =\hfill & \left(\text{USS}\text{.AABB}\text{.min}\text{.y}+\mathit{\text{j}}\mathrm{\text{Δ}}\mathit{\text{y}}-\text{ray}\text{.o}\text{.y}\right)\text{/ray}\text{.d}\text{.y}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}\mathit{\text{y}}=\left(\text{USS}\text{.AABB}.\mathrm{\text{max}}.\mathit{\text{y}}-\text{USS}\text{.AABB}\text{.}\mathrm{\text{min}}.\mathit{\text{y}}\right)/\mathit{\text{N}},\hfill \\ {\mathit{\text{t}}}_{\mathit{\text{z}},\mathit{\text{k}}}\hfill & =\hfill & \left(\text{USS}\text{.AABB}\text{.min}\text{.z}+\mathit{\text{k}}\mathrm{\text{Δ}}\mathit{\text{z}}-\text{ray}\text{.o}\text{.z}\right)\text{/ray}\text{.d}\text{.z}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}\mathit{\text{z}}=\left(\text{USS}\text{.AABB}.\mathrm{\text{max}}.\mathit{\text{z}}-\text{USS}\text{.AABB}\text{.}\mathrm{\text{min}}.\mathit{\text{z}}\right)/\mathit{\text{N}}.\hfill \end{array}$

The indices i, j and k are the indices to the USS 3D grid. The starting values of the indices and sign of the index increment depend on the sign of the coordinates of the ray direction. For example: for the set of planes parallel to the YZ plane, the starting value, the increment for index i, the limit of the index, and the coordinate value at the point of intersection are as follows:

$\begin{array}{ccc}\left(\mathit{\text{i}},\mathrm{\text{Δ}}\mathit{\text{i}},{\mathit{\text{i}}}_{\mathit{\text{l}}\mathit{\text{i}}\mathit{\text{m}}\mathit{\text{i}}\mathit{\text{t}}},\mathit{\text{x}}\right)=\{\begin{array}{cc}\left(\mathit{\text{N}},-1,-1,\text{USS}\text{.AABB}\text{.max}\text{.x}\right)& \text{if ray}\text{.d}\text{.x}<0\\ \left(0,1,\mathit{\text{N}},\text{USS}\text{.AABB}\text{.min}\text{.x}\right)& \text{otherwise}\text{.}\end{array}& & \left(11.16\right)\end{array}$

The formula of ray intersection parameters is not only simple, but also lets us write incremental expressions for the parameters of the next intersection along the ray. The incremental expressions along the axes are respectively:

$\begin{array}{lll}{\mathit{\text{t}}}_{\mathit{\text{x}},\mathit{\text{i}}+\mathrm{\text{Δ}}\mathit{\text{i}}}\hfill & =\hfill & {\mathit{\text{t}}}_{\mathit{\text{x}},\mathit{\text{i}}}+\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{x}}}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{x}}}=\left(\mathrm{\text{Δ}}\mathit{\text{i}}\mathrm{\text{Δ}}\mathit{\text{x}}\right)/\text{ray}\text{.d}\text{.x}\hfill \\ {\mathit{\text{t}}}_{\mathit{\text{y}},\mathit{\text{j}}+\mathrm{\text{Δ}}\mathit{\text{j}}}\hfill & =\hfill & {\mathit{\text{t}}}_{\mathit{\text{y}},\mathit{\text{j}}}+\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{y}}}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{y}}}=\left(\mathrm{\text{Δ}}\mathit{\text{j}}\mathrm{\text{Δ}}\mathit{\text{y}}\right)/\text{ray}\text{.d}\text{.y}\hfill \\ {\mathit{\text{t}}}_{\mathit{\text{z}},\mathit{\text{k}}+\mathrm{\text{Δ}}\mathit{\text{k}}}\hfill & =\hfill & {\mathit{\text{t}}}_{\mathit{\text{z}},\mathit{\text{k}}}+\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{z}}}\hfill \\ \hfill & \hfill & \text{where }\mathrm{\text{Δ}}{\mathit{\text{t}}}_{\mathit{\text{z}}}=\left(\mathrm{\text{Δ}}\mathit{\text{k}}\mathrm{\text{Δ}}\mathit{\text{z}}\right)/\text{ray}\text{.d}\text{.z}\hfill \end{array}$

Figure 11.3 shows the intersection points produced by this process in the 2D case. For the ray whose origin is inside the volume of the USS, the starting index values must be adjusted to guarantee the ray parameter values are greater than zero. There are multiple ways of making this adjustment. An iterative method for this adjustment is given in Listing 11.3.

tx,i = (USS.AABB.min.x + i Δx - ray.o.x)/ray.d.x

ty,j = (USS.AABB.min.y + j Δy - ray.o.y)/ray.d.y

tz,k = (USS.AABB.min.z + k Δz - ray.o.z)/ray.d.z

while (tx,i < 0 && ty,j < 0 && tz,k < 0)

  if (tx,i < ty,j && tx,i) < tz,j)

  tx,i += Δtx

  i+=Δi

  else if (ty,j < tz,k)

  ty,j += Δty

  j += Δj

 else

  tz,k+= Δtz

  k += Δk

 endif

endwhile

LISTING 11.3: Code to take into account the rays originating inside the USS bounds.

Given the start t and index values we can now write a USS-based ray–scene intersection algorithm as shown in Listing 11.4:

function USS_ray_traverse()

  OUTPUT t

{

  while(i ≠ ilimit && j ≠ jlimit j ≠ jlimit && k ≠ klimit

  [t, intersectFlag] = ray_object_list_intersect(ray,objectList[i][j][k]);

  if (intersectFlag == TRUE && pointInsideVoxel(ray.o+t*ray.d,i,j,k))

   // intesection found

  return TRUE;

  endif

  if (tx,i < ty,j && tx,i < tz,j)

   tx,i += Δtx

   i +=Δi

  else if (t(y,j) < tz,k)

   ty,j += Δty

   j += Δj

  else

   tz,k += Δtz

   k +=Δk

  endif

 endwhile

 returnFALSE

}

LISTING 11.4: An incremental algorithm for ray–USS traversal.

The algorithm uses the ray_object_list_intersect function where the ray is intersected with each of the objects in the list and returns the nearest point of intersection if there is one or more intersections. In USS an object is likely to be part of multiple voxels. So the intersection of the ray with the object list associated with the voxel may generate an intersection point outside the voxel, and thus could generate an erroneous nearest point detection. This problem is avoided by checking if the point is inside the current voxel or not. Like AABB, the USS voxel boundaries are parallel to the axis plane. So the point-inside-outside test is easily done by checking the coordinates of the point against the coordinates of voxel bounds.

The resolution of the USS grid could be specified by the user directly or indirectly as average triangle density per voxel. The finer the resolution, the faster the ray–intersection, but this requires longer pre-computation time, a larger amount of memory for the USS storage, and also adds computational overhead for ray-tracing empty parts of the scene. The performance of USS is good in scenes with homogeneous object distribution. However, USS performance can be poor if the extents of the objects in the scene are large and overlap. It also performs poorly for tracing scenes where small but high-triangle-count objects are located, in scenes with large extent, for example a teapot in a football field. To avoid this latter problem, nested USS structures have been proposed. In the next section we describe a hierarchical structure-based ray tracing acceleration technique that better handles this problem.

11.1.3.4 BVH-Based Acceleration Scheme

A bounding volume hierarchy (BVH) scheme uses a recursive partitioning technique to partition objects in the scene into two parts, and continues partitioning till the partitions contain one or a predefined number of objects (see Figure 11.4 for an example). The process creates a binary tree of object lists. As the name suggests, the bounding volume of the objects plays an important role in partitioning the objects. Like in USS, an often-chosen bounding volume is AABB. However, instead of dividing the scene uniformly along each axis as done in USS, the bounding volume is divided hierarchically by a plane parallel to any of the three axial planes. Each division creates two distinct sets of objects. The elements of each set belong to one side of the axial plane of subdivision. The sidedness of the object with respect to the dividing plane is determined by the sidedness of a candidate point of the object. For example, for a scene composed of triangles, the sidedness of a triangle is determined by the sidedness of the centroid point of the triangle. The plane of subdivision is chosen based on a certain heuristic. For example, a simple heuristic could be to choose the plane perpendicular to the longest side of the AABB, and to choose the location of the plane that divides the longest side into two equal halves or to choose the location of the plane such that exactly half of the objects are on one side and the remaining half on the other. After the subdivision, the bounding volume of each subdivision is computed. The two bounding volumes make two nodes of the hierarchy, and the nodes are recursively subdivided till the number of objects in the children nodes reaches a certain predefined minimum. The resulting hierarchy is a binary tree whose every node is a bounding volume. The BVH creation must precede any ray–scene intersection test. A content of the BVH node and sample pseudo-code for BVH creation is given in Listing 11.5.

BVHnode

{

  AABB

  objectList or // No object list for intermediate nodes.

  partitionPlane // No Partition plane for leaf node.

}

function BHVCreate() // recursive function.

  INPUT root: BVHnode

  list : list of objects

  OUTPUT left, right: BVHnode

{

  if(list. size < thresholdSize)

  {

  left= right = null

  objectList= list

  return

  }

  P: Choose subdivision plane // based on certain heuristic

  leftList= empty

  rightList= empty

  for each object in objectlist

  if (object is on left of P) insert(object,leftList)

  else insert(object,rightList)

  endfor

  leftAABB= computeAABB(leftList)

  rightAABB= computeAABB(rightList)

  leftNode= (leftAABB,leftList)

  rightNode= (rightAABB,rightList)

  root.left= leftNode

  root.right= rightNode

  root.partitionPlane= P

  BVHcreate(leftNode,leftList)

  BVHcreate(rightNode,rightList)

}

LISTING 11.5: BVH creation algorithm.

The partitioning used during BVH creation may be used to sort the object list (partition sorting), and replace the object list in the leaf node by an index to the object and the number of objects in the list. In a BVH-based ray acceleration scheme, ray–scene intersection starts with the intersection of the ray with the root AABB of the BVH tree. If there is an intersection, the ray is recursively intersected with the AABBs of the children nodes till the recursion is complete. The order in which children nodes are intersected depends on the ray direction and the plane of partition associated with the node. The recursion is continued even if intersection is found. It is terminated when the recursion stack is empty. In the case that an intersection is found with the leaf node, the t of the ray is set and is used in the ray–AABB intersection test to verify if the point of intersection is already closer. A sample pseudo-code for ray–BVH intersection is given in Listing 11.6.

function rayBHV() // recursive function

  INPUT ray : Ray,

  node : BVHnode.

  OUTPUT ray.t : nearest point of intersection if any,

  hit : true/false.

{

  if node == null return

  if (intersect(node.AABB,ray)

  {

  if (node.isLeaf)

  {

  t = intersect(ray,node.objectList)

  if (t < ray.t) ray.t = t

  }

  else

  {

   if (ray.direction[partitionPlane.axis] < 0)

   rayBVH(node.rightNode)

   else rayBVH(node.leftNode)

  }

 }

}

LISTING 11.6: Ray-BVH intersection-finding algorithm.

In ray–BVH intersection, ray–AABB intersection plays a dominant role, and hence must be efficiently implemented. While the computational effort required for BVH construction is comparable to USS structure construction, the ray–BVH outperforms ray–USS intersection in most scenes. The performance improvement may vary depending on the heuristic used for choosing the node partition plane. Partitioning by surface area heuristics (SAH) is claimed to provide improved performance as compared to dividing the node by spatial size or by number of objects. However, optimal partition finding is still an open research problem.

11.1.4 Ray Tracing for Rendering

As stated in the introduction of this book, ray tracing is used when high quality rendering is in demand. In ray tracing-based rendering, rays are traced from a virtual camera through every pixel of the rendered image to query the scene color. A general ray tracing algorithm looks something like the code shown in Listing 11.7.

function rayTraceRendering()

  INPUT camera, scene

  OUTPUT image

{

  for row = 1 to rows

  for col = 1 to cols

   ray = getRay(row,col,camera)

   image[row][col] = getColor(ray,scene)

  endfor

  endfor

}

LISTING 11.7: Fundamental ray-tracing algorithm.

The first function, getRay, computes the ray through every pixel of the image array. How the ray is computed depends on the type of camera used; for the standard pinhole camera used in rasterization-based rendering as discussed earlier, the process is simple. We first compute the rays in camera space, where the image window is parallel to the XY plane and located near units away from the coordinate origin 0. Then we transform the ray from camera space to world space using inverse camera matrix $\left({\mathit{\text{M}}}_{\text{camera}}^{-1}\right)$. For simplicity we assume that the window is centered around the Z-axis, and the image origin is at the Bottom-Left corner of the image window. Then the ray through the center of the pixel (col,row) is computed as follows:

$\begin{array}{lll}\text{ray}\text{.o}\hfill & =\hfill & {\mathit{\text{M}}}_{\text{camera}}^{-1}{\left(0,0,0\right)}^{\mathit{\text{T}}}\hfill \\ \text{ray}\text{.d}\hfill & =\hfill & {\mathit{\text{M}}}_{\text{camera}}^{-1}{\left(\mathit{\text{w}}\left(\frac{\left(\mathit{\text{c}}\mathit{\text{o}}\mathit{\text{l}}+0.5\right)}{\mathit{\text{c}}\mathit{\text{o}}\mathit{\text{l}}\mathit{\text{s}}}-0.5\right),\mathit{\text{h}}\left(\frac{\left(\mathit{\text{r}}\mathit{\text{o}}\mathit{\text{w}}+0.5\right)}{\mathit{\text{r}}\mathit{\text{o}}\mathit{\text{w}}\mathit{\text{s}}}-0.5\right),-\mathit{\text{n}}\mathit{\text{e}}\mathit{\text{a}}\mathit{\text{r}}\right)}^{\mathit{\text{T}}}\hfill \end{array}$

where Mcamera = PrspM, w and h are the width and height of the image window, and the image resolution is cols × rows. The next function, getColor, gets the color of the light coming along the ray towards the camera. This color is assigned to the pixel. The exact computation carried out in this function distinguishes one ray-traced rendering method from the other. The getColor function may simply return the diffuse color of the object, or may evaluate the direct lighting equation, as explained in Chapter 6, to get the color at the point and return it. The color may be modulated with the texture queried from the texture map associated with the surface. This process is called ray-casting-based rendering. It creates images similar to those created in simple rasterization based rendering. As the rendering using rasterization hardware produces images at a much faster rate, it is uncommon to use ray casting for the same purpose. Most ray-tracing-based rendering methods normally include some form of global illumination computation to get color in their getColor method, and a ray-tracing-based rendering method is distinguished from the others based on what is done in its getColor method. Independent of the methods, the first step in getColor is computing the nearest visible point along the ray originating at the camera, and then computing the color at that visible point. We detail the exact lighting computation technique used in two popular ray-tracing-based rendering methods: classical ray tracing and Monte Carlo ray tracing.

11.1.5 Classical Ray Tracing

Classical ray tracing for rendering was introduced to rendering literature in 1980 by Whitted [41]. This method, like ray-casting-based rendering, computes direct lighting and texturing at the point ray surface intersection, but extends it to include shadows, specular inter-reflections between multiple objects in the scene, and also may include support of transparency. The algorithm given in Listing 11.8 details the method.

function getColor()

  INPUT ray, scene

{

  (t, object, intersectFlag) = raySceneIntersect(ray,scene)

  if (intersectFlag==FALSE)return backgroundColor

  color = black

  for i=1 to #Lights

  shadowRay = computeShadowRay(ray,t,object,scene,lights[i])

  if (inShadow(t,ray,scene) == FALSE)

   color += computeDirectLight(t,ray,scene.lights[i])

  endif

  endfor

  if (isReflective(object)) // Interreflection support

   newRay = reflect(ray,t,object)

   color += object.specularReflectance * getColor(newRay,scene)

  endif

  if (isRefractive(object)) // transparency support

   newRay = refract(ray,t,object)

   color += object.transparency * getColor(newRay,scene)

  endif

  return color

}

LISTING 11.8: Algorithm for pixel color computation in classical ray-tracing.

As we see from the algorithm, shadow, mirror reflection and transparency are handled naturally by creating additional rays and tracing them. Unlike in rasterization-based rendering, no complicated algorithmic effort is required to support these features. The additional rays are termed secondary rays to distinguish them from the rays originating at the camera, which are called primary rays. For all the secondary rays, the origin is the point of intersection of the primary ray with the scene, but the directions differ.

Shadow ray directions for point light sources are computed by taking the vector difference of the shadow ray origin from the position of the light source. This approach is also extended to area light sources by dividing the surface of the area light source into smaller area patches, and computing the shadow ray directions as the vector differences of the shadow ray origin from the center of (or from randomly chosen points on) the patches. A surface point is under shadow if an object appears between the point and the light source, which translates into finding the ti of a shadow ray–object intersection and checking for the inequality 0 < ti < 1. This ray scene object intersection test, and checking the value of t is done inside function inShadow. Note that inShadow function does not require finding the nearest point of intersection. As soon as an intersection satisfying the condition 0 < t < 1 is found, there is no more need to continue with the intersection of the shadow ray with the other objects of the scene. That is why sometimes it is preferable to write an inShadow algorithm separately from the standard nearest point-finding rayScene intersection algorithm. Furthermore, an object shadowing a point is very likely to shadow the points in the neighborhood. So caching the nearest shadow object information and first checking for the intersection of the shadow ray for the neighboring primary rays with the cached object has been shown to accelerate shadow computation.

Reflection and refraction of light at the nearest point along the primary ray is handled, respectively, by setting the direction of the secondary ray to be the mirror reflection and the refraction of the primary ray and making recursive calls to getColor. The recursion stops when the ray does not intersect any object in the scene, or hits an object that is non-reflective and non-refractive. In a highly reflective/refractive closed environment the algorithm may get into an infinite recursive loop. Though such scenes are uncommon, most algorithms introduce a safety feature by keeping the count of recursive calls, and stopping recursion after the recursion count reaches a certain predefined maximum (often set at 5).

Classical ray tracing is capable of creating accurate rendering with global illumination due to inter-reflection and inter-refraction of light in scenes with mirror-like reflectors and fully or partly transparent objects. Classical ray tracing has been extended to support inter-reflection in scenes with diffuse and translucent surfaces.

11.1.6 Path Tracing

Path tracing is a Monte Carlo-based ray tracing method that renders global illumination effects in scenes with arbitrary surface properties. A path here is a collection of connected rays, starting with primary rays. The rays are recursively generated from the primary rays. So, in effect, path tracing may be considered a simple extension of classical ray tracing algorithm. The difference is that the secondary reflection rays are not restricted to mirror-like surfaces, and refracted rays are not restricted only to transparent surfaces. The ray generation from arbitrary surfaces is done by Monte Carlo direction sampling. Unlike in classical ray tracing, instead of generating a primary ray from the camera through the center of the pixel, in path tracing multiple primary rays are normally generating through random positions in the pixel. So a path tracing algorithm for a scene containing opaque surfaces has the form given in Listing 11.9 (see Figure 11.5).

function pathTraceRendering()

INPUT camera, scene

OUTPUT image

{

  for row = 1 to rows

 for col = 1 to cols

  for i = 1 to N

   ray = getPathRay(row,col,camera)

   image[row][col] = getPathColor(ray,scene)/N

  endfor

 endfor

}

function getPathColor()

  INPUT ray, scene

{

  (t, object, intersectFlag) = raySceneIntersect(ray,scene)

  if (intersectFlag==FALSE) returnbackgroundColor

  color = black

  for i=1 to #Lights

 shadowRay = computeShadowRay(ray,t,object,scene,lights[i])

 if (inShadow(t,ray,scene) == FALSE)

  color += computeDirectLight(t,ray,scene.lights[i])

 endif

  endfor

  if (isReflective(object)) // Interreflection support

  (newRay, factor) = sampleHemisphere(ray,t,object)

  color += factor * getPathColor(newRay,scene)

  endif

  return color

}

LISTING 11.9: Path tracing algorithm.

As mentioned in the early part of this section, we can see that path-tracing-based rendering and classical rendering are very similar. The difference is in computing the secondary reflection ray direction, which is done by Monte Carlo sampling of the hemisphere. It generates a random direction on the upper hemisphere of the surface, and computes light along this direction. There exists a number of different ways of sampling this direction: uniform sampling, cosine importance sampling, BRDF importance sampling, etc. Associated with the sampling is a factor that is used to multiply the patch color returned by the getPathColor function. This factor depends on how the sample is generated, and on the BRDF of the surface. Independent of the type of sampling, every Monte Carlo sampling method, also referred to as a stochastic sampling method, depends on a uniform random sampler that samples floating point range [0,1] uniformly. Many programming language implementations provide us with such a sampler function as a part of their respective libraries. Let us call this function rand(). Using this function we will discuss here a couple of hemisphere samplers and the associated multiplication factors.

Uniform Sampling of Hemisphere: A hemisphere is half of a sphere, and can have infinitely many orientations. The orientation of a hemisphere over a surface point depends on the normal at the surface. For direction sampling on an arbitrarily oriented hemisphere, we first make a uniform sampling of a canonical hemisphere, i.e., a hemisphere of unit radius around Z-axis, and then rotate the direction to match the actual orientation of the hemisphere. A sample vector from the center through the canonical hemisphere surface can be described by spherical coordinates (θ,φ), where θ is the angle around the Z-axis, and φ is the angle the projection of the direction vector makes with the X-axis. Once the sample vector is generated it is transformed using the TBN matrix. Thus the uniformly sampled direction on this canonical hemisphere is computed as in Listing 11.10.

function uniformHemisphereSample1()

  INPUT object, p, ray

  OUTPUT direction, factor

{

  // Uniform sampling the canonical hemisphere

  θ = arccos(rand())

  φ = 2π rand()

  sample_d = (sin(θ) cos(φ), sin(θ) sin(φ), cos(θ))T

  // Let T , B, N are the unit tangent, bitangent, and

  // normal vectors at the object point p

  (direction,factor) = ([T B N] * sample_d, 2πcos(θ) (object.brdf(ray.d,

sample_d))) // Rotation

}

LISTING 11.10: Algorithm for uniformly sampling an arbitrarily oriented unit hemisphere.

Another option would be to sample the whole sphere, and discard all those vectors that make an angle greater than 90 degrees with the normal to the surface at the point of interest. In fact the angle checking can be done by checking for the sign of the dot product of the sample vector with the normal, and accepting only those sample vectors that have a positive dot product. As half of the directions are discarded this method may not be considered as efficient, but one advantage of this approach is that no TBN matrix transformation is required. So we also give here (see Listing 11.11) a method based on uniform sampling of the sphere.

function uniformHemisphereSample2()

  INPUT object, p, ray

  OUTPUT direction, factor

{

  Let N be the unit normal vectors at the object point p

  while(TRUE)

  θ = arccos(1‒2*rand())

  phi = 2π rand()

  sample_d = (sin(θ) cos(ϕ), sin(θ) sin(ϕ), cosθ))T

  // Uniform sampling the canonical hemisphere

  if (dot (N,sample_d) > 0)

   (direction,factor) = (sample_d, 2π dot(N,sample_d)object.brdf(ray.d,

sample_d)) // Rotation

  endif

 endwhile

}

LISTING 11.11: Rejection-based algorithm for uniformly sampling an arbitrarily oriented unit hemisphere.

Cosine-Based Importance Sampling of Hemisphere: In this approach, sampled directions have a cosine θ distribution on the hemisphere around the normal, which means the density of the samples are maximum closer to the normal, and the density falls according to cosine function away from the normal. Cosine-sampled directions are preferred over the uniformsampled direction in Monte Carlo lighting computation. It is mostly because the color contributions brought in by the rays traced along directions away from the normal are reduced by a factor of cos θ. So it is considered better to sample according to cosine distribution and not reduce by a factor, than to attenuate the color contribution by a factor of cos θ. As the distribution is very much dependent on the angle the sampled direction makes with the normal vector, we cannot use a full sphere sampling method here. Listing 11.12 shows a popular cosine sampling method.

function cosineHemisphereSample()

  INPUT object, p, ray

  OUTPUT direction, factor

  {

  // Uniform sampling the canonical hemisphere

  $\mathit{\text{θ}}=\mathit{\text{arcsin}}\left(\sqrt{\text{rand}\left(\right)}\right)$

  φ = 2π rand()

  sample_d = (sinθ) cost(ø), sin(θ) sin(ø), cosθ))T

  // Let T , B, N are the unit tangent, bitangent, and

  // normal vectors at the object point p

  (direction,factor) = ([T B N] * sample_d, π object.brdf(ray.d,sample_d)) // Rotation

}

LISTING 11.12: Algorithm for cosine importance sampling a hemisphere.

11.2 Multi-Pass Algorithms

Ray-tracing-based lighting computation methods discussed in the above section are view-dependent computations. That means if the view is changed all the lighting-related computation must be repeated. At times it is preferable to compute the lighting distribution in the scene first (global lighting computation pass) and store it in a data structure. Then, access the data structure during the ray casting or rasterization-based rendering time (rendering pass) to compute the color for the pixel. Such approaches create faster rendering. However, global lighting pass can be very time consuming. This approach is preferable when a faster walk through renderings with global illumination is a requirement, for example, in game rendering. Literature abounds with methods to carry out the global lighting computation pass. We describe here two methods: the photon tracing and radiosity methods.

11.2.1 Photon Tracing

Photon tracing is a two-pass rendering technique. Both the passes mostly use ray tracing technique for ray-scene intersection. The first pass of the technique simulates photon behavior of light in which a large number of photons are emitted from the light source and propagated in the scene till they are either completely absorbed, or lost in the space. Each photon is emitted from certain locations of the light source, with a starting direction. The positions are randomly sampled from the surface of the light source and the direction for the sampled photon is randomly sampled from the hemisphere around the normal to the surface of that position. The photon carries with it a power spectrum that is representative of the light source, and is propagated in the scene. During the propagation the photon travels along a ray determined by the propagation direction. If the photon hits a surface during its propagation then it loses part of its power due to the surface absorption, and is reflected with attenuated power. The photon continues its propagation along the reflected direction and continues with its interactions till it no more hits any surface (lost in space) or its power is attenuated significantly. Each hit event of the photons is stored in a data structure associated with the hit surface or with the whole scene. In the second pass, the actual rendering pass, rays are traced from the eye through the view window, and are intersected with the surface. A circular neighborhood around the point of intersection is searched to find the photon hits and a weighted sum of the power spectrum of the points modulated with the surface BRDF is reflected towards the eye.

Photon tracing technique is a relatively simple, but expensive global illumination computation technique. The quality of the rendering of photon-traced scenes normally depends on the number of photons traced: the larger the number, the better the quality. Lighting computation in a relatively simple scene requires tracing hundreds of millions of photons, and even with such high numbers it is not uncommon to see noise in the parts that are indirectly visible to the light source. A lot of recent research has been devoted to increasing the efficiency and quality of photon-traced rendering.

11.2.2 Radiosity

The radiosity method was introduced in 1984 with the goal to compute inter-reflection of light in diffuse scenes. We have described the term radiosity (conventional symbol B) in Section 6.2, as an alternate term for irradiance and exitance. In the context of global lighting computation the term radiosity refers to the method used to compute equilibrium irradiance in a scene. In the radiosity method the propagation of light is modeled as a linear system of equations relating the radiosity of each surface in the scene to the radiosity of all other surfaces in the scene. The equilibrium radiosity over the surfaces are computed by solving this linear system. In the next section we develop this linear system and describe the various solution methods for solving this system. Finally we describe methods for realistic rendering of scenes with a precomputed radiosity solution.

The radiosity method is based on two important assumptions:

1. It is possible to break down the surfaces of any scene into smaller patches in such a way that radiosity is uniform over the patch. So it is assumed that surfaces of the scene have been broken down into a number of such uniform radiosity patches.
2. The surface patches are flat and are diffusely reflecting in nature.

11.2.3 Concept of Form Factor

In the radiosity method every surface patch plays the role of an emitter and a receiver. Of all the patches only a few belong to actual light sources and are inherently emitters. But all other patches emit the reflected light, and hence are secondary emitters. A part of the light flux arriving on a patch from all other patches on the scene is reflected and the rest of the light is lost in absorption. We derive here a relation between the flux emitted by the i-th emitter patch, Фi and the flux received at the receiver patch, j, from the expression derived for lighting due to area light source in Chapter 6. The expression for the irradiance at a differential area around a point on the receiver patch j is

$\begin{array}{ccc}{\mathit{\text{E}}}_{\mathit{\text{j}}}={\mathit{\text{L}}}_{\mathit{\text{i}}}{\displaystyle \underset{{\mathit{\text{A}}}_{\mathit{\text{i}}}}{\int }\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2& & \left(11.17\right)\end{array}$

where angle (θi) (θj) is the angle formed by the line connecting the two patches dAi and dAj and the normal at patch dAi(dAj), VdAi,dAj and RdAi,dAj are respectively the visibility and distance between the two differential patches. The visibility term is binary in nature, and takes value one or zero depending on whether the differential areas dAi, dAj are visible to each other or not. We assumed that the radiance of the patches is constant over the surface patch, and surfaces are diffuse. Under these conditions the flux and radiance are related by expression L = Φ/(πΑ). So we will use this relation to replace Li, the radiance of patch i, by Фi Next we compute an expression for the flux received by the total surface j due to light emitted from surface i. Here we use the relation between irradiance and flux: dФ = EdA and integrate the differential flux over the whole area to get the following equation.

$\begin{array}{ccc}{\mathrm{\text{Φ}}}_{\mathit{\text{i}}\to \mathit{\text{j}}}={\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{j}}}}{\mathit{\text{E}}}_{\mathit{\text{j}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}={\mathrm{\text{Φ}}}_{\mathit{\text{i}}}/\left(\mathit{\text{π}}{\mathit{\text{A}}}_{\mathit{\text{i}}}\right){\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{j}}}}{\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{i}}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2}}& & \left(11.18\right)\end{array}$

Now we can write an expression of the fraction of the flux emitted by patch i that reached patch j as

$\begin{array}{cc}{\mathit{\text{F}}}_{\mathit{\text{i}}\to \mathit{\text{j}}}={\mathrm{\text{Φ}}}_{\mathit{\text{i}}\to \mathit{\text{j}}}/{\mathrm{\text{Φ}}}_{\mathit{\text{i}}}=1/\left(\mathit{\text{π}}{\mathit{\text{A}}}_{\mathit{\text{i}}}\right){\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}{\displaystyle {\int }_{{\mathit{\text{A}}}_{{}_{\mathit{\text{i}}}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2& \left(11.19\right)\end{array}$

The expression of this fraction contains terms that depend on the geometry and orientation and is independent of any lighting characteristics of the surface. This fraction is called form-factor, and plays an important role in radiosity computation.

We said earlier that every patch in a scene is an emitter, be it primary or secondary. So if patch i is emitting towards patch j then patch j also must be emitting towards patch i. Using derivations as above we can also write the expression for fraction Fj→i, the fraction of flux received by patch i due to emission from patch j. And the expression is:

$\begin{array}{cc}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}=1/\left(\mathit{\text{π}}{\mathit{\text{A}}}_{\mathit{\text{j}}}\right){\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{i}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}}}{\displaystyle {\int }_{{\mathit{\text{A}}}_{{}_{\mathit{\text{j}}}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2& \left(11.20\right)\end{array}$

We notice that both the fractions are similar except for the area term in the denominator of the right hand side. That means the fractions are related to each other and the relation is: AiFi→j = AjFj→i. This relation becomes useful in the derivation of the radiosity transport equation, and also is useful during form factor computation because if we know Fi→j then we can get Fj→i and vice versa by simply applying the relation.

11.2.4 Flux Transport Equation and Radiosity Transport Equation

In this section we derive an expression for total flux leaving any patch i. The total flux leaving a patch has two components: emission term, and a reflection term. Emission is a property of the surface. If the surface is a light source then it has nonzero emission, zero otherwise. We assume that ${\mathrm{\text{Φ}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}$, the emitting flux for each patch, is known a priori. The reflection term ${\mathrm{\text{Φ}}}_{\mathit{\text{i}}}^{\mathit{\text{r}}}$ is due to the reflection of the flux received by the patch. In the paragraph above, we defined form factor to express the fractional amount of flux received from elsewhere. Part of that received light is reflection. The fraction reflected is also a material property known as reflectance ρ, and is also known a priori for all patches. We use this information to write an expression of flux leaving any patch i in terms of flux leaving from all other patches in the scene as

$\begin{array}{cc}{\mathrm{\text{Φ}}}_{\mathit{\text{i}}}={\mathrm{\text{Φ}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}+{\mathrm{\text{Φ}}}_{\mathit{\text{i}}}^{\mathit{\text{r}}}={\mathrm{\text{Φ}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}+{\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle \sum _{\mathit{\text{j}}=1}^{\mathit{\text{N}}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}{\mathrm{\text{Φ}}}_{\mathit{\text{j}}}}& \left(11.21\right)\end{array}$

This equation is called flux transport equation.

In a uniformly emitting surface the flux is related to its radiosity by the equation Φ = BA. We now use this relation to substitute flux with radiosity to derive a transport equation for radiosity.

$\begin{array}{cc}{\mathit{\text{B}}}_{\mathit{\text{i}}}{\mathit{\text{A}}}_{\mathit{\text{i}}}={\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}{\mathit{\text{A}}}_{\mathit{\text{i}}}+{\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle \sum _{\mathit{\text{j}}=1}^{\mathit{\text{N}}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}{\mathit{\text{B}}}_{\mathit{\text{j}}}{\mathit{\text{A}}}_{\mathit{\text{j}}}.}& \left(11.22\right)\end{array}$

Dividing both sides by Ai we get

$\begin{array}{cc}{\mathit{\text{B}}}_{\mathit{\text{i}}}={\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}+{\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle \sum _{\mathit{\text{j}}=1}^{\mathit{\text{N}}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}{\mathit{\text{B}}}_{\mathit{\text{j}}}{\mathit{\text{A}}}_{\mathit{\text{j}}}/{\mathit{\text{A}}}_{\mathit{\text{i}}}.}& \left(11.23\right)\end{array}$

Next we use the relation between form factors to replace Fj→iAj/Ai by Fi→j to get:

$\begin{array}{cc}{\mathit{\text{B}}}_{\mathit{\text{i}}}={\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}+{\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle \sum _{\mathit{\text{j}}=1}^{\mathit{\text{N}}}{\mathit{\text{F}}}_{\mathit{\text{i}}\to \mathit{\text{j}}}{\mathit{\text{B}}}_{\mathit{\text{j}}}.}& \left(11.24\right)\end{array}$

This equation is called the radiosity transport equation. Note that both the transport equations are very similar. The only difference is the order in which the patch indices are specified in each equation’s form factor term. In the flux transport equation the order is from patch j to patch i, and in the radiosity transport equation it is the other way around.

Both the transport equations are written for a patch i and are valid for every patch i from 1 to N. So if we write expressions for all the is, then we get a system of linear equations that expresses flux of every patch in the scene in terms of flux of all other patches. We noticed earlier that form factor Fi→j (or Fj→i) depends on the orientation and location of the patches in the scene, and hence can be computed independent of the lighting in the scene. So we should be able to solve this linear system to get the equilibrium flux in the scene due to inter-reflection of the emitted light. Before we describe the radiosity solution methods, we describe a couple of methods for computing form factors.

11.2.4.1 Computation of Form Factor

Many methods have been proposed to compute form factor. Analytical methods of computing form factor between two visible patches have been proposed. However, in real scenes most patches are partly or fully occluded from each other. So visibility must be resolved before using any analytical method. As visibility computation is the most expensive part of the computation, and analytical methods are only useful for completely visible pairs of patches, the use of analytical methods is not popular. So here we discuss a couple of numerical methods that account for visibility during the computation. We mentioned earlier that we can take advantage of the fact that if we have computed Fi→j then Fj→i can be computed by simple arithmetic. So this principle itself re-duces the computation cost by half. We describe a couple of methods for the actual form factor computation.

Form Factor by Monte Carlo Quadrature: The form factor between two patches is a four-dimensional integration. So computation of this term by Monte Carlo quadrature proceeds as shown in Listing 11.13.

Fi→j = 0

for N times

 Let pi be a randomly sampled point on Patch i

 Let ni be the normal to the patch at pi

 Let pj be arandomly sampled point on Patch j

 Let nj be the normal to the patch at pj

 d = pj - pi

 R = |d|

 shadowRay = (pi,d)

 V = inShadow(ray, scene) ? 0:1

 ΔF = dot(d/R, nj) dot(-d/R, nj) (V/πR2)

 Fi→j += (ΔF/N)

endfor

Fj→i = Fi→j (Ai/Aj)

LISTING 11.13: Algorithm for computing form factor between two patches using Monte Carlo sampling.

As we notice here, we compute visibility between two points using the inShadow function described in our ray tracing section. So the algorithm is relatively straightforward. However, we must note that visibility by ray tracing may require full ray-scene intersection, and the form factor must be computed between every pair of the patches, so this method can be expensive.

Hemisphere/Hemicube Method: This is a class of methods that computes form factor between a patch and the rest of the patches in the scene. They are approximate in nature, because the computation is carried out at only one point (mostly the center point) of the patch. So the form factor of the patch is approximated to be the same as the form factor of the differential surface around the point and rest of the scene. As the receiver surface is a differential patch the form factor equation is simplified from a double integration to a single integration.

$\begin{array}{cc}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}\approx {\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}}}={\displaystyle {\int }_{{\mathit{\text{A}}}_{\mathit{\text{j}}}}{\displaystyle \int \mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2}}& \left(11.25\right)\end{array}$

Furthermore, patches emitting towards another must all lie on the upper hemisphere of the receiving patch. So if we are interested in computing form factor between all those emitting patches, then we should in principle be integrating only over the surface of the upper hemisphere. We write here an expression for such a form factor.

$\begin{array}{cc}{\mathit{\text{F}}}_{{\mathcal{\mathscr{H}}}_{\mathit{\text{j}}}\to \mathit{\text{i}}}\approx {\displaystyle {\int }_{{\mathcal{\mathscr{H}}}_{\mathit{\text{j}}}}{\displaystyle \int \mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{i}}}{\mathit{\text{V}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{ω}}}_{\mathit{\text{j}}}}\mathit{\text{d}}{\mathit{\text{ω}}}_{\mathit{\text{j}}}\approx {\displaystyle \sum _{\mathit{\text{k}}}\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{k}}}{\mathit{\text{V}}}_{\mathit{\text{i}},\mathit{\text{k}}}\mathrm{\text{Δ}}{\mathit{\text{ω}}}_{\mathit{\text{k}}}={\displaystyle \sum _{\mathit{\text{k}}}{\mathit{\text{V}}}_{\mathit{\text{i}},\mathit{\text{k}}}\mathit{\text{f}}\mathit{\text{a}}\mathit{\text{c}}\mathit{\text{t}}\mathit{\text{o}}{\mathit{\text{r}}}_{\mathit{\text{k}}}}}}}& \left(11.26\right)\end{array}$

where ℋj represents area of patch j visible to patch i through the unit hemisphere over patch i, and dωj replaces $\mathit{\text{cos}}{\mathit{\text{θ}}}_{\mathit{\text{j}}}\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}/{\mathit{\text{R}}}_{\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{i}}},\mathit{\text{d}}{\mathit{\text{A}}}_{\mathit{\text{j}}}}^2$. For a unit hemisphere, R=1 and cosθj=1 as well. So the differential areas over the hemisphere itself represent differential solid angles. However, it is not so when we replace hemisphere by hemicube in the numerical method discussed below. The interpretation of the visibility term is slightly different here: it represents whether the patch j is visible through the differential solid angle. The summation terms in the equations approximate the integration over ℋj by the summation of discrete subdivisions of the hemisphere, Vi,k is the visibility of patch j to patch i through the solid angle k. The term factork represents the analytical form factor of the portion of hemispherical surface area subtended by the k-th discrete solid angle. It is Δωj times the cosine of the direction through the k-th solid angle, and can be computed a priori for specific discretizations. This approach of form factor computation requires that we find the visibility of patch j through each discrete subdivision of the hemisphere. Instead of finding the discrete subdivisions through which a patch j is visible, if we find the patch j visible through each discrete subdivision k of the hemisphere then this latter becomes the already discussed problem of nearest object finding along a ray. And we can write an expression for computing only a fraction of the form factor for a part of the patch j visible through the discrete solid angle as:

$\begin{array}{cc}\mathrm{\text{Δ}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}\approx \mathrm{\text{Δ}}{\mathit{\text{F}}}_{{\mathit{\text{ℋ}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}}=\mathit{\text{f}}\mathit{\text{a}}\mathit{\text{c}}\mathit{\text{t}}\mathit{\text{o}}{\mathit{\text{r}}}_{\mathit{\text{k}}}& \left(11.27\right)\end{array}$

and compute the form-factor for the whole patch j as ${\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}={\displaystyle \sum \mathit{\text{δ}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}}$. Using these formulations we can write the algorithm for computing the form factor between the patches as given in Listing 11.14.

for all j

 Fj→i = 0

endfor

Let Hi be the unit hemisphere around p_i the center of patch i.

Let N be the number of discrete subdivision of the hemisphere

for k = 1 to N

 // assuming N is the number of discrete subdivision of the hemisphere

 // Let dk be direction from the center of patch i through the center of the k— th solid angle ray = (pi ,dk)

 j = nearestObject(ray,scene)

 Fj→i += factor_k

endfor

LISTING 11.14: Algorithm for computing form factor between a patch and all other patches using a method based on projection on hemisphere.

Notice that in this algorithm we have replaced the visibility function from the form factor equation by nearest object, that returns the patch visible along the direction through the solid angle. By doing so we are able to update the form factor between the receiver patch and all those patches contributing to the receiver patch. By discretizing the hemisphere sufficiently we can get a good-enough approximation of the form factors.

A variation of the hemisphere-based method is the Hemicube method, in which the hemicube replaces the hemisphere. Hemicube covers the whole hemisphere of directions, but has the advantage that its five faces are flat rectangles. So one may use hardware-based visibility computation by projecting the whole scene on each of the hemicube faces, and rendering patch id’s on the discrete hemicube pixels. Thanks to hardware acceleration techniques, this method used to be much faster compared to ray tracing, and hence was one of the earliest acceleration methods proposed for form factor computation.

11.2.5 Solution of Radiosity System

We now have all the components required to solve the linear system of radiosity transport and flux transport equation. We can write a matrix form MA = B where A represents the vector of unknown equilibrium flux or radiosity terms, B represents the vector of emitting flux or radiosity terms, and M represents a square matrix with terms containing known reflectance and computed form-factors. The size of the vectors and matrix depends on the number of patches in the scene. Solutions for the system A = M-1B can be computed by a matrix inversion followed by a matrix vector multiplication. For a vector size N, the computational complexity of this approach is O(N3). Because of this complexity, a radiosity/flux system solution for a reasonably complex scene (> 10000 patches) is impractical. One must note here that patches must be small, to satisfy the uniform brightness requirement. So even for a smallish scene the number of patches reaching or exceeding > 10000 is common. An alternate approach to radiosity/flux solution is to use an iterative method for solving linear systems. Any of the three well-known iteration methods, Jacobi method, Gauss-Seidel method and Southwell iteration, may be used. Independent of the exact method used, all the iterative methods start with an initial guess (see Listing 11.15) for radiosity/flux distribution, update radiosity of all patches based on the previous guess, and repeat until converged.

for patch i=1 to N // Initialize

 ${\mathit{\text{B}}}_{\mathit{\text{i}}}={\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}$

endfor

LISTING 11.15: Initialization for gathering based method.

We first describe Jacobi method (see Listing 11.16) and Gauss-Seidel method (see Listing 11.17)-based algorithms for solving the radiosity system. Same approaches may be applied to solving flux systems as well.

// Jacobi Iteration : Radiosity Gathering Method

while (not converged)

  for patch i=1 to N

  $Bnew_i = B_i^e + 
ho_i sum_j B_j F_{i 
ightarrow j}$

  endfor

  for patch i=1 to N

  Bi = Bnewi

  endfor

endwhile

LISTING 11.16: Jacobi-iteration-based method for computing equilibrium radiosity.

// Gauss‒Seidel Iteration: Radiosity Gathering Method

while (not converged)

  for patch i=1 to N

  ${\mathit{\text{B}}}_{\mathit{\text{i}}}={\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}+{\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle {\sum }_{\mathit{\text{j}}}{\mathit{\text{B}}}_{\mathit{\text{j}}}{\mathit{\text{F}}}_{\mathit{\text{i}}\to \mathit{\text{j}}}}$

 endfor

endwhile

LISTING 11.17: Gauss-Seidel-iteration-based method for computing equilibrium radiosity.

One may notice that the Jacobi and Gauss-Seidel methods are small variations of each other. The Gauss-Seidel method converges relatively faster, and does not require any intermediate vector Bnew. So of the two it is the preferred method.

The third iterative algorithm (see Listing 11.18), like photon tracing, simulates the natural process of light transport. The iteration method distributes light from the brightest patches in the scene. To start with, the emitter patches are the brightest. After light from emitters is distributed, the next brightest patch is chosen and its light that has not been distributed (shot) so far is distributed, and the process is continued till the equilibrium is reached.

// Southwell Iteration: Radiosity Shooting Method

for patch i=1 to N // Initialize

  B_unshot_i = ${\mathit{\text{B}}}_{\mathit{\text{i}}}^{\mathit{\text{e}}}$

endfor

while (not converged)

  // Let j be the patch with highest B_unshot

  for patch i=1 to N

  $\mathrm{\text{Δ}}\mathit{\text{B}}={\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle {\sum }_{\mathit{\text{j}}}}$ B_unshot_j Fi→j

  Bi += ΔB

  B_unshot_i += ΔB

  endfor

  B_unshot_j = 0

endwhile

LISTING 11.18: Southwell-iteration-based method for computing equilibrium radiosity.

This algorithm requires an additional vector B_unshot to keep track of the radiosity that is yet to be shot or distributed, and the total radiosity. Every time the light distribution from a surface patch, j, is completed its B_unshot is set to zero. This avoids erroneous repeated distribution of the same light from the patches.

All these iterative algorithms require an equilibrium condition check. The equilibrium condition is assumed to have been reached when the incremental flux on all the patches is reduced below a small threshold. In practice, the iteration is continued for a fixed number of iterations. In such a case, choosing a shooting-based method improves the accuracy of overall lighting computation. In addition to the accuracy-related criterion, this algorithm is also preferred because of a couple of other reasons: less storage requirement and adaptive form factor computation. Note that the form-factor matrix required N2 storage and a similar order of computation effort. However, we may notice that the iterative solution proceeds with one patch (the brightest one) at a time. So one may choose to compute only the form-factors from this patch to all other patches using a hemisphere or hemicube-based method, and avoid the N2 computation and storage. However, this approach requires a small modification to the ΔB computation step of the shooting algorithm. The modified step is

$\begin{array}{cc}\mathrm{\text{Δ}}\mathit{\text{B}}={\mathit{\text{ρ}}}_{\mathit{\text{i}}}{\displaystyle \sum _{\mathit{\text{j}}=1}{\mathit{\text{B}}}_{\mathit{\text{j}}}{\mathit{\text{F}}}_{\mathit{\text{j}}\to \mathit{\text{i}}}{\mathit{\text{A}}}_{\mathit{\text{j}}}/{\mathit{\text{A}}}_{\mathit{\text{i}}}.}& \left(11.28\right)\end{array}$

11.2.5.1 Rendering from Radiosity Solution

The radiosity method computes equilibrium radiosity over every patch of the scene independent of whether the patch is visible or not. For creating rendered images of the scenes, a rendering step must follow a radiosity computation. The goal in this step is to create a realistic-looking visualization of the scene from the computed solution. The simplest approach is to use the patch radiosity directly for color. Radiosity from a uniformly diffuse surface patch is related to its radiance by the relation: B = πL. This radiance can be used as the color for any point on the patch visible to the virtual camera through a pixel. Such a rendering unfortunately creates a faceted appearance.

A simple approach to get rid of this faceted appearance is to compute radiosity at the vertices of the patches by area weighted averaging of the radiosity values of the connecting patches, followed by Gouraud interpolated rendering. Though the resulting appearance is smooth, it may still suffer from inaccuracy in meshing and will result in blurring of highlights, and leakage of shadow and light. In the following section we discuss a gathering-based approach for creating high-quality rendering.

Final Gathering: In final-gathering-based rendering, radiance at the surface point visible though a pixel is computed by gathering light from the hemisphere surrounding the point. Either of the direct lighting computation algorithms, the one described in Section 2.2.3 for computing from an environment light source, or the one described in Section 2.2.2, can be adapted for use in final gathering.

When using the environment lighting-based approach, the sampling process remains unchanged. Each sampled ray from the point of interest is traced to find the visible surface patch. The radiance from the visible patches along the sample rays are averaged to get the estimated radiance from the point.

The radiosity method computes outgoing equilibrium radiosity at every surface patch in the scene. Thus, every patch in the scene may be considered as a potential light source. The method for direct light computation due to an area light source can be used to compute light at any point on the scene due to each patch. Thus, one can use an area light-based gathering approach for rendering realistic looking images.

One disadvantage to using each and every patch for lighting computation during the gathering process is that in a complex scene only a small fraction of the scene patches are visible from any point in the scene. The contribution of invisible surface patches to the illumination of the point is zero and thus any computational effort spent in gathering light from such patches is wasteful. A hemi-cube-based method removes this problem. In this method a virtual unit hemi-cube is set up over the point of interest. Using the hardware Z-buffer rendering method, the scene is projected onto the faces of this hemicube. The Z-buffer algorithm eliminates invisible surfaces. Each pixel on the hemicube represents a small piece of a visible surface patch and hence can be considered as a light source illuminating the surface point at the center of the hemi-cube.

If the pixel is sufficiently small then the direct lighting due to the pixel can be approximated by BjFpixel where Bj is the radiosity of the surface patch j projected on the pixel, ${\mathit{\text{F}}}_{\mathit{\text{p}}\mathit{\text{i}}\mathit{\text{x}}\mathit{\text{e}}\mathit{\text{l}}}=\frac{\mathit{\text{cos}}{\mathit{\text{θ}}}_1\mathit{\text{cos}}{\mathit{\text{θ}}}_2}{\mathit{\text{π}}{\mathit{\text{r}}}_{12}^2}\mathrm{\text{Δ}}\mathit{\text{A}},\mathrm{\text{Δ}}\mathit{\text{A}}$ is the area of the pixel and is mostly the same for all hemicube pixels.

The cosine terms and the r term in the Fpixel expression depend on the hemicube face on which the pixel lies and the coordinates of the pixel. For the top face whose pixel coordinate is (x,y,1), ${\mathit{\text{F}}}_{\mathit{\text{p}}\mathit{\text{i}}\mathit{\text{x}}\mathit{\text{e}}\mathit{\text{l}}}=\frac{1}{\mathit{\text{π}}{\left({\mathit{\text{x}}}^2+{\mathit{\text{y}}}^2+1\right)}^2}\mathrm{\text{Δ}}\mathit{\text{A}}$. For the side faces, the pixel coordinate is (±1,y,z) and ${\mathit{\text{F}}}_{\mathit{\text{p}}\mathit{\text{i}}\mathit{\text{x}}\mathit{\text{e}}\mathit{\text{l}}}=\frac{\mathit{\text{z}}}{\mathit{\text{π}}{\left(1+{\mathit{\text{y}}}^2+{\mathit{\text{z}}}^2\right)}^2}\mathrm{\text{Δ}}\mathit{\text{A}}$. For both the front and back side the expression with pixel coordinate (x,±1,z) is $\mathrm{\text{Δ}}{\mathit{\text{F}}}_{\mathit{\text{p}}\mathit{\text{i}}\mathit{\text{x}}\mathit{\text{e}}\mathit{\text{l}}}=\frac{\mathit{\text{z}}}{\mathit{\text{π}}{\left({\mathit{\text{x}}}^2+1+{\mathit{\text{z}}}^2\right)}^2}\mathrm{\text{Δ}}\mathit{\text{A}}$. In the gathering-based computation, only ra-diosity, Bj, is projection dependent. The rest of the terms in the equation can be pre-computed and multiplied at the time of final gathering. The final gathering method can be adapted to create noise-free rendering in the photon-tracing-based method. Unlike in the radiosity method, the photon-tracing-based method does not compute equilibrium radiosity. It stores the photon information at the hit points. So radiosity must be computed every time at the surface point visible through the hemicube pixel.