Diamond-square algorithm help

[kronikel]

I am trying to create landscapes/terrain using the diamond-square algorithm.
I have a working program that just uses the midpoint displacement but of course it produces objects that are too "squareish", so I want to implement the diamond step.

I have the concept down of the diamond step, for each square you take the midpoint of each side and randomly shift its height up or down.
My problem is only setting each point ONE TIME.
I end up with floating squares and squares that don't match up because each point should be set once, but I can't find a way to know if its been set or not.

[David Makin]

It happens that I just used the diamond-square algorithm myself for a 1024*1024 array of 3 floats but easily adaptable to say a single integer array or whatever, so here you are (for 3 channels @1024*1024):

Code:

    float vals[1024][1024][3]; // probably needs to be a global or at least malloc'ed
    float r,g,b,rx,gx,bx,rn,gn,bn;
    rx = gx = bx = -2.0f;
    rn = gn = bn = 2.0f;
    vals[0][0][0] = vals[0][0][1] = vals[0][0][2] = 0.0f;
    UInt32 x, y, dx = 512;
    UInt32 sx = 1024;
    float o = 0.353553390593274f, f = 0.707106781186548f/65535.0f;
    do
    {
        y = dx;
        do
        {
            x = dx;
            do
            {
                r = vals[x][y][0] = 0.25f*(vals[x-dx][y-dx][0]+vals[(x+dx)&1023][y-dx][0]+vals[x-dx][(y+dx)&1023][0]+vals[(x+dx)&1023][(y+dx)&1023][0]) + (mtr->MTrand()&65535)*f - o;
                g = vals[x][y][1] = 0.25f*(vals[x-dx][y-dx][1]+vals[(x+dx)&1023][y-dx][1]+vals[x-dx][(y+dx)&1023][1]+vals[(x+dx)&1023][(y+dx)&1023][1]) + (mtr->MTrand()&65535)*f - o;
                b = vals[x][y][2] = 0.25f*(vals[x-dx][y-dx][2]+vals[(x+dx)&1023][y-dx][2]+vals[x-dx][(y+dx)&1023][2]+vals[(x+dx)&1023][(y+dx)&1023][2]) + (mtr->MTrand()&65535)*f - o;
                if (r>rx)
                {
                    rx = r;
                }
                if (r<rn)
                {
                    rn = r;
                }
                if (g>gx)
                {
                    gx = g;
                }
                if (g<gn)
                {
                    gn = g;
                }
                if (b>bx)
                {
                    bx = b;
                }
                if (b<bn)
                {
                    bn = b;
                }
            } while ((x+=sx)<1024);
        } while ((y+=sx)<1024);
        y = 0;
        do
        {
            x = dx;
            do
            {
                r = vals[x][y][0] = 0.25f*(vals[x][(y-dx)&1023][0]+vals[x-dx][y][0]+vals[(x+dx)&1023][y][0]+vals[x][(y+dx)&1023][0]) + (mtr->MTrand()&65535)*f - o;
                g = vals[x][y][1] = 0.25f*(vals[x][(y-dx)&1023][1]+vals[x-dx][y][1]+vals[(x+dx)&1023][y][1]+vals[x][(y+dx)&1023][1]) + (mtr->MTrand()&65535)*f - o;
                b = vals[x][y][2] = 0.25f*(vals[x][(y-dx)&1023][2]+vals[x-dx][y][2]+vals[(x+dx)&1023][y][2]+vals[x][(y+dx)&1023][2]) + (mtr->MTrand()&65535)*f - o;
                if (r>rx)
                {
                    rx = r;
                }
                if (r<rn)
                {
                    rn = r;
                }
                if (g>gx)
                {
                    gx = g;
                }
                if (g<gn)
                {
                    gn = g;
                }
                if (b>bx)
                {
                    bx = b;
                }
                if (b<bn)
                {
                    bn = b;
                }
            } while ((x+=sx)<1024);
            if ((y += dx)<1024)
            {
                x = 0;
                do
                {
                    r = vals[x][y][0] = 0.25f*(vals[x][y-dx][0]+vals[(x-dx)&1023][y][0]+vals[(x+dx)&1023][y][0]+vals[x][(y+dx)&1023][0]) + (mtr->MTrand()&65535)*f - o;
                    g = vals[x][y][1] = 0.25f*(vals[x][y-dx][1]+vals[(x-dx)&1023][y][1]+vals[(x+dx)&1023][y][1]+vals[x][(y+dx)&1023][1]) + (mtr->MTrand()&65535)*f - o;
                    b = vals[x][y][2] = 0.25f*(vals[x][y-dx][2]+vals[(x-dx)&1023][y][2]+vals[(x+dx)&1023][y][2]+vals[x][(y+dx)&1023][2]) + (mtr->MTrand()&65535)*f - o;
                    if (r>rx)
                    {
                        rx = r;
                    }
                    if (r<rn)
                    {
                        rn = r;
                    }
                    if (g>gx)
                    {
                        gx = g;
                    }
                    if (g<gn)
                    {
                        gn = g;
                    }
                    if (b>bx)
                    {
                        bx = b;
                    }
                    if (b<bn)
                    {
                        bn = b;
                    }
                } while ((x+=sx)<1024);
            }
        } while ((y+=dx)<1024);
        dx >>= 1;
        sx >>= 1;
        f *= 0.707106781186548f;
        o *= 0.707106781186548f;
    } while (dx);
