c++ - Midpoint Displacement Algorithm

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).

Edit (years later): Python version https://gist.github.com/dcousens/5573724#file-mdp-py

Answer

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