How to become an author
Search
Write a publication
Pull to refresh

Простой способ построения изотропных псевдослучайных рядов точек синего шума

Reading time20 min
Views3K
Algorithms*Mathematics*
Translation
Original author: Martin Roberts
В данной статье я описываю простой способ создания ряда точек, основанных на псевдослучайном ряде с низким расхождением, демонстрирующем улучшенные изотропные свойства синего шума. Он обеспечивает высокую скорость схождения с минимальными артефактами алиасинга.


Рисунок 1. Первые 100, 200, 500, 1000, 2000 и 5000 точек выборки из предлагаемого прогрессивного нестохастического ряда точек (уравнение 11), демонстрирующие почти изотропные характеристики синего шума с быстрой сходимостью QMC и сниженным количеством артефактов. Ряд основан на новой простой псевдослучайной последовательности с низким расхождением $R_2$.

Введение


Псевдослучайные ряды с низким расхождением применяются для создания распределений, являющихся менее регулярными, чем решётки, но более регулярными, чем простые случайные распределения (см. рисунок 2). Они играют важную роль во множестве областей численных вычислений, в том числе в физике, финансах, а в последнее время — и в компьютерной графике.


Рисунок 2. Сравнение регулярной решётки (слева) с тремя разными псевдослучайными функциями (посередине) и простым случайным распределением (справа). Заметьте, что псевдослучайные распределения выглядят менее регулярными, чем решётка, но имеют не так много скоплений и разрежений точек, как случайное распределение.

Классической статьёй о том, как использование псевдослучайного сэмплирования вместо равномерного сэмплирования может улучшить 2D-изображения, является Stochastic Sampling in Computer Graphics [Cook, 1986]. За последние десятилетия псевдослучайное сэмплирование стало важным инструментом 3D-рендеринга. Рендеринг (и сэмплеры) могут использовать для вычисления или ЦП, или GPU. Превосходным примером их использования в 3D-рендеринге является Anti-aliased Low Discrepancy Samplers for Monte Carlo Estimators in Physically Based Rendering [Perrier, 2108]. В этой статье автор рассказывает о задаче 3D-рендеринга (см. рисунок 3) следующим образом:

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

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

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

В последнее время предпринимались попытки объединить преимущества сглаживания флуктуационного сэмплирования с высокой скоростью схождения (низкой дисперсией), характерной для псевдослучайных последовательностей с низким расхождением [Christensen, 2018]. С той же самой целью я представил новый способ построения бесконечного ряда точек, являющийся чрезвычайно простым, быстрым в вычислении, прогрессивным, полностью детерминированным, который демонстрирует быструю сходимость (низкую дисперсию) с нужным изотропным спектром синего шума, а следовательно, приводит к минимальным артефактам.

Итак, приступаем…

Рисунок 3. Рендеринг в 3D (слева). В процессе сэмплирования часто возникают нежелательные артефакты алиасинга (справа). Источник :[Perrier, 2018]

Краткое содержание


  • простые флуктуационные (конечные) множества точек
  • флуктуационные псевдослучайные множества точек
  • флуктуации (бесконечных) псевдослучайных рядов
  • сравнения сходимости QMC
  • вывод
  • пример кода

Простое флуктуационное сэмплирование


Наиболее простыми флуктуационными (возмущёнными) множествами точек являются флуктуационные сетки. В нашем случае пространство сэмплирования подразделено на сетку из $m \times m$ квадратных областей. Вместо простого выбора точки в центре каждой ячейки выбирается любая (псевдо-)случайная точка внутри каждой ячейки. Преимущества флуктуационного сэмплирования сетки заключаются в простоте понимания и реализации, а также в чрезвычайной скорости вычисления. Однако оно имеет очевидное ограничение — $N$ должен быть идеальным квадратом.

Хотя прямоугольное типологическое сэмплирование может снизить строгость этого ограничения так, чтобы $N$ был произведением двух соответствующих целых чисел, я опишу другой подход, который информативен, очень изящен и намного менее известен. Представим без потери обобщённости, что пространство сэмплирования является единичным квадратом $[0,1)^2$.

Мы определим решётку Фибоначчи как решётку из $N$ точек, где:

$\pmb{p}_i(x,y) = \left( \frac{i-0.5}{N}, \left\{ \frac{i}{\phi} \right\}\; \right); \quad \textrm{for } \; i=1,2,3,…, N. \tag{1}$


где $\phi \simeq 1.618033988$ — золотое сечение, а $\{ x \}$ обозначает дробную часть $x$. Учтите, что при вычислениях этот оператор дробности обычно записывается и реализуется как $x \; \%1$ (произносится как «mod 1»), и поскольку этот пост в основном направлен на вычислительные области применения, в дальнейшем мы будем использовать эту запись.

