31 July 2019

Применение интегрирования Монте-Карло в рендеринге

Working with 3D-graphicsAlgorithmsMathematics
Translation
Original author: Chubak Bidpaa·
Все мы изучали в курсе математики численные методы. Это такие методы, как интегрирование, интерполяция, ряды и так далее. Существует два вида числовых методов: детерминированные и рандомизированные.

Типичный детерминированный метод интегрирования функции $f$ в интервале $[a, b]$ выглядит так: мы берём $n + 1$ равномерно расположенных в интервале точек $t_0 = a, t_1 = a + \frac{b - a }{n}, \ldots, t_n - b$, вычисляем $f$ в средней точке $\frac{t_i + t_{i + 1}}{2}$ каждого из интервалов, определяемых этими точками, суммируем результаты и умножаем на ширину каждого интервала $\frac{b -a}{b}$. Для достаточно непрерывных функций $f$ при увеличении $n$ результат будет сходиться к верному значению.


Вероятностный метод, или метод Монте-Карло для вычисления, или, если точнее, приблизительной оценки интеграла $f$ в интервале $[a, b]$, выглядит так: пусть $X_1, \ldots, X_n$ — случайно выбранные точки в интервале $[a, b]$. Тогда $Y = (b - a) \frac{1}{n}\sum_{i = 1}^{n}f(X_i)$ — это случайное значение, среднее которого является интегралом $\int_{[a,b]}f$. Для реализации метода мы используем генератор случайных чисел, генерирующий $n$ точек в интервале $[a, b]$, вычисляем в каждой $f$, усредняем результаты и умножаем на $b-a$. Это даёт нам приблизительное значение интеграла, как показано на рисунке ниже. $\int_{-1}^{1}\sqrt{1 - x^2} dx$ с 20 сэмплами аппроксимирует верный результат, равный $\frac{\pi}{2}$.


Разумеется, каждый раз, когда мы будем вычислять такое приблизительное значение, то будем получать разный результат. Дисперсия этих значений зависит от формы функции $f$. Если мы генерируем случайные точки $x_i$ неравномерно, то нам необходимо слегка изменить формулу. Но благодаря использованию неравномерного распределения точек мы получаем огромное преимущество: заставив неравномерное распределение отдавать предпочтение точкам $x_i$, где $f(x)$ велика, мы можем значительно снизить дисперсию приблизительных значений. Такое принцип неравномерной дискретизации называется выборкой по значимости.


Так как за последние десятилетия в методиках рендеринга произошёл масштабный переход от детерминированных к рандомизированным подходам, мы изучим рандомизируемые подходы, применяемые для решения уравнений рендеринга. Для этого мы используем случайные величины, математическое ожидание и дисперсию. Мы имеем дело с дискретными значениями, потому что компьютеры дискретны по своей сути. Непрерывные величины имеют дело с функцией плотности вероятности, но в статье мы не будем её рассматривать. Мы поговорим о функции распределения масс (probability mass function). PMF обладает двумя свойствами:

  1. Для каждого $s \in S$ существует $p(s) \geq 0$.
  2. $\sum_{s \in S}p(s) = 1$

Первое свойство называется «неотрицательностью». Второе называется «нормальностью». Интуитивно понятно, что $S$ представляет собой множество результатов некоторого эксперимента, а $p(s)$ — это результат вероятности $s$, член $S$. Исход — это подмножество пространства вероятностей. Вероятность исхода является суммой PMF элементов этого исхода, поскольку

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


Случайная переменная — это функция, обычно обозначаемая заглавной буквой, ставящая в соответствие пространству вероятностей вещественные числа:

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


Учтите, что функция $X$ — это не переменная, а функция с вещественными значениями. Она также не является случайной, $X(s)$ — это отдельное вещественное число для любого результата $s\in S$.

Случайная переменная используется для определения исходов. Например, множество результата $s$, для которого $X(s) = 1$, то есть если ht и th — это множество строк, обозначающих «орлы» или «решки», то

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


и

$ = {ht, th} $



это исход с вероятностью $\frac{1}{2}$. Запишем это как $Pr\{X=1\} = \frac{1}{2}$. Мы используем предикат $X=1$ как укороченную запись для исхода, определяемого предикатом.

Давайте взглянем на фрагмент кода, симулирующий эксперимент, описанный представленными выше формулами:

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

Здесь мы обозначаем как ranb() булеву функцию, которая возвращает true в половине случаев. Как она связана с нашей абстракцией? Представьте множество $S$ всех возможных выполнений программы, объявив два выполнения одинаковыми значениями, возвращаемыми ranb, попарно идентичными. Это значит, что существует четыре возможных выполнений программы, в которых два вызова ranb() возвращают TT, TF, FT и FF. По своему опыту мы можем сказать, что эти четыре выполнения равновероятны, то есть каждое встречается примерно в четверти случаев.

Теперь аналогия становится понятнее. Множество возможных выполнений программы и связанные с ними вероятности — это пространство вероятностей. Переменные программы, зависящие от вызовов ranb, — это случайные переменные. Надеюсь, теперь вам всё понятно.

Давайте обсудим ожидаемое значение, также называемое средним. По сути это сумма произведения PMF и случайной переменной:

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


Представьте, что h — это «орлы», а t — «решки». Мы уже рассмотрели ht и th. Также существуют hh и tt. Поэтому ожидаемое значение будет следующим:

$ 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} $


Вы можете задаться вопросом, откуда взялся $X$. Здесь я имел в виду, что мы должны назначать значение $X$ самостоятельно. В данном случае мы присвоили h значение 1, а t значение 0. $X(hh)$ равно 2, потому что в ней содержится 2 $h$.

