Процедурная гидрология: динамическая симуляция рек и озёр

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

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

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

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

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

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

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


Эта система способна быстро генерировать очень реалистично выглядящий рельеф с гидрологией. Это видео было отрендерено в реальном времени. Система способна генерировать бесконечное множество таких ландшафтов

Пояснение: я не геолог, поэтому создавал систему на основе имеющихся у меня знаний.

Концепция гидрологии

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

• Миграция рек и потоков
• Естественные водопады
• Образование каньонов
• Набухание почвы и поймы

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

• Рельеф влияет на движение воды
• Эрозия и осаждение пород влияют на рельеф

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

• В движущемся потоке вода ведёт себя иначе
• В стоячем бассейне вода ведёт себя иначе

Примечание: я часто буду упоминать потоки (streams) и бассейны (pools). В модели делается допущение, что это двухмерные явления большого масштаба. Они очень сильно снижают сложность модели.

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

Упрощённая гидрологическая модель

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

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

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

Карты потоков и бассейнов

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

Карта бассейнов описывает неподвижную воду на поверхности (лужи, пруды, озёра, океаны). Она хранит глубину воды в соответствующей позиции карты.

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

Рельеф с гидрологическим воздействием. Слой воды в рендере взят из гидрологических карт.


Скомбинированные гидрологические карты. Светло-синяя — карта потоков, тёмно-синяя — карта бассейнов.

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

Примечание: карта потоков отображается после пропускания через функцию преобразования (ease-in-out) и рендерится на рельефе на основании этого. Для получения более тонких/чётких потоков (или для игнорирования нижних значений на широких плоских областях), это отображение можно модифицировать или задать ему пороговые значения.

Вода как частица

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

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

Примечание: понятие частицы более подробно объяснено в предыдущем посте [перевод на Хабре] (и бесконечном количестве других ресурсов).

Гидрологический цикл и взаимодействие карт

Карты взаимодействуют друг с другом через гидрологический цикл. Гидрологический цикл состоит из следующих этапов:

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

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


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

Реализация

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

Примечание: я буду показывать только относящиеся к делу части кода. Более подробную информацию можно найти в репозитории на Github. Все соответствующие части кода находятся в файле «water.h».

Класс частиц

Структура частицы Drop идентична структуре из предыдущей системы. Descend и flood теперь являются членами структуры, потому что они действуют одновременно только на одну частицу.

struct Drop{
  
  //... constructors

  int index;                         //Flat Array Index
  glm::vec2 pos;                     //2D Position
  glm::vec2 speed = glm::vec2(0.0);
  double volume = 1.0;
  double sediment = 0.0;

  //... parameters
  const double volumeFactor = 100.0; //"Water Deposition Rate"

  //Hydrological Cycle Functions
  void descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale);
  void flood(double* h, double* pool, glm::ivec2 dim);
};

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

Алгоритм спуска

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

void Drop::descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale){

  glm::ivec2 ipos; 

  while(volume > minVol){

    ipos = pos; //Initial Position
    int ind = ipos.x*dim.y+ipos.y; //Flat Array Index

    //Register Position
    track[ind] = true;

    //...
  }
};

Набор параметров модифицируется картами потоков и бассейнов:

//...
  //Effective Parameter Set
  double effD = depositionRate;
  double effF = friction*(1.0-0.5*stream[ind]);
  double effR = evapRate*(1.0-0.2*stream[ind]);
//...

Примечание: я выяснил, что такое изменение параметров хорошо работает.

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

//... nind is the next position after moving the particle
  
  //Out-Of-Bounds
  if(!glm::all(glm::greaterThanEqual(pos, glm::vec2(0))) ||
     !glm::all(glm::lessThan((glm::ivec2)pos, dim))){
       volume = 0.0;
       break;
  }

  //Slow-Down
  if(stream[nind] > 0.5 && length(acc) < 0.01)
    break;

  //Enter Pool
  if(pool[nind] > 0.0)
    break;

//...

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

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

Алгоритм затопления

Частица с оставшимся объёмом может выполнять затопление из своей текущей позиции. Это происходит, если она прекращает спускаться (отсутствует ускорение) или попадает в уже существующий бассейн.

Алгоритм затопления переносит объём частицы в повысившийся уровень воды, изменяя карту бассейна. Техника заключается в инкрементном повышении уровня воды на доли объёма частицы при помощи «проверочной плоскости». При повышении уровня воды объём частицы уменьшается.


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

