Computer Graphics & Geometry

One Method to Reduce Complexity of Matrix Radiosity Algorithm

Victor A. Debelov, Igor M. Sevastyanov
Institute of Computational Mathematics 
and Mathematical Geophysics
(former Computing Center) SB RAS,



The given report is devoted to an original approach for reducing computational cost (memory and time) of the matrix radiosity algorithm. Such investigations are still quite actual because the algorithm requires significant amounts of memory and time. From other side it began to be applied in practice due to increasing power and memory of modern desktop computers. Recent efforts of computer graphics specialists were done in the following main directions.

All those approaches try to process the whole scene geometry in that or another way, and thus they do not simplify the original problem dramatically. Somewhat aside the method of imposters lies. It bases on geometric subdivision of an initial scene, but farther parts of a scene do not take part in calculations of light energy equilibrium, therefore it is lacking of photorealism. The given here method of geometric decomposition � it is a quite new modification of the matrix radiosity algorithm. The scene is divided into parts (note, that it is done formally, i.e., without preliminary geometric analysis). Each part is processed independently on each other. Next step � exchange of energy between parts. It is easy to see that an independence of parts allows: Note that in the given report we consider the just geometric decomposition of a scene in opposite to "decomposition" used in classical problems of computational mathematics. We show feasibility of the given approach for the case, when the scene is divided into parts by plane(s), really geometric bisection. The formulation and numerical experiments are presented too. Keywords: synthesis of photorealistic images, radiosity equation, geometric decomposition of a scene, collocation method, parallel computing.
  1. Introduction.

      Firstly we will spend some time to introduction in a problem of solution to the radiosity equation. Next chapter is devoted to reformulation of the equation with respect to goals of the given report.

    1.1. Radiosity Equation.

      In the given text we consider the following Fredholm integral equation: (1)

    where the kernel  

    S- a scene, piecewise smooth surface, Le(x)- emission of light energy in 3D point , k(x)- diffuse reflectance function, V(x,y)- visibility between two points, i.e., it is 1 if open straight-line segment (x,y) has no intersections with S, and 0 otherwise, Nx,Ny- normals to a scene surface in corresponding points, L(x)- unknown radiosity. The problem is to determine L(x) for all points of S. The more detailed problem statement can be found in [1,2]. The main goal of the given report is to show the feasibility of the method of geometric bisection of a scene by an arbitrary plane. We expose our idea using collocation method

      1.2. Collocation method.

    Here we present the method in brief in limits necessary to support a presentation. Let us consider the reflectance equation (1) anew. A scene is divided into a set of nonintersecting finite elements Then we choose the basis: An approximating solution is found as: in a fixed collection of scene points Really the problem is reduced to a solution of a system of linear algebraic equations with respect to unknowns . In matrix notation it looks like: . As a rule a matrix of a system becomes very huge if we want to obtain images of high quality. Thus main efforts of investigations were directed to fighting with the size of  F.

    1.3. Hierarchical Radiosity.

      A hierarchical method [3] is a modification of a classical radiosity method; it was created in order to overcome a dramatic growth of a number of finite elements, and as a consequence, a huge order of a linear system. Hierarchical method fights against such disaster of radiosity methods via restraining the amount of finite elements in acceptable limits. Some hierarchical algorithms proceed in a top-down manner, by limiting the subdivision of a scene surface to what is necessary. Algorithms of hierarchical radiosity were derived from the N-body problem [4]. They first deal with a coarse subdivision and then step by step refine it in order to catch such scene particularities as shadow boundaries, penumbral areas, etc. Note that a hierarchical algorithm does not reduce the complexity of the initial set of finite elements. A hierarchy of links between finite elements of different subdivision levels replaces the matrix. Indeed hierarchical algorithms do not improve the situation while an initial number of finite elements is already large.

    1.4. Cluster Radiosity.

      The cluster method [5] is an evolution of the previous one. It was developed in order to reduce a computational cost by the grouping finite elements into clusters and approximating the light energy transfer between them. The very detailed and representative material on cluster radiosity method one can find in [6,7]. We should notice that cluster methods try nevertheless to process the whole scene. Several modifications assume so-called preprocessing step to prepare clusters. Anyway a final solution requires considerations of all scene clusters. Cluster methods allow creating hierarchies of links as top-down (similar to hierarchical methods) as down-top via grouping of small finite elements into clusters. Thus the growth of the size of matrix is restrained. It is too problematic to select scene clusters automatically, although some advances in this way are observed, e.g., see [7]. Obviously one could thoroughly prepare clusters in advance; nevertheless the problem of light energy transfer between clusters still persists because the decreasing of information leads to the loss of small details of image.

    1.5. Progressive Radiosity
         In the recent times the techniques related to calculations of global illumination began to be used in an         industry of electronic entertainments and as a standard tool of conventional 3D editors like 3D Studio. For example, the well-known computer game Quake allows precomputing of global illumination for game levels. It is camera-independence feature that permits to do such calculations only once. One more good example is the usage of radiosity in Lightscape3.2 [8] of Discreet Company (we mean plug-in for 3D Studio). Here we can find a different approach ("interactively driven") for the solution to the radiosity equation, which is called progressive radiosity. Roughly it may be characterized as follows.
  1. At initial step a scene surface is divided into quite big finite elements. Further they may be subdivided automatically into smaller ones in regions, where the significant difference of a light intensity between adjacent elements is found, e.g., across shadow boundaries.

  2. The light energy is transferred from each emitting element to all scene elements. Really it means that we compute a matrix row, which corresponds to emitting element. It is the stage of direct illumination.

  3. The progressive radiosity algorithm examines all surface elements in order to find the one, which has maximum radiosity. It is considered as an emitter. This step (indirect illumination) is then iteratively repeated.

  4. The process ends when a specified part of light energy is distributed among scene elements, i.e., the system reaches the equilibrium state.

