UnityTerrainErosionGPU

Hydraulic and thermal erosion implement in Unity using compute shaders.

This is an example implementation of hydraulic and thermal erosion with shallow water equations. My initial motivation was to implement a game mechanic like in the From Dust game.

Disclaimer

This project is still in progress. Hydraulic erosion requires a bit more parameter tweaking and revisiting the actual implementation. But overall it works.

Demo

Just run Main Scene and use the mouse to draw modify terrain/water.

How it works

Explanation is still in progress...

Data / notation

We first need to list all the data we will operate with during the simulation. Since we are using a grid-based simulation we will need some per-cell information. Assuming we have a grid of size $W \times H$ and current cell coordinates are .

Simulation parameters
- - grid width
- - grid height
- $\Delta t$ - simulation time delta
- - gravity
- - cell size along axis
- - cell size along axis
- - pipe length
- - pip cross-section area
- - evaporation rate
- - rain rate
- Hydraulic erosion parameters
  - - sediment capacity
  - - soil suspension rate
  - - soil deposition rate
  - $l_{max}$ - erosion depth limit
- Thermal erosion parameters
  - - thermal erosion rate
  - $K^s_\alpha$ - talus angle scale
  - $K^b_\alpha$ - talus angle bias
Grid-values
- $h^{terrain}_{i,j}$ - height of terrain in cell , must be positive
- $h^{water}_{i,j}$ - height of water in cell , must be positive
- $s_{i,j}$ - suspended sediment amount in cell
- $r_{i,j}$ - terrain hardness in cell
- $f^{left}_{i,j}$ , $f^{right}_{i,j}$ , $f^{up}_{i,j}$ , $f^{down}_{i,j}$ - water flow (flux) in each direction in cell
- $p^{left}_{i,j}$ , $p^{right}_{i,j}$ , $p^{up}_{i,j}$ , $p^{down}_{i,j}$ - terrain mass flow (flux) in each direction in cell
- $v^x_{i,j}$ - water velocity in direction in cell
- $v^y_{i,j}$ - water velocity in direction in cell

Step 1. Water sources

We need to increase the water amount from water sources (e.g. rain). Right now the brush-drawing is omitted for clarity, but it happens in this step.

$h^{water}_{i,j} := h^{water}_{i,j} + \Delta t K_e$

Step 2. Water flow computation

Water flow is proportional to the height difference in each cell because of pressure emerging from the height of water. Thus, to compute flow in each direction of cell we need to calculate the height differences:

$\Delta h^{left}_{i,j} := (h^{terrain}_{i,j} + h^{water}_{i,j}) - (h^{terrain}_{i-1,j} + h^{water}_{i-1,j})$

And repeat for other directions. Then we need to compute the outgoing amount of water from the current cell to the neighbors (in each direction) which is proportional to the volume diefference - height difference $\Delta h^{left}_{i,j}$ multiplied by cell area ( $l_x \times l_y$ ):

$f^{left}_{i,j} := k \cdot max\Big(0,f^{left}_{i,j} + \frac{\Delta t g l_x l_y h^{left}_{i,j}}{l} \Big)$

and also repeat the computation for each direction. - is a special scaling factor to prevent situations when total outflow in 4 directions is higher than water volume in cell (since each direction is computed independently). Thus is defined to split the scale the outflow accordingly - cell water volume divided by total outflow per single step:

$k = min\Big(1, \frac{h^{water}_{i,j} l_x l_y}{\Delta t \big[f^{left}_{i,j} + f^{right}_{i,j} + f^{up}_{i,j} +f^{down}_{i,j} \big]}\Big)$

If outflow will exceed the total volume, the fraction will be less than 1 thus flow will be reduced.

Boundaries

If we need the water to bounce off the walls we need to disable outgouing flow at the boundaries:

$f^{left}_{0,j} := 0$

$f^{right}_{W-1,j} := 0$

$f^{down}_{i,0} := 0$

$f^{up}_{i,H-1} := 0$

At the end of this step, we got all the outgouing flow computed at each cell so we now need to use this flow information to adjust the height of the water.

Step 3. Appling the water flow

The water height in each should increase by the total amount of incoming flow and decrease by the amount of total outgoing flow.

$f^{out} = f^{left}_{i,j} + f^{right}_{i,j} + f^{up}_{i,j} +f^{down}_{i,j}$

