Hydrologie procédurale: simulation dynamique des rivières et des lacs

Remarque: le code source complet du projet est publié sur Github [ ici ]. Le référentiel contient également des informations détaillées sur la façon de lire et d'utiliser le code.

Explication: je ne suis pas géologue, j'ai donc créé le système en fonction de mes connaissances.

Concept d'hydrologie

• Migration fluviale et fluviale
• Cascades naturelles
• Formation de Canyon
• Gonflement du sol et de la plaine inondable

• Le relief affecte le mouvement de l'eau
• L'érosion et la sédimentation affectent le terrain

• Dans un cours d' eau en mouvement , l' eau se comporte différemment.
• Dans une piscine debout , l' eau se comporte différemment

Modèle hydrologique simplifié

Remarque: n'oubliez pas qu'ils sont modélisés comme des systèmes 2D.

Cartes de flux et de piscine

Remarque: les cartes de flux et de pool sont des tableaux de la même taille que la carte de hauteur.

L'eau en tant que particule

Cycle hydrologique et interaction cartographique

• Création d'une particule sur le terrain
• (.. ).
• , .
• , .
• , ( ) .
• .

la mise en oeuvre

Classe de particules

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

Algorithme de descente

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

Remarque: j'ai découvert qu'un tel changement de paramètre fonctionne 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;

//...

Remarque: la condition de débordement réinitialise également le volume afin que les particules ne passent pas à l'algorithme d'inondation.

Algorithme d'inondation

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

    //...

Remarque: pour plus de simplicité, un algorithme de remplissage à huit voies est utilisé ici . À l'avenir, il sera possible de le mettre en œuvre plus efficacement.

    //...

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

    //...

Remarque: lorsque le niveau d'eau diminue en raison d'une fuite, j'ai découvert que cela fonctionne mieux avec un faible taux de fuite. De plus, l'élimination d'une partie des roches sédimentaires facilite la mise en œuvre.

    //...

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

Enveloppe d'érosion

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

};

Remarque: la carte de hauteur est initialisée à l'aide du bruit 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--;
}

//...

Remarque: les particules pénètrent rarement dans les piscines et les quittent plus d'une ou deux fois avant de s'évaporer complètement pendant les étapes de descente, mais j'ai ajouté cela juste au cas où.

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

}

Des arbres

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

Remarque: le modèle d'arbre se trouve dans le fichier "végétation.h" et dans la fonction "World :: grow ()".

Autres détails

résultats

Comparaison des effets de rétrécissement

• Érosion à base de particules (érosion de base) qui reçoit des cartes de débit et de pool. Les pools affectent toujours la génération
• L'érosion de base avec des paramètres modifiés par des cartes hydrologiques (en combinaison avec l'érosion)
• Érosion combinée avec des paramètres modifiés par les cartes hydrologiques et affectant l'érosion par les arbres

Impact de la piscine

Remarque: j'ai obtenu plusieurs graines, dans lesquelles trois lacs se sont coulés les uns dans les autres de manière séquentielle, mais je n'ai pas voulu passer beaucoup de temps à trouver la bonne pour cet article. J'ai déjà généré trop de photos et de vidéos.

Vitesse de simulation

Remarque: peu de temps après la publication de cet article, j'ai fait un petit changement dans la méthode de calcul des normales de surface, ce qui a augmenté la vitesse de simulation. L'effet de ceci est très perceptible pendant la simulation, mais il est difficile de le comparer en raison de grandes variations dans le temps de lancement et d'inondation. Selon mes estimations, la vitesse de simulation a doublé.

Plus de belles vidéos

Conclusion et travaux pour l'avenir