Problem with my Diamond Square algorithm

[mbob61]

I imagine you are all sick to death with hearing about the dreaded "spikes in my diamond square" problem, but it is sadly a problem I am having.
I've looked around for various solutions and I'm pretty sure my problem stems from my random function but I'm not sure.

All the relevant code is below.

Code:

        private float RandGenerator(Random r, float Max, float Min)
        {
            float number = Min + (float)r.NextDouble() * (Max - Min);          
            return number;
        }

        //private float randomMinus(Random P)
        //{
        //    Random Q = new Random();
        //    float y = (float)P.NextDouble();
        //    float x = (float)Q.NextDouble() * (y < 0.5 ? 1 : -1);

        //    return x;
        //}        

        private float[][] generateGrid(int size, float maxHeight = 40.0f, float minHeight = 0.0f, float noise = 0.0f, int seed = 0)
        {
            int s = (size - 1);          
            float heightRange = 0;

            float[][] grid = new float[size][];          

            for (int i = 0; i < size; i++)
            {
                grid[i] = new float[size];
            }

            gridRandom = (seed == 0 ? new Random() : new Random(seed));

            grid[0][0] = 0;
            grid[s][0] = 0;
            grid[0][s] = 0;
            grid[s][s] = 0;

            float topLeft,
                topRight,
                bottomLeft,
                bottomRight,
                topMiddle,
                leftMiddle,
                rightMiddle,
                bottomMiddle,
                centre;

            for (int i = s; i > 1; i /= 2)
            {
                float divisor = (float)i / (s);
                //Reduces the range each time.
                //The divisor will become a smaller decimnal with every step.
                heightRange = (maxHeight - minHeight) * noise * divisor;

                //Diamond Step
                //This finds the centre of the square. This will be used as one corner of the diamond
                for (int y = 0; y < s; y += i)
                {
                    for (int x = 0; x < s; x += i)
                    {
                        //Assigns grid values for the 4 corners of the square.
                        topLeft = grid[x][y];
                        topRight = grid[x + i][y];
                        bottomLeft = grid[x][y + i];
                        bottomRight = grid[x + i][y + i];

                        //Finds the centre piece by averaging the 4 corner pieces specified above.
                        grid[x + (i / 2)][y + (i / 2)] = ((topLeft + topRight + bottomLeft + bottomRight) / 4) +
                        RandGenerator(gridRandom, -heightRange, heightRange);
                    }
                }

                //Square Step
                //This finds the centre point of the diamonds (or edge values of the square) produced from the diamond step.
                for (int y = 0; y < s; y += i)
                {
                    for (int x = 0; x < s; x += i)
                    {
                        //Assigns grid values for the 4 corners and centre of the square.
                        //These will be used to form the diamonds.
                        topLeft = grid[x][y];
                        topRight = grid[x + i][y];
                        bottomLeft = grid[x][y + i];
                        bottomRight = grid[x + i][y + i];
                        centre = grid[x + (i / 2)][y + (i / 2)];

                        topMiddle = y <= 0 ? (topLeft + topRight + centre) / 3.0f
                                           : (topLeft + topRight + centre + grid[x + (i / 2)][y - (i / 2)]) / 4.0f;

                        leftMiddle = x <= 0 ? (topLeft + bottomLeft + centre) / 3.0f
                                            : (topLeft + bottomLeft + centre + grid[x - (i / 2)][y + (i / 2)]) / 4.0f;

                        rightMiddle = x >= s - i ? (topRight + bottomRight + centre) / 3.0f
                                                 : (topRight + bottomRight + centre + grid[x + i + (i / 2)][y + (i / 2)]) / 4.0f;

                        bottomMiddle = y >= s - i ? (bottomLeft + bottomRight + centre) / 3.0f
                                                  : (bottomLeft + bottomRight + centre + grid[x + (i / 2)][y + i + (i / 2)]) / 4.0f;

                        grid[x + (i / 2)][y] = topMiddle + RandGenerator(gridRandom, -heightRange, heightRange);
                        grid[x][y + (i / 2)] = leftMiddle + RandGenerator(gridRandom, -heightRange, heightRange);
                        grid[x + i][y + (i / 2)] = rightMiddle + RandGenerator(gridRandom, -heightRange, heightRange);
                        grid[x + (i / 2)][y + i] = bottomMiddle + RandGenerator(gridRandom, -heightRange, heightRange);
                    }
                }

            }
            return grid;
        }

