Хабрахабр

[Перевод] Теория вероятностей для физически точного рендеринга

image

Введение

В рендеринге часто используется вычисление многомерных определённых интегралов: например, для определения видимости пространственных источников освещения (area light), светимости, доходящей до области пикселя, светимости, поступающей за период времени и облучения, поступающего через полусферу точки поверхности. Вычисление этих интегралов обычно выполняется при помощи интегрирования Монте-Карло, в котором интеграл заменяется ожиданием стохастического эксперимента.

Это будет сделано с практической точки зрения — предполагается, что читатель не сильно знаком с теорией вероятностей, но всё равно хочет разрабатывать эффективные и корректные алгоритмы рендеринга.
В этой статье я подробно расскажу о базовом процессе интегрирования Монте-Карло, а также о нескольких техниках, позволяющих снизить дисперсию методики.

Определённые интегралы

Определённый интеграл — это интеграл вида $\int_^{b}f(x)dx$, где $[a,b]$ — это отрезок (или область), $x$ — скаляр, а $f(x)$ — функция, которая может быть вычислена для любой точки отрезка. Как написано в Википедии, определённый интеграл — это площадь со знаком на плоскости $x$, ограниченная графиком $f$, осью $x$ и вертикальными линиями $x=a$ и $x=b$ (Рисунок 1a).

Эта концепция логично распространяется и на более высокое количество измерений: для определённого двойного интеграла площадь со знаком становится объёмом со знаком (Рисунок 1b), а в общем для определённых многократных интегралов он становится многомерным объёмом со знаком.

Рисунок 1: примеры определённых интегралов.

В иных случаях аналитическое решение невозможно, например, когда нам нужно узнать объём части айсберга, находящейся над водой (Рисунок 1c). В некоторых случаях площадь можно определить аналитически, например, для $f(x)=2$: на отрезке $[a,b]$ площадь будет равна $2(b-a)$. В таких случаях $f(x)$ часто можно определить сэмплированием.

Численное интегрирование

Мы можем приблизительно вычислить площадь сложных интегралов с помощью численного интегрирования. Один из примеров — это сумма Римана. Эта сумма вычисляется разделением области на правильные фигуры (например, прямоугольники), которые вместе образуют область, похожую на истинную область. Сумма Римана задаётся следующим образом:

$\tag{1}S=\sum_{i=1}^{n}f(x_i)\Delta x_i$

$n$ — это количество подинтервалов, а $\Delta x_i=\frac{b-a}{n}$ — ширина одного подинтервала. Для каждого интервала $i$ мы сэмплируем $f$ в одной точке $x_i$ внутри подинтервала (на Рисунке 2 эта точка находится в начале подинтервала).

Рисунок 2: сумма Римана.

Стоит заметить, что при увеличении $n$ сумма Римана сходится к действительному значению интеграла:

$\tag{2}\int_{a}^{b}f(x)dx=\lim_{||\Delta x||\to 0}\sum_{i=1}^{n}f(x_i)\Delta x_i$

Сумма Римана может использоваться и при больших размерностях (Рисунок 3). Однако здесь мы сталкиваемся с проблемой: для функции с двумя параметрами количество подинтервалов должно быть намного больше, если мы хотим достичь разрешения, сравнимого с тем, которое применялось в двухмерном случае. Это явление называют проклятием размерностей, и в более высоких размерностях оно усугубляется.

Рисунок 3: сумма Римана для двойного интеграла.

Сейчас мы оценим точность суммы Римана для следующей функции (мы намеренно выбрали сложную функцию):

$\tag{3}f(x)=\left |\sin \left (\frac{1}{2}x+\frac{\pi}{2} \right )\tan \frac{x}{27}+\sin \left (\frac{3}{5}x^2 \right )+\frac{4}{x+\pi+1}-1\right |$

График функции на отрезке $[-2.5,2.5]$ показан ниже. Для справки мы вычислили в Wolfram Alpha определённый интеграл $\int_{-2.5}^{2.5}f(x)$, получив площадь $3.12970$. График справа показывает точность численного интегрирования с помощью суммы Римана для увеличивающегося $n$.

Рисунок 4: график функции и точность суммы Римана. Даже при небольших $n$ мы получаем достаточно точный результат.

При $n=100$ погрешность составляет $~3\times10^{-4}$. Чтобы получить представление о точности, приведём числа: при $n=50$ погрешность составляет $~2\times10^{-3}$. Следующий порядок величин получается при $n=200$.

Дополнительную информацию о суммах Римана можно прочитать на следующих ресурсах:

Монте-Карло (1)

