Rendering: Monte Carlo Integration II Bernhard Kerbl Research - - PowerPoint PPT Presentation

Rendering: Monte Carlo Integration II

Bernhard Kerbl

Research Division of Computer Graphics Institute of Visual Computing & Human-Centered Technology TU Wien, Austria

With slides based on material by Jaakko Lehtinen, used with permission

න

Integrating the cosine-weighted radiance 𝑀𝑗(𝑦, 𝜕) at a point 𝑦 Integral of the light function

• ver the hemisphere, w.r.t.

direction/solid angle at 𝜕 Let’s find a solution!

How do we integrate over the hemisphere? How do we do it smartly?

Today’s Goal

Rendering – Monte Carlo Integration I 2

Sampling a Unit Disk

Imagine we have a disk-shaped surface with radius 𝑠 = 1 that registers incoming light (color) from directional light sources As an exercise, we want to approximate the total incoming light over the disk’s surface area We integrate over an area of size 𝜌 We will use the Monte Carlo integral for that

Rendering – Monte Carlo Integration II 3

2 2 r = 1

Uniformly Sampling the Unit Disk

If we can manage to uniformly sample the disk, then we can compute the Monte Carlo integral as a simple average × 𝜌 By drawing uniform samples in 𝑦 and 𝑧, we cannot cover the area precisely Inscribed square: information lost Circumscribed square: unnecessary samples

Rendering – Monte Carlo Integration II 4

Uniformly Sampling the Unit Disk

If we can manage to uniformly sample the disk, then we can compute the Monte Carlo integral as a simple average By drawing uniform samples in 𝑦 and 𝑧, we cannot cover the area precisely Inscribed square: information lost Circumscribed square: unnecessary samples

Rendering – Monte Carlo Integration II 5

Uniformly Sampling the Unit Disk

If we can manage to uniformly sample the disk, then we can compute the Monte Carlo integral as a simple average By drawing uniform samples in 𝑦 and 𝑧, we cannot cover the area precisely Inscribed square: information lost Circumscribed square: unnecessary samples

Rendering – Monte Carlo Integration II 6

If we can manage to uniformly sample the disk, then we can compute the Monte Carlo integral as a simple average By drawing uniform samples in 𝑦 and 𝑧, we cannot cover the area precisely Inscribed square: information lost Circumscribed square: unnecessary samples

This is actually somewhat ok!

Uniformly Sampling the Unit Disk

Rendering – Monte Carlo Integration II 7

Rejection Sampling

Requires a PDF 𝑞𝑌(𝑦) and a constant 𝑑 such that 𝑔 𝑦 < 𝑑𝑞𝑌(𝑦) Draw 𝜊𝑗 and 𝑌𝑗 from their respective distributions. If the point (𝑌𝑗, 𝜊𝑗𝑑𝑞𝑌 𝑌𝑗 ) lies under 𝑔 𝑦 , then the sample is accepted

loop forever: sample 𝑌 from 𝑞𝑌’s distribution if 𝜊𝑗 ⋅ 𝑑𝑞𝑌 𝑌𝑗 < 𝑔 𝑌𝑗 then return 𝑌𝑗

Rendering – Monte Carlo Integration II 8

Rejection Sampling

Requires a PDF 𝑞𝑌(𝑦) and a constant 𝑑 such that 𝑔 𝑦 < 𝑑𝑞𝑌(𝑦) Draw 𝜊𝑗 and 𝑌𝑗 from their respective distributions. If the point (𝑌𝑗, 𝜊𝑗𝑑𝑞𝑌 𝑌𝑗 ) lies under 𝑔 𝑦 , then the sample is accepted

loop forever: sample 𝑌 from 𝑞𝑌’s distribution if 𝜊𝑗 ⋅ 𝑑𝑞𝑌 𝑌𝑗 < 𝑔 𝑌𝑗 then return 𝑌𝑗

Unit disk: 𝑔 𝑦, 𝑧 = ൝1 𝑗𝑔 ( 𝑦2 + 𝑧2 ≤ 1) 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 , 𝑞 𝑦, 𝑧 = 1

4 , 𝑑 = 4

Rendering – Monte Carlo Integration II 9

Back to the Unit Disk

We do not want to waste samples if we can avoid it Instead, find a way to generate uniform samples on the disk Second attempt: draw from 2D polar coordinates

Polar coordinates defined by radius 𝑠 ∈ [0,1) and angle 𝜄 ∈ [0,2𝜌) Transformation to cartesian coordinates: 𝑦 = 𝑠 sin 𝜄 y = 𝑠 cos 𝜄

Rendering – Monte Carlo Integration II 10

Uniformly Sampling the Unit Disk?

Convert two 𝜊 to ranges 0, 1 , [0,2𝜌) for polar coordinates Convert to cartesian coordinates

Rendering – Monte Carlo Integration II 11

