Optimization of Illumination through Windows for MCRT Ildar Valiev1, Dmitry Zhdanov2, Sergey Ershov1, Alexey Voloboy1, Vadim Sokolov1, Vladimir Galaktionov1 piv@gin.keldysh.ru|ddzhdanov@mail.ru|sergey_65@mail.ru| voloboy@gin.keldysh.ru|sokolov1969@gmail.com|vlgal@gin.keldysh.ru 1 The Keldysh Institute of Applied Mathematics Russian Academy of Sciences, Moscow, Russia 2 ITMO University, St. Petersburg, Russia The paper presents an improvement of Monte-Carlo ray tracing which changes ray emission from a light source to accelerate convergence i.e. reduce the noise remained after the given run time. It is mainly intended for interior scenes illuminated from outside (e.g. skylight) through windows or other holes. The rays from light sources are generated so that they are directed to these windows. In other words, the number of rays is increased for directions that contribute to the camera image. It is shown that the proposed method allows calculating image with desirable quality several times faster. Keywords: ray tracing, lighting simulation, Monte-Carlo, optimal PDF. Below we shall describe how it must be adapted to the more 1. Introduction complex cases, when The basic Monte-Carlo ray tracing (MCRT) generates or  Light source is not parallel but it is skylight with given transforms rays at random with probability distribution determined goniogram; locally i.e. in case of emission by the light source properties and in  There are several windows that can overlap. case of surface or particle scattering by Bidirectional scattering The examples shown in Fig. 1 and 2 can be implemented by distribution function (BSDF) in the hit point. letting the user to point the windows and then mask the distribution This approach is not ultimately optimal, e.g. BSDF may send of ray origin by 1 inside or 0 outside their projection. This rays away from the virtual observer or camera. Moreover, in case approach, however, may introduce error. For example, imagine a of a parallel light source there is an ambiguity with the ray origin: house with “French windows” stayed on a white sand lawn. The formally for a parallel light the ray origin must be chosen uniformly sun rays that hit the ground near the window can reflect into it and in an infinite domain which is technically impossible. Therefore in contribute to the illumination of the interior. Meanwhile, the light practice the ray origin is chosen over the projection of scene onto emits only towards the transparent windows thus no rays are the plane orthogonal to the ray direction (see Fig. 1). directed to the opaque ground (Fig. 2). Fig. 2. Side view of a scene illuminated by a parallel light source. Rays are directed to the windows only and not to the ground. These rays (shown dashed) can reflect inside the room too but are missed. Fig. 1. Top view of a scene illuminated by a parallel light This situation can be cured if the “emission mask” is more source. Rays from the light are emitted to the scene only. sophisticated and is not 0 for some opaque areas. On the contrary, Rays that miss the window do not contribute to the the mask is not 1 for some transparent areas even they do not open interior view. into the interior. The resulting Probability density function or PDF of the Therefore the ray emission is determined not only by the light emitted ray origin (and/or direction) cannot be simply postulated source properties but also by the geometry of scene. The natural but must be calculated from the results of simulation. The method extension of that approach is when we render an interior of a room then belongs to a wide class of algorithms of an “optimal PDF” for illuminated from outside through a window (or an open door). In MCRT. this case the rays sent to the opaque walls are “lost” while illumination of the interior is created only by the rays that run 2. Relation to previous work through the window. The ray origin must then be chosen in a projection of that window onto the plane orthogonal to the Usually the quantities of interest for optical simulation, e.g. the illumination direction. luminance of the hit point, are usually the integrals of a product, This is the basic idea behind optimization of light emission. e.g. the point luminance is a convolution of BSDF with Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). illumination. Therefore the optimal probability distribution of the whole ray paths. The main ideas were suggested in [10] where scattered ray direction in MCRT must be that product. This is besides all it was proved that the relatively simple “balance usually termed as importance sampling. The problem is that at one heuristic” estimator is close to optimal, i.e. other estimators of the multiplier of that product (illumination in the above example) is the class considered cannot decrease the variance significantly. result of MCRT itself and is not known in advance and accurately. Recently Sbert et al [11] claimed this is not always so. They Since long attempts are made to improve efficiency of MCRT analyzed the conditions in the Veach’s theorem and found that by making the PDF of scattered rays as close as possible to that sometimes an estimator can be constructed significantly better than product using various approximations and heuristic. Say, surface the “balance heuristic” one. luminance under even direct illumination by an area light is an integral over that light surface. To improve efficiency of its 3. Light entrance through the windows calculation Shirley et al [1] suggested to separate the fastest varying terms of the integrand and to use them for constructing a PDF. 3.1 Simple transparent window The famous works [2], [3] improved efficiency of the backward For a parallel light source the ray origin is distributed MCRT by scattering the camera ray according to not BSDF only uniformly. On the other hand rays that miss the window are useless but its product with the illumination of the point. In [2] photon maps because they do not create interior illumination. Therefore the ray were used to estimate the angular distribution of illuminance. In [3] origin is chosen uniformly within projection of the window onto this was made with a 5D tree covering the scene which accumulated the plane orthogonal to the light direction. illumination brought by various reflected backward MCRT rays (so In case of triangulated scene geometry the triangles that there was no a separate forward MCRT phase). constitute the window glass are used. First we choose one of them There are two main problems with the method. Firstly the with probability proportional to the energy flow through that optimal distribution must now be stored if not in each hit point but triangle and then choose the ray origin uniformly within the in each voxel. This requires a lot of memory. Various triangle. approximations are used to reduce it, e.g. the space resolution is Although it is possible that the window is open or does not have quite coarse. And this results in the second problem: we have only glass inside. In this case there are no triangles inside the window to a rather rough estimation of the incoming radiation distribution. It direct rays to. The user then can select triangles that form the thus deviates from the target product and the difference must be window’s frame. After that a convex hull is constructed from them compensated by the scale of ray energy. and the rays are directed to its triangles. In case light source is not At last, we may face a deadlock with the methods like [3] where parallel the calculation is a bit more sophisticated. Say, sky light, the illumination is taken from the same MCRT rays which were like a parallel light, is an infinitely distant, but it has a smooth shot with the accumulated PDF. For example, if by some time no goniogram. rays were shot in some angular cone the estimated illumination Here we first choose direction of emission 𝒗 and after that we from it is 0. And when generating a next ray the sampling procedure choose triangle t with probability does not send it to that “black” region, so illumination from it 𝑺𝒕 |(𝒗 ⋅ 𝒏𝒕 )| remains 0 and so on. The simplest remedy is to use a mixed strategy 𝑷𝒕 = (1) ∑𝒕 𝑺𝒕 |(𝒗 ⋅ 𝒏𝒕 )| when some fraction of rays are scattered by BSDF without usage of the accumulated PDF. The difficulty though is how to find the where St is area of the t-th triangle and nt is its normal. Finally the ratio. If it is too small the narrow peaks of illumination do not ray origin 𝒙 is chosen uniformly in that triangle and translated contribute to PDF and are only taken from the small fraction of rays outside the scene domain. Direction 𝒗 is chosen with probability scattered by BSDF which creates strong noise. If the ratio is large density this kills the very benefit (acceleration) of the method. ∑𝑡 𝑆𝑡 |(𝒗 ⋅ 𝒏𝑡 )|𝑔(𝒗) 𝑝(𝒗) = This basic idea was then exploited by many researchers [4]–[8]. ∑𝑡 𝑆𝑡 ∫|(𝒗 ⋅ 𝒏𝑡 )|𝑔(𝒗)𝑑2 𝒗 While [5] uses hemispherical PDF of the form rather similar to where g is intensity of the emission goniogram. Jensen’s [2], in [6] the authors use approximation by a sum of Let us assume that sky goniogram is tabulated on a rectangular Gaussians. Its advantage is that a product of two Gaussian bells is grid (𝜗, 𝜑), and bilinearly interpolates inside cells. also a Gaussian and can be sampled efficiently. In [8] wavelets are used instead. Some authors utilize spherical harmonics whose We first choose goniogram cell [𝜗𝑖 , 𝜗𝑖+1 ] × [𝜑𝑗 , 𝜑𝑗+1 ] with product also admits analytic treatment. Also other approaches are probability used up to neural networks and machine learning which are another ∑𝑡 𝑆𝑡 ∫𝑐𝑒𝑙𝑙 |(𝒗 ⋅ 𝒏𝑡 )|𝑔(𝒗)𝑑2 𝒗 way to optimize the ray direction PDF on the base of accumulated 𝑖,𝑗 (2) 𝑃𝑖,𝑗 = statistic [7]. ∑𝑡 𝑆𝑡 ∫|(𝒗 ⋅ 𝒏𝑡 )|𝑔(𝒗)𝑑2 𝒗 All the above methods alter PDF in one point to reduce and then choose within this cell a direction with probability density variance. A radical solution is the Metropolis Light transport where, unlike MCRT which traces ray segment-by segment, we ∑𝒕 𝑺𝒕 |(𝒗 ⋅ 𝒏𝒕 )|𝒈(𝒗) 𝒑𝒊,𝒋 (𝒗) = (3) choose at random the whole path from light source to camera [9]. ∑𝒕 𝑺𝒕 ∫𝒄𝒆𝒍𝒍 |(𝒗 ⋅ 𝒏𝒕 )|𝒈(𝒗)𝒅𝟐 𝒗 𝒊,𝒋 But while the above problem of MCRT is eliminated other problems emerge whose analysis is out of our scope. After that we choose triangle t with probability (1). Less radical solution is still using MCRT but applying different If the cell is small so that the emission directions for all its 4 approaches (both for generation of scattering direction and for vertices are to the same side of triangle normal then the integral collecting illuminance) and generating several ray paths. The pixel over cell is luminance is then taken as a weighted sum of luminance each of them would brought. This is usually termed as multiple importance ∫ |(𝒗 ⋅ 𝒏𝒕 )|𝑔(𝒗)𝑑2 𝒗 = |(𝑭𝑖,𝑗 ⋅ 𝒏𝒕 )| strategy. This method does not optimize PDF in hit points 𝑐𝑒𝑙𝑙𝑖,𝑗 separately but instead it mixes contributions of several different where 4. Adaptive optimal PDF As seen from (5) the above approach masks the own spatially 𝑭𝑖,𝑗 ≡ ∫ 𝒗𝑔(𝒗)𝑑2 𝒗 uniform light source probability density 𝑔(𝒗) with mask which is 𝑐𝑒𝑙𝑙𝑖,𝑗 1 if the ray intersect any window triangle and 0 otherwise. This depends only on the goniogram cell (but not triangle). Therefore simplest alteration is not the optimal. Let us consider a more for a goniogram cell it suffices to calculate three values flexible approach when the origin and direction of the emitted ray (𝑥) (𝑦) (𝑧) 𝐹𝑖,𝑗 , 𝐹𝑖,𝑗 , 𝐹𝑖,𝑗 and then for most of triangles the flow through them are chosen with arbitrary probability density 𝜌(𝒙; 𝒗). is calculated as a dot product of this vector with triangle normal. After the light source ray (𝒙; 𝒗) is chosen its further fate is still For those few triangles which are illuminated from different sides stochastic at each hit point. So its contribution to the luminance of by the goniogram cell the energy flow must be calculated directly. pixel p is random function 𝐿𝑝 (𝒙; 𝒗). Naturally if one changes the Therefore in the above method the probability density of the distribution of rays it must be compensated by the change of the ray ray origin and direction (𝒙, 𝒗) is contribution to keep the expectation. Therefore 𝑝(𝒙; 𝒗) ∑𝑡 𝜒𝑡 (𝒙; 𝒗)𝑔(𝒗) 𝐿𝑝 (𝒙; 𝒗) = ℒ𝑝 (𝒙; 𝒗) 𝑝(𝒙, 𝒗) = 𝜌(𝒙; 𝒗) ∫ ∑𝑡 𝜒𝑡 (𝒙; 𝒗)𝑔(𝒗) 𝑑2 𝒗𝑑2 𝒙 ∑𝑡 𝜒𝑡 (𝒙; 𝒗)𝑔(𝒗) (4) where ℒ𝑝 (𝒙; 𝒗) is contribution for the correct density (5). = Assuming the camera ray path is deterministic, the noise in ∑𝑡 𝑆𝑡 ∫|(𝒗 ⋅ 𝒏𝑡 )|𝑔(𝒗)𝑑2 𝒗 pixel p is entirely due to the forward MCRT part and can be where 𝜒𝑡 (𝒙; 𝒗) = 1 if the ray (𝒙, 𝒗) intersects triangle t and 0 calculated as otherwise. 2 3.2 Multiple windows and overlap ∫〈𝐿2𝑝 〉(𝒙; 𝒗)𝜌(𝒙; 𝒗)𝑑2 𝒙𝑑2 𝒗 − (∫〈𝐿𝑝 〉(𝒙; 𝒗)𝜌(𝒙; 𝒗)𝑑2 𝒙𝑑2 𝒗) The proposed above method produces correct energy flow where 〈⋅〉 means averaging over all random fates of the ray emitted through a surface composed by triangles (or other facets) whose from light source with origin 𝒙 and direction 𝒗. The total error is projections do not overlap. If projections of N triangles overlap the just the sum over all pixels. numerator in (4) increases for this area and this causes a distortion, Varying the probability density by 𝛿𝜌 changes this total error see Fig. 3. by 𝑝2 (𝒙; 𝒗) − ∫ ∑〈ℒ𝑝2 〉(𝒙; 𝒗) 𝛿𝜌(𝒙; 𝒗)𝑑2 𝒙𝑑2 𝒗 (6) 𝜌2 (𝒙; 𝒗) 𝑝 An optimal 𝜌(𝒙; 𝒗) makes the noise minimal and in this extremum the change must vanish for any admissible 𝛿𝜌. The variation of normalized density is by definition an arbitrary function for which ∫ 𝛿𝜌(𝒙; 𝒗)𝑑2 𝒙𝑑2 𝒗 = 0. Comparing with (6) we thus conclude that it must be (the common derivation in Monte- Carlo methods) 𝑝2 (𝒙; 𝒗) ∑〈ℒ𝑝2 〉(𝒙; 𝒗) 2 = const 𝜌 (𝒙; 𝒗) 𝑝 Fig. 3. Top view of a scene with multiple windows. The rays shown with thick arrows intersect two windows. or Dashed arrows show rays directed to windows that were shadowed by an opaque surface before entering the room. 𝜌(𝒙; 𝒗) = 𝐶𝑝(𝒙; 𝒗)√∑〈ℒ𝑝2 〉(𝒙; 𝒗) (7) 𝑝 The obvious correction is to change (4) to where C is the normalization constant. ∑𝑡 𝑁 −1 (𝒙; 𝒗)𝜒𝑡 (𝒙; 𝒗) 𝑔(𝒗) The 〈ℒ𝑝2 〉(𝒙; 𝒗) is by definition independent on how we emit 𝑝(𝒙, 𝒗) = (5) ∫ ∑𝑡 𝑁 −1 (𝒙; 𝒗)𝜒𝑡 (𝒙; 𝒗)𝑔(𝒗) 𝑑2 𝒗𝑑2 𝒙 the ray (𝒙; 𝒗). Before ray tracing we create a 4D mesh in that (𝒙; 𝒗). Then when a ray is emitted we find which cell the (𝒙; 𝒗) where 𝑁(𝒙; 𝒗) = ∑𝑡 𝜒𝑡 (𝒙; 𝒗) is the number of triangles the ray (𝒙; 𝒗) intersects. belongs to and add ℒ𝑝2 (𝒙; 𝒗) to the value accumulated in that cell. As ray tracing continues we accumulate better and better This correcting factor is easily handled with the “rejection method”: we first choose ray direction and origin according to the estimation of 〈ℒ𝑝2 〉(𝒙; 𝒗) in each cell and begin using it to optimize probability density (4), i.e. actually with (1), (2) and (3). Then once emission density. The ray origin and direction are now chosen with cycles over all triangles of all windows and check if the ray density (7). As ray tracing continues and the accumulated estimate intersects them. This is done for isolated triangle, i.e. we do not care of 〈ℒ𝑝2 〉(𝒙; 𝒗) changes, the density of ray emission changes too. In if there are other triangles “in front of it”. Then with probability other words it is now time-dependent. This however does not create 1 − 𝑁 −1 (𝒙; 𝒗) the ray is rejected, i.e. the direction and origin are a feedback loop because the accumulated value does not depend on chosen anew, and this continues until the ray is accepted. The how we choose the ray start. Gradually, as the accumulated resulting density will then be (4) then it provides correct estimate of 〈ℒ𝑝2 〉(𝒙; 𝒗) improves, the density of ray start converges illumination. to its limiting optimal form. 5. Results beneath the door. An example scene is oversimplified model of two rooms separated by a wall with a door. Skylight illuminates the right room through the window and light penetrates into the left room through the gap beneath the door only (Fig. 4). Camera was in the left room and looks towards the door. Fig. 4. Layout of the “Two rooms” scene. Fig. 7. “Two rooms” scene. Generation of ray direction Figures 5–7 present the camera images for “Two rooms” scene from the light source uses an adaptive PDF calculated during the same time interval without light emission optimization (Fig. 5), for light emission towards the gap beneath One can see that the least noise is in the second method but it the door, i.e. algorithm described in p. 3 (Fig. 6), and for adaptive underestimates luminance outside of the bright area on the floor, light PDF (Fig. 7). We also calculated average illuminance and for example, in the wall area marked by red rectangular. Besides RMS over the red rectangle presented in Fig. 5-7. The calculated the gap in it looks entirely black because camera rays sent to it go values are in Table 1. to the right room which is entirely black because illumination goes to the gap only. Average RMS Another example scene is shown in Fig. 8. The house is Without optimization (Fig. 5) 0.49 0.62 illuminated with sky and sun light. There is a lot of windows Light emission toward the gap (Fig. 0.23 0.11 through which it can penetrate into the interior of the room 6) observed by the camera. Therefore the second method (explicit specification of the windows in the scene) is inconvenient here. Adaptive light PDF (Fig. 7) 0.47 0.27 Table 1. Average illuminance and RMS over the red rectangle in figures 5-7 (cd/m2). Camera position indoor Fig. 8. Layout of the “House” scene. Fig. 5. “Two rooms” scene. Light emission to the whole scene domain, without optimization. Figures 9–10 present the images calculated during the same time interval without light emission optimization (Fig. 9) and for adaptive light PDF (Fig. 10). It is seen from these figures that level of noise is noticeably low for the adaptive light PDF method. Also visually the same quality of image was obtained about 4 times faster for the adaptive light PDF. 6. Conclusion We suggested two methods that allow to improve ray tracing for scenes where the area of interest is illuminated from outside through windows or small holes and maybe along a complex path. Both methods reduce noise or, which is the same, allow to achieve its target level faster. The first one, which requires marking explicitly those holes in the scene geometry, gives better noise Fig. 6. “Two rooms” scene. Light emission to the gap reduction but in some cases it may underestimate illuminance as it [4] D. Burke, A. Ghosh and W. Heidrich. Bidirectional can be seen from Table 1. The second method constructs an optimal Importance Sampling for Direct Illumination. Eurographics PDF of light emission which minimizes the noise. It is slower but Symposium on Rendering (2005) gives a correct illuminance. Nevertheless the second method [5] Heinrich Hey and Werner Purgathofer. Importance sampling accelerates image generation in times in comparison to light with hemispherical particle footprints. In Proceedings of the emission without optimization to the whole domain. 18th Spring Conference on Computer Graphics (SCCG '02). ACM, New York, NY, USA, 107-114. (2002) DOI http://dx.doi.org/10.1145/584458.584476 [6] Sebastian Herholz, Oskar Elek, Jiˇrí Vorba, Hendrik Lensch, Jaroslav Kˇrivánek. Product Importance Sampling for Light Transport Path Guiding. Eurographics Symposium on Rendering 2016, 35(4) (2016) [7] Vorba J., Karlík O., Šik M., Ritschel T., Kˇrivánek J.: On-line learning of parametric mixture models for light transport simulation. ACM Trans. Graph. (Proc. of SIGGRAPH) 33 (2014). [8] Petrik Clarberg, Wojciech Jarosz, Tomas Akenine-Möller, and Henrik Wann Jensen. 2005. Wavelet importance sampling: efficiently evaluating products of complex functions. ACM Trans. Graph. 24, 3 (July 2005), 1166-1175. DOI: https://doi.org/10.1145/1073204.1073328 [9] Veach E., Guibas L. J. Metropolis Light Transport. // Fig. 9. Camera image for the “House” scene. Light SIGGRAPH ’97 Proceedings of the 24th annual conference on emission to the whole scene domain. Computer graphics and interactive techniques, 65–76 (1997) [10] Veach E., Guibas L. J. Optimally combining sampling techniques for Monte Carlo rendering. ACM SIGGRAPH 1995, pp. 419–428 [11] Mateu Sbert, Vlastimil Havran and Laszlo Szirmay-Kalos. Multiple importance sampling revisited: breaking the bounds. EURASIP Journal on Advances in Signal Processing (2018) 2018:15 https://doi.org/10.1186/s13634-018-0531-2 Fig. 10. Camera image for the “House” scene. Light emission by the adaptive light PDF. 7. Acknowledgements This work was supported by RFBR grants 18-01-00569 and 19-01-00435. 8. References [1] Peter Shirley, Changyaw Wang, and Kurt Zimmerman. 1996. Monte Carlo techniques for direct lighting calculations. ACM Trans. Graph. 15(1), 1–36 (1996). DOI http://dx.doi.org/10.1145/226150.226151 [2] Jensen, H.W. Importance driven path tracing using the photon map. In Rendering Techniques (proc. EGWR) (1995) [3] Lafortune E.P., Willems Y.D. A 5D Tree to Reduce the Variance of Monte Carlo Ray Tracing. In: Hanrahan P.M., Purgathofer W. (eds) Rendering Techniques ’95. Eurographics. Springer, Vienna (1995). DOI https://doi.org/10.1007/978-3-7091-9430-0_2