Midpoint Displacement in One Dimension

For the creation of the 1 dimensional height map, similar to a horizon, a very simple algorithm is used:

    def generateLine(length)

      # Due to this algorithm being simple for
      # the articles sake, length should be 
      # constrained to the powers of 2. 
      #
      # As we need a midpoint, however, length
      # should be defined as (2^n)+1.

      # Create an array which describes a line.
      # Set the default values to 0American Standard Code for Information Interchange 
      line = Array.new(length, 0)

      # Define the range of the terrain values.
      # For this example we shall take the range
      # of values to be -63 to 63, with the midpoint
      # at 0, out default value. 
      # Range is then 128.
      range = 128;

      # Work out the number of line segments
      # levels there are in the line. As each line 
      # segment level is defined by deviding by two 
      # until length is 1, the number of segments
      # is the log_2 of the length.
      noSegments = Math.log(length-1, 2)

      #Iterate through the levels
      (1 .. noSegments).each{ |level|

        # Work out the line segment length so you
        # can properly address the offset.
        segLength = (length/2**level)

        # Work out the number of line segments
        # within this level
        noSegL = (2**level)-1

        # Iterate through the line segments and 
        # apply your random function.
        (1 .. noSegL).each{ |segOffset|

          # If value is not zero, skip over it, done on a previous
          # level
          if( line[segLength*segOffset] == 0 )

            # Make sure the current value of the line
            # is directly midway between its two parents.
            line[segLength*segOffset] = (line[segLength*(segOffset-1)] + line[segLength*(segOffset+1)])/2

            # Apply random function to value
            line[segLength*segOffset] = randomFunction(line[segLength*segOffset], level, range)

          end
        }
      } 

      return line
    end

Now as you can see the most important part of that algorithm is that which I have purposely missed, randomFunction. This is where you, quite obviously, define how you want your heights defined. A good way is to simply use a bounded random function, where the bounds are defined by the line segment level you are currently within. For example, a function like:

    def randomFunction(value, level, range)

      # Roughness constant 0 < h < 1
      # As h approachs 0, the level dependent
      # bounds grow and tend to 1 for all values
      # of level
      h = 0.8

      # Define bounds in terms of level and a roughness
      # function
      multiplier = 2 ** (-h * level-1)

      # Perform random function, bound it, and then apply 
      # multiplier
      value += (rand() * (range) - (range / 2)) * multiplier

      # Return
      return value
    end

Would offset the value of the line a random amount which is bounded by a function dependent on the line segment level and a roughness coefficient. As this function is the heart of this algorithm and the 2 dimensional algorithm, I will go into some detail on the use of the roughness coefficient. If the roughness coefficient is set to 1, then the multiplier acts the same as : . Decreasing the value of h flattens out the function meaning the bounds are less constrictive and allowing for much rougher terrain. Here is a plot of the multiplier when h=0.8 and h=0.2, and a plot of the generated lines when those constraints are used. X axis is equal to line segment level.

Roughness function

1D heightmap

As you can see, the roughness coefficient makes a massive difference to the outputted numbers. For those interested in recreating the plot, I piped the output of the above code into a text file, and then used GNUPlot to make the images.

Extending into 2 dimensions - The diamond square algorithm

To extend this algorithm into the second dimension, we have to imagine the terrain data as a two dimensional array of numbers. Each number represents a height in this two dimensional field, and so each column and row can be treated similarly to above.

To split up a two dimensional array in a self-similar way we must use squares, or diamonds, which are analogous to the line segments of the 1 dimensional algorithm. Then, rather than using the ends of a line segment to work out a base height the corners of the square, or diamond, are used. For example:

Like the line segment algorithm, the mathematical 'meat' is the same random function as before. The complexity comes in managing which indexes are part of a diamond or a square. So, for example, here is a code segment which works out indexes of a square, depending on level and location, and applies the random function:

    # Get the corners of an arbitrary square and perform operation
    # on center.
    #
    # @param  array           Array to use
    # @param  topleftIndexX   X index of top left of square
    # @param  topleftIndexY   Y index of top left of square
    # @param  length          Length of the square 
    # @param  level           Level into the calculation
    def processSquare(array, topleftIndexX, topleftIndexY, length, level, range, h)

      # Get coordinates of the corners of the square
      toprightIndexX    = topleftIndexX
      toprightIndexY    = topleftIndexY + length + ((level == 0) ? -1 : 0)

      bottomleftIndexX  = topleftIndexX + length + ((level == 0) ? -1 : 0)
      bottomleftIndexY  = topleftIndexY

      bottomrightIndexX = topleftIndexX + length + ((level == 0) ? -1 : 0)
      bottomrightIndexY = topleftIndexY + length + ((level == 0) ? -1 : 0)

      middleX           = topleftIndexX + (length)/2
      middleY           = topleftIndexY + (length)/2

      # Get values
      topleftValue      = array[topleftIndexX][topleftIndexY]
      toprightValue     = array[toprightIndexX][toprightIndexY]
      bottomleftValue   = array[bottomleftIndexX][bottomleftIndexY]
      bottomrightValue  = array[bottomrightIndexX][bottomrightIndexY]

      # Get average
      average = (topleftValue + toprightValue + bottomleftValue + bottomrightValue)/4

      # Set new value
      array[middleX][middleY] = average + calculateOffset(level, range, h)
    end

