# 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

• In the pictures, it looks like each of the dots is at the same height. Are the corners at the same height as well? – deft_code May 18 '11 at 16:06
• For what it's worth, your images look very pretty. :) – ChrisE May 19 '11 at 20:39
• scrand(): I'm not entirely sure that I understand--is this supposed to return a signed character on the interval (-r/2,r/2]? The dimples, to me anyways, seem kind of like the result of overflow. The adjacent areas seem to suddenly shoot black, and then climb back up to white. Overall, it also looks like there is sharp banding, again suggesting to me that you are overflowing, perhaps. Would you mind doing the math in a higher precision (say, integer), and then clamping the values to the range [0,256] or [-128,127]? – ChrisE May 19 '11 at 20:55
• The problem was solved below, it was because I was bounds checking against the range of the random value, not its actual value. The scrand() was a temporary 'noise' function ideally returning [-128, 127] – deceleratedcaviar May 19 '11 at 20:59
• Ah, cool! Glad to hear it's working now. – ChrisE May 19 '11 at 21:03

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

• This was definetly the better way to go about it, but no cigar as of yet, it appears my image still has some 'quirks', I updated the original post. – deceleratedcaviar May 18 '11 at 22:56
• not sure but line 93 looks wrong, 64/i might need to be 64>>i (as you half the effect each octave) – Richard Fabian May 19 '11 at 10:47
• Ahhh!! Thank you so much, I can't believe that, I knew it'd be something stupid for that second problem. Loved your makeshift TGA code, sorry, I should have saved you the trouble and put up the header. – deceleratedcaviar May 19 '11 at 20:49

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.

• Some good points, obviously I am just trying to get the algorithm correct, the implementation is far from ideal, I don't even clean up memory at this point :P. – deceleratedcaviar May 19 '11 at 3:25
• At this point, the scaling over passes is 64 / i, obviously I will change that later, but it doesn't really explain the dimples effect currently experienced. :S – deceleratedcaviar May 19 '11 at 3:33

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;

• I'm very aware on how to clean up memory, but seeing as this was an algorithm prototype only, and would be re-written to use vectors, I just wanted to make sure it was correct. – deceleratedcaviar May 19 '11 at 20:51