void sampleUnitDisk() { std::default_random_engine r_rand_eng(0xdecaf); std::default_random_engine theta_rand_eng(0xcaffe); std::uniform_real_distribution<double> uniform_dist(0.0, 1.0); for (int i = 0; i < NUM_SAMPLES; i++) { auto r = uniform_dist(r_rand_eng); auto theta = uniform_dist(theta_rand_eng) * 2 * M_PI; auto x = r * sin(theta); auto y = r * cos(theta); samples2D[i] = std::make_pair(x, y); } }

We successfully sampled the unit disk in the proper range However, the distribution is not uniform with respect to the area Samples clump together at center Averaging those samples will give us a skewed result for the integral!

Clumping

• 1
• 0,5

0,5 1

• 1
• 0,5

0,5 1 Rendering – Monte Carlo Integration II 12

Uniformly Sampling the Unit Disk: A Solution

The area of a disk is proportional to 𝑠2, times a constant factor 𝜌 If we see the disk as concentric rings of width Δ𝑠, the 𝑘 inner rings up to radius 𝑠

𝑘 = 𝑘Δ𝑠 should contain 𝑠𝑘 𝑠 2

𝑂 out of 𝑂 total samples Conversely, the 𝑗𝑢ℎ sample should lie in the ring at radius 𝑠𝑗 = 𝑠

𝑗 𝑂

Since 𝜊 is uniform in [0, 1), we can switch 𝑘

𝑂 for 𝜊 to get 𝑠𝑗 = 𝑠 𝜊𝑗

Rendering – Monte Carlo Integration II 13

Uniformly Sampling the Unit Disk: A Solution

It works, and it is not even a bad way to arrive at the correct solution However, for more complex scenarios, we might struggle to find the solution so easily With the tools we introduced earlier, we can formalize this process for arbitrary setups

Rendering – Monte Carlo Integration II 14

• 1
• 0,5

0,5 1

• 1
• 0,5

0,5 1

• 1
• 0,5

0,5 1

• 1
• 0,5

0,5 1

Let’s transform a regular grid from polar to cartesian coordinates

Polar To Cartesian Coordinates

Rendering – Monte Carlo Integration II 15

𝑌 = 𝑠 cos(𝜄), 𝑍 = 𝑠 sin(𝜄) 𝑠 𝜄

Take 100k samples, transform and see which box they end up in

First Attempt to Learn the PDF

Rendering – Monte Carlo Integration II 16

𝑌 = 𝜊1 cos(2𝜌𝜊2), 𝑍 = 𝜊1 sin(2𝜌𝜊2) ξ1, ξ2

Knowing the PDF

If we know the effect of a transformation 𝑈 on the PDF, we can

Use it in the Monte Carlo integral to weight our samples, or Compensate to get a uniform sampling method after transformation

Rendering – Monte Carlo Integration II 17

𝐽𝑜𝑞𝑣𝑢 (𝜊1, 𝜊2) 𝐷𝑏𝑠𝑢𝑓𝑡𝑗𝑏𝑜 (𝑦, 𝑧) 𝑄𝑝𝑚𝑏𝑠 (r, 𝜄)

𝑈(𝑠, 𝜄)

Computing the PDF after a Transformation

Assume a random variable 𝐵 and a bijective transformation 𝑈 that yields another variable 𝐶 = 𝑈 𝐵 Bijectivity dictates that 𝑐 = 𝑈(𝑏) must be either monotonically increasing or decreasing with 𝑏 This implies that there is a unique 𝐶i for every 𝐵i, and vice versa In this case, the CDFs for the two variables fulfill 𝑄𝐶 𝑈(𝑏) = 𝑄

𝐵(𝑏)

Rendering – Monte Carlo Integration II 18

Computing the PDF after a Transformation

If 𝑐 = 𝑈(𝑏) and 𝑐 increases with 𝑏, we have:

𝑒𝑄𝐶(𝑐) 𝑒𝑏

= 𝑒𝑄𝐵(𝑏)

𝑒𝑏

If 𝑐 decreases with 𝑏 (e.g. 𝑐 = −𝑏), we have: −

𝑒𝑄𝐶(𝑐) 𝑒𝑏

= 𝑒𝑄𝐵(𝑏)

𝑒𝑏

Since 𝑞𝐶 is the non-negative derivative of 𝑄𝐶, we can rewrite as: 𝑞𝐶 𝑐 𝑒𝑐 𝑒𝑏 = 𝑞𝐵 𝑏 , 𝑞𝐶 𝑐 = 𝑒𝑐 𝑒𝑏

−1

𝑞𝐵 𝑏

Rendering – Monte Carlo Integration II 19

𝑉𝑡𝑗𝑜𝑕: 𝑒𝑄

𝑌 𝑦

𝑒𝑧 = 𝑞𝑌 𝑦 𝑒𝑦 𝑒𝑧

