# Generating terrain with hydraulic erosion

**URL:** <https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948>\
**Category:** Gallery\
**Created:** [May 5, 2019, 7:55pm UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948 "2019-05-05T19:55:38Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jb4x](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jb4x/32/789_2.png) [@jb4x](https://discourse.processing.org/u/jb4x)\
**Post date:** [May 5, 2019, 7:55pm UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948/1 "2019-05-05T19:55:38Z")

</div>

Hi,

After watching that video: [https://www.youtube.com/watch?v=eaXk97ujbPQ](https://www.youtube.com/watch?v=eaXk97ujbPQ)

And reading that paper: [https://www.firespark.de/resources/downloads/implementation%20of%20a%20methode%20for%20hydraulic%20erosion.pdf](https://www.firespark.de/resources/downloads/implementation%20of%20a%20methode%20for%20hydraulic%20erosion.pdf)

I decided to take a go at generated terrain using hydraulic erosion.

The code below definitively gives some results although they are not completely as expected.  
Here some before/after pictures:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/2/2698936ca8488b6aa6a8f06c46c48b7bec412472.jpeg)

To play with the code, just hit play, zoom a bit with the middle mouse button to see the terrain and hit space to erode the terrain a bit.

```auto
import peasy.*;

PeasyCam cam;
PRNG prng;
float terrainHeights[][];
PShape terrain;  
ErosionSystem es;

final float scaleFactorXY = 0.1;
final float scaleFactorZ = 50;
final int nbOfX = 512;
final int nbOfY = 512;
//final float xOffset = scaleFactorXY * nbOfX / 2.;
//final float yOffset = scaleFactorXY * nbOfY / 2.;

final float xOffset = 0;
final float yOffset = 0;

void settings() {
  size(1000, 1000, P3D);
}

void setup() {
  cam = new PeasyCam(this, 400);
  perspective(PI/3.0, width/height, 0.1, 200);

  es = new ErosionSystem();
  prng = new PRNG();
  terrainHeights = new float[nbOfX][nbOfY];

  float noiseScale = 0.004;
  for (int x = 0; x < nbOfX; x++) {
    for (int y = 0; y < nbOfY; y++) {
      terrainHeights[x][y] = prng.mixPerlin(noiseScale * x, noiseScale * y, 8, 2, 0.4);
    }
  }

  //createTerrainShape();
}

void draw() {
  lights();
  background(20);
  //shape(terrain);

  for (int x = 0; x < nbOfX - 1; x++) {

    beginShape(TRIANGLE_STRIP);
    fill(150);
    noStroke();

    for (int y = 0; y < nbOfY - 1; y++) {
      vertex((scaleFactorXY * x) - xOffset, (scaleFactorXY * y) - yOffset, scaleFactorZ * terrainHeights[x][y]);
      vertex((scaleFactorXY * (x+1)) - xOffset, (scaleFactorXY * y) - yOffset, scaleFactorZ * terrainHeights[x+1][y]);
    }

    endShape();
  }

  //es.showDrops(scaleFactorXY, scaleFactorZ);
}

void keyPressed() {
  if (key == ' ') {
    for (int i = 0; i < 1; i++) {
      es.erode(terrainHeights, 100000);
    }
    println("Done eroding");
    //createTerrainShape();
  }
}

void createTerrainShape() {
  terrain = createShape();
  for (int x = 0; x < nbOfX - 1; x++) {

    terrain.beginShape(TRIANGLE);
    terrain.fill(150);
    terrain.noStroke();

    for (int y = 0; y < nbOfY - 1; y++) {
      terrain.vertex((scaleFactorXY * x) - xOffset, (scaleFactorXY * y) - yOffset, scaleFactorZ * terrainHeights[x][y]);
      terrain.vertex((scaleFactorXY * (x+1)) - xOffset, (scaleFactorXY * y) - yOffset, scaleFactorZ * terrainHeights[x+1][y]);
      terrain.vertex((scaleFactorXY * x) - xOffset, (scaleFactorXY * (y + 1)) - yOffset, scaleFactorZ * terrainHeights[x][y+1]);

      terrain.vertex((scaleFactorXY * x) - xOffset, (scaleFactorXY * (y + 1)) - yOffset, scaleFactorZ * terrainHeights[x][y+1]);
      terrain.vertex((scaleFactorXY * (x+1)) - xOffset, (scaleFactorXY * y) - yOffset, scaleFactorZ * terrainHeights[x+1][y]);
      terrain.vertex((scaleFactorXY * (x+1)) - xOffset, (scaleFactorXY * (y + 1)) - yOffset, scaleFactorZ * terrainHeights[x+1][y+1]);
    }

    terrain.endShape();
  }

  println("Done creating terrain");
}

class ErosionSystem {
  private float initialDropSpeed;
  private float initialDropWaterVolume;
  private PRNG prng;
  private long prngSeed;
  private float dropsInertia; // Resistance to change direction - 0: no resistance; 1: full resistance
  private int maxDropStep;
  private float erosionRadius;
  private float gravity; // Affect the speed of the drops
  private float evaporationFactor; // Affect the speed of evaporation (1: full evaporation - 2: no evaporation) ---- TO DO: Change that factor based on the map localisation  
  private float depositFactor; // The percentage of sediment dropped when the max sediment capacity is reached ---- TO DO: Change that factor based on the map localisation  
  private float erosionFactor; // The percentage of sediment taken when the sediment quantity of the drop is lower than the max sediment capacity ---- TO DO: Change that factor based on the map localisation  
  private float sedimentCapacityFactor; // Influence the sediment capacity of a drop. The higher the value the higher the drop will be able to carry sediment in the same condition ---- TO DO: Change that factor based on the map localisation  
  private ArrayList<DropPath> dp; // Used for debuging: Draw the path of the drops.
  private boolean debugMode;

  /*
   * Basic constructor
   */
  ErosionSystem() {
    prng = new PRNG();
    init();
  }

  /*
   * Constructor with seed to get predictable results
   */
  ErosionSystem(long seed) {
    prng = new PRNG(seed);
    init();
  }

  /*
   * Initialize all the variables
   */
  void init() {
    dp = new ArrayList<DropPath>();
    debugMode = false;
    prngSeed = prng.getInitSeed();
    dropsInertia = 0.2;
    maxDropStep = 30;
    initialDropSpeed = 1.;
    initialDropWaterVolume = 1.;
    erosionRadius = 6;
    gravity = 4.;
    evaporationFactor = 0.01;
    depositFactor = 0.7;
    erosionFactor = 0.4;
    sedimentCapacityFactor = 4;
  }

  /*
   * This method is used to modify and shape an existing 3D terrain by simulating the erosion of rain
   * It works directly on the given data to avoid useless copy of huge map
   *
   * @param mapToErode A 2D float matrices containing he height of all the vertices of the terrain
   * @param nbOfIteration The number of drops to simulate
   */
  void erode(float[][] mapToErode, int nbOfIteration) {
    int maxSizeX = mapToErode.length;
    int maxSizeY = mapToErode[0].length;
    int nbOfSimulatedSteps = 0;

    for (int k = 0; k < nbOfIteration; k++) {
      if (debugMode) {
        dp.add(new DropPath());
      }

      // Create a new drop at a random position
      Drop drop = new Drop(prng.random(0, maxSizeX - 2), prng.random(0, maxSizeY - 2), initialDropSpeed, initialDropWaterVolume);
      drop.setInertia(dropsInertia);

      while (true) { // The exit conditions will be taken care of with break instructions
        // Get the height and gradients of the drop
        DropInfo dropInfo = getDropInfo(mapToErode, drop);
        if (debugMode) {
          dp.get(dp.size() - 1).add(dropInfo.x, dropInfo.y, dropInfo.h);
        }

        // Change the direction of the drop based on the gradients
        PVector newDir = new PVector(dropInfo.gradX, dropInfo.gradY);
        newDir.normalize();
        drop.move(newDir);

        // Check if the drop is still in the map and if not stop the simulation
        if (dropIsOutOfBound(mapToErode.length, mapToErode[0].length, drop)) {
          break;
        } 

        // Get the new height and gradients of the drop
        DropInfo newDropInfo = getDropInfo(mapToErode, drop);

        // Get Height difference between the old position and the new position
        // Negative value means current position is lower than old position
        float heightDiff = newDropInfo.h - dropInfo.h;

        // Update the sedimenet capcity of the drop
        drop.updateMaxSedimentCapacity(heightDiff, sedimentCapacityFactor);

        // Get the amount of sediment to erode or deposit
        float sedimentAmount = drop.updateSedimentQty(heightDiff, depositFactor, erosionFactor);

        if (sedimentAmount > 0) { // Sediment is deposited
          int i = dropInfo.i;
          int j = dropInfo.j;
          float u = dropInfo.u;
          float v = dropInfo.v;

          mapToErode[i][j] += (1 - u) * (1 - v) * sedimentAmount;
          mapToErode[i + 1][j] += (u ) * (1 - v) * sedimentAmount;
          mapToErode[i][j + 1] += (1 - u) * (v ) * sedimentAmount;
          mapToErode[i + 1][j + 1] += (u ) * (v ) * sedimentAmount;
        } else { // Sediment is taken
          erode(mapToErode, dropInfo.x, dropInfo.y, sedimentAmount, erosionRadius);
        }

        drop.updateSpeed(heightDiff, gravity);
        drop.updateWaterVolume(evaporationFactor);

        // Check if we don't go over the maximum number of steps allowed
        nbOfSimulatedSteps++;
        if (nbOfSimulatedSteps >= maxDropStep) {
          break;
        }
      }
    }
  }

  /*
   * Erode the map around a drop position. The effect of the erosion decrease linearly with the distance from the center
   *
   * @param mapToErode A 2D float matrices containing he height of all the vertices of the terrain
   * @param fromX The x position of the drop (the center of erosion)
   * @param fromY The y position of the drop (the center of erosion)
   * @param amount The amount of sediment to take from the map
   * @param radius The radius around the drop position affected by the erosion
   */
  private void erode(float[][] mapToErode, float fromX, float fromY, float amount, float radius) {
    final int offset = ceil(radius);
    final int maxNbOfVerticesToCheck = (2 * offset + 1) * (2 * offset + 1);
    final float radiusSquare = radius * radius;

    // The following arrays will contain 
    // - vertexCoordinate: the valid coordinates of the effected vertices
    // - vertexWeight: The weights of the affected vertices
    // Since some values can be negative (because the brush is on the side of the map) that array might not be used entirely
    // The realVertexNb variable keep track of how much of the vertices is actually concerned and thus how deep in the previous arrays we need to go
    int vertexCoordinate[][] = new int[maxNbOfVerticesToCheck][2];
    float vertexWeight[] = new float[maxNbOfVerticesToCheck];
    int realVertexNb = 0;
    float totalWeight = 0; // Use to normalize the values

    // Get top left grid node coordinate
    int i = (int) fromX;
    int j = (int) fromY;

    for (int x = i - offset; x <= i + offset; x++) {
      for (int y = j - offset; y <= j + offset; y++) {
        if (x >= 0 && y >= 0 && x < mapToErode.length && y < mapToErode[0].length) { // If the coordinates are valid
          float w = radiusSquare - ( (x - fromX) * (x - fromX) + (y - fromY) * (y - fromY) );

          if (w > 0) { // If the point is within the radius
            vertexCoordinate[realVertexNb][0] = x;
            vertexCoordinate[realVertexNb][1] = y;
            vertexWeight[realVertexNb] = w;
            totalWeight += vertexWeight[realVertexNb];
            realVertexNb++;
          }
        }
      }
    }

    for (int k = 0; k < realVertexNb; k++) {
      int gridX = vertexCoordinate[k][0];
      int gridY = vertexCoordinate[k][1];
      float w = vertexWeight[k];
      mapToErode[gridX][gridY] += amount * (w / totalWeight);
    }
  }

  /*
   * Check if a drop is outside of the map
   *
   * @param sizeX the x dimension of the map
   * @param sizeY the y dimension of the map
   * @param drop the drop to check
   *
   * @Return true if the drop is outside of the terrain, false otherwise
   */
  private boolean dropIsOutOfBound(float sizeX, float sizeY, Drop drop) {
    float x = drop.getPos().x;
    float y = drop.getPos().y;

    return (x < 0 || y <0 || x > sizeX - 2 || y > sizeY - 2);
  }

  /*
   * Change the inertia value of the drop
   * The Drop object is making sure that the value is within the correct range [0,1]
   */
  void setInertia(float value) {
    dropsInertia = value;
  }

  /*
   * Change the maximum of step a drop can do before not beeing simulated anymore
   * Any value less than 1 will still result in at least one simulated step
   */
  void setMaxDropStep(int value) {
    maxDropStep = value;
  }

  /*
   * Change the initial speed of a drop
   * Make sure it is in the range [0,infinity[
   */
  void setInitialDropSpeed(float value) {
    if (value < 0) {
      initialDropSpeed = 0;
    }
    initialDropSpeed = value;
  }

  /*
   * Change the initial water volume of a drop
   * Make sure it is in the range [0,infinity[
   */
  void setInitialDropWaterVolume(float value) {
    if (value < 0) {
      initialDropWaterVolume = 0;
    }
    initialDropWaterVolume = value;
  }

  /*
   * Change the initial water volume of a drop
   * Make sure it is in the range [1,infinity[
   */
  void setErosionRadius(float value) {
    if (value < 1) {
      erosionRadius = 1;
    }
    erosionRadius = value;
  }

  /*
   * Change the gravity of the system
   * Make sure it is in the range [0,infinity[
   */
  void setGravity(float value) {
    if (value < 0) {
      gravity = 0;
    }
    gravity = value;
  }

  /*
   * Change the evaporation factor of the system
   *
   * Make sure it is in the range [0,1]
   */
  void setEvaporationFactor(float value) {
    if (value < 0) {
      evaporationFactor = 0;
    }
    if (value > 1) {
      evaporationFactor = 1;
    }
    evaporationFactor = value;
  }

  /*
   * Change the deposit factor of a drop
   * Make sure that the value is in range [0,1]
   */
  void setDepositFactor(float value) {
    if (value < 0) {
      depositFactor = 0;
      return;
    }

    depositFactor = value;
  }

  /*
   * Change the erosion factor of a drop
   * Make sure that the value is in range [0,1]
   */
  void setErosionFactor(float value) {
    if (value < 0) {
      depositFactor = 0;
      return;
    }

    depositFactor = value;
  }

  /*
   * Change the sediment capacity factor of a drop
   * Make sure that the value is in range [0,infinity[
   */
  void setSedimentCapacityFactor(float value) {
    if (value < 0) {
      sedimentCapacityFactor = 0;
      return;
    }

    sedimentCapacityFactor = value;
  }

  /*
   * Display the drops as sphere
   * 
   * @Param scaleXY The scale to apply to the x and y axis
   * @Param scaleZ The scale to apply to the z axis
   */
  void showDrops(float scaleXY, float scaleZ) {
    for (int i = 0; i < dp.size(); i++) {
      dp.get(i).show(scaleXY, scaleZ);
    }
  }

  /*
   * When debugMode is true the drops are saved
   * 
   * @Param value True to activate debugMode, false otherwise
   */
  void setDebug(boolean value) {
    debugMode = value;
  }
}

class Drop {
  private PVector pos;
  private PVector dir;
  private float inertia;
  private float sedimentQty;
  private float maxSedimentCapacity;
  private float minSedimentCapacity;
  private float speed;
  private float waterVolume;

  Drop(float x, float y, float initialSpeed, float initialWaterVolume) {
    pos = new PVector(x, y);
    dir = new PVector(1, 0);
    inertia = 0.5;
    sedimentQty = 0;
    maxSedimentCapacity = 1.;
    minSedimentCapacity = 0.01;
    speed = initialSpeed;
    waterVolume = initialWaterVolume;
  }

  /*
   * Modify the direction of the drop based on a new target direction
   *
   * @Param newDir a unit vector that points toward the desired direction
   */
  private void changeDir(final PVector newDir) {
    dir = dir.mult(inertia).add(newDir.mult(1 - inertia)).normalize();
  }

  /*
   * Modify the position of the drop based on a new target direction
   *
   * @Param newDir a unit vector that points toward the desired direction
   */
  void move(final PVector newDir) {
    changeDir(newDir);
    pos.add(dir);
  }

  /*
   * Update the max sediment capacity of the drop based on the height difference between its old position and its new position
   *
   * @Param heightDiff The height difference between its old position and its new position
   */
  void updateMaxSedimentCapacity(float heightDiff, float sedimentCapacityFactor) {
    maxSedimentCapacity = - heightDiff * speed * waterVolume * sedimentCapacityFactor;

    if (maxSedimentCapacity < minSedimentCapacity) {
      maxSedimentCapacity = minSedimentCapacity;
    }
  }

  /*
   * Update the sediment quantity of the drop
   * It wil deposit or take sediment depending on the situation
   *
   * @Param heightDiff The height difference between its old position and its new position
   *
   * Return The amount of sediment to take or deposit (return a negative value when sediment is taken and a positive value when sediment is dropped)
   */
  float updateSedimentQty(float heightDiff, float depositFactor, float erosionFactor) {
    float sedimentAmount = 0; // Positive if sediment is taken from the terrain, negative if sediment is deposited 

    if (sedimentQty > maxSedimentCapacity || heightDiff > 0) { // Depositon
      if (heightDiff > 0) {
        sedimentAmount = -min(sedimentQty, heightDiff);
      } else {
        sedimentAmount = -(sedimentQty - maxSedimentCapacity) * depositFactor;
      }
    } else { // Erosion
      sedimentAmount = min((maxSedimentCapacity - sedimentQty) * erosionFactor, -heightDiff);
    }

    sedimentQty += sedimentAmount;
    return -sedimentAmount;
  }

  /*
   * Update the speed of the drop
   *
   * @Param heightDiff The height difference between its old position and its new position
   * @Param gravity The gravity force (Somehow)
   */
  void updateSpeed(float heightDiff, float gravity) {
    float temp = speed * speed + heightDiff * gravity;
    if (temp < 0) {
      speed = 0;
    } else {
      speed = sqrt(temp);
    }
  }

  /*
   * Update the water volume of the drop
   *
   * @Param heightDiff The height difference between its old position and its new position
   * @Param gravity The gravity force (Somehow)
   */
  void updateWaterVolume(float evaporationFactor) {
    waterVolume *= (1 - evaporationFactor);
  }

  /*
   * @Return the position of the drop as a PVector
   */
  PVector getPos() {
    return pos;
  }

  /*
   * Change the inertia value of the drops
   * Make sure that the value is in range [0,1]
   */
  void setInertia(float value) {
    if (value > 1) {
      inertia = 1;
      return;
    } 

    if (value < 0) {
      inertia = 0;
      return;
    }

    inertia = value;
  }

  /*
   * Change the minimum sediment capacity of a drop
   * Make sure that the value is in range [0,infinity[
   */
  void setMinSedimentCapacity(float value) {     
    if (value < 0) {
      minSedimentCapacity = 0;
      return;
    }

    minSedimentCapacity = value;
  }
}

class DropInfo {
  public float gradX;
  public float gradY;
  public float h;
  public float x;
  public float y;
  public int i; // The x coordinate of the closet top left grid point
  public int j; // The y coordinate of the closet top left grid point
  public float u; // The x offset from the closet top left grid point
  public float v; // The y offset from the closet top left grid point

  DropInfo(float x, float y, int i, int j, float h, float gradX, float gradY) {
    this.x = x;
    this.y = y;
    this.i = i;
    this.j = j;
    this.h = h;
    this.gradX = gradX;
    this.gradY = gradY;
  }
}

DropInfo getDropInfo(final float[][] terrain, final Drop drop) {
  // Get drop position
  float x = drop.getPos().x;
  float y = drop.getPos().y;

  // Get the closest top left vertex
  int i = (int)x;
  int j = (int)y;

  // Get the offset of that drop from the closet top left vertex
  float u = x - (int)x;
  float v = y - (int)y;

  // Get the height of the surrounding vertices
  float heightTL = terrain[i][j];
  float heightTR = terrain[i + 1][j];
  float heightBL = terrain[i][j + 1];
  float heightBR = terrain[i + 1][j + 1];

  // Get the height of the drop by linear interpolation
  float h = (1-v) * ( (1-u) * heightTL + u * heightTR ) + v * ( (1-u) * heightBL + u * heightBR );

  // Compute X gradient by linear interpolation
  // Positive gradients = rigth direction
  float gradX = (heightTL - heightTR) * (1 - v) + (heightBL - heightBR) * v;

  // Compute Y gradient by linear interpolation
  // Positive gradients = down direction
  float gradY = (heightTL - heightBL) * (1 - u) + (heightTR - heightBR) * u;

  return new DropInfo(x, y, i, j, h, gradX, gradY);
}

class DropPath {
  ArrayList<PVector> pos;

  DropPath() {
    pos = new ArrayList<PVector>();
  }

  void add(float x, float y, float z) {
    pos.add(new PVector(x, y, z));
  }

  void show(float scaleXY, float scaleZ) {
    pushMatrix();
    translate(scaleXY * pos.get(0).x, scaleXY * pos.get(0).y, scaleZ * pos.get(0).z);
    noStroke();
    fill(200, 20, 20);
    sphere(0.05);
    popMatrix();

    for (int i = 1; i < pos.size(); i++) {
      pushMatrix();
      translate(scaleXY * pos.get(i).x, scaleXY * pos.get(i).y, scaleZ * pos.get(i).z);
      noStroke();
      fill(200, 20, 20);
      sphere(0.05);
      popMatrix();
      
      stroke(200, 20, 20);
      strokeWeight(0.01);
      line(scaleXY * pos.get(i-1).x, scaleXY * pos.get(i-1).y, scaleZ * pos.get(i-1).z, scaleXY * pos.get(i).x, scaleXY * pos.get(i).y, scaleZ * pos.get(i).z);
    }
  }
}

class PRNG {
  private long initSeed;
  private long seed;
  private long a;
  private long c;
  private long m32;
  private int[] p;
  private PVector gradients[];
  private final float perlinBound = 1. / sqrt(2);

  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************

  /*
   * All code in this section was written by Jeremy Douglass
   * You can see the original thread here: https://discourse.processing.org/t/cross-language-random-numbers/3474
   */

  /*
   * Constructor
   * 
   * Defaults to current system time in milliseconds.
   */
  PRNG() {
    this(System.currentTimeMillis());
  }

  /* 
   * @param seed starting seed value 
   */
  PRNG(long seed) {
    this.a = 1664525;
    this.c = 1013904223;
    this.m32 = 0xFFFFFFFFL;
    this.randomSeed(seed);
    initPerlin();
  }
  
  
  /*
   * Updates seed based on past value, returns new value.
   * 
   * @return current seed
   */
   long getInitSeed() {
     return this.initSeed;
   }

  /*
   * Updates seed based on past value, returns new value.
   * 
   * @return unsigned 32bit value stored in a long
   */
  long nextLong() {
    this.seed = this.seed * a + c & m32;
    return this.seed;
  }

  /**
   * Returns next long wrapped to the max integer value.
   * Implemented so that calls to nextInt or nextLong stay in sync,
   * the next seed internally is always updated to the next long.
   * 
   * @return positive 
   */
  int nextInt() {
    return (int)(this.nextLong()%Integer.MAX_VALUE);
  }

  /**
   * Set the current seed.
   * 
   * @param newSeed the new seed value
   */
  void randomSeed(long newSeed) {
    this.seed = newSeed;
    this.initSeed = newSeed;
  }

  /**
   * @return a random float between 0 and 1
   */
  float random() {
    return random(0, 1);
  }

  /**
   * @param max maximum value to return
   * @return a random float between 0 and max
   */
  float random(float max) {
    return random(0, max);
  }  

  /**
   * @param min minimum value to return
   * @param max maximum value to return
   * @return a random float in the specified range (min, max)
   */
  float random(float min, float max) {
    return map(this.nextInt(), 0, Integer.MAX_VALUE, min, max);
  }

  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************
  // ***********************************************************************************************************************************************  

  /*
   * I built the perlin noise on top of the work above
   */

  /*
   * @return a random int between 0 included and maxInt excluded
   */
  int nextInt(int maxInt) {
    return nextInt() % maxInt;
  }

  /*
   * Initialize the parameters of the perlin noise
   * 
   * 16 unit gradients evenly spaced are created
   * Those gradients will be place on a grid randomly
   *
   * A permutation array p is randomly created
   * This array allows to randomly place the gradients on the grid 
   */
  private void initPerlin() {
    // Create 16 gradients evenly spaced around the unit circle
    gradients = new PVector[16];
    float stepSize = TWO_PI / 16;
    for (int i = 0; i < 16; i++) {
      gradients[i] = PVector.fromAngle(i * stepSize);
    }

    // Create a pseudo-random permutation array that will be used to randomly place the gradients on the grid
    int[] partA = new int[256]; // Needs to be a multiple of to give the same chance to all the gradients;
    int[] partB = new int[256]; // Needs to be a multiple of to give the same chance to all the gradients;

    for (int i = 0; i < 256; i++) {
      partA[i] = i;
      partB[i] = i;
    }

    shuffleArray(partA);
    shuffleArray(partB);

    // p consists of 2 parts to avoid beeing out of bounds values in the perlin code
    p = new int[512];
    for (int i = 0; i < 256; i++) {
      p[i] = partA[i];
      p[i + 256] = partB[i];
      ;
    }
  }

  /*
   * @param x x coordinate of the perlin noise
   * @param y y coordinate of the perlin noise
   * @return a random float between 0 and 1
   */
  float perlin(float x, float y) {
    // Get the top left grid corner coordinates
    int i = (int)x & 255; // = x % 256
    int j = (int)y & 255; // = y % 256

    // Get the gradients used in the 4 surrounding grid corners
    int gradTL = ((p[p[p[i] + j]]) & 15); // = % 16 - Top left gradient (in processing coordonate system)
    int gradTR = ((p[p[p[incr(i)] + j]]) & 15); // = % 16 - Top right gradient (in processing coordonate system)
    int gradBL = ((p[p[p[i] + incr(j)]]) & 15); // = % 16 - bottom left gradient (in processing coordonate system)
    int gradBR = ((p[p[p[incr(i)] + incr(j)]]) & 15); // = % 16 - Top left gradient (in processing coordonate system)

    // get the delta from the top left grid corner
    float u = x - (int)x;
    float v = y - (int)y;

    // Compute the dot products between the gradients at the 4 surrounding grid corners and the 4 vectors starting at the corners and heading towards the point (x, y)
    float dotProdTL = dotProd(gradients[gradTL], u, v);
    float dotProdTR = dotProd(gradients[gradTR], u - 1, v);
    float dotProdBL = dotProd(gradients[gradBL], u, v - 1);
    float dotProdBR = dotProd(gradients[gradBR], u - 1, v - 1);

    // Interpolate the values on the x-axis
    float topAverage = lerpValues(dotProdTL, dotProdTR, u); 
    float bottomAverage = lerpValues(dotProdBL, dotProdBR, u);

    // Interpolate and the y-axis and return the result
    return map(lerpValues(topAverage, bottomAverage, v), -perlinBound, perlinBound, 0, 1);
  }

  /*
   * @param x x coordinate of the perlin noise
   * @param y y coordinate of the perlin noise
   * @param nbOFreq The number of frequencies to add on top of another
   * @param freqFactor The amount by which to change the next frequency
   * @param ampFactor The amount by which to change the next amplitude
   *
   * @return a random float between 0 and 1
   */
  float mixPerlin(float x, float y, int nbOFreq, float freqFactor, float ampFactor) {
    float result = 0;
    float freq = 1;
    float amp = 1;
    final float maxValue = (1 - pow(ampFactor, nbOFreq)) / ( 1 - ampFactor);

    for (int i = 0; i < nbOFreq; i++) {
      result += perlin(x * freq, y * freq) * amp;
      amp *= ampFactor;
      freq *= freqFactor;
    }

    return result/maxValue;
  }

  /*
   * @param x x coordinate of the perlin grid point
   * @param y y coordinate of the perlin grid point
   *
   * @return The gradient used at that perlin grid point
   */
  PVector getGrad(int x, int y) {
    int i = x & 255; // = x % 256
    int j = y & 255; // = y % 256

    int gradTL = ((p[p[p[i]+j]]) & 15);
    return gradients[gradTL];
  }

  /*
   * @param a The start point
   * @param b The end point
   * @param t The value to which interpolate 
   *
   * @return The interpolation at t between a and b. The function used is not linear in order to get continuous second derivative 
   */
  private float lerpValues(float a, float b, float t) {
    return a + fade(t) * (b - a);
  }

  /*
   * @param t The value at which evaluate the function
   *
   * @return The value of the function 6 t^5 - 10 t^4 + 15 t^3 
   */
  private float fade(float t) {
    return t * t * t * ( t * ( t * 6 - 15 ) + 10 );
  }

  /*
   * @param v The first vector
   * @param x The x componant of the second vector
   * @param y The y componant of the second vector
   *
   * @return the dot product of the 2 vectors
   */
  private float dotProd(PVector v, float x, float y) {
    return v.x * x + v.y * y;
  }

  /*
   * Increase a value by 1 and make sure it does not go above 255
   * if so return 0 instead
   */
  private int incr(int a) {
    return a > 254 ? 0 : a + 1;
  }

  /*
   * Function written by PhiLho here: https://stackoverflow.com/questions/1519736/random-shuffling-of-an-array
   * Shuffle the passed array randomly
   */
  private void shuffleArray(int[] ar)
  {
    for (int i = ar.length - 1; i > 0; i--)
    {
      int index = nextInt(i + 1);
      // Simple swap
      int a = ar[index];
      ar[index] = ar[i];
      ar[i] = a;
    }
  }
}

```

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [May 6, 2019, 1:17am UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948/2 "2019-05-06T01:17:36Z")

</div>

Hi @jb4x,

Funny, I experimented different kind of hydraulic erosion algorithms with Processing about a month ago, including the one you’re referring to.

> [@jb4x](#):
>
> The code below definitively gives some results although they are not completely as expected

The results look decent to me, what exactly puzzle you with these outputs ? Sebastian Lague (author of the video you posted) 's [implementation](https://github.com/SebLague/Hydraulic-Erosion) could be useful if you want to to some tweaking. For comparison my version based on his code looked like this:

 ![22](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/a/aa4f5b53af73292cb6282a7c1c18d583683cd0b7.jpeg)

I also ran across this [great article](https://github.com/dandrino/terrain-erosion-3-ways) from Daniel Andrino on github where he explores different ways to mimic hydraulic erosion (Machine Learning, BFS and simulation). I found his numpy-based simulation to be quite similar to Sebastian Lague’s work

 ![58](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/c/c66252c5d1a2d8a84aaa6d5823b598690aa5771f.jpeg)

Lastly, if you have a thing for beautilful maps and color gradient you might want to check this [paper](http://people.csail.mit.edu/bkph/papers/Hill-Shading.pdf) from Berthold Horn for hill shading. Here below some of my notes:

 ![](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/1/1bee52ef1b9285eac17609f16928e2a447ee0ca5.png)  
 ![](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/7/7a637723ffd2631a1a6faf9d5ae232bc4f4c05eb.png)

---

<div class="post-metadata">

**Author:** ![jb4x](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jb4x/32/789_2.png) [@jb4x](https://discourse.processing.org/u/jb4x)\
**Post date:** [May 6, 2019, 6:31am UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948/3 "2019-05-06T06:31:35Z")

</div>

Hi @solub,

Thanks for all the readings, I’ll definitely dig into it when I have the time! That seems to be good stuff 😛

Your version looks gorgeous too, I can’t figure out the right parameters to get those nice creases you get down hill.

I tried several parameters but can’t manage to get them. I also tried copying Sebastian Lague parameters but it doesn’t work.

I wonder if I don’t have a weird behavior in my code but it is quite hard to debug…

---

<div class="post-metadata">

**Author:** ![jeremydouglass](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jeremydouglass/32/20_2.png) [@jeremydouglass](https://discourse.processing.org/u/jeremydouglass)\
**Post date:** [May 7, 2019, 4:53pm UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948/4 "2019-05-07T16:53:41Z")

</div>

Thanks for sharing this, @jb4x – and @solub

@jb4x, is there any possibility that you aren’t seeing certain details because you need a higher resolution total space?

---

<div class="post-metadata">

**Author:** ![jb4x](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jb4x/32/789_2.png) [@jb4x](https://discourse.processing.org/u/jb4x)\
**Post date:** [May 7, 2019, 4:57pm UTC](https://discourse.processing.org/t/generating-terrain-with-hydraulic-erosion/10948/5 "2019-05-07T16:57:04Z")

</div>

Hi @jeremydouglass,

I also thought of that and tried to increase my grid resolution to 1024\*1024 (instead of 512\*512) but I still can’t get those creases.

Also I noticed that my version is more noisy than the one of @solub and Sebastian Lague so I have probably something a bit off in my code.