Стоит заметить, что решётка Фибоначчи — это очень изящный способ размещения любого количества точек в единичном квадрате чрезвычайно регулярным образом. Например, на рисунке 2a показана стандартная решётка Фибоначчи для N=41 точек. При флуктуационном сэмплировании вместо выбора в качестве сэмплируемых точек центров каждого диска мы можем просто выбрать точку $\tilde{p_i}$ внутри каждого из $N$ дисков (см. рисунок 4b). Учитывая, что типичное расстояние между точками равно $\simeq 0.92/\sqrt{N}$,

$\tilde{\pmb{p}_i}(x,y) = \pmb{p}_i(x,y) + \pmb{\epsilon}_i (\delta_0/ \sqrt{N}),\quad \textrm{where } \delta \simeq 0.46 \tag{2}$


где $\epsilon_i(r)$ означает, что для $i$-того члена случайный сэмпл берётся из центрированной в нуле окружности с радиусом $r$.

Следует заметить, что простой прямой способ равномерного взятия сэмпла из окружности радиусом $r$ заключается в том, чтобы сначала взять два равномерных сэмпла $u_1$ и $u_2$, каждый из которых находится в интервале $[0,1)$, а затем применить следующее преобразование:

$\pmb{\epsilon}_i(r) = \left( r \sqrt{u_1} \cos(2\pi u_2), r \sqrt{u_1} \sin(2 \pi u_2) \right) \tag{3}$


Разумеется, можно снизить степень флуктуаций, просто выбрав меньшее значение $\delta$.


Рисунок 4-1. (a) Каноническая решётка Фибоначчи; (b) флуктуационная решётка Фибоначчи.

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


Рисунок 4-2. (a) Каноническая решётка Фибоначчи; (b) флуктуационная решётка Фибоначчи; и (с) случайные сэмплы.

Если $N$ известно заранее, то это простой, но эффективный способ создания изотропного множества точек.

Интересно то, что несмотря на интуитивное ощущение наличия между соседними точками постоянного расстояния, спектр Фурье для флуктуационного множества точек Фибоначчи не демонстрирует спектра классического синего шума. Тем не менее, эта флуктуационная версия очень хорошо подходит во многих областях применения и контекстах!

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

Более того, несмотря на чрезвычайно высокую скорость вычисления флуктуационной сетки и флуктуационной решётки Фибоначчи, они обе страдают от упорядочивания точек «строка за строкой». Часто более предпочтительны способы сэмпирования, при которых последовательные точки не находятся поблизости, а скорее «равномерно распределены» по всему пространству сэмплирования. Это свойство часто называют «прогрессивным сэмплированием».

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

Флуктуационные псевдослучайные множества точек


Псевдослучайный ряд с низкой расходимостью — это (бесконечный) ряд точек, распределённых более равномерно, чем простое случайное распределение, но менее регулярно, чем решётки. Такие ряды полезны во многих областях применения, например, при интегрировании методом Монте-Карло. Существует множество типов псевдослучайных рядов, в том числе Вейла/Кронекера, ван дер Корпута/Холтона, Нидеррайтера и Соболя. Я подробно рассказывал о них, в том числе и о новом, в своём предыдущем посте: The unreasonable effectiveness of quasirandom sequences. В этом посте мы рассмотрим только один из описанных рядов — ряд $R_2$. Для сравнения мы также рассмотрим ряд Холтона и ряд Соболя, поскольку сегодня их чаще всего используют в QMC.

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

$\pmb{R_2}$: ряд $R_2$ — это рекурсивный ряд Кронекера/Вейла/Рихтмайера, представленный в моём предыдущем посте. Он основан на обобщении золотого сечения, а потому может рассматриваться как естественный способ построения обобщённой решётки Фибоначчи, соединяющейся слева направо и сверху вниз. Он задаётся следующим образом:

$\pmb{t}_n = n \pmb{\alpha} \; \% 1, \quad n=1,2,3,… \tag{4}$


$\textrm{where} \quad \pmb{\alpha} =\left( \frac{1}{\varphi}, \frac{1}{\varphi^2} \right) \;\; \textrm{and} \; \varphi\ \textrm{is the unique positive root of } x^3=x+1.$


Стоит заметить, что $\varphi = 1.324717957244746 \ldots$, часто называемая пластической константой, имеет множество красивых свойств (также см. здесь).

Ряд Холтона: $d$-мерный ряд Холтона является независимыми от $d$ рядами обратного основания ван дер Корпута с покомпонентным сочленением. В этом посте мы будем придерживаться традиционной схемы с выбором в качестве оснований первых $d$ простых чисел. Следовательно, в этом посте я всегда буду использовать {2,3}-ряд Холтона.

