SebastionErosion.cpp

    // ported C# to Godot C++: "Coding Adventure: Hydraulic Erosion" by Sebastian Lague
    class SebastionErosion {
       public:
        struct HeightAndGradient {
           public:
            float height;
            float gradientX;
            float gradientY;
        };

        int seed = 0;
        // [Range (2, 8)]
        int erosionRadius = 3;
        // [Range (0, 1)]
        float inertia = .05f;              // At zero, water will instantly change direction to flow downhill. At 1, water will never change direction.
        float sedimentCapacityFactor = 4;  // Multiplier for how much sediment a droplet can carry
        float minSedimentCapacity = .01f;  // Used to prevent carry capacity getting too close to zero on flatter terrain
                                           // [Range (0, 1)]
        float erodeSpeed = .3f;
        // [Range (0, 1)]
        float depositSpeed = .3f;
        // [Range (0, 1)]

        float evaporateSpeed = .01f;
        float gravity = 4;
        int maxDropletLifetime = 30;

        float initialWaterVolume = 1;
        float initialSpeed = 1;

        // Indices and weights of erosion brush precomputed for every node
        std::vector<std::vector<int>> erosionBrushIndices;

        std::vector<std::vector<float>> erosionBrushWeights;

        Ref<RandomNumberGenerator> prng;

        int currentSeed = 0;
        int currentErosionRadius = 0;
        int currentMapSize = 0;

        // Initialization creates a System.Random object and precomputes indices and weights of erosion brush
        void Initialize(int mapSize, bool resetSeed) {
            // if (resetSeed || prng == null || currentSeed != seed) {
            //     prng = new System.Random(seed);
            //     currentSeed = seed;
            // }

            prng.instantiate();
            prng->set_seed(currentSeed);

            // if (erosionBrushIndices == null || currentErosionRadius != erosionRadius || currentMapSize != mapSize) {
            InitializeBrushIndices(mapSize, erosionRadius);

            currentErosionRadius = erosionRadius;
            currentMapSize = mapSize;
            // }
        }

        void Erode(PackedFloat32Array& map, int mapSize, int numIterations = 1, bool resetSeed = false) {
            Initialize(mapSize, resetSeed);

            for (int iteration = 0; iteration < numIterations; iteration++) {
                // Create water droplet at random point on map

                float posX = prng->randf_range(0, mapSize - 1);
                float posY = prng->randf_range(0, mapSize - 1);

                float dirX = 0;
                float dirY = 0;
                float speed = initialSpeed;
                float water = initialWaterVolume;
                float sediment = 0;

                for (int lifetime = 0; lifetime < maxDropletLifetime; lifetime++) {
                    int nodeX = posX;
                    int nodeY = posY;
                    int dropletIndex = nodeY * mapSize + nodeX;
                    // Calculate droplet's offset inside the cell (0,0) = at NW node, (1,1) = at SE node
                    float cellOffsetX = posX - nodeX;
                    float cellOffsetY = posY - nodeY;

                    // Calculate droplet's height and direction of flow with bilinear interpolation of surrounding heights
                    HeightAndGradient heightAndGradient = CalculateHeightAndGradient(map, mapSize, posX, posY);

                    // Update the droplet's direction and position (move position 1 unit regardless of speed)
                    dirX = (dirX * inertia - heightAndGradient.gradientX * (1 - inertia));
                    dirY = (dirY * inertia - heightAndGradient.gradientY * (1 - inertia));
                    // Normalize direction

                    float len = sqrt(dirX * dirX + dirY * dirY);
                    if (len != 0) {
                        dirX /= len;
                        dirY /= len;
                    }
                    posX += dirX;
                    posY += dirY;

                    // Stop simulating droplet if it's not moving or has flowed over edge of map
                    if ((dirX == 0 && dirY == 0) || posX < 0 || posX >= mapSize - 1 || posY < 0 || posY >= mapSize - 1) {
                        break;
                    }

                    // Find the droplet's new height and calculate the deltaHeight
                    float newHeight = CalculateHeightAndGradient(map, mapSize, posX, posY).height;
                    float deltaHeight = newHeight - heightAndGradient.height;

                    // Calculate the droplet's sediment capacity (higher when moving fast down a slope and contains lots of water)
                    float sedimentCapacity = MAX(-deltaHeight * speed * water * sedimentCapacityFactor, minSedimentCapacity);

                    // If carrying more sediment than capacity, or if flowing uphill:
                    if (sediment > sedimentCapacity || deltaHeight > 0) {
                        // If moving uphill (deltaHeight > 0) try fill up to the current height, otherwise deposit a fraction of the excess sediment
                        float amountToDeposit = (deltaHeight > 0) ? MIN(deltaHeight, sediment) : (sediment - sedimentCapacity) * depositSpeed;
                        sediment -= amountToDeposit;

                        // Add the sediment to the four nodes of the current cell using bilinear interpolation
                        // Deposition is not distributed over a radius (like erosion) so that it can fill small pits
                        map[dropletIndex] += amountToDeposit * (1 - cellOffsetX) * (1 - cellOffsetY);
                        map[dropletIndex + 1] += amountToDeposit * cellOffsetX * (1 - cellOffsetY);
                        map[dropletIndex + mapSize] += amountToDeposit * (1 - cellOffsetX) * cellOffsetY;
                        map[dropletIndex + mapSize + 1] += amountToDeposit * cellOffsetX * cellOffsetY;

                    } else {
                        // Erode a fraction of the droplet's current carry capacity.
                        // Clamp the erosion to the change in height so that it doesn't dig a hole in the terrain behind the droplet
                        float amountToErode = MIN((sedimentCapacity - sediment) * erodeSpeed, -deltaHeight);

                        // Use erosion brush to erode from all nodes inside the droplet's erosion radius
                        for (int brushPointIndex = 0; brushPointIndex < erosionBrushIndices[dropletIndex].size(); brushPointIndex++) {
                            int nodeIndex = erosionBrushIndices[dropletIndex][brushPointIndex];
                            float weighedErodeAmount = amountToErode * erosionBrushWeights[dropletIndex][brushPointIndex];
                            float deltaSediment = (map[nodeIndex] < weighedErodeAmount) ? map[nodeIndex] : weighedErodeAmount;
                            map[nodeIndex] -= deltaSediment;
                            sediment += deltaSediment;
                        }
                    }

                    // Update droplet's speed and water content
                    speed = sqrt(speed * speed + deltaHeight * gravity);
                    water *= (1 - evaporateSpeed);
                }
            }
        }

        HeightAndGradient CalculateHeightAndGradient(PackedFloat32Array& nodes, int mapSize, float posX, float posY) {
            int coordX = posX;
            int coordY = posY;

            // Calculate droplet's offset inside the cell (0,0) = at NW node, (1,1) = at SE node
            float x = posX - coordX;
            float y = posY - coordY;

            // Calculate heights of the four nodes of the droplet's cell
            int nodeIndexNW = coordY * mapSize + coordX;
            float heightNW = nodes[nodeIndexNW];
            float heightNE = nodes[nodeIndexNW + 1];
            float heightSW = nodes[nodeIndexNW + mapSize];
            float heightSE = nodes[nodeIndexNW + mapSize + 1];

            // Calculate droplet's direction of flow with bilinear interpolation of height difference along the edges
            float gradientX = (heightNE - heightNW) * (1 - y) + (heightSE - heightSW) * y;
            float gradientY = (heightSW - heightNW) * (1 - x) + (heightSE - heightNE) * x;

            // Calculate height with bilinear interpolation of the heights of the nodes of the cell
            float height = heightNW * (1 - x) * (1 - y) + heightNE * x * (1 - y) + heightSW * (1 - x) * y + heightSE * x * y;

            // return new HeightAndGradient(){height = height, gradientX = gradientX, gradientY = gradientY};
            return HeightAndGradient{height, gradientX, gradientY};
        }

        void InitializeBrushIndices(int mapSize, int radius) {
            // erosionBrushIndices = new int[mapSize * mapSize][];
            // erosionBrushWeights = new float[mapSize * mapSize][];

            erosionBrushIndices.clear();
            erosionBrushIndices.resize(mapSize * mapSize);
            erosionBrushWeights.clear();
            erosionBrushWeights.resize(mapSize * mapSize);

            std::vector<int> xOffsets(radius * radius * 4);
            std::vector<int> yOffsets(radius * radius * 4);
            std::vector<float> weights(radius * radius * 4);

            float weightSum = 0;
            int addIndex = 0;

            for (int i = 0; i < erosionBrushIndices.size(); i++) {
                int centreX = i % mapSize;
                int centreY = i / mapSize;

                if (centreY <= radius || centreY >= mapSize - radius || centreX <= radius + 1 || centreX >= mapSize - radius) {
                    weightSum = 0;
                    addIndex = 0;
                    for (int y = -radius; y <= radius; y++) {
                        for (int x = -radius; x <= radius; x++) {
                            float sqrDst = x * x + y * y;
                            if (sqrDst < radius * radius) {
                                int coordX = centreX + x;
                                int coordY = centreY + y;

                                if (coordX >= 0 && coordX < mapSize && coordY >= 0 && coordY < mapSize) {
                                    float weight = 1 - sqrt(sqrDst) / radius;
                                    weightSum += weight;
                                    weights[addIndex] = weight;
                                    xOffsets[addIndex] = x;
                                    yOffsets[addIndex] = y;
                                    addIndex++;
                                }
                            }
                        }
                    }
                }

                int numEntries = addIndex;
                // erosionBrushIndices[i] = new int[numEntries];
                erosionBrushIndices[i].resize(numEntries);  // warning this might not clear the vars (not sure if a problem)

                // erosionBrushWeights[i] = new float[numEntries];
                erosionBrushWeights[i].resize(numEntries);  // warning this might not clear the vars (not sure if a problem)

                for (int j = 0; j < numEntries; j++) {
                    erosionBrushIndices[i][j] = (yOffsets[j] + centreY) * mapSize + xOffsets[j] + centreX;
                    erosionBrushWeights[i][j] = weights[j] / weightSum;
                }
            }
        }
    };