Problem with diamond-square algorithm

[kronikel]

I decided to start working on my terrain generator again, and so here I am back with another problem.
It's easy to see what is going wrong, but I can't figure out how it is doing this.
There will be a point here and there that gets set either way too high or way too low which creates huge cone shapes.
There are always huge cones sticking out of the top -

And it can be seen very clearly when I set the smoothness up -

No matter how smooth it is there are still these cones sticking into and out of the terrain.

Here is how I generate the heights -

Code:

SideLength = (2^x) + 1
while SideLength >= 2
{
    HalfSide = SideLength / 2;

    //diamond step
    for (int x = 0; x < Size - 1; x = x + SideLength)
        {
             for (int y = 0; y < Size - 1; y = y + SideLength)
             {
                   Heights[x + HalfSide, y + HalfSide] = (Heights[x, y] + Heights[x + SideLength, y] + 
                       Heights[x, y + SideLength] + Heights[x + SideLength, y + SideLength]) / 4 +- Randomness;
              }
        }
    
    //square step
    for (int x = 0; x < Size - 1; x = x + SideLength)
    {
         for (int y = 0; y < Size - 1; y = y + SideLength)
         {
              if (y != 0)
              Heights[x + HalfSide, y] = (Heights[x, y] + Heights[x + SideLength, y] + 
                  Heights[x + HalfSide, y + HalfSide] + Heights[x + HalfSide, y - HalfSide]) / 4 +- Randomness;

              if (x != 0)
              Heights[x, y + HalfSide].y = (Heights[x, y].y + Heights[x, y + SideLength].y + 
                  Heights[x + HalfSide, y + HalfSide].y + Heights[x - HalfSide, y + HalfSide].y) / 4 +- Randomness;
         }
    }
    SideLength = SideLength / 2;
    Randomness = Randomness / Roughness; 
    }
}

[David Makin]

Assuming I'm correct that "randomness" there is just a simple variable then it's really a "scale" and you need to change the "+-randomness" in the diamond and square functions to "+randomness*random()" where random() is a random generator returning -1 to +1.
The "correct" value of "roughness" is sqrt(2) if I remember correctly - better implemented as 1/sqrt(2) and using *roughness instead of divide.

Also one trick to adjust the roughness to give a more realistic effect for landscapes is to compute the average base height value first then use an extra multiplier based on this i.e. instead of just adding randomness*random() compute the base height first and then add randomness*random()*(value based on height difference from zero, useful range to return from say 0.5 for height zero to 1.0 for all height differences >a given magnitude).

[kronikel]

What I posted was just a rough summary of my code, but to generate the random part I use -
Average - r + ((float)(random.NextDouble() * r) * 2)

And if it helps any here is the actual code

Code:

private void CreateHeights()
        {
            if (cbUseLand.Checked == false) 
                return;
            int
                Size = Convert.ToInt32(System.Math.Pow(2, int.Parse(tbDetail.Text)) + 1),
                SideLength = Size - 1,
                d = 1025 / (Size - 1),
                HalfSide;
            Heights = new Point3D[Size, Size];
            float
                r = float.Parse(tbHeight.Text),
                Roughness = float.Parse(RoughnessBox.Text);

            for (int x = 0; x < Size; x++)
                for (int y = 0; y < Size; y++)
                    Heights[x, y] = Make3DPoint(x * d, 740, y * d);

            while (SideLength >= 2)
            {
                HalfSide = SideLength / 2;

                for (int x = 0; x < Size - 1; x = x + SideLength)
                {
                    for (int y = 0; y < Size - 1; y = y + SideLength)
                    {
                        Heights[x + HalfSide, y + HalfSide].y =
                          (Heights[x, y].y +
                          Heights[x + SideLength, y].y +
                          Heights[x, y + SideLength].y +
                          Heights[x + SideLength, y + SideLength].y) / 4 - r + ((float)(random.NextDouble() * r) * 2);
                    }
                }

                for (int x = 0; x < Size - 1; x = x + SideLength)
                {
                    for (int y = 0; y < Size - 1; y = y + SideLength)
                    {
                        if (y != 0)
                            Heights[x + HalfSide, y].y = (Heights[x, y].y + Heights[x + SideLength, y].y + Heights[x + HalfSide, y + HalfSide].y + Heights[x + HalfSide, y - HalfSide].y) / 4 - r + ((float)(random.NextDouble() * r) * 2); 
                        if (x != 0)
                            Heights[x, y + HalfSide].y = (Heights[x, y].y + Heights[x, y + SideLength].y + Heights[x + HalfSide, y + HalfSide].y + Heights[x - HalfSide, y + HalfSide].y) / 4 - r + ((float)(random.NextDouble() * r) * 2);
                    }
                }
                SideLength = SideLength / 2;
                r = r / Roughness;
            }
        }

[David Makin]

