Chapter 11Fractals Objects of Fractional DimensionMorris Firebaugh© William C. Brown Publishers, Inc. 
We have seen that one of the best measures of visual realism in computer graphics is the extent to which shading models and rendering algorithms incorporate the physics of light and its interactions with objects. Yet, in spite of the successes of ray tracing and radiosity techniques, experts in computer graphics can consistently detect the artificial nature of synthetic images (i.e., computergenerated images). Even novices in computer graphics often sense something artificial about synthetic images, with the objects themselves providing clues on the rendering technique. Glass balls, teapots, and checkerboards suggest ray tracing while walls, tables, and lamp shades hint at radiosity. It is rare, indeed, to find a synthetic image of natural objects that is at all convincing in its visual realism. That is, the representation and rendering techniques we have discussed to date are extremely effective in displaying objects that are humanmade, symmetrical, or geometric primitives, but they fail miserably in generating realistic images of such common objects as trees, clouds, and mountains. Why is that?
The answer is both simple and profound, and it is deeply rooted in the concept of dimensionality itself. The simple answer is given by Yale mathematician Benoit Mandelbrot in the introductory quotation who observes that most natural objects do not conform to simple geometric shapes. The more profound answer is that the geometry best describing most natural objects is not traditional 3D Euclidean geometry, but rather what Mandelbrot describes as the fractal geometry of nature. Mandelbrot coined the term fractal to describe both natural and mathematical objects that share the property of having fractional dimensionality, a term we define shortly. The concept of dimensionality plays a key role in distinguishing fractal objects from those 1D, 2D, and 3D objects we have studied so far. For this reason, the study of fractals as fractionaldimension objects is the natural concluding chapter of our NDimensional approach to computer graphics.
Fractal geometry is another confirmation of the principle
of model authenticity introduced in previous chapters. This principle
holds that realistic simulation of physical systems requires recognition
and accurate simulation of the physics governing the behavior of those systems.
The computer graphics task of rendering realistic images of natural scenes
requires, according to this principle, that we model such scenes using the
appropriate geometryfractal geometry. Fractal geometry is one aspect of
the emerging field of chaos and nonlinear dynamics, a revolutionary new
area of physics and mathematics research. The remarkably realistic renderings
of scenes based on fractal geometry confirm the principle of model authenticity.
Perhaps the most noteworthy aspect of fractal geometry is the enormous range of objects and systems it is capable of describing. A partial list includes: mountain landscapes, clouds, rainfall patterns, snowflakes, frost patterns, island formations, sea coasts, river basin patterns, tree branching, arteries, veins, bronchioles, highenergy physics particle jets, percolation patterns, moon craters, galaxy structure, Brownian motion of molecules, structure of magnetic domains, radio noise, stock market fluctuations, and the mountain landscape in the movie Star Trek II: Wrath of Khan.
The science of chaos and fractal geometry is a classical example of an interdisciplinary field. The following areas have both contributed toand profited greatly fromthe study of fractal geometry:
Chaos theory and fractal geometry have also contributed significantly to the following fields:
The unifying and central role of fractal geometry is illustrated
in Figure 11.1.
Fractal geometry and chaos theory play a central and unifying role in many disciplines.
Fractal objects share certain characteristics which distinguish them from more traditional objects defined by Euclidean geometry. These include:
Let's examine some of these properties in more detail.
Our indoctrination with traditional Euclidean geometry
makes the integer classification of ndimensional systems
seem like the only natural and logical one. We all know intuitively what
1D, 2D, and 3D objects are. A line or smooth curve is a 1D object; a plane,
polygon, or Bézier patch is a 2D object ; and a sphere, cylinder,
or superquadric is a 3D object. However, Figure 11.2 suggests that life
may not always be so simple.






