Random walk radiosity with infinite path length1

PII:

Comput. & Graphics, Vol. 22, No. 2±3, pp. 161±166, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0097-8493(98)00002-8 0097-8493/98 $19.00 + 0.00

WSCG'97

RANDOM WALK RADIOSITY WITH INFINITE PATH LENGTH{ MATEU SBERT Institut d'InformaÁtica i Aplicacions, Universitat de Girona, Lluis Santalo s/n, E 17003 Girona, Spain e-mail: [email protected] AbstractÐThe error in an unbiased Monte Carlo method is characterized by the variance. By knowing the variance of di€erent Monte Carlo estimators for Radiosity (and also their cost) we should be able to obtain the most ecient of them. This paper gives the variances for two such estimators, the shooting and gathering in®nite path length random walk estimators. This completes a previous work of the author on ®nite path length estimators. # 1998 Elsevier Science Ltd. All rights reserved. Key words: rendering, radiosity, Monte Carlo, random walk.

1. INTRODUCTION

We study in this paper the variance of two Monte Carlo discrete random walk estimators for radiosity, shooting and gathering random walk with in®nite path length. The study of the variance is important because it gives the expected square error for an unbiased estimator. Thus, by knowing the variance of di€erent Monte Carlo estimators for radiosity (and also their cost) we should be able to obtain the most ecient. Gathering random walk proceeds sending paths from the patches of interest to gather energy when a source is hit. Shooting random walk shoots paths carrying energy from the sources, to update the visited patches. We consider in both cases that the random walk proceeds according to the discrete form factor matrix (that is, patch-to-patch form factors). In®nite path length estimators are such that the path never ends. Obviously, in a simulation, one has to cut o€ the path, obtaining a biased estimator. In this case the expected square error is the variance plus the squared bias. However, if we can assure a small bias, the variance will give a good approximation of the square error. The solution obtained with a random walk estimator will converge to the solution of the Radiosity system (see for instance [1, 2] for random walk solutions of a system of equations). Also, when the size of the patches decreases, it will converge to the one found with Particle Tracing [3], which uses point-to-point form factors. Path-tracing [4], and even distributed ray-tracing [5, 6] can be considered as the limiting case of gathering random walk for the non-discrete case (without the shadow ray). Bidirectional ray-tracing [7, 8] is a mixture of { Extended version of the paper ``Variance of two in®nite path length random walk radiosity estimators'' presented at WSCG'97, Plzen, February 1997. 161

non-discrete shooting and gathering. Refs [9] and [10] can be seen as a breadth-®rst approach to a shooting random walk estimator, which in turn would be the depth-®rst approach. As far as the author knows, only two papers have addressed to date the variance of random walk estimators for radiosity. Ref. [11] gives a bound for a shooting estimator, and Ref. [12] gives variances for a wide family of shooting and gathering estimators. The purpose of this paper is to enlarge the work in Ref. [12] addressing the in®nite path length estimators. The structure of this paper is the following. In Section 2 we will give previous results on ®nite path length estimators. In Sections 3 and 4 the variance of the shooting and gathering in®nite path length estimators will be obtained. In Section 5 results will be given to support the theoretical ®ndings, and our conclusions and future research will be presented in Section 6. 2. PREVIOUS WORK

In Ref. [11] a bound was given for a shooting ®nite path length estimator. This bound was intended more to study the complexity than to give a realistic approximation of error. In Ref. [12] di€erent shooting and gathering estimators with ®nite expected path length were studied, and their variances given. Three estimators for shooting, FT/ 1 ÿ Ri, FT/Ri and FT and three for gathering Es/ 1 ÿ Rs, Es/Rs and Es were analyzed. In Table 1 the characteristics of the di€erent estimators are given. The shooting estimators, ®rst de®ned in Ref. [11], are as follows: A path (or particle) starts from source s with probability proportional to the power of the source. It advances from patch to patch according to the

162

M. Sbert Table 1. Di€erent random walk estimators (the meaning of the di€erent quantities is in Table 2) Patch scored Shooting FT 1ÿRi FT Ri

Variance  ÿ bi     bi FAT Ri i R1i ‡ 2xi ÿ bi   bi RAi Fi T …1 ‡ 2Ri xi † ÿ bi

Last

bi

All but last

FT Gathering

All

Es 1ÿRs



FT Ri Ai …1ÿRi †

Last

Ri pi

Es Ss …1ÿR bis ÿ b2i s†

Es Rs

All but last

Ri pi

2 s Ss Es ‡2b Rs bis ÿ bi

Es

All

form factor probabilities (simulated with a cosinus distribution [9]). On each hit patch i, a die-survive (or absorption-re¯ection) test is done according to the probabilities (1 ÿ Ri, Ri), where Ri is the re¯ectivity of the patch. The patch radiosity gets updated ever (FT estimator), only when the path survives (FT/Ri estimator) or only when the path dies (FT/ 1 ÿ Ri estimator). The name of the estimators re¯ect the amount by which the received power (and hence the radiosity) of the patch increases. This is for a single path simulation, with N particles or paths each one carries the total power divided by N. As for the gathering estimators, a path begins at a patch i with probability pi, and the same transition and die-survive probabilities as for the shooting case are used. On each hit source the radiosity of the patch i origin of the path gets updated ever (Es estimator), only when the path survives (Es/Rs estimator) or only when the path dies (Es/1 ÿ Rs estimator). Again the name of the estimators re¯ect the amount by which the received radiosity of the patch increases (after dividing by pi). With N particles or paths it should be divided by N. The best estimators for each case (that is, the ones with minimum variance) were found to be the FT and Es, which agrees with intuition, as both update each patch in a path (this advantage can be shown to overbalance the positive covariances). Two in®nite path length estimators were also characterized in Ref. [12], but we were only able to give bounds for their variances. These two estimators correspond in practice to biased estimators, because the path is cut o€ under some criterion (having reached a predetermined length, or being the accumulated re¯ectivities inferior to some threshold). In the next two sections we will reexamine those estimators and ®nd their variances.

3. AN INFINITE PATH LENGTH SHOOTING ESTIMATOR

Let us ®rst consider what the expected value of any unbiased Monte Carlo estimator should be for the incoming power on a patch. Let us suppose that

Ri pi

Ss …Es ‡ 2bs †bis ÿ b2i

the initial power of source s is Fs, fi is the incoming power on patch i, Fkl denotes the form factor from patch k to patch l, and Rk denotes the re¯ectance of patch k. Then we have, by developing the power system in Neumann series and dropping the zero order term: X XX fi ˆ Fs Fsi ‡ Fs Fsh Rh Fhi s

‡

s

XXX s

h

h

Fs Fsh Rh Fhj Rj Fji ‡




j

This can be expressed as: …2† …3† fi ˆ f…1† i ‡ fi ‡ fi ‡




where f(1) f(2) i =SsFsFsi, i =SsShFsFshRhFhi, (3) fi =SsShSjFsFshRhFhjRjFji and so on. That is, represents the power arrived directly from the f(1) i sources, f(2) represents the power arrived after one i bounce, as so on. Let us now consider the following simulation. An in®nite length path g starts from source s with probability ps=Fs/FT (that is, according to the power of the source), and from here on it evolves according to the transition probabilities given by the form factors. For instance, from s it will go to patch j with probability Fsj. The random variables (2) (3) f(1) i , fi , fi , . . . are de®ned in the following way: All of those random variables are initially null. If the path g arrives at patch i at length l, and if s, h1, h2, . . ., hl ÿ 1, i is the trajectory the path has folTable 2. Meaning of the di€erent quantities appearing in Table 1; the sux i means for patch i, sux s indexes the sources Ei bi FT Ai Ri xi bis pi

Emissivity Re¯ected radiosity = BiÿEi Total power in scene Area Re¯ectivity Received power (or radiosity) due to self-emitted unit power (or emittance) Re¯ected radiosity on i due to source s Probability for a path to begin at i

Random walk radiosity with in®nite path length

lowed, then the value of f(l) i is set to Rh1 Rh2 . . . Rhl ÿ 1 FT. Let us de®ne also a new random variable fi as:

X X s

h

Rh FT

Fs Fsh Fhi ˆ f…2† i FT

and so on. Then, we have ^ …2† ^ …1† ^ …2† E…f^ i † ˆ E…f^ …1† i ‡ fi ‡ . . .† ˆ E…fi † ‡ E…fi † ‡ . . . …2† ˆ f…1† i ‡ fi ‡ . . . ˆ fi

Thus the random variable f(1) is a centered estimai tor for the power arrived to patch i after l bounces, and the sum of all this family of estimators gives a new centered estimator fi which corresponds to the total incoming power arrived to patch i after any number of bounces. Our aim now is to obtain the variance for this estimator. We can decompose Var(fi) in the following way ^ …2† Var…f^ i † ˆ Var…f^ …1† i ‡ fi ‡ . . .†

…1†

1n
We ®nd ®rst the cross terms, which are not null because a path can arrive at length n on patch i and again at length m: XX XX X ^ …m† E…f^ …n† ... ... Rh1 . . . Rhnÿ1 FT
i fi † ˆ s

h1

hnÿ1 hn‡1

hmÿ1

Rh1 . . . Rhnÿ1 Ri Rhn‡1 . . . Rhmÿ1 FT
Fs Fsh . . . Fhnÿ1 i Fihn‡1 . . . Fhmÿ1 i FT 1 XX X ˆ ... R2h1 . . . R2hnÿ1 s

h1

hnÿ1

Fs Fsh1 . . . Fhnÿ1 i
X X Ri FT ... Rhn‡1 . . . Rhmÿ1 . . . hn‡1

i

i

1n

i

n
i

where ci is the incoming energy in the same environment having changed all the re¯ectivities by their square and xi is the expected value of the incoming power (or radiosity) on patch i due to a unit power (or unit emittance) on the same patch i. xi is equal to Iiÿ1/Ri, where Ii is the self-importance of patch i (for a de®nition of importance, see [13]). Also we have X Fs …1† E…f^ i…1†2 † ˆ …FT †2 Fsi ˆ FT f…1† i ˆ FT ci F T s (1) because f(1) i =ci . XX Fs E…f^ i…2†2 † ˆ …Rh FT †2 Fsh Fhi ˆ FT c…2† i FT s h

and so on. Then we ®nally obtain Var…f^ i † ˆ FT …ci…1† ‡ c…2† i ‡ . . .† ‡ 2Ri FT ci xi ÿ f2i ˆ ci FT …1 ‡ 2Ri xi † ÿ f2i where ci is the incoming energy in the same environment having changed all the re¯ectivities by their square. For the radiosity estimator BÃi=Ei+Ri/Aifi we have

2 ^ …2† ^ 2 ˆ E……f^ …1† i ‡ fi ‡ . . .† † ÿ …E…fi ††

ˆ E…f^ i…1†2 † ‡ E…f^ i…2†2 † ‡ . . . X 2 ‡2 E…f^ i…n† f^ …m† i † ÿ fi

n
ˆ 2Ri FT ci xi

The expected values are X Fs E…f^ …1† FT Fsi ˆ f…1† i †ˆ i FT s E…f^ …2† i †ˆ

emittance) on the same patch i. Then ! X …n† X …mÿn† X X E…f^ …n† f^ …m† † ˆ 2Ri FT c x 2 1n

^ …2† ^ …3† f^ i ˆ f^ …1† i ‡ fi ‡ fi ‡ . . .

163

hmÿ1

Fihn‡1 . . . Fhmÿ1 i
…mÿn† ˆc…n† i Ri FT xi

where c(n) is the incoming energy after n bounces in i the same environment having changed all the re¯ecÿ n) tivities by their square and x(m is the expected i value of the incoming power (or radiosity) on patch i after m ÿ n bounces due to a unit power (or unit

Var…B^ i † ˆ ˆ

R2i …c FT …1 ‡ 2Ri xi † ÿ f2i † A2i i b i FT …1 ‡ 2Ri xi † ÿ b2i Ai

…2†

where bi is the re¯ected radiosity (or radiosity due to incoming energy) in the same environment having changed all re¯ectivities by their squares. The variance of this estimator can be shown to be lower than the one for the FT estimator [12]. But this does not necessarily mean that this estimator is better. First, the cost of the in®nite estimator is determined by how we cut o€ the in®nite path. This can be done on a predetermined length, or alternatively when the accumulated re¯ectivity (that is, the product of re¯ectivities along a path) is less than a preestablished threshold. Second, this procedure imposes a bias on the solution which should be taken into account when comparing errors. If we want no bias, we can use Russian Roulette, which consists simply in switching to some ®nite path length estimator (such as the ones in Table 1) to distribute the left energy. This will however increase the cost considerably. We can alternatively consider acceptable this small percentage of undistributed energy (the variance accounts for the noise in the image, the bias for the undistributed energy). Also, if this threshold is small enough we can be con®-

164

M. Sbert

dent that the variance of the resulting biased estimator is close to the variance of the in®nite length estimator. 4. AN INFINITE PATH LENGTH GATHERING ESTIMATOR

Let us ®rst consider what the expected value of any unbiased Monte Carlo estimator should be for the radiosity of a patch. Let us suppose that the emittance of source s is Es, bi is the re¯ected radiosity, or radiosity of patch i due to the received power (that is, bi=BiÿEi, and so for a non-emitter patch, it equals the total radiosity), Fkl denotes the form factor from patch k to patch l, and Rk denotes the re¯ectance of patch k. Then we have, by developing the radiosity system in Neumann series (dropping the zero order term): X XX bi ˆRi Es Fis ‡ Ri Es Fih Rh Fhs s

‡ Ri

XXX h

j

s

h

s

Es Fih Rh Fhj Rj Fjs ‡




This can be expressed as: …2† …3† bi ˆ b…1† i ‡ bi ‡ bi ‡




where b(1) b(2) i +RiSsEsFis, i =RiSsShEsFihRhFhs, b(3) =R S S S E F R F R F i s h j s ih h hj j js and so on. That is, i b(1) represents the radiosity due to direct illuminai represents the radiosity after one tion, b(2) i bounce, and so on. It is also useful to de®ne the following quantities: …2† bis ˆ b…1† is ‡ bis ‡




b(1) is represents the radiosity due to direct illumination from source s, b(2) represents the radiosity is after one bounce from source s, and so on. It is clear that: X bi ˆ bis

Now let us ®nd the expected value of those random variables. Applying the de®nition of expected value, and remembering that the probability of selecting patch i is pi, the probability of landing on source s just after leaving patch i is Fis, we have E…b^…1† i †ˆ

^…2† ^…3† b^i ˆ b^…1† i ‡ bi ‡ bi ‡




s

Ri

Es pi Fis ˆ b…1† i pi

Now, to go from patch i to a source s in a two length path we can pass through any patch h, so we have XX Es E…b^…2† Ri Rh pi Fih Fhs ˆ b…2† i † ˆ i pi s h and so on. Then, we have ^…2† ^…1† ^…2† E…b^i † ˆ E…b^…1† i ‡ bi ‡


† ˆ E…bi † ‡ E…bi † ‡


…2† ˆ b…1† i ‡ bi ‡


ˆ bi

So it is clear that the random variable bÃ(l) i is a centered estimator for the radiosity due to the power arrived on patch i after l bounces, and the sum of all this family of estimators gives a new centered estimator bÃi which corresponds to the total radiosity of patch i due to the power arrived after any number of bounces. Our aim now is, as before, to obtain the variance for this estimator, which we will do decomposing Var(bÃi) in the same way as in Ã(m) Equation (1). The terms of the form E(bÃ(n) i bi ) are not null, because if a path arrives at length n on source s it can also arrive later at source s' at length m, and we ®nd them as in the previous section: E…b^i…n† b^i…m† † ˆ

s

Let us now consider the following simulation. A in®nite length path g starts from patch i with probability pi (this probability can be considered as the initial or emitted importance of the patch), and from here on it evolves according to the transition probabilities given by the form factors. For instance, from i it will go to patch j with probability Fij. Let us de®ne now the ranÃ(2) Ã(3) dom variables bÃ(1) i , bi , bi , . . . in the following way: All of those random variables are initially null. If the path g happens to arrive at source s at length l, and if i, h1, h2, . . ., hl ÿ 1, s is the trajectory the path has followed, then the value of bÃ(l) is set i to RiRh1Rh2 . . . Rhl ÿ 1 Es/pi. Let us de®ne also a new random variable bÃi as:

X

ˆ

XXX s

s0

...

h1

XX

...

hnÿ1 hn‡1

X

Ri Rh1 . . .

hmÿ1

Rhnÿ1

Es
Ri Rh1 . . . Rhnÿ1 Rs Rhn‡1 . . . pi

Rhmÿ1

Es0
pi Fih1 . . . Fhnÿ1s Fshn‡1 . . . Fhmÿ1s0 pi

XX s

h1

Fhnÿ1s


...

X hnÿ1

R2i R2h1 . . . R2hnÿ1

XXX s0

Es pi Fih1 . . . pi

Rs Rhn‡1 . . .

hn‡1 hmÿ1

Rhmÿ1 Fshn‡1 . . . Fhmÿ1s0

Es0 pi

1 X …n† X …mÿn† 1 X …n† …mÿn† b b 0 ˆ b b pi s is s0 ss pi s is s where b(n) is the radiosity due to the incoming is energy after n bounces in the same environment having changed all the re¯ectivities by their square.

Random walk radiosity with in®nite path length

165

and so on. Then we obtain 1X 1X Var…b^i † ˆ Es …bis…1† ‡ b…2† b bs ÿ b2i is ‡


† ‡ 2 pi s pi s is ˆ

1X …Es ‡ 2bs †bis ÿ b2i pi s

where bis is the radiosity due to the incoming energy in the same environment having changed all the re¯ectivities by their square. For the radiosity our estimator is simply bÃi+Ei, and as Ei is a constant we have 1X Var…b^i ‡ Ei † ˆ …Es ‡ 2bs †bis ÿ b2i …3† pi s

Fig. 1. Numbering the patches in the test scene. Patches 1±9, with re¯ectance 0.3, are not shown. Patch 4 is the emitter. Patches 10±18 have re¯ectance 0.4, 19 to 27 0.5 and so on.

This variance can be proven to be less than the one for the Es estimator [12]. But again here we can repeat the considerations of the previous section on which estimator is better. Both shooting and gathering variances found are for a single path. With N independent paths, the variances must be divided by N. 5. RESULTS

Then X X n
1n

ˆ

1X pi

s

! ^…m† E…b^…n† i bi †

ˆ

1 X X …n† X …mÿn† b b pi s 1n is n
bis bs

and also E…b^i…1†2 † ˆ

X s

Ri

Es pi

2

pi Fis ˆ

X Es s

pi

b…1† is

 X X X Es …2† Es 2 …2†2 ^ E…bi † ˆ Ri Rh pi Fih Fhs ˆ b pi pi is s s h

Presented in Fig. 2 are some experiments performed on a very simple scene, a cubical enclosure with each face divided into nine equal size patches (Fig. 1), the re¯ectivities of the faces being 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 respectively, and a source with emissivity 1 in the middle of the ®rst face, in patch 4. Thus patches 1±9 receive no direct lighting and have re¯ectivity 0.3, patches 10±18 re¯ectivity 0.4, and so on. For this scene we computed a reference solution with a ®nite path length estimator and 108 paths. This provided us with the bi values. A reference solution was also computed for the same scene with the re¯ectivities squared, that is, each re¯ectivity was substituted by its square. This provided us with the bis values. Then 100 runs of 104 paths each for both in®nite estimators were computed (taking for gathering pi=Ai/AT, the fraction of total area), and used to obtain the square errors, and thus an

Fig. 2. Comparison of the expected variances (Q) and the experimentally obtained square errors for the (a) shooting and (b) gathering in®nite path length estimator, for the 54 equal area patches of a cube (on x axis), with face re¯ectivities 0.3, 0.4, 0.5, 0.6, 0.7 and 0.8. A source with emittance 1 is in the middle of the ®rst face.

166

M. Sbert

estimated value of the variances for a single path. The in®nite paths were cut o€ when the product of re¯ectivities along the path was <0.001. The formulae for the variances are the formulae 2 and 3, with the approximation xi=0. Figure 2 shows that the obtained results are in concordance with the theoretically expected ones. Although the scene used in the test has no occlusions, it should be noted that the variance of a patch radiosity does not depend on whether it is due to direct or indirect illumination. 6. CONCLUSIONS

We have given here the variances for the in®nite path length random walk estimators for radiosity. These formulae complete the study in Ref. [12] for ®nite path length estimators. A study of the relative eciencies of both kind of estimators remains to be done, and also the generalization of the results to the R-G-B case. AcknowledgementsÐThis project has been funded in part with Grant No. TIC 95-0614-C03-03 of the Spanish Government. Many thanks to my supervisor, Xavier Pueyo, and to Joaquim GelabertoÂ. REFERENCES

1. Hammersley, J. M. and Handscomb, D. C., Monte Carlo Methods, Methuen & Co, London, 1975. 2. Rubinstein, R. Y., Simulation and the Monte Carlo Method. John Wiley & Sons, New York, 1981.

3. Pattanaik, S. N. and Mudur, S. P., Computation of global illumination by Monte Carlo simulation of the particle model of light. In Proceedings of the Third Eurographics Workshop on Rendering, 71. Bristol, U.K., May 1992. 4. Kajiya, J. T., The rendering equation. Computer Graphics (SIGGRAPH'86 Proc.), 1986, 20, 143±150. 5. Cook, R. L., Porter, T. and Carpenter, L., Distributed ray tracing. Computer Graphics (SIGGRAPH'84 Proc.), 1984, 18, 137±145. 6. Ward, G. J., Rubinstein, F. M. and Clear, R. D., A ray tracing solution for di€use interre¯ection. Computer Graphics (SIGGRAPH'88 Proc.), 1988, 22, 85±92. 7. Veach, E. and Guibas, L., Optimally combining sampling rendering techniques for Monte Carlo rendering. Computer Graphics (SIGGRAPH'95 Proc.), 1995, 29, 419±428. 8. Lafortune, E. and Willems, Y., Bi-directional path tracing. In Proceedings of CompuGraphics, 145. Alvor, Portugal, December 1993. 9. Shirley, P., A ray tracing method for illumination calculation in di€use specular scenes. In Proceedings of Graphics Interface '90, 205. Toronto, Ontario, May 1990. 10. Feda, M. and Purgathofer, W., Progressive ray re®nement for Monte Carlo radiosity. In Proceedings of the Fourth Eurographics Workshop on Rendering. Paris, June 1993, pp. 15±26. 11. Shirley, P., Time complexity of Monte Carlo radiosity. In Proceedings of EUROGRAPHICS'91. Vienna, September 1991, pp. 459±465. 12. Sbert, M., Error and complexity of random walk Monte Carlo radiosity. IEEE Transactions on Visualization and Computer Graphics, 1997, 3, 23±38. 13. Pattanaik, S. N. and Mudur, S. P., The potential equation and importance in illumination computation. Computer Graphics Forum, 1993, 12, 131±136.