Diamond Square Algorithm Implementation Problem

Im working on an implementation of the Diamon Square algortihm, and I can't quite seem to get it to look right.

I'm working off this tutorial http://www.lighthouse3d.com/opengl/terrain/index.php?mpd2 which might help explain some of my code.

The function that I know I need to change somehow is below.

It simply takes the min and max X/Z coordinates for the current square being subdivided (which index into an array of Y-Values where these min's and max's already have a Y value computed) and a displacement.

I'm pretty sure the step where I calculate a new midpoint is correct (setting E) but the part that I know is wrong is where I set the midpoint of each side (F G H I)

The forumla seems like it should work, but on high levels of recursion I see almost veritcal drops instead of a smoother transition:


An alternate approach I tried that -almost- worked was just doing:

   yMap[xMin][zMidpnt] = (A + C ) / 2;  // F
   yMap[xMidpnt][zMin] = (A + B) / 2 ;  // G
   yMap[xMax][zMidpnt] = (B + D)/ 2;  // H
   yMap[xMidpnt][zMax] = (C + D) / 2 ;  // I

Which produces a more better landscape but isn't really performing the square step of the algorithm correctly (I think) because it does not add any random displacement to these midpoints.

Any help on how to fix this would be appreciated.

void genLandscapeRec(int xMin, int xMax, int zMin, int zMax, float disp) {
   float nDisp =  (float)(disp * (float)pow(2, -(double)roughnes)); // New displacement
   int xMidpnt = (xMin+xMax) / 2;
   int zMidpnt = (zMin+zMax) / 2;
   GLfloat A = yMap[xMin][zMin];
   GLfloat B = yMap[xMax][zMin];
   GLfloat C = yMap[xMin][zMax];
   GLfloat D = yMap[xMax][zMax];
   GLfloat E = (A + B + C + D) / 4 + randDisp(disp); // Midpoint height
   yMap[xMidpnt][zMidpnt] = E; // Set E
   yMap[xMin][zMidpnt] = (A + C + E) / 3 + randDisp(disp);  // F
   yMap[xMidpnt][zMin] = (A + B + E) / 3 + randDisp(disp);  // G
   yMap[xMax][zMidpnt] = (B + D + E) / 3 + randDisp(disp);  // H
   yMap[xMidpnt][zMax] = (C + D + E) / 3 + randDisp(disp);  // I
   if ((xMidpnt - xMin) > 1) { // Subdivide if new midpoint length will be > 1
      genLandscapeRec(xMin, xMidpnt, zMin, zMidpnt, nDisp);
      genLandscapeRec(xMidpnt, xMax, zMin, zMidpnt, nDisp);
      genLandscapeRec(xMin, xMidpnt, zMidpnt, zMax, nDisp);
      genLandscapeRec(xMidpnt, xMax, zMidpnt, zMax, nDisp);
float randDisp(float disp) {
    return (((float)rand()) / ((float)RAND_MAX)) * (disp + (disp/3)) - disp/3;

Open in new window

Who is Participating?
I wear a lot of hats...

"The solutions and answers provided on Experts Exchange have been extremely helpful to me over the last few years. I wear a lot of hats - Developer, Database Administrator, Help Desk, etc., so I know a lot of things but not a lot about one thing. Experts Exchange gives me answers from people who do know a lot about one thing, in a easy to use platform." -Todd S.

iMr_KiAuthor Commented:
The following globals and initialization function may also be helpful:
// fractal landscape globals
static GLfloat recDepth = 9;
static GLfloat roughnes = 1;
static GLfloat yMap[1025][1025];
void genLandscape(void) {
   int dMax = (int)pow(2, recDepth);
   int dMin = 0;
   float disp = 3;
   yMap[dMin][dMax] = 0.0;
   yMap[dMin][dMax] = 0.0;
   yMap[dMax][dMin] = 0.0;
   yMap[dMax][dMax] = 0.0;
   genLandscapeRec(dMin, dMax, dMin, dMax, disp);

Open in new window

hey iMr_Ki,

probably the randDisp-function is returning wrong (too high) values. check the values it is returning with an assert-ation

iMr_KiAuthor Commented:
randDisp works fine, its at the bottom of the code snippet, simply returns floats in the range of the parameter passed (with a preference on positive values).
Cloud Class® Course: Microsoft Windows 7 Basic

This introductory course to Windows 7 environment will teach you about working with the Windows operating system. You will learn about basic functions including start menu; the desktop; managing files, folders, and libraries.

I sort of got this working better....

The issue was the test for "doneness"....

this test
   if ((xMidpnt - xMin) > 1) { // Subdivide if new midpoint length will be > 1

didn't return values in the expected range

it returned large intervals.... and it was never less than zero.


ie, xMin was always positive.. and so was xMidpiont....

the sides of the square were monotonically increasing... so in the interval 1 to 10 for instnace... xMid is 5, xMin is 1....

but midpoint never crosses 1.. so interval is big...


iMr_KiAuthor Commented:
Jaggedness was being caused by overwriting edge midpoints that had already been computed.

Solution would be to code snippet to something like:

   yMap[xMin][zMidpnt] = (A + C + E) / 3 + randDisp(disp);  // F
   yMap[xMidpnt][zMin] = (A + B + E) / 3 + randDisp(disp);  // G
   yMap[xMax][zMidpnt] = (B + D + E) / 3 + randDisp(disp);  // H
   yMap[xMidpnt][zMax] = (C + D + E) / 3 + randDisp(disp);  // I

Checking to see if midpoints had been previously set.

Experts Exchange Solution brought to you by

Your issues matter to us.

Facing a tech roadblock? Get the help and guidance you need from experienced professionals who care. Ask your question anytime, anywhere, with no hassle.

Start your 7-day free trial
It's more than this solution.Get answers and train to solve all your tech problems - anytime, anywhere.Try it for free Edge Out The Competitionfor your dream job with proven skills and certifications.Get started today Stand Outas the employee with proven skills.Start learning today for free Move Your Career Forwardwith certification training in the latest technologies.Start your trial today
Game Programming

From novice to tech pro — start learning today.

Question has a verified solution.

Are you are experiencing a similar issue? Get a personalized answer when you ask a related question.

Have a better answer? Share it in a comment.