GPU Terrains

Introduction

My first thermal erosion algorithm has a lot of good features, and it’s faster than other ones, but I’m not totally satisfied by the timing performance, so I’ve developed a new thermal erosion algorithm that aims to be as fast as possible and to give good visual results as my first one.

Some technical info

In order to have a shorter execution time I’ve abandoned the reversed structure based on two phases of the previous implementations.

Now all the fragments can compute how much terrain has to go out and how much terrain has to come in without any data computed by the neighborhoods, this allows to reduce the execution times more than 50% (because the first phase disappear), as I’m going to show in the next benchmarks post.

Of course the new model is more empirical than my first one, and requires some assumptions that can give results that are less realistic, from a simulation point of view, but the results are good enough from the visual view, so this algorithm is a good compromise between performance and realism.

Erosion in action

As my first algorithm, also the new one continue the erosion for any number of iterations:

The previous group of images shows the same terrain eroded by 100, 200, 400 and 800 iterations. More iterations means a smoother and lower (on average) terrain.

The next image shows how the new algorithm erodes a (blurred) cube of rock:

The top left panel is the result of 50 iterations, the following ones (in clockwise) are the result of 100, 200 and 400 iterations. As it’s possible to notice, the new algorithm gives results similar to ones obtained by the first and third algorithms, but the eroded terrain is smoother like the one eroded by the second algorithm.

The new algorithm has just one parameter, c, that affects the level of erosion of the algorithm, as showed in the following image:

A lower value of c (0.25 in the top left panel) means more erosion, so a lower terrain, instead a bigger value (1.0 in the bottom right panel) means less erosion and an higher terrain.

Finally I want to compare the result of the erosion of my two algorithms:

On the left (bottom and top panels) there’s the same image of a terrain eroded by 100 iterations of my first algorithm, on the right there’s the terrain eroded by 50 iterations (top panel) and by 100 iterations (bottom panel).

I’ve compared a terrain eroded by 100 iterations with a terrain eroded by the half of iterations because I think they are more similar (on average height) than the two terrains eroded by 100 iterations and this is a good thing because that means I can get results similar to the ones I can get from the first algorithm using half the iterations, so shortening the times any more.

Introduction

I’ve developed a thermal erosion algorithm from scratch, reusing some ideas from the two previous algorithm I’ve implemented and mixing them with some ideas I had working on this project.

A new algorithm is born

I’ve developed a new algorithm because I was not totally satisfied by the visual results of the first two algorithms I’ve implemented, in order to understand what I mean, look at the following image:

This is the result of 200 iterations of the first thermal erosion algorithm on a (blurred) cube of rock, does this look realistic? Is it possible to see a side of a mountain so smooth in nature? I guess no!

The second algorithm gave worse visual results, even if i’s faster than the first one.

So I’ve decided to develop a new algorithm.

Some technical info

In my algorithm I’ve kept the reversed structure used for the other algorithms, so also my algorithm is based on two phases on each step, but my algorithm doesn’t compute the minimum or maximum difference among heights as the other algorithms do, so it’s faster (about 7-10% compared to the first one).

My algorithm uses the talus angle “T” and the “c” constant as they are used in the first algorithm, but I’ve introduced another value “F“, called erosion Factor, that allows me to do some funny things, as I’m going to show you later.

Furthermore it needs just 3 components of a fragment (as the second algorithm), so I guess it’s possible to make it faster using RGB textures instead of RGBA ones that are required by the first algorithm.

New features and differences with other algorithms

The first important feature is that I can generate terrains that are really similar to ones eroded by the first algorithm, look at this image for an example:

The first pyramid has been obtained by 200 iterations of the first erosion algorithm, the second pyramid has been obtained by 200 iterations of my algorithm, it’s really similar to the first one, instead, the last pyramid, has been obtained by 250 iterations of my algorithm, it’s almost the same of the first one.

Another important feature is that the first erosion algorithm can alter the terrain just for a limited number of iterations (about 300), instead my algorithm continue to level the terrain for thousands of iterations as it’s possible to notice in the following image:

The first two squares (from the top) are the results of 200 and 1000 iterations of the first erosion algorithm, while the other two squares are the results of 200 and 1000 iterations of my algorithm.

But the most important feature is that I can obtain different eroded terrains just varying the new parameter F as showed in the following image:

These pyramids are the results of 200 iterations of my erosion algorithm on a cube of rock (not blurred this time), they differ just for the F value (it’s doubled each time) and they look much more realistic, in my opinion!

What does this mean in practice? Here an example:

Changing the T and c values I can obtain just two different kind of eroded terrains using the first algorithm

These two images are the results of 200 iterations of erosion, the first one is based on the values: T = 0.005 - c= 0.5, while the second one on the values: T = 0.003 - C = 0.1, so just two different kind of eroded terrains and they are not so different!

Instead my algorithm can obtain different terrains just varying the F parameter, furthermore it can generate terrains that are quite different changing the other parameters, look at these images:

The first group of images is obtained using the default values for T and c (T = 0.005 and c = 0.05), but varying the F value from 1 to 8 (doubling the value each time), instead the second group of images is generated using the values T = 0.003 - c = 0.1 and varying the F value as in the first group. Visual results are much better than the first algorithm, aren’t they?

For all these reasons I think that my algorithm is better than the other ones I’ve implemented.

Introduction

I’ve developed another thermal erosion algorithm originally proposed by Benes, Marak, Simek and Slavik using a fragment shader.

Also this time I’ve used the reversed approach showed in the first erosion algorithm.

The original algorithm

The algorithm is executed examining each cell of the map and its 8 neighbours according to a Moore scheme like showed in the following picture:

For each neighbour it’s computed the difference between the processed cell and the neighbour:

d[i] = h - h[i];

the minimum positive difference is stored in d_min, and the sum of all the positive differences (this number is n) is stored in d_tot.

Now it’s possible to update all the n cells (where d[i] is bigger than 0.0) using this formula:

h[i] = h[i] + h_out * (d[i] / d_tot);

where h_out is the amount of material redistributed from the central cell (h) to its neighbours, this value is computed by the following formula:

h_out =  n * d_min / (n + 1.0);

so the central cell is updated removing h_out from its height:

h = h - h_out;

My reversed version

My implementation uses two phases to perform the erosion as in the previous algorithm:

in the first phase the shader inspects all the neighbours according to the Moore scheme so to compute d_min, d_tot and n, then, using these values, it computes the amount of material to redistribute h_out. Finally, the values: h, h_out and d_tot are stored in the first three components of the fragment (R, G, B) so to be used in the second phase.

the second phase (another execution of the shader) uses the values stored in all the neighbours to compute how much material has to be moved to the current fragment using the formula:

h_in = h_in + h_out[i] * d[i] / d_tot[i]);

Finally the height is updated summing the two values to the original one:

h = h + h_in - h_out;

A graphical example

The following image shows the erosion performed on a cube of rock, the 4 squares are the results of 50, 100, 150 and 200 steps of erosion:

As it’s possible to notice, this thermal erosion is less effective then the first one, it acts more like a blur.

Introduction

I’ve developed the thermal erosion algorithm originally proposed by Musgrave, Kolb and Mace using a fragment shader.

I’ve altered the algorithm in order to get over the limits of the shaders, in fact a fragment shader can’t scatter data, it can just gather (a shader can read all the fragments, but it can write just the current one).

The original algorithm

The algorithm is executed examining each cell of the map and its 8 neighbours according to a Moore scheme like showed in the following picture:

For each neighbour it’s computed the difference between the processed cell and the neighbour:

d[i] = h - h[i];

the maximum positive difference is stored in d_max, and the sum of all the positive differences that are bigger than T (this numer is n), the talus angle, is stored in d_tot.

Now it’s possible to update all the n cells (where d[i] is bigger than T) using this formula:

h[i] = h[i] + c * (d_max - T)  * (d[i] / d_tot);

