PatientZero 3 июн в 09:52

Симуляция жизни частиц в браузере на WebGPU

Средний

18 мин

3.8K

ФизикаМатематика*Разработка игр*Обработка изображений*

Туториал

Перевод

Автор оригинала: Nikita Lisitsa

Я люблю физические симуляции, а в особенности симуляции частиц. Обычно я реализую что-то на основе традиционной физики, но недавно наткнулся на забавную нефизическую модель, которая может демонстрировать поведение, напоминающее жизнь.

Я написал на C++ прототип для собственного движка, а потом решил, что будет интересно попробовать запустить его в браузере при помощи WebGPU API. Он заработал на удивление хорошо, позволяя создавать подобные симуляции:

В посте я расскажу, как он устроен внутри.

Модель Particle Life

В Интернете есть достаточно много ресурсов (в основном, это видео), на которых рассказывается об этой модели, например, 1, 2 и 3, а в ещё одном говорится о конкретных формулах, использованных в симуляции. Из-за достаточно претенциозного названия, состоящего из часто встречающихся слов, поиск по «particle life» не вернёт никаких релевантных результатов (за исключением видео на YouTube), поэтому мне не удалось найти первоисточник этой модели. В одном видео утверждается, что источником вдохновения для создания этой модели стали микробиологические исследования, что кажется правдоподобным, учитывая то, насколько органичными и мобильными выглядят симуляции.

Эта модель, несмотря на своё название, никак не связана с конвеевской игрой «Жизнь». Она больше похожа на Particle Lenia, однако правила гораздо проще.

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

Это может казаться странным и противоречащим физике, но

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

Живые существа потребляют пищу и кислород (а иногда и другие химические вещества), превращая их в энергию, которую они могут тратить на перемещение и поиск других потребляемых ресурсов. Асимметричные силы Particle Life имитируют это, пропуская цикл преобразования веществ в энергию и просто добавляют в систему новую энергию.

Кроме того, это позволяет частицам преследовать друг друга — например, в описанном выше сценарии частица Б постоянно преследует частицу А, которая, в свою очередь, убегает от частицы Б, что в какой-то степени имитирует взаимоотношения хищника и жертвы из биологии. Однако в этой модели частицы не создаются и не уничтожаются, так что частица Б никогда не «съест» частицу А.

Для простоты и добавления согласованности в симуляцию частицы обычно разбиваются на множество типов с таблицей, определяющей способ взаимодействия для каждой пары типов частиц. Лично я использовал модель, предложенную в видео:

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

В виде формул линейная зависимость силы от расстояния может быть записана так:

$F(r) = A\cdot \max\left(0, \frac{|r-R|}{R}\right)$

Где A — величина силы (отрицательная в случае отталкивания), R — радиус действия силы а r — расстояние между двумя частицами.

Вот график силы как функции от расстояния, где A=5 и R=10:

А вот соответствующая потенциальная энергия, если вам нравится физика:

В Desmos можно поэкспериментировать с этими значениями.

Если сложить силу столкновения и, например, силу притяжения, то мы получим такой график силы от расстояния:

И такую потенциальную энергию:

Они тоже есть на Desmos.

Как видите, в случае столкновения + притяжения на некоем расстоянии есть локальный минимум, где сила равна нулю. Это будет идеальной точкой для частиц, то есть они по возможности будут придерживаться такого расстояния.

Обратите внимание, что сила (а значит, и энергия) не увеличиваются до бесконечности, когда расстояние близко к нулю, как это происходит с типичным потенциалом $\frac{1}{r}$ в случае гравитации или электромагнетизма. Это значит, что сила столкновения имеет конечную величину и её могут преодолеть другие силы, если они достаточно большие или радиус достаточно мал. Например, на этом скриншоте несколько частиц слились в одну:

Впрочем, это нас устраивает, ведь это создаёт интересные поведения и не позволяет нашим силам слишком сильно увеличиваться (что обычно происходит в симуляциях потенциала с $\frac{1}{r}$ ).

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

Почему WebGPU?

Стоит объяснить, почему я перенёс симуляцию на GPU. Это типичная хорошо распараллеливаемая задача, для которых и предназначены GPU. Но почему именно WebGPU?

