О сэмплировании дисков Пуассона

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

На полностью противоположном конце спектра в 2007 году Роберт Бриджсон опубликовал статью на одну страницу, которая набрала почти 1000 цитирований и понимается полностью менее чем за 10 минут. В ней предложено простое решение задачи, часто возникающей в компьютерной графике и симуляциях: случайное размещение объектов так, чтобы они не оказывались слишком близко друг к другу.

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

Алгоритм Бриджсона

Пусть желаемое минимальное расстояние между точками равно r, а работа ведётся в d-мерном пространстве. Алгоритм Бриджсона выглядит так:

1. Разбить пространство на сетку с длиной стороны ячейки r/√d. Это гарантирует, что в каждой ячейке сетки может находиться не более одной точки.
2. Инициализировать список active одной случайной точкой, выбранной равномерно из пространства.
3. Пока active не пуст:
   1. Выбрать элемент p из active равномерно случайным образом.
   2. Равномерно сэмплировать кольцо (annulus) с центром в p, внутренним радиусом r и внешним радиусом 2r, не более k раз. Если найдена валидная точка (используя сетку для эффективной проверки коллизий), добавить её в active и выбрать новую p. Если за k попыток валидная точка не найдена, удалить p из active. Бриджсон рекомендует k = 30.

Самый простой способ равномерно сэмплировать кольцо — сгенерировать случайный единичный вектор v⃗ ∈ ℝᵈ и число x, выбранное равномерно из интервала [1/2ᵈ, 1), а затем итоговая точка вычисляется как 2rx^(1/d)·v⃗. Несколько лет назад по этой теме было выпущено видео, объясняющее, почему это работает. В двух измерениях выбор единичного вектора эквивалентен выбору угла θ ∈ [0, 2π). В более высоких измерениях можно нормализовать вектор, каждая компонента которого взята из нормального распределения.

Улучшения

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

Начнём с двумерного улучшения. Рассмотрим момент, когда алгоритм размещает точку p, а затем сэмплирует её кольцо, получая новую точку q. Точку p будем называть родителем q. В отношении между этими точками содержится ценная информация. Когда неизбежно придётся сэмплировать кольцо с центром в q, целый диапазон углов можно не рассматривать, потому что точки в этом диапазоне окажутся слишком близко к p.

Визуальная интуиция здесь проста, но перевод её в формулу — утомительное тригонометрическое упражнение. Опуская детали, конус, образованный этими граничными линиями, центрирован по углу α, а его ширина равна 2β, где α — это угол между направлениями от q к p, а β зависит от расстояния между p и q и минимума двух арккосинусов, определяющих пересечение конуса с внутренней и внешней окружностями кольца.

Единственная интересная деталь этой формулы — минимум в выражении для β. Он учитывает то, что либо внутренняя, либо внешняя окружность кольца может ограничивать конус — в зависимости от расстояния между p и q. Минимум нужен, чтобы гарантировать, что пересечение конуса и кольца целиком содержится в окружности. Граничные точки конуса «перескакивают» с внешней окружности на внутреннюю, когда расстояние между точками пересекает значение √3·r.

Реализация этой оптимизации требует только хранения родителя каждой точки. Тогда можно вычислить углы конуса и генерировать угол θ для следующей точки вне этого конуса.

Сравнение показывает, что при использовании родительской оптимизации алгоритм Бриджсона генерирует заметно больше точек за то же число итераций, чем без неё (эксперимент проводился на сетке со стороной ℓ = 100 при r = 1, по 100 испытаний на точку данных).

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

Вместо изменения способа выбора угла к следующей точке, второе улучшение меняет способ выбора расстояния до неё. Рассмотрим распределение расстояний от каждой точки кольца до его центра. Его функция распределения (CDF) пропорциональна x^d на интервале [r, 2r] (объяснение — в упомянутом видео). А что если заменить показатель степени на некоторую константу c, отличную от d? Тогда точки можно смещать ближе к центру или дальше от него. При c ≠ 0 точная CDF задаётся кусочной формулой, равной нулю при x ≤ r, единице при x > 2r и определённым выражением через x^c на промежутке между ними.

Меняя c, можно наблюдать, как меняются CDF и вид 500 случайных точек в кольце. Напомним, что c = 2 даёт равномерное распределение. Ползунок допускает и значение c = 0, хотя формула для F₀(x) в этом случае не определена из-за деления на ноль — она заменяется пределом при c → 0, который выражается через логарифмы по основанию 2.