As a rule the process of distributing of light energy is done with simultaneous displaying of the current result on the screen (the result after current iteration). It allows user to interrupt it if an image satisfies him. A remarkable investigation is demonstrated in a work [7]. The author combines the advantages of as cluster as progressive radiosity algorithms. However his algorithm takes into account the whole scene also. Obviously we say about a step of energy distribution only. The idea of this algorithm is too attractive because of it demonstrates the techniques how to approach to the required image step by step without computing all elements of the matrix. Actually in each step of the algorithm a scene is considered as consisting from 2 parts: 1) element, emitting the light into scene; 2) all other finite elements.

  1.6. Imposters

The radiosity method allows displaying on a screen all changes in scene illumination, thus it is quite suitable in such applications as urban navigations of cars [9]. Geometry of a city is nonvarying but its illumination depends on daytime, weather, etc. We may explain it in the following way. Let geometry be divided into blocks (similar to 3D cubes), each block is considered as a separate scene to which the radiosity method is applied. In order to simulate farther views imposters were used. Each imposter is a precalculated image mapped onto a corresponding wall of an imaginary cube. Obviously the imposters algorithm requires preliminary calculations. We do not consider here the advances of this approach but idea of imposters show us one of the ways to omit one part of a scene while calculating light energy equilibrium in other scene part. 

2. Geometric Decomposition


Fig.1. The idea of scene bisection 

All the methods mentioned above are derived from the classic matrix radiosity method. Their driving formulas are constructed starting from consideration of the full matrix of form-factors. Really the efforts were done to reduce amount of form-factors that should be actually computed. As a fact the matrix is not used explicitly. We decided to return to the classic matrix radiosity algorithm. We consider the fixed subdivision of a given scene into a set of finite elements. In other words the algorithm we suggest does not change matrix elements while solution thus preserving accuracy provided for by particular subdivision. This feature differences it from other radiosity algorithms because it preserves all scene details. Let us consider a fictitious plane in a scene, which obviously divides a scene into two parts � geometric bisection of a scene. Further we consider both scene parts separately, and independently on each other. Now we compute equilibrium of light energy in each part (independently!). A fictitious plane is a means to exchange light energy between scene parts. The process of energy redistribution is repeated cyclically. It is not important which methods of energy distribution are used for different scene parts. It means that for first scene part we could apply hierarchical algorithm, and for second part the progressive one. Moreover the choice of a method could depend on geometric and reflectance features of scene objects in the corresponding parts. Besides, we could use different radiosity algorithm in a different iteration steps. Honestly speaking it is the main direction of our investigations. In the given report we will illustrate our approach basing on the simple matrix radiosity algorithm. By the application of the approach to scene parts we could decrease the computational cost applying scene bisections recursively. Thus the main advantage of the given algorithm is a possibility to represent a calculation of radiosity equation for a complex scene as a sequence of more simple tasks. Obviously the approach should be a base of a parallel algorithm to solve the radiosity problem.

2.1. Derivation of driving formulas.

  Let us consider the mathematical formulations of the algorithm of the geometric bisection. A scene surface is divided into two parts:

Write the equation (1) in the following form: 


Then we rewrite it as a system 




Thus the system has a view: 

The approximating solution is given by the following iteration scheme: 



Terms and

define the exchange of light energy between scene parts, separated by fictitious plane.

2.2. Implementation

After division of a scene into two parts each finite element belongs to one of parts. Using appropriate enumeration we can write: and. Classic matrix radiosity algorithm requires computing the following table of form factors: , or the same in a short form:

The suggested algorithm work as follows.

  1. Computing of the table F11 for the scene part D1
  2. Computing of the table F22 for the scene part D2.
  3. Computing of the table F12, which is used while energy transfer from D2 to D1.
  4. Computing of the table F21, which is used while energy transfer from D1 to D2.
  5. Calculation of equilibrium of light energy in scene part D1.
  6. Calculation of equilibrium of light energy in scene part D2.
  7. Fulfill energy transfer from the scene part D1 to D2, using matrix  F21.
  8. Fulfill energy transfer from the scene part D2 to D1 using matrix F12.
  9. Single iteration of the algorithm consists of steps 5 � 8.
