Chapter 11

Fractals -

Objects of Fractional Dimension

Morris Firebaugh

© William C. Brown Publishers, Inc.

 

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 humanly impossible to understand a single word of it; without these, one is wandering about in a dark labyrinth.

Galileo

 

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. ... Mathematicians have disdained this challenge, however, and have increasingly chosen to flee from nature by devising theories unrelated to anything we can see or feel. ...


Responding to this challenge, I conceived and developed a new geometry of nature and implemented its use in a number of diverse fields.

Benoit Mandelbrot

 

What is new is the ability to start with an actual image and find the fractals that will imitate it to any desired degree of accuracy. Since our method includes a compact way of representing these fractals, we end up with a highly compressed data set for reconstructing the original image.

Michael Barnsley and Alan Sloan



Fractals play a very special role in computer graphics. They form the basis for some of the most beautiful and exquisite images ever generated by humans or machines. Computer graphics provided the most significant tool used in the discovery of fractal geometry. Because many natural objects are fractal in nature, fractal geometry gives us new tools for modeling natural scenes. With the solution of the inverse problem, fractal representations provide an image compression technique for solving one of the most troublesome problems of computer graphics--the huge size of graphics image files.

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., computer-generated 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 human-made, 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 fractional-dimension objects is the natural concluding chapter of our N-Dimensional 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 geometry--fractal 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.
 
 

The Fractal Geometry of Nature

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, high-energy 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 to--and profited greatly from--the 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.
 


Figure 11.1

Fractal geometry and chaos theory play a central and unifying role in many disciplines.


 

Properties of Fractals

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.
 

Fractional Dimensions

Our indoctrination with traditional Euclidean geometry makes the integer classification of n-dimensional 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.

 


(a)


(b)


(c)


(d)


(e)


(f)

Figure 11.2

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 high-altitude 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 smaller-scale features. The L3 answer generated by the surveyors' meter sticks represents an even finer-grained 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.3

Relationship 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 r1 = 1

2D: N r2 = 1

3D: N r3 = 1



This may be generalized for dimension, D, as:

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:

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&hellip; for the von Koch snowflake.
 



(a) 



(b) 



(c) 

Figure 11.4

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 sub-structure 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.

 

 



(a)



(b)

Figure 11.5

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:

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 steep-walled 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.
 

Classification of Fractals

From the preceding discussion, several properties of fractals provide the basis for a classification system, including dimensionality, self-similarity, 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.


Figured 11.6

Classification scheme for fractals based on the role of chance.



Deterministic Fractals

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.
 

Stochastic Fractals

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.
 

Generating Fractals

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.

Categories of Fractal Algorithms

The fractals of interest in this chapter can all be generated by the following classes of algorithms:

Linear Replacement Mapping

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.


Linear Replacement Mapping Algorithm

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:

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.


Pascal Program Koch


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 "ping-ponging" 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.




(a) Nit = 3 

(b) Nit = 4


(c) Nit = 5

(d) Nit = 6

Figure 11.7

Triadic von Koch curves ("snowflake") for iterations 3 - 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; 



and



procedure MapLine (var outfile: RealFile; x1, y1, x2, y2: real);
{Procedure to successively refine the input line into new output.}
const
Nout = 9;
var
i, j: integer;
xp, yp: array[1..Nout] of real;
begin
dx := (x2 - x1) / 4;
dy := (y2 - y1) / 4;
xp[1] := x1;
yp[1] := y1;
xp[2] := xp[1] + dx;
yp[2] := yp[1] + dy;
xp[3] := xp[2] + dy;
yp[3] := yp[2] - dx;
xp[4] := xp[3] + dx;
yp[4] := yp[3] + dy;
xp[5] := xp[4] - dy;
yp[5] := yp[4] + dx;
xp[6] := xp[5] - dy;
yp[6] := yp[5] + dx;
xp[7] := xp[6] + dx;
yp[7] := yp[6] + dy;
xp[8] := xp[7] + dy;
yp[8] := yp[7] - dx;
xp[9] := x2;
yp[9] := y2;
for i := 1 to Nout do
write(outfile, xp[i], yp[i]);
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.

 


(a) Nit = 3 


(b) Nit = 4 

Figure 11.8

Quadric von Koch fractal with further refinement.




Note that the deterministic fractals produced by replacement mappings display all of the advertised properties of fractals: self-similarity, complex structures arising from simple rules, a non-integer dimensionality, and production by iteration.
 

Iterated Function Systems

Iterated function systems (IFSs) have a number of fascinating properties which make them particularly valuable tools in fractal geometry. These properties include:

  • IFSs generate fractals through a repeated application of a set of contractive affine transformations.


  • IFSs can generate both rigid, geometric shapes and realistic, life-like forms. We illustrate these two extremes with the Sierpinski triangle and the Black Spleenwort fern.


  • IFSs can generate fractals using either deterministic procedures or random (chaotic) processes. Remarkably, these two different algorithms generate the same fractal!


  • IFSs demonstrate the concept of strange attractors. This concept has emerged from the new science of dynamic systems and chaos.


  • IFSs provide the basis for a highly efficient image compression system developed by Michael Barnsley.

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 Wi					 [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.


Iterated Function System Algorithm

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).



Sierpinski Triangle

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 upper-right triangle. Below we list a Pascal program for implementing the IFS algorithm.
 

Pascal Program IFS


program IFS;
{Program to compute fractals using IFSs.}
{Ref.: Barnsley, "Fractals Everywhere , p. 88}
const
pixdim = 120;
type
pic = array[1..pixdim, 1..pixdim] of Boolean;
vec = array[1..4] of real;
dimvec = array[1..4] of integer;
var
s, t: pic;
a, b, c, d, e, f, p: vec;
x, y: dimvec;
i, j, k, dpix: integer;
box: rect;

{********** 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 half-diagonal 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.
 

Table 11.1

Coefficients of IFS map for Sierpinski Triangle Fractal

Wi 

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 

where pixdim = size of original object in pixels.
 

Figure 11.9 demonstrates the iterative refinement produced by the IFS algorithm applied to the object triangle.

 

(a)

(b)

(c)

(d) 


(e) 


(f) 

Figure 11.9

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, half-size 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 list-map. 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)


(b)


(c)


(d)


(e)


(f)

Figure 11.10

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)


(b)


(c)


(d)


(e)


(f)

Figure 11.11

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 non-linear 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.
 

Black Spleenwort Fern

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:

  • Rather than mapping the whole object window to a new image window with each iteration, only a single point is mapped by each transformation.
  • The particular transformation used to map the point is selected randomly according to an a priori probability specified for each affine transformation of the set. This probability is indicated as "p" in Table 11.2.

Table 11.2

Coefficients of IFS map for Fern Fractal

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.


Random Iteration IFS Algorithm

1. Define set of contractive, affine transformation coefficients, (a...f,p)i
2. Repeat for n iterations

    2.1 Pick transformation i randomly with probability pi
    2.2 Apply X´= X Wi
    2.3 Plot point (xi´, yi´)
    2.4 Set X = X´.



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.
 


(a) n = 2000 


(b) n = 20,000 


(c) n = 200,000 

Figure 11.12

Black Spleenwort Fern after varying numbers of iterations, n.




Pascal Program Fern


program Fern;
{Program to generate fractal fern.}
{Ref.: Barnsley, "Fractals Everywhere"}

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 self-similarity. That is, the main frond has numerous subfronds, each of which have numerous sub-subfronds, and so on, with the daughter object closely resembling (but not duplicating) the parent object.
 


Figure 11.13

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