Computing the PDF after a Transformation

Let’s interpret 𝑞𝐶 𝑐 =

𝑒𝑐 𝑒𝑏 −1

𝑞𝐵 𝑏 It is the probability density of 𝐵, multiplied by

𝑒𝑐 𝑒𝑏 −1 𝑒𝑐 𝑒𝑏 −1

has two intuitive interpretations: the change in sampling density at point 𝑏 if we transform 𝑏 by 𝑈

• r,

the inverse change in the volume of an infinitesimal hypercube at point 𝑏 if we transform 𝑏 by 𝑈

Rendering – Monte Carlo Integration II 20

Multidimensional Transformations

If we try to apply the above to the unit disk, we fail at 𝑦 = 𝑠 sin 𝜄 We can’t evaluate

𝑒𝑦 𝑒𝑠 −1

: the transformation that produces one target variable is dependent on both input variables and vice-versa We cannot compute the change in the PDF between individual variables, we must take them all into account simultaneously It’s matrix time

Rendering – Monte Carlo Integration II 21

Multidimensional Transformations

We write the set of 𝑂 values from a multidimensional variable Ԧ 𝐵 as a vector Ԧ 𝑏 and the 𝑂 outputs of transformation 𝑈 as a vector 𝑐: Ԧ 𝑏 = 𝑏1 ⋮ 𝑏𝑂 , 𝑐 = 𝑐1 ⋮ 𝑐𝑂 = 𝑈

1( Ԧ

𝑏) ⋮ 𝑈𝑂( Ԧ 𝑏) = 𝑈( Ԧ 𝑏) Instead of quantifying the change in volume incurred by 𝑈 𝑏 ,

𝑒𝑈 𝑏 𝑒𝑏

, our goal is now to quantify the change incurred by 𝑈 Ԧ 𝑏

Rendering – Monte Carlo Integration II 22

The Jacobian Matrix

For a transformation 𝑐 = 𝑈( Ԧ 𝑏), we can define the Jacobian matrix that contains all 𝑐

𝑘, 𝑏𝑗 combinations of partial differentials

𝐾𝑈( Ԧ 𝑏) =

𝜖𝑐1 𝜖𝑏1

⋯

𝜖𝑐1 𝜖𝑏𝑂

⋮ ⋱ ⋮

𝜖𝑐𝑁 𝜖𝑏1

⋯

𝜖𝑐𝑁 𝜖𝑏𝑂

If we consider Ԧ 𝐵’s domain as a space with 𝑂 axes, 𝐾𝑈( Ԧ 𝑏) gives the change of the edges of an infinitesimal hypercube from Ԧ 𝑏 to 𝑈( Ԧ 𝑏)

Rendering – Monte Carlo Integration II 23

The Jacobian Matrix, Visualized

Change of edges of an infinitesimal hypercube with extent 1,1

𝐾𝑈( Ԧ 𝑏) = 𝜖𝑐1 𝜖𝑏1 ⋯ 𝜖𝑐1 𝜖𝑏𝑂 ⋮ ⋱ ⋮ 𝜖𝑐𝑂 𝜖𝑏1 ⋯ 𝜖𝑐𝑂 𝜖𝑏𝑂

Rendering – Monte Carlo Integration II 24

𝑐 Ԧ 𝑏

1 𝜖𝑐1 𝜖𝑏1 𝜖𝑐2 𝜖𝑏1 1 𝜖𝑐1 𝜖𝑏2 𝜖𝑐2 𝜖𝑏2

The Jacobian

The columns of a square matrix 𝑁 can be interpreted as the natural base vectors of a space 1 ⋮ , 1 ⋮ if they were transformed by 𝑁 The determinant 𝑁 of 𝑁 computes the volume

• f a parallelepiped spanned by these vectors[3]

𝐾𝑈 , called the Jacobian of 𝑈, gives the volume change at Ԧ 𝑏 by 𝑈

Rendering – Monte Carlo Integration II 25

𝑐

𝑲𝑼 𝒃

Computing the PDF of a Transformation

Let’s try polar coordinates again: 𝑦 𝑧 = 𝑈 𝑠 𝜄 = 𝑠 sin 𝜄 𝑠 cos 𝜄

𝜖𝑈 𝑠 𝜄 𝜖 𝑠 𝜄

= 𝐾𝑈 =

𝜖𝑦 𝜖𝑠 𝜖𝑦 𝜖𝜄 𝜖𝑧 𝜖𝑠 𝜖𝑧 𝜖𝜄

= cos 𝜄 −𝑠 sin 𝜄 sin 𝜄 𝑠 cos 𝜄 = r 𝑞 𝑦, 𝑧 =

𝑞(𝑠,𝜄) 𝑠

