Diamond Square Algorithm - How to fix horrible diamond shaped artefacts and spikes?

Question

I was attempting to follow this article explaining the implementation of a wrapping version of the diamond square algorithm.

I tried running it with dimensions of 128x128 and a "feature size" of 16, like in the article, but the resultant terrain is this monstrosity:

I read through past questions on this topic and am aware that there can be some issues with somewhat noticeable diamond or square shaped artifacts when using this algorithm, but this is definitely beyond that.

Here is my code:

private float GetHeightAtPoint(int x, int z)
{
    if (x < 0) { x = width + x; }
    if (z < 0) { z = height + z; }

    return fractal[(x % (width - 1)) + (z % (height - 1)) * width];
}

private void SetHeightAtPoint(int x, int z, float value)
{
    if (x < 0) { x = width + x; }
    if (z < 0) { z = height + z; }

    fractal[(x % (width - 1)) + (z % (height - 1)) * width] = value;
}

private void SquareStep(int x, int z, int size, float value)
{
    int halfSize = size / 2;

    float a = GetHeightAtPoint(x - halfSize, z - halfSize);
    float b = GetHeightAtPoint(x + halfSize, z - halfSize);
    float c = GetHeightAtPoint(x - halfSize, z + halfSize);
    float d = GetHeightAtPoint(x + halfSize, z + halfSize);

    float meanHeight = (a + b + c + d) / 4.0f;
    float adjustedHeight = meanHeight + value;

    SetHeightAtPoint(x, z, adjustedHeight);
}

private void DiamondStep(int x, int z, int size, float value)
{
    int halfSize = size / 2;

    float a = GetHeightAtPoint(x - halfSize, z - halfSize);
    float b = GetHeightAtPoint(x + halfSize, z - halfSize);
    float c = GetHeightAtPoint(x - halfSize, z + halfSize);
    float d = GetHeightAtPoint(x + halfSize, z + halfSize);

    float meanHeight = (a + b + c + d) / 4.0f;
    float adjustedHeight = meanHeight + value;

    SetHeightAtPoint(x, z, adjustedHeight);
}

private void DiamondSquare(int stepSize, float scale)
{
    int halfStep = stepSize / 2;
    for (int z = halfStep; z < height + halfStep; z += stepSize)
    {
        for (int x = halfStep; x < width + halfStep; x += stepSize)
        {
            SquareStep(x, z, stepSize, Random.Range(-1.0f, 1.0f) * scale);
        }
    }

    for (int z = 0; z < height; z += stepSize)
    {
        for (int x = 0; x < width; x += stepSize)
        {
            DiamondStep(x + halfStep, z, stepSize, Random.Range(-1.0f, 1.0f) * scale);
            DiamondStep(x, z + halfStep, stepSize, Random.Range(-1.0f, 1.0f) * scale);
        }
    }
}

/*
 * Generate a fractal map using diamond square algorithm and store in fractal[]
 */
private void GenerateFractal(int featureSize)
{
    fractal = new float[width * height];

    for (int z = 0; z < height; z += featureSize)
    {
        for (int x = 0; x < height; x += featureSize)
        {
            SetHeightAtPoint(x, z, Random.Range(0.0f, maxHeight));
        }
    }

    int sampleSize = featureSize;
    float randomScale = 1.0f;

    while (sampleSize > 1)
    {
        DiamondSquare(sampleSize, randomScale);

        sampleSize = sampleSize / 2;
        randomScale = randomScale / 2.0f;
    }
}

Could anyone give me any guidance here? One thing that was a bit suspect was that I had to make slightly different wrapping array access functions in order to make this run as I was getting negative values. You can see that my functions have the if statements checking for negatives, whereas the original implementation does not.

Thanks in advance

Answer 1

The main issue is your code for DiamondStep is the same as for SquareStep. As a result, it is populating the uninitialized target point using other uninitialized points. It should be using points in a diamond around the target:

float a = GetHeightAtPoint(x - halfSize, z);
float b = GetHeightAtPoint(x + halfSize, z);
float c = GetHeightAtPoint(x, z + halfSize);
float d = GetHeightAtPoint(x, z - halfSize);

As a minor issue, the article you linked wraps the coordinates using x & (width - 1). This assumes width is a power of 2, and works for positive or negative coordinates. If you want to do the same with the modulus operator it should be x % width, not x % (width - 1).