Monte Carlo Integration Application in Rendering

We all studied numerical methods in the course of mathematics. These are methods such as integration, interpolation, series, and so on. There are two types of numerical methods: deterministic and randomized.

Typical deterministic function integration function $$$f$in the range $$$[a, b]$It looks like this: we take $$$n + 1$evenly spaced points $$$t_0 = a, t_1 = a + \ frac {b - a} {n}, \ ldots, t_n - b$calculate $$$f$at midpoint $$$\ frac {t_i + t_ {i + 1}} {2}$of each of the intervals defined by these points, summarize the results and multiply by the width of each interval $$$\ frac {b -a} {b}$. For sufficiently continuous functions $$$f$with increasing $$$n$the result will converge to the correct value.


The probabilistic method, or the Monte Carlo method for calculating, or, more precisely, an approximate estimate of the integral $$$f$in the range $$$[a, b]$looks like this: let $$$X_1, \ ldots, X_n$- randomly selected points in the interval $$$[a, b]$. Then $$$Y = (b - a) \ frac {1} {n} \ sum_ {i = 1} ^ {n} f (X_i)$Is a random value whose average is an integral $$$\ int _ {[a, b]} f$. To implement the method, we use a random number generator that generates $$$n$points in the interval $$$[a, b]$we calculate in each $$$f$, average the results and multiply by $$$b-a$. This gives us the approximate value of the integral, as shown in the figure below. $$$\ int _ {- 1} ^ {1} \ sqrt {1 - x ^ 2} dx$with 20 samples approximates the correct result equal to $$$\ frac {\ pi} {2}$.
')

Of course, every time we calculate such an approximate value, we will get a different result. The variance of these values ​​depends on the shape of the function. $$$f$. If we generate random points $$$x_i$unevenly, then we need to slightly change the formula. But thanks to the use of uneven distribution of points, we get a huge advantage: forcing uneven distribution to give preference to points $$$x_i$where $$$f (x)$large, we can significantly reduce the variance of the approximate values. This principle of non-uniform sampling is called sampling by significance .


Since over the past decades, a large-scale transition from deterministic to randomized approaches has taken place in rendering techniques, we will study the randomized approaches used to solve rendering equations. For this we use random variables, mathematical expectation and variance. We are dealing with discrete values, because computers are discrete in nature. Continuous quantities deal with the probability density function , but in the article we will not consider it. We’ll talk about the probability mass function. PMF has two properties:

1. For each $$$s \ in S$exists $$$p (s) \ geq 0$.
2. $$$\ sum_ {s \ in S} p (s) = 1$

The first property is called non-negativity. The second is called "normality." Intuitively, that $$$S$represents the set of results of some experiment, and $$$p (s)$Is the result of probability $$$s$member $$$S$. The outcome is a subset of the probability space. The probability of an outcome is the sum of the PMF elements of this outcome, since

$$

$$Pr \ {E \} = \ sum_ {s \ in S} p (s)$$

A random variable is a function, usually denoted by a capital letter, which sets real numbers in the probability space:

$$

$$X: S \ rightarrow \ boldsymbol {R}.$$

Note that the function $$$X$- This is not a variable, but a function with real values. She is also not random , $$$X (s)$Is a separate real number for any result $$$s \ in S$.

A random variable is used to determine outcomes. For example, many results $$$s$, for which $$$X (s) = 1$, that is, if ht and th are many lines denoting “eagles” or “tails,” then

$$

$$E = {s \ in S: X (s) = 1}$$

and

$$

$$= {ht, th}$$

it is an outcome with probability $$$\ frac {1} {2}$. We write it as $$$Pr \ {X = 1 \} = \ frac {1} {2}$ . We use the predicate $$$X = 1$as a shortened entry for the outcome determined by the predicate.

Let's take a look at a piece of code simulating an experiment described by the formulas presented above:

headcount = 0 if (randb()): // first coin flip headcount++ if (randb()): // second coin flip headcount++ return headcount 

Here we denote by ranb() Boolean function that returns true in half the cases. How is it related to our abstraction? Imagine a lot $$$S$all possible executions of the program, declaring two executions the same values ​​returned by ranb , pairwise identical. This means that there are four possible executions of the program in which two ranb() calls return TT, TF, FT, and FF. From our own experience, we can say that these four accomplishments are equally probable, that is, each occurs in about a quarter of cases.

Now the analogy is becoming clearer. The many possible executions of a program and their associated probabilities are the probability space. Program variables that depend on ranb calls are random variables. I hope you understand everything now.

Let's discuss the expected value, also called the average. This is essentially the sum of the product of PMF and a random variable:

$$

$$E [X] = \ sum_ {s \ in S} p (s) X (s)$$

Imagine that h are “eagles” and t are “tails”. We have already covered ht and th. There are also hh and tt. Therefore, the expected value will be as follows:

$$

$$E [X] = p (hh) X (hh) + p (ht) X (ht) + p (th) X (th) + p (tt) X (tt)$$

$$

$$= \ frac {1} {4}. 2 + \ frac {1} {4}. 1 + \ frac {1} {4}. 1 + \ frac {1} {4} .0$$

$$

$$= 1 \ text {QED}$$

You may wonder where it came from $$$X$. Here I meant that we should assign meaning $$$X$by yourself. In this case, we assigned h to 1, and t to 0. $$$X (hh)$equals 2 because it contains 2 $$$h$.

Let's talk about distribution. The probability distribution is a function that gives the probabilities of various outcomes of an event.

When we say that a random variable $$$X$has a distribution $$$f$then should indicate $$$X \ sim f$.

Scattering values ​​accumulated around $$$X$is called its dispersion and is defined as follows:

$$