При рендеринге почти никакие (а может быть, вообще никакие?) интегралы не являются однократными. Это значит, что мы быстро столкнёмся с проклятием размерностей. Кроме того, сэмплирование функции на равных интервалах подвержено недостаточности выборки и искажениям: мы можем пропустить важные значения функции или получить непредусмотренные взаимные помехи между сэмплируемой функцией и паттерном сэмплирования (Рисунок 5).

Рисунок 5: искажения приводят к потере частей сэмплируемой функции (красная) и, в этом случае, к полностью неверной интерпретации функции.

Аналогично сумме Римана, в ней тоже используется сэмплирование функции во множестве точек, но в отличие от детерминированного паттерна суммы Римана мы применяем фундаментально недетерминированный ингредиент: случайные числа. Эти проблемы решаются при помощи техники, называемой интегрированием Монте-Карло.

Интегрирование Монте-Карло основано на следующем наблюдении: интеграл можно заменить ожиданием стохастического эксперимента:

$\tag{4}\int_{a}^{b}f(x)dx=(b-a)E \left [f(X) \right ]\approx \frac{b-a}{n}\sum_{i=1}^{n}f(X)$

Другими словами, мы сэмплируем функцию $n$ раз в случайных точках в пределах отрезка (обозначаемого заглавной $X$), усредняем сэмплы и умножаем на ширину отрезка (для одномерной функции). Как и в случае суммы Римана, при стремлении $n$ к бесконечности среднее значение сэмплов сходится к ожиданию, то есть к истинному значению интеграла.

Немного теории вероятностей

Здесь важно понимать каждое из используемых понятий. Давайте начнём с ожидания: это значение, ожидаемое для одного сэмпла. Учтите, что это не обязательно возможное значение, что может показаться нелогичным. Например, когда мы бросаем кубик, ожидание равно $3.5$ — среднему всех возможных исходов: $(1+2+3+4+5+6)/6=21/6=3.5$.

Это может показаться очевидным, но для интегрирования Монте-Карло нам требуются равномерно распределённые случайные числа, т.е. Второе понятие — это случайные числа. Подробнее об этом мы поговорим позже. каждое значение должно иметь равную вероятность генерации.

Даже когда мы возьмём небольшое количество чисел, ожидаемое среднее значение, а также ожидание каждого отдельного сэмпла должно быть одинаковым. Третье понятие — это отклонение и связанная с ним дисперсия. Отклонение — это разность между ожиданием и исходом эксперимента: $X-E(X)$. Однако при вычислении уравнения 4 мы редко будем получать такое значение.

На практике это отклонение имеет интересное распределение:

Это график нормального распределения, или распределения Гаусса: он показывает, что не все отклонения одинаково вероятны. На самом деле, примерно 68,2% сэмплов находится в интервале $-1\sigma..1\sigma$, где $\sigma$ (сигма) — это стандартное отклонение. Стандартное отклонение можно описать двумя способами:

  • Стандартное отклонение — это мера изменчивости данных.
  • 95% точек данных находятся в пределах $2\sigma$ от среднего.

Для определения стандартного отклонения существует два метода:

  1. Стандартное отклонение $\sigma=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left (X_i-E\left [X\right ]\right )^2}$: его можно вычислить, если есть дискретное распределение вероятностей и известно ожидание $E[X]$. Это истинно для кубиков, у которых $X={1,2,3,4,5,6}$ и $E[X]=3.5$. Подставив числа, получим $\sigma=1.71$.
  2. Также стандартное отклонение сэмплов можно вычислить как $\sigma=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}\left (X_i-X\right )^2}$. Подробнее об этом можно прочитать в Википедии.

Если <img src="https://habrastorage.org/getpro/habr/formulas/f26/e25/4a7/f26e254a766b3149249186fac74821bd.svg" alt="$\sigma=1. Проверка: правильно ли это? Мы знаем, что ${2,3,4,5}$ удовлетворяют этому критерию, а $1$ и $6$ — нет. 71$" data-tex="inline"/>, то мы заявляем, что 68,2% сэмплов находятся в пределах 1,71 от 3,5. Если бы наш кубик мог выдавать любое значение в интервале $[1..  Четыре из шести — это 66,7%. 6]$, то мы бы получили ровно 68,2%.

Вместо стандартного отклонения часто используют связанное понятие дисперсии, которое определяется как $Var\left [X\right ]=\sigma^2$. Поскольку используется квадрат, дисперсия всегда положительна, что помогает в вычислениях.

Монте-Карло (2)

Выше мы приближённо вычислили уравнение 3 при помощи суммы Римана. Теперь мы повторим этот эксперимент с интегрированием Монте-Карло. Вспомним, что интегрирование Монте-Карло задаётся следующим образом:

$\tag{5}\int_{a}^{b}f(x)dx=(b-a)E \left [f(X) \right ]\approx \frac{b-a}{n}\sum_{i=1}^{n}f(X)$

Переведём это в код на C:

double sum = 0;
for( int i = 0; i < n; i++ ) sum += f( Rand( 5 ) - 2.5 );
sum = (sum * 5.0) / (double)n;

Результат для значений от $n=2$ до $n=200$ показан на графике ниже. Из него можно предположить, что интегрирование Монте-Карло проявляет себя гораздо хуже, чем сумма Римана. Более внимательное изучение погрешности говорит, что при $n=200$ средняя погрешность суммы Римана равна $0.0002$, а у Монте-Карло — $0.13$.

Рисунок 6: погрешность Монте-Карло при 2..200 сэмплах.

Показанное ниже уравнение является развёрнутой версией использованного выше, получающее два параметра: В более высоких размерностях эта разница снижается, но не устраняется полностью.

1y\right )\cos\left (2. $f(x,y)=\left |\sin\left( \frac{1}{2}x+\frac{\pi}{2}\right )\tan \frac{x}{27}+\sin \left ( \frac{1}{6}x^2\right )+\frac{4}{x+\pi+1}-1\right |\left |\sin\left ( 1. 3x\right )\right | (6)$

Рисунок 7: график представленного выше уравнения.

5,2. В области определения $x∈[-2. 5,2. 5],y∈[-2. 8685$. 5]$" data-tex="inline"/> объём, ограниченный этой функцией и плоскостью $xy$, равен $6. 043$. При $n=400$ (20×20 сэмплов) погрешность суммы Римана равна $0. 33$. При том же количестве сэмплов средняя погрешность интегрирования Монте-Карло составляет <img src="https://habrastorage.org/getpro/habr/formulas/15e/d37/a0c/15ed37a0c0a7bdea1cb69a1b403d6d45.svg" alt="$0. Чтобы разобраться в этой проблеме, мы изучим хорошо известную технику снижения дисперсии для интегрирования Монте-Карло под названием «стратификация». Это лучше, чем предыдущий результат, но разница всё равно значительна.

Рисунок 8: влияние стратификации; a) сэмплы с плохим распределением; b) сэмплы с равномерным распределением.

На Рисунке 8a для сэмплирования функции используются восемь случайных чисел. Стратификация повышает равномерность случайных чисел. На Рисунке 8b показано влияние стратификации: область определения разделена на восемь страт, и в каждой страте выбирается случайная позиция, что улучшает равномерность. Так как каждое число выбирается случайным образом, часто они оказываются неравномерно распределёнными по области определения.

На Рисунке 9a показан график результатов со стратификацией и без неё. Влияние на дисперсию достаточно очевидно. При $n=10$ средняя погрешность для 8 страт равна <img src="https://habrastorage.org/getpro/habr/formulas/bf6/317/a9a/bf6317a9a463403679fb5e8f3282509a.svg" alt="$0. На Рисунке 9b показана погрешность приблизительного значения. 07$" data-tex="inline"/>, а для 200 страт она снижается до $0. 05$; для 20 страт — $0.  Исходя из этих результатов кажется, что стоит использовать большое количество страт. 002$. Во-первых, количество сэмплов всегда должно быть кратно количеству страт; во-вторых, как и в сумме Римана, стратификация страдает от проклятья размерностей. Однако стратификация обладает недостатками, усиливающимися при увеличении количества страт.

Рисунок 9: стратификация и дисперсия: a) приблизительное значение при количестве сэмплов от n=2 до n=200; b) отклонение.

Выборка по значимости

В предыдущих разделах мы сэмплировали уравнения равномерно. Расширение интегрирующей функции Монте-Карло позволяет нам изменить ситуацию:

$\tag{7}\int_{a}^{b}f(x)dx=(b-a)E \left [f(X) \right ]\approx \frac{b-a}{n}\sum_{i=1}^{n}\frac{f(X)}{p(X)}$

Здесь $p(X)$ — это функция распределения вероятностей (probability density function, pdf): она определяет относительную вероятность того, что случайная величина примет определённое значение.

1$" data-tex="inline"/>, pdf равняется 1 (Рисунок 10a), и это значит, что каждое значение имеет одинаковую вероятность выбора. Для равномерной случайной величины в интервале $0.. 5]$, то получим вероятность в <img src="https://habrastorage.org/getpro/habr/formulas/f68/2d0/19e/f682d019e8fd1f11da97e75f9863c232.svg" alt="$0. Если мы проинтегрируем эту функцию по <img src="https://habrastorage.org/getpro/habr/formulas/a83/cde/062/a83cde0624e801266fa07e48d5532188.svg" alt="$[0,0. Для $X>\frac{1}{2}$ мы, очевидно, получим ту же вероятность. 5$" data-tex="inline"/> того, что $X<\frac{1}{2}$.