, or 𝑞 𝑠, 𝜄 = 𝑠 𝑞 𝑦, 𝑧 : a uniform density in 𝑦, 𝑧 must be proportional to 𝑠 in (𝑠, 𝜄)

Rendering – Monte Carlo Integration II 26

Compare PDFs After Transformation

Rendering – Monte Carlo Integration II 27

Measured Using 𝑞 𝑠, 𝜄 =

1 2𝜌 and 𝑞 𝑦, 𝑧 = 𝑞(𝑠,𝜄) 𝑠

How To Visualize the Computed PDF Conversion?

Python code, for the interested

Rendering – Monte Carlo Integration II 28

num_bins = 20 xlist = np.linspace(-1.0, 1.0, num_bins) ylist = np.linspace(-1.0, 1.0, num_bins) bin_volume = 4 / (num_bins*num_bins) X, Y = np.meshgrid(xlist, ylist) X += 1/num_bins Y += 1/num_bins r = np.sqrt(X*X+Y*Y) uniform_dist = 1/(2 * np.pi) pdf_transform = 1/r Z = uniform_dist * pdf_transform * bin_volume Z[r>1] = 0 fig,ax = plt.subplots(1,1) cp = plt.pcolor(X, Y, Z,cmap='viridis',norm=matplotlib.colors.LogNorm()) fig.colorbar(cp) # Add a colorbar to a plot plt.show()

Sampling Joint PDFs Correctly

For independent variables, the joint PDF 𝑞(𝑦, 𝑧, … ) is 𝑞𝑌 𝑦 𝑞𝑍 𝑧 … In general, this is an assumption that we should not rely on Furthermore, after a transformation, only the joint PDF is known The proper way to sample multiple variables 𝑌, 𝑍 is to compute

the marginal density function 𝑞𝑌 𝑦 of one the conditional density function 𝑞𝑍 𝑧|𝑦 of the other

Rendering – Monte Carlo Integration II 29

Marginal and Conditional Density Function

Assume we have obtained the joint PDF 𝑞 𝑦, 𝑧 of variables 𝑌, 𝑍 with ranges [𝑏𝑌, 𝑐𝑌) and [𝑏𝑍, 𝑐𝑍) In a 2D domain with 𝑌, 𝑍 we can think of 𝑞𝑌 𝑦 as the average density of 𝑞 𝑦, 𝑧 at a given 𝑦 over all possible values 𝑧 We can obtain the marginal density function for one of them by integrating out all the others, e.g.: 𝑞𝑌 𝑦 = ׬

𝑏𝑍 𝑐𝑍 𝑞 𝑦, 𝑧 𝑒𝑧

We can then find 𝑞 𝑧|𝑦 =

𝑞(𝑦,𝑧) 𝑞𝑌(𝑦)

Rendering – Monte Carlo Integration II 30

Adding More Variables

What to do for multiple variables, e.g. 𝑌, 𝑍 and 𝑎?

Find first marginal density 𝑞𝑌 𝑦 = ׬

𝑏𝑎 𝑐𝑎 ׬ 𝑏𝑍 𝑐𝑍 𝑞 𝑦, 𝑧, 𝑨 𝑒𝑧 𝑒𝑨

Find first conditional density 𝑞𝑌 𝑧, 𝑨|𝑦 = 𝑞 𝑦,𝑧,𝑨

𝑞𝑌 𝑦

Find second marginal density 𝑞𝑍 𝑧|𝑦 = ׬

𝑏𝑎 𝑐𝑎 𝑞 𝑦, 𝑧, 𝑨 𝑒𝑨

Find second conditional density 𝑞𝑌 𝑨|𝑦, 𝑧 = 𝑞 𝑧,𝑨|𝑦

𝑞𝑍 𝑧|𝑦

Integrate + invert first marginal, first and second conditional densities Sample each of them Extend ad libitum to even more variables

Rendering – Monte Carlo Integration II 31

Sampling the Unit Disk: The Formal Solution

The size of the sampling domain in cartesian coordinates is 𝜌 Since we want uniform sampling and sample probabilities must integrate to 1, the PDF in cartesian coordinates is 𝑞 𝑦, 𝑧 =

1 𝜌

We know that 𝑞 𝑠, 𝜄 = 𝑠 𝑞 𝑦, 𝑧 , so we want 𝑞 𝑠, 𝜄 =

𝑠 𝜌

𝑞𝑆 𝑠 = ׬

2𝜌 𝑞 𝑠, 𝜄 𝑒𝜄 = 2𝑠 and 𝑞 𝜄|𝑠 = 𝑞(𝑠,𝜄) 𝑞𝑆(𝑠) = 1 2𝜌

Rendering – Monte Carlo Integration II 32

Sampling the Unit Disk: The Formal Solution

If we draw samples for our 𝑠, 𝜄 with the above PDFs, we get uniform distribution in (𝑦, 𝑧) after applying transformation 𝑈