$$\ boldsymbol {Var} [X] = E [(X - \ bar {X}) ^ 2]$$

Where $$$\ bar {X}$Is average $$$X$.

$$$\ sqrt {\ boldsymbol {Var}}$called standard deviation . Random variables $$$X$and $$$Y$are called independent if:

$$

$$Pr \ {X = x \ text {and} Y = y \} = Pr \ {X = x \}. Pr \ {Y = y \}$$

Important properties of independent random variables:

1. $$$E [XY] = E [X] E [Y]$
2. $$$\ boldsymbol {Var} [X + Y] = \ boldsymbol {Var} [X] + \ boldsymbol {Var} [Y]$

When I started with a story about probability, I compared continuous and discrete probabilities. We examined discrete probability. Now let's talk about the difference between continuous and discrete probabilities:

1. Values ​​are continuous. That is, the numbers are infinite.
2. Some aspects of analysis require mathematical subtleties such as measurability .
3. Our probability space will be infinite. Instead of PMF, we should use the probability density function (PDF).

PDF Properties:

1. For each $$$s \ in S$we have $$$p (s) \ geq 0$
2. $$$\ int_ {s \ in S} p (s) = 1$

But if the distribution $$$S$evenly , then the pdf is defined like this:


With continuous probability $$$E [X]$defined as follows:

$$

$$E [X]: = \ int_ {s \ in S} p (s) X (s)$$

Now compare the definitions of PMF and PDF:

$$

$$\ mathbb {PMF} \ rightarrow p_y (t) = Pr \ {Y = t \} \ text {for} t \ in T$$

$$

$$\ mathbb {PDF} \ rightarrow Pr \ {a \ leq X \ leq b \} = \ int_a ^ bp (r) dr$$

In the case of continuous probability, random variables are better called random points . Because if $$$S$Is the probability space, and $$$Y: S \ rightarrow T$displayed in a different space than $$$\ mathbb {R}$then we should call $$$Y$random point , not a random variable. The concept of probability density is applicable here, because we can say that for any $$$U \ subset T$we have:


Now let's apply what we have learned to the sphere. The sphere has three coordinates: latitude, longitude, and complement of latitude. We use longitude and latitude addition only in $$$\ mathbb {R} ^ 2$, two-dimensional Cartesian coordinates applied to a random variable $$$S$turn her into $$$S ^ 2$. We get the following detail:

$$

$$Y: [0, 1] \ times [0, 1] \ rightarrow S ^ 2: (u, v) \ rightarrow (\ cos (2 \ pi u) \ sin (\ pi v), \ cos (\ pi v) \ sin (2 \ pi u) sin (\ pi v))$$

We start with a uniform probability density $$$p$at $$$[0, 1] \ times [0, 1]$, or $$$p (u, v) = 1$. Look at the uniform probability density formula above. For convenience, we will write $$$(x, y, z) = Y (u, v)$.

We have an intuitive understanding that if you select points evenly and randomly in a unit square and use $$$f$to convert them to points on a unit sphere, they will accumulate next to the pole. This means that the obtained probability density in $$$T$will not be uniform. This is shown in the figure below.


Now we will discuss ways to approximate the expected value of a continuous random variable and its application to determine the integrals. This is important because in rendering we need to determine the value of the reflectivity integral :

$$

$$L ^ {ref} (P, \ omega_o) = \ int _ {\ omega_i \ in S _ {+} ^ {2}} L (P, - \ omega_i) f_s (P, \ omega_i, \ omega_0) \ omega_i. \ boldsymbol {n} d \ omega_i,$$

for various values $$$P$and $$$\ omega_0$. Value $$$\ omega$Is the direction of the incident light. Code generating a random number uniformly distributed in the interval $$$[0, 1]$and taking the square root, creates a value in the range from 0 to 1. If we use PDF for it, since this is a uniform value, then the expected value will be equal $$$\ frac {2} {3}$. Also this value is the average value $$$f (x) = \ sqrt {x}$in this interval. What does this mean?

Consider Theorem 3.48 from the book Computer Graphics: Principles and Practice. She says that if $$$f: [a, b] \ rightarrow \ mathbb {R}$is a function with real values, and $$$X \ sim \ boldsymbol {U} (a, b)$is a uniform random variable in the interval $$$[a, b]$then $$$(b-a) f (x)$Is a random variable whose expected value has the form:

$$

$$E [(b-a) f (x)] = \ int_a ^ b f (x) dx.$$

What does this tell us? This means that you can use a randomized algorithm to calculate the value of the integral if we execute the code many times and average the results .

In the general case, we get a certain value $$$C$, as in the integral shown above, which needs to be determined, and some randomized algorithm that returns an approximate value $$$C$. Such a random variable for a quantity is called an estimator . An estimator is considered to be distortion- free if its expected value is $$$C$. In the general case, estimators without distortions are preferable to distortions.

We have already discussed discrete and continuous probabilities. But there is a third type, which is called mixed probabilities and is used in rendering. Such probabilities arise due to pulses in the distribution functions of bidirectional scattering, or pulses caused by point sources of illumination. Such probabilities are defined in a continuous set, for example, in the interval $$$[0, 1]$but not strictly defined by the PDF function. Consider the following program:

 if uniform(0, 1) > 0.6 : return 0.3 else : return uniform(0, 1) 

In sixty percent of the cases, the program will return 0.3, and in the remaining 40%, it will return a value evenly distributed in $$$[0, 1]$. The return value is a random variable with a probability mass of 0.6 at 0.3, and its PDF at all other points is specified as $$$d (x) = 0.4$. We must define the pdf as:


In general, a mixed variable random variable is a variable for which there are a finite set of points in the PDF definition area, and vice versa, uniformly distributed points where the PMF is not defined.