Hello All!

Today we’re going to implement a very popular algorithm, the Diamond Square algorithm. Simplifying all the theory about, it’s a fractal algorithm, because it uses recursive behavior (depends on results previously calculated same way) and random elements. Note that was in middle-80’s when computers were enough powerful (at least for the street-man) to be used for study themes like fractals, Chaos Theory, Quantum Physics and all Mathematica involving random values. So, begin to appear nice algorithms (and nice results as well) applied to games, terrain generation, flows simulation, aerodynamics, etc.

Ok, let’s go: Diamond Square algorithm it’s based in the Middle-point Displacement algorithm. This algorithm is so easy as, from a set of points geometrically located, create a new point in set’s geometrical middle with the average value of the set, plus a random value based in the difference between values of points. Yet the Middle-point Displacement it’s good enough, as the set it’s always squared, it generates ‘square-shaped’ random artifacts. A solution for that, at 2D implementation, it’s an intermediate variation calculating at same iteration the top-down-left-right points as well, as shown in diagram :

Seems very easy to implement, but there are some factors that can block the first attempts. Some of them :

• Recursive implementations sometimes uses values that are yet not calculated, so new points are calculated erroneous.
• Border points at the diamond got 3, not 4, points to get the average. As this factor can allow wrapping of the result, to not control it can as well give incorrect results.
• To ensure existence of a midpoint to calculate until end of iterations, number of points must be a power of 2, plus 1. So, there must be exactly 3, 5, 9, 17, …, 513, … points by side.
• Random value to add must be in the order of the iteration. You can add some ‘disorder’ raising this value, but raising it too much can give you a complete random result, with no relation between points.

Coding. First of all, we must set the initial conditions. There will be initial random values for the corners of our squared array W.

      private const int _HIGHRANGE = 1000;  // max high difference will be this value * Entrophy    

    // prepare work array
    private void Initialize(int[,] W)
    {
        // Let's fill in corners
        W[ 0,0] = rnd.Next() % _HIGHRANGE;
        W [_MAXINDEX, 0] = rnd.Next() % _HIGHRANGE;
        W [0, _MAXINDEX] = rnd.Next() % _HIGHRANGE;
        W [_MAXINDEX, _MAXINDEX] = rnd.Next() % _HIGHRANGE;
    }

As did it, we prepare a function to give us a value between (-range,range), and adding some disorder (Entrophy S)

     private Random rnd = new Random();    // generic random object

    // gets a random int in range * entrophy
    public int r(int range,float S)
    {
        return (int)((rnd.NextDouble()*2 - 1)*range*S);
    }

 

Now, there’s no more than the main loop itself. As a fractal got a recursive nature, it can be realized as a recursive function, but anyway it is more clear to realize and explain it via an iterative way: we will run until all posible iterations (this will give us the positions of all midpoints in all cases, that it’s, all 2^n points). At each iteration, we calculate the average of the corners and put them in the middlepoint (diamond step). At each iteration as well, we will calculate the top-down-left-right points (square step). To avoid index-out-of-array and as well incorrect average problems, we will search for the values to average and the denominator to use at each point as well.

    
  private const int _MAXINDEX = 512;      // Map will go from 0..MAX_INDEX (included). So number of cells 
                                           // will be MAX_INDEX+1. Must be 2^n to ensure midpoint at each iteration. 

 // main loop of diamond square algorithm (S is entrophy)
    private void DoDiamondSquare(int[,] W, float S)
    {
        int hs, x, y;
        int A, B, C, D, M, n;

        // Lets iterate through side size until size is too small
        for (int it = _MAXINDEX; it > 1; it/=2) 
        {
            hs = it /2;
            
            //Midpoints
            for (y = hs; y < _MAXINDEX; y += it)
                for (x = hs; x < _MAXINDEX; x += it) 
                {
                    // each square corner
                    A = W [x - hs, y - hs];
                    B = W [x - hs, y + hs];
                    C = W [x + hs, y - hs];
                    D = W [x + hs, y + hs];

                    // center point will be average plus a random in-range value
                    W [x, y] = ((A + B + C + D) / 4) + r (hs,S);
                }

            // Going through each square point                
            for (y = 0; y < _MAXINDEX+1; y += hs)
                for (x = y % it == 0? hs : 0; x < _MAXINDEX+1; x += it) // getting offset of x in function of y 
                {      
                    M = n = 0; // Sum and denominator

                    // this way we can calculate border points
                    try { M += W [x + hs, y]; n++;} catch(Exception){}
                    try { M += W [x - hs, y]; n++;} catch(Exception){}
                    try { M += W [x, y + hs]; n++;} catch(Exception){}
                    try { M += W [x, y - hs]; n++;} catch(Exception){}

                    // lets average sum plus random value
                    W [x, y] = M / n + r (hs,S)/2;
                }
        }
    }

