Примечание: полный исходный код проекта, а также пояснения о его использовании и чтении можно найти на Github [здесь].

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

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

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

Результаты рендерятся при помощи урезанной версии моего движка Homebrew OpenGl Engine, который я модифицировал для рендеринга 2D-массива точек в качестве карты высот. Урезанную версию движка намного проще понять, если вас интересует изучение OpenGL на C++.

Реализация

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

Процесс эрозии на основе частиц выполняется очень просто:

• Мы создаём частицу в случайной точке поверхности
• Она перемещается/скользит по поверхности, используя стандартные классические механики (о них будет сказано ниже)
• Выполняем перенос вещества/осадочных пород между поверхностью и частицей (это тоже объясняется ниже)
• Испаряем часть частицы
• Если частица находится за пределами карты или слишком мала, то уничтожаем её
• Повторяем процесс с нужным вам количеством частиц.

Примечание: параметры системы объяснены в соответствующем разделе. Самым важным является фактор шага времени dt, который пропорционально масштабирует все параметры. Это позволяет нам увеличить частоту симуляции (ценой повышения шума), не меняя относительный масштаб параметров. Это можно увидеть в показанном ниже коде.

Частицы

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

struct Particle{
  //Construct particle at position _pos
  Particle(glm::vec2 _pos){ pos = _pos; }

  glm::vec2 pos;
  glm::vec2 speed = glm::vec2(0.0); //Initialize to 0

  float volume = 1.0;   //Total particle volume
  float sediment = 0.0; //Fraction of volume that is sediment!
};

Примечание: на случай если вы не знакомы с библиотекой GLM: я использую её для выполнения векторных операций.

Частица имеет позицию и скорость, определяющую то, как она перемещается. Кроме того, она имеет объём и долю, определяющую, какая часть объёма является осадочными породами.

Движение: классическая механика

Движение частицы по поверхности симулируется при помощи классической механики. Если вкратце, позиция x частицы изменяется скоростью v, изменяемой ускорением a.

Примечание: выделенные жирным буквы означают, что величина является вектором.

Ещё мы знаем, что сила равна массе, умноженной на ускорение:


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

Следовательно, мы можем сказать, что ускорение a пропорционально вектору нормали поверхности n, разделённому на массу частицы.


где k — константа пропорциональности, а m — масса частицы. Если масса равна объёму, умноженному на плотность, то полную систему движения частицы мы получим с помощью классической механики:

//... particle "drop" was spawned above at random position

glm::ivec2 ipos = drop.pos; //Floored Droplet Initial Position
glm::vec3 n = surfaceNormal(ipos.x, ipos.y);  //Surface Normal

//Accelerate particle using classical mechanics
drop.speed += dt*glm::vec2(n.x, n.z)/(drop.volume*density);
drop.pos   += dt*drop.speed;
drop.speed *= (1.0-dt*friction);  //Friction Factor

//...

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

Процесс образования осадочных пород: массообмен

Процесс образования осадочных пород физически происходит как перенос осадочных пород с земли на частицу и обратно в точке расположения частицы ("массообмен").

В химических технологиях перенос массы (т.е. изменение массы/осадочных пород с течением времени) между двумя фазами (в данном случае это поверхность земли и капля) обычно описывается при помощи коэффициентов массообмена.

Массообмен пропорционален разности между концентрацией c и равновесной концентрацией c_eq:


где k — константа пропорциональности (коэффициент массообмена). Эту разность между равновесной и действительной концентрацией часто называют «движущей силой».

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

Коэффициент массообмена можно интерпретировать различными способами:

• Как частоту перехода между фазами (здесь это «скорость отложения»)
• Скорость, с которой концентрация капли стремится к равновесной концентрации

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

//...

//Compute Equilibrium Sediment Content
float c_eq = drop.volume*glm::length(drop.speed)*(heightmap[ipos.x][ipos.y]-heightmap[(int)drop.pos.x][(int)drop.pos.y]);

if(c_eq < 0.0) c_eq = 0.0;

//Compute Capacity Difference ("Driving Force")
float cdiff = c_eq - drop.sediment;

//Perform the Mass Transfer!
drop.sediment += dt*depositionRate*cdiff;
heightmap[ipos.x][ipos.y] -= dt*drop.volume*depositionRate*cdiff;

//...

Примечание: изменение концентрации внутри частицы полностью описывается уравнением массообмена. Изменение в карте высот дополнительно умножается на объём частиц, потому что мы меняем его пропорционально не концентрации, а массе (концентрация умножить на объём).

Другие аспекты

Карта высот инициализируется многослойным шумом Перлина со случайным seed.

В конце каждого шага времени частица теряет немного массы в соответствии со скоростью испарения:

//...

drop.volume *= (1.0-dt*evapRate);

//...

Этот процесс повторяется для тысяч частиц, создаваемых в случайных местах и симулируемых по отдельности (в моём случае вычисления производятся последовательно в ЦП).

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

Результаты

Готовый код процесса эрозии составляет примерно 20 строк без комментариев.

Вот сравнение «до и после» для 10 выборок. Симуляция генерирует очень красивые хребты на возвышениях, по сторонам некоторых хребтов откладываются осадочные породы, что приводит к созданию красивых плато.


Выборка из десяти сравнений «до и после». Написанный мной шейдер получает на входе два цвета. Также я реализовал наложение теней, зависящий от расстояния туман и затенение по Фонгу.

Вот ещё десять (разных) выборок просто с результатами:


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

По разному выбирая seed при инициализации карты, можно создавать управляемые разные результаты.

Время симуляции

Время симуляции напрямую связано со сроком жизни частиц и количеством симулируемых частиц.

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

• Трение и инерция: скорость перемещения частиц
• Размер сетки: вероятность того, что частицы выпадут с карты
• Скорость испарения: скорость исчезновения частиц

При параметрах по умолчанию на симуляцию отдельной частицы требуется по 10-100 миллисекунд, что в результате даёт 10-20 секунд на симуляцию 200 000 частиц.

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

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

Оптимизировать код можно снижением затрат на повторное создание меша или повышением скорости одного шага времени для отдельной частицы.

Работа на будущее

Вскоре я вставлю этот код в свой проект Territory в качестве основы для генерации карты высот в генераторе рельефа.

Я уже написал фреймворк для упрощённой симуляции динамики жидкости в климатических системах, но он ещё не готов к публикации. Эта система получает климатические паттерны из карты высот. Это будет темой для будущего поста о физически точных процедурных климатических системах!

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

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

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