; Now we normalise the ranges created:
    r = 255.0f/(rx-rn);
    g = 255.0f/(gx-gn);
    b = 255.0f/(bx-bn);
    y = 0;
    UInt8 *dest=(UInt8*)texMem;
    do
    {
        x = 0;
        do
        {
            *dest++ = (vals[x][y][0]-rn)*r;
            *dest++ = (vals[x][y][1]-gn)*g;
            *dest++ = (vals[x][y][2]-bn)*b;
            *dest++ = 255;
        } while (++x<1024);
    } while (++y<1024);

Ignoring the /65535.0 in the calculation of "f":
"o" should be half the value of f.
The larger the value of f (and o) then the sharper the "landscape" - in fact anything >1/sqrt(2) for f is effectively introducing true randomness beyond fBm - using 1.0 is really true randomness, obviously smaller values of f (and o) will give smoother landscape/plasma.
Note that MTRand is simply an implementation of the Mersenne Twister returning 32-bit randoms.
If you make "o" larger than 0.5*f then landscape at lower "heights" will be rougher and if you make "o" smaller than 0.5*f then landscape at higher "heights" will be rougher.

[kronikel]

I really couldn't understand any part of that o.O
Partly because the way we are doing things is so different and partly because I'm not good with C.

I created a working version where every time a new vertex is made it checks in an array of previously made vertexes to see if that one has already had a height set to it. If it has, it uses the height that was set before, and if not it generates a random height.
But this is the most inefficient thing imaginable, its far too slow.
I'm pretty sure I'm doing this all wrong even though I'm getting good results.
I think you're supposed to work through the squares from left to right, top to bottom, which I'm not doing.

[David Makin]

OK, well I thought after I posted that in fact there are 2 different ways of doing it - either the straightforward method (as I posted in C) where you're creating the complete array or the more tricky method where you need the algorithm done on an individual pixel basis (e.g. as a GLSL routine).
If you don't follow C, here is a pseido-code version that I hope you can follow - this is "off-the-top-of-my-head" so I may make some silly mistakes:

For an area 1024*1024:

float sharpness = 1.0/sqrt(2.0)

