Несмещной пат-трейсер (18 стр)

Имбирная Ведьмочка
Спасибо за ответы, интересно наблюдать за развитием, желаю успехов!

Одна эпоха рендерится таким образом:
1. Сначала пускается фиксированное количество путей от источников.
2. Каждый раз, когда путь от источника попадает в поверхность — создаётся коннектор между точкой попадания и объективом камеры.
3. Ждём, когда все пути от источников закончатся. Завершённые пути собираются в пул.
4. Для каждого пикселя экрана, который принадлежит этой эпохе (каждая эпоха рендерит только часть пикселей), пускаем по одному пути от камеры.
5. Каждый раз, когда путь от камеры попадает в отражающую поверхность — из пула от-источников рандомно выбирается фиксированное число путей, и для каждого пути, для каждой вершины этого пути создаётся коннектор между этой вершиной и точкой попадания.
6. Когда путь от камеры попадает в излучающую поверхность, её цвет напрямую записывается в пиксель.

Этот процесс регулируется тремя параметрами:
— число пикселей в одной эпохе, N,
— размер пула от-источников в эпохе, R,
— количество выборок из пула на каждый путь от камеры, k.

В контексте мультипл-импортанса, в этом процессе, нам становится важно учесть эти параметры при вычислении вероятностей различных техник. Если посчитать общее количество сэмплов, сгенерированных разными методами в течение одной эпохи, то у нас выйдет:
— R сэмплов с одной вершиной от камеры,
— N сэмплов с нулём вершин от источника,
— k*N сэмплов по всем остальным схемам.
Соответственно, когда мы рассчитываем относительные вероятности альтернатив (для вычисления веса сэмпла в мультипл-импортансе), нам следует учесть соотношения между этими параметрами, а именно:
— В соотношении \(P_{1 \leftrightarrow n} / P_{2 \leftrightarrow n-1}\), помимо обычных коэффициентов для геометрии и материала, у нас дополнительно появится \(k \ N / R\).
— В соотношении \(P_{n \leftrightarrow 1} / P_{n+1 \leftrightarrow 0}\) дополнительно появится \(1/k\).

(картинки с примерами пока что рендерятся)

Влияние учёта k (количество коннекторов на пиксель) на развесовку можно увидеть в этой паре изображений.

Рендер с connections_per_lens_sample = 32, где этот параметр учитывается в развесовке:

Рендер с теми же параметрами, но без учёта в развесовке:

_

Влияние учёта R/N (соотношение свето-сэмплов к пиксель-сэмплам) можно рассмотреть при двух ситуациях — когда это соотношение больше единицы, и когда оно меньше.

light_samples_per_epoch / pixels_per_epoch = 1/32, с учётом в развесовке:

Без учёта:

_

light_samples_per_epoch / pixels_per_epoch = 1024, с учётом в развесовке:

Без учёта:

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

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

Имбирная Ведьмочка
> каждая эпоха рендерит только часть пикселей
Эпоха это bucket? или это последовательность пикселей из общего буфера? или это вообще рандомные пиксели?
> Сначала пускается фиксированное количество путей от источников
Это для всех эпох или для текущей? Если для текущей, то наверное будет различимая разница в шуме, например если эпоха выглядит как отдельный bucket, то на границе разных bucket по идеи будет видимая граница?

Ruslan
> Эпоха это bucket? или это последовательность пикселей из общего буфера? или это
> вообще рандомные пиксели?
Понятия не имею, что за бакет.

В первом приближении, одна эпоха — это один независимый сэмпл. Эпоха может сэмплить все пиксели сразу — тогда она становится по смыслу эквивалентной "фрейму".

Контекст эпохи в коде выглядит так:

struct Epoch {
    std::shared_mutex mutex;
    Arena<PathNode> path_node_arena;
    Arena<Connector> connector_arena;
    AlignedBuffer light_sample_buffer;
    size_t epoch_number;
    LightTriangleDistribution light_triangles_distribution;
    std::atomic<uint64_t> use_count;
    std::atomic<int> phase_one_progress;
    std::atomic<int> completed_light_paths_count;
    std::atomic<int> phase_two_progress;
};

path_node_arena и connector_arena — это арена-аллокаторы. По ходу трассировки путей, временные объекты аллоцируются из этих арен, а удаление происходит через быстрый сброс арены при инициализации эпохи. (типичный арена-аллокатор, короче)

