Ian Ashdown is president of byHeart Software Limited. He has a B.App.Sc. in electrical engineering from the University of British Columbia and 13 years of software engineering experience, and is currently Research and Development Manager for Ledalite Architectural Products Inc., a commercial lighting fixture manufacturer. He can be reached at 620 Ballantree Road, West Vancouver, B.C., Canada V7S 1W3, tel. (604) 922-6148, fax. (604) 987-7621 or CompuServe (72060,2420).
Introduction
Say "realistic computer-generated images" and most of us think of various ray-tracing techniques. True realism in computer-generated imagery, however, requires an accurate rendition of diffuse reflections, color bleeding between surfaces, penumbrae along shadow boundaries, and detailed shading within shadows. These are effects that ray-tracing techniques can generally model only by tracing enormous numbers of rays, with rendition times extending from hours to days on even the fastest desktop computers and workstations.If rendition time is an issue (and when is it not?), there is an alternative — radiosity. While radiosity in many ways complements ray tracing, radiosity accurately models area light sources, diffuse reflections, color bleeding, and realistic shadows. Ray tracing works best for modeling point light sources, specular reflections, and refraction effects. Radiosity used in combination with ray-tracing, or alone, can create stunningly realistic images in reasonable amounts of time.
Admittedly, the radiosity method has its limitations. It requires that all surfaces be ideal diffuse reflectors. Each surface must be subdivided into meshes of small planar polygons (usually rectangles and triangles) called patches, with each patch further subdivided into polygonal elements. (If these polygonal subdivisions are not carefully chosen, aliasing artifacts may appear in the rendered image.) The method typically needs several megabytes of RAM to serve as its working memory, and often requires as much CPU time as ray tracing to model a complex scene. By itself, the radiosity method cannot efficiently model specular reflections or transparent surfaces. Finally, it is strictly applicable to closed environments only, where every ray of light must eventually intersect a surface.
Despite these limitations, the radiosity method is a practical and useful image-synthesis technique, particularly for scenes such as architectural interiors where most surfaces are diffuse reflectors, light sources may include windows and skylights, and shadows must be faithfully rendered. If necessary, ray-tracing techniques can be used to model the relatively few specular surfaces typically present in such environments.
Radiosity has two distinct advantages over ray tracing. First, it generates view-independent solutions. That is, the method precalculates the luminance (or photometric brightness) of every surface element in the environment. Once it does so, a view of the environment can be calculated with relatively little effort for any synthetic camera position. Second, progressive radiosity, one form of the radiosity method, is iterative. It first estimates the luminance of every surface element, then progressively improves its estimates until they converge to a solution. This means that a geometrically-correct view of the environment can be displayed almost immediately, while the element luminances can be updated and refined after each iteration.
The radiosity method offers a powerful capability that ray-tracing techniques cannot — interactive walk-throughs of the environment. Given for example a 3-D architectural description of a building interior for which the element luminances have been precalculated, a user can visually "walk" through the rooms using a mouse or 3-D pointing device. Today's workstations are capable of displaying only coarse polygon renditions at interactive rates. Tomorrow's virtual-reality machines will be capable of much more, and will undoubtedly rely on the radiosity method to synthesize realistic images in real-time.
Like ray tracing, the radiosity method can be successfully implemented on a desktop microcomputer — a 486 system running Windows 3.x is more than adequate. Regrettably, it is not possible to present a complete implementation within the confines of this article, or to discuss the many improvements to the original method that have been developed in the past few years. What is presented here is the basic theory of radiosity and C code for a very simple radiosity-method demonstration program. For readers interested in implementing their own radiosity-based render, a comprehensive bibliography in included on the Code Disk.
Radiosity Theory
The radiosity method was introduced to the computer graphics community by Goral et al. (1984), who based their work on radiative heat-transfer theory and thermal engineering (e.g., Siegel and Howell, 1981). Their paper was quite technical and relied on the terminology of thermal engineering: radiative heat-transfer equations, finite-element analysis, double-area integrals, Stoke's Theorem, and so forth. Fortunately, the method can be described in a much more intuitive but still rigorous manner.Consider the room with a table and a ceiling light shown in Figure 1. All of the surfaces are ideal diffuse (or Lambertian) reflectors, which reflect an incident beam of light illuminating the surface equally in all directions. (If the surface is self-luminous, it emits as well as reflects light equally in all directions.) Each surface is subdivided into a mesh of rectangular elements.
The only self-luminous elements in the room are those of the ceiling light. Light (or luminous flux) emitted by these elements will directly illuminate the floor, walls, and table top. The sides of the table will be in shadow, and the ceiling elements will not be directly illuminated. Depending on the reflectance of each element, some light will be reflected back into the room from the directly illuminated elements and the remainder absorbed. The reflected light will illuminate other surface elements, including those of the table sides and ceiling. Again, some of this light will be reflected back into the room and the remainder absorbed. This process repeats until all of the light is eventually absorbed. It can be extended to colored surfaces by considering the reflectance of each surface in three or more spectral bands (e.g., red, green, and blue) and independently calculating the light (now called spectral radiant flux) reflected between elements in each band.
The process is clearly iterative, and can be applied without modification to any closed environment. If a record is kept of how much light each element reflects and/or emits, its luminance is given by
L = F/pA
where F is the total amount of flux emitted and/or reflected by the element, and A is its area. Having the element luminances, you can use a 3-D graphics package to display a view of the room from any position. As far as the graphics package is concerned, the room is no more than a collection of small polygons. (Most computer-graphics textbooks offer discussions of the algorithms needed to display 3-D polygons, e.g., view transformation, clipping, hidden-surface removal, backface culling, scan conversion, and Gouraud shading. These algorithms are required by, but not part of, the radiosity method.) Since the luminance of each element (polygon) is independent of the position and orientation of the synthetic camera, the solution is view-independent.
The original approach (now referred to as the full radiosity method) presented by Goral et al. (1984), expressed the above process, using finite-element analysis techniques, as a system of linear equations in an n-by-n matrix, where n is the total number of elements. The element luminances were then determined by solving the system of equations. Unfortunately, this method requires O(n2) memory and processing time. An environment with a mere 1,000 elements requires at least 4MB of memory and hours of processing time.
Progressive Radiosity
A better approach — the progressive radiosity method, introduced by Cohen, et al. (1988) — reduces memory requirements and improves processing time significantly. You define three variables for each element Ei: the element luminance (Li), the total amount of reflected and/or emitted flux (Fitotal) and the amount of flux that has been received by the element but not yet reflected (Fiunsent). In addition, you have several constants: the flux initially emitted by self-luminous elements (Fiemit), the element area (Ai) and the user-specified minimum amount of flux (Fmin). Elements which are not self-luminous have a Fiemit value of zero, and the process terminates when the quantity of unsent flux for every element is less than Fmin . Figure 2 shows the process in pseudocode.Progressive radiosity method differs from the full radiosity method in that it only requires O(n) memory. Instead of the minimum 4MB required by full radiosity, the progressive radiosity method can render an image of a 1,000 element environment using less than 100KB of memory. It also has the advantage of being explicitly iterative — the estimated element luminances are available for display purposes after each pass through the DO-WHILE loop. While the progressive radiosity method is no faster than full radiosity in computing the final element luminances, a useful image is available in seconds. (Selecting the element with the greatest unsent flux maximizes the rate of convergence.)
Apart from the somewhat cryptic line send Fiunsent to other elements, the above pseudocode is very easy to implement.
Form Factors and Hemi-Cubes
Figure 3 shows the simplest case of two elements E1 and E2, where E1 sends its (reflected or emitted) flux into space and E2 intercepts a portion of it. The fraction of flux received by E2 to that leaving E1 is called the form factor F12. Given F1unsent and the reflectance R2 of element E2, the amount of flux received by E2 is F1unsent F12, and the amount of flux subsequently reflected by E2 is F1unsent F12 R2.Determining F12 for two arbitrary elements is not an easy task. While analytic formulae have been developed for a number of simple geometries such as parallel and perpendicular rectangles and circular disks, (e.g., Siegel and Howell, 1981), most instances require solving a double-integral equation using contour integration. The problem is compounded for occluded environments, where E2 may be partially obscured by an intervening element E3 when viewed from E1 and thus have an apparently complex geometry.
Determining form factors in occluded environments using contour integration is a numerically complex and computationally demanding problem. Fortunately, approximate form factors can be computed to any desired degree of precision using the hemi-cube method introduced by Cohen and Greenberg (1985).
Figure 4 shows two elements E1 and E2, where the sending element E1 is much smaller than E2 and has been placed at the center of an imaginary unit hemi-cube (half a cube, similar to a hemisphere). A third element E3 has been placed behind E2. The size, shape, location, and orientation of E3 is such that E2 exactly obscures it as seen from E1. The outlines of both E2 and E3 as seen from E1 are projected onto the surface of the hemi-cube. Since E2 exactly obscures E3, the two outlines are identical. This also means that the two form factors F12 and F13 are identical, thus demonstrating that the distance of an element Ex from E1 doesn't matter when calculating the form factor F1x. What does matter is the projection of Ex onto the hemi-cube surrounding E1.
It is relatively easy to calculate the projection of a 3-D polygon onto a flat surface such as a hemi-cube face — this is exactly what you do when you perform a perspective projection of a polygon onto the viewing plane in a 3-D graphics program. Again, the necessary algorithms are discussed at length in most computer-graphics textbooks.
Suppose that you subdivide each face of the hemi-cube into an array of square elements (called patches) and scan convert the projection of the 3-D polygon onto each face. That is, you consider each patch of a hemi-cube face as being equivalent to a pixel in a frame buffer and discretize the projection. The delta form factor DF of a 3-D polygon which projects onto exactly one patch of a hemi-cube face is then very easy to calculate. For the top face, it is (Cohen and Greenberg 1985)
DF = DA/p (x2 + y2 + 1)2
and for any side face
DF = DAz/p (y2 + z2 + 1)2
where DA is the patch area. These formulae are strictly accurate only when the size of the elements is much smaller than the distance between them. In practice, the top face of the hemi-cube is usually subdivided into 50-by-50 to 100-by-100 patches.
Form factors are additive. That is, if the projection of a polygon covers a number of hemi-cube patches on one or more faces, the polygon from factor is simply the sum of the delta form factors of the covered patches. The delta form factors are constant for any given hemi-cube, and can therefore be calculated during program initialization and stored in memory. (The symmetry of the hemi-cube is such that, if you subdivide the top face into n-by-n patches and the side faces into n-by-n/2 patches, only 3n/8 delta form factors need to be stored.)
The progressive radiosity method pseudocode line send Fiunsent to other elements can thus be expanded as shown in Figure 5, where Rj is the reflectance of element Ej. (Note that placing the hemi-cube over an element involves positioning the hemi-cube origin at the geometric center of the element and aligning the hemi-cube z axis with that of the element normal. The x and y axes can be arbitrarily oriented with respect to the global coordinates.)
Projecting polygonal elements onto the hemi-cube consumes most of the CPU time required by the radiosity method. Recker et al. (1990) presented a "modified hemi-cube" for progressive-radiosity calculations which is on average twice as efficient as the original hemi-cube method. Further improvements can be achieved through careful analysis of the polygon clipping, hidden-surface removal, and scan-conversion algorithms used by the 3-D graphics package.
Substructuring Techniques
One problem with the hemi-cube technique outlined above is that a complex environment may have ten thousand or more elements. This implies that each iteration of the progressive radiosity method must calculate as many as 100 million form factors (although hidden-surface removal will likely reduce this number by an order of magnitude).Even in simple environments, a large number of elements are generally required to determine the fine surface-luminance details, particularly at shadow boundaries. The hemi-cube method also requires that the size of the sending element be small compared to the distance to the receiving element. If either of these conditions is violated, the outlines of the polygonal elements may be visible in the rendered images.
This problem can be overcome by using a hierarchy of elements. Cohen et al. (1986) proposed a two-level substructuring technique in which surfaces are subdivided into large polygonal patches (not to be confused with hemi-cube patches) and each patch further subdivided into elements. For the progressive-radiosity method (Cohen et al., 1988), luminous flux is sent from each patch (not element) to all other elements. The result is a marked decrease in the number of form factors to be calculated.
Various other substructuring techniques have been proposed in the literature (see the bibliography). Of these, the rapid hierarchical-radiosity algorithm proposed in Hanrahan et al. (1991) is particularly noteworthy. It is an optimized form of the progressive-radiosity method with an integral adaptive-substructuring technique. The algorithm automatically subdivides the surfaces into a hierarchy of patches and elements that minimizes the number of required form factor calculations while maintaining a user-specified error criterion.
Adaptive-substructuring techniques relieve the user of having to subdivide surfaces by hand within a 3-D CAD program. However, there are a number of issues which most of these techniques do not address, yet which affect the accuracy of the final image. They are discussed at length in Baum et al. (1991), an excellent paper which also describes an intelligent substructuring program for preprocessing AutoCad files.
Ray-Traced Form Factors
Another problem with the hemi-cube method is that it is prone to aliasing when determining form factors in complex images. The problem can be alleviated by increasing the number of patches on the hemi-cube. A more efficient technique proposed by Wallace et al. (1989) uses ray tracing to determine form factors by numerical integration. Unfortunately, the details of the technique are beyond the scope of this article.Shirley (1991) presented a probabilistic ray-tracing technique which implicitly calculates form factors. The luminous flux reflected from an ideal diffuse surface has a cosine distribution about the surface normal, in accordance with Lambert's Cosine Law. Using the progressive-radiosity method, the unsent flux from a sending element is divided into equally-sized "packets" and assigned to rays that are sent in random directions according to a cosine density distribution. The flux in each packet is transferred to the element that intercepts the ray. While not particularly efficient, the technique is very easy to implement, as demonstrated by RAY_RAD (see Listing 1) .
RAY_RAD — A Simple Radiosity Renderer
RAY_RAD is arguably the simplest demonstration of the progressive-radiosity method that can be written. It achieves this dubious honor by avoiding 3-D graphics algorithms wherever possible. For example, RAY_RAD considers only one environment: an empty room with a single ceiling light. Since every surface is fully visible to every other surface, there is no need to perform hidden-surface calculations. All surfaces are aligned with the coordinate system axes, and so 3-D view transformation calculations can be simplified. The program displays its results as tables of element luminances; therefore, there is no need to specify a camera position and orientation or to perform 3-D clipping, polygon scan conversion, or shading calculations. Following a program outline suggested by Shirley (1991), RAY_RAD uses probabilistic ray tracing to implicitly determine its form factors, and so does not need to project surface elements onto a hemi-cube. The result is a program which demonstrates the progressive radiosity method and no more. If you are contemplating writing your own radiosity renderer, please do not use RAY_RAD as a starting point!In RAY_RAD, main essentially follows the progressive-radiosity method pseudocode presented in Figure 2 and Figure 5, except that it loops through every element rather than selecting the one with the greatest unsent luminous flux. This makes it easier to display the results (by calling display_surface) after each iteration. send_flux calls select_ray to select a ray with a random origin (to minimize aliasing effects) within the sending element and a random direction selected according to a cosine distribution in the element's vertical plane, then calls find_element to determine which element the ray intercepts before updating the total and unsent flux values of the intercepting element. The functions local_to_global and global_to_local transform the ray coordinates between the global and local (element) coordinate systems.
Conclusions
In less than a decade, the radiosity method has gone from an academic curiosity to a useful computer-graphics technique. Research into improving the method and its associated sub-structuring techniques is continuing. Given its capabilities, it is surprising that the method has not received more attention outside the academic community. Most of the original academic papers have been published in the annual ACM SIGGRAPH conference proceedings, which are not widely available. Some of the more interesting papers are available only as university technical reports.Hopefully, this situation will change soon. Watt (1990) offers a very accessible discussion of full and progressive radiosity, while Chen (1991) has released C source code for a simple but effective progressive-radiosity renderer into the public domain. (It requires a 3-D graphics library to render polygons, but is otherwise complete.) Readers interested in developing their own renderers are referred to the annotated bibliography (included with the source code listings on the Code Disk). As for myself, I am currently writing a book on the radiosity method. Interested readers are invited to write, call, or fax me with any questions or suggestions.
References
Arvo, J., ed., 1991. Graphic Gems II. Academic Press, San Diego. CA.Baum, D.R., S. Mann, K.P. Smith and J.M. Winget. 1991. "Making Radiosity Usable: Automatic Preprocessing and Meshing Techniques for the Generation of Accurate Radiosity Solutions." Computer Graphics 25:4 (Proc. ACM SIGGRAPH '91), 51-60.
Chen, S.E. 1991. "Implementing Progressive Radiosity with User-Provided Polygon Display Routines." Arvo 295-298, 583-597.
Cohen, M.F. and D.P. Greenberg. 1985. "The Hemi-Cube: A Radiosity Solution for Complex Environments." Computer Graphics 19:3 (Proc. ACM SIGGRAPH '85), 31-40.
Cohen, M.F., D.P. Greenberg, D.S. Immel and P.J. Brock. 1986. "An Efficient Radiosity Approach for Realistic Image Synthesis." IEEE Computer Graphics and Application 6:3, 26-35.
Cohen, M.F., S.E. Chen, J.R. Wallace and D.P. Greenberg. 1988. "A Progressive Refinement Approach to Fast Radiosity Image Generation." Computer Graphics 22:4 (Proc. ACM SIGGRAPH '88), 75-84.
Goral, C. M., K.E. Torrance, D.P. Greenberg and B. Battaile. 1984. "Modeling the Interaction of Light Between Diffuse Surfaces." Computer Graphics 18:3 (Proc. ACM SIGGRAPH '84), 213-222.
Hanrahan, P., D. Salzman and L. Aupperle. 1991. "A Rapid Hierarchical Radiosity Algorithm." Computer Graphics 24:4 (Proc. ACM SIGGRAPH '91), 197-206.
Recker, R.J., D.W. George and D.P. Greenberg. 1990. "Acceleration Techniques for Progressive Refinement Radiosity", Computer Graphics 24:2 (1990 Symposium on Interactive 3D Graphics), 59-66.
Shirley, P. 1991. "Radiosity Via Ray Tracing." Arvo 306- 310.
Siegel, R. and J.R. Howell. 1981. Thermal Radiation Heat Transfer. Hemisphere Publishing, Washington DC.
Wallace, J.R., K.A. Elmquist and E.A. Haines. 1989. "A Ray Tracing Algorithm for Progressive Radiosity." Computer Graphics 23:3 (Proc. ACM SIGGRAPH '89), 315-324.
Watt, A. 1990. Fundamentals of Three-Dimensional Computer Graphics. Addison-Wesley, Reading, MA.