Scape: 3. Procedural extensions

January 7, 2012 Graphics

This is the third article in the series on Scape, a GPU-based terrain editor, picking up where I left off in the previous article on procedural noise techniques. In this article, I’ll discuss two novel extensions to the more common procedural terrain generation algorithms covered in my previous post.

Panoramic canyon terrain in Scape

If you don’t know what Scape is, please see my first article on Scape, explaining the project and its rendering pipeline, together with video showing what it’s capable of in real time.  Furthermore, as this third article builds on concepts and code explained in the second post on Scape, I would recommend reading that before continuing to read this, if you haven’t done so already. To understand how the procedural algorithms presented here fit into Scape’s brush-based editing pipeline, see the fourth article. The fifth article provides a summary of the project, together with all links to the project’s research, articles, source code and binaries.

The algorithms discussed here have been implemented as Cg pixel shaders for use in Scape’s internal GPU-based brush pipeline. But for the purpose of this article, I’ve extracted those and built an FX Composer project around them for you to play with. The project files can downloaded from the Downloads section below.

Swiss turbulence

Please compare the two images below. (Click on them to enlarge). The one on the left is generated using the standard ridged turbulence function. It’s not bad, but notice how uniform and contextless the terrain looks. It’s full of features and details, but it doesn’t resemble how mountains typically look after being battered by wind, rain and changes in temperature for a long time.

Normally, lots of material on the surface of mountains gets displaced from the peaks to the valleys over time, carving gulleys on its slopes and smoothing the valleys below. And that’s what you miss in the image on the left. There are ways to actually simulate erosion on a heightfield, but they tend to be quite slow as they require many iterations. But you can try to fake it quite successfully, as demonstrated by the image on the right.

The heightfield in the image on right was created in Scape in only a few minutes using a brush based on a novel turbulence algorithm. This algorithm is still based on combining different Perlin noise octaves, but it tries to apply those in a very specific way to create results resembling glacier-eroded mountain ranges. And as I couldn’t come up with a ‘scientific’ name for this specific look, I chose to name it after a country that has similar features in its landscape.

float swissTurbulence(float2 p, float seed, int octaves,
                      float lacunarity = 2.0, float gain = 0.5,
                      float warp = 0.15>)
{
     float sum = 0;
     float freq = 1.0, amp = 1.0;
     float2 dsum = float2(0,0);
     for(int i=0; i < octaves; i++)
     {
         float3 n = perlinNoiseDeriv((p + warp * dsum)*freq, seed + i);
         sum += amp * (1 - abs(n.x));
         dsum += amp * n.yz * -n.x;
         freq *= lacunarity;
         amp *= gain * saturate(sum);
    }
    return sum;
}

Please notice that without the code in red, this algorithm is equivalent to the ridged turbulence function. What the additional pieces of code exactly do is described in the next couple of paragraphs. The perlinNoiseDeriv function is similar to the perlinNoisePseudoDeriv function described in the previous article, but this function actually returns the correct analytic Perlin noise gradient (i.e., the derivatives dn/dx and dn/dy) in n.yz, instead of, well, something else. The value in n.x is again equal to what the (only) direct output of the perlinNoise function would be.

The noise octaves are scaled and summed into sum while iterating over them from coarse to fine, just like in the original turbulence function. But because I wanted to base this algorithm on the ridged turbulence look, I’m summing 1-abs(n.x) instead of n.x. This effectively replaces a call to the equivalent ridgedNoise call. So, there’s nothing really new thus far.

dsum is a sum weighted similarly to sum. But it sums scaled noise gradients instead of the noise values. The 2D dsum is used to do a lookup for finer octaves at a location closer to the nearest ridge, based on the local gradient information from coarser levels. This causes the features on the slopes to get elongated. The scale of this distortion is controlled by warp. Note that applying the chain rule reveals that the gradient for 1-abs(n.x) is actually n.xy*-sign(n.x) and not n.yz*-n.x. But to prevent artifacts caused by the discontinuity in the sign function in the former, I used the latter.