Ряд Соболя: ряд Соболя имеет множество превосходных свойств низкого расхождения, в том числе высокие скорости схождения. В этом посте построение ряда Соболя (и выбор параметров) будет выполняться при помощи Mathematica версии 10. В документации Mathematica говорится, что реализация Соболя является «закрытыми, проприетарными и полностью оптимизированными генераторами из библиотеки MKL Intel».

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

На рисунке 5 показаны эти три ряда с низким расхождением вместе с их Фурье-преобразованиями и радиальными спектральными функциями. Мы можем это увидеть:

  • Точки $R_2$ очень равномерно распределены. Для многих областей применения такое распределение идеально, но для других областей оно «слишком равномерно».
  • во всех рядах с низким расхождением заметны чёткие паттерны, которые часто приводят к нежелательным артефактам, возникающим на отрендеренных при помощи QMC изображениях.
  • спектры Фурье для всех рядов с низким расхождением состоят из множества пространственно регулярных пиков
  • ни один из псевдослучайных множеств точек не является изотропным
  • радиальная спектральная функция ряда $R_2$ состоит только из дискретного количества пиков.
  • радиальные спектральные функции рядов Холтона и Соболя непрерывны, но не гладки.

То есть, артефакты рендеринга являются непосредственным результатом (i) дискретных ярких пиков в спектрах Фурье и/или (ii) дискретных пиков в радиальных спектральных функциях.

Краткое примечание о том, почему здесь может быть полезен термин «синий шум».

В общем смысле распределение $R_2$ можно назвать синим шумом, потому что его низкие (синие) частоты очевидно подавлены. Однако он состоит только из дискретных пиков, поэтому, возможно, лучше называть его «дискретным синим шумом» (или не называть его синим шумом вообще). В более общем контексте для «синего шума» требуются (i) изотропные спектры Фурье и (ii) радиальная спектральная функция, плоская как у белого шума, но с полностью подавленными низкими частотами. Превосходный обзор различных спектров синего шума см. в [Heck, 2012]

Кроме того, некоторые авторы [Perrier, 2018] замечают, что в большинстве способов сэмплирования диском Пуассона, например, в best candidate Митчелла, используется не строго синий шум, потому что их спектральные функции малы, но не равны нулю при наименьших частотах. В этом смысле лучше называть его «аппроксимированным синим шумом». Также замечают, что разница между почти нулём и полным нулём может сильно отрицательно влиять на скорость схождения.


Рисунок 5. Множества точек (n=2000), основанные на разных псевдослучайных рядах: (i) $R_2$, (ii) Холтона; и (iii) Соболя, со (iv) случайным рядом для сравнения. Спектральные функции ни одного из псевдослучайных рядов не походят ни на случайный белый шум, ни на классический синий шум. Все они имеют заметные и регулярные дискретные пики, и ни один из них не является изотропным — два классических показателя того, что они будут создавать визуальные артефакты.

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

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

Интуитивно понятно, что этот критический уровень флуктуаций будет тесно связан с типичным расстоянием (поделённым на $1/\sqrt{N}$) между соседними точками. На рисунке 6a показано минимальное расстояние между соседними точками для рядов $R_2$, Холтона и Соболя. Стоит заметить, что ряд $R_2$ является единственным рядом с низким расхождением, при котором минимальное расстояние между соседними точками снижается со скоростью всего $1/\sqrt{N}$. Для рядов Холтона и Соболя расстояне падает с гораздо большими скоростями, близкими к $1/N$. Аналогично, на рисунке 6b показано, что среднее расстояние между соседними точками $R_2$ не менее $\simeq \delta_0/ \sqrt{N}$ для $\delta_0 /\simeq 0.76$ при всех значений $N$. Эти показатели дают понять, что идеальной будет флуктуация $R_2$ в пределах радиуса в половину этого значения. Эмпирические проверки подтверждают это.


Рисунки 6a, 6b. Метрика ближайшего соседа для ряда $R_2$ (синий); Холтона (оранжевый); и Соболя (зелёный). Все показанные расстояния являются абсолютными расстояниями между ближайшими соседями, умноженными на $\sqrt{N}$. Заметьте, что ряд $R_2$ является единственным, у которого минимальное расстояние падает со скоростью всего $1/\sqrt{N}$, у остальных оно падает со скоростью примерно $1/N$. Также он имеет наибольшее среднее расстояние для всех значений $n$.

При помощи этого подхода был эмпирически определён критический радиус флуктуации для $R_2$, Холтона и Соболя, показанные в таблице ниже. Соответствующие им распределения точек показаны на рисунке 5. Заметьте, что флуктуация для ряда Соболя должна снижаться гораздо медленнее чем $\sim 1/\sqrt{N}$, чтобы маскировать его очень отчётливые артефакты, поэтому я не рекомендую его для любой области применения, но я всё равно представил его для полноты изучения.


