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

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

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

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

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

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

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

(1)S=∑i=1nf(xi)Δxi
n — это количество подинтервалов, а Δxi=b−an — ширина одного подинтервала. Для каждого интервала i мы сэмплируем f в одной точке xi внутри подинтервала.

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

(2)∫abf(x)dx=lim||Δx||→0∑i=1nf(xi)Δxi
Сумма Римана может использоваться и при больших размерностях. Однако здесь мы сталкиваемся с проблемой: для функции с двумя параметрами количество подинтервалов должно быть намного больше, если мы хотим достичь разрешения, сравнимого с тем, которое применялось в двухмерном случае. Это явление называют проклятием размерностей, и в более высоких размерностях оно усугубляется.

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

(3)f(x)=|sin⁡(12x+π2)tan⁡x27+sin⁡(35×2)+4x+π+1−1|
График функции на отрезке [−2.5,2.5] показан ниже. Для справки мы вычислили в Wolfram Alpha определённый интеграл ∫−2.52.5f(x), получив площадь 3.12970. График справа показывает точность численного интегрирования с помощью суммы Римана для увеличивающегося n.

Чтобы получить представление о точности, приведём числа: при n=50 погрешность составляет  2×10−3. При n=100 погрешность составляет  3×10−4. Следующий порядок величин получается при n=200.

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

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

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

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

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

(4)∫abf(x)dx=(b−a)E[f(X)]≈b−an∑i=1nf(X)
Другими словами, мы сэмплируем функцию n раз в случайных точках в пределах отрезка (обозначаемого заглавной X), усредняем сэмплы и умножаем на ширину отрезка (для одномерной функции). Как и в случае суммы Римана, при стремлении n к бесконечности среднее значение сэмплов сходится к ожиданию, то есть к истинному значению интеграла.

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

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

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

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

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

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

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

  1. Стандартное отклонение σ=1n∑i=1n(Xi−E[X])2: его можно вычислить, если есть дискретное распределение вероятностей и известно ожидание E[X]. Это истинно для кубиков, у которых X=1,2,3,4,5,6 и E[X]=3.5. Подставив числа, получим σ=1.71.
  2. Также стандартное отклонение сэмплов можно вычислить как σ=1n−1∑i=1n(Xi−X)2. Подробнее об этом можно прочитать в Википедии.

Проверка: правильно ли это? Если 

σ=1.71, то мы заявляем, что 68,2% сэмплов находятся в пределах 1,71 от 3,5. Мы знаем, что  2,3,4,5 удовлетворяют этому критерию, а 1 и  6 — нет. Четыре из шести — это 66,7%. Если бы наш кубик мог выдавать любое значение в интервале  [1..6], то мы бы получили ровно 68,2%.
Вместо стандартного отклонения часто используют связанное понятие дисперсии, которое определяется как Var[X]=σ2. Поскольку используется квадрат, дисперсия всегда положительна, что помогает в вычислениях.

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

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

(5)∫abf(x)dx=(b−a)E[f(X)]≈b−an∑i=1nf(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.

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

f(x,y)=|sin⁡(12x+π2)tan⁡x27+sin⁡(16×2)+4x+π+1−1||sin⁡(1.1y)cos⁡(2.3x)|(6)

В области определения x∈[−2.5,2.5],y∈[−2.5,2.5] объём, ограниченный этой функцией и плоскостью xy, равен 6.8685. При n=400 (20×20 сэмплов) погрешность суммы Римана равна 0.043. При том же количестве сэмплов средняя погрешность интегрирования Монте-Карло составляет 0.33. Это лучше, чем предыдущий результат, но разница всё равно значительна. Чтобы разобраться в этой проблеме, мы изучим хорошо известную технику снижения дисперсии для интегрирования Монте-Карло под названием «стратификация».

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

Влияние на дисперсию достаточно очевидно. При n=10 средняя погрешность для 8 страт равна 0.05; для 20 страт — 0.07, а для 200 страт она снижается до 0.002. Исходя из этих результатов кажется, что стоит использовать большое количество страт. Однако стратификация обладает недостатками, усиливающимися при увеличении количества страт. Во-первых, количество сэмплов всегда должно быть кратно количеству страт; во-вторых, как и в сумме Римана, стратификация страдает от проклятья размерностей.