To smoothen the valleys, the algorithm also uses an amp calculation that is slightly different from the original turbulence function. Here, the amp is not only multiplied by gain for each octave, but also by the (clamped version of the) intermediate values of the sum of the noise octaves. This causes the amplitude for finer details to fade out quickly near valleys, but remain quite potent near peaks. To compensate for the on-average decrease in amplitude, the gain parameter should be a bit larger than ‘normal’ (for example, 0.6 instead of 0.5).

As an additional step, the p input can be slightly distorted before passing it on to the swissTurbulence function, offsetting p by a random 2D vector to add more variety to the terrain, distorting (too) smoothly shaped slope and ridge features. The random 2D vector itself can be the result of another two calls to the some turbulence-type function: once for each component in p. Unlike the actual final turbulence function, the turbulence functions used for the distortion should only require a few (e.g. 4) octaves or so, as more wouldn’t be (as) beneficial to this effect. Furthermore, the distortion can also have its own set of scale, gain and lacunarity parameters to make the whole effect more controllable.

What follows is the Cg implementation of the perlinNoiseDeriv function. The code that is identical to that of the perlinNoise function is shown in black. The additional code required to analytically calculate the local gradient is shown in red. In the turbulence function above, it’s used to find the direction towards the nearest ridge. But obviously, it could be used for other purposes as well, like analytically calculating the procedural surface’s normal directions, for example.

float3 perlinNoiseDeriv(float2 p, float seed)
{
    // Calculate 2D integer coordinates i and fraction p.
    float2 i = floor(p);
    float2 f = p - i;

    // Get weights from the coordinate fraction
    float2 w = f * f * f * (f * (f * 6 - 15) + 10); // 6f^5 - 15f^4 + 10f^3
    float4 w4 = float4(1, w.x, w.y, w.x * w.y);

    // Get the derivative dw/df
    float2 dw = f * f * (f * (f * 30 - 60) + 30); // 30f^4 - 60f^3 + 30f^2

    // Get the derivative d(w*f)/df
    float2 dwp = f * f * f * (f * (f * 36 - 75) + 40); // 36f^5 - 75f^4 + 40f^3

    // Get the four randomly permutated indices from the noise lattice nearest to
    // p and offset these numbers with the seed number.
    float4 perm = tex2D(samplerPerlinPerm2D, i / 256) + seed;

    // Permutate the four offseted indices again and get the 2D gradient for each
    // of the four permutated coordinates-seed pairs.
    float4 g1 = tex2D(samplerPerlinGrad2D, perm.xy) * 2 - 1;
    float4 g2 = tex2D(samplerPerlinGrad2D, perm.zw) * 2 - 1;

    // Evaluate the four lattice gradients at p
    float a = dot(g1.xy, f);
    float b = dot(g2.xy, f + float2(-1,  0));
    float c = dot(g1.zw, f + float2( 0, -1));
    float d = dot(g2.zw, f + float2(-1, -1));

    // Bi-linearly blend between the gradients, using w4 as blend factors.
    float4 grads = float4(a, b - a, c - a, a - b - c + d);
    float n = dot(grads, w4);

    // Calculate the derivatives dn/dx and dn/dy
    float dx = (g1.x + (g1.z-g1.x)*w.y) + ((g2.y-g1.y)*f.y - g2.x +
               ((g1.y-g2.y-g1.w+g2.w)*f.y + g2.x + g1.w - g2.z - g2.w)*w.y)*
               dw.x + ((g2.x-g1.x) + (g1.x-g2.x-g1.z+g2.z)*w.y)*dwp.x;
    float dy = (g1.y + (g2.y-g1.y)*w.x) + ((g1.z-g1.x)*f.x - g1.w + ((g1.x-
               g2.x-g1.z+g2.z)*f.x + g2.x + g1.w - g2.z - g2.w)*w.x)*dw.y +
               ((g1.w-g1.y) + (g1.y-g2.y-g1.w+g2.w)*w.x)*dwp.y;

    // Return the noise value, roughly normalized in the range [-1, 1]
    // Also return the pseudo dn/dx and dn/dy, scaled by the same factor
    return float3(n, dx, dy) * 1.5;
}

