Computers & Graphics 26 (2002) 351–354
Technical Section
Optimal ray shooting in Monte Carlo radiosity Alex Brusia, Mateu Sbertb,*, Philippe Bekaertc, Robert Toblerd, Werner Purgathofere a
b
ECCSA, Barcelona, Spain " Departament d’Inform. i. Matematica, Institut d’Informatica i Aplicacions, Universitat Aplicada de Girona, Avgda. LI. Santalo, s/n, 17071 Girona, Spain c K.U. Leuven Belgium d VRVis Center Vienna, Austria e Vienna University of Technology, Vienna, Austria
Abstract In this paper we deal with the Monte Carlo radiosity problem of how to aim more rays towards where they are needed more: small patches and high-reflectivity regions. We derive the optimal ray shooting probabilities when taking the mean-square-error as the error measure. The standard case of using the form factors is shown to be optimal in the particular case of equal relation between reflectivity and square root of patch area throughout the scene. We also present an adaptive algorithm to compute optimal probabilities. Examples showing optimal direct illumination distribution are given, exhibiting concordance with the theoretically predicted results. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Radiosity; Monte Carlo; Direct illumination; Variance
1. Introduction Power has usually been expanded in Monte Carlo radiosity using the form factors probabilities. Either in breadth-first [9,1,4] or in depth-first (random-walk) Monte Carlo radiosity [6,2,8] the distribution of light has been done according to the cosine of the angle with the normal to the patch, a probability density distribution which mimics the form factors. In [3] a new strategy was introduced to send more rays towards complex regions of the scene, where the degree of complexity was considered mainly to be the number of objects found in the region. The hemisphere was divided into spherical triangles and weights based on the number of visible objects were obtained in a preprocess for each triangle. After that rays were distributed using these weights. This heuristic gave good results for small objects, while the overall appearance of the image did not suffer any *Corresponding author. E-mail addresses:
[email protected] (A. Brusi), mateu@ ima.udg.es (M. Sbert).
penalty. We will address a similar strategy here but from a different perspective. Instead of defining a heuristic, we will consider an error measure, the expected value of the mean square error (MSE), and find the probabilities that optimize it. We will thus show that the use of form factors is only optimal in a very particular case, when the relation of reflectivity to square root of patch area is equal for all patches. We will also give the optimal source selection probability for the case with more than one source. Finally, in the results section we apply the method for the computation of direct illumination for some test scenes and show the concordance with predicted results.
2. Mean square error and optimal distribution probabilities Without losing generality we restrict ourselves here to the expansion of direct illumination. Posterior shooting steps have identical treatment, although the technique
0097-8493/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 0 2 ) 0 0 0 6 1 - 4
A. Brusi et al. / Computers & Graphics 26 (2002) 351–354
352
will be most meaningful for important secondary sources. 2.1. One source case We consider the expected value of the MSE in direct illumination as the error measure to optimize. Given a diffuse source s with power Fs we send the light from it towards patch i with probability psi : The estimator for the reflected radiosity Bi is given by Ri Fs Fsi ; B# i ¼ psi Ai
ð1Þ
where Fsi is the form factor from source s to patch i and Ri and Ai are the reflectivity and area of patch i: It is an unbiased estimator, EðB# i Þ ¼
Ri Fs Fsi Ri Fs Fsi psi ¼ ¼ Bi : psi Ai Ai
ð7Þ
2.2. Optimal gain We want to obtain the quotient of the optimal MSE, MSEopt ; i.e. the one using probabilities (5), and the error obtained when using the standard probabilities, i.e. the form factors. Using formula (6) for the numerator and obtaining the denominator by substituting psi with Fsi in (4) we obtain after some algebra pffiffiffiffiffipffiffiffiffiffi P EðMSEopt Þ 2 ioj ðRi Rj = Ai Aj ÞFsi Fsj : ð8Þ ¼ P R2i R2j EðMSEÞ ioj ð Ai þ Aj ÞFsi Fsj Also, we can obtain the percent gain
ð2Þ
On the other hand, the variance of the estimator B# i is Ri Fs Fsi 2 Ri Fs Fsi 2 psi VarðB# i Þ ¼ psi Ai Ai R2i Fsi2 F2s 1 ¼ 1 : ð3Þ psi A2i Area averaging over all patches, we obtain the expected value of the MSE of the direct illumination: 1 X R2i Fsi2 F2s 1 1 ; ð4Þ EðMSEÞ ¼ AT i psi Ai where AT is the total area. This is the quantity to optimize. Using Lagrange multipliers to optimize (4) P with the constraint i psi 1 ¼ 0; we obtain as optimal values for psi Ri Fsi psi p pffiffiffiffiffi : Ai
where G is the normalization factor X Rk Fsk pffiffiffiffiffiffi : G¼ Ak k
ð5Þ
It is immediate to see that in the case when all areas and reflectivities are equal, the optimal probabilities correspond to the form factors. However, this is not a necessary condition to obtain the same (MS) error. The necessary and sufficient condition is from (5) when all reflectivities arepproportional to the square root of the ffiffiffiffiffi area, i.e. Ri p Ai for all i: Observe also that when Ai -0 probability (5) also decreases to zero. This is because for decreasing area the form factor pffiffiffiffiffi becomes proportional to the area, and thus psi -Ri Ai : The optimal MSE value is obtained by substituting probabilities (5) into (4) ! 1 X Ri Fsi F2s RF pffiffiffiffiffi G pi ffiffiffiffiffisi ; EðMSEÞ ¼ ð6Þ AT i Ai Ai
DðMSEÞ EðMSEopt Þ ¼1 EðMSEÞ EðMSEÞ pffiffiffiffiffi pffiffiffiffiffi 2 P ioj ðRi = Ai Rj = Aj Þ Fsi Fsj ¼ P : 2 2 ioj ðRi =Ai þ Rj =Aj ÞFsi Fsj
ð9Þ
pffiffiffiffiffi We see again from this formula that when Ri p Ai for all i the gain is null. Also, from the numerator in (9), the more disequilibrium the scene has, the more the gain, where by disequilibrium we mean small patches with high reflectivity and bigger ones with low reflectivity. 2.3. Many sources case Suppose we have 1; y; nS sources each with power Fs : We have to distribute the N rays available optimally between the different sources. Here we consider that the cost is the same for each ray, independent of the source. The total reflected radiosity for patch i; Bi ; is the sum of the reflected radiosities due to the different sources X Bi ¼ Bis : ð10Þ s
On the other hand, the different estimators for the Bis quantities are independent because they use different rays each. Thus we can write, using Ns rays for source s: X Ai X 1 VarðB# is Þ EðMSEÞ ¼ AT s Ns i X 1 X Ai VarðB# is Þ ¼ Ns i AT s X 1 EðMSEs Þ ð11Þ ¼ Ns s P with the restriction s Ns ¼ N: Again we can optimize, and obtain the solution pffiffiffiffiffiffiffiffiffiffiffiffiffi Ns p MSEs : ð12Þ
A. Brusi et al. / Computers & Graphics 26 (2002) 351–354
Now we consider the optimal MSE obtained for each s (6) vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u uX Ri Fsi F2 Ri Fsi s t pffiffiffiffiffi Gs pffiffiffiffiffi Ns p i Ai Ai sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q ffiffiffiffiffiffi X Ri Fsi pffiffiffiffiffi Gs ¼ Fs Gs2 ¼ Fs Gs : E Fs ð13Þ i Ai This solution differs by the factor Gs from the standard source probability selection. 2.4. An adaptive algorithm to use the optimal probabilities To use the probabilities given in (5) in an algorithm means to know the form factors. This makes it apparently impractical. However, we can devise an adaptive algorithm using hemispherical probabilities (see for instance [5]). We divide the hemisphere into equal cosine weighted regions. Each region is sampled initially with equal p weight, ffiffiffiffiffi and counters are used for accumulating the Ri = Ai value found for each ray. At the end of this preliminary step the counter for each region will give the new weights for sampling. Equally, the weight of each source can be determined adaptively. We first distribute the rays according to the power Fs and determine Gs by averaging all contributions for all rays. That is, if rsi are the number of rays that go to patch i from source s and rs the number of rays from source s then 1 X Ri rsi pffiffiffiffiffi ð14Þ Gs E rs i Ai
353
was found to be around 0.75, and 0.65 and 0.63 for (b) and (c). This strange ‘‘extra’’ gain is due to the fact that just stratification in itself induces an important gain for this scene. If we then compare with the stratified case (with equal weight to each hemispherical region) the gain is as expected. The number of hemispherical regions considered was 32 ð4 8Þ; but similar results were also obtained for just 16 regions. Also, the number of rays to gain information is not critical, obtaining similar results for the range 103 –106 : To test now for the impact on the visible quality of the images obtained we have used the table and chairs scene (Fig. 1). The top image corresponds to the non-optimal case, the bottom one to the optimal case. Fig. 2 shows a close-up of the table. Both optimal and non-optimal were obtained casting 2 105 lines. The number of regions was also 32. To find out the average gain, 50 runs were made against a reference scene giving an average gain of 12% (0.88 MSE ratio) for the whole scene, while 13.6% for table, chairs and floor and 9.5% for the walls. The reflectivity of the walls was 0.7, and 0.3 everywhere else. To test for small patches, we divided the top of the table more densely. The results expected were an improvement on regions of higher reflectivity (walls) and smaller area (table), and a decreased quality on low reflectivity and bigger area patches (floor). Now, although the increase in quality in the walls is not obvious from the images in Fig. 1 (except for the upper row of patches near the ceiling), a reduction of noise is clearly visible on the table top in Fig. 2, while the loss in precision in the floor is attenuated. In this test scene, a pure stratification strategy only gave a 2.2% improvement.
because rsi =ri approximates Fsi : The weights are progressively refined with each new step, and the direct illumination radiosity solutions are combined in the usual way using weights proportional to the number of rays cast [7]. Although this is a safe combination, it is probably not the best one because the last steps use more information than the first ones. Other combinations could be tried giving more weight to the latter iterations.
3. Results We have analyzed the optimal probabilities (5) with different test scenes. The first one is a cube divided into 54 equal patches, with one source in the middle of a face. The reflectivities for the five additional faces were 0:3; 0:4; 0:5; 0:6; 0:7 in case (a), 0:1; 0:3; 0:5; 0:7; 0:9 for case (b) and 0:2; 0:8; 0:2; 0:5; 0:2 for (c). The theoretically optimal MSE ratio was found to be 0:92; 0:73 and 0.70, respectively. But very interestingly the gain in the simulations was higher, especially for case (a), which
Fig. 1. Comparison of the optimal solution (b) with standard solution (a) for the test scene with 2 105 rays cast in each.
354
A. Brusi et al. / Computers & Graphics 26 (2002) 351–354
trian–Spanish Acciones Integradas HU 1998-0015 and HU2000-0011 and Catalan-Flemish Joint Action ABM/ acs/ACI98-19.
References
Fig. 2. Top view of the test scene with standard (a) and optimal (b) solutions with 2 105 rays cast in each.
Although more important improvements can be obtained when small patches have higher reflectivity, this simple test scene suggests its use whenever there is an important area disequilibrium in the patches.
4. Conclusions and future research We have presented in this paper a novel strategy in Monte Carlo radiosity to send more rays where they are more needed, that is, patches with high reflectivity and/ or small area. Our approach is based on the optimization of the mean square error. Although the results obtained have been in good correspondence with visual quality, future research will be addressed towards optimizing further error measures and looking for their impact on the image.
Acknowledgements This project has been funded in part with Grant TIC98-0586-C03-02 of the spanish Government, Aus-
[1] Feda M, Purgathofer W. Progressive ray refinement for Monte Carlo radiosity. Fourth Eurographics workshop on rendering, Paris, France, June 1993. p. 15–26. Eurographics Technical Report Series EG 93 RW. [2] Keller A. Quasi-Monte Carlo radiosity. Proceedings of the Seventh Eurographics Workshop on Rendering, Rendering techniques’96. Wien: Springer, 1996. p. 101–10. [3] Jolivet V, Plemenos D, Sbert M. A pyramidal hemisphere subdivision radiosity. Formal definition and improvements. WSCG’2000, Plzen, Czech Republic, 2000. p. 228–35. [4] Neumann L. Monte Carlo radiosity. Computing 1995;55(1):23–42. [5] Neumann A, Neumann L, Bekaert P, Willems Y, Purgathofer W. Importance-driven stochastic ray radiosity. Proceedings of the Seventh Eurographics Workshop on Rendering, Rendering techniques ’96. Wien: Springer, 1996. p. 111–22. [6] Pattanaik SN, Mudur SP. Computation of global illumination by Monte Carlo simulation of the particle model of light. Proceedings of the Third Eurographics Workshop on Rendering, Bristol, UK, May 1992. p. 71–83. [7] Sbert M, P!erez F, Pueyo X. Global Monte Carlo. A Progressive Solution. Rendering Techniques ’95, Proceedings of the Sixth Eurographics Workshop on Rendering, Wien: Springer, 1995. p. 231–9. [8] Sbert M. Error and complexity of random walk Monte Carlo radiosity. IEEE Transactions on Visualization and Computer Graphics 1997;3(1):23–38. [9] Shirley P. A ray tracing method for illumination calculation in diffuse-specular scenes. Proceedings of Graphics Interface ’90, Toronto, Ont., May 1990. p. 205–12.