Если вы читали мои статьи раньше, то знаете, что я обожаю WebGPU. Я больше десятка лет работал с OpenGL, и меня от него тошнит. Это очень старый, ужасно несогласованный API с кучей легаси-решений, который перестал развиваться примерно восемь лет назад. Однако Vulkan, который должен был стать альтернативным современным API, требует тысячу строк подготовительного кода, просто чтобы показать на экране треугольник, даже без каких-то буферов вершин и атрибутов. Не поймите меня неправильно, Vulkan довольно крут, но у меня просто нет столько свободного времени, чтобы работать с ним в своих хобби-проектах.

Эту нишу заполнил для меня WebGPU. Это современный и чистый API, и в то же время он достаточно лаконичный. Мне кажется, что такой и должна быть графика реального времени. Я уже использовал его в нескольких десктопных проектах (через wgpu-native), в том числе для полного трассировщика лучей Монте-Карло и симулятора воды над рельефом, и применяю его в текущем основном проекте — градостроительном симуляторе средневековой деревни. Мне очень нравится этот API.

Более того, будучи современным API, он поддерживает такие возможности, как compute-шейдеры и атомарные операции, которые, как мы увидим ниже, очень важны для этой симуляции. Насколько я знаю, WebGL их не поддерживает.

Есть и ещё один приятный бонус — он работает в браузере! Или, по крайней мере, в некоторых экспериментальных версиях некоторых браузеров с некоторыми флагами, но ситуация постепенно будет улучшаться.

Цикл симуляции

Итак, как же нам симулировать частицы с помощью WebGPU? Сам цикл симуляции должен быть достаточно простым:

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

По сути, это симплектический метод Эйлера, который легко написать и который обладает удобными теоретическими свойствами.

Я храню частицы в одном большом буфере GPU. Частица — это простая struct такого вида:

struct Particle
{
    // Позиция
    x : f32,
    y : f32,

    // Скорость
    vx : f32,
    vy : f32,

    // Тип
    species : f32,
}

Я не храню позицию и скорость в виде vec2f просто для соответствия требованиям выравнивания, а тип частицы храню как f32, просто чтобы можно было загружать частицы из JavaScript как один большой Float32Array.

Итак, частица — это всего пять float, то есть 20 байтов. Теоретически, мы можем просто выполнять по одному compute-шейдеру на каждую частицу, которая циклически обходит все частицы, вычисляет все силы и обновляет скорость и позицию. Однако вычисление сил быстро становится узким местом: это операция , потому что нам нужно попарно вычислить все силы между всеми частицами. Это сильно ограничивает количество частиц, которые мы можем обработать. Моя первая реализация на CPU могла выдерживать всего 4096 частиц, даже с многопоточностью.

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

@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;

@group(1) @binding(0) var<uniform> simulationOptions : SimulationOptions;

@compute @workgroup_size(64)
fn particleAdvance(@builtin(global_invocation_id) id : vec3u)
{
    // Защита от чтения из записи за концом массива массива
    if (id.x >= arrayLength(&particles)) {
        return;
    }

    let width = simulationOptions.right - simulationOptions.left;
    let height = simulationOptions.top - simulationOptions.bottom;

    var particle = particles[id.x];

    // Применяем трение
    particle.vx *= simulationOptions.friction;
    particle.vy *= simulationOptions.friction;

    // Перемещаем частицу с учётом скорости
    particle.x += particle.vx * simulationOptions.dt;
    particle.y += particle.vy * simulationOptions.dt;

    // Выполняем коллизиии с границами
    if (particle.x < simulationOptions.left) {
        particle.x = simulationOptions.left;
        particle.vx *= -1.0;
    }

    if (particle.x > simulationOptions.right) {
        particle.x = simulationOptions.right;
        particle.vx *= -1.0;
    }

    if (particle.y < simulationOptions.bottom) {
        particle.y = simulationOptions.bottom;
        particle.vy *= -1.0;
    }

    if (particle.y > simulationOptions.top) {
        particle.y = simulationOptions.top;
        particle.vy *= -1.0;
    }

    // Записываем новое состояние частицы
    particles[id.x] = particle;
}

(Мой реальный код также обрабатывает мышь и зацикленные границы, но принцип остаётся тем же).

simulationOptions.friction выглядит примерно как $\exp(-\text{friction}\cdot\Delta t)$ ; причины этого я объяснял в другой статье.

Вычисляем взаимодействия

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