And … that’s all! Really, if you pass to that functions the work array W (size must be 2^n+1, remember) you will get a fractal map of the initial values.

The most of the fractals got a very big resemblance with Nature phenomena. This algorithm it’s, in fact, one of the most popular ways to generate random terrains, understanding result as a high map. Anyway, it can be as well be used as an air density map, or used to search for nature-looking streams. We will show this mapping our map to different color palettes:

 

    private const int _MAXPALETTES = 5;   // palette set number
    public static Color[][] Palette = new Color[_MAXPALETTES][]; // palette set

      // here we generate the palette set
     private void CreatePalettes()
     {
         int i;
 
         // Moduled gray
         Palette[0] = new Color[_HIGHRANGE];
         for (i=0;i<1000;i++)
             Palette[0][i] = Color.FromArgb(0xff,i%255,i%255,i%255);
         
         // Earth like
         Palette[1] = new Color[_HIGHRANGE];
         for (i=0;i<100;i++) // deep sea
             Palette[1][i] = Color.FromArgb(0x0f,0,0,rnd.Next()%40+5);
 
         for(i=100;i<300;i++) // sea
             Palette[1][i] = Color.FromArgb(0x0f,0,0,i/2-50);
         
         for(i=300;i<500;i++) // country-side
             Palette[1][i] = Color.FromArgb(0x0f,i/2-150,i-300,0);
 
         for(i=500;i<700;i++) // forest
             Palette[1][i] = Color.FromArgb(0x0f,0,355-i/2,0);
 
         for (i=700;i<950;i++) // mountain
             Palette[1][i] = Color.FromArgb(0x0f,i-700,i-700,i-700);
 
         for (i=950;i<1000;i++) // high mountain
             Palette[1][i] = Color.FromArgb(0xff,i-800,i-800,i-800);
 
         // Sky
         Palette[2] = new Color[_HIGHRANGE];
         for (i=0;i<1000;i++)
             Palette[2][i] = Color.FromArgb(0xff,(i/4)+5,(i/4)+5,155);
 
         // Volcano
         Palette[3] = new Color[_HIGHRANGE];
         for (i=0;i<1000;i++)
             Palette[3][i] = Color.Black;
 
         for (i=0;i<5;i++) // some lava rivers...
         {
             Palette[3][1 + 150*i] = Palette[3][13 + 150*i] = Color.DarkRed;
             Palette[3][2 + 150*i] = Palette[3][12 + 150*i] = Color.Red;
             Palette[3][3 + 150*i] = Palette[3][11 + 150*i] = Color.Red;
             Palette[3][4 + 150*i] = Palette[3][10 + 150*i] = Color.Orange;
             Palette[3][5 + 150*i] = Palette[3][9  + 150*i] = Color.Yellow;
             Palette[3][6 + 150*i] = Palette[3][8  + 150*i] = Color.LightYellow;
             Palette[3][7 + 150*i] = Color.White;
         }
 
         // Storm
         Palette[4] = new Color[_HIGHRANGE];
         for (i=0;i<1000;i++)
             Palette[4][i] = Color.FromArgb(0xff,30,30,(i/5)+50);
 
         for (i=0;i<4;i++) // some lightings...
         {
             Palette[4][1 + 150*i] = Palette[4][13 + 150*i] = Color.DeepSkyBlue;
             Palette[4][2 + 150*i] = Palette[4][12 + 150*i] = Color.SkyBlue;
             Palette[4][3 + 150*i] = Palette[4][11 + 150*i] = Color.Blue;
             Palette[4][4 + 150*i] = Palette[4][10 + 150*i] = Color.LightSteelBlue;
             Palette[4][5 + 150*i] = Palette[4][9  + 150*i] = Color.LightSteelBlue;
             Palette[4][6 + 150*i] = Palette[4][8  + 150*i] = Color.LightCyan;
             Palette[4][7 + 150*i] = Color.White;
         }
     } 