Рисунок 10: распределения вероятностей. a) постоянная pdf, при которой каждый сэмп имеет равную вероятность выбора; b) pdf, при которой сэмплы ниже 0.5 имеют бОльшую вероятность выбора.

В данном случае вероятность генерации числа меньше $\frac{1}{2}$ равна 70%. На Рисунке 10b показана другая pdf. Это можно реализовать при помощи следующего фрагмента кода:

float SamplePdf()
{ if (Rand() < 0.7f) return Rand( 0.5f ); else return Rand( 0.5f ) + 0.5f;
}

Эта pdf задаётся следующим образом:

4,if x<\frac{1}{2}\\0. $\tag{8}p(x)=\left\{\begin{matrix}1. 6,otherwise\end{matrix}\right.$

Числа $1.4$ и $0.6$ отражают необходимость того, что вероятность $x<\frac{1}{2}$ была равна 70%. При интегрировании pdf по $[0..\frac{1}{2}]$ даёт $1.4\times\frac{1}{2}$, а $0.6\times\frac{1}{2}$ равно $0.3$. Это иллюстрирует важное требование для всех pdf в целом: результат интегрирования pdf должен быть равен 1. Ещё одно требование заключается в том, что $p(x)$ не может быть равна нулю, если $f(x)$ не равна нулю: это бы значило, что части $f$ имеют нулевую вероятность сэмплирования, что, очевидно, влияет на значение.

Несколько подсказок, позволяющих уяснить концепцию pdf:

  • Одно значение pdf не описывает вероятность: следовательно, локально pdf может быть больше 1 (например, как в только что рассмотренной pdf).
  • Однако интеграл по области определения pdf является вероятностью, а значит, что интегрирование pdf даёт 1.

Одно значение можно интерпретировать как относительную возможность появления конкретного значения.

В случае нормального распределения эта случайная величина является отклонением от среднего. Стоит учесть, что нормальное распределение — это функция распределения вероятностей: оно даёт нам вероятность того, что какая-то случайная величина будет находиться в определённом интервале. Как и любая приличная pdf, результат интегрирования нормального распределения равно 1.

Следовательно, уравнение 7 позволяет нам выполнять неравномерное сэмплирование. Оно компенсирует его делением каждого сэмпла на относительную вероятность его выбора. Важность этого показана на Рисунке 11a. График функции демонстрирует значительный интервал, в котором её значение равно $0$. Сэмплирование в этой области бессмысленно: к сумме ничего не прибавляется, мы просто делим на большее число. Вспомните айсберг с Рисунка 1c: нет смысла сэмплировать высоту в большой площади вокруг айсберга.

Рисунок 11: pdf для функции с нулевыми значениями.

Заметьте, что эта pdf на самом деле равна нулю для интервала значений. Pdf, использующая это знание о функции, показана на Рисунке 11b. Мы можем расширить эту идею за пределы нулевых значений. Это не превращает её в неправильную pdf: в некоторых местах функция равна нулю. На самом деле, идеальная pdf пропорциональна сэмплируемой функции. Сэмплы лучше всего тратить на те места, в которых функция имеет значимые значения. Ещё более качественная pdf показана на Рисунке 12b. Очень хорошая pdf для нашей функции показана на Рисунке 12a. В обоих случаях мы не должны забывать нормализовать её, чтобы интеграл равнялся 1.

Рисунок 12: улучшенные pdf для функции на Рисунке Figure 11.

Pdf на Рисунке 12 ставят перед нами две задачи:

  1. как создать такие pdf;
  2. как сэмплировать такие pdf?

Ответ на оба вопроса одинаков: нам этого делать не нужно. Во многих случаях функция, которую мы хотим интегрировать, неизвестна, и единственный способ определения мест, в которых она значима — это её сэмплирование, и именно для этого нам и нужна pdf; классическая ситуация «курицы и яйца».

В таких случаях очень приблизительная pdf часто лучше, чем никакой pdf. Однако в других случаях у нас есть примерное представление о том, где функция может давать более высокие значения или нулевые значения.

Пара сэмплов дают представление о форме функции, и на основании неё мы нацеливаем последующие сэмплы в те места, где ожидаем высоких значений, которые используем для улучшения pdf, и так далее. Также у нас может иметься возможность создания pdf «на лету».

Серьёзной сложностью является построение pdf. В следующей статье мы применим эти концепции к реализации рендеринга. Мы исследуем несколько случаев, когда pdf помогают в сэмплировании.

Показать больше

Похожие публикации

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»