Чтобы использовать это, мы применим типичное решение с пространственным хэшированием/биннингом: создадим квадратную сетку из ячеек/бинов размером 32×32, рассортируем частицы по этим бинам, а затем будем вычислять взаимодействия только между соседними бинами. Звучит просто, и на CPU это реализуется практически тривиально (например, на C++ можно создать 2D-массив, хранящий std::vector<particle>). Однако для создания подобных структур данных на GPU требуется много времени — например, мы не можем распределить память GPU непосредственно из шейдера!

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

Для этого мы воспользуемся типичным трёхэтапным процессом:

Этап 1: вычисляем количество частиц в каждом бине
Этап 2: параллельно считаем частичную сумму для вычисления смещения каждого бина
Этап 3: помещаем каждую частицу в соответствующий бин

На этапах 1 и 3 активно используются атомарные операции шейдеров. Этап 1 достаточно прост — выполняем compute-шейдер для каждой частицы, вычисляем индекс бина (линеаризованный) для этой частицы и увеличиваем размер этого бина:

@group(0) @binding(0) var<storage, read> particles : array<Particle>;

@group(1) @binding(0) var<uniform> simulationOptions : SimulationOptions;

@group(2) @binding(0) var<storage, read_write> binSize : array<atomic<u32>>;

@compute @workgroup_size(64)
fn fillBinSize(@builtin(global_invocation_id) id : vec3u)
{
    if (id.x >= arrayLength(&particles)) {
        return;
    }

    // Считываем данные частицы
    let particle = particles[id.x];

    // Вычисляем линеаризованный индекс бина
    let binIndex = getBinInfo(vec2f(particle.x, particle.y),
        simulationOptions).binIndex;

    // Увеличиваем размер бина
    atomicAdd(&binSize[binIndex + 1], 1u);
}

Так как это выполняется параллельно множеством вызовов шейдеров, здесь нужно использовать атомарные операции. Я выполняю инкремент массива по индексу binIndex + 1, оставляя нулевой элемент со значением 0, потому чтоблагодаря этому binOffset[i + 1] - binOffset[i] всегда даёт нам размер i-того бина. Стоит учитывать, что этом случае размер массива binSize на 1 больше, чем реальное количество бинов.

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

fn getBinInfo(position : vec2f, simulationOptions : SimulationOptions) -> BinInfo
{
    let binSize = simulationOptions.binSize;
    let gridSize = simulationOptions.gridSize;

    let binId = vec2i(
        clamp(i32(floor((position.x - simulationOptions.left) / binSize)), 0, gridSize.x - 1),
        clamp(i32(floor((position.y - simulationOptions.bottom) / binSize)), 0, gridSize.y - 1)
    );

    let binIndex = binId.y * gridSize.x + binId.x;

    return BinInfo(binId, binIndex);
}

(В реальном коде gridSize не передаётся как uniform, а вычисляется на месте, потому что сначала я написал код так, а потом не стал менять его).

На этапе 2 нам нужно преобразовать массив размеров бинов в массив смещений бинов. Так как частицы из одного бина в памяти расположены по соседству, бин с каким-то индексом начинается, когда мы уже поместили в массив все частицы из бинов с меньшими индексами. Это значит, что смещение O[i] i-того бина оказывается просто суммой размеров S[j] всех предшествующих бинов: $O[i] = \sum\limits_{j=0}^{i-1} S[j]$ . В псевдокоде это выглядит так:

offset[0] = 0;
for (int i = 1; i < binCount; ++i)
    binOffset[i] = binOffset[i - 1] + binSize[i];

Это называется вычислением частичных сумм (prefix sums). Разумеется, это не параллельный алгоритм. Чтобы использовать всю мощь GPU, нужно применить для этих вычислений параллельный алгоритм. Это называется задачей параллельных частичных сумм; о ней мы поговорим в следующем разделе.

Пока допустим, что мы каким-то образом уже вычислили эти смещения. На этапе 3 нам нужно поместить («отсортировать») частицы в новый массив в зависимости о того, к какому бину принадлежит частица. Это можно сделать, снова запустив compute-шейдер для каждой частицы. Так как мы хотим поместить разные частицы в разные индексы нового массива, нужно снова воспользоваться атомарными операциями:

@group(0) @binding(0) var<storage, read> source : array<Particle>;
@group(0) @binding(1) var<storage, read_write> destination : array<Particle>;
@group(0) @binding(2) var<storage, read> binOffset : array<u32>;
@group(0) @binding(3) var<storage, read_write> binSize : array<atomic<u32>>;

@group(1) @binding(0) var<uniform> simulationOptions : SimulationOptions;