Where calculateOffset is the random function in this application. The diamond calculation algorithm is very similar and looks like this:

    # Get the edges of an arbitrary diamond and perform operation
    # on center
    #
    # @param  array           Array to use
    # @param  topIndexX       X index of top of diamond
    # @param  topIndexY       Y index of top of diamond
    # @param  length          Length of diamond
    # @param  level           Level into the calculation
    def processDiamond(array, topIndexX, topIndexY, arraylength, level, range, h)

      # Adjust
      arraylength -= 1
      length = arraylength/(2 ** level)

      #Get coordinates of the diamond
      rightIndexX   = topIndexX + length/2
      rightIndexY   = (topIndexY == length) ? length/2 : topIndexY + length/2 

      leftIndexX    = topIndexX + length/2
      leftIndexY    = (topIndexY == 0) ? arraylength - length/2 : topIndexY - length/2

      bottomIndexX  = (topIndexX + length/2 == arraylength) ? length/2 : topIndexX + length
      bottomIndexY  = topIndexY

      middleX       = topIndexX + length/2
      middleY       = topIndexY

      # Get values
      topValue      = array[topIndexX][topIndexY]
      rightValue    = array[rightIndexX][rightIndexY]
      bottomValue   = array[bottomIndexX][bottomIndexY]
      leftValue     = array[leftIndexX][leftIndexY]

      # Get average
      average = (topValue + rightValue + bottomValue + leftValue)/4

      # Set new value
      array[middleX][middleY] = average + calculateOffset(level, range, h)

      # Wraps
      if(middleX == arraylength)
        array[0][middleY] = array[middleX][middleY]
      end
      if(middleY == 0)
        array[middleX][arraylength] = array[middleX][middleY]
      end
    end

The only difference with the above snippet is the different indices it retrieves, and that it must handle wrap around for some of the edges.

So currently, we can create arbitrary diamonds and squares within a 2-dimensional array and assign a fuzzy average of the edges according the h value. Now all we need is some code to manage traversing through the levels of iterations and through the diamonds and squares themselves. Here is my solution:

    # The main control loop for the algorithm.
    #
    # @param  lengthExp       Length exponent
    # @param  range           Value range
    # @param  h               Roughness constant
    def generateTerrain(lengthExp, range, h)

      length = (2 ** lengthExp) + 1

      array = Array.new
      array = createArray(array, length)

      #Go through Levels (irerative recursion)
      (0 .. lengthExp - 1).each{ |level|

        # Iterator for the Square part of the algorithm
        # Will go through the x-axis coords
        (0 .. (2 ** level) -1 ).each { |sqx|

          # Y axis coords
          (0 .. (2 ** level) -1).each { |sqy|

            gap = length/2 ** level
            x = (0 + (gap*sqx))
            y = (0 + (gap*sqy))

            processSquare(array, x, y, gap, level, range, h)
          }
        }

        # Iterator for the diamond part of the algorithm
        (0 ... (2 ** (level+1))).each { |dix|

          # Offset in the number of points on the y-axis. Dependant
          # on if x iteration is even or odd.
          offset = (dix.even?) ? 1 : 2
          (0 ... (2 ** (level+1)/2)).each { |diy|

            gap = (length/2 ** (level+1))
            ygap = 2 * gap

            x = (0 + (gap*dix))
            if (dix.even?)

              y = 0 + (ygap*diy)
            else

              y = gap + (ygap*diy)
            end
            processDiamond(array, x, y, length, level, range, h)
          }
        }
      }
      return array
    end

And this gives us our array with its height map hidden inside. Using a library like RMagick we can output images like the ones shown above. To create the gray scale image, the following code was used:

  image = Image.new(array.length, array.length)

  (0 ... array.length).each { |x|
    (0 ... array[x].length).each { |y|
      val = array[x][y] * (2**9)
      # Create greyscale image
      image.pixel_color(x, y, Pixel.new(val, val, val, val))
    }
  }
  image.display

Which just takes the value in the array, and multiplies it by 512 which gives the values a range of . This gives us the gaseous image that has been generated above.