Rendering – Monte Carlo Integration II 33

𝑠 𝜄

Sampling the Unit Disk: The Formal Solution

Integrate marginal and conditional PDFs and invert: we get the same solution as before:

𝑠 = P

𝑆 −1 𝜊1 =

𝜊1 𝜄 = 𝑄Θ

−1 𝜊2 = 2𝜌𝜊2

𝑞 𝜄|𝑠 is constant: no matter what radius we are looking at, all positions on a ring of that radius (angle) should be equally likely Final integral: 𝑆𝐻𝐶𝑢𝑝𝑢𝑏𝑚 =

𝜌 𝑂 σ𝑗=1 𝑂

𝑆𝐻𝐶(𝑆𝑗 𝑡𝑗𝑜 Θ𝑗, 𝑆𝑗 𝑑𝑝𝑡 Θ𝑗)

Rendering – Monte Carlo Integration II 34

Moving on to the Hemisphere

This took as a while, but we have seen all the formal procedures We only need to switch from integrating planar area to points 𝜕

• n hemisphere surface (i.e., vectors 𝑦, 𝑧, 𝑨 with length 1)

Use spherical coordinates and bijective 𝑈 from (𝑠, 𝜄, 𝜚) to 𝑦, 𝑧, 𝑨 : 𝑦 = 𝑠 sin 𝜄 cos 𝜚 𝑧 = 𝑠 sin 𝜄 sin 𝜚 𝑨 = 𝑠 cos 𝜄

Rendering – Monte Carlo Integration II 35

Deriving Integration Over Hemisphere

Each direction 𝜕 represents an infinitesimal surface area portion 𝑒𝜕 How do we integrate a function 𝑔(𝜕) with differential 𝑒𝜕? Integration over points on hemisphere surface 𝜕, w.r.t. (𝜄, 𝜚)

Rendering – Monte Carlo Integration II 36

𝜄 𝜚

𝜕 𝑒𝜕

𝑜

Deriving Integration Over Hemisphere

We assume a planar surface with an upright facing normal 𝑜 We use the integral intervals 𝜄 ∈ 0,

𝜌 2 , 𝜚 ∈ [0, 2𝜌)

I.e., a curve from perpendicular to parallel for 𝜄, a ring for 𝜚

Rendering – Monte Carlo Integration II 37

2𝜌 𝜌 2

𝑜

Deriving Integration Over Hemisphere

We can split the surface along 𝜄 into ribbons of width Δ𝜄 → 𝑒𝜄 The upper edge of the ribbon is slightly shorter than the lower If we keep adding more and more ribbons, this difference vanishes

Rendering – Monte Carlo Integration II 38

Δ𝜄

Deriving Integration Over Hemisphere

As a ribbon’s width goes to 𝑒𝜄, its area becomes its length times 𝑒𝜄 We can find this length by projecting the ribbon to the ground Using trigonometry, we find the length of a ribbon is 2𝜌sin 𝜄

Rendering – Monte Carlo Integration II 39

Deriving Integration Over Hemisphere

As a ribbon’s width goes to 𝑒𝜄, its area becomes its length times 𝑒𝜄 We can find this length by projecting the ribbon to the ground Using trigonometry, we find the length of a ribbon is 2𝜌sin 𝜄

Rendering – Monte Carlo Integration II 40

sin 𝜄

𝜄

𝑜

The length of a ribbon spans the entire interval 𝜚 ∈ [0, 2𝜌) Convert the length to an integral over 𝑒𝜚: 2𝜌sin 𝜄 = ׬

2𝜌 sin 𝜄 𝑒𝜚

The final integral: ׬

Ω 𝑔 𝜕 𝑒𝜕 = ׬

𝜌 2 ׬

2𝜌 𝑔(𝜕) sin 𝜄 𝑒𝜚 𝑒𝜄

Deriving Integration Over Hemisphere

Rendering – Monte Carlo Integration II 41

𝜄 𝜚

𝜕 𝑒𝜕

𝑜

Deriving PDF for Hemisphere Sampling

Integral of 𝑔(𝜕) over area Δ𝜕 = ׬

Δ𝜕 𝑔(𝜕) 𝑒𝜕

Integral of 𝑔(𝜕) w.r.t. (𝑒𝜄, 𝑒𝜚) = ׬

Δ𝜄 ׬ Δ𝜚 𝑔 𝜕 sin 𝜄 𝑒𝜚 𝑒𝜄

Integration domain and 𝑔 𝜕 are identical, thus: 𝑒𝜕 = sin 𝜄 𝑒𝜚 𝑒𝜄 𝜕 → 𝜄, 𝜚 is bijective, we have 𝑞 𝜄, 𝜚 d𝜄 𝑒𝜚 = 𝑞 𝜕 𝑒𝜕 and: 𝑞 𝜄, 𝜚 = sin 𝜄 𝑞 𝜕