На каждом шаге мы выполняем заливку (flood-fill) из позиции частицы (т.е. рекурсивно проверяем позиции соседей), добавляя все позиции, расположенные выше исходной плоскости (текущего уровня воды) и ниже проверочной плоскости во «множество затопления». Это та область рельефа, которая является частью бассейна.

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

void Drop::flood(double* height, double* pool, glm::ivec2 dim){

  index = (int)pos.x*dim.y + (int)pos.y;
  double plane = height[index] + pool[index];  //Testing Plane
  double initialplane = plane;                 //Water Level

  //Flood Set
  std::vector<int> set;
  int fail = 10; //Just in case...

  //Iterate while particle still has volume
  while(volume > minVol && fail){

    set.clear();
    bool tried[dim.x*dim.y] = {false};

    //Lowest Drain
    int drain;
    bool drainfound = false;

    //Recursive Flood-Fill Function
    std::function<void(int)> fill = [&](int i){

      //Out of Bounds
      if(i/dim.y >= dim.x || i/dim.y < 0) return;
      if(i%dim.y >= dim.y || i%dim.y < 0) return;

      //Position has been tried
      if(tried[i]) return;
      tried[i] = true;

      //Wall / Boundary of the Pool
      if(plane < height[i] + pool[i]) return;

      //Drainage Point
      if(initialplane > height[i] + pool[i]){

        //No Drain yet
        if(!drainfound)
          drain = i;

        //Lower Drain
        else if( pool[drain] + height[drain] < pool[i] + height[i] )
          drain = i;

        drainfound = true;
        return; //No need to flood from here
      }

      //Part of the Pool
      set.push_back(i);
      fill(i+dim.y);    //Fill Neighbors
      fill(i-dim.y);
      fill(i+1);
      fill(i-1);
      fill(i+dim.y+1);  //Diagonals (Improves Drainage)
      fill(i-dim.y-1);
      fill(i+dim.y-1);
      fill(i-dim.y+1);
    };

    //Perform Flood
    fill(index);

    //...

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

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

Если найдена точка утечки, мы перемещаем частицу (и её объём «переполнения») к точке утечки, чтобы она снова могла начать спускаться. Затем уровень воды спускается до высоты точки утечки.

    //...

    //Drainage Point
    if(drainfound){

      //Set the Particle Position
      pos = glm::vec2(drain/dim.y, drain%dim.y);

      //Set the New Waterlevel (Slowly)
      double drainage = 0.001;
      plane = (1.0-drainage)*initialplane + drainage*(height[drain] + pool[drain]);

      //Compute the New Height
      for(auto& s: set) //Iterate over Set
        pool[s] = (plane > height[s])?(plane-height[s]):0.0;

      //Remove some sediment
      sediment *= 0.1;
      break;
    }

    //...

Примечание: при снижении уровня воды вследствие утечки я выяснил, что лучше всего это работает при малом коэффициенте утечки. Кроме того, помогает выполнению устранение части осадочных пород.

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

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

    //...

    //Get Volume under Plane
    double tVol = 0.0;
    for(auto& s: set)
      tVol += volumeFactor*(plane - (height[s]+pool[s]));

    //We can partially fill this volume
    if(tVol <= volume && initialplane < plane){

      //Raise water level to plane height
      for(auto& s: set)
        pool[s] = plane - height[s];

      //Adjust Drop Volume
      volume -= tVol;
      tVol = 0.0;
    }

    //Plane was too high and we couldn't fill it
    else fail--;

    //Adjust Planes
    float approach = 0.5;
    initialplane = (plane > initialplane)?plane:initialplane;
    plane += approach*(volume-tVol)/(double)set.size()/volumeFactor;
  }

  //Couldn't place the volume (for some reason)- so ignore this drop.
  if(fail == 0)
    volume = 0.0;

} //End of Flood Algorithm

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

Обёртка эрозии

Класс мира содержит все три карты в виде обычных массивов:

class World {

public:
  void generate();            //Initialize
  void erode(int cycles);     //Erode with N Particles

  //...

  double heightmap[256*256] = {0.0};
  double waterstream[256*256] = {0.0};
  double waterpool[256*256] = {0.0};

};

Примечание: карта высот инициализируется при помощи шума Перлина.

Каждый шаг гидрологии для отдельной частицы состоит из следующего:

//...