Jordan turbulence

The heightfield shown in the next two images is, again, fully procedurally generated by mixing Perlin noise functions, without the use of any iterative simulation techniques. I thought it kind of looked like some of the quite eroded mountains in Jordan, hence the name.

Just like the swissTurbulence function, this algorithm uses the perlinNoiseDeriv function to warp and scale noise from finer octaves. Yet, it generates a completely different type of terrain, and exposes many more parameters.

float jordanTurbulence(float2 p, float seed, int octaves, float lacunarity = 2.0,
                       float gain1 = 0.8, float gain = 0.5,
                       float warp0 = 0.4, float warp = 0.35,
                       float damp0 = 1.0, float damp = 0.8,
                       float damp_scale = 1.0)
{
    float3 n = perlinNoiseDeriv(p, seed);
    float3 n2 = n * n.x;
    float sum = n2.x;
    float2 dsum_warp = warp0*n2.yz;
    float2 dsum_damp = damp0*n2.yz;

    float amp = gain1;
    float freq = lacunarity;
    float damped_amp = amp * gain;

    for(int i=1; i < octaves; i++)
    {
        n = perlinNoiseDeriv(p * freq + dsum_warp.xy, seed + i / 256.0);
        n2 = n * n.x;
        sum += damped_amp * n2.x;
        dsum_warp += warp * n2.yz;
        dsum_damp += damp * n2.yz;
        freq *= lacunarity;
        amp *= gain;
        damped_amp = amp * (1-damp_scale/(1+dot(dsum_damp,dsum_damp)));
    }
    return sum;
}

What is quite different here from other turbulence function implementations, is that the coarsest octave is calculated before the actual for loop. This allows the effects from the largest features on the look of smaller features to be scaled independently from the effects of smaller features on even smaller features. Of course, an even more flexible way to implement this would be to provide controllable constants for all individual octaves as uniform arrays (or constant buffers), saving GPU cycles and offering more control at the same time.

As always, sum accumulates the weighed sum of noise octave values. But this time it sums over n2.x, the square of n.x. Squaring the noise produces a slightly billowy effect, causing the jordanTurbulence to a have an even more rounded look in comparison to the sharply-edged swissTurbulence.  Note that when warp0, warp, damp0 and damp are all zero, dsum_warp and dsum_damp remain zero as well. In that case, the algorithm is basically equivalent to a standard turbulence function that sums the squares of perlinNoise octaves. The image on the right shows what that looks like.

According to the chain rule, the gradient of (n.x)2 is 2.0*n.x*n.yz, or 2.0*n2.yz. Hence, dsum_damp is accumulating a scaled version of the gradients. This dsum_damp is then used to smoothly scale down the next octave’s amplitude when the coarser features are relatively flat. This creates an effect that somewhat resembles the effects of thermal erosion.

Unlike the swissTurbulence function, the jordanTurbulence function dampens the amplitude of finer octave based on gradients, not on heights. That means that more horizontal areas near peaks will be treated the same as the more horizontal areas near valleys, allowing relatively flat areas without finer features to form at any height. damp0 controls the damping effect of the coarsest octave on finer octaves. damp controls the damping from finer octaves on even finer octaves. How much the amplitude can be damped by the gradients is controlled by damp_scale.

The dsum_warp variable is calculated quite similarly to dsum_damp, but serves a different purpose. It is used to create kind of gulley-like downhill features, mimicking one of the effects of fluvial erosion. Therefore, it deserved its own tweakable set of parameters. As a side-effect, it also seems to ‘squash’ features of all sizes. warp0 controls the amount of warping the coarsest octave can cause in all other octaves. Similarly, warp controls the warping effect of mid-range features on even smaller features.

All the above steps apply similar concepts as the swissTurbulence function but do it in a slightly different way, resulting in a completely different type of landscape. The same holds for pre-warping p to bend, stretch and compress horizontally using two additional calls to turbulence, as suggested as an optional additional step in swissTurbulence function. When applied correctly to jordanTurbulence‘s input, it helps to hide the ‘average distance’ between peaks, but also stretch features to the point they start to look like strata.