Rendering – Monte Carlo Integration II 42

Δ𝜕 Δ𝜄

Deriving PDF for Hemisphere Sampling, Formal

Target distribution in 𝜕, which is (𝑦, 𝑧, 𝑨) with x2 + 𝑧2 + 𝑨2 = 1 The transformation 𝑈 from (𝑠, 𝜄, 𝜚) to 𝑦, 𝑧, 𝑨 : 𝑦 = 𝑠 sin 𝜄 cos 𝜚 𝑧 = 𝑠 sin 𝜄 sin 𝜚 𝑨 = 𝑠 cos 𝜄 The Jacobian of the transformation 𝑈 gives 𝐾𝑈 = 𝑠2 sin 𝜄 Hence, we have 𝑞 𝑠, 𝜄, 𝜚 = 𝑠2 sin 𝜄 𝑞(𝑦, 𝑧, 𝑨) = 1 sin 𝜄 𝑞(𝜕)

Rendering – Monte Carlo Integration II 43

Uniformly Sampling the Unit Hemisphere

The domain, i.e., the unit hemisphere surface area, is 2𝜌. Uniformly sampling the domain over 𝜕 implies 𝑞 𝜕 =

1 2𝜌

Hence, since 𝑞 𝜄, 𝜚 = 𝑠2 sin 𝜄 𝑞 𝜕 , we want 𝑞 𝜄, 𝜚 =

sin 𝜄 2𝜌

Marginal density 𝑞Θ(𝜄): ׬

2𝜌 𝑞 𝜄, 𝜚 𝑒𝜚 = sin 𝜄

Conditional density 𝑞 𝜚|𝜄 : 𝑞 𝜄,𝜚

𝑞Θ(𝜄) = 1 2𝜌

Rendering – Monte Carlo Integration II 44

Uniformly Sampling the Unit Hemisphere – Complete

Antiderivative of 𝑞Θ(𝜄): ׬ sin 𝜄 𝑒𝜄 = 1 − cos 𝜄 (added constant 1) Antiderivative of 𝑞 𝜚|𝜄 : ׬

1 2𝜌 𝑒𝜚 = 𝜚 2𝜌

Invert them to get 𝜄 = cos−1 𝜊1 (cos is symmetric), 𝜚 = 2𝜌𝜊2 Apply transformation 𝑈 on (𝜄, 𝜚) to obtain uniformly distributed 𝜕 Done!

Rendering – Monte Carlo Integration II 45

Arbitrary Hemispheres

The orientation of sampled hemispheres depends on surface normal Use the tangent, bitangent and normal vectors of the intersection 𝜕𝑥𝑝𝑠𝑚𝑒 = 𝑦 ⋅ 𝑢𝑏𝑜𝑕𝑓𝑜𝑢 + 𝑧 ⋅ 𝑐𝑗𝑢𝑏𝑜𝑕𝑓𝑜𝑢 + 𝑨 ⋅ 𝑜𝑝𝑠𝑛𝑏𝑚 Or, if available, use ToWorld transformation methods

Rendering – Monte Carlo Integration II 46

Applications for Uniform Hemisphere Sampling

Diffuse lighting based on last lecture’s insights with constant 𝑔

𝑠

What do we use for 𝑀𝑗? Full rendering equation: next time This time: ambient occlusion, direct lighting

Rendering – Monte Carlo Integration II 47

Consider all unblocked directions around 𝑦 as indirect light sources, integrate 𝑊 𝑦, 𝑦 + 𝛽𝜕

cos 𝜄 𝜌

• ver directions 𝜕 around the normal

Limit ray length to 𝛽, return free if no intersection closer than 𝛽

Ambient Occlusion

Rendering – Monte Carlo Integration II 48

Fine geometric details on objects are accentuated by the absence of ambient light due to the shadows cast by close-by geometry Integrate over directions 𝜕 on the unit hemisphere defined by point 𝑦, normal 𝑜 𝑊 𝑦, 𝑦 + 𝛽𝜕 = ቊ1 𝑗𝑔 𝑦 → 𝑦 + 𝛽𝜕 𝑔𝑠𝑓𝑓 0 𝑗𝑔 𝑦 → (𝑦 + 𝛽𝜕) 𝑐𝑚𝑝𝑑𝑙𝑓𝑒

Ambient Occlusion

Rendering – Monte Carlo Integration II 49

Cosine-Weighted Importance Sampling

Recap: variance is low if the sampling function mimics the signal We use 𝑔

𝑠 = 1 𝜌 for ambient occlusion, therefore the contribution of

signal samples varies mostly with cos 𝜄 It would be best to apply importance sampling: use a sampling strategy for 𝜕, such that 𝑞 𝜕 ∝ cos 𝜄 We have gone through all the necessary steps. Try to solve this formally with the inversion method as an exercise!

