Note that this is a special case of Equation 11.6 for which c = 0. Under this mapping, all points in the complex z plane fall into one of the following three categories:F(z) = z^{2} [11.12]
2. Request initial point, z1, from the user
3. Repeat until done
3.2 Draw and label line from z_{n} to z_{n+1}
program Complex_Map;
{Program to compute and plot}
{Complex mappings.}
const
scale = 100.0;
size = 400;
type
complex = record
r: real;
i: real
end;
var
Nit, xp, yp: integer;
sd4, sd2: integer;
z, znew: complex;
circle: rect;
out: text;
{**********Prod *************}
procedure prod (a, b: complex; var c: complex);
{Does complex multiplication: c = a bullet b}
begin
c.r := a.r * b.r  a.i * b.i;
c.i := a.r * b.i + a.i * b.r;
end;
begin
rewrite(out, 'Ch11:Complex.dat');
showtext;
writeln('Initial Value of z: (Zr,Zi)?');
ReadLn(z.r, z.i);
showdrawing;
Nit := 1;
ForeColor(cyanColor);
sd4 := size div 4;
sd2 := size div 2;
SetRect(circle, sd4, sd4, 3 * sd4, 3 * sd4);
FrameOval(circle);
ForeColor(blueColor);
drawline(0, sd2, size, sd2);
drawline(sd2, 0, sd2, size);
PenSize(2, 2);
xp := round(z.r * scale + size div 2);
yp := round(z.i * scale + size div 2);
ForeColor(redColor);
WriteLn(' Nit Re(z) Imag(z)');
WriteLn(out, ' Nit Re(z) Imag(z)');
writeln(Nit : 3, ' ', z.r, z.i);
writeln(out, Nit : 3, ' ', z.r, z.i);
moveto(xp, yp);
WriteDraw(Nit : 3);
moveto(xp, yp);
repeat {Iterative loop}
prod(z, z, znew);
z.r := znew.r;
z.i := znew.i;
xp := round(z.r * scale + size div 2);
yp := round(z.i * scale + size div 2);
lineto(xp, yp);
Nit := Nit + 1;
write(Nit : 3, ' ', z.r, z.i);
write(out, Nit : 3, ' ', z.r, z.i);
WriteDraw(Nit : 3);
MoveTo(xp, yp);
readln;
until button;
end.
Output of Complex_Map demonstrates
graphically the dynamics of the complex mapping and the effects of strange
attractors on the trajectories of the points under this transformation.
First, let's examine the motion of a point with modulus < 1. After ten
mappings, the initial point is pulled into the strange attractor at (0,0).
Nit Re(z) Imag(z)
1 7.0e1 5.0e1
2 2.4e1 7.0e1
3 4.3e1 3.4e1
4 7.4e2 2.9e1
5 7.9e2 4.3e2
6 4.4e3 6.8e3
7 2.7e5 6.0e5
8 2.8e9 3.2e9
9 2.4e18 1.8e17
10 3.2e34 8.8e35
11 0.0e+0 0.0e+0
Trajectory of mapping for initial point of modulus > 1.0.
Note effect of strange attractor at •.
1 9.0e1 6.0e1
2 4.5e1 1.1e+0
3 9.6e1 9.7e1
4 1.6e2 1.9e+0
5 3.5e+0 5.9e2
6 1.2e+1 4.1e1
7 1.5e+2 1.0e+1
8 2.3e+4 3.1e+3
9 5.2e+8 1.4e+8
10 2.5e+17 1.5e+17
11 3.9e+34 7.2e+34
12 INF INF
Figure 11.16 indicates the trajectory of a point with modulus > 1. Now the point flees toward • and reaches it after eleven transformations, according to Complex_Map.
Finally, let's examine the chaotic behavior of points in the Julia set itself. These are points lying on the circle of radius one which is the boundary between the basin of attraction of the (0,0) attractor and that of the attractor at • .
Figure 11.17 illustrates a very important feature of chaotic
systemsthe extreme sensitivity of such systems to initial conditions.
In principle, the trajectory of the point should wander through the Julia
set circle indefinitely. In practice, the limited precision of real numbers
prevents this ideal behavior, and the trajectory falls off the Julia set
into the strange attractor. Many natural systems, such as the threebody
gravitational problem, turbulent flow, and the weather, are chaotic systems
subject to extreme sensitivity to initial conditions. This sensitivity
is captured in the parable of the butterfly effectthe disturbance
caused by a butterfly fluttering its wings in China can, in principle,
propagate into a hurricane in the Caribbean Sea!
In the above discussion, the unit circle is shown to be the Julia set for the complex transformation, F(z)= z^{2}. While the example nicely illustrates the concepts of iterative mapping and domains of attraction, the Julia set for this transformation is less than sensational. However, the addition of a single, complex constant to the transformation produces a family of Julia sets, most of which are interesting and some of which are quite dazzling.
The interesting quadratic Julia sets arise for various
choices of c in the transformation:
F(z) = z^{2} + c [11.13]
There are three algorithms for plotting the Julia sets corresponding
to the boundaries of the strange attractor basins at 0 and •. We can identify
these as:
First we discuss each algorithm, then list a Pascal implementation
for the forward mapping algorithm, and finally show output generated by
this program.
(a) 
(b) 
Figure 11.17Initial trajectory of point located on the Julia set (modulus = 1). The initial point is (0.9230769, 0.3846153). Nit Re(z) Imag(z) 1 9.2e1 3.8e1 2 7.0e1 7.1e1 3 8.4e3 1.0e+0 4 1.0e+0 1.7e2 5 1.0e+0 3.3e2 6 1.0e+0 6.7e2 7 9.9e1 1.3e1 8 9.6e1 2.6e1 
However, after about twenty iterations, the trajectory begins to fall off the chaotic Julia set and spiral into the strange attractor at (0,0). By point 30 it has visibly disappeared, and the following computed points document its "last gasp." 30 1.3e3 3.3e3 31 9.2e6 8.4e6 32 1.3e11 1.5e10 33 2.4e20 3.9e21 34 5.5e40 1.8e40 35 0.0e+0 0.0e+0 
The backward mapping algorithm is suggested by an examination of Figures 11.1511.17. Note that all points z _{ n+1 } on a given trajectory are farther from the unit circle boundary (Julia set) than the corresponding z _{ n } point. This suggests that, if we reverse the mapping direction, the Julia set itself will act as an attractor for the trajectory. This, in fact, is the case. That is, by selecting any initial point on the z plane and mapping it backwards, it will be pulled onto the Julia set within the first 1030 iterations and spend the rest of its iteration life hopping from point to point on the set.
Since forward iteration is defined as z_{n+1}
= (z_{n})^{2}+ c, to reverse the direction of motion along
the trajectory, we must write:
z_{n+1} = [11.14]
Note that there are two roots for each iteration of the backward
mapping algorithm. That is, there are two possible precursor points which
can map to any given point under the inverse transformation. The backward
mapping algorithm may be summarized as:
1. Start with an arbitrary point, z 1 , on the complex plane.
2. Repeat for i = 1 to 30 {This gets trajectory onto Julia
set.}
2.2 Select (+) or () root randomly with equal probability
2.3 Set zi = zi+1.
3. Repeat until set is outlined distinctly
{Julia set will "grow in" point by point.}
3.2 Select (+) or () root randomly with equal probability
3.3 Plot point zi+1
3.4 Set zi = zi+1.
The forward mapping algorithm iterates on each point in the vicinity of the Julia set to see if it flees to • (set pixel white) or lies on the Julia set or basin of the zero attractor (set pixel to color). This results in what Devaney calls the filledin Julia set. The algorithm for forward mapping is a relatively simple extension of the Complex_Map algorithm.
1. Request starting parameters
1.2 Input (Cr,Ci), the complex constant, c.
2. For each pixel in the 2 ¥
2 complex plane
2.1.2 IF (zn+1)2 > 10, exit 2.1 loop (pixel = white)
2.1.3 Set zn = zn+1 .
The boundary scanning algorithm combines the forward mapping
algorithm with the formal definition of Julia sets as the boundary between
the two basins of attraction. The basic idea is to color a pixel as part
of the Julia set if it does not flee to • but an adjacent pixel does. This
algorithm can be summarized as:
1. Select an M ¥ N pixel grid to scan a 2 ¥ 2 complex plane.
2. For each (m,n) pixel on this grid
2.1.2 IF zn+1 > 2, color pixel white and exit 2.1 loop.
2.2.2 IF at least one of these four escapes, color (m,n) black
2.2.3 IF all four do not escape, color (m,n) white.
Of the three algorithms, the first and last display the
true Julia set, that is, the boundary between two basins of attraction.
The second algorithm for filledin Julia sets produces more visually appealing
images and is implemented below.
program Julia;
{Program to compute and plot Julia set.}
const
scale = 0.01;
R = 10;
type
complex = record
r: real;
i: real
end;
var
i, j, k, n, row, col, Nit: integer;
x, y: real;
z, znew, c: complex;
done, gone: Boolean;
procedure prod (a, b: complex; var c:complex);
{Does complex multiplication: c = a bullet b}
begin
c.r := a.r * b.r  a.i * b.i;
c.i := a.r * b.i + a.i * b.r;
end;
procedure add (a, b: complex; var c:complex);
{Does complex addition: c = a + b}
begin
c.r := a.r + b.r;
c.i := a.i + b.i;
end;
procedure plot (c, r: integer);
{Procedure to pixel (c.r).}
begin
moveto(c, r);
lineto(c, r);
end;
begin
showtext;
writeLn('How many iterations at each point?');
readLn(Nit);
writeLn('Value of C: (Cr,Ci)?');
ReadLn(c.r, c.i);
showdrawing;
ForeColor(blackColor);
n := 0;
for col := 1 to 400 do
for row := 1 to 400 do
begin
z.r := (col  200) * scale;
z.i := (row  200) * scale;
repeat
prod(z, z, znew);
add(znew, c, z);
gone := (z.r * z.r + z.i * z.i > R);
n := n + 1;
done := (n > Nit);
until done or gone or button;
if done then
plot(col, row);
n := 0;
end;
end.