THIS IS INSIDE THE UPDATE METHOD

            if (keyboard.GetPressedKeys().Length > 0)
            {
                if (!keyDown)
                {
                    if (keyboard.IsKeyDown(Keys.Up))
                    {
                        roughness++;
                    }
                    if (keyboard.IsKeyDown(Keys.Down))
                    {
                        roughness--;
                    }
                    if (keyboard.IsKeyDown(Keys.Enter))
                    {
                        seed = initialRandom.Next(100);
                    }                    

                    //else if (keyState.IsKeyDown(Keys.Space))
                    //    doRotate = !doRotate;
                    keyDown = true;
                    generateTerrain();                          
                }
            }
            else keyDown = false;
    }
}

http://imgur.com/htXV0

I've been tinkering around with various randoms but can't seem to get the result i should be having. Perhaps my problem is elsewhere?
Hopefully its something really simple I'm missing. Apologies if i've missed something really simple, I'm not a great coder and I'm new to this concept but I'm trying my best

Cheers,  Mike

[fractower]

The diamond square algorithm is a clever way to produce a linear combination of wavelets. Unfortunately the wavelet basis function has a spike. Kronikel made a nice animation of how the basis functions develop in the following thread.


http://www.fractalforums.com/programming/problem-with-diamond-square-algorithm/


The thread is quite long, but you might find it useful.

[mbob61]

Cheers, I'll give that a look.
I had a thought about running through the first few iterations and taking another average, this time using the height of the point as well. Hopefully this should help flatten out certain values if they are too much taller than the rest.
I know this isn't the nicest way but hopefully it might fix the problem.

Mike

[darylj]

