=Paper=
{{Paper
|id=Vol-2485/paper8
|storemode=property
|title=Optimization of Illumination through Windows for MCRT
|pdfUrl=https://ceur-ws.org/Vol-2485/paper8.pdf
|volume=Vol-2485
|authors=Ildar Valiev,Dmitry Zhdanov,Sergey Ershov,Alexey Voloboy,Vadim Sokolov,Vladimir Galaktionov
}}
==Optimization of Illumination through Windows for MCRT==
https://ceur-ws.org/Vol-2485/paper8.pdf
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