FX Composer

The accompanying downloadable NVIDIA FX Composer project contains both the functions from the previous article and this article. The functions’ outputs are visualized as procedural textures and are generated in real time. The implemented DX9/DX10 FX and (OpenGL) CgFX effects expose all parameters as user-controllable sliders.

jordanTurbulence pixel shader in FX Composer

I hope these algorithms will inspire others to start experimenting with mixing noise functions in novel ways themselves as well, as there’s probably still a lot to be gained in this field. In the next (to be written) article, I will cover the details of Scape’s user-controllable brush pipeline that internally uses the procedural algorithms to do the actual editing on a heightfield.

Downloads

  • FX Composer 2.5 project, extending the project from the previous article with the two new algorithms as FX and CgFX effects.  (ZIP)

Comments (18)

Michael
January 13, 2012

Hello Giliam,

I realy need to thank you for putting this stuff online. Your work is very inspiring and it motivated me to start looking into shaders. I always dreamed of having my own planet-engine abd after doing it on the CPU with limited success, I will now use a fresh approach based on your work.

Wish you all the best,

Michael

Arianna
January 17, 2012

Remarkable information! I have been looking for something such as this for some time now. Many thanks!

Nick Anderson
March 31, 2012

Hi Gilliam,

Thank you very much for making this code available, and for presenting it in such a pleasant readable format!

Warm Regards,
Nick

Nicolaas
August 30, 2012

So I implemented your swiss turbulence in my procedural planet renderer, check out http://irio.co.uk/sites/cirio/?p=788

If you have any ideas on how to implement stuff like, say, more realistic rivers / valleys, let me know =)

Giliam
August 30, 2012

Hi Nicolaas. That’s really awesome! And thanks for sharing! Yes, getting realistically flowing rivers can be quite tricky. But perhaps the people at http://en.spaceengine.org/forum/21-589-1 (who also added the procedural erosion algorithms to their planet generator) might have some ideas. Good luck!

Nicolaas
September 5, 2012

hehe ok thanks. Rivers & valleys can wait, I added some new screenshots with an improved erosion thingy & new cloud systems : http://irio.co.uk/sites/cirio/?p=800 . I think I should start focusing on getting out a preview code, perhaps minecraft-alike game with automatic updates.. so thinking about real-time modification of the landscape, but it’s pretty hard when you don’t have time and have to do other stuff :P

Anyway, thanks for the inspiration, really what I needed just now. Not too many people working with this stuff =)

Michael
January 4, 2013

Hi Im not sure what saturate is supposed to do in the amp *= gain * saturate(sum); For swissTurbulence. Could you please provide some info.

Thanks!

Giliam
January 5, 2013

Sure. As you probably know, saturate(x) is simply a shorthand for clamp(x, 0, 1) (or, x<1?(x>0?x:0):1 ), and it’s purpose is to prevent amp from potentially going through the roof. As n.x is between -1 and 1, the first sum can already be 2.0. The saturate makes sure that amp is always decreasing (as long as gain < 1). There are probably better ways of doing something similar, but this one is simple, robust and efficient. Good luck!

Michael
January 6, 2013

Giliam-

Thanks for the info. For those of us that are playing in oldschool cpu land ie renderman and the like, would you be willing to post a version of youre perlinNoiseDeriv, Swiss and Jordan code that is just straight up C++? Id also love to see a 3d or 4d version.

All the best.

-Michael

Giliam
January 6, 2013