And result will be :

 

I don’t know you, but I get yet amazed for that images!

You got an online version of this code here (diamond square algorithm) (PD: images got a lower resolution due server optimization work, but you can use any image format)

Hope you had enjoyed. In next posts we will raise up that high maps & raytrace them, and so, we will as well try to navigate through them…

Here full code:

DiamondSquareAlgorithm.aspx

 <%@ Page Language="C#" Inherits="OnlineRepository.DiamondSquareAlgorithm" %>
<!DOCTYPE html>
<html>
<head runat="server">
    <title>DiamondSquareAlgorithm</title>
    <style type="text/css">
        p{text-align:center;font-variants: small-caps; font-name=Lucida Console; color:white;}
    </style>
</head>
<body bgcolor="#00000F" >
    <form id="form1" runat="server">
    <asp:ScriptManager runat="server" id="ScriptManager1">
    </asp:ScriptManager>
    <p>DIAMONDSQUARE ALGORITHM</p>
    <table align="center" >
    <tr>
    <td><asp:Label runat="server" AssociatedControlID="tbEntrophy" Text="Entrophy (0-1 for controled high limits):" ForeColor="white" Font-size="small" /></td>
    <td><asp:TextBox id="tbEntrophy" runat="server" Text="1" Columns="6" Style="width:150pts;text-align: right" ToolTip="Entrophy (bigger, more disorded &amp; bumped)" /></td>
    </tr>
    <tr>
    <td><asp:Label runat="server" AssociatedControlID="ddPalette" Text="Applied palette" ForeColor="white" Font-size="small" /></td>
    <td><asp:DropDownList id="ddPalette" runat="server" >
        <asp:ListItem Text="Gray module" Value="0" Selected="true"></asp:ListItem>
        <asp:ListItem Text="Earth" Value="1" Selected="false"></asp:ListItem>
        <asp:ListItem Text="Sky" Value="2" Selected="false"></asp:ListItem>        
        <asp:ListItem Text="Volcano" Value="3" Selected="false"></asp:ListItem>
        <asp:ListItem Text="Storm" Value="4" Selected="false"></asp:ListItem>
        </asp:DropDownList> </td>
    </tr>
    </table>
    <asp:UpdatePanel runat="server" id="UpdatePanel1">
        <ContentTemplate>
            <p>
            <asp:Button id="btGenerate" runat="server" Text="Generate field" OnClick="btGenerate_Click" /><br/>
            <asp:Label runat="server" Text="Any image yet generated" id="InfoLabel" font-size="small" font-name="verdana" ></asp:Label><br/>
            <asp:Image runat="server" id="UHighMap" ImageAlign="AbsMiddle" Width="513" Height="513" BorderStyle="None" />
            </p>
        </ContentTemplate>
    </asp:UpdatePanel>
    </form>
</body>
</html>

DiamondSquareAlgorithm.aspx.cs

 public partial class DiamondSquareAlgorithm : System.Web.UI.Page
{
    private const int _MAXINDEX = 512;      // Map will go from 0..MAX_INDEX (included). So number of cells 
                                          // will be MAX_INDEX+1. Must be 2^n to ensure midpoint at each iteration.
                
    private const int _HIGHRANGE = 1000;  // max high difference will be this value * Entrophy

    private Random rnd = new Random();    // generic random object

    private const int _MAXPALETTES = 5;   // palette set number
    public static Color[][] Palette = new Color[_MAXPALETTES][]; // palette set