In the given implementation the iterative process assumes simultaneous display of current state of calculations (current state of an image), which allows user to stop iterations interactively after he gets necessary quality of an image. It is obvious that steps 1 � 4 can be done while preprocessing stage as the most of algorithms do. Calculations of both form factors tables of different scene parts in steps 1 and 2 can be done in different processors in parallel. Steps 5 and 6 (the light distribution inside each scene part) can be done in parallel too. Many of parallel computer architectures allow parallel processing of steps 1, 2, 3, and 4. And parallel processing of steps 7 and 8 is possible also. The computational requirements (an order of a system of linear equations, memory, and time) of the bisection algorithm in each its step are less than ones of the classic matrix radiosity algorithm. Moreover the given algorithm gains even if it is fulfilled in a single processor, because it requires at least less memory in each step.

2.3. Experiments

  In the presented test scenes we use finite elements � 3D rectangles. We do not apply any interpolation in order to smooth pictures; we try to show the algorithm behavior using coarse finite elements. Many illustrations of the given report are omitted in the text because of the absence of color possibilities in the proceedings.

2.3.1. "Constructive Wood"

  The idea of this test scene was taken from [1]. The scene consists of several two-sided diffuse planes, the only one of which is emitting light; see fig. 1. If this scene would be rendered by conventional ray tracing algorithm (not Monte Carlo), then its image would contain only four vertical bright strips � direct observation of a light source. The geometry of finite elements is shown in fig. 3. 


Fig. 2. Top view of scene geometry 


Fig. 3. Wireframe view of scene finite elements 

We used vertical plane, which divided the scene into two symmetrical parts as shown in fig. 1. The result shown in fig. 4 confirms that our algorithm correctly processes diffuse interreflections.


Fig. 4. Image of "constructive wood"

2.3.2. Shadows


Fig. 5. Geometry of a test scene (rough subdivision) 

The second test scene should demonstrate the work of the bisection algorithm iteration by iteration. As it is shown in fig. 5 the scene is divided by plane into two parts, which have different subdivision into finite elements. Only left scene part contains a light source (see fig. 6, 7) and a column, which are reasons of a shadow. Below we present visual results of a series of iterations. Fig. 6 demonstrates that right part of a scene is pure dark. It is obvious result because initially this scene part has no light sources. Figures 6, 7, and 8 show how energy transfer between parts influences on the current image from iteration to iteration. The shown three iterations demonstrate convergence of the image to the image rendered by classic matrix radiosity algorithm (see fig. 9).


Fig. 6. After step 6 on first iteration


Fig. 7. After step 6 on second iteration 


Fig. 8. After step 6 on third iteration 

In a conclusion let us look at times spent in different steps of the algorithm. First table is devoted to steps that are fulfilled at preprocessing stage. And the second one � times of steps of iterations.


Fig. 9. Classic matrix radiosity

Scene \ Step 1 2 3 4
Constr.wood 3 4 9 9 26 23 9
Shadow 10 11 40 40 101 91 40
Where: is time while sequential processing of steps; - time while steps 1 and 2 are done in parallel; - full parallelism. It is expected that is time of classic matrix algorithm nevertheless the algorithm of geometric bisection requires less memory.
Scene \ Steps 5 6 7 8
Constr.wood 16 16 0..5 0.5 33 17
Shadow 40 56 2 1 99 59
The given table shows that the bisection algorithm gains time almost twice with respect to the classic sequential matrix algorithm.

  3. Conclusion

Let us note the advantages of the algorithm of geometric bisection of a scene while solution of the radiosity equation.

. The presented numerical experiments have shown the feasibility and usefulness of the bisection radiosity algorithm. While sequential processing it reduces memory requirements. It looks similar to cluster algorithms except the fact that clusters (scene parts) are created automatically. Possibility to assemble a scene from precomputed parts.


The given work was supported by Russian Foundation for Basic Research under grants № 99-01-00577 and № 99-07-90422ск. 


  1. Cohen, M.F., Wallace, J.R. Radiosity and Realistic Image Synthesis. Academic Press, 1993, 379p.

  2. Heckbert, P.S. Finite Element Methods for Radiosity. In "Global Illumination" course, SIGGRAPH'93. See also

  3. Hanrahan, P., Salzman, D., Aupperle, L. A rapid hierarchical radiosity algorithm. Computer Graphics 25, 4 (July 1991), p. 197-206.

  4. Appel, A.A. An efficient program for many-body simulation. SIAM J. Sci. Stat. Computing 6(1), 1985, p.85-103.

  5. Willmott, A.J., Heckbert, P.S., Garland, M. Face Cluster Radiosity. Eurographics Workshop on Rendering, Granada, Spain, June 1999.

  6. Willmott, A.J., Heckbert, P.S. An Empirical Comparison of Radiosity Algorithms. Technical Report CMU-CS-97-115. See also
  7. Willmot, A.J. Hierarchical Radiosity with Multiresolution Meshes. Ph.D. thesis proposal. 206 p. Draft 17/6/2000.

  9. Sillion, F.X., Drettakis, G., Bodelet, B. Efficient Impostor Manipulation for Real-Time Visualization of Urban Scenery Proceedings of Eurographics'97, Budapest, Hungary, September 4-8, 1997.

Computer Graphics & Geometry