//Spawn Particle
glm::vec2 newpos = glm::vec2(rand()%(int)dim.x, rand()%(int)dim.y);
Drop drop(newpos);

int spill = 5;
while(drop.volume > drop.minVol && spill != 0){

  drop.descend(heightmap, waterstream, waterpool, track, dim, scale);

  if(drop.volume > drop.minVol)
    drop.flood(heightmap, waterpool, dim);

  spill--;
}

//...

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

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

Функция эрозии оборачивает этот код и выполняет шаги гидрологии для N частиц, напрямую изменяя карту потоков:

void World::erode(int N){

  //Track the Movement of all Particles
  bool track[dim.x*dim.y] = {false};

  //Simulate N Particles
  for(int i = 0; i < N; i++){
   
    //... simulate individual particle

  }

  //Update Path
  double lrate = 0.01;  //Adaptation Rate
  for(int i = 0; i < dim.x*dim.y; i++)
    waterstream[i] = (1.0-lrate)*waterstream[i] + lrate*((track[i])?1.0:0.0);

}

Здесь массив track передаётся функции спуска. Я выяснил, что одновременное отслеживание движения нескольких частиц и соответствующие изменения обеспечивают для карты потоков улученные результаты. Скорость адаптации (adaptation rate) определяет, насколько быстро удаляется старая информация.

Деревья

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

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

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

//... descend function
double effD = depositionRate*max(0.0, 1.0-treedensity[ind]);
//...

Деревья погибают, если они находятся в бассейне, или поток под ними слишком мощен. Кроме того, они имеют случайную вероятность гибели.


Благодаря наложению теней и картам нормалей даже очень простые спрайты деревьев делают рельеф красивее.

Примечание: модель дерева можно найти в файле «vegetation.h» и в функции «World::grow()».

Другие подробности

Результаты визуализируются при помощи самодельной обёртки OpenGL, которая выложена [здесь].

Результаты

Частицы могут создаваться на карте согласно любому нужному вам распределению. В своих демонстрациях я создавал их равномерным распределением на картах размером 256×256.

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

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

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


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

Сравнение воздействия сужения потоков

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

Я симулировал одинаковый рельеф три раза:

• Эрозия на основе частиц (базовая эрозия), получающая карты потоков и бассейнов. Бассейны всё равно влияют на генерацию
• Базовая эрозия с параметрами, изменяемыми гидрологическими картами (в сочетании с эрозией)
• Скомбинированная эрозия с параметрами, изменяемыми гидрологическими картами и влияющими на эрозию деревьями

Примечание: это довольно хаотичная система, и тип проявляющейся гидрологии сильно зависит от рельефа. Сложно найти «очень показательный» пример рельефа.


Рендер рельефа базовой системы


Рендер рельефа скомбинированной системы


Рендер рельефа системы с деревьями

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

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

Другие эффекты более заметны при непосредственном наблюдении за рельефом.


Гидрологическая карта базовой системы


Гидрологическая карта скомбинированной системы


Гидрологическая карта системы с деревьями

Примечание: эти результаты были сгенерированы ровно за 60 секунд времени симуляции.

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


Пример записи того, как расположение деревьев может помогать сохранять местоположение потоков. Это тот же рельеф, что и раньше, со всеми активированными эффектами.

Воздействие бассейнов

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


Пример видео образования бассейнов на более грубом исходном рельефе с бОльшим перепадом высот. Верхнее озеро физически расположено над нижним и сливает получаемую воду напрямую в нижнее озеро.

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

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


Спустя ещё одну минуту генерации сами собой появились несколько новых бассейнов.


При более крутом угле разница высот бассейнов заметнее.


Гидрологическая карта чётко показывает, что центральный бассейн сливается в нижний бассейн.

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

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

Скорость симуляции

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

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

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

Все видео в этой статье записаны в реальном времени, то есть в целом симуляция выполняется быстро.

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

Ещё красивые видео

Симуляция гидрологии на неровном вертикальном рельефе.


Немного более гладкий рельеф. Некоторые озёра немного получаются немного колеблющимися.

Симуляция гидрологии на плоском гладком рельефе.

Последующие видео были записаны после улучшения скорости.


Ещё более плоский и плавный рельеф.


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


Более неровный рельеф с формированием рек.

Я могу продолжать так вечно.

Вывод и работа на будущее

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

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

Было бы интересно добавить разные типы почвы в виде карты почв, из которой можно было бы брать параметры эрозии.

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

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

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