Julia sets (filledin) by forward mapping algorithm.
Program Julia produces the output for the c values shown in Figure 11.18. A number of interesting features are apparent in this figure. Perhaps the most striking is the rich variety of fractal objects "hidden" in the simple quadratic function, F(z) = z^{2} + c. The approximate selfsimilarity appears in each figurelobes have sublobes which have subsublobes, branches have subbranches which have subsubbranches, and so on. Finally, Julia sets range from formal mathematical patterns, (a) and (b), through suggestive natural shapes, (c) and (d), to the "monsters," (e)  (g), that distressed early mathematicians.
A curious side effect of program Julia
is the apparent dependence of the final fractal shape upon the number
of iterations, Nit,
as indicated in (e)  (g). Conversion of the complex type to double precision
yields identical results, thus eliminating roundoff error as the source
of the side effect. This suggests that perhaps the effect is real, that
is, that the true Julia set is a disconnected set of points rather than
the connected set shown in (e). To verify this assumption, we color encoded
the "flight time," that is, the number of iterations required for the point
to flee to infinity. The results of this run are shown in Figure 11.19.
Julia set for c = (0.4,0.6) using colorcoding, shown as shading, of flight times.
Another method of understanding flight times is through
3D encoding of the iterative map. Figure 11.20 shows the Julia set after
100 iterations as black islands surrounded by an ocean whose depth is inversely
proportional to the number of iterations required to flee to infinity.
Julia set with flight times encoded as 3D inverse depths.
Figures 11.18(e)  (g) correspond to slicing Figure 11.20 at Nit = 10, 20, and 50 respectively, and it is quite obvious that the filledin Julia set is strongly dependent on the number of iterations.
Note that parameters used to represent Julia sets are
the complex z numbers themselves. As we move on to the Mandelbrot
set, the parameters used to represent the set are the complex c numbers.
The Mandelbrot set has been described as the most complex
mathematical object ever discovered. It is generated on the complex plane
as the set of points (mr,mi)
which, when successively mapped with the function,
z_{n+1} =(z_{n})^{2} µ [11.15]
do not cause the complex number z_{n} to fly off to infinity. Equation 11.14 is the form in which Benoit Mandelbrot originally defined this set, but it is clear from comparison of this equation to Equation 11.6 that m = c of our previous formulation.
The algorithm for generating the Mandelbrot set is a minor
modification of the forward mapping algorithm for Julia sets. The main
difference is that the key parameter used to represent the set is the parameter
m (equivalent to c in the Julia set algorithm).
We can summarize this algorithm as:
1. Select a window in complex m space and a viewport in screen pixel space and the appropriate scale and offset parameters for mapping window to viewport.
2. For all pixels in the viewport
2.2 Set z1 = (0,0)
2.3 Repeat
Compute zn+1 =(z_{n})^{2} µ
IF (z_{n+1})^{2}> 10 then gone = true
n = n+1
until gone or n > 100.
2.4 IF not gone then plot pixel corresponding to (mr, mi)
2.5 Set n = 0.
Implementing this algorithm produces the stark, beetlelike
Mandelbrot set representing those complex m
values for which the attractor at • in complex zspace does not pull the
point to infinity under iterated mapping by the function F(z) = z2  m.
A relatively simple modification to the algorithm encodes the "flight time"
that it takes those pixels not in the Mandelbrot set to flee to •. This
code can be used to colorcode points surrounding the Mandelbrot set. Such
color coding provides additional insight into the dynamics of iterative
mapping trajectories and adds greatly to the visual interest of the graphics.
Below we present a Pascal implementation of a colorcoded Mandelbrot set
generator.
{Program to compute and plot the Mandelbrot set.
const
Nit = 100;
scale = 0.005;
R = 10;
type
complex = record
r: real;
i: real
end;
var
i, j, k, n, row, col: integer;
x, y: real;
z, znew, c: complex;
done, gone: Boolean;
procedure prod (a, b: complex; var c: complex);
{Does complex multiplication: c = a bullet b}
begin
c.r := a.r * b.r  a.i * b.i;
c.i := a.r * b.i + a.i * b.r;
end;
procedure sub (a, b: complex; var c: complex);
{Does complex subtraction: c = a  b}
begin
c.r := a.r  b.r;
c.i := a.i  b.i;
end;
procedure plot (c, r, n: integer);
{Procedure to pixel (c.r) in color code n.}
begin
case n of
0..4:
ForeColor(blueColor);
5:
ForeColor(cyanColor);
6:
ForeColor(greenColor);
7:
ForeColor(magentaColor);
8..11:
ForeColor(redColor);
12..20:
ForeColor(yellowColor);
21..99:
ForeColor(whiteColor);
100:
ForeColor(blackColor);
otherwise
ForeColor(blackColor);
end;
moveto(c, r);
lineto(c, r);
r := 400  r;
moveto(c, r);
lineto(c, r);
end;
begin
for col := 1 to 400 do
for row := 1 to 200 do
begin
z.r := 0.0;
z.i := 0.0;
c.r := (col  100) * scale;
c.i := (200  row) * scale;
repeat
n := n + 1;
prod(z, z, znew);
sub(znew, c, z);
done := (n > Nit);
gone := (z.r * z.r + z.i * z.i > R);
until done or gone or button;
plot(col, row, n);
n := 0;
end;
end.
The Mandelbrot Set (black) with colorcoded flight times.
The output of program Mandelbrot is shown in Figure 11.21. The labeling was added to the figure by an image processing program.
The fascination inspired by the Mandelbrot set stems from several sources. This set exhibits all of the previously discussed features distinguishing fractals, including an infinite complexity which emerges as the researcher explores selected regions in finer and finer detail. With a program no more complex than program Mandelbrot, anyone can soon be discovering new regions of the set which no human has ever seen before. Such a voyage of discovery is documented in the captivating video animation, Nothing but Zooms, produced by the Cornell University Supercomputer Center.
Let's explore the Mandelbrot set by successively magnifying
regions near the boundary that look interesting. Using the more conventional
mapping notation, F(z) = z ^{2} +
c, we can examine a smaller region of the complex c plane defined by diagonally
opposite corners, c1 and c2,
as shown in the following figures. (Note that the convention using c positive
along the real axis effectively mirrors Figure 11.21 about the imaginary
axis.)
Segment of the Mandelbrot set with corners at
c1 = (0.64167, 0.6930) and c2= (0.453,0.511).
Segment of the Mandelbrot set with corners at
c1 = (0..5963,0.5821) and c2 = (0.5550,0.5390).
Segment of the Mandelbrot set with corners at
c1 = (0.5624,0.5623) and c2 = (0.5572,0.5570).
In just three successive magnifications we are already examining, in Figure 11.24, a region ~7¥10^{6} the size of the original Figure 11.21. This exploration may be continued indefinitely and soon reveals elegant fractal patterns never observed before.
Visualizing the Mandelbrot set in 3D provides additional
insight into the fractal structure of this object. Figure 11.25 shows the
same region of Figure 11.23 as a shaded relief map illuminated from the
upper right. The complete Mandelbrot set is shown as a 3D object in Figure
11.26. This representation is helpful in visualizing the set itself as
a region of stable trajectories and all points off the set as "sliding
off to infinity."
Segment of the Mandelbrot set with corners at
c1 = (0.5624,0.5623) and c2 = (0.5572,0.5570)
displayed as a shaded relief map.
Mandelbrot set rendered as a 3D "island" of stability. The greater the depth off shore, the faster the iterated z point flees towards infinity.
To conclude the discussion of complex plane mapping, we demonstrate the relationship of the Mandelbrot set to Julia sets with Figure 11.27. Recall that Julia sets are represented on the z plane while the Mandelbrot set uses the complex number, c, as its basis. Different Julia sets are distinguished by different values of c. This implies that a Mandelbrot set displayed as an m ¥ n pixel image is really an index or catalog of m ¥ n Julia sets for which the complex c serves as the key.
To examine elements of this catalog, we can select various c values on the Mandelbrot set and compute the corresponding Julia sets. Figure 11.27 illustrates this procedure, indicating a number of points on the Mandelbrot set and their corresponding Julia sets (compressed to twentyfive percent of their full linear scale). An interesting aspect of this figure is that some regions show a rather intuitive similarity in structures (e.g., filamentary regions of the Mandelbrot set giving filamentary Julia sets) while other regions show surprising differences in structure.
This concludes our discussion of fractals generated by iterative mapping on the complex plane. We have concentrated on the Mandelbrot set and Julia sets of the mapping F(z) = z^{2} + c for two reasons. First is the key role they played in the historical development of fractal geometry. Second is their intrinsic beauty, mathematical simplicity, and awesome complexity.
The reader should not, however, get the impression that this introduction to complex plane mapping exhausts the topic. The most obvious extension is to study the Julia sets of other complex analytic functions. Two classes of functions which have been studied extensively include:
Study of the dynamic behavior of iterated complex mappings
(i.e., trajectories) continues to yield information on the location
of strange attractors, their basins of attraction, and splendid fractal
images. Research results emerging from these studies have been extended
to understanding the stability of dynamical physical systems such as the
solar system.
The Mandelbrot set as catalog of Julia sets. The Julia set in the lower, righthand corner of the figure, for example, corresponds to (c_{r},c_{i}) = (0.13,0.64).
The single characteristic shared by all of the fractals discussed to this point is that they are deterministic. That is, running the fractal generating algorithm again will result in a structure identical with that generated in the first run. However, one of the distinguishing features of many natural objects (e.g., coastlines, clouds, and mountains) is that they are statistically selfsimilar but unique. In order to examine this vary large class of fractal objects we must introduce the concept of stochastic (or random) fractals.
The distinguishing feature of stochastic fractals is that the dominant behavior of the fractal structure depends on random processes. Thus, no two coastlines, river basins, or mountain skylines will ever be identical, although certain small segments may closely resemble each other.
Stochastic fractals are particularly significant for the following reasons. First, landscapes generated as stochastic fractal images are so convincing in their depiction of natural scenes as to leave no doubt in the mind of the observer that they accurately represent nature. The logical implication of this observation is that natural objects such as mountains, river basins, and the cratered surface of the moon must obey fractal geometry. A wealth of recent research in physics, chemistry, meteorology, and related scientific disciplines has confirmed the fractal nature of a great variety of natural objects.
This is particularly important for computer graphics because of the principle of model authenticity. Model authenticity requires that realistic images of computergenerated (i.e., artificial) objects must be based on models which are consistent with "the real thing." This means that the most realistic renderings of natural scenes can be achieved using stochastic fractal models.
The second significant aspect of stochastic fractals is that they have provided the first real commercial application of fractals. As we have already indicated, several science fiction movies have already used fractal landscapes and firestorms as the natural background for their imaginary worlds. As fractal algorithms improve and fractal ASIC hardware appears, this trend will most certainly accelerate.
The two examples of stochastic fractals we consider are
Brownian motion and fractal landscapes. The first example represents one
of the earliest random processes observed in physics, and the second illustrates
the basic algorithms used in generating the imaginary worlds of the movie
industry.
In 1828, Robert Brown, a botanist, noted an irregular zigzag motion of pollen grains floating on the surface of liquids. This effect, now known as Brownian motion, was fully explained in one of Einstein's famous 1905 papers. Einstein's theory and Jean Perrin's exhaustive set of measurements of the motion of particles of various sizes in different fluids firmly established the kinetic theory of matter.
Brownian motion is an important example of the general class of problems called the randomwalk problem. The randomwalk problem, in its simplest 1D form, involves the motion resulting from the flip of a coin. If it's heads, the coinflipper takes one step ahead; if it's tails, s/he steps back one step. The resulting motion is called a randomwalk. In the case of the Brownian motion of floating particles (2D) or particles suspended in fluids (3D), the motion is caused by forces exerted on the particles by collisions with molecules of the fluid. The statistical fluctuations in the rate of collisions from various directions result in unbalanced forces causing the particles to skitter across the surface of the liquid.
To illustrate Brownian motion in 1D we present two algorithms,
the first a randomwalk algorithm and the second a midpoint displacement
algorithm. Both algorithms are expressed in terms of a displacement along
one axis, d(t), as a function of time, t.
1. Set d(0) = 0 for time t = 0.
2. For t = 1 to 256
2.1 Flip a coin with a balanced random number generator, ran
If heads (i.e., ran > 0)
d(t) = d(t1) + step
else
d(t) = d(t1)  step
2.2 Plot (t, d(t)).
The Pascal implementation of this algorithm is virtually
verbatim, with the only extensions being the balanced random number generator,
ran, and the pen size of 2¥2
pixels to mask the intrinsic jitter built into the algorithm. This implementation
produces a trace of the 1D Brownian motion in which the displacement
is plotted along the yaxis and time along the xaxis. A
selected example of the output of Random_Walk is
shown in Figure 11.28.
program Random_Walk;
{Program to generate Brownian trace}
{by RandomWalk algorithm.}
const
xmax = 256;
ymax = 200;
step = 2;
var
x, d, t: integer;
i, xp, yp: integer;
flip: real;
time: DateTimeRec;
{********** ran **********}
function ran: real;
{Function to return a random real number}
{on range 1 < ran < +1.}
begin
ran := random / 32768
end;
begin {Main program}
pensize(2, 2);
d := 0;
{Get "Random" randomized.}
gettime(time);
randSeed := time.Minute * 60 + time.Second;
for i := 1 to xmax do
begin
flip := ran;
if flip > 0 then
x := step
else
x := step;
d := d + x;
xp := i;
yp := ymax div 2  d;
drawline(xp, yp, xp, yp);
end;
end.
Mandelbrot has observed that the displacement trace in 1D Brownian motion is reminiscent of the skyline of a mountainous scene. An interesting feature of Random_Walk is that to get an "interesting" profile, like that shown in Figure 11.28, the user must filter out a number of "dull" runs showing much less relief. This sorting for appealing output may appear to be cheating yet is similar to what sightseers do on vacations. Millions of visitors each year all travel thousands of miles through relatively dull countryside to marvel at the prominent relief of the Teton mountains of Wyoming. Both Nature and Random_Walk produce more dull than interesting scenes.
The second algorithm, random midpoint displacement (RMD), grew out of Norbert Wiener's studies of Brownian motion in the 1920s. The basic concept is similar to the linear replacement mapping algorithm we used to generate von Koch fractals, with a random displacement of the center point of the line substituted for the deterministic replacement pattern of the von Koch fractal. This algorithm can be summarized as follows.
1. Initialize parameters
1.2 Set number of iterations, Nit = log2(Tmax)
1.3 Set variance sigma = Tmax/2 {arbitrary choice of scale factor.}
2. Repeat for n = 1 to Nit {i.e., subdivide line until we hit pixel level}
2.1 Sigma = sigma/2 {reduce variance in half}
2.2 Repeat for all line segments at this level
2.2.1 Divide line segment in half, creating two new line segments
2.2.2 Displace midpoint by ran¥sigma, where 1 < ran ? +1.
2.2.3 Plot displaced point.
The Pascal program, MidPoint,
is a straightforward implementation of this algorithm in which the main
program handles the bookkeeping and the procedure RMD
subdivides each line sent to it and applies the random midpoint displacement.
To conserve space we have also deleted the function ran and its
initialization.
program MidPoint;
{Program to generate Brownian trace by}
{method of Random Midpoint Displacement.}
const
xmax = 256;
ymax = 200;
type
image = array[0..xmax] of real;
var
trace: image;
n, m, Nit, d, Mmax: integer;
sigma: real;
time: DateTimeRec;
{*********** RMD ********}
procedure RMD (var trace: image;
i, j: integer; sigma: real);
{Function to apply random displacement}
var
k: integer;
dy, ave: real;
begin
k := (j  i) div 2;
dy := sigma * ran;
ave := (trace[i] + trace[j]) / 2;
trace[i + k] := ave + dy;
moveto(i + k, ymax div 2  round(trace[i + k]));
lineto(i + k, ymax div 2  round(trace[i + k]));
end;
begin {Main Program}
pensize(1, 2);
sigma := xmax div 2;
Nit := round(ln(xmax) / ln(2));
d := xmax;
trace[1] := 0;
trace[xmax] := 0;
for n := 1 to Nit do
begin
sigma := sigma / 2;
Mmax := round(exp((n  1) * ln(2)));
for m := 1 to Mmax do
RMD(trace, (m  1) * d, m * d, sigma);
d := d div 2
end;
end.
Brownian trace generated by Random Midpoint Displacement program, MidPoint.
While MidPoint effectively implements the algorithm shown, it efficiently bypasses the problem associated with bisecting and storing lines by taking advantage of the fact that, at each level, lines are defined by the end points which have been plotted at that level. These end points are fixed, and successive midpoint bisections operate only on the intervals between the fixed end points. By successively bisecting the trace down to the pixel level and plotting each midpoint as it is generated, MidPoint generates the Brownian trace without resorting to any line plotting.
Mandelbrot has generalized the technique of random midpoint displacement to what he calls fractional Brownian motion (fBm) and applied fBm techniques to the study of river discharges, coastlines, distribution of natural resources, and mountain terrain. One of the leading advantages of the random midpoint displacement algorithm over simpler techniques like the randomwalk method is that it is more easily extensible to generating fractal landscapes. For this reason it was the first and has remained the favorite algorithm for creating artificial scenes.
The simple random midpoint displacement
algorithm presented here suffers from flaws which prevent precise simulation
of Brownian motion. The most obvious of these is the fixed nature of each
midpoint displacement after it is generated. While this flaw is not apparent
in 1D Brownian traces, it causes "creases" in 3D Brownian landscapes. Richard
Voss has indicated a refinement, successive random additions, which treats
all points equivalently and eliminates the crease problem .
(a) (b)
Random midpoint displacement algorithm for generating fractal landscapes. Midpoints are displaced randomly in sign and magnitude in the vertical directionhere 2 and 3 upwards and 1 downwards. The process is repeated on each of the four resulting triangles.
The random midpoint displacement algorithm is readily extended to generate fractal mountain landscapes. The basic idea is to start with a planar 2D figure such as a quadrilateral or triangle and successively subdivide each side, randomly displacing its midpoint as shown in Figure 11.30.
The RMD algorithm suggested by Figure
11.30 is a 3D extension of the 1D RMD algorithm for Brownian motion. It
can be summarized as follows.
1. Initialize starting parameters
1.1 Define initial triangle as three (x,y) pairs in the z = 0 plane
1.2 Define initial variance as s = <L>/2, where <L> = average side length
1.3 Define light source direction
1.4 Define camera position
2. Repeat until desired resolution is obtained
2.1 Reduce s = s/2
2.2 For all triangles at this level of refinement
2.2.1 Compute midpoints of each side
2.2.2 Displace each midpoint by Dz = ran¥s
3. Project resulting 3D model onto a 2D image
3.1 Use back face removal to eliminate invisible surfaces
3.2 Sort visible triangles by distance from camera
3.3 Project triangles onto image with Painter's algorithm
3.4 Shade each triangle with simple shading model
As might be expected, the generation
of the 3D model of the fractal landscape requires less computation than
rendering the model. One additional problem introduced by the 3D extension
of the 1D RMD algorithm is the increased complexity of bookkeeping. The
complexity is not apparent in the first iteration transforming the planar
triangle (a) into the fourtriangle profile of (b). The problem arises
at the next level of bisection, however, when step two of the algorithm
is applied to each of the four resulting triangles in (b). If the algorithm
is first applied to the leftmost of the four triangles of (b) and then,
independently, to the center triangle, in general a tear or "cave" in the
surface will develop. This results due to the fact that the midpoint of
the boundary edge has undergone two independent displacements. This unwanted
sideeffect violates the physical constraint that the surface must be continuous.
(a) N_{it} = 1 
(b) N_{it} = 3 
(c) N_{it} = 5 
(d) N_{it}= 7 
(e) N_{it} = 7 (shaded) 
(f) N_{it} = 7 (shaded, with sea level = 0) 
Generation of fractal mountain landscape by RMD algorithm.
Two possible solutions to this problem include:
To achieve visual realism of fractal mountain landscapes, this simple RMD algorithm is usually augmented with a color scheme typically encoded by the zvalue of displacement. Triangles with the highest elevations may be colored white to simulate snow; the next lower elevation band is given some mix of gray/brown/yellow to simulate rock structure; the band beneath that may be shaded various colors of green to indicate forest foliage; and the lowest levels may be colored a uniform blue to simulate a lake at the bottom of the valley. Such a heightencoded color scheme achieves an impressive degree of visual realism.
In Figure 11.31 we present the output of one simple RMD program.
Figures (a)  (d) represent wire frame output for the
odd levels of refinement from level 1 through level 7 which is approaching
the pixel resolution of the screen. The camera viewing angle was selected
to be 30° with respect to the z = 0 plane. Figure (e) also used 7 levels
of iteration, but was rendered as a shaded model. An interesting feature
of this particular program is a "seal level" control which was set at z
= 0 in Figure (f). This mode allows the generation and study of fractal
coastlines and islands.
We have already indicated that the first major application of fractal geometry was the generation of artificial landscapes for the imaginary worlds of the movie industry. There are two additional applications which should prove to be of even more immediate interest and value to computer graphics users.
Fractal design tools open up new opportunities for designers
and artists to produce complex and random patterns with more appeal than
output from conventional geometry programs. Fractal image compression techniques
promise to resolve many of the problems associated with the enormous storage
demands imposed by graphical images. Products in both of these applications
areas are commercially available.
Many programs have been written to generate and explore particular types of fractals. The public domain MacFractal generates Sierpinski triangles, Brownian traces, and fractal mountain landscapes. The Beauty of Fractals Laboratory provides elegant tools for investigating the Julia and Mandelbrot sets. To the best of our knowledge, however, no single program has appeared which generates all of the standard fractal types.
Several commercially available painting programs now incorporate
generators for certain fractal shapes. We illustrate here three examples
 a Brownian trace generator and two fractal tree generators.
The program PixelPaint Professional contains a
fractal line option with a user controlled "randomness" parameter shown
in Figure 11.32. With the randomness set in the range 90%95%, the line
tool generates a good approximation to a 2D random walk of a Brownian particle
(shown in the sample window).
(a) 
(b) 
(c) 

(d) 
(e) 
Fractal trees generated using Kid Pix.
Turning the randomness down produces a less jagged line more closely resembling a 1D Brownian displacement trace. Although the design process rarely calls for a single fractal line, this tool is useful in designing irregular natural phenomena such as surfaces cracks and lightning. In Figure 11.33 we show a lightning stroke scene constructed with the fractal line generator.
(a) 
(b) 
(c) 
(d) 
(e) 
Fractal trees generated using PixelPaint program.
Fractal plant generation is the subject of intense research which is yielding interesting results on plant morphology and realistic images. The charming children's program, Kid Pix, contains a stochastic fractal tree generator capable of producing interesting effects. Upon selection of the tree icon, each click of the mouse produces a random tree of varying sizes and shapes at the cursor position as shown in the top row of Figure 11.34. The user has two control options: the color and tree size category. The option key switches the generator from standard to large size trees shown in the bottom row.
As Figure 11.34 indicates, even simple fractal tree algorithms can produce trees of great variability suggestive of natural objects. The primary limitation of the fractal trees of Kid Pix is the lack of leaves. The following program simulates leaves by terminating each terminal branch with a branch cluster.
An undocumented feature of the regular PixelPaint program includes an elegant fractal tree generator. Upon selecting Special Effects + pencil icon + option key, each click of the mouse button produces a fractal tree such as the ones shown in Figure 11.35.
Note the interesting 3D effect generated by the combination of outlining the left side of each branch with black and the sequential generations of the branches. In addition, this algorithm varies the length of branches randomly, enhancing the realism of the resulting tree images. This algorithm fails the model authenticity test primarily through the uniform color of trunk, branches, and leaves. In Figures 11.35(d) we attempted to overcome this problem by dragging the trunk portions of a brown fractal tree onto a green fractal tree.
Figure 11.36 shows how fairly realistic scenes can be designed by generating fractal trees on a smoothly shaded background. The illusion of a 3D scene is created through use of the Painter's Algorithm to creating a background and add trees sequentially from the background to the foreground. This technique, in combination with the intrinsic 3D appearance of individual trees, creates an impression of 3D structure in the scene.
Fall scene designed using PixelPaint fractal trees.
These examples illustrate the design potential of simple stochastic fractal tools. The obvious extension of these techniques is to combine in a single paint/drawing program the following stochastic fractal forms:
Such a program would provide the designer of natural scenes
the same power and flexibility that sophisticated CAD systems provide for
the design of conventional geometric objects.
As the integration of the graphics areas of animation, video, and multimedia progresses, effective image compression techniques are assuming increasing importance. The fractal image compression technique invented by Michael Barnsley is probably the most significant application of fractals in computer graphics. The history of the evolution of this technique is a fascinating example of harnessing abstract mathematical concepts to solve difficult technological problems.
To illustrate the potential for image compression using fractal techniques, consider the fractal fern shown in Figure 11.13 and the set of twentyeight iterated function system (IFS) coefficients shown in Table 11.2 used to generate the fern image. Assume your assignment is to transmit the image of the fern to a colleague overseas. To keep your colleague's options for further image processing or electronic publishing open, you decide on magnetic media (disk) or, if both sites support it, email. Then the question arises: What format provides for most efficient transmission?
One option is to send the image in PICT or TIFF image file format. Depending on the desired resolution, this option would require between 32 Kb (for the image shown) to several Mb. A second option, suggested by Table 11.2, is to send the 28 real coefficients from which the fractal image was generated. This corresponds to 112 bytes of information. Particularly if you are using email and a slow modem, the choice is obvious. That is, the IFS fractal encoding of this image achieves an image compression ratio of between 285 and 10,000 or more.
But, you protest, those twentyeight numbers by themselves are not sufficient for generating the image. This is a proper objection, and the complete solution is to send along the rules for interpreting the coefficients. This means sending along the Pascal program, Fern.p, at a cost of some 1,919 bytes. Even sending along this more general program (capable of interpreting other IFS fractal images as well), you are achieving an image compression ratio of over 1,500 at the higher resolution.
Even more impressive results emerge as we examine how the image compression ratio depends on resolution. For PICT and TIFF files, doubling the resolution corresponds to quadrupling the image file size. However, the IFS coefficients are capable, in principle, of infinite resolution. So as the trend towards higher and higher resolution graphics images continues, the advantage of fractal image compression techniques increases.
So the advantages of sending the fractal fern image in
compressed IFS format is beyond question. However, the image you wish to
transmit to your colleague may not be a fern, tree, or Sierpinski triangle,
all of whose coefficients are tabulated in Barnsley's book, Fractals
Everywhere. Suppose the image you wish to transmit is a captured video
frame, a scanned image of your grandmother, or an elegant design that you
created with a paint program. What IFS coefficients do you email to your
colleague? The answer to this question involves solving what Barnsley defines
as the inverse problem.
When Barnsley began his research into fractals, the predominant research approach was to postulate a new linear replacement map, IFS, complex map, or stochastic process and then investigate the nature of the resulting fractal image. This research approach is represented as the top row of Figure 11.37.
Barnsley formulated (and later solved) what he called
the inverse problem:
The inverse problem can be conceptualized as the bottom
row of Figure 11.37.



(a) The fractal research agenda



(b) The inverse problem
Traditional fractal research vs. the inverse problem.
As Figure 11.37 indicates, the traditional fractal research agenda can be stated: Given an algorithm, what is the image? The inverse problem can be stated: Given an image, what is the algorithm?
The Collage Theorem of Barnsley is the key to solving the inverse problem. The Collage Theorem says, in effect, that an image may be represented to any desired degree of accuracy by a union of contractive affine transformations of itself. This is quite a mouthful, so let's review what some of the terms mean and then illustrate it with examples.
Affine transformations include the scaling, rotations,
and shears we studied as homogeneous, 2D transformations. Contractive means
that any two points on the transformed image are closer together than they
are on the original image. As the first example of the Collage Theorem,
assume we wish to generate the IFS code for a unit square as shown in Figure
11.38. The four affine transformations, W1,
W2, W3,
and W4 shown in Table 11.3 map
the original unit square into the four labeled rectangles which serve as
the collage for representing the original square.
(a) 
(b) 
Collage Theorem used to encode a square as a 4transformation IFS.
The coefficients in Table 11.3 can be written down from a careful examination of Figure 11.38(a), interpreted using equation 11.4. For example, transformation W3, corresponding to rectangle 3, requires scaling the unit square by 0.25 in x and 0.5 in y and shifting it 0.5 units in y (coefficients a, d, and f). Since its area is one eighth of a unit, we assign its probability 0.13. The coefficients for the other three contractive affine transformations may be similarly read off by examination of the collage (a).
Running program IFS with the coefficients of Table 11.3 will, indeed, generate the unit square, our original image. Thus we see how the Collage Theorem provides the critical tool for mapping an image into the coefficients of the IFS algorithm.
The collage of a square from the four rectangles of Figure 11.38 provides a "perfect fit," with no overlap nor uncovered areas of the target image. In comparison with natural objects like clouds, mountains, and trees, this example appears a bit contrived. How could the Collage Theorem be applied to generate the fractal algorithm for an irregular natural object like a tree or leaf?
Barnsley proposes the Collage algorithm as an interactive
geometric modeling tool for using the Collage Theorem to find the IFS code
for an image. The algorithm defines the steps necessary to transform a
target image, T, into a set of contractive affine transformations,
{Wn, pn:
n = 1,2,Ö,N}, which specify the attractor image,
A, which closely represents T.
1. Process the target image, T, to generate a closed polygon approximation.
2. Repeat for n = 1 to N until the collage, T´, adequately represents the target T
2.1 Duplicate image T to produce image Dn
2.2 Contract and transform Dn until it fits an unfilled niche in T
2.2.1 Drag Dn into the niche in T
2.2.2 Record the transformation, Wn = {an,bn,cn,dn,en,fn}
2.2.3 Compute and record pn = area(Dn)/area(T)
The collage, T´, is formally defined as the union of the Wn transformations:
[11.16]
where
Wn = nth tile of the collage.
Barnsley provides a measure, h(T,T´) of the "goodness of fit" between the collage, T´, and the original image, T. If h(T,T´) is small, the collage is a good fit to the target, and the Collage Theorem guarantees that the attractor, A, defined by the set Wn will also give a good fit to T. The goal of the Collage algorithm is the achieve the maximum goodness of fit with the minimum number of collage tiles. This corresponds to simultaneously minimizing both h(T,T´) and N.
Let's use an image of a real natural object to illustrate the Collage algorithm. Figure 11.39(a) is the scanned image of a maple leaf from the author's own backyard. Figure 11.39(b) represents the output of three image manipulation processes:
1. Convert (a) to a 2color level image (i.e., all black interior),
2. Apply an image outline tracer,
3. Convert the bitmapped outline to a series of polygon sections and group them.
The output of step (3) is a polygonal object that is easily manipulated by any drawing program. Steps 2.2 and 2.2.1 of the Collage Algorithm were performed using the mouse and various image manipulation options of a drawing program to produce the collage of Figure 11.39(c). Steps 2.2.22.2.3 were omitted for the purpose of this illustration but could be accomplished by careful, quantitative use of standard drawing programs.
An undocumented feature of the drawing program displayed a visual "goodness of fit" indicator shown in Figure 11.39(d). First, the five tiles representing T´ were grouped, removed from the target, and colored black. Next, the target leaf, T, was colored green. Finally, the black collage was dragged over the green target, and, unexpectedly, the area in common between the two images turned pink. Thus, with N fixed at five, we visually maximized the goodness of fit by maximizing the pink area. This corresponds to minimizing h(T,T´) by minimizing the sum of black and green areas.
A little reflection reveals that the Collage Theorem is both mysterious and powerful. Note that in neither the procedure outlined for building the square collage of Figure 11.38 nor that for the leaf collage of Figure 11.39 is there any specification of the shape parameters of the original image. That is, the strange attractor, A, defined by the IFS coefficients of Table 11.3 "knows" that it is a square, and Wn recorded for the leaf transformation leading to Figure 11.39(c) "know" that they should regenerate not just any leaf, but this particular leaf. Both of these IFS codes have encoded the shape of the original object without recording a single absolute coordinate.
The power of the Collage Theorem is exploited in the commercial image compression system we turn to next.
The two illustrations shown for solving the inverse problem using the Collage Algorithm range from a perfect solution for the square to a messy first approximation for the leaf. Using more care, a drawing program user could undoubtedly achieve a better collage for the leaf than the one shown. In addition, using a few more tiles to fill the vacant niches in T would significantly improve the goodness of fit.
However, the reader must sense that applying the Collage
Algorithm manually as we did in these two examples is not an efficient
image compression technique. The question then becomes: Can the process
be automated? That is, is it possible to build a system which takes raw
images as input and produces a set of IFS codes from which the raw image
can be reconstituted? The answer, happily, is yes!
(a) 
(b) 
(c) 
(d) 
The Collage Algorithm. Generating IFS code for a leaf.
Barnsley conceived a system for automatically encoding images as IFS coefficients, and he and Alan Sloan, a mathematics colleague at the Georgia Institute of Technology, have established Iterated Systems, Inc. for developing and marketing image compression products using this technology. While the details are proprietary, the outline of the process are given in the Barnsley references above.
The fractal encoding image compression process is summarized in the following steps.
In addition to the tremendous compression ratios, the fractal image compression technique has other advantages. First is its stability. As the collage becomes more accurate, the resulting IFS code better represents the original image, but the code does not need to be exact in order to get acceptable image encoding. Second, the IFS code is robust. That is, small perturbations in the code will not seriously damage the image.
What problems remain with the fractal image compression
technique? The principle problems are associated with the computationally
intensive nature of the encoding and decoding phases. Barnsley and Sloan's
early results indicated that complex color images require about 100 hours
to encode and 30 minutes to decode using software on a dual68020 system.
This would appear to rule our fractal image compression as a feasible technique,
at least for personal workstations.
A highly successful approach to solving the bottlenecks for specific applications in computer science is to work out the solution in software and, after the algorithm is tested and proven, convert the algorithm to firmware. This firmware takes the form of Application Specific Integrated Circuits (ASICs). ASICs can be considered as special processors dedicated to solving specific tasks at high speeds. The fractal image encoding task is a natural candidate for conversion to ASICs.
Iterated Systems has developed a fractal image compression board which they market under the name P.OEM™ ("Pictures for OEMs"). The P.OEM™ board consists of eight fractaltransform ASICs and an Intel i960 RISC processor. Image playback (decoding) is performed by improved decompression software modules. Both the P.OEM™ board and image decompression software is compatible with 80286 and higher PCs.
The features and performance specifications of the P.OEM™ system are summarized below.
Iterated Systems offers Fractal Factory™ Video Player
for decompressing video tape sequences and playing them back on 286/386/486
VGA systems. It also offers Fractal Factory™ Slide Projector for
displaying up to sixty high definition slides stored on a single 1.2 Mb
disk. For users wishing to use the Fractal Factory™ playback modules
but not willing to buy the P.OEMô image compression board, Iterated Systems
provides image compression services. To give some indication of the quality
of the fractal image compression/decompression process, we show as scanned
image of a promotional photograph from Iterated Systems, Inc. in Figure
11.40. This image was scanned at 150 dots/inch as a TIFF file and cropped
slightly to reduce the file size.
Scanned image of Iterated Systems photograph indicating the compression ratio possible with fractal image compression. The original photograph indicates no image degradation due to the compression/decompression cycle.
From these specifications it is apparent that fractal
image compression is an established technology with great commercial potential.
It has the promise of reducing image storage and retrieval tasks to a time
and space scale similar to that required for text documents. This technological
advance has serious implications for all aspects of multimedia, with particular
relevance to the integration of video information on personal workstations.
Perhaps the most important lesson to be learned from this discussion is
that even the most esoteric mathematical formalism such as IFSs and the
Collage Theorem can yield practical technological products in a very short
time span.
As objects of fractional dimension, fractals provide a fascinating conclusion to our discussion of Ndimensional graphics. Computer graphics played a key role in the discovery of fractal geometry, and fractal geometry has returned the favor by providing important tools for the design of natural objects and the compression of graphics images. Fractal geometry illustrates again the importance of the principle of model authenticity to computer graphics. Since the structure of most natural objects is fractal, synthetic images of such objects must incorporate fractal geometry models in order to achieve visual realism. The movie industry has recognized this principle and now routinely uses fractal geometry as the basis for its imaginary landscapes.
Fractal geometry is a manifestation of the dynamic behavior of chaotic systems, a huge class of natural systems ranging from the neural network activity of the brain to motion of galaxies. Iterated function systems (IFSs) were shown to be capable of representing objects of arbitrary shape, and examples were presented for geometric shapes such as squares and natural shapes such as leaves. Finally, the solution of the IFS inverse problem was demonstrated as an effective technique for fractal image compression. Additional readings on fractal geometry are listed in this reference.