I've been working on a code for something like the diamond-square algorithm, and I've read through the thread started by Kronikel, the webpage mentioned there by DarkBeam (http://www.gameprogrammer.com/fractal.html), and looked at some papers by Mandelbrot on fBm and the ones by Fournier et al, Gavin Miller, and a couple of others. I haven't yet digested Mandelbrot's criticism of diamond-square from the technical correspondence section of the August 1982 volume of Communications of the ACM, or looked at the reply by Fournier et al, but I think that's probably worth doing. A paper by Lau et al (1995) on Self-Similar Traffic Generation has a section (2.3) which I think does a really nice job of explaining how the midpoint displacement (i.e., the random offset) calculation is supposed to work---a lot better job than the description on the above website!

My reason for looking at all of this is because I wanted to make sure I was calculating the random offset correctly. And I think the tent peg-effect actually does come from doing this incorrectly, because I'm only getting tent pegs when I use incorrect parameters in my model.

The correct way to calculate the random offset values is stated by Miller on pg 40 of his paper, and shown in more detail in the Lau paper I mentioned. At each iteration, you're supposed to draw randomly from a Normal distribution with mean zero and standard deviation k2^(-iH), where k is your initial scale factor and H is like the Hurst exponent in fBm. Therefore, you initially sample from a Normal distribution with standard deviation k, then you divide by 2^H, and then by 2^H again, etc., at subsequent iterations.

The thing is, H has to be between 0 and 1. With H=0, the standard deviation stays the same at each iteration (division by 1), and you get a surface based on totally uncorrelated random Gaussian noise. With H=0.5, you get a "Brownian" surface (if it's okay to call it that: Mandelbrot may roll in his grave). And with H=1, you divide your standard deviation by 2 each time you halve your length-scale. This surface is strongly correlated, since the random offset at each subsequent midpoint can't deviate from the average by more than half of what it could in the previous iteration.

The problem with the tent pole-effect, I think, is that you might be effectively using a value of H that's greater than 1, and therefore smoothing the surface too much with subsequent iterations. That's why the first set of values dominates the final image: the subsequent points are too close to the average because you haven't introduced enough randomness---and it only gets smoother with each iteration. At least this is what happens when I set H>1 in my code.

I hope this helps. Incidentally, the images that Gavin Miller generated with the diamond-square algorithm (Figure 2) also had the tent pole-effect, but I don't see the value of H he was using anywhere in the paper.

[Tglad]

That sounds like the problem is simply that mbob61 needs to use a normal distribution rather than uniform distribution. He is already using H=1 since standard deviation is proportional to the width of each square (incidentally, it would be better to scale the heightRange by sqrt 2 for the diamond part, since it is operating on a larger width. (Also better to divide by 4.0f than 4)).
mbob61 if you can't get rid of the tent poles, another way to get a similar effect is to sum since waves together, the amplitude of each sine wave is proportional to the wave length, but you need to do include more waves with smaller wavelength... to discuss more. ... Or you could do Veronoi noise (Worley noise) terrain.

[David Makin]

Here's my DS for UF - this does the DS on a pixel-by-pixel basis as necessary for standard UF rendering (and of course using shaders).
I've posted it as the UPR with formulas so I think if you put the whole in a UPR file it should work OK if you then load "Fractal1" from the UPR.
Below the first (at base desktop resolution with just a 1*1 wrapped area) is a second at max resolution (31 depths of DS and max map size 2^31 DS squares).
I was going to convert this to shader code but just never got around to it - I know the UF version has noticeable artifice but most is from the poor random generator as you'll see simply by noting large changes in style as you change the seed.

I don't think this has tent-pole spikes,  but I haven't seen exactly this version in 3D yet

Code:

Fractal1 {
::7N9T+jn2lNVTPqNMUw7Ix/Brcvg/I+ri8lKEHq0eavWpVWJOgXSsTdMtL/77LQMgUvlMjfz8
  evxuLZby2+vveFCl95enp6wdISF6v+28JjmqQnc+jnyGhWg6tXdpJDZuimkr1nnMV7t/xhez
  e2H2RobJktUMhubP6nPBp4ZQyuXPKeLf7vmFqw4bXoEbJc48Y1uDgHuk/M69e/nndpBbYncL
  lNzSrWv6W/cr/bsjZfMYq+ht58xU8SotCFHtN+8VDBjRDu8pYrZ4Sf2PanmWvawOO6DHvXuL
  kdJDdjSxUatuWCN9GKrmj1SoT4oB7xgBvRrlUGjueVXMBaZvV9g9L/c5z+M6SNncNnNxuOUn
  v3FsDwidv3OEDtv/7L2kbzluhq5CBXTXNVZ3U+OTFa8jW3Is5JgUfkshWz8HTnspRweJWSwC
  piQUCOfWCgz5aNU68XNWIExbYc4nO7U2c1BDqPM5bd3jZQxpOoVDxgb9q4l8/T17DObClcjO
  bGO4LTxwQ3mLN9VlG/t3OY/GsJY7hx5Yy26BibqNNEjwUA+jSxs9W2QQ+Qr7LYgai9xkRVLY
  KJdBlyZL4EuST5a2sOFORhjULEUmuQIKixoKtkIK4S8tY9OlUqUyHUQYeHvWp4CcdBXrL4CY
  7ykPtnRqfUiWoxyCO86YBXK4KIyWwrf1dGhpoFXY8iUC4aFR/QKhcBHGcWtC/i7qHLGCk/q6
  HapJPIYEBT/kQ/0fiQKlUCcZZ55wr5DkulMp8YhOfv6f0zqONA==
}

DiamondSquare.ufm:testSquare {
; 0                   1                   2                   3
; b       *1          b        *1         b        *1         b
;
;
;          0          1         2         3         4
;*2       *3         *2        *3        *2        *3        *2
;
; 4                   5                   6                   7
;          5          6         7         8         9
; b       *1          b        *1         b        *1         b
;
;
;         10         11        12        13        14
;*2       *3         *2        *3        *2        *3        *2
;
; 8                   9                  10                  11
;         15         16        17        18        19
; b       *1          b        *1         b        *1         b
;
;
;         20         21        22        23        24
;*2       *3         *2        *3        *2        *3        *2
;
;12                  13                  14                  15
;
; b        *1         b        *1         b        *1         b
::jx0cIgn2FnVbvNytR4vXg+fg5Caskstud1bWXypDOXSuPUg0C0rB9DuuAUa5ql27uUmkrsVQ
  +x3ZGy99135GcB9AObtcHO8Ze7ZGKvPVtln+t/5/EjJztM9Nznsa2tlPKZbYBlPYEiI45rxf
  jrpFHEcL+J4z3IvFenXENPPSlNCfaM++icrMlNSuRee443tZ+kRXHJOYTO/aUSQEZu0SYIOV
  xtskbCXdb9j5J3MbZjnfCOiRBTntciWwTH91HkPJSHP+vEONoWoTVCJz47fGhiBhgFe9X7QM
  8z9i63+A82ZTD+PtQ7lhjLdIHLhJ+Q+xSQiPF10xdXzHuv5DZVPEze6tBOoxID8J25MPWF5R
  yYvUnqlCtwTDIFt7J9AOBbaPPz7aGPTuJCjnVn1R/zNimRbiooZ4KcRvBnqUH+2uKryKdWw7
  28g7hST4SW1CooZgNt0tQlN58QN2pzFNwODHanRvlMzJzLl+uaIxqDJUmcDgj/DS+uHN94Jj
  gQ893u54N3d756biOf2tjbLGsStgwDoonH+cCvopwLcCPbAhRoeHYYLKXy7/vfz9nHGABgZL
  /9ZVA0aggwWwNsHcn3U45tgbfhX2U4lOhXQC3T0VooJw7dFP+3cs8NH7+GcPX1U9X5U/yhV/
  6S1vqr6XXq+VdV/baq+34U/qeq/zHXq8+/P4tfx+aRqR87PbmIGHlgpyJOck4CRuIxYwuexJ
  8V6qUNzKVz5OHfDt9yqKq1oXVzLVF9rrapx//V60ywJP75evSf34Lr+qlhXqLvfZAn5LqM0p
  TH8aFkdKOMolK/DqYtGE9yO8wYWLY8HURdDYsoTKlHGLaBjvsF/fZyX+yntETT1lw1HmE3r7
  ZjCj2VCfyiEkZBePAmsGBij+VP2a1SfZY9Sdam329Ce3FPTXfnY3VOUCjGhLY6yJP0Ywlv66
  Yuxy+mvh9wbDHz+tfjNqaln2sBm5B/0J6TjZ/7aUhSiT9MfwZYpBNIF52xviDgCTgu81Xb2x
  TFjhQCEin4eimYriCvvslxWqMhSSfD+DsatlKcBqtcZqqgGhOmTaNSEzLSdDVbl2UBcCv6Hl
  8MVeE7jPUw1C2ja+hDioXhyAzInT+qrWWOc6BumnxIj0hyd8DWpCF7VxvPj9j4bel7V+jDHn
  2nTkJRBn5fgjzrN3HlSQtjwRELzFGmNRAKwC2ALVcUkykTFTpV1CjKtgOUVMtCcwTbGTw/99
  sj80CBKyM2e5RQlI+WMZhXtXw41i0ZzzbsjHl2E26JrLRDEYZGFTlPl9RlHmO7sjSA1+mu6Z
  Z4sJw/bqrPBKC9b/kqoWF+tOLcxVLWPf1i1Tq/YPvw/Egn6oQzTTL3o3n5ucDEPtWhO3wkGG
  kkoY7U5WtKNFuzWHdt9Etx/hbjUg2Dd0MsJwlj4RRGmSDhoMFCc/RCR9O66g6Rh24COwb5Mj
  MTmy1QSReuQPl9BQLycjlnvT00Jbcu1OqLGk2nI4BGCpwq1bh5dqstQKGzqI3LIavIXbPmL/
  DO3wgZLmg/oSxguOkKsVlN9sTn7dK7nMGBw0AaF9igzOTwBvO4z+lPw+ZoOTGL3xpkaHs5sU
  +JhurChDFdXebIqO3JcBGCF5qi9JMAw5kWOopLpL/VssJPC95dUIgAD/EbdQQwE8HM6OxGyB
  iXMDS9QaB8C4y9JCNRKUD2OaDFHClHArEyFgiEyN79UOETcHsDgmiUFbTl57dpka0T0RdtPL
  0e6ttSDrJtAidEJd02OIS7zzA8pFRF+cLnTEyvR8iuNBfXCLhnesxx85gWPYQ8DQydMfn1HS
  zV6MeK7Rl++u6zvdqYY1igJLWHA17Qa7juyew5UdkGfeiPZorZW6pdp7UeLb0mNVZvjxMl9K
  VUZ6CBud7K0croPuOAmEk7yLZchHpCxHxTv29B6UmvTL4G4ohgOoQF0cdoSVfmUPP2rrVQD3
  L5K/7xgdVoNiBKIYpIbC9JykvwxFAuJ/quDaK7HS457xQKu47l9wm3hljZvpWsYZ4+N+q6Lg
  4plWAmT6g4iuqDYThKoH5QdPcming8AQzlMdluRoyCaXHLAQabx4YGs2fXTbomvD8+bFWnBz
  RpqFyT+NotaToEAi8G8vFpRuuAbBCXsYoMGrFnZArF0rSurXMIVEbd9njpCBowCPcvxQ6oFS
  nxI2dKWjC2RdP7+mHO9V+Bco5QaPUC2Rr3MJNZ+7NWSrpSC+UTl8xDidy4TIYpsdD7dBwxUY
  8Tq47myT3r0Q1aGVAEFNExhhUlEyxzLy2CBMIpar0eJxLaQGbs3ifyC3sZXaczmlJg4VveMc
  M0jITzf0DEzUozP4wQGPIUWUWSCIKTpxKWt0VPDrOQBVZ3M077Objroi2t3Y9Au2eo8zB0W1
  RKsnIe59CbzKaKQXYa7cfOl9c9mbCbPkbTJhNFGgLSBpKZFZeOCEDICzFOPGMNkRG5pRMgYw
  JhEHuUyuBC4EaNpBRcDAB8QWUbt7ZgiTNb7T5WR4QV7CIc+TYHpWKPBi7pCjh8IJCeUfvFE9
  pOZ1IGy0BOiMjjHnofBzHKc/bKrwPDXqnUobblukn+BtSReiIlwQsh4X+OQy5wdTyyesZjUI
  viDlwnkNbEMu9p1xpmOkHDroBrDwAQQx3C3tx77Krl84Abz3g6ZgC06iDup/x3Ln6IceXChA
  DvP5k7PZgjeiuLdP+pPirC33x0nca6VBXFGs6q1hhrXtcx6BIrKpnwx8heK79JquGy2y+znR
  HtQfWZPjLQP1eROWJ136pKwYNYxBTXgevgprdjssVYsfKG4y/IQDxA/Rh/SltuX4s22xvYwx
  xP5NAsQMVeP0YlyipA/BtiGbj3opZ1Q2dMl3LIG6S2hcutQX9EMPlnHz5IsAXgRlJKdAd0Vr
  DvKlmKwVbPKVFY6f+l9uVFFB6f70/ahBYE0nKPN33CBdPPYmP3czQx4ZUhKc7I7pUn/QLeoQ
  q7lzSNfbsLiyF9dsDoNj84IXnXAi7gM2Oqp2EwUSE9S7nJru8LdodWNuavo9ipdSb/hutnuj
  cLQ0ZvmHJxYL4NEPJQ0irfiuLCUtb1yd2+X3yXqT/J7YB4jhOXT52hhaOzS9JYaY60oWH20+
  zI5K0bOlbMErwR6gDvCkUfHkNBIHU4c0YfkBvLufyX6aSCnpeWFFwZ1lsVxPgLt6Y6mS24UB
  LlO4+BrtKVqPWV/1Q1IU9B8L5iyj5aL2cqfczqLEtDc/LEnickmNqxE2GXZCESzAyIZO1KwP
  Hu77IQgfDTDQgDbxkC+GcAZy2bwjXiLfxUdxF3OE9NeZAjdw5DMiUBl2UPEBQp31n9dsvHmT
  yUGVck8YPBEM1LPC6ZBA+C4AdN0vwpuxoGoyCqhIOJBciA23Cdc2B8O/MKQ/ionolbxOi1Lf
  Xbo1ffneB77/CMnqCoL=
}

mmf.ucl:MMFa-for3D {
;
; Probably the most complex colouring ever !
; This colouring only really applies to the
; MMF "3D" fractal types and will not produce
; very interesting results for other fractals.
;
::cqgJBhn2LtMzLxcsiXuUQBlzMvUStCFsVhiSNxc0Q5q0kXuSJ10Ss0cKBs0lkZJ5kKQZVyXf
  dTBjdRhkzPn8LtoMzLdl4lLA6SwEWD==
}

Code:

Fractal2 {
fractal:
  title="Fractal2" width=928 height=696 layers=1
  credits="Dave Makin;12/11/2012;D J Makin;12/20/2011;Dave Makin;10/5/\
  2011;David Makin;6/15/2008;Frederik Slijkerman;7/23/2002"
layer:
  caption="Background" opacity=100 method=multipass
mapping:
  center=2.883899937206161385/2.2345097056422888 magn=4.8116135E8
formula:
  maxiter=100 percheck=off filename="DiamondSquare.ufm"
  entry="testSquare" p_depth=31 p_rand=31 p_sharp=0.70710678118655
  p_seed=22 p_scale=0.35 p_fast=yes
inside:
  transfer=none
outside:
  transfer=linear repeat=no filename="mmf.ucl" entry="MMFa-for3D"
gradient:
  smooth=yes rotation=1 index=0 color=8463872 index=253 color=15892593
  index=256 color=11466239 index=260 color=3289716 index=270
  color=3778876 index=285 color=4885604 index=299 color=4618637
  index=314 color=4896907 index=328 color=4765810 index=340
  color=3313824 index=354 color=6008197 index=367 color=5933480
  index=383 color=11707844 index=391 color=11316394 index=399
  color=16777215
opacity:
  smooth=no index=0 opacity=255
}