На рисунке 7 видно, что флуктуационный R2, флуктуационный Холтон и даже флуктуационный Соболь создают распределения точек изотропного синего шума.


Рисунок 7. Флуктуационные множества точек (n=2000), основанные на различных псевдослучайных рядах: (i) $R_2$, (ii) Холтона; и (iii) Соболя, со (iv) случайным для сравнения. Заметьте, что преобразование Фурье для всех псевдослучайных версий теперь изотропно, с небольшими, но наблюдаемыми признаками синего шума. Также заметьте, что флуктуация для ряда Соболя должна быть гораздо больше, чтобы замаскировать очень крупные артефакты. Стоит заметить, что радиальная спектральная функция флуктуационного R2 и флуктуационного Холтона кажется похожей на рандомизированные и флуктуационные сетки Оуэна (не показаны), но их спектры Фурье очень отличаются.

Стоит заметить, что для $R_2$ мы можем определить критический уровень флуктуации иным образом — как максимально возможный уровень флуктуаций, который обеспечивает нулевое значение радиальной спектральной функции для самых низких частот. Эмпирически это определение даёт нам для $R_2$ похожее значение $\delta_0 \simeq 0.76$, но соответствующие уровни для Холтона и Соболя будут гораздо ниже указанных в таблице, и должны падать со скоростью $1/n$, а не $1/\sqrt{n}$, как в случае $R_2$.

В связи с желательными спектрами Фурье часто требуется построение множеств точек с хорошим разделением между соседними точками. Ниже показаны минимальное и среднее расстояние для n=500 точек при различных методах. Помеченные звёздочкой данные взяты из статьи «Progressive Multi-jittered Sample Sequences» [Christensen, 2018], а жирным шрифтом указаны новые результаты. (Стандартные данные для Холтона и Соболя, а также общие данные для случайного рядаThe ). Данные представлены в порядке убывания среднего расстояния.


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

В противовес расхождению, метрики разделения точек основаны на длине. Интуитивно понятно, что распределения с хорошим разделением точек обычно обладают низким расхождением, и наоборот. Поэтому на практике эти понятия очень близки и часто (произвольно и некорректно) используются как взаимозаменяемые. Однако иногда могут возникать тончайшие различия между этими отличающимися понятиями. Оптимальный ряд точек с низким расхождением может иметь хорошее, но неоптимальное разделение точек. И наоборот, ряды с идеальным разделением точек могут иметь хорошее, но неидеальное расхождение. Иными словами, «дисковые сэмплы best candidate имеют хорошее распределение пространства между точками, но неравномерное распределение в целом» [Christensen, 2018]. Это может помочь нам в более глубоком интуитивном понимании того, почему ряд Соболя имеет превосходные скорости схождения, но только неплохое разделение точек; и почему аппроксимированное дисковое сэмплирование может иметь почти оптимальные характеристики синего шума, но только неплохие скорости схождения.

Также стоит заметить, что если судить только по метрике минимального/среднего разделения точек, то покажется, что флуктуационные ряды Холтона на практике бесполезны и не отличимы от случайного шума. Однако из рисунка 4 чётко видно, что флуктуационный ряд Холтона демонстрирует меньше скоплений и пробелов, чем случайный шум. Следовательно, каждый ряд имеет свои достоинства и недостатки, а отдельные метрики только дают понимание об узком аспекте общих характеристик. (Ещё одна ключевая метрика, не указанная в этой таблице — время, необходимое для вычисления этих множеств точек. Большинство классических распределений синего шума основаны на методе принятия/отбрасывания или методик наподобие имитации отжига, а поэтому довольно медленны по сравнению с методами, вычисляющими точки напрямую.)

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

Флуктуационные бесконечные ряды точек


Коэффициент $1/\sqrt{N}$ в уравнении (3) означает, что это выражение можно использовать только для конечных множеств точек.

Во многих случаях $N$ заранее неизвестно.

В ситуациях, когда известен аппроксимированный интервал $N$ (например, $100 < N <1000$), на удивление хорошие результаты можно получить, присвоив размеру флуктуации фиксированное значение $N_{\textrm{avg}}$ в пределах этого интервала. В нашем примере хорошо подойдёт значение $N_{\textrm{avg}} \simeq 400$.

Однако такой аппроксимированный метод не полностью устраивает присутствующих среди нас пуристов! Если мы хотим подвергнуть флуктуации бесконечный ряд, то нам нужно полностью избавиться от этой зависимости от $N$. К счастью, оказывается, что существует простая подстановка, позволяющая нам выполнять флуктуацию бесконечного ряда с неизменно хорошими результатами. Требуемая подстановка выглядит так:

$\frac{1}{\sqrt{N}} \mapsto\frac{1}{2\sqrt{i-i_0}}; \quad \textrm{where } i_0 = 0.700$


То есть размер флуктуации для $i$-той точки теперь является строго убывающей функцией от $i$.

Следовательно, уравнение для флуктуации бесконечного ряда $R_2$ можно упростить:

$\tilde{\pmb{p}_i} (x,y) \rightarrow \pmb{p}_i(x,y) + \epsilon \left( \frac{k}{2 \sqrt{i-i_0}}\right ) ; \quad \textrm{for } i = 1,2,3,…,N; \quad k \simeq 0.38, \; i_0 = 0.700 \tag{6}$


Радиус флуктуации для каждого из бесконечных псевдослучайных рядов:


Распределение точек этого флуктуационного бесконечного ряда показан на рисунке 1 в самом начале статьи. Более того, как следует из рисунка 8, эта адаптация из флуктуации множества точек в флуктуацию бесконечного ряда практически никак визуально не влияет на качество!


Рисунок 8. Сравнение использования флуктуации замкнутого конечного множества точек (сверху) с флуктуацией открытого бесконечного ряда точек (внизу). Верхние ячейки — это три разных множества точек с N=100, N=300 и N=1000. Для сравнения: в ячейках внизу представлены первые 100, 300 и 1000 сэмплируемых точек одного бесконечного ряда точек. Заметьте, что качество соответствующих распределений визуально почти неразличимо.

Вычисление флуктуации


Форма флуктуации

Выше мы обсуждали только радиальную флуктуацию. Эмпирические тесты показывают, что визуальная разница в качестве почти неразличима, если флуктуация просто основана на квадрате с длиной стороны $\sqrt{\pi}r$, а не на окружности радиуса $r$. Хотя это имеет незначительное отрицательное влияние на спектры Фурье и некоторые метрики разделения точек, я не обнаружил измеримой разницы в скоростях схождения. Поэтому ради простоты и/или ускорения вычислений флуктуации можно упростить до

$\pmb{\epsilon}_i(L) = \sqrt{\pi}L (u_1,u_2); \quad (u_1,u_2) \sim U[0,1)^2 \tag{7}$


Размер флуктуаций

В отличие от рандомизации и перемешивания, которые по сути являются процессами вида «всё или ничего», мы можем управлять величиной применяемой флуктуации. Следовательно, мы можем ввести настраиваемый параметр $\lambda$, определяющий степень флуктуации. То есть подвергнутые флуктуации точки можно вычислить так:

$\tilde{\pmb{p}_i} (x,y) \rightarrow \pmb{p}_i(x,y) + \lambda \frac{\delta}{2\sqrt{N}} \epsilon_i(x,y) ; \quad \textrm{for } i = 1,2,3,…,N; \quad \delta_0 \simeq 0.76, \; \lambda >0\; \tag{8}$


На рисунке 9 показан ряд $R_2$ для N= 600, подвергнутый флуктуации на разных уровнях. Заметьте, что $\lambda = 0$ соответствует исходному множеству точек без флуктуаций; при $\lambda = 1$ мы получаем критический уровень флуктуаций; а при $\lambda >2$ мы по сути имеем белый шум. В областях применения, требующих псевдослучайного размещения объектов, «не слишком равномерного и не слишком случайного», наиболее визуально приятными результаты получаются при $\lambda \simeq 50\%$.


Рисунок 9. Изменение относительного размера флуктуаций, от 0% до 100%. Где 100% представляет собой критический порог между исходным дискретным спектром Фурье и белым шумом. В общем случае наиболее визуально приятными оказываются значения в окрестности $\lambda =50\%$. Меньшие значения $\lambda$ означают пониженное расхождение, улучшенное разделение точек и ускорение схождения, но при этом увеличивается количество артефактов алиасинга.

Нестохастические флуктуации


Это последний фрагмент картины. Вас не должно удивлять, что случайные функции, используемые для применения флуктуаций, не обязаны обладать криптографическим качеством. Флуктуационное сэмплирование должно отвечать только довольно слабым критериям, поскольку оно играет вторичную роль в структуре распределения точек. Следовательно, там, где требуется простота и/или высокая скорость вычислений, можно (и нужно?) просто заменить всю рандомизацию самыми любимыми высокоскоростными функциями хеширования, выдающих значения в интервале $[0,1)$.

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

К счастью, на вопрос можно уверенно ответить положительно!

Насколько я знаю, пока нет опубликованных способов явного и детерминированного создания изотропных рядов точек синего шума. Все известные мне способы (для конечных множеств точек или открытых рядов) требуют перемешивания, перестановки, методов принятия/отбрасывания и/или итеративной релаксации.

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