Hi Michael. Sadly, I currently don’t have any C++ implementations to post. But I hope the following tips will make it easier for you to try it out yourself. The perlinNoiseDeriv function is actually an ordinary 3D noise function which only supports discrete steps in its z coordinate (it’s second parameter, that is, which is used in the turbulence functions as a seed, thus effectively sampling different XY planes for different octaves). Moreover, it also calculates and returns the noise’s dN/dx and dN/dy derivates at the same time. But you can do the same with any noise generator using forward differencing. That is, n.x = sampler(x,y,z), n.y = (sampler(x+e,y,z) – n.x)/e, n.z = (sampler(x,y+e,z) – n.x)/e, where sampler(x,y,z) is some 3d noise function and e is a very small positive number. The two turbulence function should be very easy to rewrite to C++. (Obviously, one difference is that vector classes are typically named vector2, vector3 and vector4 in C++, while Cg denotes these as float2, float3 and float4, respectively). Once you’ve got that far, it should be very straighforward to extend into 4D, although I have no idea if it would still produce similar results. Regards, Giliam

David
May 2, 2014

“I hope these algorithms will inspire others to start experimenting with mixing noise functions in novel ways themselves as well, as there’s probably still a lot to be gained in this field.”

You are absolutely correct. The future of textures is procedural, as all texture surfaces can be mathematically described. This terrain is a-mazing, and I’m wondering if it could be incorporated into the Unity3D engine. Standard pillowy puffs for terrain is dead, and this has to be the future. Please let me know if you can.

Jim
October 25, 2014

Thanks for these excellent algorithms, but how do you calculate the derivatives for these recursive functions for the surface normals for lighting calculations? Particularly with the Jordan turbulence function I’m struggling to analytically calculate a derivative.

Giliam
October 28, 2014

Hi Jim. Never tried it, but it should be possible to calculate the x and y derivatives for all helper variables while iterating through the octaves and build it up somewhat similar to the sum itself. But note that as the ‘sum’ is already a function of both the noise and it’s x and y derivative, the x and y derivative of ‘sum’ (i.e. what you’re looking for) will depend on the second derivative of the noise as well. That is, it will be a function of n, dn/dx, dn/dy, ddn/ddx, ddn/ddy, ddn/dxdy. So, although possible to calculate it analytically, a numerical approach (e.g. central differencing using four samples of the Jordan function as-is) is a lot easier to implement and likely to be more efficient for most purposes. Hope this helps!

John
March 13, 2015

Hey Giliam, this is some great work! Very inspiring. Do you have any tips on how to get your swiss turbulence working with 3D noise and 3D derivatives specifically for use around a sphere? It seems like I need to calculate the local up vector for each vertex and apply the step 2 stretch with a Vector3 dsum based on that. I was able to replicate your other steps and its looking pretty cool already, but step 2 seems like the most important part to achieve that eroded look. Thanks for your help!

Giliam
March 14, 2015

Hi John. Thanks! Assuming you already have a way to calculate the 3D derivative of the local 3D noise position p (either analytically or by using central differences), I think you could replace “dsum += amp * n.yz * -n.x;” by “dsum += amp * n.yzw * -n.x;” where n is assumed to now be float4(n, dx, dy, dz). The 3d dsum vector could then be used to apply the warping of p for the next octave. Ideally, you would want to do all of this using a pair of polar coordinates centered around p, but it’s probably also just good enough to do a direct (and much simpler) reprojection onto the sphere in 3d. In that case, you can replace “(p + warp * dsum)*freq” by “normalize(projectOnTangentSurface(p, warp * dsum))*freq” where projectOnTangentSurface(float3 p, float3 delta) equals (1-dot(p, delta))*p + delta. The call to normalize() may perhaps even be overkill, so you can try it without as well. The above assumes that the sphere is centered at the origin and p is on the unit sphere. The final 3d coordinate would then be “sphereCenterPos + p * sphereRadius * (1 + someTweakedConstant*sum)”. I never actually tried this, so no guarantees whatsoever :-). And in any case, this stuff is almost as much an art as it is a science, so be sure to try out and compare some variations on the theme. Hope this helps and be sure to let me know how it works out!

ManuelLog
March 24, 2015

Superb Website, Continue the beneficial work. Thanks for your time!

Kevin Roast
March 1, 2017

Thanks for making this work available – very interesting and well explained. Helped me write me this: http://www.kevs3d.co.uk/dev/shaders/terrain3.html

Robin
March 22, 2019

Amazing! and just what I’ve been looking for!
Thank you!

Leave a comment