    // gets a random int in range * entrophy
    public int r(int range,float S)
    {
        return (int)((rnd.NextDouble()*2 - 1)*range*S);
    }

    // main load event. Let's create palettes
    public void Page_Load(object sender, EventArgs args)
    {
        // First time page load, lets generate palettes  
        if (!Page.IsPostBack) 
            CreatePalettes ();
    }

    // here we generate the palette set
    private void CreatePalettes()
    {
        int i;

        // Moduled gray
        Palette[0] = new Color[_HIGHRANGE];
        for (i=0;i<1000;i++)
            Palette[0][i] = Color.FromArgb(0xff,i%255,i%255,i%255);
        
        // Earth like
        Palette[1] = new Color[_HIGHRANGE];
        for (i=0;i<100;i++) // deep sea
            Palette[1][i] = Color.FromArgb(0x0f,0,0,rnd.Next()%40+5);

        for(i=100;i<300;i++) // sea
            Palette[1][i] = Color.FromArgb(0x0f,0,0,i/2-50);
        
        for(i=300;i<500;i++) // country-side
            Palette[1][i] = Color.FromArgb(0x0f,i/2-150,i-300,0);

        for(i=500;i<700;i++) // forest
            Palette[1][i] = Color.FromArgb(0x0f,0,355-i/2,0);

        for (i=700;i<950;i++) // mountain
            Palette[1][i] = Color.FromArgb(0x0f,i-700,i-700,i-700);

        for (i=950;i<1000;i++) // high mountain
            Palette[1][i] = Color.FromArgb(0xff,i-800,i-800,i-800);

        // Sky
        Palette[2] = new Color[_HIGHRANGE];
        for (i=0;i<1000;i++)
            Palette[2][i] = Color.FromArgb(0xff,(i/4)+5,(i/4)+5,155);

        // Volcano
        Palette[3] = new Color[_HIGHRANGE];
        for (i=0;i<1000;i++)
            Palette[3][i] = Color.Black;

        for (i=0;i<5;i++) // some lava rivers...
        {
            Palette[3][1 + 150*i] = Palette[3][13 + 150*i] = Color.DarkRed;
            Palette[3][2 + 150*i] = Palette[3][12 + 150*i] = Color.Red;
            Palette[3][3 + 150*i] = Palette[3][11 + 150*i] = Color.Red;
            Palette[3][4 + 150*i] = Palette[3][10 + 150*i] = Color.Orange;
            Palette[3][5 + 150*i] = Palette[3][9  + 150*i] = Color.Yellow;
            Palette[3][6 + 150*i] = Palette[3][8  + 150*i] = Color.LightYellow;
            Palette[3][7 + 150*i] = Color.White;
        }

        // Storm
        Palette[4] = new Color[_HIGHRANGE];
        for (i=0;i<1000;i++)
            Palette[4][i] = Color.FromArgb(0xff,30,30,(i/5)+50);

        for (i=0;i<4;i++) // some lightings...
        {
            Palette[4][1 + 150*i] = Palette[4][13 + 150*i] = Color.DeepSkyBlue;
            Palette[4][2 + 150*i] = Palette[4][12 + 150*i] = Color.SkyBlue;
            Palette[4][3 + 150*i] = Palette[4][11 + 150*i] = Color.Blue;
            Palette[4][4 + 150*i] = Palette[4][10 + 150*i] = Color.LightSteelBlue;
            Palette[4][5 + 150*i] = Palette[4][9  + 150*i] = Color.LightSteelBlue;
            Palette[4][6 + 150*i] = Palette[4][8  + 150*i] = Color.LightCyan;
            Palette[4][7 + 150*i] = Color.White;
        }
    }

    // prepare work array
    private void Initialize(int[,] W)
    {
        // Let's fill in corners
        W[ 0,0] = rnd.Next() % _HIGHRANGE;
        W [_MAXINDEX, 0] = rnd.Next() % _HIGHRANGE;
        W [0, _MAXINDEX] = rnd.Next() % _HIGHRANGE;
        W [_MAXINDEX, _MAXINDEX] = rnd.Next() % _HIGHRANGE;
    }

