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,
E-mail:

• Abstracts
• 1. Introduction
  • 1.1. Radiosity Equation.
  • 1.2. Collocation method.
  • 1.3. Hierarchical Radiosity.
  • 1.4. Cluster Radiosity.
  • 1.5. Progressive Radiosity
  • 1.6. Imposters
• 2. Geometric Decomposition
  • 2.1. Derivation of driving formulas.
  • 2.2. Implementation
    • 2.3.1. "Constructive Wood"
    • 2.3.2. Shadows
• 3. Conclusion
• Acknowledgements
• References

Abstract  

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.

• Adaptive scene subdivision (hierarchical radiosity methods) was elaborated in order to prevent the growth of the number of finite elements in acceptable limits.

• Grouping of scene objects � cluster methods, when the light energy exchange is estimated between groups of scene elements (e.g., polygons) but between individual elements.

• To solve more simple problems in each part.

• To process each part in parallel (coarse grain parallelism).

• To process each part by different modifications of radiosity algorithms.

• Possibility of recursive subdivision to gain a desirable rate of computational cost.

• Possibility to assemble a scene from precomputed parts.

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, L_e(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, N_x,N_y- 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

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.

  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.

  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 

Denote

Thus the system has a view: 

. 

The approximating solution is given by the following iteration scheme: 

, 

and

  , . 

Terms and

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

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 F₁₁ for the scene part D₁
2. Computing of the table F₂₂ for the scene part D₂.
3. Computing of the table F₁₂, which is used while energy transfer from D₂ to D₁.
4. Computing of the table F₂₁, which is used while energy transfer from D₁ to D₂.
5. Calculation of equilibrium of light energy in scene part D₁.
6. Calculation of equilibrium of light energy in scene part D₂.
7. Fulfill energy transfer from the scene part D₁ to D₂, using matrix  F₂₁.
8. Fulfill energy transfer from the scene part D₂ to D₁ using matrix F₁₂.
9. Single iteration of the algorithm consists of steps 5 � 8.

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.

  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"

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

Scene \ Steps	5	6	7	8
Constr.wood	16	16	0..5	0.5	33	17
Shadow	40	56	2	1	99	59

  3. Conclusion

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

• Formal subdivision of a scene into parts without preliminary analysis of its geometry that could be quite time-consuming process.

• To solve more simple problems in each scene part.

• To process each part in parallel.

• To process each part by different modifications of radiosity algorithms.

• Possibility of recursive subdivision of a scene to gain a desirable rate of computational cost.

• Possibility to assemble a scene from precomputed parts. Really it means that we need to compute only sub-matrices and

  Acknowledgements

References

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 http://www.cs.cmu.edu/afs/cs/user/ph/www/femrad93.ps.gz.
   
   
3. Hanrahan, P., Salzman, D., Aupperle, L. A rapid hierarchical radiosity algorithm. Computer Graphics 25, 4 (July 1991), p. 197-206.
   
   
1. Appel, A.A. An efficient program for many-body simulation. SIAM J. Sci. Stat. Computing 6(1), 1985, p.85-103.
   
   
2. Willmott, A.J., Heckbert, P.S., Garland, M. Face Cluster Radiosity. Eurographics Workshop on Rendering, Granada, Spain, June 1999.
   
   
3. Willmott, A.J., Heckbert, P.S. An Empirical Comparison of Radiosity Algorithms. Technical Report CMU-CS-97-115. See also http://www.cs.cmu.edu/~radiosity/emprad-tr.html
4. Willmot, A.J. Hierarchical Radiosity with Multiresolution Meshes. Ph.D. thesis proposal. 206 p. Draft 17/6/2000. http://sulfuric.graphics.cs.cmu.edu/~ajw/thesis/ajw-thesis.pdf
5. http://www.lightscape.com/main/
   
   
6. 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