RANDOM FRACTAL FORGERIES

Richard F. Voss
IDM Thomas J. Watson Research Center
Yorktown Heights, NY 10598 USA

ABSTRACT

Mandelbrot's fractal geometry provides both a description and a mathematical model for
many of the seemingly complex shapes found in nature. Such shapes often possess a re-
markable invariance under changes of magnification. This statistical self-similarity may
be characterized by a fractal dimension, a number that agrees with our intuitive notion
of dimension but need not be an integer. In section I, a series of computer generated and
rendered random fractal shapes provide a visual introduction to the concepts of fractal
geometry. These complex images, with details on all scales, are the result of the simplest
rules of fractal geometry. Their success as forgeries of the natural world has played an
important role in the rapid establishment of fractal geometry as a new scientific discipline
and exciting graphic technique. Section IT presents a brief mathematical characterization
of these forgeries, as variations on Mandelbrot's fractional Brownian motion. The im-
portant concepts of fractal dimension and exact and statistical self-similarity and self-
affinity will be reviewed. Finally, section m will discuss independent cuts, Fourier
filtering, midpoint displacement, successive random additions, and the Weierstrass-
Mandelbrot random function as specific generating algorithms.

I. VISUAL INTRODUCTION TO FRACTAL FORGERIES

Figure C 1, like all of the color illustrations accompanying this article (but collected else-
where in this volume), is an example of a fractal forgery. It is actually composed of three
separate forgeries which will be discussed in more detail below: the fractal landscape in
the foreground, Fig. C2, its fractally distributed craters, Fig. C3, and the generalization
of Brownian motion onto a sphere rising in the background, Fig. C4. The objects re-
presented in the images exist only as arrays of pseudo-random numbers in a computer.
Yet, the resemblance to the real world is unmistakable. My success with these forgeries
and the extent to which they capture some of the essential characteristics of nature is due
to a new branch of mathematics, fractal geometry [1], as conceived and developed by
Benoit Mandelbrot.

NATO ASI Series. Vol. F17

Fundamental Algorithms for Computer Graphics
Edited by R. A. Earnshaw
© Springer-Verlag Berlin Heidelberg 1985
 806

According to Galileo (1623),

"Philosophy is written in this grand book - I mean universe - which stands continuously
open to our gaze, but it cannot be understood unless one first learns to comprehend the
language in which it is written. It is written in the language of mathematics, and its
characters are triangles, circles and other geometrical figures, without which it·is hu-
manly impossible to understand a single word of it; without these, one is wandering
about in a dark labyrinth" .

Galileo was, of course, stating one of the basic tenets of modern science. In order to un-
derstand or simulate nature, one must first comprehend its language. To scientists, na-
ture's languages are mathematics and geometry presents the specific structures for
understanding and manipulating shapes. Galileo was, however, wrong in terms of the
dialect used by much of nature. The inability of the triangles and circles of Euclidean
geometry (the actual dialect to which he referred) to accurately describe the natural
world has been unappreciated until recently. In retrospect it is obvious that, "clouds are
not spheres, mountains are not cones, coastlines are not circles, and bark is not smooth,
nor does lightning travel in a straight line" [Mandelbrot, ref. 1]. Fortunately, there is now
a geometry that is appropriate for the irregular shapes of the real world.

mathematical monsters: some early fractals

The turn of the century coincided roughly with an upheaval in the world of mathematics.
Minds conceived of strange monsters seemingly without counterpart in nature: Cantor
sets, Weierstrass functions, Peano curves, curves without derivatives and lines that could
fill space. Having once discovered these monsters (and congratulated themselves on a
creativity superior to nature), mathematicians banished the beasts, mostly sight unseen,
to a mathematical zoo. They could imagine no use for, nor interest in, their creations by
natural scientists. Nature, however, was not so easily outdone. As shown by Benoit
Mandelbrot, a mathematician at IBM and Harvard, many of these "monsters" do, in fact,
have counterparts in the real world. illuminated by computer graphics, they have finally
been recognized as some of the basic structures in the language of nature's irregular
shapes: the fractal geometry of nature.

Fractals (a word coined by Mandelbrot in 1975) have blossomed tremendously in the

past few years and have helped reconnect pure mathematics research with both the na-
tural sciences and computing. Within the last 5-10 years fractal geometry and its con-
cepts have become central issues in most of the natural sciences: physics, chemistry,
biology, geology, meteorology, and materials science. At the same time, fractals are of
 807

interest to graphic designers and filmmakers for their ability to create new and exciting
shapes and artificial but realistic worlds. Computer graphics has played an important role
in the development and rapid acceptance of fractal geometry as a valid new discipline.
The computer rendering of fractal shapes leaves no doubt of their relevance to nature.
Conversely, fractal geometry now plays a central role in the realistic rendering and
modelling of natural phenomena in computer graphics.

fractals: self-similarity and dimension

What then are fractals? How are they different from the usual Euclidean shapes? The
answer can be illustrated with one of the early mathematical monsters: the von Koch
snowflake curve (first proposed around 1904). Figure 5 illustrates an iterative or recur-
sive procedure for constructing a fractal curve. A simple line segment is divided into
thirds and the middle third is replaced by two segments forming part of an equilateral
triangle. At the next stage in the construction each of these 4 segments is replaced by 4
new segments with length 1/3 of their parent according to the original pattern. This
procedure, repeated over and over, yields the beautiful von Koch curve shown at the top
of Fig. 5 and demonstrates that the iteration of a very simple rule can produce seemingly
complex shapes with some highly unusual properties. At each stage in its construction
the length of the curve increases by a factor of 4/3. Thus, the limiting curve crams an
infinite length into a finite area of the plane without intersecting itself. The curve has
detail on all length scales. Indeed, the closer one looks, the more powerful the microscope
one uses, the more detail one finds. More important, the curve possesses an exact self-
similarity. Each small portion, when magnified, can reproduce exactly a larger portion.
The curve is said to be invariant under changes of scale.

This property of self-similarity or scaling is one of the central concepts of fractal geom-
etry. It is closely connected with our intuitive notion of dimension. An object normally
considered as one-dimensional, a line segment, for example, also possesses a similar
scaling property. It can be divided into N identical parts each of which is scaled down
by the ratio r= liN from the whole. Similarly, a two-dimensional object, such as a square
in the plane, can be divided into N self-similar parts each of which is scaled down by a
factor r= 1IN I 12. A three-dimensional object like a solid cube may be divided into N
little cubes each of which is scaled down by a ratio r=1/N I /3. With self-similarity the
generalization to fractal dimension is straightforward. A D-dimensional self-similar ob-
ject can be divided into N smaller copies of itself each of which is scaled down by a factor
r where r= 1/N liD or
 808

(1.1)

Conversely, given a self-similar object of N parts scaled by a ratio r from the whole, its
fractal or similarity dimension is given by

D = 10g(N)/ 10g(1/r). (1.2)

