=Paper=
{{Paper
|id=Vol-2763/CPT2020_paper_s6-2
|storemode=property
|title=Direct calculation of the optimal weight for MIS
|pdfUrl=https://ceur-ws.org/Vol-2763/CPT2020_paper_s6-2.pdf
|volume=Vol-2763
|authors=Sergey Ershov,Alexey Voloboy,Dmitry Zhdanov,Vladimir Frolov,Andrey Zhdanov
}}
==Direct calculation of the optimal weight for MIS==
https://ceur-ws.org/Vol-2763/CPT2020_paper_s6-2.pdf
Direct calculation of the optimal weight for MIS
S.V. Ershov1, A.G. Voloboy1, D.D. Zhdanov2, A.D. Zhdanov2, V.A. Frolov1
voloboy@gin.keldysh.ru
1
Keldysh Institute for Applied Mathematics, Moscow, Russia
2
National Research University of Information Technologies, Mechanics, and Optics, St. Petersburg, Russia
A Monte-Carlo ray tracing is nowadays standard approach for lighting simulation and generation of realistic images. A widely
used method for noise reduction in Monte-Carlo ray tracing is combing different means of sampling, known as Multiple Importance
Sampling (MIS). For bi-directional Monte-Carlo ray tracing with photon maps (BDPM) the join paths are obtained by merging
camera and light sub-paths. Since several light paths are checked against the same camera path and vice versa, the join paths
obtained are not statistically independent. Thus the noise in this method does not obey the laws which are correct in simple classic
Monte-Carlo with independent samples. And, correspondingly, the MIS weights that minimize that noise must also be calculated
differently. In this paper we calculate these weights for a simple model scene directly minimizing the noise of calculation. This is a
pure direct numerical minimization that does not involve any doubtful hypothesis or approximations. We show that the weights
obtained are qualitatively different from those calculated from classic “balance heuristic” for Monte-Carlo with independent samples.
They depend on the scene distance, but not only on scattering properties of the surfaces and the distribution of light source emission.
Keywords: Monte-Carlo ray tracing, bi-directional ray tracing, photon maps, reduction of noise, multiple importance sampling,
weights.
1. Introduction In this paper we calculate the optimal weights
directly, i.e. find the minimum of the sample variance of
A powerful method of solution of the rendering the pixel value. This is performed for a simple model
equations is Monte Carlo ray tracing (MCRT). It is scene. We demonstrate that these weights depend on the
widely used in calculation of the global illumination [1, geometry of the scene and on the number of light and
2]. Its main problem is noise, and it strongly depends on camera rays per iteration, while the known MIS
the method of generation of random points. Therefore formulae from [3, 4] include only the BDFs and
there were and are a lot of papers devoted to the optimal distribution of light source emission.
choice of the probability distribution of ray scattering
[3–9]. One of the powerful approaches here is the so- 2. BDPM and weights in it
called Multiple Importance Sampling (MIS). Its idea is
The basic idea of BDPM is that we trace several
that we generate several random samples (rays)
camera and several light rays. Then for each pair of
according to different “strategies” i.e. probability
light + camera paths, we try to merge them in a join
distributions and then sum with weights their
trajectory that connects light and camera. If they do join
contributions to image luminance.
we increment the accumulated luminance. Then the next
The mathematics behind that was produced in the
pair is processed. After all light rays had been checked
famous thesis by E Veach [3] where the theorem was
against all camera rays they all are discarded and new
proved about several simple schemes of weight
sets of rays are generated etc. Generation of the sets of
calculation. It was proved there that the resulting noise
rays and then cycling over all pairs constitute one
is close to its minimal value. This theorem applies to the
iteration of the process. The luminance calculated in
classic MCRT method when successive random points
different iterations is statistically independent.
are absolutely independent.
So, an iteration which uses 𝑁𝑁𝐵𝐵 camera rays (through
Lighting simulation meanwhile frequently uses not
given pixel!!) and 𝑁𝑁𝐹𝐹 light rays (for all pixels!)
that simple MCRT but more advanced methods like bi-
increases accumulated luminance of a pixel by
directional Monte-Carlo path tracing (BDPT), bi- 𝑁𝑁𝐹𝐹 𝑁𝑁𝐵𝐵
(1)
directional Monte-Carlo ray tracing with photon maps 1 1
(BDPM) [2], their combination termed sometimes 𝛥𝛥𝛥𝛥 = � � 𝐶𝐶𝑖𝑖,𝑗𝑗
𝑁𝑁𝐹𝐹 𝑁𝑁𝐵𝐵
BDCM [8, 9] etc. Here the successive trajectories are 𝑖𝑖=1 𝑗𝑗=1
not quite independent, for example, in the BDPM the where 𝐶𝐶𝑖𝑖,𝑗𝑗 is the contribution from the pair of 𝑖𝑖-th light
same forward path is “merged” with all the backward and 𝑗𝑗-th camera rays. Similarly to it can be written as
paths. Therefore the resulting joined full trajectories 𝐶𝐶 = � � 𝑤𝑤
(𝑐𝑐) (𝑙𝑙)
𝐸𝐸 (𝑐𝑐) (𝑥𝑥⃗ )𝐸𝐸 (𝑙𝑙) (𝑥𝑥⃗ )𝐾𝐾(𝑥𝑥⃗
(𝑐𝑐) (2)
𝑚𝑚+𝑛𝑛,𝑚𝑚 𝑚𝑚 𝑛𝑛 𝑚𝑚
have common ”tail” and thus are not independent. 𝑚𝑚 𝑛𝑛
As a result, the noise in these methods follows other (𝑙𝑙) (𝑐𝑐)
− 𝑥𝑥⃗𝑚𝑚 )𝑓𝑓(⋯ ; 𝑥𝑥⃗𝑚𝑚 )
rules than in the simple or classic MCRT [6]. Therefore where 𝑚𝑚 cycles over all camera path vertices and 𝑛𝑛
the weights that minimize this noise are likely different cycles over all light path vertices, 𝐾𝐾 is the integration
from those which minimize the noise functional in the (𝑐𝑐)
classic MCRT. We shall prove it for the example of a kernel and 𝑓𝑓is BDF in luminance units at the point 𝑥𝑥⃗𝑚𝑚 .
very simple model scene. In this scene the noise level is Like in [3, 4], 𝑤𝑤𝑘𝑘,𝑚𝑚 is the weight for junction at the 𝑚𝑚-
dictated by geometric factors, but not by object optical th camera vertex when the join path of 𝑘𝑘 vertices (i.e.
properties in form of bi-directional distribution function the light half of the join path has 𝑘𝑘 − 𝑚𝑚 vertices). It
(BDF), while the Veach formulae [3, 4] relate weights must be a function of that full path such that
to the BDFs along the ray path.
Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY
4.0)
𝑘𝑘−1
plane 3
� 𝑤𝑤𝑘𝑘,𝑚𝑚 = 1
plane 1
plane 2
𝑚𝑚=0
Notice that in principle we have different sets of
weights for joint paths of different total length 𝑘𝑘 (this is
obvious because they are functions of 𝑘𝑘 arguments).
L2 L3
3. Direct calculation of optimal weights
It follows from (1) and (2) that the increment of the Fig. 1: The model scene
pixel luminance from one algorithm iteration is also
linear in weights: it is a sum of weights times some The spatial distribution of illumination of plane 3 is
random functions: 𝐼𝐼(𝑥𝑥, 𝑦𝑦). Since transmittance is the Lambert the angular
∞ 𝑘𝑘 distribution of incident light is irrelevant.
1
𝛥𝛥𝛥𝛥 = � � � 𝑤𝑤𝑘𝑘,𝑚𝑚 (𝜁𝜁⃗𝑖𝑖,𝑗𝑗 )𝐶𝐶𝑘𝑘,𝑚𝑚 (𝜁𝜁⃗𝑖𝑖,𝑗𝑗 )
𝑁𝑁𝐹𝐹 𝑁𝑁𝐵𝐵 Join paths and weights
𝑘𝑘=1 𝑚𝑚=0 𝑖𝑖,𝑗𝑗
where 𝑖𝑖 enumerates light rays in one iteration, 𝑗𝑗 This scene has no reflection and all paths connecting
enumerates camera rays (through this pixel) in one the camera and light source are qualitatively the same
iteration and 𝜁𝜁⃗ is the join path from them. 𝐶𝐶𝑘𝑘,𝑚𝑚 (𝜁𝜁⃗𝑖𝑖,𝑗𝑗 ) is camera → 𝑥𝑥⃗1 → 𝑥𝑥⃗2 → 𝑥𝑥⃗3 → light
the contribution from the this pair (𝑖𝑖, 𝑗𝑗) constrained to
conditions where 𝑥𝑥⃗1 = 0 �⃗ is fixed and the segment between lights
1. The paths merge into join path source and plane 3 is ignored because does not affect
2. They merge at the 𝑚𝑚-th camera vertex the path contribution. Therefore the full path is
3. The join path has 𝑘𝑘 vertices completely described by its two variable vertices, 𝑥𝑥⃗2
If these conditions are not satisfied 𝐶𝐶𝑘𝑘,𝑚𝑚 vanish. This and 𝑥𝑥⃗3 .
resolves ambiguity how to define the join path if the The camera and light rays can meet at planes 1, 2
camera and light halves did not merge. and 3 whose contributions are taken with weights
The sets of join paths from different iterations are 𝑤𝑤0 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ), 𝑤𝑤1 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) and 𝑤𝑤2 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ). Because of
independent. So for this linear form the average over normalization 𝑤𝑤0 + 𝑤𝑤1 + 𝑤𝑤2 it suffices to calculate 𝑤𝑤0
iterations (= the limiting luminance) is also linear in and 𝑤𝑤1 .
weights. The mean square of luminance value If camera and light ray meet at plane 𝑚𝑚, it is
calculated in one iteration is a quadratic form in ambiguous whether 𝑥𝑥⃗𝑚𝑚 is camera or light hit (they can
weights whose “coefficients” are averages that can be differ by integration kernel radius). We choose camera
calculated in ray tracing. hit then.
Therefore we can find those weights that minimize
the variance of the pixel luminance, i.e. the noise. These Calculation of contribution
are the optimal weights. Its direct calculation that does (𝑐𝑐) (𝑐𝑐) (𝑐𝑐)
Camera path is (𝑥𝑥⃗1 , 𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) and light path is
not include an approximations and hypothesis is (𝑙𝑙) (𝑙𝑙) (𝑙𝑙) (𝑐𝑐)
regrettably very expensive numerically. So we shall (𝑥𝑥⃗3 , 𝑥𝑥⃗2 , 𝑥𝑥⃗1 ) where 𝑥𝑥⃗𝑚𝑚 is the hit point of camera ray
(𝑙𝑙)
perform it for a simple model scene, but even this at the 𝑚𝑚-th plane, 𝑥𝑥⃗𝑚𝑚 is the hit point of light ray at the
example will give us some important conclusions. 𝑚𝑚-th plane (light ray goes from plane 3 to plane 2 then
(𝑐𝑐) �⃗ is fixed. As said above we do
to plane 1), and 𝑥𝑥⃗1 = 0
4. Simple model scene and calculations for it
not consider the light ray before it hits the plane 3; just
(𝑙𝑙)
Scene layout we start the ray by choosing the point 𝑥𝑥⃗3 at random.
The contribution of these two sub-paths is
To simplify our calculations we use a model scene (𝑙𝑙) (𝑙𝑙) (𝑐𝑐) (𝑙𝑙) (𝑐𝑐)
where all join paths have the same (and small) length. It 𝐶𝐶 = 𝑤𝑤0 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 )𝐸𝐸 (𝑐𝑐) (𝑥𝑥⃗1 )𝐸𝐸 (𝑙𝑙) (𝑥𝑥⃗1 )𝐾𝐾(𝑥𝑥⃗1
(𝑙𝑙)
consists of 3 parallel planes with diffuse transparency; − 𝑥𝑥⃗1 )𝑓𝑓1 (⋯ )
planes 1 and 3 have the Lambert BDF. BDF of the (𝑐𝑐) (𝑙𝑙) (𝑐𝑐) (𝑙𝑙) (𝑐𝑐) (𝑙𝑙)
+𝑤𝑤1 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 )𝐸𝐸 (𝑐𝑐) (𝑥𝑥⃗2 )𝐸𝐸 (𝑙𝑙) (𝑥𝑥⃗2 )𝐾𝐾(𝑥𝑥⃗2 − 𝑥𝑥⃗2 )𝑓𝑓2 (⋯ )
middle plane 2 is arbitrary and can be made very sharp (𝑐𝑐) (𝑐𝑐) (𝑐𝑐) (𝑙𝑙) (𝑐𝑐)
+ 𝑤𝑤2 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 )𝐸𝐸 (𝑐𝑐) (𝑥𝑥⃗3 )𝐸𝐸 (𝑙𝑙) (𝑥𝑥⃗3 )𝐾𝐾(𝑥𝑥⃗3
(when direction of a transmitted ray is close to that of (𝑙𝑙)
the incident ray). The planes are orthogonal to Oz and − 𝑥𝑥⃗3 )𝑓𝑓3 (⋯ )
are positioned at 𝑧𝑧 = 0, 𝑧𝑧 = 𝐿𝐿2 and 𝑧𝑧 = 𝐿𝐿2 + 𝐿𝐿3 where 𝑓𝑓1 = 𝑓𝑓2 = 𝜋𝜋 −1 while
respectively. Camera looks at plane 1 at normal 𝜗𝜗2
1 − 2 1
direction. We consider a single pixel such that the 𝑓𝑓2 = 𝑒𝑒 2𝛽𝛽
2𝜋𝜋𝛽𝛽 2 cos𝛾𝛾
camera ray hits plane 1 at (0, 0, 0) point. The rightmost where 𝜗𝜗 is the angle between the incident and scattered
plane 3 is illuminated by light source from right side, rays, 𝛾𝛾 is the angle between the scattered ray and the
see Fig. 1. normal and 𝛽𝛽 is the width. As to the integration kernel,
we use the simplest one:
1 1, |𝑥𝑥⃗| ≤ 𝑅𝑅
𝐾𝐾(𝑥𝑥⃗) = �
𝜋𝜋𝑅𝑅2 0, |𝑥𝑥⃗| > 𝑅𝑅
where 𝑅𝑅 is integration radius. It is small. 𝓒𝓒 = �𝐴𝐴�|𝑊𝑊⟩ + 𝐵𝐵
(𝑐𝑐) (𝑐𝑐) (𝑐𝑐)
Denoting the camera path as 𝜉𝜉⃗ ≡ (𝑥𝑥⃗1 , 𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) 𝓒𝓒2 = ⟨𝑊𝑊|𝐷𝐷|𝑊𝑊⟩ + 2�𝐹𝐹�|𝑊𝑊⟩ + 𝐵𝐵2
(𝑙𝑙) (𝑙𝑙) (𝑙𝑙)
and light path as 𝜂𝜂⃗ ≡ (𝑥𝑥⃗3 , 𝑥𝑥⃗2 , 𝑥𝑥⃗1 ) we can write where the overbar denotes the average, and
2
𝐷𝐷 ≡ |𝐴𝐴⟩⟨𝐴𝐴|
𝐶𝐶(𝜉𝜉⃗, 𝜂𝜂⃗) = � 𝑤𝑤𝑚𝑚 (𝑥𝑥⃗2 (𝜉𝜉⃗, 𝜂𝜂⃗), 𝑥𝑥⃗3 (𝜉𝜉⃗, 𝜂𝜂⃗))𝐶𝐶𝑚𝑚 (𝜉𝜉⃗, 𝜂𝜂⃗)
⟨𝐹𝐹| ≡ 𝐵𝐵⟨𝐴𝐴|
𝑚𝑚=0
where The mathematical expectation of pixel luminance is
(𝑐𝑐) (𝑙𝑙) (𝑐𝑐) (𝑙𝑙)
𝐶𝐶𝑚𝑚 ≡ 𝐸𝐸 (𝑐𝑐) (𝑥𝑥⃗𝑚𝑚 )𝐸𝐸 (𝑙𝑙) (𝑥𝑥⃗𝑚𝑚 )𝐾𝐾(𝑥𝑥⃗𝑚𝑚 − 𝑥𝑥⃗𝑚𝑚 )𝑓𝑓𝑚𝑚+1 (⋯ ) the same for all weights, thus the limiting average
are independent from weights. Notice that usually in �𝐴𝐴� = 0. But for a finite number of iterations, when
MCRT the ray energies can be calculated convergence is incomplete, the sample average can be
deterministically from the path. But even if they are slightly depending on weights. Thus the sample average
random this does not affect our derivation, just �𝐴𝐴� ≠ 0 and while calculating the sample variance over
𝐶𝐶𝑚𝑚 becomes “more random”. a finite number of iterations we must account for this
Pixel luminance calculated during 𝑀𝑀 iterations, each dependence. This sample variance over 𝑀𝑀 iterations is
of which uses 𝑁𝑁𝐹𝐹 light rays and 𝑁𝑁𝐵𝐵 camera rays, is thus
𝑀𝑀
1 1
𝐿𝐿 = � 𝓒𝓒𝑠𝑠 𝑉𝑉𝑉𝑉𝑉𝑉(𝐿𝐿) = �⟨𝑊𝑊|𝐷𝐷|𝑊𝑊⟩ + 2�𝐹𝐹�|𝑊𝑊⟩ + 𝐵𝐵2
𝑀𝑀 𝑀𝑀
𝑠𝑠=1 2
𝑁𝑁𝐹𝐹 𝑁𝑁𝐵𝐵 (3) − ��𝐴𝐴�|𝑊𝑊⟩ + 𝐵𝐵� �
1 1 so the weights which minimize it satisfy
𝓒𝓒𝑠𝑠 ≡ � � 𝐶𝐶(𝜉𝜉⃗𝑠𝑠,𝑗𝑗 , 𝜂𝜂⃗𝑠𝑠,𝑖𝑖 )
𝑁𝑁𝐵𝐵 𝑁𝑁𝐹𝐹 �𝐷𝐷 − �𝐴𝐴��𝐴𝐴��|𝑊𝑊⟩ = −|𝐹𝐹⟩ + ⟨𝐵𝐵⟩�𝐴𝐴�
𝑖𝑖=1 𝑗𝑗=1
The contribution from each iteration (3) is random which is just a system of simultaneous linear equations.
variable and contributions from different iteration. However numerical experiments shown that the
Therefore the variance of calculated luminance is solution can be rather ragged. To improve the situation a
1 regularization term can be added to the minimization
𝑉𝑉𝑉𝑉𝑉𝑉(𝐿𝐿) = 𝑉𝑉𝑉𝑉𝑉𝑉(𝓒𝓒) equation which is a penalty for high gradients.
𝑀𝑀
𝑉𝑉𝑉𝑉𝑉𝑉(𝓒𝓒) = ⟨𝓒𝓒2 ⟩ − ⟨𝓒𝓒⟩2
The averages are over the ray ensembles. They can 5. Results
be approximately estimated from the sum over We performed the calculations for the case when
iterations (the usual practice called sample mean and • plane positions 𝐿𝐿2 = 1, 𝐿𝐿3 = 3,
sample variance). • BDF of plane 2 has width 𝛽𝛽 = 3∘
Tabulated weights • illumination density 𝐼𝐼 is
1 𝑟𝑟𝐿𝐿𝐿𝐿 2
For numerical calculations let us subdivide the � � = 3333 𝑟𝑟3 ≤ 𝑎𝑎
whole admissible area in (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) space in cells. The 𝐼𝐼(𝑟𝑟3 ) = �300 𝑎𝑎
1 𝑎𝑎 < 𝑟𝑟3 < 𝑟𝑟𝐿𝐿𝐿𝐿
weight is constant 𝑤𝑤𝑚𝑚,𝛼𝛼 within cell. Since formally 𝑥𝑥⃗2 0 𝑟𝑟3 ≥ 𝑟𝑟𝐿𝐿𝐿𝐿
and 𝑥𝑥⃗3 can be infinite we take some finite area and where 𝑟𝑟𝐿𝐿𝐿𝐿 = 20 is the radius of illuminated area and
subdivide it in a usual way, unbounded space outside it 𝑎𝑎 = 0.02 is the “aperture” (radius) of its bright central
constituting “the last cell”. Let 𝜒𝜒𝛼𝛼 (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) be 1 inside part
the 𝛼𝛼-th cell and 0 outside it. Then the contribution • integration radius 𝑅𝑅 = 0.003
from the 𝑠𝑠-th iteration (3) becomes • the number of rays 𝑁𝑁𝐵𝐵 = 100, 𝑁𝑁𝐹𝐹 = 107
(𝑚𝑚) (𝑚𝑚) (2)
𝐴𝐴𝛼𝛼,𝑠𝑠 ≡ 𝓒𝓒𝛼𝛼,𝑠𝑠 − 𝓒𝓒𝛼𝛼,𝑠𝑠 The scene is axisymmetrical. Therefore all functions
(2) of (𝑥𝑥⃗2 , 𝑥𝑥⃗3 ) actually depend on 3, not 4 variables:
𝐵𝐵𝑠𝑠 ≡ � 𝓒𝓒𝛼𝛼,𝑠𝑠
(𝑟𝑟2 , 𝑟𝑟3 , 𝜑𝜑) where 𝜑𝜑 is the angle between vectors 𝑥𝑥⃗2 and
𝛼𝛼
𝑁𝑁𝐹𝐹 𝑁𝑁𝐵𝐵 𝑥𝑥⃗3 (notice that 𝜑𝜑 and 2𝜋𝜋 − 𝜑𝜑 give the same result!).
(𝑚𝑚) 1 1 BTW one can prove that the optimal weights are
𝓒𝓒𝛼𝛼,𝑠𝑠 ≡ � � 𝜒𝜒𝛼𝛼 (𝜉𝜉⃗𝑠𝑠,𝑗𝑗 , 𝜂𝜂⃗𝑠𝑠,𝑖𝑖 )𝐶𝐶𝑚𝑚 (𝜉𝜉⃗𝑠𝑠,𝑗𝑗 , 𝜂𝜂⃗𝑠𝑠,𝑖𝑖 )
𝑁𝑁𝐵𝐵 𝑁𝑁𝐹𝐹 independent from 𝜑𝜑.
𝑖𝑖=1 𝑗𝑗=1
As said above rays fill the area with 𝑟𝑟2 and 𝑟𝑟3 up to
Combining the tables that relate to the weights 𝑤𝑤0
infinity. So we chose a finite area for each, now 0 ≤
and 𝑤𝑤1 into single array: 1
𝑤𝑤0 𝑟𝑟2 ≤ , 0 ≤ 𝑟𝑟3 ≤ 0.05, subdivided it into equal cells and
𝑊𝑊 ≡ �𝑤𝑤 � 2
1 then added the last cell which completes to the whole
(0) 1
𝐹𝐹𝑠𝑠 infinite domain, e.g. < 𝑟𝑟2 < ∞.
𝐹𝐹𝑠𝑠 ≡ � (1) � 2
𝐹𝐹𝑠𝑠 Trial calculations were done for several numbers of
(0)
𝐴𝐴𝑠𝑠 cells. It happened that although the calculated weights
𝐴𝐴𝑠𝑠 ≡ � (1) � differ, the noise level is nearly the same (as it is
𝐴𝐴𝑠𝑠
common for optimization). Since the weights are not
we can write
needed per se, but only the noise reduction by them, we
𝓒𝓒𝑠𝑠 = ⟨𝐴𝐴𝑠𝑠 || 𝑊𝑊⟩ + 𝐵𝐵𝑠𝑠
can use as small cells as enough to saturate the noise
Then, since the values for different 𝑠𝑠 are
level.
independent the average value and the mean square are
It happened that it was enough 1 cell in 𝜑𝜑 (i.e. characteristics.
weights actually do not depend on it!), 2 cells in 𝑟𝑟2 (0 ≤ Meanwhile in the “balance heuristic” or “power
1 1
𝑟𝑟2 ≤ and < 𝑟𝑟2 < ∞) and 26 cells in 𝑟𝑟3 (the first 25 of heuristic” [3, 4] this function is known in advance. Very
2 2
roughly, it calculates the weight from the ratio of BDF
size 0.002 and the last 0.05 < 𝑟𝑟3 < ∞).
at junction point to the sum of BDFs at all the vertices
The noise was calculated for the following cases.
of the join path. We therefore conclude that the
The calculation results are shown in Table 1:
balance/power heuristic, derived for the usual MCRT, is
1. rays meet at plane 1 only
not truly optimal for BDPM because there the
2. rays meet at plane 2 only
“samples” (join paths) are correlated because use the
3. rays meet at plane 3 only
same light and/or camera path several times.
4. rays meet at plane 3 only
5. optimal weights are used Acknowledgments
Table 1. The calculation results
case 𝒘𝒘𝟎𝟎 𝒘𝒘𝟏𝟏 𝒘𝒘𝟐𝟐 𝑳𝑳×105 RMS,% The study was carried out within the framework of
1 1 0 0 26.484 208% the RFBR grants 18-01-00569, 18-31-20032 and 20-01-
2 0 1 0 26.352 254% 00547.
3 0 0 1 25.997 146%
1, 𝑟𝑟3 ≤ 𝑎𝑎 0, 𝑟𝑟3 ≤ 𝑎𝑎 References:
4 0 � � 26.352 61%
0, 𝑟𝑟3 > 𝑎𝑎 1, 𝑟𝑟3 > 𝑎𝑎
[1] Matt Pharr and Greg Humphreys. 2010. Physically Based
5 0 optimal 25.961 54%
Rendering, Second Edition: From Theory to
Implementation (2nd ed.). Morgan Kaufmann Publishers
Optimal weights were calculated from statistic Inc., San Francisco, CA, USA.
accumulated in 10000 iterations (Fig. 2). Optimization [2] H. W. Jensen, Global illumination using photon maps, in
was constrained to 𝑤𝑤0 = 0. Proceedings of the Eurographics Workshop on Rendering
Techniques '96, (London, UK, UK), pp. 21–30, Springer-
Verlag, 1996.
[3] Eric Veach. A dissertation: Robust Monte-Carlo methods
for light transport simulation, 1997.
[4] Jiri Vorba. Bidirectional photon mapping. In Proceedings
of CESCG 2011: The 15th Central European Seminar on
Computer Graphics, Prague, 2011.
[5] I. Georgiev, J. Křivánek, T. Davidovič, and Ph. Slusallek.
2012. Light transport simulation with vertex connection
and merging. ACM Trans. Graph. 31, 6, Article 192
(November 2012)
[6] S. Ershov, D. Zhdanov, and A. Voloboy. Estimation of
noise in calculation of scattering medium luminance by
MCRT. Mathematica Montisnigri, XLV: 60–73, 2019.
[7] S. Ershov, D. Zhdanov, A. Voloboy, M. Sorokin. Treating
diffuse elements as quasi-specular to reduce noise in bi-
directional ray tracing // Keldysh Institute Preprints.
Fig. 2: The optimal weight 𝑤𝑤1 as a function of 𝑟𝑟2 ; 2018. No. 122. 30 p. doi:10.20948/prepr-2018-122-e
𝑤𝑤2 = 1 − 𝑤𝑤1 ; 𝑤𝑤0 was constrained to 0 [8] S. Popov, R. Ramamoorthi, F. Durand, and G. Drettakis,
Probabilistic Connections for Bidirectional Path Tracing,
6. Conclusions Computer Graphics Forum, 2015.
We see that even in case of a direct optimization [9] N. Dodik, Implementing probabilistic connections for
bidirectional path tracing in the Mitsuba Renderer, Sept.
(which gives the best result without false minima,
2017.
approximations etc.) the gain is moderate; it is about
3fold as compared to the best “fixed BDD” strategy. About the authors
This is not bad because 3fold in noise is equivalent to a
9fold increase of speed. Sergey V. Ershov, PhD, senior researcher in the Keldysh
Institute of Applied Math of RAS. E-mail:
At qualitative level we see that the optimal weights
measure@spp.keldysh.ru
are not local i.e. we cannot calculate the weight (which Alexey G. Voloboy, Doctor of Science in physics and
is as we remember a function of the vertices of join mathematics, leading researcher in the Keldysh Institute of
path) from that path only. Indeed, in the above Applied Math of RAS. E-mail: voloboy@gin.keldysh.ru
calculation illumination of the rightmost plane was Dmitry D. Zhdanov, PhD, associate professor in the
3333 times lower for 𝑟𝑟3 > 𝑎𝑎. Let us compare it with the Information Technologies, Mechanics and Optics (ITMO)
case of uniform illumination. In this case the optimal University. E-mail: ddzhdanov@mail.ru
weight is very close to 𝑤𝑤2 = 1 (for all paths), so for a Andrey D. Zhdanov, PhD student in the Information
join path with |𝒙𝒙3 | ≤ 𝑎𝑎 the optimal weight is different Technologies, Mechanics and Optics (ITMO) University. E-
mail: adzhdanov@itmo.ru
for the uniform and not uniform illumination. In other Vladimir A. Frolov, PhD, senior researcher in the Keldysh
words, the weight for this path depends on illumination Institute of Applied Math of RAS. E-mail: vova@frolov.pp.ru
outside it.
Surely the optimal weight is still a function of the
join path but this function depends on the global scene