Rendering – Monte Carlo Integration II 50

Smart Cosine-Weighted Importance Sampling

Malley’s method: uniformly pick 𝑦, 𝑧 samples on the unit disk Project them to the hemisphere surface 𝑨 = 1 − 𝑦2 − 𝑧2 Done! Your samples are now distributed with 𝑞 𝜕 ∝ cos 𝜄 Why does this work? Try to come up with your own proof!

Rendering – Monte Carlo Integration II 51

We can use Monte Carlo integration to compute the direct illumination from light sources in the scene at a point 𝑦 Naïve version: sample unit hemisphere uniformly, hoping to hit light sources Check closest hit for each direction 𝜕 Use 𝑀𝑗 𝑦, 𝜕 = ൝𝑀𝑓

𝑚

𝑗𝑔 ℎ𝑗𝑢 𝑏 𝑚𝑗𝑕ℎ𝑢 𝑚 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓

𝑦 𝑜

Finding Light

Rendering – Monte Carlo Integration II 52

𝑚

Sampling a Light Source, Revisited

A special kind of importance sampling: integrate over light sources! Use 𝑊 𝑦, 𝑧 = ቊ1 𝑗𝑔 𝑞𝑏𝑢ℎ 𝑔𝑠𝑝𝑛 𝑦 𝑢𝑝 𝑧 𝑗𝑡 𝑣𝑜𝑐𝑚𝑝𝑑𝑙𝑓𝑒 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 Pick 𝑧 on light, e.g. uniformly

cos 𝜄𝑧 𝑒𝐵𝑧 𝑠2

: change in volume of infinitesimal 2D hypercube at 𝑧 projected onto 𝑦’s hemisphere

Rendering – Monte Carlo Integration II 53

Sampling a Light Source, Revisited

For inclusion of simple material color and BRDF, use 𝑔

𝑠 = 𝜍 𝜌

Works extremely well if the light occupies only a small solid angle

Rendering – Monte Carlo Integration II 54 100 samples per pixel, hemisphere sampling 100 samples per pixel, light source sampling

Hemisphere Sampling Illustrated

Rendering – Monte Carlo Integration II 55

Light Source Sampling Illustrated

Rendering – Monte Carlo Integration II 56

Monte Carlo Integration

But where does the actual integration step happen? In the basic case, directly in the main sampling loop for each pixel!

Static scene: Samples through pixels 𝑞 always hit the same point 𝑦 Once 𝑦 has been hit, the sampling of its hemisphere follows PDF Return sampled values (colors), weighted by the corresponding PDF

Use 𝑂 samples for 𝑞, sum color values weighted by PDF, average:

Rendering – Monte Carlo Integration II 57

𝐺𝑂 = 1 𝑂 ෍

𝑗=1 𝑂 𝑀𝑗(𝑦, 𝜕)

𝑞(𝜕)

Monte Carlo Integration as Loop Over Pixel Samples

We achieve integration around 𝑦 with multiple samples through 𝑞

A bit wasteful, but is a general, valid solution We will see in a second why this is convenient

Weight returned values by PDF, sum up and divide by 𝑂

Rendering – Monte Carlo Integration II 58

𝑦 𝑞

Monte Carlo Integration as Loop Over Pixel Samples

Given: camera, pixel p, scene, pdf rgb = {0,0,0} for i in [0, N) do ray = rayThroughPixel(camera, p) x = findIntersection(ray, scene) color = getIntegratorValue(x) rgb += color end for rgb /= N function getIntegratorValue(x)

• mega = getDirectionOnHemisphere(x, pdf)

Li = evaluateLight(x, omega, scene) return Li / pdf(omega)

Rendering – Monte Carlo Integration II 59

1 𝑂 ෍

𝑗=1 𝑂 𝑀𝑗(𝑦, 𝜕)

𝑞(𝜕)

Stratified Sampling

With random uniform sampling, we can get unlucky

E.g. all samples clump in a corner If we don’t know anything of the integrand, we want a relatively uniform sampling

Not regular, though, because of alias patterns!

Subdivide domain into non-overlapping regions (e.g. a regular grid). Each region is called a stratum Pick new random sample in stratum with lowest number of samples

Rendering – Monte Carlo Integration II 60

Monte Carlo Integration for Pixels

In the first lecture, we used supersampling to fight off aliasing Pixels are another instance where we use Monte Carlo integration Choosing samples within pixels is an instance of stratified sampling Uniform 2D distribution, average over a pixel rectangle: box filter

We will see more advanced methods for filtering in future lectures If we didn’t use at least one sample per pixel, we would leave holes

Rendering – Monte Carlo Integration II 61

Monte Carlo Integration for Pixels