@compute @workgroup_size(64)
fn sortParticles(@builtin(global_invocation_id) id : vec3u)
{
    if (id.x >= arrayLength(&source)) {
        return;
    }

    // Считываем данные частицы
    let particle = particles[id.x];

    // Вычисляем линеаризованный индекс бина
    let binIndex = getBinInfo(vec2f(particle.x, particle.y),
        simulationOptions).binIndex;

    // Атомарно вычисляем индекс этой частицы
    // в соответствующем бине
    let indexInBin = atomicAdd(&binSize[binIndex], 1);

    // Записываем частицу в отсортированный массив
    let newParticleIndex = binOffset[binIndex] + indexInBin;
    destination[newParticleIndex] = particle;
}

Здесь перед запуском шейдера нам нужно обнулить массив binSize. Повторюсь, поскольку множество шейдеров выполняется параллельно, для вычисления индекса частицы в бине мы используем атомарные операции.

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

@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@group(0) @binding(1) var<storage, read> binOffset : array<u32>;
@group(0) @binding(2) var<storage, read> forces : array<Force>;

@group(1) @binding(0) var<uniform> simulationOptions : SimulationOptions;

@compute @workgroup_size(64)
fn computeForces(@builtin(global_invocation_id) id : vec3u)
{
    if (id.x >= arrayLength(&particles)) {
        return;
    }

    // Считываем данные частицы
    var particle = particles[id.x];
    let type = u32(particle.type);

    // Вычисляем бин, содержащий эту частицу
    let binInfo = getBinInfo(vec2f(particle.x, particle.y), simulationOptions);

    // Вычисляем диапазон соседних бинов для итераций
    var binXMin = binInfo.binId.x - 1;
    var binYMin = binInfo.binId.y - 1;
    var binXMax = binInfo.binId.x + 1;
    var binYMax = binInfo.binId.y + 1;

    // Учитываем границы сетки
    binXMin = max(0, binXMin);
    binYMin = max(0, binYMin);
    binXMax = min(binInfo.gridSize.x - 1, binXMax);
    binYMax = min(binInfo.gridSize.y - 1, binYMax);

    let particlePosition = vec2f(particle.x, particle.y);

    var totalForce = vec2f(0.0, 0.0);

    // Итеративно обходим соседние бины
    for (var binX = binXMin; binX <= binXMax; binX += 1) {
        for (var binY = binYMin; binY <= binYMax; binY += 1) {
            // Вычисляем линеаризованный индекс бина
            let binIndex = binY * binInfo.gridSize.x + binX;

            // Находим интервал частиц из этого бина
            let binStart = binOffset[binIndex];
            let binEnd = binOffset[binIndex + 1];

            // Итеративно обходим частицы из этого бина
            for (var j = binStart; j < binEnd; j += 1) {
                // Предотвращаем взаимодействие частиц с собой
                if (j == id.x) {
                    continue;
                }

                let other = particlesSource[j];
                let otherType = u32(other.type);

                // Применяем силу между частицами

                let force = forces[type * u32(simulationOptions.typeCount) + otherType];

                var r = vec2f(other.x, other.y) - particlePosition;
                let d = length(r);

                if (d > 0.0) {
                    let n = r / d;

                    totalForce += force.strength * max(0.0, 1.0 - d / force.radius) * n;
                    totalForce -= force.collisionStrength * max(0.0, 1.0 - d / force.collisionRadius) * n;
                }
            }
        }
    }

    // Обновляем скорость, считая массу единичной
    particle.vx += totalForce.x * simulationOptions.dt;
    particle.vy += totalForce.y * simulationOptions.dt;

    particles[id.x] = particle;
}

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

forces — это массив размера (где M — количество типов частиц), в котором хранятся параметры сил, действующих между частицами.

Параллельная частичная сумма

Параллельное вычисление частичных сумм — хорошо изученная задача. Я использовал простейшее решение, названное в Википедии Algorithm 1, а в статье GPU Gems 3 — naive parallel scan.

Иллюстрация из Википедии наглядно всё объясняет:

Основная идея заключается в выполнении множества прогонов по всему массиву $\log_2(N)$ раз (где N — размер массива), каждый раз прибавляя x[i-step] к x[i], где step — это степень двойки. В псевдокоде на итерации k мы делаем следующее:

step = (1 << k);
parallel for (int i = step; i < array.size; ++i)
    array[i] += array[i - step];

