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.
The basic Monte-Carlo ray tracing (MCRT) generates or transforms rays at random with probability distribution determined locally i.e. in case of emission by the light source properties and in case of surface or particle scattering by Bidirectional scattering distribution function (BSDF) in the hit point.
This approach is not ultimately optimal, e.g. BSDF may send rays away from the virtual observer or camera. Moreover, in case of a parallel light source there is an ambiguity with the ray origin: formally for a parallel light the ray origin must be chosen uniformly in an infinite domain which is technically impossible. Therefore in practice the ray origin is chosen over the projection of scene onto the plane orthogonal to the ray direction (see Fig. 1).
Below we shall describe how it must be adapted to the more complex cases, when
Light source is not parallel but it is skylight with given goniogram;
There are several windows that can overlap.
The examples shown in Fig. 1 and 2 can be implemented by letting the user to point the windows and then mask the distribution of ray origin by 1 inside or 0 outside their projection. This approach, however, may introduce error. For example, imagine a house with “French windows” stayed on a white sand lawn. The sun rays that hit the ground near the window can reflect into it and contribute to the illumination of the interior. Meanwhile, the light emits only towards the transparent windows thus no rays are directed to the opaque ground (Fig. 2).
Fig. 1. Top view of a scene illuminated by a parallel light source. Rays from the light are emitted to the scene only.
Rays that miss the window do not contribute to the
interior view.
Therefore the ray emission is determined not only by the light source properties but also by the geometry of scene. The natural extension of that approach is when we render an interior of a room illuminated from outside through a window (or an open door). In this case the rays sent to the opaque walls are “lost” while illumination of the interior is created only by the rays that run through the window. The ray origin must then be chosen in a projection of that window onto the plane orthogonal to the illumination direction.
This is the basic idea behind optimization of light emission. 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.
This situation can be cured if the “emission mask” is more sophisticated and is not 0 for some opaque areas. On the contrary, the mask is not 1 for some transparent areas even they do not open into the interior.
The resulting Probability density function or PDF of the emitted ray origin (and/or direction) cannot be simply postulated but must be calculated from the results of simulation. The method then belongs to a wide class of algorithms of an “optimal PDF” for MCRT.
Usually the quantities of interest for optical simulation, e.g. the luminance of the hit point, are usually the integrals of a product, 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 scattered ray direction in MCRT must be that product. This is usually termed as importance sampling. The problem is that at one multiplier of that product (illumination in the above example) is the result of MCRT itself and is not known in advance and accurately.
Since long attempts are made to improve efficiency of MCRT by making the PDF of scattered rays as close as possible to that product using various approximations and heuristic. Say, surface luminance under even direct illumination by an area light is an integral over that light surface. To improve efficiency of its calculation Shirley et al [1] suggested to separate the fastest varying terms of the integrand and to use them for constructing a PDF.
The famous works [2], [3] improved efficiency of the backward MCRT by scattering the camera ray according to not BSDF only but its product with the illumination of the point. In [2] photon maps were used to estimate the angular distribution of illuminance. In [3] this was made with a 5D tree covering the scene which accumulated illumination brought by various reflected backward MCRT rays (so there was no a separate forward MCRT phase).
There are two main problems with the method. Firstly the optimal distribution must now be stored if not in each hit point but in each voxel. This requires a lot of memory. Various approximations are used to reduce it, e.g. the space resolution is quite coarse. And this results in the second problem: we have only a rather rough estimation of the incoming radiation distribution. It thus deviates from the target product and the difference must be compensated by the scale of ray energy.
At last, we may face a deadlock with the methods like [3] where the illumination is taken from the same MCRT rays which were shot with the accumulated PDF. For example, if by some time no rays were shot in some angular cone the estimated illumination from it is 0. And when generating a next ray the sampling procedure 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 when some fraction of rays are scattered by BSDF without usage of the accumulated PDF. The difficulty though is how to find the ratio. If it is too small the narrow peaks of illumination do not contribute to PDF and are only taken from the small fraction of rays scattered by BSDF which creates strong noise. If the ratio is large this kills the very benefit (acceleration) of the method.
This basic idea was then exploited by many researchers [4]–[8].
While [5] uses hemispherical PDF of the form rather similar to Jensen’s [2], in [6] the authors use approximation by a sum of Gaussians. Its advantage is that a product of two Gaussian bells is also a Gaussian and can be sampled efficiently. In [8] wavelets are used instead. Some authors utilize spherical harmonics whose product also admits analytic treatment. Also other approaches are used up to neural networks and machine learning which are another way to optimize the ray direction PDF on the base of accumulated statistic [7].
All the above methods alter PDF in one point to reduce variance. A radical solution is the Metropolis Light transport where, unlike MCRT which traces ray segment-by segment, we 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.
Less radical solution is still using MCRT but applying different
approaches (both for generation of scattering direction and for
collecting illuminance) and generating several ray paths. The pixel
luminance is then taken as a weighted sum of luminance each of
them would brought. This is usually termed as multiple importance
strategy. This method does not optimize PDF in hit points
separately but instead it mixes contributions of several different
(
Recently Sbert et al [11] claimed this is not always so. They analyzed the conditions in the Veach’s theorem and found that sometimes an estimator can be constructed significantly better than the “balance heuristic” one.
For a parallel light source the ray origin is distributed uniformly. On the other hand rays that miss the window are useless because they do not create interior illumination. Therefore the ray origin is chosen uniformly within projection of the window onto the plane orthogonal to the light direction.
In case of triangulated scene geometry the triangles that constitute the window glass are used. First we choose one of them with probability proportional to the energy flow through that triangle and then choose the ray origin uniformly within the triangle.
Although it is possible that the window is open or does not have glass inside. In this case there are no triangles inside the window to direct rays to. The user then can select triangles that form the window’s frame. After that a convex hull is constructed from them and the rays are directed to its triangles. In case light source is not parallel the calculation is a bit more sophisticated. Say, sky light, like a parallel light, is an infinitely distant, but it has a smooth goniogram.
Here we first choose direction of emission and after that we choose triangle t with probability =
|( ⋅ )| ∑ |( ⋅ )| where St is area of the t-th triangle and nt is its normal. Finally the ray origin is chosen uniformly in that triangle and translated outside the scene domain. Direction is chosen with probability density ( ) =
∑ |( ⋅ )| ( ) where g is intensity of the emission goniogram.
Let us assume that sky goniogram is tabulated on a rectangular grid ( , ), and bilinearly interpolates inside cells.
We first choose goniogram cell [ , +1 ] × [ , +1] with probability , = ∑ ∫ ,
|( ⋅ )| ( ) 2 and then choose within this cell a direction with probability density , ( ) =
∑ |( ⋅ )| ( )
∑ ∫ , |( ⋅ )| ( )
After that we choose triangle t with probability (
If the cell is small so that the emission directions for all its 4 vertices are to the same side of triangle normal then the integral over cell is
∫ |( ⋅ )| ( ) 2 = |( , ⋅ )| depends only on the goniogram cell (but not triangle). Therefore for a goniogram
cell it suffices to calculate three values , ( ) , , ( ) , ,
( ) and then for most of triangles the flow through them is calculated as a dot product of this vector with triangle normal. For those few triangles which are illuminated from different sides by the goniogram cell the energy flow must be calculated directly.
Therefore in the above method the probability density of the ray origin and direction ( , ) is ( , ) = =
∑ ( ; ) ( ) ∫∑ ( ; ) ( ) 2 2
∑ ( ; ) ( ) where ( ; ) = 1 if the ray ( , ) intersects triangle t and 0 calculated as 3.2
The proposed above method produces correct energy flow
through a surface composed by triangles (or other facets) whose
projections do not overlap. If projections of N triangles overlap the
numerator in (
As seen from (
After the light source ray ( ; ) is chosen its further fate is still
stochastic at each hit point. So its contribution to the luminance of
pixel p is random function ( ; ). Naturally if one changes the
distribution of rays it must be compensated by the change of the ray
contribution to keep the expectation. Therefore
( ; ) = ℒ ( ; )
( ; )
( ; )
where ℒ ( ; ) is contribution for the correct density (
Assuming the camera ray path is deterministic, the noise in pixel p is entirely due to the forward MCRT part and can be ∫〈 2 〉( ; ) ( ; ) 2 2 − (∫〈 〉( ; ) ( ; ) 2 2 ) where 〈⋅〉 means averaging over all random fates of the ray emitted from light source with origin and direction . The total error is just the sum over all pixels.
Varying the probability density by changes this total error − ∫ ∑〈ℒ2〉( ; )
( ; ) 2 2 2( ; ) 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 (
This correcting factor is easily handled with the “rejection
method”: we first choose ray direction and origin according to the
probability density (
will then be (
The 〈ℒ2〉( ; ) is by definition independent on how we emit the ray ( ; ). Before ray tracing we create a 4D mesh in that ( ; ). Then when a ray is emitted we find which cell the ( ; ) belongs to and add ℒ2( ; ) to the value accumulated in that cell.
As ray tracing continues we accumulate better and better
estimation of 〈ℒ2〉( ; ) in each cell and begin using it to optimize
emission density. The ray origin and direction are now chosen with
density (
we choose the ray start. Gradually, as the accumulated estimate of 〈ℒ2〉( ; ) improves, the density of ray start converges to its limiting optimal form.
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.
Figures 5–7 present the camera images for “Two rooms” scene calculated during the same time interval without light emission optimization (Fig. 5), for light emission towards the gap beneath the door, i.e. algorithm described in p. 3 (Fig. 6), and for adaptive light PDF (Fig. 7). We also calculated average illuminance and RMS over the red rectangle presented in Fig. 5-7. The calculated values are in Table 1.
Average 0.49 0.23 0.47
RMS 0.62 0.11 0.27
Fig. 5. “Two rooms” scene. Light emission to the whole scene domain, without optimization. beneath the door.
Fig. 7. “Two rooms” scene. Generation of ray direction from the light source uses an adaptive PDF
One can see that the least noise is in the second method but it underestimates luminance outside of the bright area on the floor, for example, in the wall area marked by red rectangular. Besides the gap in it looks entirely black because camera rays sent to it go to the right room which is entirely black because illumination goes to the gap only.
Another example scene is shown in Fig. 8. The house is illuminated with sky and sun light. There is a lot of windows through which it can penetrate into the interior of the room observed by the camera. Therefore the second method (explicit specification of the windows in the scene) is inconvenient here.
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.
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 reduction but in some cases it may underestimate illuminance as it can be seen from Table 1. The second method constructs an optimal PDF of light emission which minimizes the noise. It is slower but gives a correct illuminance. Nevertheless the second method accelerates image generation in times in comparison to light emission without optimization to the whole domain.
Fig. 9. Camera image for the “House” scene. Light
emission to the whole scene domain.
Fig. 10. Camera image for the “House” scene. Light
emission by the adaptive light PDF.
This work was supported by RFBR grants 18-01-00569 and 19-01-00435.