light_sample_buffer — это массив, в который сохраняются завершённые пути от источника.

epoch_number — это номер эпохи по порядку.

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

use_count — это тупо счётчик ссылок. Каждый ray query, находящийся в процессе обработки — одна ссылка. Генераторы первичных луч-запросов от камеры и от источников ещё добавляют по ссылке, пока не закончат работу. Когда use_count.load(std::memory_order_acquire) == 0 — это значит, что вся работа по эпохе завершена, и объект теперь можно сбросить и переинициализировать по-новой.

Остальные переменные управляют работой первичных генераторов. Когда контекст инициализируется — phase_one_progress, completed_light_paths_count и phase_two_progress все зануляются.

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

В реализации, phase_xxx_progress — это количество уже выданных "билетов", и заодно номер следующего "билета", подлежащего выдаче. Когда генератор запрашивает билеты — он производит int iq = phase_xxx_progress.fetch_add(batch_size, std::memory_order_relaxed) и, соответственно, получает билеты с номерами от iq до iq + batch_size. Дополнительно при этом проверяется, что номера билетов не превосходят заранее обозначенного максимума — все номера выше этого предела считаются за "фиктивные билеты" и отбрасываются без обработки.

phase_one_progress — это билетный счётчик для начального генератора путей от источников. На каждый вытянутый билет:
1. Выбирается рандомный треугольник согласно light_triangles_distribution текущей эпохи.
2. Выбирается рандомная точка на этом треугольнике и рандомное направление.
3. Аллоцируется и заполняется PathNode для сгенерированной вершины. В этот PathNode записывается и номер вытянутого билета.
4. Выбранная точка, направление и указатель на PathNode заворачиваются в ray query отправляются на конвеер.

Для каждого такого запроса, когда конвеер выдаст результат:
1. Для полученной точки пересечения (если есть) аллоцируется и заполняется новый PathNode.
2. Производим рассчёт коннектора между точкой попадания и центром камеры. Если он проходит базовые проверки (значение не нулевое, точка попадает в кадр) — аллоцируем Connector и делаем для него вторичный ray query.
3. С учётом свойств материала, бросаем монетку на факт отражения и выбираем направление.
4а. Если отражение есть — делаем новый ray query и отправляем обратно на конвейер трассироваться дальше.
4б. Если отражения не случилось (луч поглотился) — то указатель на последний PathNode в пути сохраняем в light_sample_buffer по индексу, равному номеру билета, из которого этот путь был изначально заспавнен; и делаем completed_light_paths_count.fetch_add(1, std::memory_order_release).

completed_light_paths_count — это, соответственно, количество путей от источников, которые закончили трассировку и готовы к совмещению с путями от камеры. Второй генератор, прежде чем начать тянуть свои билеты, сначала ждёт completed_light_paths_count.load(std::memory_order_acquire) >= light_samples_per_epoch.

phase_two_progress — это билетный счётчик для генератора путей от камеры. Он на каждый вытянутый билет:
1. Детерминированно вычисляет координаты экранного пикселя из номера билета.
2. Для указанного пикселя генерируется направление луча (один пиксель покрывает целую небольшую область изображения).
3. Аллоцируется и заполняется PathNode для первой вершины, в него записываются координаты выбранного экранного пикселя.
4. Отправляется ray query на конвейер.

Результаты запросов для путей от камеры обрабатываются следующим образом:
1. Для полученной точки пересечения (если есть) аллоцируется и заполняется новый PathNode.
2. Проверяем материал в точке попадания на предмет светимости. Если есть — сразу считаем результат и записываем во фреймбуффер.
3. Из пула завершённых путей от источников выбираем несколько рандомных представителей. Для каждого выбранного пути от источника:
    3>. Для каждой вершины в выбранном пути (итерация по односвязному списку из PathNode):
        3>>. Делаем предрассчёт коннектора, делаем проверки, если они проходятся — аллоцируем Connector и отправляем вторичный ray query.
4. С учётом свойств материала, бросаем монетку на факт отражения и выбираем направление.
5а. Если отражение есть — делаем новый ray query и отправляем обратно на конвейер трассироваться дальше.
5б. Если отражения нет — ничего не делаем, путь закончился.