In that case the code looks correct, so all I can think is that there's something inherently wrong with "random.NextDouble() " either in terms of it not actually returning randoms from zero to one or it having sufficient non-randomness to cause the peaks - e.g. maybe it's returning values from -1 to +1 rather than 0 to 1 ?

On the subject of diamond-square landscapes I'm working on a (shader-friendly) implementation that will allow distance-estimated ray stepping because it can return the height very quickly for any point - I'm still improving it so that it's randomness extends infinitely rather than being restricted to one square, obviously to get the optimum minimum distance estimation is a little tricky especially as I've got to relearn any of the maths I once knew with respect to what's essentially statistics in this case

[kronikel]

That sounds pretty interesting.
But I'm pretty certain that the randomness is working.
I have tried using Random.Next(-r, r)
and also Random.Next() / (Int32.MaxValue - 1)
Which is another way of generating a number from 0 to 1, but all give me identical results.

I have also made a version of this in delphi where I have lots of experience in generating random numbers, yet still I get identical results.

There has to be something wrong with which points I'm using to average or which points I'm setting or something.
One thing to note is that the points that are messing up aren't in random places, they form a grid.

Edit:
This is very interesting...
When I take out the randomness completely, so all it does is take the average of the 4 surrounding points, I still get these strange "dimples" popping up everywhere.

[David Makin]

Maybe try viewing the heights as colours in a 2D plane just to check that it's not your 3D rendering that's causing the problem ?

[fractower]

The output seems to match what one would expect from the code. It appears to be designed to write to each point only once with an exponentially decaying randomness. The early points are essentially random. When the random factor decays, the algorithm defaults to an approximate solution to a Laplacian with a few point sources (early points). This solution produces a 1/r type behavior relative to the point sources which is essentially what you are seeing.

[kronikel]

That's a good idea David, that's actually how I first wrote the algorithm. So here's what I get -

And another example -


And to further rule out it being the way I render the points I completely made a new method for rendering the points. It is a couple lines long and as basic as it can get. Just drawing wire from point to point -


@fractower - The randomness can't be the issue because when I take out the randomness completely, so it simply sets a height based on the average of the 4 corner heights, I get identical results.

[knighty]

What value of "roughtness" are you using? Maybe it's too high? I'm getting the same effect as you when using too high values. A value around 2 (better when less than 2) should give good results. When roughtness is bigger than two the spikes begins to appear.

Fractower's explanation is right: The algorithm acts somehow as a (multigrid) solver for laplace equation.

Reminds me another terrain generation algorithm: take wihte noise and filter it (using FFT). How is it called?...

[knighty]

It also looks like there is something wrong with your random value. Its mean sould be 0. Otherwise the spikes appear.

[David Makin]

Reminds me another terrain generation algorithm: take wihte noise and filter it (using FFT). How is it called?...

FFT ? - Spectral Synthesis ? - too slow even for pre-computation at runtime unless you have fairly small grids i.e. really requires that the resulting heights are permanently stored for use.
Perlin noise is fast enough for runtime pre-computation and provides pretty decent results - in fact for pre-computed grids Perlin noise is considerably better than diamond-square unless you considerably enhance the d-s algorithm but if you need realtime computation of height (y) given a location (x,z) then on either GPU or CPU a good d-s implementation should be considerably faster - even if the d-s is somewhat enhanced though maybe not to the point of being as "good" aesthetically as Perlin noise.

[kronikel]

I will say it again, the randomness is not the issue at all.
Roughness is not the issue either.
I have used 2 along with many other numbers but no matter how smooth I make the terrain there will still be large dips.

I made a whole completely new algorithm exactly copying this one - http://www.intelegance.net/content/codebase/DiamondSquare.pde
And I still get the dips regardless of roughness/randomness.
I worked out the algorithm by hand and the numbers I got were the same as what the code gives me.
You can easily see by hand how the cones come up, but I can't find a way to keep it from happening.

[lycium]

I think Knighty nailed it: you seem to be adding strictly positive random contributions, instead of having a mean of 0.

[kronikel]

Ok this is getting a little frustrating.
IT HAS NOTHING TO DO WITH RANDOMNESS, and my code returns a number between -r and +r, with a mean of 0.

If you use

Code:

    r = AnyNumberAtAll;
    for (int a = 0; a < 20000; a++)
            {
                f += -r + (float)(random.NextDouble() * (r * 2));
            }
            Write((f / 20000).ToString());

It will average together 20,000 random numbers, and the average comes out to be almost exactly 0 every time, with a max of r and a minimum of -r.

[fractower]

I agree that the problem is not the random number generator. The issue is that the first few points calculated dominate the terrain preventing a regression to the mean. They act and look like tent poles. On the positive side it produces interesting terrain, but on the bad side it produces the spike artifacts.

One solution is to go back an recalculate the initial spikes once their influence has been established. For example run the same alogrithm in parallel but delayed a couple of iterations like a madrigal. Each location ends up being touched twice, but I think it will eliminate the tent poles.