Итерация начинается с i = step , потому что в противном случае в array[i - step] мы попытаемся считать адрес до начала массива. На самом деле, этот код не будет работать, потому что перезапись array[i] влияет на то, что будут считывать другие вызовы (с другими значениями i), поэтому необходимо использовать два массива: один для чтения, другой для записи. Примерно так:

parallel for (int i = 0; i < input.size; ++i) {
    if (i < step)
        output[i] = input[i];
    else
        output[i] = input[i] + input[i - step];
}

Для вычисления частичной суммы мы создаём временный массив и играем в пинг-понг входным и временным массивами $\log_2(N)$ раз; при этом step получает значения увеличивающихся степеней двойки. На первой итерации (k = 0, step = 1) массив binSize становится входным, а временный массив — выходным. На следующей итерации (k = 1, step = 2) временный массив входной, а binSize — выходной и так далее. С точки зрения WebGPU, это требует создания двух групп привязок (bind group) с одинаковой структурой: одна выполняет чтение из массива binSize и записывает во временный массив, другая выполняет обратную операцию.

При такой методике получающаяся частичная сумма будет в первом или втором массиве, в зависимости от чётности количества выполненных шагов, то есть чётности $\log_2(N)$ . Это довольно неудобно, поэтому я всегда округляю это до ближайшего чётного числа: благодаря этому результат всегда будет находиться в одном и том же массиве binSize, а лишняя итерация частичных сумм просто будет автоматически копировать один буфер в другой (посмотрите, что произойдёт с кодом, когда step >= array.size).

Объяснение получилось долгим, но реализация достаточно проста:

@group(0) @binding(0) var<storage, read> input : array<u32>;
@group(0) @binding(1) var<storage, read_write> output : array<u32>;
@group(0) @binding(2) var<uniform> stepSize : u32;

@compute @workgroup_size(64)
fn prefixSumStep(@builtin(global_invocation_id) id : vec3u)
{
    if (id.x >= arrayLength(&input)) {
        return;
    }

    if (id.x < stepSize) {
        output[id.x] = input[id.x];
    } else {
        output[id.x] = input[id.x - stepSize] + input[id.x];
    }
}

Чтобы поддерживать различные значения stepSize и одновременно умещать все операции в едином проходе compute, я создал буфер со всеми нужными мне значениями stepSize и использую его как источник для uniform-значения stepSize со всеми динамическими смещениями в вызове GPUComputePassEncoder.setBindGroup.

Итак, при $\log_2(N)$ итераций (округлёнными до чётного числа) мы выполняем этот compute-шейдер (по одному вызову на бин; или, в моём случае, binCount + 1 вызовов по изложенным выше причинам), передавая туда и обратно данные между двумя буферами. Вот, как это выглядит на JavaScript:

prefixSumIterations = Math.ceil(Math.ceil(Math.log2(binCount + 1)) / 2) * 2;

...

binningComputePass.setPipeline(binPrefixSumPipeline);
for (var i = 0; i < prefixSumIterations; ++i) {
    binningComputePass.setBindGroup(0, binPrefixSumBindGroup[i % 2], [i * 256]);
    binningComputePass.dispatchWorkgroups(Math.ceil((binCount + 1) / 64));
}

(Динамическое смещение i * 256 используется из-за требования кратности динамических смещений 256 байтам. Количество рабочих групп вычисляется как Math.ceil((binCount + 1) / 64), потому что шейдер использует рабочую группу размером 64 и я запускаю их для массива размером binCount + 1).

Можно было бы использовать более сложный алгоритм вычисления частичных сумм, но это не требуется: например, даже при 10000 бинов все эти вычисления сортировки частиц и частичных сумм занимают на моём GeForce GTX 1060 около 0,1 мс. Вся производительность всегда тратится на вычисление сил между частицами, так что оптимизация этапа разбиения на бины не даёт никакой выгоды.

Рендеринг

Я рендерю частицы в виде квадратов, которые во фрагментном шейдере преобразуются в круги. Для каждой частицы у нас есть 2 треугольника, то есть 6 или 4 вершины (с дублированием вершин или без него), хранящие данные одной частицы. Это довольно неудобно: мне бы не хотелось дублировать данные частиц и использовать для этого инстансинг. Поэтому мои вершины не имеют никаких атрибутов и непосредственно считывают массив частиц как буфер хранения только для чтения на основании ID вершин:

