Accelerated radiosity method for complex environments

Comput & Graphics Vo[. IJ,. No. I. pp. 6 5 - ' 1 . 1990 Printed in Gr~at 8drain.

0097-8493/90 $ 3 0 0 ~" .00 ~ 1990 Pergamon Press plc

Eurographics '89 Award Paper

ACCELERATED RADIOSITY METHOD FOR COMPLEX ENVIRONMENTS HAU XU, QUN-SHENG PENG and YOU-DONG LIANG National Laboratory of CAD & CG, Zhejiang University. Hangzhou, China Abstract--As form-factor calculation costs about 90% of the computing time when applying the radiosity approach for realistic image synthesis, it is of great significanceto reduce the required computation. An accelerated radiosity algorithm for general complex environments, based on environment localization and the directional form-factorconcept, is presented in this paper. First, we subdivide the object space into many regions. Objects contained in each region are adjacent to each other and pose more global illumination effects to their neighbours. Next, form-factors are calculated in each local environment. The radiant light energy transfer between different regions is evaluated at their common boundaries. Directional form-factors are introduced to simulate the interaction of light between local environments and between nondiffuse surfaces. Comparisons are made to existing algorithms. Statistic results and theoretical analysis show that the new algorithm is much faster than previous ones. The technique is especially fit for interactive design and animation sequence generation since modification to the shape or location of objects usually happens in local environments. because of its overwhelming computational expense. Wallace's approach seems more clever by calculating the global illumination effects of the environment in two passes. The first pass accounts for the diffuse component of the light energy leaving each surface with extension to include the effects of specular to diffuse reflection and transmission between surfaces, while the second pass deals with the specular component and generates the desired image. A so-called "image" method is adopted in its implementation to simulate the planar mirror surfaces which exert light energy transfer via specular reflections. The performance of this method may degenerate when multiple mirror reflections are inspired and the computational expense may possibly go up by exponential order. This paper presents an accelerated radiosity approach for general complex environments based on environment localization and the directional formfactor concept. Computational expense of the new approach for diffuse environments is near linear instead of square, while additional expense for nondiffuse environment is trivial. Furthermore, our approach has great application potential for CAD and animation. Interactive design is an important feature of CAD. Traditionally, once a modification is made, the whole display process will have to be repeated no matter how trivial the modification is. The radiosity method, due to its view-independence, lends itself to interactive design. Views of a scene in different directions may be readily generated without much computation. Unfortunately, any geometric modification in the environment necessitates the recalculation of all form-factors, which is expensive. Environment localization may promote the interactivity of a design cycle, and quick display can be expected since form-factor calculation is limited within a few local environments in which modifications are made. Animation is another potential application area of our technique. Animation and realistic image synthesis are two main branches of computer graphics. Up to

I. I N T R O D U C T I O N

The radiosity method is one of the major rendering techniques for realisticimage synthesis.Since Goral et al.[l] introduced the method from thermal engineering to computer graphics in 1984, much research work has been done on it:Cohen et al. [2 ] presented a hemicube algorithm in 1985 and a substructuring technique in 1986; Immel et al.[4] extended the radiosity approach to nondiffuse environment in 1986; Wallace et M.[5] proposed a two-pass solution in 1987 which integrates the radiosity method and ray-tracing for e@cient rendering of general environments. One of the key issues in applying the radiosity method for rendering is the form-factor calculation. The form-factor of patch i to patchj specifiesthe fraction of the energy radiated from patch i which lands on patch j. Obviously, the form-factor calculation takes the complexity of O(NZ), N representing the total number of patches in the environment. Statisticsshow that it costs more than 90% of the computational expense of the radiosity method. In order to achieve a fast radiosity solution, approaches must be found to reduce the computation required by form-factor calculations. Although the substructuring technique may generate images of high quality without solving the whole diffuse environment on a fine subdivision basis, it is still excessive when the environment gets more and more complex. Light energy radiation in a nondiffuse environment is spatial and spectral dependent. It is recognized that lmmel's method for solving the light energy distribution among the whole environment by spatial and spectral sampling over each surface patch is impractical

This paper is one of the award winning papers from EUROGRAPHICS "89. the annual Conference of the European Association for Computer Graphics. The paper is published here in a revised form. with permission of North Holland Publishing Company (Amsterdam). the Publisher of EUROGRAPHICS "89 Conference Proceedings. 65

HAU XU, QuN-SHENG PENG. and YOU-DONG LIANG

66

date, they have been developing almost independently. Animation requires quick response while realistic image synthesis is usually time-consuming. Our technique is suitable for animation too, because form-factor calculation needs to be performed only for the local environments containing moving objects. 2. ENVIRONMENT LOCALIZATION

Previous radiosity approaches require the whole environment to be projected onto the hemi=cube of each surface patch. Hidden surfaces are then identified and removed using a Z-buffer algorithm, thus introducing computation expense of O ( N 2). As it is well known, light energy falls off as the inverse square of the distance the light travels. Thus, an object poses more global illumination effects to its neighbouring objects than to those located far away. With this observation, we can divide the whole space into some subspaces. Objects contained in each subspace are located near to each other. These objects along with the boundary of each subspace form a local environment. Form-factors are calculated in each local environment and then integrated to obtain a global solution. Light energy transfer between different regions is accounted for at their common boundaries. The boundaries of each subspace are assumed to be transparent, causing neither deflection nor reduction when light crosses them. Assume that an environment is divided into n local environments and each consists of p object patches and q boundary patches. Then, form-factor calculation for all patches within a local environment takes time of O((p + q)Z) and O ( n , ( p + q)2) for the whole en= vironment. If we keep p, q invariant by increasing n when the environment to be rendered turns more and more complex, then our algorithm appears to be a lin= ear approach. 3. THEORETICAL F U N D A M E N T A L S

0 = angle between surface normal N and c/i.; p = the bidirectional reflectance of the surface. Bidirectional reflectance / transmittance p(c/o.,, c/i.) represents the physical reflection properties of the surface (see Fig. I ). It describes the relationship between an incoming intensity I~. (c/i.) and the outgoing intensity lo.t (c/o.t) and can be expressed by four portions corresponding to: diffuse reflection, nondiffuse reflection, diffuse transmission and nondiffuse transmission[5]: P(

c/OUt, c/in )

kapa .~ k,p,(c/out, c/i.) =

~

(reflection),

bs

k,apta + k,,p,,( c/ou,, c/i.) (transmission). a"

bt

Solving Eq. ( 1 ) analytically is difficult. Numerical techniques, however, provide efficient means to approach a satisfactory solution. We approximate the light energy leaving point O on the surface with a number "of discrete directional intensities, each describing the average intensity of patch i along an outgoing direction d within a small solid angle ~2a. Thus:

I a= Ed + ~ f~ p(c/d,c/I)Iin(c/in)COSOdo.) l=I 1 = E a + ~ P(C/a, C/t) I- I