Samples are randomly jittered in each stratum Ergo, we don’t hit the same 𝑦 with each pixel sample (𝑞𝑦, 𝑞𝑧) inside pixel 𝑞 We just add two random variables! Instead of integrating over a hemisphere, we are integrating over all surface points visible to a pixel and their respective hemispheres The box filter over pixel color samples implements uniform integral

Rendering – Monte Carlo Integration II 62

න

𝑞𝑧

න

𝑞𝑦

𝑀𝑓(𝑦, 𝑤) 𝑒𝑞𝑦 𝑒𝑞𝑧 𝑤1 𝑤2 𝑦1 𝑦2 𝑞

The Box Filter for a Pixel (Monte Carlo Integration)

Given: camera, pixel coordinates (pixel_x, pixel_y), scene rgb = {0,0,0} for i in [0, N) do px = pixel_x + uniform_random_sample() // 𝜊𝑗 = 𝑄𝜊 𝜊𝑗 = 𝑄𝜊

−1(𝜊𝑗)

py = pixel_y + uniform_random_sample() ray = rayThroughPixelPos(camera, px, py)

x = findIntersection(ray, scene)

rgb += getIntegratorValue(x) end for rgb /= N

Rendering – Monte Carlo Integration II 63

Nothing else necessary: offsets in x and y are uniformly distributed, they are independent, the domain size of each pixel is 1x1=1 and so we don’t have to change the integration at all.

Low-Discrepancy Series

Stratified sampling is a special case of low-discrepancy sampling[4] Replace the built-in RNG with a sample generation algorithm that sacrifices randomness for good spatial distribution Great for:

Faster convergence (reduction of noise) GPUs (they don’t love random numbers) Transparent and portable sampling behavior

Rendering – Monte Carlo Integration II 64

Halton Sequences

For 𝑜-D, pick 𝑜 different, co-prime bases (e.g. 2,3) Based on radical inverse: an integer 𝑏 can be written with base 𝑐 and 𝑒 ∈ [0, 𝑐 − 1] as 𝑏 = σ𝑗=1

𝑛 𝑒𝑗(𝑏)𝑐𝑗−1

For 𝑐 = 2, 𝑒 ∈ [0,1]: this is binary representation of integers: 13 = 0 × 20 + 1 × 21 + 1 × 22 + 1 × 23 Radical inverse with 𝑛 digits: Φ 𝑏 = σ𝑗=1

𝑛 𝑒𝑗(𝑏) 𝑐𝑗

Rendering – Monte Carlo Integration II 65

Default RNG 2,3 Halton Sequence

Jheald, Wikipedia, Halton Sequence – “Own work”

Halton Sequences

Start at any integer value for each of the 𝑜 variables

Need data structure to represent 𝑏 = σ𝑗=1

𝑛 𝑒𝑗(𝑏)𝑐𝑗−1

Can be written for arbitrary 𝑐 with some bit fiddling

For each new 𝑜-D sample, increment all 𝑜 integers Compute their radical inverse with their proper base Using same sequence for all pixels → patterns

Rendering – Monte Carlo Integration II 66

Default RNG 2,3 Halton Sequence

Jheald, Wikipedia, Halton Sequence – “Own work”

Halton Sequences

Large primes will behave similarly → patterns Repetitive patterns should always be avoided! Option 1: use different starting values for different instances (e.g. for each pixel) Option 2: instead of incrementing 𝑒𝑗, cycle through random permutations of [0, 𝑐 − 1] in each instance

E.g.: 𝑐 = 3, 0, 2, 1 instead of 0, 1, 2

Rendering – Monte Carlo Integration II 67

29,31 Halton Sequence 29,31 Halton Scrambled

Moving Forward

Get comfortable with all approaches to integration and sampling

The mental image of “area under the curve” eventually collapses (infinite-dimensional integral coming up next!) Transformations between sample domains may be non-trivial Once you grasp the underlying concepts, applying the math is easy

We have seen simpler explanations for the most important parts

Uniformly sampling a hemisphere Cosine-weighted sampling of a hemisphere

Rendering – Monte Carlo Integration II 68

References and Further Reading

Slide set based mostly on chapter 13 of Physically Based Rendering: From Theory to Implementation [1] Steven Strogatz, Infinite Powers: How Calculus Reveals the Secrets of the Universe [2] Video, Why “probability of 0” does not mean “impossible” | Probabilities of probabilities, part 2: https://www.youtube.com/watch?v=ZA4JkHKZM50 [3] Video, The determinant | Essence of linear algebra, chapter 6: https://www.youtube.com/watch?v=Ip3X9LOh2dk [4] SIGGRAPH 2012 Course: Advanced (Quasi-) Monte Carlo Methods for Image Synthesis, https://sites.google.com/site/qmcrendering/ [5] Wikipedia, Van der Corput Sequence, https://en.wikipedia.org/wiki/Van_der_Corput_sequence

Rendering – Monte Carlo Integration II 69