Для создаваемых по ходу дела коннекторов обработка простая:
1. По результату запроса смотрим — есть ли прямая видимость между концами коннектора или нет. Если нету — бросаем и ничего не делаем.
2. Если в коннекторе есть ссылка на PathNode со стороны камеры — то заходим туда и вытаскиваем координаты пикселя, откуда этот путь вылез. Если пути от камеры нет — то вычисляем координаты пикселя исходя из направления коннектора.
3. Записываем значение в коннекторе во фреймбуффер.

Собственно, «детерминированно вычисляет координаты экранного пикселя из номера билета» — тут, на самом деле, есть простор для фантазии. В первых версиях, когда я только реализовал систему — "билета" назначались в соответствие пикселям тупо по порядку — билет 0 идёт на пиксель 0 (по порядку лежания в памяти), билет 1 на пиксель 1, и так далее до компайл-тайм константы pixels_per_epoch. Следующая эпоха продолжит подсчёт — билет 0 на пиксель pixels_per_epoch, билет 1 на пиксель (pixels_per_epoch+1), и так далее. И так пока не будет заполнено всё изображение, после чего счёт начинается с начала:

constexpr int epochs_per_frame =
    (width * height + pixels_per_epoch - 1) / pixels_per_epoch;
<...>
int pixel_index = (epoch_number % epochs_per_frame) * pixels_per_epoch + ticket_number;
int ix = pixel_index % width;
int iy = pixel_index / width;
if (iy >= height) break;

pixel:  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
ep 0:   0  1  2  3  4  5  6  7
ep 1:                           0  1  2  3  4  5  6  7
ep 2:                                                   0  1  2  3  4 | 5..7 discarded
ep 3:   0  1  2  3  4  5  6  7
ep 4:                           0  1  2  3  4  5  6  7
ep 5:                                                   0  1  2  3  4 | 5..7 discarded
...

Полоски, на самом деле, становятся заметы только если прочие параметры настроены криво. Например, если взять light_samples_per_epoch = 1 — по факту, это будет рассчёт прямого освещения для одного точечного источника в рамках целой полоски, и в этом случае границы между эпохами — где у каждой будет свой набор из light_samples_per_epoch путей — будут чётко выделяться на фоне гладких поверхностей. Если же брать что-то более разумное — например, light_samples_per_epoch = pixels_per_epoch — то тогда границы между эпохами становятся незаметными на фоне общего шума.

Тем не менее, в текущей версии, я сделал более равномерное распределение пикселов одной эпохи по экрану:

Вот как выглядит рендер, если вернуть полосочное распределение и поставить говно параметры:

constexpr int pixels_per_epoch = 0x1000;
constexpr int light_samples_per_epoch = 1;
int pixel_index = (epoch_number % epochs_per_frame) * pixels_per_epoch + ticket_number;

Если поставить параметры по-нормальному, однако, то полоски пропадают, даже если их не перемешивать:

constexpr int pixels_per_epoch = 0x1000;
constexpr int light_samples_per_epoch = pixels_per_epoch;
int pixel_index = (epoch_number % epochs_per_frame) * pixels_per_epoch + ticket_number;

А вот так, для контекста, выглядят говёные параметры в текущей версии разбилечивания:

Имбирная Ведьмочка
> Понятия не имею, что за бакет

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

Наверно, дам небольшой математический ликбез о том, что вообще такое "мультипл-импортанс Монте-Карло". Не знаю, правда, кому он назначается — Смайл скорее всего и так его знает, а Алекс, например, как обычно запутается в закорючках и опять ничего не поймёт.

В базовой форме, метод Монте-Карло — это способ численно посчитать интеграл вида

    \(\displaystyle I = \int_\mathbb{D} f[\mathbf{x}] \ \mathrm{d}\mathbf{x}\),

где \(\mathbb{D}\) — это область интегрирования, \(\mathbf{x}\) — координаты точки на этой области, а \(f[\mathbf{x}]\) — это, собственно, подынтегральная функция.

Способ заключается в том, чтобы задать на \(\mathbb{D}\) некое случайное распределение \(\mu\), после чего сделать выборку из этого распределения \(\mathbf{x} \hookleftarrow \mu\) и посчитать:

    \(\displaystyle J[\mathbf{x}] = \frac{f[\mathbf{x}]}{p_\mu[\mathbf{x}]}\),

где \(p_\mu[\mathbf{x}]\) — это плотность распределения \(\mu\) в точке \(\mathbf{x}\).