@group(0) @binding(0) var<storage, read> particles : array<Particle>;
@group(0) @binding(1) var<storage, read> species : array<Species>;
@group(1) @binding(0) var<uniform> camera : Camera;

struct VertexOut
{
    @builtin(position) position : vec4f,
    @location(0) offset : vec2f,
    @location(1) color : vec4f,
}

const offsets = array<vec2f, 6>(
    vec2f(-1.0, -1.0),
    vec2f( 1.0, -1.0),
    vec2f(-1.0,  1.0),
    vec2f(-1.0,  1.0),
    vec2f( 1.0, -1.0),
    vec2f( 1.0,  1.0),
);

@vertex
fn vertexCircle(@builtin(vertex_index) id : u32) -> VertexOut
{
    let particle = particles[id / 6u];
    let offset = offsets[id % 6u] * 1.5;
    let position = vec2f(particle.x, particle.y) + offset;
    return VertexOut(
        vec4f((position - camera.center) / camera.extent, 0.0, 1.0),
        offset,
        species[u32(particle.species)].color
    );
}

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

@fragment
fn fragmentCircle(in : VertexOut) -> @location(0) vec4f
{
    let alpha = clamp(camera.pixelsPerUnit * (1.0 - length(in.offset)) + 0.5, 0.0, 1.0);
    return in.color * vec4f(1.0, 1.0, 1.0, alpha);
}

Вот, как это выглядит:

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

@vertex
fn vertexGlow(@builtin(vertex_index) id : u32) -> VertexOut
{
    let particle = particles[id / 6u];
    let offset = offsets[id % 6u];
    let position = vec2f(particle.x, particle.y) + 12.0 * offset;
    return CircleVertexOut(
        vec4f((position - camera.center) / camera.extent, 0.0, 1.0),
        offset,
        species[u32(particle.species)].color
    );
}

@fragment
fn fragmentGlow(in : VertexOut) -> @location(0) vec4f
{
    let l = length(in.offset);
    let alpha = exp(- 6.0 * l * l) / 64.0;
    return in.color * vec4f(1.0, 1.0, 1.0, alpha);
}

Эффект сияния, усиленный для наглядности

Всё это рендерится при помощи простого аддитивного смешения (а не обычного over blending!) в HDR-текстуру rgba16float, а затем композитингом выводится на экран (на самом деле, на HTML canvas) с использованием тональной компрессии ACES и дизеринга синего шума, чтобы скрыть полосы, возникающие из-за тёмных цветов на краях сияния. Я экспериментировал с различными способами смешения и кривыми тональной компрессии и выбрал наиболее приятное на мой взгляд сочетание.

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

@vertex
fn vertexPoint(@builtin(vertex_index) id : u32) -> VertexOut
{
    let particle = particles[id / 6u];
    let offset = 2.0 * offsets[id % 6u] / camera.pixelsPerUnit;
    let position = vec2f(particle.x, particle.y) + offset;
    return CircleVertexOut(
        vec4f((position - camera.center) / camera.extent, 0.0, 1.0),
        offset,
        species[u32(particle.species)].color
    );
}

Затем фрагментный шейдер вычисляет площадь пересечения круга и текущего пикселя. Так как вычисление пересечений круга и квадрата — достаточно нетривиальная задача, я аппроксимизирую круг при помощи квадрата, вычисляю пересечения (которые сводятся к min/max по каждой из координат) и умножаю на $\frac{\pi}{4}$ , чтобы компенсировать разницу соотношений площадей квадрата и круга:

@fragment
fn fragmentPoint(in : CircleVertexOut) -> @location(0) vec4f
{
    let s = vec2f(camera.pixelsPerUnit);
    let d = max(vec2f(0.0), min(in.offset * s + 0.5, s) - max(in.offset * s - 0.5, - s));
    let alpha = (PI / 4.0) * d.x * d.y;
    return vec4f(in.color.rgb, in.color.a * alpha);
}

Это похоже на хак, но работает: изображение почти не мерцает, а переход между двумя режимами рендеринга практически незаметен.

Галерея

Я уже давал ссылку в начале статьи, но если вам не хочется возвращаться, то к симуляции можно перейти здесь. При каждой перезагрузке страницы создаётся новая случайная система, или же можно просто нажимать на «Randomize». Можно сохранять текущие правила в файл JSON или загружать собственные правила. Поделиться случайно сгенерированным правилом можно, нажав на «Copy link».

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

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