The fractal dimension, unlike the more familiar notion of Euclidean dimension, need not
be an integer. Any segment of the von Koch curve is composed of 4 sub-segments each
of which is scaled down by a factor of 1/3 from its parent. Its fractal dimension D =
log(4)/log(3) or about 1.26.... This non-integer dimension, greater than one but less
than two, reflects the unusual properties of the curve. It somehow fills more of space than
a simple line (D= 1), but less than a Euclidean area of the plane (D=2). Mandelbrot[l]
gives many variations of the von Koch construction. As D increases from 1 toward 2 the
resulting "curves" progress from being "line-like" to "filling" much of the plane. Indeed,
the limit D-2 gives a Peano or "space-filling" curve. Although the von Koch curve has
a fractal dimension D = log(4)/log(3), it remains a "curve" with a topological dimension
of one. The removal of a single point cuts the curve in two pieces.

statistical self-simllarity

Although objects in nature rarely exhibit exact self-similarity, they do often possess a
related property, statistical self-similarity. Perhaps the simplest example is a coastline.
As with the von Koch curve, the closer one looks at a coastline, the more detail one sees.
In a measurement of the length of a coastline, the more carefully one follows the smaller
wiggles, the longer it becomes. A drive along a coast highway is shorter than a walk along
the beach. Moreover, each small section of a coastline looks like (but not exactly like)
a larger portion. This concept of statistical self-similarity can also be quantified in terms
of fractal dimension as above. When using a ruler of size r to measure a coastline's
length, the total length equals the ruler size r x the number of steps of size r, N(r), taken
in tracing the coast
LENGTH = r x N(r).
As with the snowflake, N(r) varies on the average as 1/rD and

LENGTH ex: r x l/rD = 1/rD-t. (1.3)

The variation of apparent coastline length with ruler size has been studied by Richardson
as summarized in [1]. Real coastlines can, in fact, be characterized by fractal dimensions
D of about 1.15 to 1.25, close to the log(4)/log(3) of the von Koch curve.
 809

The property that objects can look statistically similar while at the same time different in
detail at different length scales, is the central feature of all of the random fractal forgeries
presented here. Figure 6 shows a series of views, at increasing magnification, of a com-
puter generated fractal landscape with just such a statistically self-similar coastline. The
upper left image shows a section of the landscape. A small portion, outlined by a black
box, is then magnified to fill the second frame. In each succeeding frame a portion about
the coastline, indicated by a small box, is magnified to fill the next frame. The total
magnification is given above each frame. This zooming process can be continued indef-
initely. In Fig 6. the magnification is carried to a total of more than 16 million and the
final frame repeats the starting view. Each frame is clearly a magnification of a portion
of the previous, yet each frame also looks like a different part of the same landscape at
the same magnification. This is the statistical self-similarity of a random fractal. It is
random in the sense that (unlike the von Koch curve) a large scale view is insufficient to
predict the exact details of a magnified view. The way in which the detail varies as one
changes length scale is once again characterized by a fractal dimension. The irregular
fractal surface is more than a simple surface (D=2) and the sample shown in Fig. 6 has
D=2.2. The fractal dimension of the coastline is one less than that of the surface itself.
Here the coastline D= 1.2.

Mandelbrot landscapes

For the idealized Mandelbrot fractal landscape of Fig 6. the statistical self-similarity
extends from arbitrarily large to arbitrarily small scales. Actual landscapes, on the other
hand, can be statistically self-similar only over a finite (but often quite large) range of
distances. The largest variations may be limited by the size of the planet or the force of
gravity (the materials may not be strong enough to support arbitrarily high mountains).
The smallest scales may be limited by the smoothing of erosion, the basic grain size of the
rock and sand or, at the very least, by the atomic nature of the particles. The mathemat-
ical ideal is an approximation to the real world (or vice versa to some mathematicians).
Mandelbrot's fractal geometry, however, remains by far the best approximation and, to
date, the most widely successful mathematical model.

As discussed below, almost all algorithms for generating the fractal landscapes effectively
add random irregularities to the surface at smaller and smaller scales similar to the proc-
ess of adding smaller and smaller line segments to the von Koch curve. Once again, the
fractal dimension determines the relative amounts of detail or irregularities at different
distance scales. Surfaces with a larger D seem rougher. This effect is illustrated in Fig.
C7 which shows the "same" surface but with different fractal dimensions. Figure C7(a)
 810

has a relatively low fractal dimension for the surface, D = 2.15, that is appropriate for
much of the earth. An increase in the fractal dimension to D=2.5 is shown in Fig C7(b).
A further increase to D=2.8 in Fig. C7(c) gives an unrealistically rough surface for an
actual landscape. Fractal geometry specifies only the relative height variations of the
landscape at different length scales. For the samples given here the color changes were
based on height above water level and local slope.

The image in Fig. C7(a) does not, however, resemble many of the valleys on the Earth's
surface which are significantly eroded. A flat bottomed basin can, however, be approxi-
mated mathematically by simply taking the height variations and scaling them by a
power-law. This process is illustrated in Fig. C7(d) where the original landscape of Fig.
C7(a) is "cubed". The effect of such a power greater than 1 is to flatten the lower ele-
vations near near the water emphasizing the peaks. Scaling with a power less than one
has the opposite effect of flattening the peaks while increasing the steepness near the
water. A cube-root processing of Fig. C7(a) is shown in Fig. C7(e). This forgery gives
the impression of river erosion into an otherwise fairly smooth plain.

fractally distnDuted craters

In order to give the impression of a lunar landscape, as shown in the foreground of Fig.
1, it is necessary to add craters to the "virgin" fractal landscape of Fig C7(a). Each
crater is circular with a height profile similar to the effect of dropping marbles in mud.
The trick in achieving realism is to use the proper distribution of crater sizes. For the
moon, the actual distribution is fractal or power-law, with many more small craters than
large ones. Specifically, the number of craters having an area, a, greater than some
number A, N(a>A) varies as 1/A. The computer simulation is shown in Fig. C3. This
is the same surface that forms the foreground of Fig. Cl.

fractal planet: Brownian motion on a sphere

Rather than add increasing detail at progressively smaller scales, the rising fractal planet
in Fig. CI was generated by a different algorithm that is closer to an actual model for the
evolution of the earth's surface. Its surface is a generalization of Brownian motion or a
random walk on a sphere. A random walk is the sum of many independent steps, and the
sphere's surface becomes the result of many independent surface displacements or faults.
Each fault encircles the sphere in a random direction and divides it into two hemispheres
which are displaced relative to each other in height. Inside the computer, the surface of
the sphere is mapped onto a rectangular array similar to a flat projection map of the earth.
Figure C8(a) shows the surface variations, as represented by regions of different color,
 811

both for the full and the flattened sphere after only 10 random faults. After 60 faults in
Fig. C8(b), the flattened sphere gives the impression of a cubist painting. With 750
faults in Fig. C8(c), the sphere has developed a localized land mass and the effects of an
individual fault are becoming lost. After more than 10,000 faults and the addition of
polar caps in Fig C8(d), the sphere's surface becomes a plausible planet with D=2.5.
Mapping the data of Fig. C8(d) back onto a sphere gives the fractal planet in Fig. 1 or
the enlarged view of the rotated sphere shown in Fig. C8(e). Remarkably, the same
sphere is also a good forgery of a much more arid generic planet. Fig. C8(e) shows the
same data but with different coloring and without water.

fractal flakes and clouds