Давайте поговорим о распределении. Распределение вероятностей — это функция, дающая вероятности различных исходов события.

Когда мы говорим, что случайная переменная $X$ имеет распределение $f$, то должны обозначить $X \sim f$.

Рассеяние значений, скопившихся вокруг $X$, называется её дисперсией и определяется следующим образом:

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


Где $\bar{X}$ — это среднее $X$.

$\sqrt{\boldsymbol{Var}}$ называется стандартным отклонением. Случайные переменные $X$ и $Y$ называются независимыми, если:

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


Важные свойства независимых случайных переменных:

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

Когда я начал с рассказа о вероятности, то сравнивал непрерывную и дискретную вероятности. Мы рассмотрели дискретную вероятность. Теперь поговорим о разнице между непрерывной и дискретной вероятностями:

  1. Значения непрерывны. То есть числа бесконечны.
  2. Некоторые аспекты анализа требуют таких математических тонкостей, как измеряемость.
  3. Наше пространство вероятностей будет бесконечным. Вместо PMF мы должны использовать функцию плотности вероятностей (PDF).

Свойства PDF:

  1. Для каждого $s \in S$ у нас есть $p(s) \geq 0$
  2. $\int_{s\in S}p(s) = 1$

Но если распределение $S$ равномерно, то PDF определяется так:

image

При непрерывной вероятности $E[X]$ определяется следующим образом:

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


Теперь сравним определения PMF и 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 $


В случае непрерывной вероятности случайные величины лучше называть случайными точками. Потому что если $S$ — пространство вероятностей, а $Y : S \rightarrow T$ отображается в другое пространство, отличающееся от $\mathbb{R}$, тогда мы должны назвать $Y$ случайной точкой, а не случайной величиной. Понятие плотности вероятностей применимо здесь, потому что можно сказать, что для любого $U \subset T$ мы имеем:

image

Теперь давайте применим то, что мы узнали, к сфере. Сфера имеет три координаты: широту, долготу и дополнение широты. Долготу и дополнение широты мы используем только в $\mathbb{R}^2$, двухмерные декартовы координаты, применённые к случайной величине $S$, превращают её в $S^2$. Получаем следующую детализацию:

$ 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)) $


Мы начинаем с равномерной плотности вероятностей $p$ при $[0, 1] \times [0, 1]$, или $p(u, v) = 1$. Посмотрите выше формулу плотности равномерной вероятности. Для удобства мы запишем $(x, y, z) = Y(u, v)$.

У нас есть интуитивное понимание, что если выбирать точки равномерно и случайно в единичном квадрате и использовать $f$ для преобразования их в точки на единичной сфере, то они будут скапливаться рядом с полюсом. Это означает, что полученная плотность вероятностей в $T$ не будет равномерной. Это показано на рисунке ниже.


Теперь мы обсудим способы приблизительного определения ожидаемого значения непрерывной случайной величины и его применения для определения интегралов. Это важно, потому что в рендеринге нам нужно определять значение интеграла отражающей способности:

$ 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, $


для различных значений $P$ и $\omega_0$. Значение $\omega$ — это направление падающего света. Код, генерирующий случайное число, равномерно распределённое в интервале $[0, 1]$ и берущий квадратный корень, создаёт значение в интервале от 0 до 1. Если мы используем для него PDF, поскольку это равномерное значение, то ожидаемое значение будет равно $\frac{2}{3}$. Также это значение является средним значением $f(x) = \sqrt{x}$ в этом интервале. Что это означает?

Рассмотрим теорему 3.48 из книги «Computer Graphics: Principles and Practice». Она гласит, что если $f : [a, b] \rightarrow \mathbb{R}$ является функцией с вещественными значениями, а $X \sim \boldsymbol{U}(a, b)$ является равномерной случайной величиной в интервале $[a, b]$, то $(b-a)f(x)$ — это случайная величина, ожидаемое значение которой имеет вид:

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


Что это нам говорит? Это значит, что можно использовать рандомизированный алгоритм для вычисления значения интеграла, если мы достаточно много раз выполним код и усредним результаты.

В общем случае мы получим некую величину $C$, как в показанном выше интеграле, которую нужно определить, и некий рандомизированный алгоритм, возвращающий приблизительное значение $C$. Такая случайная переменная для величины называется эстиматором. Считается, что эстиматор без искажений, если его ожидаемое значение равно $C$. В общем случае эстиматоры без искажений предпочтительнее, чем с искажениями.

Мы уже обсудили дискретные и непрерывные вероятности. Но существует и третий тип, который называется смешанными вероятностями и используется в рендеринге. Такие вероятности возникают вследствие импульсов в функциях распределения двунаправленного рассеяния, или импульсов, вызванных точечными источниками освещения. Такие вероятности определены в непрерывном множестве, например, в интервале $[0, 1]$, но не определены строго функцией PDF. Рассмотрим такую программу:

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

В шестидесяти процентах случаев программа будет возвращать 0.3, а в оставшихся 40% она будет возвращать значение, равномерно распределённое в $[0, 1]$. Возвращаемое значение — это случайная переменная, имеющая при 0.3 массу вероятности 0.6, а его PDF во всех других точках задаётся как $d(x) = 0.4$. Мы должны определить PDF как:

image

В целом, случайная переменная со смешанной вероятностью — это такая переменная, для которой существует конечное множество точек в области определения PDF, и наоборот, равномерно распределённые точки, где PMF не определена.
Tags:метод монте-карлочисленные методыдискретная математикарендеринг графики
Hubs: Working with 3D-graphics Algorithms Mathematics
+13
3.4k 48
Comments 1