Мы не только можем заменить случайные флуктуации полностью детерминированной функцией, но и структура предлагаемой формулы невероятно близка к структуре самого ряда $R_2$. (Пока вы не спросили — как бы ни хотел я использовать для флуктуаций ряд с низким расхождением, все эмпирические тесты пока демонстрируют, что флуктуации должны быть намного «менее регулярными», так как каждая точка флуктуации должна быть более независимой от предыдущей точки, чем это оказывается в псевдослучайных рядах.)

Общеизвестно, что Вейл продемонстрировал следующее: для любого иррационального числа $\alpha$ дробная часть кратных величин $x$ имеет равномерное распределение. Очевидно, что это образует базис ряда $R_2$. Менее известно то, что он также показал, что для почти любого числа $\beta>1$ дробная часть степеней $\beta$ распределена равномерно (см. fractional part of powers).

То есть вместо выбора случайного числа из интервала $[0,1)$ мы можем выбрать константу $\beta>1$, а затем использовать $i$-тый член этого ряда. В частности, для нашего случая, вместо $u_1$ и $u_2$, являющихся случайными переменными из интервала $[0,1)$, мы зададим $\pmb{u} = (u_1,u_2)$ следующим образом:

$\pmb{u_i} = (\beta_1^i \; \% 1, \beta_2^i \; \%1) \quad \textrm{where } (\beta_1,\beta_2) = \left( \frac{3}{2}, \frac{4}{3} \right). \tag{9}$


Если вы не работаете с ПО, автоматически обрабатывающим произвольные длины integer (например, с Mathematica), то обязательно столкнётесь с проблемами переполнения/исчезновения значащих разрядов для $n>70$, если будете вычислять его непосредственно через выражение frac = np.pow(1.5,n) %1.

Следовательно, несмотря на то, что теоретически $\beta_1,\beta_2$ могут быть любыми вещественными числами больше 1, я выбрал 3/2 и 4/3, потому что:

  • они дают хорошие эмпирические результаты, и
  • они имеют вид $1+1/b$ для целочисленного $b$, степени которого легко вычисляются.

Базовый код (который можно значительно ускорить, если кэшировать промежуточные результаты), реализующий это, представлен в следующем разделе.

Тесты интегрирования QMC


