Midpoint Displacement Algorithm

Question

This question has come mainly out of sheer desperation, after spending several hours trying to figure out the problem.

If you cast your eyes to the picture above, you should see that my midpoint displacement algorithm algorithm is working (somewhat) successfully; in producing a somewhat coherent noise pattern.

However, it is leaving behind a black dotted grid on the image, and I have no idea why. I can foresee this being an issue in the math, but I just cannot see it; nor was this pointed out in any online resources as a possible issue; so any help will be appreciated as to hunting down this bug.

unsigned char** mdp(unsigned char** base, unsigned base_n, unsigned char r) {
    size_t n = (2 * base_n) - 1;

    unsigned char** map = new unsigned char*[n];
    for (unsigned i = 0; i < n; ++i) map[i] = new unsigned char[n];

    // Resize
    // 1 0 1
    // 0 0 0
    // 1 0 1
    for (size_t i = 0; i < n; i += 2) {
        for (size_t j = !(i % 2 == 0); j < n; j += 2) {
            map[i][j] = base[i / 2][j / 2];
        }
    }

    // Diamond algorithm
    // 0 0 0
    // 0 X 0
    // 0 0 0
    for (size_t i = 1; i < n; i += 2) {
        for (size_t j = 1; j < n; j += 2) {
            unsigned char& map_ij = map[i][j];

            unsigned char a = map[i - 1][j - 1];
            unsigned char b = map[i - 1][j + 1];
            unsigned char c = map[i + 1][j - 1];
            unsigned char d = map[i + 1][j + 1];
            map_ij = (a + b + c + d) / 4;

            unsigned char rv = std::rand() % r;
            if (map_ij + r < 255) map_ij += rv; // EDIT: <-- thanks! the bug! `map_ij + rv`, not `r`
            else map_ij = 255;
        }
    }

    // Square algorithm
    // 0 1 0
    // 1 0 1
    // 0 1 0
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = (i % 2 == 0); j < n; j += 2) {
            unsigned char& map_ij = map[i][j];

            // get surrounding values
            unsigned char a = 0, b = a, c = a, d = a;
            if (i != 0) a = map[i - 1][j];
            if (j != 0) b = map[i][j - 1];
            if (j + 1 != n) c = map[i][j + 1];
            if (i + 1 != n) d = map[i + 1][j];

            // average calculation
            if (i == 0) map_ij = (b + c + d) / 3;
            else if (j == 0) map_ij = (a + c + d) / 3;
            else if (j + 1 == n) map_ij = (a + b + d) / 3;
            else if (i + 1 == n) map_ij = (a + b + c) / 3;
            else map_ij = (a + b + c + d) / 4;

            unsigned char rv = std::rand() % r;
            if (map_ij + r < 255) map_ij += rv;
            else map_ij = 255;
        }

    }

    return map;
}

If you have any tips or resources other than http://www.gameprogrammer.com/fractal.html and http://www.lighthouse3d.com/opengl/terrain/index.php?mpd2 for fractal based terrain generation, I would appreciate them as comments also.

Edit:

This is the new image, as per Fabians suggestion (ty), however it still has some strange quirks, which you should be able to see straight away (small 'dimples' everywhere).

What could be causing this strange behavior? Updated source code: http://www.pastie.org/1924223

Edit:

Many thanks to Fabian in finding the bounds check error, for those interested, here is the current solution as 512x512 png. And current source code (modified by Fabian).

Answer 1

The algorithm recursively adds a value, but the value can be positive or negative (normally +-1 / (2^octave) )

If you start at zero and only add positive values, then you can only go up, and that is why you're seeing the vertices pulled down.

try starting at 127 rather than zero for the four corners, and also try signed char (then checking your bounds both top and bottom)

EDIT

so, two more things need to change in the main (64>>i) to get the half effect at each octave, and also your output function (the one that maps the final[][] tp the imgdta[], you just need to put in

imgdta[(i * n) + j] = 128 + final[i][j];

rather than the if else block.

another thing, I'm not sure why, but your bounds check is failing (that's lines 38 and 65) if you remove the check totally, you notice some new dark blobs too, so i reckon you might need to promote to a larger type before doing the bounds check if you want the more noisy picture you get with "64/i".

ANOTHER EDIT

just found out what it was, you're comparing against 'r', not 'rv', in the bounds check. Here's the fixed code: http://pastie.org/1927076

Answer 2

Two things that jump out:

1. Do you have a compelling reason to do this in fixed-point? There's nothing wrong with it per se, and there are plenty of reasons for using it (most notably memory requirements if you're planning on going up towards a HUGE map), but I would definitely start with a floating-point version of the algorithm and convert it to fixed-point after you've got it working; that should, if nothing else, eliminate one plausible source of errors (in particular, I suspect your clamping may be causing problems, and the conditions for when to add/substract rv).
2. While it's hard to tell from the code I see, it doesn't seem like your randomized height displacement is scaling with the level, and that combined with the issue in (1) may be causing some of the problems - you shouldn't be displacing by the same amount at each level.

Oh, and one non-algorithmic thing: I strongly recommend not doing allocations in your mdp() function; pass in two different already-allocated arrays and do the iteration 'in-place', going from one to the other. If nothing else, this'll let you ping-pong back and forth as you do the layers rather than having to allocate a new array each time out.

Answer 3

Further to the above, you're currently not deleting the memory that you're allocating. To remedy this, change line 104 from:

for (unsigned i = 1; i < 6; ++i) final = mdp(final, n, 64 / i);

to

for (unsigned i = 1; i < 6; ++i) {
  signed char** new_final = mdp(final, n, 64 / i);
  for (unsigned i = 0; i < n; ++i)
    delete[] final[i];
  delete[] final;
  final = new_final;
}

and add this after writing out to the tga file:

for (unsigned i = 0; i < n; ++i)
  delete[] final[i];
delete[] final;