Ii.(¢i.)cos 0 d¢~.

(2)

t

Note that in a local environment Ii. (c/i.) may arrive at the surface from other surfaces within the local environment or from the boundary of the subspace, c/i. is just opposite to an outgoing direction of patch i. Thus,

AND FORMULATIONS

li.(c/i.) = ~ [j(C/i.)HID(O,j, c/i.)

3.1. Radiosity equations for global environments Kajiya [ 8 ] and lmmel [4 ] presented independently a rendering equation which gives a full description of global illumination effects:

Je

= ~ I y H I D ( O , j , c/i.)

jEP

where P is the set of all patches in the local environment and /7 t represents the average outgoing intensity of

&u,(c/oo,) = E(C/o°0

+ f P(C/.... c/i.)li.(c/i.)cos # dw

iEP

(l)

where

ou~,,~

~in

lou, = the outgoing intensity for the surface; E = the outgoing intensity due to emission by the surface: 9. = the sphere of incoming directions; c/out = the outgoing direction; c/i. = the incoming direction: I~. = an incoming intensity arriving at the surface from the environment:

Fig. I. Rendering equation geometry.

Accelerated radiosity method patch j within a small solid angle 12,. To obtain the average directional intensity I f of patch i, we apply an area integration to both sides of Eq. (2), i.e.,

67

thus

Ij' ' F~j/F,'

If=Z

( 5)

iEP

re = ~, -7-1 f a E:

7_, Yo

+

. p ( % , %) ~1 -7

jEP I'I

:

f~ Ja~ FtJ C°Src O J

× H I D ( i , j , ~i.)do~ dAi

= Ef + Z Zo ~rp(~,, %)Ii' -~,1 T fA~c°s0~r jEP I-I

i

where i is a boundary patch. In each local environment we have the following radiosity equations (for simplicity we will often refer to the radiosity equations for a local environment as "local equations"):

l

× H I D ( i , j , ~i,)dw dA~.

(3)

D

I f = E l + Z Z p(gz,,~k,)lf~rF~j

(6.1)

jEP I-!

Denote where i is an object patch; 1 f~ fa cos0

HID(i,j,~bl.)d,.,dAo

=

l,"

Z

IjFij/Fi '' '

(6.2)

jeP

Comparing F[~ with the standard form-factor expression, we find that they take almost the same form except their different integration areas. We define F[j as directional form-factor for it describes the fraction of energy leaving A~ which arrives on Aj in the discrete outgoing direction f~t. As ff=. is just opposite to the outgoing direction ~bt of patch i, thus, Eq. (3) can be rewritten as: D

/f = Ef + Z Z G % ( % , ¢,,)~-r~.

(4)

j E P I=l

Note that the bidirectional reflectance/transmittance of a given surface can be determined according to its actual reflection attributes. However, we still have to specify p(~bd, ~b,) for the hypothesis boundary patch. Imagine the boundaries of each subspace as perfect transparent surfaces, which do not affect the light propagation passing through. Thus, the outgoing intensity of a boundary patch in a specified direction should be the same as its corresponding incoming intensity both in magnitude and in direction: zoo,, = t ; J ,

hence according to (2)

o

loa= = Z #(~ba, fft)Ii.(¢'i.)

f,

cos 0 dco,

l- I

,'

where i is a boundary patch. Note that a single set of local equations is not of full rank and cannot be solved independently. To find the light energy distribution among the entire environment, we must solve all local equations simultaneously. Unfortunately, the number of unknowns contained in such equations may be enormous and solving them directly appears to be an unpractical approach. By close observation, we find that the directional intensities of the boundary patches of each subspace are just intermediate variables and can be eliminated from the whole equations. This is accomplished by repeatedly replacing any directional intensities of the boundary patches of subspace appearing in (6.1) with (6.2) until no one appears. Specifically, if the light.leaving patch i along direction 1 passes through a boundary patch k and reaches patch j, then

=

and F~j is accumulated within the solid angle f~. The procedure is implemented simply by accumulating F~ with combined form-factor F[kF~j/F~where k is a boundary patch. Memory can be saved with the fact that F[k, F~j will never be referred to after the combination. This is significant for the environment with many local environments. After the replacement, we have obtained a set of global mdiosity equations: N

we get I p(¢,,, %) =

0

% ÷ - ~#

l

% = -~,.

Fi~Fkj/Fk

D

/~, = E¢ + Z Z p ( % . ¢,,)/;'.F~,, j=l

1-1

i = I, 2 . . . . .

N.

(7)

f cos 0 d~ I

Denote:

The equations are reduced to the conventional form for a diffuse environment: N

F~

= f~ cos 0 d,.,/r, I

14:l-F

I;=E;+p, ~, lifo, j=|

i = 1,2 . . . . . N

(8)

68

HAU XU, QUN-SHENG PENG, a n d Y o u - D O N G LP, NG

where D

F,,= EFt,. /= I

3.2. Two-pass method for nondiffuse environment The key issue of the successful adoption of the twopass method to realistic image synthesis is the inclusion of the illumination effects caused by diffuse to specular, specular to specular, specular to diffuse light energy transfer mechanisms. This is achieved in its current implementation by employing an image method in the preprocess to calculate the additional form-factor which estimates the indirect light energy transfer between two diffuse surface patches. The image method is not an ideal approach because it turns the real environments into much more complicated ones and it is restricted to perfect planar mirrors, not suitable for transmission and fuzzy reflection. Here we present a general approach to evaluate the indirect illumination effects between two surface patches which not only costs fewer computations but also accounts for transmission and fuzzy reflection. It can be seen that directional intensity I/a in (7) consists of two components, diffuse and nondiffuse:

l/a=Di+N/a, E/a=Hi+L/a

kapa +

k,p,(g'o.t,g'i.)

a"

bs

Thus

Oi + N/a = E/a +

;;(:

,.)

k d4 b,

j=t /=l 1 l X (Dj + N~)rF o N

D

= Hi + Z Z (Dj + NI)kaoaF~J + L/a j=t I-t N

N

Di = H, + Pdka Z DjFij.

k = 1,2 . . . . .

N.

j-I

The replacement procedure is, in fact, a mathematical simulation of light propagation through the entire environment. During the process, form-factors between any two diffuse surface patches are accumulated and the accumulated part is nothing but Wallace's additional form-factors. The improved radiosity equations have therefore included the global illumination effects by specular to diffuse light transfer mechanism, which are not presented by the traditional radiosity equations. The approach can be easily extended to account for transmission and fuzzy reflection since there are no extra constraints on the bidirectional reflectance. 4. C A L C U L A T I O N O F T H E D I R E C T I O N A L FOILM-FACrORS

and so does the bidirectional reflectance:

Pi =

may form nondiffuse loops, i.e.. a portion of the light energy radiated from a nondiffuse surface patch may feed back to itself after several specular reflections or transmissions, in which case, the process of replacement may never stop. To avoid this, depth of the replacement must be controlled. After cutting off the nondiffuse remainder, we have obtained the following equations:

D

+ Z Z (Oj + N~)k,p,(¢d, Ct)~rF~j/b,. j=t I-t

In the previous section we introduced the concept of directional form-factor. Calculating these directional form-factors analyticallyisdil~cult.The hemi-cube algorithm offers a way to evaluate them. However, direction consistency must be kept between different hemi-cubes. As each hemi-cube has its own local coordinate system, the constraint m a y pose considerable computation. To facilitatethe identificationof directions, Immel et al. presented a global cube method[4]. In the following, we extend the global cube algorithm to calculate the directional form-factors. First,each side of a global cube isdivided into square cells.Each cellspecifiesa discretedirection.Since global cubes of allpatches have the same orientation,direction consistency within the whole environment is easily maintained. Each cell is further divided into pixels. Each pixel defines a small solid angle of view from the

Decomposeit: N

D,=Hi+

D

Z Z(Dj+N~)kapaF~j

(9.1)

j - I i=l

N

~

D

= L/a + N E (Dj + N)k,m(¢~, ¢t)=Fb/b,. (9.2) Analogous to the process of removing all I/a of boundary patches from the radiosity equations, we repeatedly replace Nia in (9.1) with (9.2). Note that the light energy transfer within a nondiffuse environment

Fig. 2. Directional form-factor geometry.

pixel

69

Accelerated radiosi~ method Table I. Statistics of environment complexity.

center of the cube. Surfaces within the local environment are projected onto the global cube and hidden surfaces are then removed at the pixel level (see Fig. 2). Directional form-factors may be calculated with a numerical integrating technique: F~j =

The number of local environments The number of object patches The number of boundary patches

fn HID(i,J, ~b)cos0 doJlTr I

=

~

2

84 217 129

78

3

4

5

152

163

169

102

156

108

HID(i, j, p)AF,

pEPIXEL(I)

where AFp is a delta

1

more accurate solution with less direction sampling. This is significant for complex environments.

form-factor:

and up and ~ are small solid angle and direction associated with each pixel. Unlike Immel's method, our approach is a two-pass solution. The directional form-factor concept introduced for correctly calculating the diffuse component of the light energy leaving each surface patch permits

5. IMPLEMENTATIONSAND DISCUSSIONS We implemented the algorithm on a Micro VAXII with a PS340 raster display. The flow chart in Fig. 3 gives an overall description of our implementation. We did not incorporate the substructuring technique into our algorithm, though its embodiment is straightforward. Fig. 4 demonstrates a moderately complex environment which contains 785 patches. The traditional ra-

and Meshing

1

I FormFactor Calculation

M

I

0

d i f i

I

Integrating~ [ofm-factom into the globalones

C a

Bidirectionalreflectance J

t i 0 n

Additionalform-factom l

&

A

~-.--.--~

Diffuse reflectance and emission I

n

i m a

Solution

! i 0 n

-V

l Viewing I

Rendering

I

Fig. 3. Program flow chart.

70

HAU XU, QuN-SHENG PENG, and You-DONG LIANG

Fig. 4.

Fig. 5.

Fig. 6.

71

Accelerated radiosity method Table 2. Comparison of two methods (1 = I min.). The nondiffuse environment

Form-factor calculation

The new approach

164.35

Traditional method

215.02

diosity method takes 215 minutes to calculate the formfactors for this environment while our approach takes 164 minutes to accomplish it. Table 1 and Table 2 list the statistics. Because of the limited complexity the improvement seems to be trivial. Fig. 5 demonstrates a nondiffuse environment. The environment is relatively simple and needs no localization. The computational cost for directional form-factors is 59 minutes while the replacement procedure for additional formfactor calculation takes only 61 seconds CPU time. Fig. 6 involves an environment almost the same as that in Fig. 5 except that the mirror surface is replaced by a diffuse one. Computing this environment takes 57.6 minutes. As we did not implement Wallace's image method, direct comparisons between our approach and Wallace's is not available. However, theoretical analysis predicts that Wallace's approach takes a bit less than twice of the cost, i.e., 115.2 minutes. Less computational cost is not the only benefit brought about by our approach; our method also fits for interactive design and animation applications. Note that form-factors of a static local environment remain invariant even though object movements or shape modifications are made in other local environments. Thus, only the form-factors of those local environments with varied substances need to be recalculated and then integrated for global solution. The main problem raised here is how to interpolate the patch's intensities as well as the trials between two or more keys. It calls for an elegant interpolation scheme. Developments of Computer Graphics in recent years have proven that the radiosity method is a powerful realistic image synthesis technique. However, there is still much work to be done to make it more practical. Speed is a major problem. On the other face of the computer graphics coin, interactivity and animation construct another research field. Nevertheless, their developments were almost independent from the efforts for pursuing realism; only a few works tried to relate to both fields[10]. We believe that the interpenetration of works in these two fields will greatly promote the

Integration for form-factors 20.12 0

Solution 3.5 3.02

development of Computer Graphics and this is also our future research direction.

Acknowledgments--We wish to thank Prof. T. G. Jin and members of the Computer Service Center at Zhejiang University for their support. Specialthanks go to Mr. P. Shao and Mr. M. Shao for stimulatingdiscussions.This research project is under the scheme of National Natural Science Foundation of China. REFERENCES i. C. M. Goral, K. E. Torrance, D. P. Greenberg, and B. Battaile, Modelingthe interaction of light between diffuse surfaces. Comp. Graphics 18(3), 213-222 (July 1984). 2. M. F. Cohen and D. P. Greenberg, The hemi-cube: A radiosity solution for complex environment. Comp. Graphics 19(3), 31-40 (July [985). 3. M. F. Cohen, D. P. Greenberg, and D. S. lmmel, An efficient radiosity approach for realistic image synthesis. Comp. Graphics andAppl. 6 (2), 26-35 (March 1986). 4. D. S. lmmel, M. F. Cohen, and D. P. Greenberg, A radiosity method for non-diffuse environments. Comp. Graphics 20(4), 133-142 (August 1986). 5. J. R. Wallace, M. F. Cohen, and D. P. Greenberg, A two pass solution to the rendering equation: A synthesis of my-tracingand radiosity method. Comp. Graphics21(4), 311-320 (July 1987). 6. A. Fujimoto, T. Tanaka, and K. Iwata, ARTS:Accelerated Ray-Tracing System. Comp. Graphics and Appl. 6 (3), 16-26 (April 1986). 7. Q. S. Peng, Y. N. Zhu, and Y. Liang, A fast ray-tracing algorithm using space indexing techniques. Proc. EUROGRAPH '87, Amsterdam, 11-23 (August 1987). 8. J. T. Kajiya, The rendering equation. Comp. Graphics 20(4), 143-150 (August 1986). 9. I. Navazo, D. Ayala, and P. Brunet, A geometric modeller based on the octree representation of polyhedra. Comp. Graphics Forum (May 1986). 10. A. N. Glassner, Spacetime ray-tracing for animation. Comp. Graphics andAppl. $ (3), 60-70 (March 1988). 11. E. M. Sparrow and R. D. Cess, Radiant Heat Transfer, Hemisphere Publishing Corp., Washington ( 1978). 12. R. L Cook, T. Porter, and L. Carpenter, Distributed ray tracing. Comp. Graphics 18(3), 137-145 (July 1984). 13. M. F. Cohen, S. Eric Chen, J. R. Wallace, and D. P. Greenberg, A Progressive refinement approach to fast radiosity image generation. Comp. Graphics 22(4), 75-84 (August 1988).