define array float v(1024,1024)
float f = 1.0/65535.0
float o = 0.5
int dx = 512
int sx = 1024
int x
int y
repeat
   y = dx
   repeat
      x = dx
      repeat
          v(x,y) = 0.25*(v(x-dx,y-dx) + v((x+dx)&1023, y-dx) + v(x-dx,(y+dx)&1023) + v((x+dx)&1023,(y+dx)&1023)) + (rand()&65535)*f - o
      until ((x=x+sx)>1023)
   until ((y=y+sy)>1023)
   y = 0
   repeat
     x = dx
     repeat
        v(x,y) = 0.25*(v(x,(y-dx)&1023) + v(x-dx,y) + v((x+dx)&1023,y) + v(x,(y+dx)&1023) + (rand()&65535)*f - o
     until ((x=x+sx)>1023)
     y = y + dx
     x = 0
     repeat
        v(x,y) = 0.25*(v(x,y-dx) + v((x-dx)&1023,y) + v((x+dx)&1023,y) + v(x,(y+dx)&1023) + (rand()&65535)*f - o
     untiil ((x=x+sx)>1023)
   until ((y=y+dx)>1023)
   f = f*sharpness
   sx = sx>>1
   o = o*sharpness
until ((dx = (dx>>1))==0)

Now you have an array v(1024,1024) of floats depending on the use you can adapt the above to get the max and min values of the array and then normalise the values based on those as you wish.
Not that using larger values of sharpness increases the randomness/jaggedness of the landscape/plasma and reducing the sharpness smooths the landscape/plasma. Also larger values of sharpness will increase the overall range of values and smaller values will reduce the range.
Using a sharpness of 1/sqrt(2) is the "correct" value.
Changing the offset value "o" so it is more than 0.5 will increase randomness/jaggedness of smaller "heights" and smooth the larger "heights" and changing "o" so it is less than 0.5 will increase the randomness/jaggedness of larger heights and smooth those lower.

[kronikel]

Ok I am trying to understand what is going on in that last code but a there are a couple things I don't know

what is the "&" operator doing? I saw some places that say "&1023" or "&65535"
also the way the loops end, I'm not sure what conditions are being asked for on the "until" part.
"until ((x=x+sx)>1023)"
is that "until (x + sx) > 1023"?

and I'm not sure what's going on at all with "until ((dx = (dx>>1))==0)"

[David Makin]

Ok I am trying to understand what is going on in that last code but a there are a couple things I don't know

what is the "&" operator doing? I saw some places that say "&1023" or "&65535"
also the way the loops end, I'm not sure what conditions are being asked for on the "until" part.
"until ((x=x+sx)>1023)"
is that "until (x + sx) > 1023"?

and I'm not sure what's going on at all with "until ((dx = (dx>>1))==0)"

&1023 is bitwise AND e.g. 1024&1023 is zero 1040&1023 is 16 etc.

"until ((x=x+sx)>1023)"
This is just shorthand for:

   x = x + sx
until (x>1023)

"until ((dx = (dx>>1))==0)"

is

   dx = dx>>1
until (dx==0)

(dx>>1 is dx shifted right by one i.e. essentially dx/2)

[kronikel]

I understand it all now except for the "&"
Is there a way to replace the bitwise operations with something else?
I never have understood them.

Edit:
After looking around at some more examples I've noticed some things.
The dimensions of the starting square should be (x^2) + 1 and also that the points are cycled through using a "delta" or a "step".
My method is completely different and my starting square can be any width and height, and any shape, such as a rectangle.
I don't cycle through points using a step amount, I cycle through squares and calculate the points on that square using midpoints of its parent square.
So It's hard to relate my code to other people's.
I probably shouldn't be doing it this way but I can't seem to figure out how to change it to the normal method of doing it.

[David Makin]

You don't really need to understand the bitwise &, other than that if you're &ing with a value of 2^n-1 such as 15 or 255 or 1023 or 65535 etc. then a&b is exactly the same as "mod(a,b+1)" or "a%(b+1)" or "a mod (b+1)" but several lifetimes faster
The algorithm I've used is adaptable to any rectangle 2^n by 2^m units (with care) and does not produce the nasty squarish artifice associated with the quick-and-dirty centre-sides and centre method
However if you need a per-pixel method then to be honest I'm still trying to find a method more optimum in terms of overrall speed/artifice than the old Vista-Pro triangles method as the diamond-square is decidedly tricky for that purpose especially if you need the square/rectangle to wrap correctly.

[kronikel]

I absolutely can't find an example that works when I convert it to the language I'm working with.
If I could do that, then I could probably get an understanding of the algorithm and be able to make my own.
I converted that last example you gave me into delphi and it compiles but returns an array full of 0's.

[xenodreambuie]

In Delphi, just replace & with " and ". Replace >> with shr, and << with shl, but only for positive integers.

Looking up the parent squares can be efficient if you are storing the randomly generated points after calculating them. Each point needs a flag so you know if it's been calculated already. If you are not going to do that, then you need to instead use a random function that always gives the same value for a given coordinate.

The last time I did any midpoint subdivision it was using the triangle method on an Amiga in the 80's. I don't remember if I tried any tricks to hide ridging, but I used quadratic interpolation of each midpoint for a smooth option, and raymarching to get pixel detail instead of just rendering triangles back to front.

[kronikel]

I just came up with a successful and quick method for only setting my points once.
I create a blank bitmap the size of the original square colored black.
Every time I set the height of point (x, y) I check the color of the pixel (x, y) in the bitmap.
If the color is 0 then I set a random height and color the pixel in the bitmap to match the random height.
If the color is not 0 then I set the point's height to the color of the bitmap at that point.