Logo by haltenny - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Check out the originating "3d Mandelbulb" thread here
 
*
Welcome, Guest. Please login or register. March 28, 2024, 01:09:58 PM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: Diamond-square algorithm help  (Read 1752 times)
0 Members and 1 Guest are viewing this topic.
kronikel
Navigator
*****
Posts: 79



« on: April 08, 2011, 06:14:27 AM »

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.
Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #1 on: April 08, 2011, 09:01:56 PM »

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

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kronikel
Navigator
*****
Posts: 79



« Reply #2 on: April 08, 2011, 10:55:01 PM »

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.
Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #3 on: April 09, 2011, 06:03:42 AM »

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

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kronikel
Navigator
*****
Posts: 79



« Reply #4 on: April 11, 2011, 03:58:32 AM »

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)"
Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #5 on: April 11, 2011, 04:45:34 AM »

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)
« Last Edit: April 11, 2011, 04:47:16 AM by David Makin » Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kronikel
Navigator
*****
Posts: 79



« Reply #6 on: April 11, 2011, 05:35:58 AM »

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.
« Last Edit: April 11, 2011, 09:46:03 PM by kronikel » Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #7 on: April 11, 2011, 10:33:05 PM »

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 wink
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 smiley
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.
Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kronikel
Navigator
*****
Posts: 79



« Reply #8 on: April 11, 2011, 10:48:46 PM »

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.
Logged
xenodreambuie
Conqueror
*******
Posts: 124



WWW
« Reply #9 on: April 11, 2011, 11:28:32 PM »

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

Regards, Garth
http://xenodream.com
kronikel
Navigator
*****
Posts: 79



« Reply #10 on: April 12, 2011, 02:11:06 AM »

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.
Logged
Pages: [1]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Problem with diamond-square algorithm Programming « 1 2 3 4 » kronikel 56 22744 Last post July 16, 2012, 10:50:15 AM
by setec
Problem with my Diamond Square algorithm Programming mbob61 5 3436 Last post December 12, 2012, 12:17:35 AM
by David Makin
A Loom for Weaving Diamond Fabric Images Showcase (Rate My Fractal) JoeFRAQ 0 707 Last post November 30, 2013, 07:55:34 PM
by JoeFRAQ
Square in square Images Showcase (Rate My Fractal) Caleidoscope 2 980 Last post February 11, 2016, 01:24:45 PM
by Caleidoscope
diamond-square midpoint displacement algorithm Landscape/Terrain Generation claude 0 3278 Last post May 25, 2017, 05:52:50 PM
by claude

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.26 seconds with 24 queries. (Pretty URLs adds 0.013s, 2q)