Магия заключается в том, что математическое ожидание оценки \(J\) — это искомый интеграл \(I\):

    \(\displaystyle \text{mean} \langle J \rangle = \int_\mathbb{D} p_\mu[\mathbf{x}] \ J[\mathbf{x}] \ \mathrm{d}\mathbf{x}\)

    \(\displaystyle \text{mean} \langle J \rangle = \int_\mathbb{D} p_\mu[\mathbf{x}] \ \frac{f[\mathbf{x}]}{p_\mu[\mathbf{x}]} \ \mathrm{d}\mathbf{x}\)

    \(\displaystyle \text{mean} \langle J \rangle = \int_\mathbb{D} f[\mathbf{x}] \ \mathrm{d}\mathbf{x}\)

    \(\displaystyle \text{mean} \langle J \rangle = I\)

В мультипл-импортансе, мы используем сразу несколько распределений \(\mu_1, \mu_2, ..., \mu_i, ..., \mu_N\), каждое из который генерирует свой собственный эстиматор \(J_i[\mathbf{x}] = f[\mathbf{x}] / p_{\mu_i}[\mathbf{x}]\).

Для каждой точки пространства мы задаём "развесовку" \(w_i[\mathbf{x}]\), которая говорит — с каким весом конкретно в точке \(\mathbf{x}\) конкретно i-тый эстиматор вкладывается в общий результат. То есть, мы производим своего рода интерполяцию между несколькими эстиматорами. Соответственно, чтобы это было действительно интерполяцией — в каждой точке области интегрирования, веса эстиматоров должны суммироваться к единице:

    \(\displaystyle \sum_{i=1}^N w_i[\mathbf{x}] = 1\)

