Contents
The purpose of this page is not so much to provide a mathematically complete and rigorous account of the techniques used for fractal landscape generation, but rather to provide an overview of these techniques and their application that is accessible to someone with high-school level maths. Those who have a more complete understanding of this field will notice that there are some complexities that are glossed over. If you notice one of these, please don't bother to email a correction to me - it's an intentional simplification for the purposes of clarity rather than a mistake.
The Fractal Geometry of Nature
In conventional geometry, shapes are defined by simple formulae. Thus, an ellipse might be expressed as a function of radius and angle in polar coordinates. Such formulae are simple mappings - we put in a value of the angle in, we get a value of the radius out. We put a value for the radius in, we get 4 angles out in the range 0 -> 2p.
For a fractal object, this ceases to be the case. The key to fractals is that they have recursive definitions. In other words this means that, for example, the position of a each point on a fractal path is given by the sum of an infinite series. In most cases, this series has no analytic solution, so all we can hope to do is to converge on the correct solution by calculating successive terms in the series. An example of such a series is that which describes the Mandelbrot Set:
Zn+1 = Zn + C
Z and C are complex numbers. If we begin with Z0 = 0 + 0i and C = x + y i, and we calculate successive terms, |Z| (the size of Z) either becomes larger and rushes towards infinity, or it stays close to zero, trapped in a so-called basin of attraction. If it were possible to continue to n = INF, we would either find |ZINF| to be infinite, in which case the initial point, C, lies outside the set, or to be finite, in which case C lies inside the set.
Another sort of fractal is best described in geometric terms, using a rule of replacement. Such an example is the Koch Snowflake. It is described by a procedure:
If we continue this for about 5 iterations, we get the following shape:
In the case of the Koch Snowflake (the above shows just the top - the full figure does resemble a snowflake) we have a recursive geometric description of the set, with each new approximation to the fractal a little closer to the final, infinitely complex shape.
To help us describe fractals, we use a measure called 'fractal dimension'. This is a measure of the increase in complexity at each iteration. In the case of a one-dimensional example such as the Koch Snowflake (where the fractal object is a perimeter rather than a surface), the fractal dimension is between 1.0 and 2.0. It is, if you like, a measure of roughness of the surface. We calculate it using the following formula:
Fractal Dimension = log(Length iteration n+1) / Log(Length iteration n)
The above replacement rule replaces a section of perimeter with a section 4/3 times the length. This gives us log 4 / log 3 = 1.26 as the fractal dimension of the curve.
If the fractal object is a surface, and it is displaced in a third dimension (i.e.: like a fractal landscape), then the fractal dimension is between 2.0 and 3.0.
Now we have the basic terms with which to describe fractals, we can examine some natural processes. Benoit Mandelbrot, beginning in the 1960s, observed and recorded the properties of many natural phenomena, from turbulence in fluids to the stock market. After twenty years of conducting investigations into the mathematics governing these traditionally intractable problems, he published a book called The Fractal Geometry of Nature, which describes how many natural phenomena follow the same sort of pattern as that observed in fractal functions.
Nature is not, however, fractal in the strict mathematical sense. No two trees are alike, and nor are any two mountains, but there is only one Koch Snowflake, one Mandelbrot set and so on. The best we can do is to say that natural processes have some characteristics which are similar to those of fractals. Once we have made such a connection, we can then use our measure of fractal dimension to describe the 'approximate fractal dimension' of a natural object.
Simulating Natural Objects
If we know that a particular type of object is approximately fractal in nature, and we wish to model it, we need to generate some parameterised model, incorporating features present in the real object (in particular, some sort of controlled randomness), without losing the controllability of the fractal model. The way we do this is using a 'random fractal' process. Instead of using a set geometric rule of replacement as we did with the Koch Snowflake, we generate different rules of replacement based on the values of the model's parameters, each of which differs from the others in some controlled random way.
This is exactly the process that we use to generate fractal landscapes. As soon as we introduce randomness into the equation, we can no longer talk about fractal dimension as an analytical property of the system. Instead, it becomes a statistical property, with a mean and a variance. However, we usually save ourselves the trouble, and quote mean fractal dimension.
Subdivision and Displacement
There are several schemes which can be used to inject randomness into a fractal process, but fractal behaviour is not the only requirement for a believable landscape. Indeed, there is rather poor correlation between pure fractal processes and real landscapes. This means that it is necessary to post-process and possibly combine different fractals to achieve the desired result. Anyone who has ever used a photo-editing program will know that the higher the quality of the source, the greater the range of effects that can be achieved, and the same is true of fractals.
A 'high-quality' source fractal process is one which is stationary (all locations have the same statistical properties) and isometric (has the same properties in all directions). We could achieve both of these using a technique called Fourier synthesis, but this is extremely slow compared to other methods of generation. The reason why Fourier synthesis has such ideal properties is that the frequency spectrum of a landscape generated using this method is continuously variable - which is certainly true of natural ones as well. There is, however, almost no difference between a landscape generated this way, and one generated using a technique called recursive subdivision and displacement with cubic interpolation. The largest difference is a discrepancy in generation time of 20 times or more. With this in mind, I will not examine Fourier synthesis here.
The process of recursive subdivision and displacement is as follows:
We express the landscape as a grid of heights. For each 4 grid points, we work out the height of the points between them (interpolation), and then we add random values to all of them. One stage of this is shown below.
The red points represent the initial grid, and the black points represent the new points, whose heights are calculated by interpolation. All of these points then have random values added to them.
There exist a few different ways of calculating the heights of the new points. We can either just average the heights of the old points around them (linear interpolation) or we can try and follow the curve of the landscape a little more closely, using cubic interpolation. The easiest way to do this is to calculate the partial derivatives at the old points, and use these in conjunction with the heights to construct a bicubic surface between 4 adjacent grid points. This has the advantage of continuity of gradient at the edges of the surface, meaning that there are no ugly horizontal and vertical ridges.
Controlling Fractal Dimension of the New Surface
As was mentioned previously, it is necessary to be able to control the fractal parameters of the surface. The most important of these is fractal dimension, which is a measure of how rough the landscape is. The fractal dimension of a fractal landscape must lie between 2.0 and 3.0. Outside these ranges, we cease to have fractal behaviour. For no reason other than simplicity, it is also decided that the mean fractal dimension of the landscapes we generate should be constant. It can be shown that a necessary and sufficient condition for constant mean fractal dimension is a frequency profile of:
Amplitude = 1 / frequencyb
Where beta is a constant for any given landscape. We can approximate this frequency profile by assuming that the frequency at any iteration of our subdivision scheme is given by the reciprocal of the smallest square size at that point in time.
There remains one thing to be done: we need to work out what range of values for beta will give rise to fractal dimensions between 2.0 and 3.0. As it turns out, they are related by the following formula:
Fractal Dimension = (7 - b)/2
This means that beta must be between 1.0 and 3.0. At b = 1.0, the system will only just converge (fractal dimension = 3.0), while at b = 3.0, the landscape would be very smooth indeed (fractal dimension = 2.0).
Non-Isotropic Behaviour and some Solutions
The simple subdivision algorithm already described, while quite powerful in some circumstances, does not behave quite ideally. The most noticeable artefact is that landscape features can appear somewhat square. This arises because the algorithm operates on a square grid, and the interpolation is not sufficient to remove all directional dependence. In other words, the landscape is not isotropic.
The Diamond-square algorithm is one frequently used solution to this problem. This works by subdividing the square grid into diamond-shaped sections, and the subdividing these into square-shaped sections again.
Thus, at each step, the interpolated position of one of the new points depends on exactly four other points.
Unfortunately, while it seems like a good idea on paper, if you use this in conjunction with linear interpolation, you get some numerical instabilities in the form of spires in the landscape. Using cubic interpolation, this effect is at least partially muted, and can produce landscapes with some desirable characteristics.
Getting Rid of Stationarity
As mentioned earlier, one of the desirable characteristics of a basic single fractal surface is stationarity - having the same statistical properties regardless of location. The algorithms mentioned all approximate stationarity, to a greater or lesser degree. Now, suppose you take out a map, and have a good long look at it. You will see that mountains are bumpy and valleys are flatter. You will see that the higher you get, the more ridges and protrusions appear, and the quicker the gradient fluctuates. All of these properties are non-stationary attributes, and we must find some way of mimicking them in our fractal generator, if we hope to generate interesting and varied terrain on the same map.
There are two principle ways in which this can be achieved. The first is to take advantage of the properties of multiplication, and multiply two fractal maps together, assuming that both maps are scaled in the range 0 to 1. This creates a new map, on which the roughness is proportional to the height. Hence, we have made the valleys smoother and the mountains rougher. Crucially, we can now multiply maps of different Unit Square Sizes together, and of different betas and so on, each time generating a completely different effect at the desired scale. The whole point of fractals is that they operate on many different scales at once, making this a particularly desirable attribute.
While this can be quite a subtle tool in some cases, it is still something of a hack - every effect on each scale must be added separately by the user. In the limit, this almost amounts to making your own fractal by hand - which is EXACTLY what we want to avoid doing.
By far the most portable and advanced solution is the fractal modification of crossover scale. The crossover scale of a landscape, as I define it here, is the scale at which fractal behaviour starts. If the lowest frequency element in the map had a wavelength of 256 pixel span, then the crossover scale would be 256 everywhere on the map. Suppose, however, that it were possible to make the crossover scale continuously variable. We would have some areas (those with large crossover scales) that were covered by ranges of high mountains. We would have other areas (those with small crossover scales) that were low, smooth and almost flat down to quite small scales, but without losing the fractal attributes that are characteristic of natural phenomena.
What FLG does
FLG provides a simple user interface which allows the user to specify crossover scales, multiply fractal maps together, use linear or cubic interpolation and change the distribution of the random numbers that are used to generate the landscape. It also provides a number of useful post-processing operations - such as canyonization and glaciation, and an advanced feature which enables it to produce fold mountains with fractal folding functions.
Catmull-Rom Splines, Cubic surfaces and Interpolating Subdivision
Because of the importance of splines in generating fractal landscapes, I thought I might present some additional notes on how to interpolate grid points using efficient subdivision schemes. For those of you not interested in the methods, here is a summary of the results:
Catmull Rom Cubic Coefficients
The cubic catmull-rom spline is given by:
Where y is the result vector, t is the parameter (between 0 and 1) and x(n) are the grid vectors. In other words, we have 4 points, x1...x4 and we're trying to interpolate smoothly between the middle two, x2 and x3. At t=0, the above equation evaluates to y = x2, and at t=1, it evaluates to x3. Intermediate values give points in-between.
1-Dimensional Cubic subdivision
Evaluating this at t=0.5 gives us a cubic approximation for the midpoint of x2 and x3, which is exactly what we need for subdivision:
This is the cubic approximation for the midpoint of x2 and x3.
2-Dimensional Cubic subdivision
To obtain a 2-D subdivision scheme, we just multiply the coefficient vector of the 1-dimensional scheme by a 4x4 point matrix, and then by the 1-d vector again:
This gives us an expression for the midpoint of the interior 4 points, in terms of the surrounding 16 points.
Bi-cubic resampling
Therefore, to up-sample the entire landscape by 2 (prior to fractal perturbation), all you have to do is apply (1) to the top edge and the left edge, and apply (2) to get the midpoint of each patch. If you have a repeating landscape (which is the easiest choice), then the values simply wrap around at the edges.
How do derive it all
I'm not going to do a rigorous mathematical derivation of all this, because I don't really want to write that much. However, the basic sketch goes as follows:
- Write out the taylor series expansions of f(a+h) and f(a-h)
- Re-arrange each of them in terms of f'(a)
- Combine the two equations to cancel the error terms. This gives you a 3-point derivative operator:
- Write down the generalized cubic:
- Observe that for gradient continuity at the endpoints of the x2...x3 interval, the following conditions must be true:
- At t=0, f(t) = x2
- At t=1, f(t) = x3
- At t=0, f'(t) = h2
- At t=0, f'(t) = h3
Where h2 and h3 are the gradients at x2 and x3.
- Observe that the formula from 3) can be used to replace h2 and h3 with combinations of the points x1...x4.
- Use the conditions to form a set of 4 simultaneous equations, and solve for the coefficients A, B, C and D
The end result of this process is the equation for the cubic Catmull-Rom spline that we started with.
|
