The classic Monte-Carlo ray tracing is a powerful technique for simulating almost all effects in ray optics, but it may be prohibitively slow for, for example, calculation of images seen by a lens camera. Therefore, in practice there are often used its various modifications, in particular, bi-directional stochastic ray tracing with photon maps. The well-known flaw of all stochastic methods is their noise. The noise level, that is, the root mean square of pixel brightness calculated during given time, depends, besides all, on the number of rays traced from the light source and from the camera. The choice of the optimal parameters must provide the lowest noise level in a fixed time. This article is devoted to the choice of the optimal number of rays that minimizes the noise. It is proved that this minimal noise is in the same time homogeneous over the image. We produce the formulae to calculate the optimal number of rays from several coefficients which can be obtained from a bi-directional ray tracing of several auxiliary variants. It happens that this optimum is rather wide i.e. the noise level changes with the number of rays slowly, which allows to choose it including other factors e.g. limit this number to save memory. realistic rendering, bi-directional Monte-Carlo ray tracing, photon maps, denoising GraphiCon 2021: 31st International Conference on Computer Graphics and Vision, September 27-30, 2021, Nizhny Novgorod, Russia
Nowadays light simulation is widely used in realistic computer graphics and for design of new materials and optical systems [1]. If wave effects may be neglected in simulation, then stochastic ray tracing methods are preferable. This group of methods includes the simulation of light transport using the Metropolis method [2] and stochastic ray tracing [3]. The classical forward Monte Carlo ray tracing (which starts tracing from the light source) is inefficient for image generation, and therefore it is replaced by bidirectional modifications of it [4–6]. One of the most popular methods among them is the bidirectional stochastic ray tracing with photon maps (BDPM – Bidirectional Photon Mapping) [5, 7]. A well-known flaw of all stochastic methods is that they produce noisy results. Therefore, the noise reduction problem always has primary importance. It is considered in many works, e.g., see [8–10].
The level of noise in BDPM mainly depends on the random scattering of the forward and backward rays, on the choice of the vertex for their merging (or, in other words, on the vertex of the camera ray trajectory at which luminance is estimated from photon maps), and also on the number of forward and backward rays traced in one iteration step. The majority of studies is devoted to the first two issues (e.g., [9–12]), and the number of rays got less attention. However, this is an important factor, and it often happens that the number of forward rays is already redundant and its further increase only increases the computation time but does not reduce the noise. In other cases, it may happen that the number of forward rays is indeed critical, while the number of traced backward rays is redundant. It is usually difficult to predict which fraction of the forward and backward rays is optimal: however, a good choice can speed up the computations by several times.
2021 Copyright for this paper by its authors.
In this paper, we propose a method which allows to estimate the optimal number of rays. A general rule determining noise in the BDPM was derived in [8].
One must realize that the rate of convergence of stochastic methods like BDPM is not actually determined the variance of the image calculated in one iteration (or another number of them), but the variance after the fixed time of calculations. Indeed, if the variance of one iteration is decreases twice while the time to calculate it increases tenfold, the situation obviously worsens.
Therefore, calculation of the optimum i.e. those number of camera and light paths that provide “the best image” (or the fastest convergence) we must account for both the variance from one iteration and the time spent on one iteration.
In this paper we produce the simple approximating law for that time and, combining it with the law of dependence of the variance from one iteration on the number of rays, we derive the formulae for the optimal number of rays. It happens that those distribution of camera rays per pixel that makes the noise homogeneous is just the one which makes it minimal. This is especially convenient because eliminates the problem of which spatial aggregate to treat as the “noise amplitude”.
Generally the number of camera rays traced through a pixel may depend on this pixel, i.e. vary across the image. The variance of the contribution of one iteration of BDPM can be calculated individually for each pixel and is [11, 13]
( ) ( )≡ ⟨⟨ 2⟩⟩( )− ⟨⟨ ⟩⟩2( ) ( ) ≡ ⟨⟨ ⟩2 ⟩ ( )− ⟨⟨ ⟩⟩2( ) ( )≡ ⟨⟨ ⟩2 ⟩ ( )− ⟨⟨ ⟩⟩2( ) respectively. where is the number of light paths per iteration, ( )is the number of camera paths through the pixel p per iteration, ( ( ), ( ))is the contribution to the pixel’s luminance from merging of the light path ( ) with camera path ( ), its average ⟨⟨ ⟩⟩ obviously equals the limiting (exact) pixel’s luminance ( ), and ⟨⋅⟩ and ⟨⋅⟩ denote the averaging over (the ensemble of) light and camera paths,
The variance of the image obtained during time Т is obviously where is the (average) time spent on one iteration. Therefore the value that determines convergence rate and which, therefore, must be minimized to obtain the best image, is parts: ● ● ● path. ( )≡ ( ) The total iteration time is the sum of these partial timings denoted as , , respectively.
The time spent on the first parts is obviously linear in the number of rays. For the third part, it is linear in the number of their pairs. The general functional form is then
Let’s estimate the average time per iteration. The iteration in BDPM consists of three successive
= ∑
( ) ( ),
= ,
=
(
Summing, we have = + + = + ∑ ( )) ( + ∑
( ) ( ) + ∑
( ) ( )) Usually the number of rays is much greater than 1, so we can neglect the terms −1( )и −1: ( )= ( 1 ( + ( )) + )( + ∑
is homogeneous i.e.
Usually the noise is not homogeneous over the image, i.e. it depends on the pixel. It is not very convenient because it is difficult to decide, which is better to have less noise on the floor yet stronger in the ceiling or vice versa? Therefore usually we use a single measure of noise, e.g. the maximum over the image, which better corresponds to the visual perception of “noisy” or “clear” images.
Therefore, what we need to minimize is ( ).
It happens that (it will be so should ( )be continuous instead of integers) in this case the noise Indeed, let us suppose the opposite: that in some pixel ∗ the noise is below the maximum i.e. ( ∗)< ( ). Now let us decrease (infinitesimally as now we treated it as a continuous parameter!) ( ∗). Obviously, in all the other pixels but ∗ the noise will decrease after it because of decrease of the time spent on one iteration : ( )= ( ) ≡ while ( )for these pixels won’t change. In the same time in pixel ∗ the change of noise level will be ( ∗)= − 2( ∗) 1 (
+ ( ∗)) ( ∗)+ ( ∗)
. Therefore, ( ) = ( ) and In view of our supposition ( ∗)< ( ) the maximum is achieved outside ∗ where for each pixel ( ) − ( ∗)= ( ( ) − ( ∗)) +
+ ( ∗)) ( ∗) (through those pixels where the noise is below the maximum) we can make it more optimal (smaller). In other words, an inhomogeneous distribution cannot be the optimal (since it admits improvement). Therefore, the optimal noise distribution is homogeneous.
In fact the limitation that ( )must be integer distort this simple solution, but in case ( )is large enough, the effect must be weak. One may not though hope to reach the homogeneous distribution with small number of camera rays, so it is better to have it large enough so that it be at least 3 in the “worse” pixels (recall that the roundoff error is below half unit!).
So we should estimate the number of camera rays that yields a homogeneous distribution of noise. For the sake of simplicity let us consider the case when 1 ( )is negligible (as compared to the other noise components) which for the case of BDD>0 is if not always then very frequently. Then and one can see that the homogeneous distribution of noise is achieved when ( )= + ( )) = ( + ( )) where , denote the image size in pixels and × ( 1
1 ≡ ∑ ( ), ≡ ∑
Now we have just two free parameters: and the scale factor for the number of camera rays (notice the actual average number of camera rays per pixel will deviate from it because of rounding) which must be varies to minimize the (now homogeneous!) noise. 4. The optimal and has
Substituting (
+
(
One easily finds that for the given , minimum in is achieved for
Denoting we can write this minimum as ( + which is achieved for = ( 2 − ̂ ) + √ ̂ + + ̂ + ̂ = 2 + 2 √ ̂ − + ̂
+ ̂ ̂
Obviously, this is a monotone increasing function of and thus a monotone decreasing function of
. Therefor there is no extremum in and formally the greater this value, the better. However too
many rays per pixels is bad as concerning memory etc. Therefore it is reasonable to take such that
the noise is only (1 + ), where is a small number, times the limiting noise value for → ∞ i.e.
2 + 2 √ ̂
−
̂
+ ̂
(
−√ ̂ + √ ̂ + ( ̂ (We naturally discarded the negative root) or =
̂
2 − ̂
=
)
2
− ̂
(
+ 2√ ̂ √ ̂ (√ ̂ +
2
√ ̂ )
Substituting it in (
The value of , i.e. time spent on average on one light (forward) ray is rather easy to measure. We trace some large number of FMCRT rays , measure the time spent and divide it by .
But ( )and ( )are more difficult to obtain. In principle, we can perform a series of BMCRT simulations for single pixel each, tracing a large number of rays through this pixel and then dividing the time spent by their number. This is expensive, though, because for an accurate measurement of time requires that it be at least milliseconds, while there usually is millions of pixels.
The value of ( )also can be measured straightforwardly, but this is even more difficult.
Happily, we do not need all this because the expressions (
By definition, the two former are BMCRT time and merging time when there are ( )camera rays per pixel and one FMCRT ray. The two latter are for ( )camera rays per pixel.
These values can be measured easily just by doing a few iterations of bi-directional ray tracing for some specially chosen ( )and some (more or less arbitrary) medium about, say, 10000.
Namely, let us take ( )= × ( )where the constant is chosen so that the average number
of rays per pixels be about, say, 10. Then we measure the time spent on the backward ray tracing ,
on the forward ray tracing and time spent on merging the paths . In view of (
( ) ( )= ∑
Similarly, choosing ( )= × ( )we can measure the two remaining sums.
We use the famous Cornell Box with an isotropic point light source slightly below the center of the ceiling. All scene surfaces are gray Lambert with albedo kd=50%. Both direct and indirect illumination was taken from photon maps. “Diffuse depth” i.e. the maximal allowed number of diffuse events for camera ray [11] was BDD=1. Scene image is shown in Figure 1.
There are two variants which differ in the radius of integration sphere R: “large” one equals 0.0083 of the scene size, and the “small” is 0.0015 of the scene size. Naturally, one can expect that smaller radius requires more FMCRT rays.
Calculation of ̂ and ̂ was done by the usual bi-MCRT, using radius 0.0083 of the scene size, 10000 iteration had been done with = 10000 light paths and ( )= 25 camera rays per pixel. The matrices ( )and ( )were calculated as described in [13]. The sums are presented in Table 1. Averages of the BMCRT term and cross-term for the two test variants
R = 0.0083 of scene size
Calculation of time-related terms was done like described in Section 5 for 10000 light paths and average of 50 camera paths per pixel. Obviously the pure BMCRT and FMCRT terms are independent of radius. Since measurement of time is not very accurate (because of loading by background processes, adaptive changing of CPU clock etc), the figures are rounded (Table 2). Time-related sums for the two test variants R = 0.0083 of scene size
R = 0.0015 of scene size
̂
̂
̂
̂
coefficients taken from Table 2.
Approximation of the iteration time (
One can see that approximation is quite good for large but not for small . Mainly this is real run time is larger. Finally the optimal numbers of rays are in Table 5. because (9) did not include rounding of the number of camera rays which increase ( )and thus the
REMARK. Rather large times are because we did not optimized the calculations and did not use
The number of rays for the two test variants
Obviously, the small estimates for is rather absurd since it yields less than one ray per pixel
everywhere. Small ( )results in strong rounding and the error of rounding which can reach ±0.5 is
a serious distortion when ( )is about 3. So we used = 5 as the low bound. As explained above,
there is no formal optimum i.e. the greater the better, but too large values result in unjustified memory
consumption. So for we adopted = 50 for the high bound. The distribution of ( )for low and
high bounds is presented in Figure 2. The effect of this distortion is maximal where ( )is small, i.e.
in the blue image areas, mainly the ceiling and the small box. Later we shall see that it is these areas
where the noise deviates from the mean value.
(right). The actual average is 5.6 and 54, respectively, because rounding to integers increases the
values as compared to (
As explained in Section 2, what determines the speed of convergence is the product (
One can see that the noise image roughly has three distinct areas: walls of the main box, ceiling and the small box, the noise being nearly homogeneous throughout each area. The corresponding values are printed below each image for quantitative comparison.
The areas where the noise considerably deviates from the mean level coincide with areas where
( )is below 3 or 5, see Figure 2. Meanwhile prediction of noise level (
The situation is qualitatively similar to the case of larger radius (Table 6). Again one sees the same
three distinct areas: walls of the main box, ceiling and the small box, the noise being nearly
homogeneous throughout each area. The corresponding values are printed below each image for
quantitative comparison. Again the areas where the noise level deviates from mean are those where
( )is below 2 or 3 (cf. Figure 2) and so rounding to integers deviated them from (
From comparison of the calculated variants one can see that 1. Noise is the more homogeneous the larger the (natural: because of the rounding ( ) becomes less essential) and the larger the ratio / (which looks strange).
The noise for the most of the area is not sensitive to the and . In other words, the optimum is “very wide” and deviating from it even threefold does not change the mean noise level much. This is good in many aspects because in real simulation of complex scenes it may be difficult to find the exact values of the number of rays (as they use coefficients that must be ray-traced themselves). Also time per iteration is difficult to measure as it is a random value and additionally it may have regular trends because of repeating caching etc. 3. The areas where the noise considerably deviates from the mean level coincide with areas where ( )is below 2 or 3, see Figure 2. Therefore it is reasonable to choose (which all the same does not have a precise optimum value) large enough so that ( )exceeds, say, 3 in the most of the image.
This paper operates the absolute noise, while sometimes it is better to work with the relative noise (divided by the luminance in this pixel). Our approach can be easily adapted to this case as well.
We should like to thank Elissey Birukov for his help with computations.
[9] I. Georgiev, J. Krivánek, T. Davidovic, P. Slusallek, Light transport simulation with vertex
connection and merging, ACM Trans. Graph., 31(