В мультипл-мипортансе, на каждую "итерацию" полного алгоритма, для каждого эстиматора в наборе делается независимая выборка \(\mathbf{x}_i \hookleftarrow \mu_i\). Иными словами — мы генерируем целый массив точек \(\mathbf{X} = (\mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_N) \hookleftarrow (\mu_1 \otimes \mu_2 \otimes ... \otimes \mu_N)\). После чего, мы считаем оценку интеграла для каждой из этих точек и суммируем их с локальными весами:

    \(\displaystyle J'_i[\mathbf{x}] = w_i[\mathbf{x}] \ \frac{f[\mathbf{x}]}{p_{\mu_i}[\mathbf{x}]}\)

    \(\displaystyle J_{MI}[[\mathbf{X}]] = \sum_{i=1}^N J'_i[\mathbf{x}_i] = \sum_{i=1}^N w_i[\mathbf{x}_i] \ \frac{f[\mathbf{x}_i]}{p_{\mu_i}[\mathbf{x}_i]}\)

Следует обратить внимание, что оценки комбинируются именно с локальными весами — вычисленные в разных точках для разных эстиматоров в наборе, соответственно, в этом выражении они не обязаны суммироваться в единицу. Отдельно нормализовывать их тоже не следует.

Общая оценка, при всём при этом, всё же сходится к корректному значению, что можно доказать, если явно посчитать математическое ожидание:

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \sum_{i=1}^N \text{mean} \langle J'_i \rangle\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \sum_{i=1}^N \left( \ \int_\mathbb{D} p_{\mu_i}[\mathbf{x}] \ J'_i[\mathbf{x}] \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \sum_{i=1}^N \left( \ \int_\mathbb{D} p_{\mu_i}[\mathbf{x}] \ w_i[\mathbf{x}] \ \frac{f[\mathbf{x}]}{p_{\mu_i}[\mathbf{x}]} \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \sum_{i=1}^N \left( \ \int_\mathbb{D} w_i[\mathbf{x}] \ f[\mathbf{x}] \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \int_\mathbb{D} \left( \ \sum_{i=1}^N \left( w_i[\mathbf{x}] \ f[\mathbf{x}] \right) \ \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \int_\mathbb{D} \left( \ \left( \sum_{i=1}^N w_i[\mathbf{x}] \right) \cdot f[\mathbf{x}] \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = \int_\mathbb{D} \left( \ 1 \cdot f[\mathbf{x}] \ \mathrm{d}\mathbf{x} \ \right)\)

    \(\displaystyle \text{mean} \langle J_{MI} \rangle = I\)

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

— Проделать весь K-путь, начиная от камеры, и просто случайно в конце попасть в светящийся треугольник — такой путь можно обозначить схемой \(K \ \ || \ \ 0\), что означает — K вершин, сгенерированные от камеры, плюс 0 вершин, сгенерированные от света.

— Рандомно выбрать 1 вершину среди светящихся треугольников, затем независимо проделать (K-1)-путь от камеры, и затем соединить концы прямым отрезком — это будет схема \(K-1 \ \ || \ \ 1\).

— Рандомно выбрать 1 вершину среди светящихся треугольников, пустить из неё луч и таким образом получить 2-ю вершину от света, затем независимо проделать (K-2)-путь от камеры, и затем соединить концы прямым отрезком — это будет схема \(K-2 \ \ || \ \ 2\).

— Соответственно, ещё целая группа схема \(K-n \ \ || \ \ n\) — когда мы независимо генерируем два пути, длиной n вершин от света и (K-n) вершин от камеры, и соединяем концы отрезком.

— Схема \(1 \ \ || \ \ K-1\) — когда мы (K-1)-путь от света соединяем напрямую с центром камеры.

А вот схема \(0 \ \ || \ \ K\) у меня отсутствует — она бы соответствовала ситуации, когда у камеры есть физически смоделированный объектив конечного размера, и путь от света случайно попал прямо на поверхность объектива.

Соответственно, для каждого сгенерированного K-пути, все возможные схемы — это и есть независимые эстиматоры, которые я комбинирую через мультипл-импортанс.

Какие именно брать распределения и как назначать веса — можно выбирать относительно свободно. Железных условий в этом только два:
1. Сумма весов в любой точке должна быть равна единице.
2. Любая точка, где подынтегральная функция ненулевая \(f[\mathbf{x}] \neq 0\), должна быть "достижимой" хотя бы из одного из исходных распределений.

Достаточно распространённая практика — это вычислять веса через плотности вероятности исходных распределений, с той мыслью, что если в какой-то малой области у одного из распределений плотность вероятности намного выше, чем у остальных — то соответствующий ему эстиматор будет более "детально" рассматривать конкретно эту область, и ему следует дать соответственно более высокий вес.

Конкретный метод, который я использую у себя — это брать веса пропорционально квадратам плотности вероятности:

    \(\displaystyle w_i[\mathbf{x}] = \frac{(p_{\mu_i}[\mathbf{x}])^2}{\sum_{k=1}^N (p_{\mu_k}[\mathbf{x}])^2}\)

Ну и, собственно, описанная выше модификация — это что у путей по схеме \(n || 0\) вероятность дополнительно умножается на 1 / connections_per_lens_sample, а у схемы \(1 || n\) соответственно на 1 / connections_per_lens_sample * ((float)light_samples_per_epoch / pixels_per_epoch), ввиду того, что эти пути генерируются принципиально другими методами и в другом количестве.

Имбирная Ведьмочка
> light_triangles_distribution — это дискретное распределение светящихся
> треугольников. Когда трассируется путь от источника — мы делаем выборку из
> этого распределения, получаем айди треугольника в геометрии, и затем выбираем
> из равномерного распределения по поверхности этого треугольника.
Вы могли бы показать код? Если это не секрет. Когда мне нужно было делать равномерную выборку для сэмплирования  lightMesh, то я растерезировал развертку меша в текстуру, далее собирал все не пустые тексели в буфер, и потом рандомно перемешивал эти тексели в буфере. Это работало, но как по мне  это не очень то оптимально, может подсмотрю для себя чего полезного?

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

}:+()___ [Smile]
> Почему не просто выбор треугольника и случайной точки в нем?
Треугольники то  имеют разную площадь, получится, что из самого большого треугольника например будет браться меньше сэмплов относительно его площади.
Пока писал ответ, понял, что можно продублировать в буфере индекс треугольника соответственно его площали относительно остальных, но даже так будет огромная погрешность хм.

Ruslan
> Вы могли бы показать код?
Как Смайл и сказал — выбирается треугольник по произвольному дискретному распределению, а внутри него — тупо равномерно рандомная точка.

Ruslan
> Треугольники то имеют разную площадь, получится, что из самого большого
> треугольника например будет браться меньше сэмплов относительно его площади.
Распределение между треугольниками не равномерное, в отдельном массиве для каждого треугольника явно задаём вероятность.

Для выборки из этого распределения конкретно у меня используется https://en.wikipedia.org/wiki/Alias_method, но это уже детали реализации.