and the main cell with this other formula:

h = h - (d_max - (n * d_max * T / d_tot));

The Talus angle T is a threshold that determines which slopes are affected by the erosion, instead the c constant determines how much material is eroded.

My reversed version

As I’ve mentioned in the introduction I can’t implement the original algorithm with a fragment shader because it can’t update the neighbour fragments.

My implementation uses two phases to perform the erosion:

in the first phase the shader inspects all the neighbours according to the Moore scheme so to compute d_max, d_tot and n, then these values are stored in the last three component (G, B, A) of the texel, instead the first one (R) stores the height value.

the second phase (another execution of the shader) uses the values stored in all the neighbours to compute how much material has to be moved to the current fragment using the formula:

h_in = h_in + (c * (d_max[i] - T) * d[i] / d_tot[i]);

then it’s computed how much material should be removed from the current fragment using this formula:

if(n > 0)
	h_out = d_max - (n * d_max / d_tot * T);
else
	h_out = 0.0;

Finally the height is updated summing the two values to the original one:

h = h + h_in - h_out;

A graphical example

The following image shows the erosion performed on a cube of rock, the 4 squares are the result of 50, 100, 150 and 200 steps of erosion:

I’ve got this result setting T=0.005 and c=0.5 .

Useful resources

• Jacob Olsen, Realtime procedural terrain generation

Introduction

Noise synthesis is a fast method to generate realistic terrains with fractal properties. A good way to obtain noise synthesis is to generate coherent noise. Coherent noise is a smooth pseudo-random noise that has 3 main properties:

1. the same input generates the same output
2. a small change in the input gives a small change in the output
3. a big change in the input gives a random change in the output

Perlin noise is a coherent noise generated interpolating a noise function.

Perlin Noise

A way to generate Perlin noise is using the sum of several coherent noise functions with increasing frequencies and decreasing amplitude, each function added to the final noise is called an octave because each function has double the frequency of the previous one.

The following image shows 8 octaves of coherent noise:

frequency and amplitude for each octave are computed by the following code:

freq = pow(2.0, i);
ampl = pow(persistence, i);

Where i is the current octave number (from 0 to 7 in this case) and persistence is a value that determines the amplitude variation, a bigger value (> 0.5) means a rough output, while a smaller value (< 0.5) means a smooth output.

Adding more octaves gives a more detailed noise that can be used as heightmap, this is the result of the sum of 8 octave:

Implementation notes

I’ve implemented the code described by Hugo Elias with a fragment shader, so I’m not explaining the whole code, because it’s really similar to his pseudo-code, the only different function is the random number generator I use, it’s inspired by a tutorial of the Lumina IDE.

This is the formula:

1.0 - 2.0 * fract(sin(dot((cos(vec2(x, y)) * sin(seed)), vec2(15.12345, 91.98765))) * 115309.76543);

This random number generator returns a value in the range [-1.0, 1.0], this is important for something I’m going to explain later., the seed value is passed as uniform variable.

I’ve chosen to use the cosine interpolation as interpolation function instead of the built-in mix function because it gives better visual results summing 4 or less octaves and there’s no much gain in the execution time.

I’ve added another little improvement to the original code, a factor that allows to zoom the generated map, this is a common property in the fractal world.

The following image shows the same map from several zoom levels:

Some enhancements

A fast way to generate new terrains from a starting one is to apply some mathematical function to the base map. In order to get effective visual results the values has to be in the range [-y, +x], so the values has not to be normalized before the filtering.

The following image shows some different functions applied to the terrain showed above in this post:

The top-left square has been generated using the sin function and shows a stronger contrast and more erosion, the top-right one uses the abs function and shows a stronger erosion, the bottom-left square is obtained by the sqrt function and shows interesting lines that looks like rivers, finally, the bottom-right square is generated by the composition of the sin and abs functions and shows the main features of the first two squares combined.

Useful resources

• Hugo Elias: Perlin noise
• Libnoise: glossary
• Ziggyware: Focus on terrain generation, Perlin noise