Чтобы сэмплировать радиус при произвольном значении c, применяется обратное преобразование сэмплирования. При c = 0 радиус равен r·2^x, где x — равномерная случайная величина на интервале [0, 1). В остальных случаях используется 2ry^(1/c), где y — равномерная случайная величина между 1 и 1/2^c. Границы интервала меняются местами в зависимости от знака c.

Тепловая карта показывает влияние c на количество точек, сгенерированных алгоритмом Бриджсона (тот же эксперимент: 100 испытаний на сетке со стороной ℓ = 100 при r = 1). Из неё следует, что стоит устанавливать очень отрицательное значение c или даже брать предел при c, стремящемся к минус бесконечности, — тогда каждая точка окажется на расстоянии ровно r от своего родителя. Хотя это действительно максимизирует число сгенерированных точек и создаёт более плотную упаковку, достигается это за счёт «ощущения случайности» распределения. В крайнем случае c = −∞ появляются артефакты вроде длинных цепочек точек и пробелов там, куда ограниченное расстояние не может дотянуться. Восстановить дерево генерации точек постфактум в этом случае тоже не составляет труда.

Поэтому нужен баланс между максимизацией плотности точек и сохранением случайности. Если зафиксировать некоторое 15 ≤ k ≤ 40, лучше всего работает c = −1.4 − 17/√k при использовании родительской оптимизации. Эта формула выведена эмпирически, чтобы число сгенерированных точек примерно соответствовало ожидаемому результату равномерного и максимального сэмплера дисков Пуассона (подробнее об этом — ниже). Для каждого k от 15 до 40 бинарным поиском находилось значение c, при котором окружности радиуса r/2 вокруг каждой точки покрывали бы 54,7% всей площади. Этот процент — насыщенное покрытие круглыми дисками в модели случайной последовательной адсорбции. Разумеется, это применимо только в двух измерениях, и в более высоких размерностях c потребует другой настройки.

Стипплинг

До сих пор r оставался константой, но это не обязательно. Минимальное расстояние между точками можно задавать динамически через функцию r: ℝᵈ → ℝ. Тогда после размещения точки p кольцо сэмплируется с внутренним радиусом r(p). Забавное применение этого — задать r как яркость каждого пикселя изображения, получив эффект стипплинга (точечной графики). Чёрно-белый вариант делает это напрямую, а цветной комбинирует три набора точек Пуассона — по одному на цветовой канал.

Алгоритм Бриджсона по своей природе последовательный, но существуют и другие, рассчитанные на параллельное выполнение и дающие значительный прирост производительности. Из них особенно выделяется PixelPie, который полностью выполняется на GPU. На основе этого алгоритма был создан Poisson Cam — стамплер видео в реальном времени, использующий сэмплирование дисков Пуассона. Эта работа над проектом оказалась особенно интересной, поскольку заодно пришлось освоить программирование шейдеров, Rust, алгоритмы stream compaction и, конечно, сам алгоритм PixelPie.

Максимальность и равномерность

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

Прежде чем описывать сам алгоритм, стоит выделить три момента, которые делают его лучше алгоритма Бриджсона:

1. Максимальность. После завершения работы алгоритма гарантированно невозможно разместить ещё одну точку, не нарушив свойство диска Пуассона.
2. Равномерность. Алгоритм сэмплирует равномерное распределение по всем максимальным наборам точек дисков Пуассона.
3. Детерминированность. Алгоритм не полагается на отбраковку — неудачных попыток разместить точку просто не бывает.

Максимальность и равномерность достигались и раньше многими алгоритмами, самый популярный из которых — иерархическое бросание дротиков (hierarchical dart throwing), но алгоритм Митчелла первым добивается этого без какой-либо отбраковки. Он к тому же весьма производителен: строгого бенчмаркинга не проводилось, но реализация Митчелла работает примерно с той же скоростью и генерирует примерно то же количество точек, что и реализация алгоритма Бриджсона с родительской оптимизацией при k = 20 и c = −5.2. Единственный заметный минус этого алгоритма — существенно большая сложность. Упрощённое описание выглядит так:

1. Разбить пространство на сетку с длиной стороны ячейки r/√d.
2. Пока есть место для размещения ещё одной точки:
   1. Случайно выбрать ячейку c из сетки с весом, пропорциональным оставшейся площади.
   2. Разложить c на непересекающиеся треугольники и «клинья» (chocks)¹.
   3. Выбрать случайный треугольник или клин t с весом, пропорциональным площади.
   4. Равномерно сэмплировать точку p из t и добавить её в итоговый набор.
   5. Вырезать из сетки круг радиуса r с центром в p.

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