Самые распространённые тесты сходимости, используемые для сравнения QMC — это гауссова кривая (рисунок 7a). Так как это плавная изотропная функция, случайный сэмплер сходится за $(\mathcal{O}(1/n^{0.5})$ и все псевдослучайные ряды с низким расхождением обычно сходятся примерно между $(\mathcal{O}(1/n^{0.75})$ и $(\mathcal{O}(1/n)$, причём ряды Соболя и $R_2$ сходятся немного лучше, чем флуктуационный $R_2$, который в свою очередь чуть лучше Холтона.

QMC с функциями с дискретными разрывами не очень часто используются в финансах и физике, но очень распространены в рендеринге, потому что объекты зачастую имеют резкие границы. Следовательно, тесты сходимости, задействующие функции с разрывами, демонстрируют, насколько хорошо сэмплирование будет рендерить края. То есть, насколько серьёзными будут артефакты алиасинга! Мы ввидим, что в сравнении с гауссианой сэмплер Соболя для треугольника и диска демонстрирует значительное ухудшение времени сходимости. Наиболее заметна плохая сходимость флуктуационного $R_2$ в дискретной одномерной ступенчатой функции. Это является неудачным ограничением флуктуационного $R_2$.

Стоит заметить, что источниками вдохновения для этих эмпирических тестов послужили две превосходные статьи: "Correlated Multi-jittered Sampling" [Kensler, 2013] и "Progressive multi-jittered Sampling" [Christensen, 2018].


Рисунок 10. Сравнение скоростей схождения QMC для различных функций. (i) синяя = $R_2$, (ii) зелёная = флуктуационный $R_2$, (iii) красная = Холтон, (iv) оранжевая = Соболь, и (v) серая = случайная. Также добавлены линии $\mathcal{O}(1/n)$, $\mathcal{O}(1/n^{0.75})$ и $\mathcal{O}(1/n^{0.50})$. Заметьте, что $R_2$ стабильно имеет самые лучше скорости сходимости. Для плавной гауссианы, дискретного треугольники и дискретного диска флуктуационный $R_2$ так же хорош, как Холтон. Для ступенчатой функции флуктуационный $R_2$ проявляет себя плохо, потому что в одномерной ступенчатой функции над его проекциями к несчастью доминирует флуктуация белого шума.

Выводы


Если объединить результаты предыдущих разделов, (уравнения 4, 5, 7, 8), то формула для конечного множества точек параметризированного прямоугольного флуктуационного $R_2$ имеет вид

$\begin{align} \textrm{For } i&=1,2,3,…,N\tag{10} \\ \textrm{define } \tilde{\mathbf{p_i}} &= (\alpha_1 i, \alpha_2 i ) + \frac{\lambda \delta_0 \sqrt{\pi}}{2\sqrt{N}} \left(\beta_1^i, \beta_2^i \right) \;\; \%1 \nonumber \\ & \nonumber \\ \textrm{where } \varphi &= \sqrt[3]{\frac{1}{2} + \frac{\sqrt{69}}{18} } + \sqrt[3]{\frac{1}{2} – \frac{\sqrt{69}}{18} } \simeq 1.324717957244746, \nonumber \\ \alpha &= \left( \frac{1}{\varphi},\frac{1}{\varphi^2} \right) , \; \beta = \left( \frac{3}{2}, \frac{4}{3} \right), \nonumber \\ \delta_0 &\simeq 0.76, \; \lambda >0. \nonumber \\ \end{align}$


Аналогично, если объединить результаты предыдущих разделов (уравнения 4, 6, 7, 8), то формулу бесконечного ряда точек параметризированного прямоугольного флуктуационного $R_2$ можно изящно выразить следующим образом:

$\begin{align} \textrm{For } i&=1,2,3,…\tag{11} \\ \textrm{define } \tilde{\mathbf{p_i}} &= (\alpha_1 i, \alpha_2 i ) + \frac{\lambda \delta_0 \sqrt{\pi}}{4\sqrt{i-i_0}} \left(\beta_1^i, \beta_2^i \right) \;\; \%1 & \nonumber \\ \textrm{where } \varphi &= \sqrt[3]{\frac{1}{2} + \frac{\sqrt{69}}{18} } + \sqrt[3]{\frac{1}{2} – \frac{\sqrt{69}}{18} } \simeq 1.324717957244746, \nonumber \\ \alpha &= \left( \frac{1}{\varphi},\frac{1}{\varphi^2} \right) , \; \beta = \left( \frac{3}{2}, \frac{4}{3} \right), \nonumber \\ \textrm{and } \delta_0 &\simeq 0.76, \; i_0 = 0.700, \lambda >0. \nonumber \\ \end{align}$


Простая реализация этой формулы на GLSL находится здесь [Vassvik, 2018].

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

  • настолько прост, что его можно выразить менее чем в 50 строках простого кода;
  • основан на ряде с низким расхождением ($R_2$)
  • может быть полностью детерминированным или стохастическим, в зависимости от требований. ($\epsilon$)
  • может создавать конечные множества точек или бесконечные ряды точек.
  • способен параметризироваться и изменяться в пределах от ряда с низким расхождением до простого белого шума.($\lambda_0$)

Модифицирование бесконечных псевдослучайных рядов с низким расхождением так, чтобы они демонстрировали более желательные спектры синего шума — совершенно новая область исследований. Превосходные современные методы и исследования можно посмотреть в [Christensen, 2018] и [Perrier, 2018]. Perrier описывает, как модифицировать ряды Соболя при помощи перестановок, чтобы получать очень хорошо распределённые ряды при условии, что $N$ является степенью 16, а Christensen использует стратифицированную флуктуацию. Стоит заметить, что в отличие от метода Perrier, наш метод не требует особых значений $N$. Более того, наш метод более схож с методом Christensen, также основанным на флуктуациях. Нет никаких сомнений, что описанные Perrier и Christensen методы обеспечивают лучшие результаты рендеринга, чем мой, но мой новый подход способен обеспечить в определённых обстоятельствах некоторые преимущества, что в дальнейшем может способствовать развитию современных методов сэмплирования.

Одной из важнейших областей дальнейшей работы является исследование проекций флуктуационного $R_2$ на более низкие размерности. Хотя проекции ряда $R_2$ без флуктуаций обладают низким расхождением, используемый сейчас флуктуационный процесс совершенно маскирует это из-за доминирующего влияния белого шума. Надеюсь, что более изощрённый способ реализации флуктуаций может привести к созданию лучших проекций. Работы Christensen и других исследователей дают понять, что результативным путём может стать стратификация. Ещё одной областью дальнейших исследований является расширение флуктуационного ряда на более высокие размерности. Несмотря на то, что $R_d$ естественным и почти тривиальным способом расширяется до более высоких размерностей, это является более сложной задачей, потому что проекции $R_d$ (как и у большинства псевдослучайных рядов) не имеют хорошего распределения, не говоря уже о ситуации с флуктуационными версиями! Третьей областью дальнейшей работы является определение того, как похожая флуктуационная техника может использоваться для устранения артефактов в псевдослучайной маске дизеринга, описанной в конце моей предыдущей статьи.

Несмотря на очевидный факт, что ещё многое можно улучшить, моя новая техника обеспечивает инновации, применимые практически ко всем псевдослучайным рядам с низким расхождением (с рядом $R_2$ в качестве примера):

  1. метод применения флуктуаций (постепенно снижающейся величины) к псевдослучайному ряду с низким расхождением таким образом, чтобы получившийся ряд демонстрировал изотропные и аппроксимативные характеристики синего шума и приводил к минимальным артефактам рендеринга;
  2. полностью детерминированный метод (без необходимости рандомизации, перемешивания, итераций, имитации отжига или линейного программирования) на основании псевдослучайного множества или ряда точек с низким расхождением, демонстрирующий изотропные и аппроксимативные спектры синего шума и приводил к минимальным артефактам рендеринга; и
  3. метод построения параметризированного ряда, который может плавно варьироваться от чистого псевдослучайного ряда до почти синего шума, а затем к простому белому шуму.

Пример кода


Задача этого кода — помочь в понимании математики, содержащейся в этом посте. При его написании основной целью была читаемость даже для пользователей, незнакомых с Python, поэтому он не так оптимизирован по скорости. Заметьте, что в посте я использую условность, что отсчёт индексов начинается с 1. Однако в коде на Python используется стандартная нумерация, то есть индексы массивов начинаиются с нуля!

import numpy as np

# Returns a pair of deterministic pseudo-random numbers
# based on seed i=0,1,2,...
def getU(i):
  useRadial = True # user-set parameter
  # Returns the fractional part of (1+1/x)^y
  def fractionalPowerX(x,y):
    n = x*y
    a = np.zeros(n).astype(int)
    s = np.zeros(n).astype(int)
    a[0] = 1
    for j in range(y):
      c = np.zeros(n).astype(int)
      s[0] = a[0]
      for i in range(n-1):
        z = a[i+1] + a[i] + c[i]
        s[i+1] = z%x
        c[i+1] = z/x #integer division!
      a = np.copy(s)
    f =0;
    for i in range(y):
      f += a[i]*pow(x,i-y)
    return f
  #
  u = np.zeros(2)
  v = np.zeros(2)
  v = [ fractionalPowerX(2,i+1), fractionalPowerX(3,i+1)]  
  if useRadial:
    u = [pow(v[0],0.5)*np.cos(2*np.pi * v[1]), pow(v[0],0.5)*np.sin(2*np.pi * v[1])] 
  else:
    u = [v[0], v[1] ]
  return u  
return [ fractionalPowerX(2,i+1), fractionalPowerX(3,i+1)]

# Returns the i-th term of the canonical R2 sequence
# for i = 0,1,2,...
def r2(i):
  g = 1.324717957244746 # special math constant
  a = [1.0/g, 1.0/(g*g) ]
  return [a[0]*(i+1) %1, a[1]*(1+i) %1]

# Returns the i-th term of the jittered R2 (infinite) sequence.
# for i = 0,1,2,...
def jitteredR2(i):
  lambd = 1.0 # user-set parameter
  useRandomJitter = False # user parameter option
  delta = 0.76 # empirically optimized parameter
  i0 = 0.300 # special math constant
  p = np.zeros(2)
  u = np.zeros(2)
  p = r2(i)
  if useRandomJitter: 
     u = np.random.random(2)
  else:
    u = getU(i)
  k = lambd*delta*pow(np.pi,0.5)/(4* pow(i+i0,0.5)) # only this line needs to be modified for point sequences
  j = [k*x for x in u]             # multiply array x by scalar k
  pts= [sum(x) for x in zip(p, j)] # element-wise addition of arrays p and j
  return [s% 1 for s in pts]  # element-wise %1 for s

\begin{array}{l}
\textrm{r2}: & [0.7548, 0.5698], [0.5097, 0.1396], [0.2646, 0.7095], [0.0195, 0.2793], [0.7743, 0.8492],…\\
\textrm{u}: &[0.5000, 0.3333], [0.2500, 0.7777], [0.3750, 0.3703], [0.0625, 0.1604], [0.5937, 0.2139],…\\
\text{jitteredR2}: & [0.0623, 0.7747], [0.5835, 0.3694], [0.3479, 0.7917], [0.0310, 0.3091], [0.8708, 0.8839] ,…\\
\end{array}
Tags:
Hubs:
If this publication inspired you and you want to support the author, do not hesitate to click on the button
+19
Comments1
1749
Karma
6.8
Rating
@PatientZero

Переводчик-фрилансер

Send message

Articles

Show the best of all time
Unlock dark mode

Your account

Sections

Information

Services

Support
© 2006–2024, Habr