Sierpinski Gasket . What starts as a simple set of sixteen 1D lines becomes more complex and "space filling" with each iteration. By the sixth iteration, the 1D figure has evolved into a 2D figure.
It is clear that what starts out as a 1D outline tends toward a 2D figure as each level of recursion is carried out.
Another experiment illustrating the ambiguity of our traditional concept of dimension is the following. Start with a flat sheet of white paper. It is clearly a 2D object. Now crumple it up somewhat. What is it now? Well, it is a sort of random set of somewhat curved 2D patches. Now crumple it tighter and tighter until it is about the size of a golf ball. Hold it between your thumb and forefinger and call across the room to your friend, "Let's go play golf!" She'll reply, "Fine, I see you've got a ball." But a golf ball is 3D and this is only a sheet of paper. Is it 2D or 3D?
So the concept of dimensionality must be considered in a more systematic and analytic way. The following definition of dimension D was proposed by Hausdorff early in the twentieth century. Consider the problem of measuring the coast line of England. Using a satellite photograph and a ruler with appropriate divisions, say 1 km, we would measure a length we can designate L1. To get a more accurate answer, we could use a series of highaltitude aircraft photographs and a ruler with finer divisions of say one hundred or ten m to determine a new length, L2. If this were still not adequate, we could send out teams of surveyors with meter sticks who would report back a length, L3.
Note that L1 represents the length of coast line of the major land masses, bays, and other geologic features with a scale of a few kilometers or more. However, big bays have little bays and big peninsulas have smaller peninsulas, and L2 is a better representation of such smallerscale features. The L3 answer generated by the surveyors' meter sticks represents an even finergrained measurement. In principle, this refinement could continue indefinitely down through the centimeter and millimeter scales appropriate for rocks and pebbles until we reached the millimicron scale of atoms.
How are the measurements, L1, L2, and L3 related? Obviously, L1 < L2 < L3. That is, the result obtained is determined by the scale of the ruler used to make the measurement. In other words, the size of irregular objects depends upon the number of increments used to determine the size. This is not the case for regular Euclidean objects. This scale dependence of fractal objects provides the basis for the definition of dimensionality.
Consider the relationship between the number of incremental
elements, N, and the element size, r, for the following 1D, 2D, and 3D objects.
1D
r = 1/4, N = 4
Figure 11.3Relationship between r and N, where: r = scale size of element, N = number of elements. 
2D
r = 1/4, N = 16 
3D
r = 1/4, N = 64 
Note that the relationship between the scale r and number of elements, N, may be written:
1D: N r^{1} = 1
2D: N r^{2} = 1
3D: N r^{3} = 1
This may be generalized for dimension, D, as:
N r^{D} = 1 [11.1]
The quantity D is defined as the Hausdorff dimension and may be extracted from Equation 11.1 and expressed in terms of the number of elements and their size as:
D = ln(N)/ln(1/r) [11.2]
Note that nothing in the defining Equation 11.2 requires that D assumes only integer values. For fractal objects D takes on noninteger values. Consider D for the following fractal objects.
The von Koch triangular snowflake has some very interesting
properties. The geometric generator applied iteratively to each section
of the figure is shown in 11.4(a), and the results of the first iteration
given in 11.4(c). Note that the length of the "shoreline" (perimeter)
of the original triangle is L1
= 3; after one iteration, the perimeter is L2 = (4/3)L1; after two
iterations, L3 = (4/3)2L1;
and so on. Thus the perimeter of the von Koch snowflake goes to infinity
as the number of iterations goes to infinity, even though the curve itself
is bounded by the square subtending 11.4(c). Substituting the number of
elements used to measure a side (N = 4) and the length of each (r = 1/3),
Equation 11.2 yields the dimensionality D = 1.2618… for the von
Koch snowflake.



von Koch fractal triangle generator. Applying the geometric substitution shown in (a) to each segment of the triangle (b) produces first iteration of the von Koch snowflake, shown in (c). Since N = 4 and r = 1/3, the Hausdorff dimension of the Koch triangle is D @ 1.26.
In light of the previous discussion of the length of coastline,
one might expect that the dimensionality of an object would increase as
the complexity of substructure increases. The von Koch square, shown with
its generator and the results of the first iteration in Figure 11.5, exhibits
a complexity greater than that of the triangular snowflake. Now, to characterize
the increase in length of shoreline after each iteration, we need N = 8
elements each of length r = 1/4 unit. Substituting these values into Equation
11.2 does, in fact, yield an increased dimensionality of D = 1.5.

