Hidrología de procedimiento: simulación dinámica de ríos y lagos.

Nota: el código fuente completo del proyecto se publica en Github [ aquí ]. El repositorio también contiene información detallada sobre cómo leer y usar el código.

Explicación: no soy geólogo, así que creé el sistema basado en mi conocimiento.

Concepto de hidrología

• Migración de ríos y arroyos
• Cascadas naturales
• Formación del cañón
• Hinchazón del suelo y llanuras aluviales.

• El alivio afecta el movimiento del agua.
• La erosión y la sedimentación afectan el terreno.

• En una corriente en movimiento , el agua se comporta de manera diferente.
• En una piscina de pie , el agua se comporta de manera diferente.

Modelo hidrológico simplificado

Nota: no olvide que están modelados como sistemas 2D.

Mapas de corrientes y piscinas

Nota: los mapas de flujo y agrupación son matrices del mismo tamaño que el mapa de altura.

El agua como partícula.

Ciclo hidrológico e interacción cartográfica.

• Creando una partícula en el terreno
• (.. ).
• , .
• , .
• , ( ) .
• .

Implementación

Clase de partícula

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);
};

Algoritmo de descenso

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]);
//...

Nota: descubrí que dicho cambio de parámetro funciona bien.

//... 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;

//...

Nota: la condición de desbordamiento también restablece el volumen para que las partículas no pasen al algoritmo de inundación.

Algoritmo de inundación

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);

    //...

Nota: para simplificar, aquí se usa un algoritmo de relleno de ocho vías . En el futuro, será posible implementarlo de manera más efectiva.

    //...

    //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;
    }

    //...

Nota: cuando el nivel del agua disminuye debido a las fugas, descubrí que esto funciona mejor con una tasa baja de fugas. Además, la eliminación de parte de las rocas sedimentarias ayuda a la implementación.

    //...

    //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

Envoltura de erosión

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};

};

Nota: El mapa de altura se inicializa utilizando el ruido Perlin.

//...

//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--;
}

//...

Nota: las partículas rara vez ingresan a las piscinas y las dejan más de una o dos veces antes de que se evaporen por completo durante las etapas de descenso, pero agregué esto por si acaso.

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);

}

Arboles

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

Nota: el modelo de árbol se puede encontrar en el archivo "vegetacion.h" y en la función "Mundo :: crecer ()".

Otros detalles

resultados

Comparación de los efectos del estrechamiento de los flujos.

• Erosión basada en partículas (erosión básica) que recibe mapas de flujo y piscinas. Las piscinas aún afectan la generación
• Erosión básica con parámetros modificados por mapas hidrológicos (en combinación con erosión)
• Erosión combinada con parámetros modificados por mapas hidrológicos y que afectan la erosión de los árboles.

Impacto de piscina

Nota: obtuve varias semillas, en las que tres lagos fluyeron entre sí secuencialmente, pero no quería perder mucho tiempo buscando la correcta para este artículo. Ya generé demasiadas fotos y videos.

Velocidad de simulación

Nota: poco después de la publicación de este artículo, realicé un pequeño cambio en el método de cálculo de las normales de superficie, lo que aumentó la velocidad de simulación. El efecto de esto es muy notable durante la simulación, pero es difícil compararlo debido a las grandes variaciones en el tiempo que lleva lanzarse e inundar. Según mis estimaciones, la velocidad de simulación se ha duplicado.

Más hermosos videos

Conclusión y trabajo para el futuro.