Up to now the fractal forgeries have been pretentious surfaces. The addition of irreg-
ularities on smaller and smaller scales raise their dimensions above a topological value
of 2 to a fractal value 2<D<3. It is, however, also possible to generate fractals with a
topological dimension 3 whose scaling irregularities raise the fractal dimension to
3<D<4. Rather than landscapes, which consist of random heights as a function of
2-dimensional surface position, such objects might represent a random temperature or
water vapor density as a function of 3-dimensional position. The simplest means of
viewing such constructions is with their zerosets (equivalent to the coastlines of the fractal
landscapes) as shown in Fig. C9. For a temperature distribution T(x,y,z) with fractal di-
mension 3<D<4, the zerosets are all points for which T(x,y,z)-To = 0 and have a fractal
dimension of D-l. They can be displayed by treating all points for which T(x,y,z»To
as opaque while all others are transparent. The resulting fractal flakes are shown in Fig.
C9(a) where D-l=2.5 and'Fig. C9(b) where D-l=2.2. A more accurate display is also
possible by allowing the local light scattering to vary as T(x,y,z). In this case, the flakes
are transformed into realistic looking fractal clouds as shown in Figs C9(c) and (d). The
correspondence to actual clouds is not surprising. Clouds and rain areas are the natural
shapes where fractal behavior has been experimentally verified over the widest range of
distance scales[lJ. A combination of such a fractal cloud above a cratered fractal surface
is shown in Fig. C9(e). Here light scattered by the cloud produces shadows of varying
intensity on the landscape.

scaling randomness in time: l/ff3 noises

All of the above forgeries have been based on fractal randomness in space. The com-
puter generated images show, the extent to which simple fractal constructions capture the
essence of natural shapes. Changes in time, however, have many of the same similarities
 812

at different scales as changes in space. To the physicist, unpredictable changes of any

quantity V varying in time t are known as noise. Graphical samples of typical noises V(t)
are shown in Fig. 10. To the left of each sample is a representation of its spectral
density. The spectral density, Sy(f), gives an estimate of the mean square fluctuations
at frequency f and, consequently, of the variations over a time scale of order 1/f. The
traces made by each of these noises is a fractal curve and there is a direct relationship
between the fractal dimension and the logarithmic slope of the spectral density which
will be discussed below. Indeed, it is the understanding and simulation of such noises that
forms the basis for all of the fractal forgeries presented here.

Figure 10(a) shows a white noise, the most random. It could be produced by a pseudo-
random number generator and is completely uncorrelated from point to point. Its spectral
density is a flat line, independent of frequency. Alternately, a white noise contains equal
amounts of all frequencies (like a white light). Figure 10(c) shows a Brownian motion
or a random walk, the most correlated of the three noise samples. It consists of many
more slow (low frequency) than fast (high frequency) fluctuations and its spectral den-
sity is quite steep. It varies as l/f2. Formally, the Brownian motion of Fig. 10(c) is the
integral of the white noise of Fig. lO(a). In the middle, in Fig. 10(b) is an intermediate
type of noise known as 1If noise because of the functional form of its spectral density.
In general, the term llf noise is often applied to any fluctuating quantity V(t) with
Sy(f) varying as 1/fP over many decades with 0.5</3<1.5. Although both white and
1/f2 noise are well understood in terms of mathematics and physics, 1/f noise, whose
origin remains a mystery after more than 60 years of investigation, represents the most
common type of noise found in nature.

There are no simple mathematical models that produce llf noise other than the
tautological assumption of a specific distribution of time constants. Little is also known
about the physical origins of llf, but it is found in many[2] physical systems: in almost
all electronic components from simple carbon resistors to vacuum tubes and all semi-
conducting devices; in all time standards from the most accurate atomic clocks and quartz
oscillators to the ancient hourglass; in ocean flows and the changes in yearly flood levels
of the river Nile[3] as recorded by the ancient Egyptians; in the small voltages measur-
able across nerve membranes due to sodium and potassium flow; and even in the flow
of automobiles on an expressway[4]. 1/f noise is also found in music.
 813

fractal music

One of my most exciting discoveries[5,6] was that almost all musical melodies also mimic
l/f noise. Music has the same blend of randomness and predictability that is found in
1/f noise. If one takes a music score and draws lines between successive notes of the
melody one finds a graph remarkably similar in quality to Fig. 10(b). Some of the actual
measured spectral densities for different types of music is shown in Fig. 11. There is little
to distinguish these measurements on widely different types of music from each other or
from the 1/f noise of Fig. 10(b). This type of analysis is surprisingly insensitive to the
different types of music. It emphasizes, instead, the common element in music and sug-
gests an answer to a question that has long troubled philosophers[7]. In the words of
Plato, "For when there are no words (accompanying music), it is very difficult to recog-
nize the meaning of the harmony and rhythm, or to see that any worthy object is imitated
by them". Philosophers generally agreed on the imitative nature of the arts. It seemed
obvious that painting, sculpture or drama imitated nature. But what does music imitate?
The measurements suggest that music is imitating the characteristic way our world
changes in time. Both music and l/f noise are intermediate between randomness and
predictability. Like fractal shapes there is something interesting on all (in this case, time)
scales. Even the smallest phrase reflects the whole.

It is, of course, possible to use fractals, in this case as 1/fP noises, for music as well as
landscape forgery. Fig. 12 shows samples of "music" generated from the three charac-
teristic types of "noise" shown in Fig. 10. Although none of the samples in Fig. 12 cor-
respond to a sophisticated composition of a specific type of music, Fig. 12(b) generated
from 1/f noise is the closest to real music. Such samples sound recognizably musical, but
from a foreign or unknown culture.

fractals and the return to nature

Man has always been confronted with a world filled with seemingly complex irregular
shapes and random fluctuations. In the quest for understanding, natural science has
progressed by concentrating primarily on the simplest of systems. In this process, it has
moved away from the direct experience of nature to the electronically instrumented lab-
oratory. After all, who could describe, let alone understand, the profile of a mountain or
the shape of a cloud? With fractal geometry, the quest for scientific understanding and
realistic computer graphic imagery can return to the everyday natural world. As these
forgeries show, nature's complex shapes are only seemingly so.
 814

D. MATHEMATICAL CONSIDERATIONS: fractional Brownian motion

This section presents an expository summary of the major mathematical definitions and
relations used in the random fractal forgeries as condensed from Mandelbrot[I]. When
no other reference is indicated [1] may be assumed as the source.

The most useful mathematical model for the random fractals found in nature has been the
fractional Brownian motion (mm) of Mandelbrot and Wa1lis[1,8]. It is an extension of the
central concept of Brownian motion that played an important role in both physics and
mathematics. Moreover, it forms the basis for all of the forgeries presented here and the
standard by which various generating algorithms may be compared.

A fractional Brownian motion, VH(t), is a single valued function of one variable, t (usu-
ally time). Its increments VH(t2)-VH(tt) have a Gaussian distribution with variance

(0.1)

where the brackets < and> denote averages over many samples of VH(t) and the pa-
rameter H has a value O<H<I. Such a function is both stationary and isotropic. Its mean
square increments depend only on the time difference trtt and all t's are statistically
equivalent. The special value H=I/2 gives the familiar Brownian motion with tJ.V2 a:
tJ.t.

As with the usual Brownian motion, although VH(t) is continuous, it is nowhere

differentiable. Nevertheless, many constructs have been developed (and are relevant to
the problem of light scattering from fractals) to give meaning to "derivative of fractional
Brownian motion" as fractional Gaussian noises. Such constructs are usually based on
averages of VH(t) over decreasing scales. The derivative of normal Brownian motion,
H= 1/2, corresponds to the uncorrelated white Gaussian noise of Fig. 10(a), and
Brownian motion is said to have independent increments. Formally, for any three times
such that tt<t<t2' tJ.Vt=VH(t)-VH(tt) is statistically independent of
tJ.V2=VH(t2)-VH(t) for H=1/2. For H>I/2 there is a positive correlation both for the
increments of VH(t) and its derivative fractional Gaussian noise in Figs. 10 (b) and (c).
For H<I/2 the increments are negatively correlated. Such correlations, moreover, ex-
tend to arbitrarily long time scales.

As with coastlines, VH(t) shows a statistical scaling behavior. If the time scale t is
changed by the factor r, then the increments tJ.VH change by a factor~. Formally,

(11.2)
 815

Unlike statistically self-similar coastlines, however, a VH(t) trace requires different scal-
ing factors in the two coordinates (r for t but ~ for VH) reflecting the special status of
the t coordinate. Each t can correspond to only one value of VH but any specific VH
may occur at multiple t's. Such non-uniform scaling is known as self-affinity rather than
self-similarity.

self-simiIar vs self-aff"me fractals.

The distinction between similarity and affinity is important. By way of summary[l], a

self-similar object is composed of N copies of itself (with possible translations and ro-
tations) each of which is scaled down by the ratio r in all E coordinates from the whole.
More formally, consider a set S of points at positions x = (xl, ... , xE) in Euclidean space
of dimension E. Under a similarity transform with real scaling ratio O<r<l, the set S
becomes rS with points at rx = (rxl> ... , rXE)' A bounded set S is self-similar when S is
the union of N distinct (non-overlapping) subsets each of which is congruent to rS.
Congruent means identical under translations and rotations. The fractal or similarity di-
mension of S is then given by

logN
or D = - - - (11.3)
log l/r

It is useful to note that the fractal dimension D also characterizes the covering of the set
S by E-dimensional "boxes" of linear size L. If the entire S is contained within one box
of size L max , then each of the N=l/~ subsets will fall within one box of size
L = rLmax' Thus, the number of boxes of size L, N(L), needed to cover S is given by

(11.4)

This definition of box dimension is one of the most useful methods for estimating the
fractal dimension of a given set.

One can also estimate the "volume" or "mass" of the set S by a covering with boxes of
linear size L. If one considers only distances of order L about a given point in S, one
finds a single box of size L with E-dimensional volume L E. If the distance scale about the
same point is increased to Lmax=L/r, one now finds a total of N=l/~=(Lmax/L)D
boxes of the same size boxes. Thus, the mass within a distance Lmax of some point in S,
M(Lmax) = NxM(L) = M(L)x (Lmax/L)D or

M(L) IX: LD. (11.5)

 816

The fractal dimension D, thus, also corresponds to the commonly used mass dimension
in physics. Mass dimension also has a strong connection with intuitive notions of di-
mension. The amount of material within a distance L of a point in a one-dimensional
object increases as V. For an E-dimensional object it varies as LE. The mass dimension
is another extremely useful method for estimating the fractal dimension of a given object.

The set S is also self-similar if each of the N subsets is scaled down from the whole by a
different similarity ratio rn. In this case, D is given implicitly by

N
1= L
0=1
(11.6)

which reduces to the familiar result in Eq. (11.3) when all of the rn are equal.

The set S is statistically self-similar if it is composed of N distinct subsets each of which

is scaled down by the ratio r from the original and is identical in all statistical respects to
rS. The similarity dimension is again given by Eq. (11.3). In practice, it is impossible to
verify that all moments of the distributions are identical, and claims of statistical self-
similarity are usually based on only a few moments. Moreover, a sample of a random set
(such as a coastline) is often statistically self-similar for all scaling ratios t. Its fractal di-
mension is usually estimated from the dependence of box coverings N(L) or mass M(L)
on varying L as in Eqs. (11.4) and (11.5).

Under an affine transform, on the other hand, each of the E coordinates of x may be
scaled by a different ratio (r!> ... , rE). Thus, the set S is transformed to r(S) with points
at r(x) = (rix!> ... , rExE). A bounded set S is self-affine when S is the union of N dis-
tinct (non-overlapping) subsets each of which is congruent to r(S). Similarly, S is statis-
tically self-affine when S is the union of N distinct subsets each of which is congruent in
distribution to r(S). The fractal dimension D, however, is not as easily defined as with
self-similarity.

the relation of D to H for self-aff"me fractional Brownian motion

The assignment of a fractal dimension D to a self-affine set can be illustrated with a trace
of fractional Brownian motion VH(t) from above. It may be helpful to glance again at
Fig. 10. Consider, for convenience, a trace of VH(t) covering a time span Ilt= 1 and a
vertical range 11VH= 1. VH(t) is statistically self-affine when t is scaled by rand VH is
scaled by r. H. Suppose the time span is divided into N equal intervals each with
Ilt=l/N. Each of these intervals will contain one portion of VH(t) with vertical range
 817

Il VH = IltH = l/NH. Since O<H< 1 each of these new sections will have a large vertical
to horizontal size ratio and the occupied portion of each interval will be covered by
IlVH/llt = (l/NH)/(l/N) = N/NH square boxes of linear scale L=1/N. In terms of
box dimension, as t is scaled down by a ratio r= l/N the number of square boxes covering
the trace goes from 1 to N(L) = number of intervals x boxes per interval or

(11.7)

Thus, by comparison with Eq. (11.4),

D= 2- H for a l/fP noise. (11.8)

Consequently, the trace of normal Brownian motion has D = 1.5.

The zeroset of mm is the the intersection of the trace of VH(t) with the t axis, the set of
all points such that VH(t) = O. The zeroset is a disconnected set of points with
topological dimension zero and a fractal dimension Do = D-1 = 1-H that is less than 1
but greater than O. Although the trace of VH(t) is self-affine, its zeroset is self-similar.

self-afrmity in higher E: Mandelbrot landscapes and clouds

The traces of mm such as Fig. 10(c) bear a striking resemblance to a mountainous hori-
zon. The modelling of the irregular Earth's surface as a generalization of traces of mm
was first proposed by Mandelbrot. The single variable t can be replaced by coordinates
x and y in the plane to give VH(x,y) as the surface altitude at position x,y. In this case,
the altitude variations of a hiker following any straight line path at constant speed in the
xy plane is a fractional Brownian motion. In analogy with Eq. (11.1),

(11.9)

Once again, the fractal dimension D must be greater than the topological dimension 2 of
the surface. Here,

D= 3- H for a fractal landscape. (11.10)

The intersection of a vertical plane with the surface VH(x,y) is a self-affine mm trace
with D=2-H, smaller by one than the value of Eq. (11.10). Similarly, the zeroset of
VH(x,y), its intersection with a horizontal plane, also has a fractal dimension Do=2-H.
This intersection, which produces a family of (possibly disconnected) curves, forms the
coastlines of the VH(x,y) landscape. Since the two coordinates x and yare, however,
equivalent, the coastlines of VH(x,y) are self-similar, not self-affine.
 818

This generalization of rum can continue to still higher dimensions to produce the self-
affine fractal temperature distribution T(x,y,z) as VH(x,y,z). Here, the temperature var-
iations of an observer moving at constant speed along any straight line path in space
generate a rum and the fractal dimension

D= 4- H for a fractal cloud. (11.11)

The zeroset VH(x,y,z) = constant now gives a self-similar fractal with Do=3-H.

To summarize, a statistically self-affine fractional Brownian function, V H of x = (xl>

... , xE) in E+ 1 Euclidean dimensions satisfies

(11.12)

and has a fractal dimension

D =E + 1 - H. (11.13)

The zerosets of VH(x) form a statistically self-similar fractal with dimension Do = E-H.

spectral densities for mm and the spectral exponent f3

As mentioned in section I, random functions in time V(t) are often characterized[9,10]

by their spectral densities Sy(f). If V(t) is the input to a narrow bandpass filter at fre-
quency f and bandwidth M, then Sy(f) is the mean square output V(f) divided by M,
Sy(f) = 1V(f) 12/ M. Sy(f) gives information about the time correlations of V(t). As
shown in Fig. 10, when Sy(f) increases steeply at low f, V(t) varies more slowly. If one
defines V(f,T) as the Fourier transform of a specific sample of V(t) for O<t<T,

V(f,T) = ~foTV(t)e2'ITiftdt, then Sv(f)a:T 1V(f,T) 12 as T .. 00. (11.14)

An alternate characterization of the time correlations of V(t) is given by the 2 point

autocorrelation function

Gy ( T) provides a measure of how the fluctuations at two times separated by T are related.
Gy ( T) and Sy(f) are not independent. In many cases they are related by the Wiener-
Khintchine relation[ 9,10]
 819

(11.15)

For a Gaussian white noise Sy(f) = constant and Gy('I") = ~V28('I") is completely un-
correlated. For certain simple power laws for Sy(f), Gy ( '1") can be calculated exactly.
Thus, for

(11.16)

Moreover, Gy('I") is directly related to the mean square increments of mm,

(11.17)

Roughly speaking, Sy(f) ex: l/ffJ corresponds to Gy('I") ex: 'l"l-fJ and a mm with 2H=~-1
from Eqs. (11.1) and (11.17). Thus, the statistically self-affine fractional Brownian func-
tion, VH(i), has a fractal dimension D and spectral density Sy(f) ex: l/ffJ, for the fluctu-
ations along a straight line path in any direction in E-space with

3-~
D = E+I-H = E + - - (II. IS)
2

This result agrees with other[I,S,II] "extensions" of the concepts of spectral density and
Wiener-Khintchine relation to non-stationary noises where some moments may be unde-
fined. Moreover, it provides an extremely useful connection between D, H and ~ for fi-
nite simulations. For H in the range O<H<I, E<D<E+l, and 1<~<3. The value H ~

O.S is a good choice for many natural phenomena.

Although the formal definition of fBm restricts H to the range O<H< 1, it is often useful
to consider integration and an appropriate definition of "derivative" as extending the
range of H. Thus, integration of a mm produces a new mm with H increased by 1, while
"differentiation" reduces H by 1. When H .... 1, the derivative of mm looks like a mm
with H .... O. In terms of spectral density, if V(t) has Sy(f) ex: l/f fJ then its derivative
dV /dt has spectral density f2/ffJ = l/ffJ- 2 . In terms of Eq. (11.19), differentiation de-
creases ~ by 2 and decreases H by 1.
 820

III. ALGORITHMS: approximating mm on a ("mite grid

Section I established the visual connection between many, seemingly complex, shapes in
the natural world and statistically self-similar and self-affine fractals. Section II reviewed
some of the mathematical considerations and established the relation between fractal
dimension D, the H parameter of mm, and the spectral density exponent /3. This section
discusses various methods for producing a finite sample of mm as a noise (E= 1), a
landscape (E=2), or a cloud (E=3). As with nature itself, this sample will typically be
limited by both a smallest size or resolution A and a largest scale Lmax. Consideration
will be given to sample statistics (mathematically, how well does it approximate mm?),
visual features (does it look natural?) and computation characteristics (how does com-
putation time vary with sample size? can one magnify or extend the sample beyond its
boundaries?). In general, the closer a given algorithm approximates mm, the more "re-
alistic" the resulting image looks. Although no specific consideration will be given to
image rendering, it is important to note that in most cases rendering of a fractal sample
requires far more computation than sample generation. This is particularly true for ren-
dering packages based on Euclidean shapes which have an extremely difficult time with
the "infinite" number of surface patches on a fractal surface.

Brownian motion as independent cuts

In the usual definition, Brownian motion, VB(t) = VH(t) for H= 1/2, is the integral of a
Gaussian white noise W(t),

(ill.l)

Brownian motion may also be considered as the cumulative displacement of a series of

independent jumps. A pulse at time tj causes a jump of magnitude Aj (a Gaussian ran-
dom variable) in V(t). The response to such a pulse, AjP(t-1j) has the form of a step
function with P(t)=l for t>O and P(O)=l for t<O. Thus, VB(t) can also be written as
the sum of independent cuts at random (Poisson distributed) times 1j

VB(t) = L"" AjP(t - tj). (ill. 2)

i=-oo

This latter formulation of Brownian motion is useful since it can be generalized to circles
and spheres to produce the fractal planet of Fig. CS. The time coordinate t may be re-
placed with the angular position (J on a unit circle. A Brownian motion in (J must then
 821

have the periodicity that VB(O) = VB(O + 2'11). In this case, VB(O) = It" W(O')dO'.
The P(O) correspond to half circles, P(O)=l for 0<0<'11 and P(O)=O for '11<0<2'11, and
Brownian motion on a circle becomes the summation of independent "half-circles" uni-
formly distributed in OJ,

L
00

VB(O) = AjP(O - OJ (111.3)

i=-oo

In the generalization to a sphere, VB becomes a function of r ,the position on the unit

sphere, and corresponds to the addition of random hemispheres whose positions rj are
uniformly distributed on the sphere surface,

L
00

VB(r) = AjP(r - rJ (I1I.4 )

i=-oo

Here, Per - rj) = 1 for r • rj >0 and zero otherwise. The evolution of a Brownian
sphere under this summation process is shown in Fig. C8.

A flat Brownian relief VB(x,y) can similarly be constructed from the addition of ran-
domly placed and randomly oriented faults in the plane. The profile of such a fault cor-
responds to the step function of random amplitude Aj. Such a surface will have a fractal
dimension D=2.S.

The independent addition of such step function faults is extremely expensive computa-
tionally. Each fault requires additions to roughly half of the surface elements. Moreover,
the resulting fractal always corresponds to a fBm with H= 1/2. Nevertheless, the proce-
dure is straightforward and represents, historically, the first[l] method used for produc-
ing fractal landscapes.

random cuts with H #: 1/5: CampbeU's theorem

The summation in Eq. (III.2) is a special case of Campbell's theorem (1909). Consider
a collection of independent pulses occurring at random times tj corresponding to an av-
erage rate th. Each pulse produces the profile AjP(t-tj). The random function
Vet) = LAjP(t - t), will then have a spectral density[9],

(111.5)

where P(f) is the Fourier transform of pet), P(f) = IP(t)e 27Tiftdt. For normal Brownian
motion, each excitation by the white noise Wet) in Eq. (111.1) produces the step function
 822

response P(t) with P(f) ex: l/f and Sy(f) ex: 1/f2. With a suitable power-law choice of
P(t), one can generate an approximation to fBm with any H. In fact, fBm can be put in
a form[I,8] similar to Eq. (III.l),

f(H
1
+ 1/2)
It
-Q()
(t _ t,)H-l/2W (t')dt'.

Thus, fBm can be constructed from independent pulses with response P(t) = t H-1/ 2
corresponding to P(f) ex: 1/fH +1/2 from Eq. (11.16) and Sy(f) ex: l/f 2H + 1 in agreement
with Eq. (11.18). This was the first method used for generating samples of fractional
Gaussian noise and fractional Brownian landscapes. Unless only a few "pulses" or
"faults" are acceptable, this process is, in general, too expensive computationally.

Fast Fourier Transform fdtering

Another straightforward algorithm for fBm is to directly construct a random function

with the desired spectral density ex: l/ff3. For all practical purposes, the output of a
pseudo-random number generator produces a "white noise" W(t). Filtering W(t) with
a transfer function T(f) produces an output, V(t), whose spectral density,

Sv(f) ex: 1T(f) 12Sw(f) ex: 1T(f) 12.

Thus, to generate a l/f f3 noise from a W(t) requires T(f) ex: 1/ff3/ 2 .

A continuous function of time, V(t), may be approximated by a finite sequence of N

values, V n, defined at discrete times tn = nl1t, where n runs from 0 to N-l and I1t is the
time between successive values. The discrete Fourier transform (DFT) defines V n in
terms of the complex Fourier coefficients, vm' of the series:

(111.6)

where the frequencies fm = m/Nilt for m=O to N/2 - 1. For a fBm sequence with
Sy(f) varying as l/ff3 the coefficients must satisfy

(III. 7)

The relation between f3 and D and H is given by Eq. (11.18) with E= 1. The vm may be
obtained by multiplying the Fourier coefficients of a white noise sequence by l/f/3/2 or
by directly choosing complex random variables with mean square amplitude given by Eq.
 823

(ill.7) and random phases. One possible variation sets 1vm 1 = 1/mP/ 2 and only ran-
domizes the phases. With an FFT (Fast Fourier Transform) algorithm the evaluation of
Eq. (ill.6) requires of order NlogN operations[12] to produce a series of N points. As
such it offers a significant improvement over the random cuts described above.

Since the FFT includes all possible frequency components over the range of distance
scales from Ato Lmax=NA it represents an accurate approximation to ffim. It has, how-
ever, several drawbacks. The entire sample must be computed at once and it is difficult
to vary the degree of detail across the sample or to extend the sample across its original
boundaries. In addition, the result is periodic in time, Vn = Vn+N. Such boundary
constraints become more obvious as B ... 3, H ... l, and D ... l. This effect may be reduced
by generating a longer sequence and keeping only a portion (typically 1/4 to 1/2).
Nevertheless, the FFT transform has been used to produce some of the most striking still
images of random fractals (as shown in the ~olor figures).

Although neither the random cuts nor the FFT filtering represent a recursive method of
generating random fractals, the process of adding more and more frequencies is similar
to the construction of the von Koch snowflake. Fig. 13 shows the increasing complexity
in the FFT construction of a ffim sample with H=O.8 as higher frequencies are included.

The procedure is easily extended to functions of 2 coordinates to generate a fractal sur-

face VH(x,y). VH(x,y) should have the same form for its autocorrelation function as a
ffim. Eq. (11.15) can be extended to the xy plane

where r is a point in the xy plane and 8 is_a displacement in the_xy plane. Since all di-
rectio~ in the xy plane are equivalent, S(k) depends only on 1k 1 = k = (ki + k~)1/2
an~ 18 I. For this 8 dependence.!o correspond to the T dependence of <V(t)V(t+T»,
S(k) must vary as 1/k1+P and S(k) ex: Scut(k)/k. The extra factor of k compensates for
the 2 dimensional diffenptial, 2'ITkdk, in the inte~and. For 3 coordinates, the spectral
density of VH(x,y,z), S(k) ex: Scut(k)/k2. and S(k) varying as 1/k2+P produces a 1/fP
noise for a sample VH along any line.

The corresponding fractal surface is appr-oximated on a finite N by N grid to give

VH(Xn,ym) where xn = nA and Ym = mAo The 2 dimensional complex FFT can be used
to evaluate the series
 824

where kq = q/(NA) and kr = r/(NA) are the spatial frequencies in the x and y directions.
For a fBm surface corresponding to f3 = 1 + 2H = 7 -2D, the coefficients vqr must satisfy

(111.8)

The landscapes of Fig. C7 with 2<D<3 were generated with this process.

For a 3 dimensional fractional Brownian volume, VH(x,y,z), the Fourier coefficients

vqrs will satisfy

(III. 9)

with 3<D<4. This method was used to produce the fractal flakes and clouds of Fig. C9.

random midpoint displacement

Random midpoint displacement is a recursive generating technique that was applied to

normal Brownian motion as early as the 1920's by N. Wiener. It is a natural extension
of the von Koch construction and figures in many of the fractal samples described by
Mandelbrot[I]. Its use in computer graphics has been widely popularized by Fournier,
Fussell, and Carpenter[13,14]. For H #: 1/2 or E > 1, it sacrifices mathematiCal purity
for execution speed in its approximation to mm.

Consider the approximation to a simple mm,VH(t), where the mean square increment
for points separated by a time ~t= 1 is (J2. Then, from Eq. (Ill), for points separated
by a time t, < 1VH(t)-VH(O) 12> = t 2H X(J2. If, for convenience, VH(O)=O, then the
points at t= ± 1 are chosen as samples of a Gaussian random variable with variance (J2 to
satisfy Eq. (Ill). Given these initial conditions, one defines the midpoints at

(lIll0)

where ~1 is a Gaussian random variable with zero mean and variance ~I that is deter-
mined by the condition that the increments from 0 to ± 1/2 must satisfy Eq. (11.1).

(J 2 (J 2 2
2 _(J_[ 1 _ 2 2H -
~1 = 2H
2]
•
(III.ll)
2
 825

The first term is the desired total variance from Eq. (11.1) while the second term repres-
ents the fluctuations already in VH(±I)-VH(O) due to the previous stage. As H+l, ~t
+0, D+ 1, no new fluctuations are added at smaller stages, and VH(t) remains a col-
lection of smooth line segments connecting the starting points. At the second stage,
VH(± 1/4) = 0.5[VH(0) + VH(±1/2)] + ~2 where ~2 has variance

2 ~2 2
~~ = _G_ __1 = _G_[ 1 _ 22H-2l (III. 12)
42H 4 42H

At the nth stage, the length scale has decreased to 1/2n and a random Gaussian variable
~n is added to the midpoints of the (n-l)th stage with variance

(III. 13)

As expected for a fBm, at a length scale r= 1/2n one adds randomness with mean square
variations varying as r2H.

Although' this process does produce a fractal, the result is, unfortunately, not
stationary[I,15] for all H. Once a given point at tj has been determined, its value remains
unchanged in all later stages. All additional stages change t<1j independent from t>1j
and the correlations required of mm with H¢I/2 are not present. More specifically, by
construction, the increments < 1VH(± l)-VH(O) 12> = G 2 . For a stationary process, the
same should be true of all increments with ~t= 1. However, the absence of correlation
across an earlier stage requires that
2
< 1VH(1/2) - VH( -1/2) 12> = 2< 1VH(1/2) - VH(O) 12> = 2 ~H'
2

This gives the desired result G2 only for the H= 1/2 of normal Brownian motion.

Figure 14 shows the generation of a mm sequence with H=0.8. Points generated at

different stages have different statistical properties in their neighborhoods. This often
leaves a visible trace that does not disappear as more stages are added. The effect is more
pronounced as H + 1. These visible artifacts, which are a consequence of the lack of
stationarity of the mathematical approximation, are particularly visible on fractal sur-
faces. Figure 16 shows a zenith view of such a midpoint displacement surface with
H=0.8. In the generation of a midpoint displacement surface on a square grid each step
proceeds in two stages. First the midpoints of each of the squares is determined from its
4 comers and shifted by a random element ~. This determines a new square lattice at 45
 826

degrees to the original and with lattice size 1/.'2. In the second stage, the midpoints of
the new lattice receive a random contribution smaller by 1/(.l2)H from the first stage.
This produces the new square lattice with a scale 1/2 the original. The traces of early
stages are readily visible in Fig. 16. These artifacts, which occur at all stages, can not be
eliminated by local smoothing.

In spite of its mathematical failings, the speed of midpoint displacement, and its ability
to add "detail" to an existing shape make it a useful fractal algorithm for some applica-
tions. To generate N points requires only order N operations. To extend the sequence
by just one point at the smallest scale beyond its original endpoints, however, requires
an additionallog2N operations.

successive random additions

In many respects the non-stationary artifacts of midpoint displacement are similar to the
staircase effect of aliased raster display lines. With midpoint displacement, once deter-
mined, the value at a point remains fixed. At each stage only half of the points are de-
termined more accurately. If one imagines the process of magnifying an actual object,
as the spatial resolution increases all points are determined more accurately. In terms of
the Nyquist sampling theorem, to approximate N real points requires N/2 complex fre-
quencies or N/2 sine and cosine components. When the resolution is doubled to 2N
points, the additional high frequency components alter all of the original values. Mid-
point displacement only adds the additional sine (or cosine) components. Conversely,
in reducing the resolution of an image an anti-aliasing procedure will average over
neighboring pixels. Simply keeping the same value as the center (the equivalent of mid-
point displacement) produces the objectionable staircase edge.

I call the process of adding randomness to all points at each stage of a recursive subdivi-
sion process successive random additions. This enhancement reduces many of the visible
artifacts of midpoint displacement and the generation still requires only order N oper-
ations to generate N points. The computation of the midpoints is the same as midpoint
displacement. The only difference is in the number of random elements. For a sequence
of N elements, N/2 points in the final stage had only one random addition. N/4 points
in the previous stage had 2 random additions. N/8 had 3 and so on. The series converges
to the 2N random additions for the N elements.

The zoom sequence of Fig. 6 was produced with successive random additions. The
artifacts of the square lattice are not as visible as with the midpoint displacement surface
 827

of Fig. 15. Figure 16 shows the construction of a mm sequence with H=O.8 by succes-
sive random additions.

With successive random additions, at each stage all points are treated equivalently. This
has the additional advantage that the resolution at the next stage can change by any fac-
tor r<1. For midpoint displacement r must be 1/2. Thus, given a sample of Nn points at
stage n with resolution A, stage n+ 1 with resolution rA is determined by first interpolating
the Nn+l = Nn/r new points from the old values. In practice, this can be accomplished
using either linear or spline interpolation. A random element ~n is then added to all of
the new points. At stage n with scaling ratio r<l, the ~ will have a variance

A2 ( n)2H (III. 14)

Un a: r .

When l/r is an integer, the generation of a sequence of N points requires order C(r)N
operations. The coefficient C(r) varies as Ln(1 - r)n over the number of stages.

The fractal dimension D for the generated objects is determined only by H, r can be
varied independently. Variations in r change the lacunarity of the fractal. Figure 17 shows
three samples of successive random addition surfaces (all with H=O.8, D=2.2) generated
from differing r. The zoom sequence in Fig. 6 was generated with 1/r=2. As l/r in-
creases to 4 in Fig. 17(a) and 8 in Fig. 17(b) the few characteristic resolutions at which
randomness has been added become visible and the lacunarity increases. As 1/r is de-
creased below 2 to v'2 in Fig. 17(c) the lacunarity decreases and the surfaces approach
the characteristics of the FFT samples. Empirically, a value of 1/r much smaller than 2
produces little observable change.

The free choice of a value for r is an important addition for filmmaking. The addition of
irregularities to a surface can be accomplished continuously as the resolution is slowly
increased from frame to frame.

Weierstrass-Mandelbrot random fractal function

Each stage in a midpoint displacement process increases the highest spatial frequency by
a factor l/r=2. Each stage of successive random additions increases this frequency by
1/r> 1. Thus, both are related to Mandelbrot's generalization[ 1,11] of the Weierstrass
non-differentiable function, VMW(t). Whereas the Fourier series of Eq. (III.6) involves
a linear progression of frequencies, the Weierstrass function involves a geometric
progression. In the notation of this paper,
 828

L AnrnH sin(2'ITr- nt + cf>n),

QC)

VMW(t) = (ill. 15)

n=-oo

where An is a Gaussian random variable with the same variance for all nand cf>n is a
random phase uniformly distributed on 0-2'IT. The original Weierstrass function did not
include the terms for n<O which add small spatial frequencies and large scale fluctu-
ations. As with successive random additions, the addition of a new term to the VMW sum
decreases the spatial resolution by r and adds new fluctuations with variance A2r2H. In
terms of spectral density, although VMW contains only discrete frequencies fn= 1/rn,
each component has a mean square amplitude a: r2Hn a: 1/ f~H in a bandwidth .:1f a: f.
Thus, the spectral density

amplitude 2
M
= (ill. 16)

in agreement with Eq. (11.18).

Although the VMW sum involves an infinite number of components, all practical appli-
cations introduce both low and high frequency cutoffs. For a range of distance scales
from A to Lmax only log(Lmax/A) components are needed. Frequencies much higher
than l/A average out to zero contribution over scales of size A while frequencies much
lower than l/Lmax contribute only an overall offset and slope. The VMW random fractal
function of Eq. (111.15) allows the direct computation of VMW(t) at any t from a stored
set of order log(Lmax/A) coefficients An and cf>n' Figure 18 illustrates this process for
H=O.8 as the number of terms is increased.

With the use of table lookup for sin(x), the VMW fractal function can be made extremely
fast. The fractal is represented by a only few stored coefficients and there is no penalty
for calculating the function outside its original boundary. The resolution can easily be
changed by changing the number of components used. A variable r allows changing
lacunarity independent of D. Other periodic functions can be used in place of sin(x). For
example, both midpoint displacement and successive random additions, correspond to a
triangle function. In addition, midpoint displacement sets all cf>n = constant. Extensions
to E> 1 are possible with periodic functions in all E coordinates.

acknowledgement

This work was made possible by the help, encouragement, and inspiration of Benoit
Mandelbrot.
 829

suggested references

1. Mandelbrot, B.B. The Fractal Geometry of Nature, (Freeman, New York) 1982 and
references therein. See also, Fractals: Form, Chance, and Dimension, W. H. Freeman
and Co., San Francisco (1977).
2. Voss, R F. "l/f (flicker) noise: a brief review", Proc. 32rd Annual Symposium on
Frequency Control, Atlantic City, (1979), 40-46 and references therein.
3. Mandelbrot, B. B. and Wallis, J. R "Some Long-Run Properties of Geophysical
Records", Water Resources Research 5, (1969) 321-340.
4. Musha, T. and Higuchi, H. "The 11f Fluctuation of a Traffic Current on an Ex-
pressway", Jap. J. Appl. Phys. 15, (1976), 1271-1275.
5. Voss, R F. and Clarke, J. "l/f Noise in Music: Music from 1/f Noise", J. Accous.
Soc. Am. 63, (1978),258-263.
6. Voss, RF. and Clarke, J. " 'llf noise' in music and speech", Nature 258, 317-8
(1975).
7. Gardner, M. "White and brown music, fractal curves, and one-over-f noise",
Mathematical Games column in Scientific American, April 1978 p16.
8. Mandelbrot, B.B. and Wallis, J. W. "Fractional Brownian motions, fractional
noises, and applications", SIAM review, 10 (1968) 422-437.
9. For example see: Freeman, J. J. Principles of Noise, John Wiley & Sons, Inc., New
York, (1958), Chapter 1, "Fourier Series and Integrals." or Robinson, F.N.H. Noise
and Fluctuations, Clarendon Press, Oxford, (1974).
10. A good discussion is found in Reif, F. Statistical and Thermal Physics, McGraw-Hill
Book Co., New York, (1965), Chapter 15, "Irreversible Processes and Fluctu-
ations."
11. Berry, M.V. and Lewis, Z.V. "On the Weierstrass-Mandelbrot fractal function",
Proc. R. Soc. Lond. A, 370 (1980) 459-484.
12. Cochran, W. T. et al. What is the Fast Fourier Transform? Proc. IEEE 55, (1967),
1664-1677.
13. Carpenter, L. "Computer rendering of fractal curves and surfaces", SIGGRAPH
'80 Conference Proceedings, (1980) 109.
14. Fournier, A., Fussell, D., and Carpenter, L. "Computer rendering of Stochastic
Models", Comm. of the ACM, 25, (1982) 371-384
15. Mandelbrot, B.B. "Comment on Computer rendering of Fractal Stochastic Mod-
els", Comm. of the ACM, 25, (1982) 581-584, and the response by Fournier,
Fussell, and Carpenter.
 830

Figure 5. Construction of the exactly self-similar von Koch curve. At each stage in the
construction a line segment (on the bottom) is replaced by 4 smaller segments of length
1/3 of its parent (as shown in the middle). The top figure shows the curve after 6 stages
(4096 line segments) with a fractal dimension D = log(4)/log(3) = 1.26 ... .
MAG 1.00~ +00 MAG -4.001:+00 MAG 3.20E:+01

MAG 2.56E:+02 MAG -4.10E:+OJ MAG 6.55E+0-4

MAG I.118E+07 MAG 1.001:+00

Figure 6. Zoom sequence of a statistically self-similar fractal landscape (D=2.2). Each

succeeding picture shows a blowup of the framed portion of the previous image. As the
surface is magnified, a small portion looks similar to (but not exactly the same as) a larger
portion. The total magnification corresponds to 16 million.
 831

log!

1$ r " .f\ aJ' " .l' J' J' J' D;- ~ ;J .J ;J iJ I ~

Iif; ~ J iii .l' J' J'.".,.;/,; ~ iii J "". ;;J'
14 II J iii ~ OJ OJ iii "P iii OJ J'iiI J'.l' J .J F «)

Figure lO. Samples of typical "noises", Figure 12. Samples of stochastically com-
V(t), the random variations of a quantity posed fractal music based on the different
in time. (a) White noise, the most random. types of noises shown in Fig.lO.(a) "white"
(b) 11f noise, an intermediate but very music is too random. (b)" 1/flO music is
commonly found type of fluctuation in na- the closest to actual music (and most
ture, its origin is, as yet, a mystery. (c) pleasing). (c) "Brown" or 1/f2 music is
Brownian motion or a random walk. To the too correlated.
left of each sample is a graphical represen-
tation of that noises spectral density,
Sv(f), a measurement technique for char-
acterizing the time correlations in the noise.
 832

107 107
(Il)
106 106
(c)
105 105
(c)

104 104
(d) t;::
(d)
~ 115
.j 103 .j 103
a
~
(e)
~ 102 102
-fi -fi
'0- '0-
10' 10'

10°

10-'

f (Hz) f (hz)

Pitch fluctuations from different musical Pitch fluctuations n western music

cultures (a) Medieval music up to 1300
(a) the Ba-Benzele Pygmies (b) Beethoven. 3rd Symphony
(b) traditional music of Japan (c) Debussey. piano works
(c) classical ragas of India (d) R. Strauss. ein Heldeniebe
(dl folk songs of old Russia (e) 1he Beaties. Sgt. Pepper
(e) American blues

Figure 11. Spectral density measurements of the pitch variations in various types of music
showing their common correlations as 11f noise.
 833

Freq.
used

1025
513
257
129
65
~
33
17
9
5
3

Figure 13. Increasing detail in a FFT generated sample (N=1024) of fBM for H=O.8,
fJ=2.6, and D= 1.2 as higher frequencies are included. To generate N points requires
order NlogN operations.

points
used
1025
513
257
129
65
~ 33
17

~ 9
5
3

Figure 14. Increasing detail in a midpoint displacement generated sample of fBM for
H=O.8 and D= 1.2 as the number of stages is increased. At each stage a random dis-
placement is added to each of the midpoints of between the previous values. Note that
large fluctuations in an an early stage remain as "special points" in later stages. This
sample is not stationary. To generate N points requires order N operations.
 834

Figure 15. Zenith view of a midpoint displacement surface on a square lattice for H=O.8
and D=2.2. The illumination was along one lattice direction. The non-stationary char-
acter is visible as the prominent shadows. This may be compared with Figs. 6 and 17
which wt!re also generated on a square lattice.

01"- UCI ,:uo .... o L» tA l

Figure 17. Zenith view of a successive random addition surfaces on a square lattice for
H=O.8 and D=2.2. The illumination was along one of the axis. By varying the expan-
sion factor l/r before new randomness is added, the surface /acunarity (or texture) can
be varied without changing D. (a) l/r=4. (b) 1/r=8. (c) l/r=1.41. The lacunarity
increases as l/r increases and for large l/r only a few fluctuation scales are visible. As
r-l the lacunarity approches that of the FFT surfaces.
 835

points
used
1025
513
257
129
65
~ 33
17
9
5
3

Figure 16. Increasing detail in a successive random addition sample of fBM for H=O.8
and D= 1.2 as the number of stages is increased. At each stage a random displacement
is added to each of the points, not just the midpoints. This sample is a better approxi-
mation to fBm than the midpoint displacement of Fig. 14 but the generation requires
order N operations for N points.

used

17
15
13
.e 11
>
9
7
5
3
1

Figure 18. Increasing detail as more components are included in a pseudo-random

VMW fractal function based on the Weirstrass-Mandelbrot random series. Here l/r =
2/(-'5-1) = 1.618 .... Computation of N points covering a spatial irregularity scale from
).. to Lmax involves order Nlog(Lmax/)..) operations.