by Hin Jang
Radiosity is the rate at which energy leaves a surface. A Lambertian surface is one that reflects an amount of light from a unit differential area dA proportional to the cosine of the angle between the surface normal N and the direction of the light source2. For an environment consisting of such surfaces, the radiosity equation is
|\ | B(x) = Be(x) + p(x) | G(x, x')B(x') dx' | \| = Be(x) + p(x)E(x) = Be(x) + Br(x)
where B(x), the radiosity at point x measured in energy per unit time per unit area, is the sum of emitted radiosity B e (x) and the reflected radiosity B r (x). Reflectivity, a function of wavelength, is denoted p(x). The integral is taken over the hemisphere about x. G(x, x') accounts for the relative orientation, distance and visibility of the surfaces.
cos(tx) cos(tx') G(x, x') = -------------------- V(x, x') pi |x - x'|2
V(x, x') is equal to one if point x' is visible from point x, zero otherwise. E(x) is irradiance, the amount of energy per unit area received from other surfaces.
Solving the radiosity equation involves projecting the emittance and reflectance functions onto a set of basis functions. For radiosity B(x), its approximation is B^(x), a linear combination of n basis functions {Ni}i = 1, ..., n. The approximation for reflectivity p(x) is similarly defined.
n ---- B(x) ~= B^(x) = \ BiNi(x) / ---- i=1 n ---- p(x) ~= p^(x) = \ plNl(x) / ---- l=1
where B i and p l are the coefficients of the chosen orthonormal bases. The coefficients B i are the inner product i
> 6. As such, the approximation of B(x) can be written asn ---- B(x) ~= B^(x) = \ i> Ni(x) / ---- i=1
The solution to the original radiosity equation for receiving surface i is approximated by the following linear system
n ---- B(x) = Bi = Bei + \ KijBj / ---- j=1
where
|\ |\ | | Kij = | p^(x) dx | G(x, x')Nj(x')Ni(x) dx' | | \| \|
The reflected radiosity B r (x) is
n n ---- ---- Br(x) = p^(x)E(x) = \ plNl(x) \ EiNi(x) / / ---- ---- l=1 i=1
If the radiosity equation is decomposed as follows
- irradiance from scene radiosities
The operator G for a general function f is|\ | G(f)(x) = | G(x, x')f(x') dx' | \|so that
E(x) = G(B)(x) - radiosities from irradiances and reflectivities
The operator S isS(f)(x) = p(x)f(x)so that
Br(x) = S(E)(x)
the radiosity equation becomes
B(x) = Be(x) + K(B)(x) = Be(x) + S o G(B)(x)
and the linear system can be written as
n ---- Bi = Bei + \ KijBj / ---- j=1 n n ---- ---- = Bei + \ Sik \ GkjBj / / ---- ---- k=1 j=1
where
Kij = <S o G(Nj), Ni> Gkj = <G(Nj), Nk> Sik = <S(Nk), Ni> n ---- = \ pll Nk Ni> / ---- l=1
The decomposition into operators G and S allows for efficient representation of each operator individually 3.
Hierarchical Radiosity
For n elements in an environment, a linear system of n equations must be solved to yield the radiosity solution. This requires an algorithm to compute n * n coefficients representing the interaction of light energy between each pair of elements. To avoid the enormous computational complexity, the operations can be decomposed into n blocks for a given accuracy 9. In each block, the magnitude of interaction is about the same. The approach of hierarchical radiosity subdivides, recursively, each input surface into a set of subpatches until the measure of interaction is constant across a given subpatch 7. Each node in the hierarchy represents an area of the original surface. Two nodes are linked if the interaction between them can be computed to within some predefined accuracy.
The following code fragment establishes all linkages between initial patches p and q. FormFactor( ) computes the percentage of light interaction as the integral of G(x, x'), defined eariler, with respect to the area of the receiver patch, taking into account, also, the degree of occlusion. If either form factor is larger than the estimate Fe, the patch is subdivided into four new quadrilaterials. The subdivision is stored in quadtree data structure. Subdivide( ) returns false if the patch cannot be subdivided further, in that its area is less than Ae.
void Refine(Patch *p, Patch *q, double Fe, double Ae) { double Fpq, Fqp; Fpq = FormFactor(p, q); Fqp = FormFactor(q, p); if (Fpq < Fe && Fqp < Fe) Link(p, q); else { if (Fpq > Fqp) { if (Subdivide(q, Ae)) { Refine(p, q->ne, Fe, Ae); Refine(p, q->nw, Fe, Ae); Refine(p, q->se, Fe, Ae); Refine(p, q->sw, Fe, Ae); } else Link(p, q); } else { if (Subdivide(p, Ae)) { Refine(q, p->ne, Fe, Ae); Refine(q, p->ne, Fe, Ae); Refine(q, p->ne, Fe, Ae); Refine(q, p->ne, Fe, Ae); } else Link(p, q); } } }
Once all form factors have been determined, the radiosity for each patch is calculated. Gather( ) accumlates the total amount of energy received by a patch directly and from its parent subpatches. The average brightness of each patch is stored in B and its diffuse colour in is stored in Cd. The brightness, gathered from the list of all interactions in q, is stored in Bg.
void Gather(Patch *p) { Patch *q; double Fpq; if (p != NULL) { p->Bg = 0.0; for (q = p->interactions; q != NULL; q = q->next) { Fpq = FormFactor(p, q) p->Bg += Fpq * p->Cd * q->B; } Gather(p->sw); Gather(p->se); Gather(p->nw); Gather(p->ne); } }
The disadvantage of hierarchical radiosity is that shadow boundaries tend to be jagged because subdivision occurs regularly and, therefore, does not follow the contour of the shadow. One way to remove the discontinuities is to mesh the environment along the curve. This method of discontinuity meshing, however, increases the number of surfaces in a scene, magnifying the computational complexity.
Galerkin Radiosity
The accuracy to which radiosity is computed in the previous method is dependent on surface geometry. Complex interactions among non-planar surfaces require finer subdivision, at the cost of greater computation, since the bases are assumed to be piecewise constant across the subpatches. One way to avoid the limiting feature of hierarchical radiosity is to project the functions onto a higher order basis. Galerkin radiosity is a means by which the integral equation for radiosity can be solved in terms of a basis set of non-constant functions across a surface 12. The basis set is {Tk(s, t) | k = 0, 1 ... n}, where s and t are the parametric variables of the surface and k denotes a particular function of the set. {Tk(s, t)} is an orthonormal set of basis functions where the coefficients for a given radiosity function over surface i is
Bki = i, Tk>
so that the original function can be approximated with
n ---- Bi(s, t) ~= \ Bki Tk(s, t) / ---- k=1
Applying the Galerkin method for radiosity begins with the radiosity equation of two variables
|\ |\ | | B(s, t) = Be(s, t) + p(s, t) | | G(s, t, u, v)B(u, v) dudv | | \| \|
Expanding B(u, v) in terms of the basis set {Tl(u, v)} gives
|\ |\ | | B(s, t) = Be(s, t) + p(s, t)Bl | | G(s, t, u, v)Tl(u, v) dudv | | \| \|
By taking the inner product of both sides with the kth basis set function Tk(s, t), the radiosity equation can be written as a matrix equation
---- Bki - Eki = \ Blj Kklij / ---- j, l
where Kklij, Bki, and Eki are the form factors, patch radiosities, and emittances, respectively. The radiosity solution computed by this method is a list of basis set expansion coefficients Bki for each surface i and basis function k 12. The radiance at a point (s, t) on surface i is recovered from these coefficients using
n ---- Bi(s, t) ~= \ Bki Tk(s, t) / ---- k=1
Galerkin radiosity allows direct evaluation of the radiosity equation without the need to tesselate a curved surface. Meshing is only required when two surfaces are extremely close to each other and is not needed to model variations in intensity across a surface 12.
[1] Bastos, R., M. Goslin, and H. Zhang, Efficient Rendering of Radiosity Using Textures and Bicubic Reconstruction, TR-96-025, Department of Computer Science, University of North Carolina, Chapel Hill, 1996
[2] Foley, J.D., A. van Dam, S.K. Feiner, and J.F. Hughes, Computer Graphics Principles and Practice, Second Edition, Addison-Wesley, Reading, 723-724, 1990 [3] Gershbein, R., P. Schroder, and P. Hanrahan, Textures and Radiosity: Controlling Emission and Reflection with Texture Maps, Research Report CS-TR-449-94, Department of Computer Science, Princeton University, 1993 [4] Gershbein, R., An Adaptive Gauss Method For Computing Irradiance Coefficients of Galerkin Radiosity Systems, TR-485-95, Department of Computer Science, Princeton University, 1995 [5] Goral C.M., K.E. Torrance, D.P. Greenberg, and B. Battaile, "Modeling the Interaction of Light Between Diffuse Surfaces," Computer Graphics, 18(3):213-222, July 1984 [6] Gortler, S.J., P. Schr?der, M.F. Cohen, and P. Hanrahan, "Wavelet Radiosity," Computer Graphics, SIGGRAPH 1993 Proceedings, 27(4):221-230 [7] Hanrahan, P., D. Salzman, and L. Aupperle, "A Rapid Hierarchical Radiosity Algorithm," Computer Graphics, SIGGRAPH 1991 Proceedings, 25(4):197-206 [8] Kajiya, J.T., "The Rendering Equation," Computer Graphics, SIGGRAPH 1986 Proceedings, 20(4):143-149 [9] Lischinski, D., F. Tampieri, and D.P. Greenberg, "Combining Hierarchical Radiosity and Discontinuity Meshing," Computer Graphics, SIGGRAPH 1993 Proceedings, 27(4):199-208 [10] Pellegrini, M., "Monte Carlo Approximation of Form Factors with Error Bounded a Priori," ACM Proceedings of the Eleventh Annual Symposium on Computational Geomerty, 287-296, 1995 [11] Smits, B., J. Arvo, and D. Greenberg, "A Clustering Algorithm for Radiosity in Complex Environments," Computer Graphics, SIGGRAPH 1994 Proceedings, 28(4):435-442 [12] Zatz, H.R., "Galerkin Radiosity: A Higher Order Solution for Global Illumination," Computer Graphics, SIGGRAPH 1993 Proceedings, 27(4):213-220
