Fractal images are strikingly beautiful, bright, and complex and they're everywhere these days - from art galleries to image compression to movie graphics to environmental modelling (see Figure 1). It really is amazing to see a concept from the dusty world of mathematics jump up and grab the world's attention. Fractals and associated concepts came into being just at the time that powerful computers became cheap and widely accessible, and this is no coincidence. In the past it was widely assumed that complex phenomena such as coastlines, trees, and mountain shapes could only be described by complex and difficult mathematics, which could only be understood by a select few.

This view has now been turned on its head and we can use simple equations to generate immense complexity just like that seen in nature. We can generate this complexity by using a computer to churn through millions of simple computations, hence the relationship between the emergences of fractals and computing. The reason I like fractals so much is that, together with chaos theory and nonlinear dynamics, they have brought mathematics into the "cool science" domain.

Another reason I like fractals is that although they look complex, they are simple to generate (well, simple for a computer!). When equipped with a computer and some programming skills, this new branch of mathematics is virgin territory that can be explored in an experimental and visual manner. I find that to be an exciting prospect.

In this article I will show you how simple it is to generate images of fractal landscapes complete with coastlines, sea, and relief shading, just like those illustrated in Figure 1. Fortunately the GDI+ library of the .NET framework makes it equally simple to render these images in plan view as bitmaps. Ideally these should be presented in 3D but that would introduce another level of complexity that is beyond the scope of the current article.

Generating Fractal Landscapes Making Random Numbers
First of all the .NET framework provides a Random class specifically for generating random numbers; however, for the current example we require the numbers to be distributed according to the normal (hopefully familiar bell curve) distribution, whereas Random gives us uniform or evenly distributed random numbers (i.e., each number is equally likely to occur). In a normal distribution numbers near a specified average occur more often than numbers away from the average. A normal distribution is characterised by a mean (i.e., average) and a variance, which describes how far away from the mean values may occur (see Figure 2).

For example, a large variance results in data spread widely around the average, whereas a small variance gives data tightly clustered around the average. In order to implement normally distributed random numbers in code I have created a new class called normal, which is derived from the .NET Random class, as this saves a whole lot of work. It has two private data members for mean and variance implemented as properties, which can be specified in the constructor. The base class constructor of Random is called and its Random number generator is seeded using the current date and time. Finally I have added a method called NextNormal, which gives normally distributed Random numbers by slightly modifying Random's NextDouble method.

Random Displacements
Much of the algorithm for generating a fractal landscape involves making random adjusts to the height of points on a square grid. If the random displacements were always of the same order of magnitude then as square grid became smaller and smaller, the displacements would become large in comparison to the individual squares and leading to a very jagged and unrealistic landscape. We counteract this by scaling the random displacements according to size of the smallest square. Of course we could try any strategy we like to achieve this and in doing so many weird and wonderful landscapes can be generated. The usual approach however is to generate random displacements from a normal distribution with a mean of 0 and a variance of (2h²)a where h is the side length of the square and a is a value between 0 and 1. Therefore, as h gets smaller the variance gets smaller giving proportioned random displacements. The a parameter additionally allows us to control the jaggedness of the landscape closer to 0 means jagged and closer to 1 means smoother. The fractal dimension (D) of the landscape is calculated from a as D = 3 - a.

Landscape Algorithm
The algorithm is relatively simple and consists of taking a large square and filling in the details by progressively quartering the square and adjusting heights randomly.

I have decided to take an approach that makes the coding less taxing and more efficient. Usually we would have to deal with three space coordinates (x,y,z), however if we use a square array of floats to represent the landscape then, by choosing a suitable size, the indices of the array can correspond to the (x,y) coordinates and the value stored at each array location represents the height (z). For example, suppose that I sized my array to be 9 by 9 (= 2³ + 1), then I can continuously half my array to get integer indices i.e., half of (0,8) is 4, half of (0, 4) is 2, half of (4, 8) is 6, etc. Therefore choosing my array to be of size [2ⁿ + 1, 2ⁿ + 1] where n is a positive whole number makes halving squares painless.

For simplicity I have chosen n = 9, so my array is 513 by 513 and ultimately a bitmap of size 512 by 512 will be generated. A variable called level is used to orchestrate the level of detail which needs to be inserted at any given stage and can take values between 1 and n (i.e. 9 right now). Using the level I can calculate that the current square is of length 2ⁿ/2^level-1 and that the quartered square is of length 2n/2level. For example, at level 1 the current square is of length 512 and the quartered square is of length 256.

To begin with I set the level to 1 and give values to the vertices of the current square (0,0), (0,512), (512,512), (512,0) of 0 plus a random displacement with variance (2h²)a where h is 2ⁿ/2^level-1, the current square length. Two functions are then called iteratively as the level is increased by one each time to fill in the details of the landscape.

This procedure is illustrated in Figure 3. To begin with we have a square grid with known elevations previously calculated for the vertices (i.e., the grey dots in Figure 3(a)). The center points (i.e., the black dots in Figure 3(b)) can be calculated as the average elevation of the four adjacent vertices (grey dots) plus a random displacement with variance (2h²)a where h is 2ⁿ/2^level, the quartered square length. Finally elevations for the midpoints of each side (i.e., the black dots in Figure 3(c)) can be calculated using the average elevation of the four adjacent grey dots plus a random displacement with variance (2h²)a as before.