# Diamond-square terrain generation problem

The problem is that I get these steep cliffs all over the map. It happens on the edges, when the terrain is recursively subdivided:

Here is the source:

void DiamondSquare(unsigned x1,unsigned y1,unsigned x2,unsigned y2,float range)
{
int c1 = (int)x2 - (int)x1;
int c2 = (int)y2 - (int)y1;
unsigned hx = (x2 - x1)/2;
unsigned hy = (y2 - y1)/2;
if((c1 <= 1) || (c2 <= 1))
return;

// Diamond stage
float a = m_heightmap[x1][y1];
float b = m_heightmap[x2][y1];
float c = m_heightmap[x1][y2];
float d = m_heightmap[x2][y2];
float e = (a+b+c+d) / 4 + GetRnd() * range;

m_heightmap[x1 + hx][y1 + hy] = e;

// Square stage
float f = (a + c + e + e) / 4 + GetRnd() * range;
m_heightmap[x1][y1+hy] = f;
float g = (a + b + e + e) / 4 + GetRnd() * range;
m_heightmap[x1+hx][y1] = g;
float h = (b + d + e + e) / 4 + GetRnd() * range;
m_heightmap[x2][y1+hy] = h;
float i = (c + d + e + e) / 4 + GetRnd() * range;
m_heightmap[x1+hx][y2] = i;

DiamondSquare(x1, y1, x1+hx, y1+hy, range / 2.0);   // Upper left
DiamondSquare(x1+hx, y1, x2, y1+hy, range / 2.0);   // Upper right
DiamondSquare(x1, y1+hy, x1+hx, y2, range / 2.0);   // Lower left
DiamondSquare(x1+hx, y1+hy, x2, y2, range / 2.0);       // Lower right


}

Parameters: (x1,y1),(x2,y2) - coordinates that define a region on a heightmap (default (0,0)(128,128)). range - basically max. height. (default 32)

Help would be greatly appreciated.

• Without looking hard at your code, it looks like you probably have the wrong corners in the wrong calls in the 4 recursive calls at the end. The map looks like each square is rotated/flipped before calculating the next set, thus sub-dividing the map on strange cliffs. The bottom edge of the top right square looks like it matches the right edge of the top left square, and so on. – DampeS8N Sep 20 '12 at 14:04
• I'm not sure what you mean. The center of the coordinate system is in the top left corner, x axis points to the right and y down. So in the first iteration (x1=0,y1=0),(x2=128,y2=128) and (x1+hx=64,y1+hy=64) is the center of the square. The square is thus divided into 4 subsquares: ((0,0)(64,64)),((64,0)(128,64)),((0,64)(64,128)) and ((64,64)(128,128)). Looks fine to me... – kafka Sep 20 '12 at 14:45

In each subdivision level, the "square" step relies on the results of the "diamond step". But it also factors in the diamond-step produced in the adjacent cell, which you aren't accounting for. I'd rewrite the DiamondSquare function to iterate Breadth-first, instead of depth-first as you currently have it.

Your first issue is that since you re-calculate the square edges twice, it ignores the contribution of the adjacent centerpoint. For example, in the article you reference,

P = (J + G + K + E)/4 + RAND(d)


P = (J + G + J + E)/4 + RAND(d)


i.e. it factors in the current centerpoint twice, not the adjecent centerpoint. This is why you need to go breadth-first, so that you have the previous centerpoints calculated.

Here's my code and the output: .

void DiamondSquare(unsigned x1, unsigned y1, unsigned x2, unsigned y2, float range, unsigned level) {
if (level < 1) return;

// diamonds
for (unsigned i = x1 + level; i < x2; i += level)
for (unsigned j = y1 + level; j < y2; j += level) {
float a = m_heightmap[i - level][j - level];
float b = m_heightmap[i][j - level];
float c = m_heightmap[i - level][j];
float d = m_heightmap[i][j];
float e = m_heightmap[i - level / 2][j - level / 2] = (a + b + c + d) / 4 + GetRnd() * range;
}

// squares
for (unsigned i = x1 + 2 * level; i < x2; i += level)
for (unsigned j = y1 + 2 * level; j < y2; j += level) {
float a = m_heightmap[i - level][j - level];
float b = m_heightmap[i][j - level];
float c = m_heightmap[i - level][j];
float d = m_heightmap[i][j];
float e = m_heightmap[i - level / 2][j - level / 2];

float f = m_heightmap[i - level][j - level / 2] = (a + c + e + m_heightmap[i - 3 * level / 2][j - level / 2]) / 4 + GetRnd() * range;
float g = m_heightmap[i - level / 2][j - level] = (a + b + e + m_heightmap[i - level / 2][j - 3 * level / 2]) / 4 + GetRnd() * range;
}

DiamondSquare(x1, y1, x2, y2, range / 2, level / 2);
}


• Yes, I was also thinking along the lines of breadth-first approach. These fractals are always causing me problems. It was the same with Perlin noise and L-systems. You're awesome. – kafka Sep 21 '12 at 8:17

One possibility is that you are taking a shortcut with your implementation that the algorithm on your linked page does not.

For the square stage, you are calculating the height of the points with

float f = (a + c + e + e) / 4 + GetRnd() * range;
m_heightmap[x1][y1+hy] = f;


which the page's algorithm indicates to use if you are wrapping your map. This gives the appearance that you are using the "next square over"'s height value to calculate this one's. In the simplest, first case, the central point (with height 'e') is used on both the left and right sides to calculate f.

However, the algorithm you reference has you use the actual values of the other squares/diamonds to help you calculate the value of the height of this square point. In their algorithm, the second level point is calculated with the following formula:

N = (K + A + J + F)/4 + RAND(d)


Notice the lack of duplication of a value in there?

I think you might want to try to go for using the non-wrapping versions of the formulae given, those will recurse better, I think.

F = (A + C + E)/3 + ...