    // main loop of diamond square algorithm, (S is entrophy)
    private void DoDiamondSquare(int[,] W, float S)
    {
        int hs, x, y;
        int A, B, C, D, M, n;

        // Lets iterate through side size until size is too small
        for (int it = _MAXINDEX; it > 1; it/=2) 
        {
            hs = it /2;
            
            //Midpoints
            for (y = hs; y < _MAXINDEX; y += it)
                for (x = hs; x < _MAXINDEX; x += it) 
                {
                    // each square corner
                    A = W [x - hs, y - hs];
                    B = W [x - hs, y + hs];
                    C = W [x + hs, y - hs];
                    D = W [x + hs, y + hs];

                    // center point will be average plus a random in-range value
                    W [x, y] = ((A + B + C + D) / 4) + r (hs,S);
                }

            // Going through each diamond point                
            for (y = 0; y < _MAXINDEX+1; y += hs)
                for (x = y % it == 0? hs : 0; x < _MAXINDEX+1; x += it) // getting offset of x in function of y 
                {      
                    M = n = 0; // Sum and denominatior

                    // this way we can calculate border points
                    try { M += W [x + hs, y]; n++;} catch(Exception){}
                    try { M += W [x - hs, y]; n++;} catch(Exception){}
                    try { M += W [x, y + hs]; n++;} catch(Exception){}
                    try { M += W [x, y - hs]; n++;} catch(Exception){}

                    // lets average sum plus random value
                    W [x, y] = M / n + r (hs,S)/2;
                }
        }
    }

    // dumping work array to bitmap, using palette P. We mark under-zero and over-highrange values
    private void DumpToBitmap(Color[] P, int[,] W, Bitmap B)
    {
        int x, y;
        Color c;

        for (x = 0; x < _MAXINDEX + 1; x++)
            for (y = 0; y < _MAXINDEX + 1; y++) 
            {
                if (W [x, y] == 0)            c = Color.Red;
                else
                if (W[x,y] < 0)                c = Color.DarkRed;
                else
                if (W[x,y] > _HIGHRANGE-1)    c = Color.Yellow;
                else                         c = P[W[x,y]];

                B.SetPixel (x, y, c);
            }
    }

    // functional to send dump bitmap and send it to web
    private void DumpToBitmapAndSendToWeb(int[,] W)
    {
        // Create bitmap, and fill it 
        Bitmap B = new Bitmap (_MAXINDEX+1, _MAXINDEX+1,System.Drawing.Imaging.PixelFormat.Format24bppRgb);
        DumpToBitmap (Palette[Int32.Parse(ddPalette.SelectedItem.Value)],W, B);

        using (MemoryStream m = new MemoryStream ()) 
        {
            B.Save (m, System.Drawing.Imaging.ImageFormat.Gif);
            UHighMap.ImageUrl = "data:image/Gif;base64," + Convert.ToBase64String (m.ToArray ());
        }
    }

    // an info search...
    private void SearchMinMax(int[,] W)
    {
        int  min=_HIGHRANGE    
            ,max=-_HIGHRANGE;

        for (int y=0;y<_MAXINDEX+1;y++)
        for (int x=0;x<_MAXINDEX+1;x++)
        {
            min = W[x,y] < min ? W[x,y] : min;
            max = W[x,y] > max ? W[x,y] : max;
        }

        InfoLabel.Text = "min/max:" + min.ToString() + "/" + max.ToString();
    }
    
    // let's do a field...
    public void btGenerate_Click(object sender, EventArgs args)
    {
        Stopwatch sw = Stopwatch.StartNew();

        // Work array
        int[,] W = new int[_MAXINDEX+1,_MAXINDEX+1];

        // Create fractal
        Initialize(W);
        DoDiamondSquare (W,float.Parse(tbEntrophy.Text));

        // Dump result
        DumpToBitmapAndSendToWeb(W);
        SearchMinMax(W);
    
        // Finish
        sw.Stop ();
        InfoLabel.Text += " / total ellapsed time : " + sw.ElapsedMilliseconds.ToString() + "ms";
    }        
}