outgoing flow is the same as at the previous step.

$f^{in} = f^{right}_{i-1,j} + f^{left}_{i+1,j}+ f^{up}_{i,j-1} + f^{down}_{i,j+1}$

incoming flow is the outgoing flow from neighbor cells in opposite directions.

The total volume change of the column is:

$\Delta V =\Delta t ( f^{in} - f^{out} )$

Finally apply the volume change to the water column, since we store the height and not the volume we need to divide by cell area:

$h^{water}_{i,j} := h^{water}_{i,j} + \frac{\Delta V}{l_x l_y}$

And that's it for water flow. It is also called shallow water equations using pipe model. At the end of this step, the water can slide down the terrain, create vertical waves and so on. But to apply erosion we will need to do more stuff.

Step 3.5. Compute water velocity

In further computation we will need the information about the water velocity in each cell. We can compute it using information about the water flow:

$v^{x}_{i,j} := \frac{f^{right}_{i-1,j} - f^{left}_{i,j} + f^{right}_{i,j} - f^{left}_{i+1,j}}{2}$

$v^{y}_{i,j} := \frac{f^{up}_{i,j-1} - f^{down}_{i,j} + f^{up}_{i,j} - f^{left}_{i,j+1}}{2}$

$\vec{v}_{i,j} = \{v^{x}_{i,j}; v^{y}_{i,j} \}$

the velocity in each axis is the average total flow in each pipe along that axis. So for axis we have 2 neighbor cells: and , thus we can compute the total flow for each neighbor and average across neighbors.

Note: This is only partially physically accurate. For the true velocity - this amount should be scaled by something (include pipe area and length).

Step 4. Hydraulic erosion and deposition

While water flows over terrain it takes (erodes) and transports some amount of soil. After a while, some suspended sediment will be deposited to the ground. This process is mostly defined by the sediment transport capacity of the water flow. There are many complex models regarding these processes, but we will use the simple empirical equation:

$C_{i,j}=K_c \cdot sin(\alpha_{i,j}) \cdot |\vec{v}_{i,j}|$

Step 5. Thermal Erosion

...

Project structure

Shaders
- Shaders/Erosion.compute - all computational stuff happening there in form of separate compute kernels (functions) acting like passes and responsible for different things. Look through that file if you are interested in the actual algorithm implementation.
- Shaders/Water.shader - Surface shader for rendering water plane. In vertex shader vertex positions are updated from state texture and normals are computed. It has basic lighting and alpha decay depending on depth.
- Shaders/Surface.shader - A lit shader to render the terrain surface. In vertex shader vertex positions are updated from state texture and normals are computed.
- Shaders/InitHeightmap.shader - A special shader to initialize initial state from common grayscale heightmap texture. Since state texture is a float texture and operates with values higher than 1 the original heightmap texture should be scaled. This shader is used in the special material used in Simulation.cs.
Scripts/Simulation.cs - main Monobehavior responsible for compute shader setup, dispatching computation to the GPU, texture creation, parameter sharing, and dispatching drawing.
Scripts/ChunkedPlane.cs - main Monobehavior responsible for terrain mesh creation.

References

Mei, Xing, Philippe Decaudin, and Bao-Gang Hu. "Fast hydraulic erosion simulation and visualization on GPU." 15th Pacific Conference on Computer Graphics and Applications (PG'07). IEEE, 2007. http://www.nlpr.ia.ac.cn/2007papers/gjhy/gh116.pdf
Jákó, Balázs, and Balázs Tóth. "Fast Hydraulic and Thermal Erosion on GPU." Eurographics (Short Papers). 2011. http://old.cescg.org/CESCG-2011/papers/TUBudapest-Jako-Balazs.pdf
Warning: poor paper quality along with math mistakes. But has nice ideas.
Št'ava, Ondřej, et al. "Interactive terrain modeling using hydraulic erosion." Proceedings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation. Eurographics Association, 2008. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.173.5239&rep=rep1&type=pdf

TODO:

Better explanation of the implementation
More descriptive comments in code
Quality of life:
- Better editor - camera controls and better brush controls
- Different initial state loaders (from terrain data, from 16bit textures, from .raw)
- Terrain chunks to simplify rendering of distant terrain parts

bshishov / unityterrainerosiongpu Goto Github PK