Smoothing issue with Diamond-Square algorithm

Question

I am using the diamond-square algorithm to generate random terrain. It works fine except I get these large cone shapes either sticking out of or into the terrain. The problem seems to be that every now and then a point gets set either way too high or way too low.

Here is a picture of the problem

And it can be better seen when I set the smoothness really high

And here is my 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);

        //seeding all the points
        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;
        }
    }

Answer 1

Gavin S. P. Miller gave a SIGGRAPH '86 talk about how Fournier, Fussel & Carpenter's original algorithm was fundamentally flawed. So what you're seeing is normal for any naive implementation of the Diamond Square algorithm. You will require a separate approach for smoothing, either post each Diamond-Square compound step, or as a post-process to all diamond-square iterations (or both). Miller addressed this. Weighting and box or gaussian filtering are one option; seeding the initial array to a greater degree than just the initial 4 points (i.e., replicating the resultsets of the first few steps of diamond-square either manually or using some built-in intelligence, but instead supplying unbiased values); the more initial information you give the array before increasing the detail using diamond-square, the better your results will be.

The reason appears to be in how the Square step is performed. In the Diamond step, we've taken the average of the four corners of a square to produce that square's centre. Then, in the subsequent Square step, we take the average of four orthogonally-adjacent neighbours, one of which is the square's centre point we just produced. Can you see the problem? Those original corner height values are contributing too much to the subsequent diamond-square iteration, because they are contributing both through their own influence AND through the midpoint that they created. This causes the spires (extrusive and intrusive), because locally-derived points tend more strongly toward those early points... and because (typically 3) other points do as well, this creates "circular" influences around those points, as you iterate to higher depths using Diamond-Square. So these kinds of "aliasing" issues only appear when the initial state of the array is underspecified; in fact, the artifacting that occurs can be seen as a direct geometric consequence of using only 4 points to start with.

You can do one of the following:

• Do local filtering -- generally expensive.
• Pre-seed the initial array more thoroughly -- requires some intelligence.
• Never smooth too many steps down from a given set of initial points -- which applies even if you do seed the initial array, it's all just a matter of relative depths in conjunction with your own maximum displacement parameters.

Answer 2

The actual flaw in the above algorithm is an error in conceptualization and implementation. Diamond square as an algorithm has some artifacting but this is range based artifacts. So the technical max for some pixels is higher than some other pixels. Some pixels are directly given values by the randomness while others acquire their values by the diamond and squared midpoint interpolation processes.

The error here is that you started from zero. And repeatedly added the value to the current value. This causes the range of diamond squared to start at zero and extend upwards. It must actually start at zero and go both up and down depending on the randomness. So the top range thing won't matter. But, if you don't realize this and naively implement everything as added to the value, rather than starting at zero and fluctuating from there, you will expose the hidden artifacts.

Miller's notes were right, but the flaw is generally hidden within the noise. This implementation is shows those problems. That is NOT normal. And can be fixed a few different ways. This was one of the reasons why after I extended this algorithm to remove all the memory restrictions and size restrictions and made it infinite and deterministic1, I then still switched away from the core idea here (the problems extending it to 3d and optimizing for GPUs also played a role.2

Answer 3

Instead of just smoothening with an average, you can use a 2-D median filter to take out extremes. It is simple to implement, and usually generates the desired effect with a lot of noise.

Answer 4

I believe the size of the displacement r in each iteration should be proportional to the size of the current rectangle. The logic behind this is that a fractal surface is scale invariant, so the variation in height in any rectangle should be proportional to the size of that rectangle.

In your code, the variation in height is proportional to r, so you should keep it proportional to the size of your current grid size. In other words: multiply r by the roughness before the loop and divide r by 2 in each iteration.

So, instead of

r = r / Roughness;

you should write

r = r / 2;