(b) 
von Koch fractal square generator. Applying the substitution indicated in (a) to each segment converts the square into the fractal curve shown after the first iteration in (b). Since N = 8 and r = 1/4, the dimensionality of this fractal is D = 1.5.
An intuitive interpretation of the dimensionality of fractals has been given by Richard Voss, one of Mandelbrot's colleagues at IBM, who stated, "Fractal dimension measures how the length of a coast line changes as we change the size of the ruler." We will discuss shortly algorithms for generating fractal landscapes which accurately simulate mountain scenes. An interesting result of the research on the dimensionality of natural fractal objects is the X.2 Rule. This rule states that very realistic fractal coast lines, mountains, and clouds are generated by objects with dimensionality Di, where:
Coastlines Di = 1.2
Mountains Di = 2.2
Clouds Di = 3.23.3
A natural interpretation of the fractional dimension of
coastlines involves their "area filling" character. That is, as
we move from the simple triadic Koch curve through the more complex quadric
Koch curve to the highly irregular Sierpinski curve the dimension ranges
from 1.26 through 1.5 to nearly 2. Mountain landscapes involve distorting
a smooth, 2D surface by stochastic processes acting on the third dimension.
If the disturbance is small, "rolling hills" with D ~ 2.1 result;
increasing the dimension to D > 2.5 causes a rugged landscape with jagged
peaks and steepwalled canyons. Modeling of clouds requires specifying at
each 3D point a fourth dimension such as the temperature or translucency
of the condensation at that point.
From the preceding discussion, several properties of fractals provide the basis for a classification system, including dimensionality, selfsimilarity, and generation technique. One property which cleanly distinguishes two distinct classes of fractals is the role of chance in their generation and resulting structure. Figure 11.6 summarizes this classification scheme.
Classification scheme for fractals based on the role of chance.
Deterministic fractals have structures which are fixed uniquely by the algorithm employed in their creation. That is, for a given set of parameters, a deterministic fractal generator will produce identical structures each time it is run. Chance plays no role in the final structure of the deterministic fractal.
Interestingly, identical fractals may be produced by distinct
algorithms. The Sierpinski gasket, for instance, can be generated by the
"Chaos Game" algorithm or by an iterated function system algorithm
discussed presently. Some of these algorithms do employ random number generators
within the deterministic program. The final structure of the fractal, however,
shows no indication of random processes and is a function of the control
parameters only.
The distinctive characteristic of stochastic fractals is that random processes play a central role in determining the structure of the fractal object. Phenomena such as turbulence, seacoast, and mountain formation are deterministic at a physical level but are extremely sensitive to the initial conditions, a property common to all chaotic systems. The complex interaction of the system with the initial conditions and subsequent environment results in apparently random turbulent behavior and fractal mountain landscapes.
It is virtually impossible to accurately represent the
initial condition and successfully simulate the subsequent system development
of such natural systems. A far more tractable approach has been to simulate
the fractal geometry of such objects by introducing random processes in
creating the scenes. The approach has proven successful in simulating Brownian
motion, percolation behavior, and mountain geography.
The best way to understand fractals is to study algorithms for generating the major types of fractals, build programs to implement these algorithms, and observe the graphical results. The broad categories of fractals defined in the last section can be further refined according to the general algorithmic approach used in generating the fractal.
The fractals of interest in this chapter can all be generated by the following classes of algorithms:
The basic idea of replacement mapping is that a generator function maps a given, parent structure into a new, more complex child structure. An original object is defined, and, in the first iteration, the replacement mapping is applied to each element of the original object, producing a refined child object. In the second iteration, each element of the child object is transformed into an even more refined grandchild object, and so on. Theoretically, a fractal generated by replacement mapping would require an infinite number of iterations. In practice, for graphical applications, the iteration or recursion need continue only until the structures being mapped fall to subpixel dimensions. The simplest form of replacement mapping is linear replacement mapping, the conversion of a single line into a more complex object.
The basic algorithm for linear replacement mapping is summarized below.
1. Define initial structure in terms of lines defined by end points.
2. Define replacement mapping replacing each line with refined set of lines.
3. Iterate refinement until desired level is achieved.
To implement this algorithm we add two additional conditions:
1. The algorithm should be generalized enough to map any arbitrary configuration of lines by a conversion of each line into another arbitrary set of lines.
2. The algorithm should be capable of iteration to an arbitrary level of refinement.
In order to implement the first condition, the first two steps of the algorithm can be isolated to individual subroutines. To satisfy the second condition we can use a file swapping technique to avoid limiting the problem size by the dimensions of internal data structures. The following Pascal program implements the Linear Replacement Mapping Algorithm, subject to the generality conditions.
program
Koch;
{Program to draw and iteratively refine
}
{ von Koch curves by the method of
}
{ linear replacement mapping. }
type
RealFile = file of real;
var
x1, y1, x2, y2, dx, dy: real;
x, y: integer;
i, j, Nit: integer;
infile, outfile: RealFile;
Name: string;
done: Char;
{********** Define Object ************}
procedure
DefineObject (var Name: string);
{Procedure to generate and write original
database.}
const
xmax = 512;
ymax = 342;
var
xc, yc: real;
dx, dy: real;
begin
Name := 'Triadic Koch Curve';
rewrite(outfile, 'Ch11:Koch.datA');
xc := xmax / 2;
dy := ymax / 4;
dx := 1.1547 * dy;
write(outfile, (xc  dx), 2.5 * dy);
write(outfile, xc, 0.5 * dy);
write(outfile, (xc + dx), 2.5 * dy);
write(outfile, (xc  dx), 2.5 * dy);
close(outfile);
end ;
{************** SetIO ****************}
procedure
setIO (Nit: integer; var infile, outfile: RealFile);
{Procedure to set input & output
files as a function of Nit.}
var
even: boolean;
begin {Father/Son
data files .}
even := ((Nit mod 2) = 0);
if even
then
begin
reset(infile, 'Ch11:Koch.datA');
rewrite(outfile, 'Ch11:Koch.datB');
end
else
begin
reset(infile, 'Ch11:Koch.datB');
rewrite(outfile, 'Ch11:Koch.datA');
end;
end;
{************** Refine ***************}
procedure
Refine (var outfile: RealFile; x1, y1, x2, y2: real);
{Procedure to successively refine the
input line into new output.}
const
Nout = 5;
var
i: integer;
xp, yp: array[1..Nout] of
real;
begin
dx := (x2  x1) / 3;
dy := (y2  y1) / 3;
xp[1] := x1;
yp[1] := y1;
xp[2] := x1 + dx;
yp[2] := y1 + dy;
xp[3] := xp[2] + 0.5 * dx + 0.866 *
dy;
yp[3] := yp[2]  0.866 * dx + 0.5 *
dy;
xp[4] := xp[2] + dx;
yp[4] := yp[2] + dy;
xp[5] := x2;
yp[5] := y2;
for i := 1 to 5 do
write(outfile, xp[i], yp[i]);
end;
begin {Main
program}
Nit := 0;
DefineObject(name);
repeat
hideall;
showdrawing;
setIO(Nit, infile, outfile);
{Read first point to get started.}
read(infile, x1, y1); x := round(x1);
y := round(y1);
moveto(x, y);
{Now iterate until eof(infile.}
repeat
read(infile, x2, y2);
x := round(x2);
y := round(y2);
lineto(x, y);
Refine(outfile, x1, y1, x2, y2); {Refine
line.}
x1 := x2; {Substitute last into first.}
y1 := y2;
until eof(infile);
close(infile);
close(outfile);
Nit := Nit + 1;
moveto(20, 20);
writedraw(name);
moveto(20, 30);
writedraw(' Nit = ', Nit : 1);
readln(done);
until (done
= 'd');
end
Note that the original object is defined by procedure DefineObject by writing out a set of points to a file. These points effectively define the vertices of an initial polygon. The second feature of this program is that, by "pingponging" input and output files, the complexity of the resulting fractal is limited only by the size of the storage medium and/or the compute time for a given iteration, and not by some arbitrary upper limit imposed by an internal data array.
Procedure Refine generates
the replacement shown in Figure 11.4(a). Running the program produces the
figures shown in 11.4(b) and 11.4(c) for the first two iterations. Output
for the next four iterations is shown in Figure 11.7.

(b) Nit = 4 
(c) Nit = 5 
(d) Nit
= 6 
The structure of the Koch program makes its modification
for producing other linear replacement mapped fractals a relatively simple
affair. For instance, to generate the quadric von Koch fractals of Figure
11.5, the subroutines DefineObject and MapLine which implement
steps 1 and 2 of the algorithm, respectively, are replaced by:
procedure
DefineObject (var name: string);
{Procedure to generate and write original
database.}
const
xmax = 512;
ymax = 342;
var
xc, yc: real;
d, min: real;
begin
name := 'Quadric Koch Curve';
if ymax
< xmax then
min := ymax
else
min := xmax;
rewrite(outfile, 'Ch11:Koch.datA');
xc := xmax / 2;
yc := ymax / 2;
d := min / 5;
write(outfile, (xc  d), (yc  d));
write(outfile, (xc + d), (yc  d));
write(outfile, (xc + d), (yc + d));
write(outfile, (xc  d), (yc + d));
write(outfile, (xc  d), (yc  d));
close(outfile);
end;
Substituting these two procedures and rerunning the program
produces the first two iterations shown in Figure 11.5 and the next two
iterations shown in Figure 11.8.
One of the difficulties in designing the mapping algorithm is to identify the symmetries of the problem in order to build a transformation which is independent of orientation of the original line. The LOGO turtle graphics language is the natural language for implementing such mappings. Turtle graphics commands such as FORWARD N, RIGHT q, and LEFT q, all with respect to the present turtle pointer direction, avoid the ad hoc tricks used in the MapLine procedure.
More complex space filling transformations, such as the
Sierpinski gasket, are most naturally implemented by mapping a set of
lines into a new, more complex set. This can be considered as a pattern
replacement mapping instead of a linear (line) replacement mapping.


Quadric von Koch fractal with further refinement.
Note that the deterministic fractals produced by replacement
mappings display all of the advertised properties of fractals: selfsimilarity,
complex structures arising from simple rules, a noninteger dimensionality,
and production by iteration.
Iterated function systems (IFSs) have a number of fascinating properties which make them particularly valuable tools in fractal geometry. These properties include:
Recall from earlier chapters that an affine transformation in 1D is defined as:
x´ = S x + T [11.3]
In 2D, the set of affine transformations may be written in Barnsley's notation as
[11.4]
where 1 <= i <= n = the number of sequential transformations. This may be reformulated in terms of the familiar homogeneous coordinates as
X´= X W_{i} [11.5]
where
X´= [ x´ y´ 1]
= transformed 2D coordinate vector,
X = [ x y 1]
= original 2D coordinate vector,
= contractive affine transformation.
The term contractive simply means that the effect of applying the transformation to an object is to produce a compressed image, that is, one in which any two points are closer together than they were in the original object. From the homogeneous representation of Equation 11.5, the coefficients ai, bi, ci, and di are recognized as the product of scale and rotation factors, and ei and fi as translations. The index i indicates that the first iteration of an iterated function system generally involves the simultaneous application of n distinct transformations to the present object to produce a new image. This image becomes the object of the next iteration. Each distinct sequence, Wn, generates a distinct fractal, and Barnsley tabulates the Wn for a number of different fractal objects.
1. Define an initial graphical object in pixel array, t(i,j).
2. Define set of contractive, affine transformation coefficients, (a...f)n.
3. Repeat until screen resolution reached
3.1 Repeat for all pixels of object
3.1.1 If pixel ON, repeat for all n transformations
Apply X´= X Wi
Set image, s(X´) = ON
3.2 Replace t(i,j) by s(i,j)
3.3 Plot t(i,j).
Now let's investigate the behavior of this deterministic
algorithm by applying it to several initial objects. An intuitively "natural"
object on which to iterate is the upperright triangle. Below we list a
Pascal program for implementing the IFS algorithm.
{********** Pset **************}
procedure
pset (x, y: integer);
{Procedure to plot point at pixel (x,y).}
begin
moveto(x, y);
lineto(x, y);
end;
{********* Define Object ***********}
procedure
DefineObject (var t: pic);
{Procedure to define initial graphical
object}
var
i, j: integer;
begin
{Initialize object to an upper halfdiagonal
square.}
for i :=
1 to pixdim do
for j :=
1 to pixdim do
if j <
i then
begin
t[i, j] := true;
pset(i, j);
end;
SetRect(box, 1, 1, pixdim, pixdim);
end;
{************* SetCoef ***********}
procedure
SetCoef (var a, b, c, d, e, f: vec);
{Procedure to set up coefficient of
affine transform.}
var
i: integer;
begin {Set
problem parameters for Sierpinski triangle.}
for i :=
1 to 3 do
begin
a[i] := 0.5;
b[i] := 0;
c[i] := 0;
d[i] := 0.5;
e[i] := pixdim / 2;
f[i] := 1;
end;
e[1] := 1;
f[3] := pixdim / 2;
end;
{***************Main Program***************}
begin DefineObject(t);
SetCoef(a, b, c, d, e, f);
dpix := pixdim div 2;
{Map next generation of image A(n+1)
= W(A(n))}
repeat
for i :=
1 to pixdim do
for j :=
1 to pixdim do
if t[i,
j] then
begin
for k :=
1 to 3 do
begin
x[k] := trunc(a[k] * i + b[k] * j +
e[k]);
y[k] := trunc(c[k] * i + d[k] * j +
f[k]);
s[x[k], y[k]] := true;
end;
end;
{Now plot.}
EraseRect(box);
for i :=
1 to pixdim do
for j :=
1 to pixdim do
begin
t[i, j] := s[i, j];
s[i, j] := false;
if t[i, j] then
pset(i, j);
end;
dpix := dpix div 2;
WriteLn('dpix = ', dpix);
until button
or (dpix < 1)
end.
Note that the first two steps of the IFS algorithm have
been isolated in the procedures DefineObject and SetCoef,
respectively. This permits easy extension to alternative initial objects
and sets of affine transformations. The set of three affine transformation
coefficients for producing a Sierpinski triangle are listed in Table 11.1.
W_{i} 
a 
b 
c 
d 
e 
f 
1 
0.5 
0 
0 
0.5 
1 
1 
2 
0.5 
0 
0 
0.5 
pixdim/2 
1 
3 
0.5 
0 
0 
0.5 
pixdim/2 
1 
Figure 11.9 demonstrates the iterative refinement produced
by the IFS algorithm applied to the object triangle.
(a) 
(b) 
(c) 
(d) 

(f) 
Sierpinski triangle IFS applied to triangular object.
The behavior of the Sierpinski IFS can be described as "eating away" the original solid triangle parent to produce three solid, halfsize sons and an identically sized blank space. This process is iteratively applied to each solid son at each level.
What happens if the Sierpinski triangle IFS is applied
to some initial object other than the "natural" initial triangle?
To explore the behavior of this IFS, let's select an object as "different"
from the original initial triangle as possible. One such object is the photographic
negative of the original square image which can be produced by changing
all 0s to 1s and all 1s to 0s in the original listmap. This is readily
accomplished in DefineObject by changing
the IF test from "j < i" to "j > i." Making this
single change in the IFS program results in the iteration sequence shown
in Figure 11.10.
(a) 





Iterations of the Sierpinski IFS applied to an original object which is the negative of that in Figure 11.9. Note similarity between 11.9(f) and 11.10(f).
The surprising result apparent from comparison of the final iterations in Figures 11.9 and 11.10 is that, starting from radically different initial states, the Sierpinski IFS converged to the same fractal shape. If one considers each iteration as progressive moments in time, both initial objects have moved, under the influence of the IFS, to the same fractal structure. It appears that the Sierpinski fractal triangle has some "strange attraction" for both of these initial shapes.
If these two very different objects are both attracted
to the same final state, the natural question arises: Will the Sierpinski
IFS attract any initial shape to the same fractal structure? To explore
this possibility, let's start with some whimsical initial object of a kind
entirely different from the two previous triangles. Figure 11.11 indicates
what happens to a smiley face under the influence of the Sierpinski IFS.
(a) 





Iterations of the Sierpinski IFS on an arbitrary initial object. Again, note the similarity of Figures 11.10(f) and 11.11(f).
This result confirms our speculation that the fractal Sierpinski
triangle does indeed act as a strange attractor of initial objects under
the influence of this particular IFS. This concept of strange attractor
emerged from the study of chaos and is central to the dynamic behavior of
nonlinear systems.
Strange attractors exist not only for mathematical objects
such as the Sierpinski triangle but also for natural objects such as trees
and ferns. We turn next to an outstanding example of one such natural fractal.
Barnsley has shown that the IFS specified in Table 11.2 is capable of generating a fractal pattern closely resembling the Black Spleenwort fern leaf. The algorithm for iterating the function system differs from that used above in the following ways:
Wi 
a 
b 
c 
d 
e 
f 
p 
1 
0 
0 
0 
0.16 
0 
0 
0.01 
2 
0.85 
0.04 
0.04 
0.85 
0 
1.6 
0.85 
3 
0.2 
0.26 
0.23 
0.22 
0 
1.6 
0.07 
4 
0.15 
0.28 
0.26 
0.24 
0 
0.44 
0.07 
The random iteration algorithm proposed by Barnsley for generating the fern is summarized below.
1. Define set of contractive, affine transformation coefficients,
(a...f,p)i
2. Repeat for n iterations
The main program of the Pascal implementation is given below.
The output of program Fern is shown
in Figure 11.12 for varying numbers of iterations. The early iterations
of the program seem to sprinkle points at random across the page, but soon
the hazy outline of the fern leaf appears. With an increased number of iterations
the interior points gradually fill in, and after a few hundred thousand
iterations, the leaf is essentially solid at the resolution shown.



Black Spleenwort Fern after varying numbers of iterations, n.
type
vec = array[1..4] of
real;
var
a, b, c, d, e, f, p: vec;
x, y, newx, newy, r: real;
k, xp, yp: integer;
n:LongInt;
begin
{Main Program}
{Select fractal model and read parameters.}
GetCoef(a, b, c, d, e, f, p);
{Initialize coordinates and counter.}
x := 0;
y := 0;
n := 1;
ForeColor(greenColor);
repeat
n := n + 1;
r := Rand; {Pick random r; 02r<1.}
if (r >
p[1]) and (r < (p[1] + p[2])) then
k := 2
else if
(r > (p[1] + p[2])) and (r < (p[1] + p[2] + p[3])) then
k := 3
else if
(r > (p[1] + p[2] + p[3])) then
k := 4
else
k := 1;
newx := a[k] * x + b[k] * y + e[k];
newy := c[k] * x + d[k] * y + f[k];
x := newx;
y := newy;
xp := 200 + round(75 * x);
yp := 550  round(50 * y);
pset(xp, yp);
until n=200000
end.
Note that in program Fern,
the user must supply routines SetCoef and
a random number generator, Rand. Figure 11.13
shows the output of the IFS fern program after some 8 million iterations
at a more detailed resolution. This figure demonstrates convincingly several
of the distinguishing features of fractal objects. First, the Black Spleenwort
fern exhibits an approximate selfsimilarity. That is, the main frond has
numerous subfronds, each of which have numerous subsubfronds, and so on,
with the daughter object closely resembling (but not duplicating) the parent
object.
The Black Spleenwort fern generated by the IFS algorithm.
Secondly, the fractal fern is a strong strange attractor. All of the leaves shown in Figures 11.12 and 11.13 start with point(0,0). However, it is easy to verify that the final shape of the fractal generated with the fern IFS is independent of the starting point. That is, the dynamics of the fern IFS attracts any starting point onto the identical final fractal fern. Finally, the similarity of the IFS fern fractal to natural fern leaf structures is so striking that it strongly suggests that natural ferns must exhibit fractal geometry.
Note that, although the sequence of points leading to the fractal fern (trajectory) is stochastic, the fractal itself is deterministic with a structure determined by the stra