Julia Sets

Julia sets arise from combining the concept of strange attractors with iterative mapping on the complex plane. The Julia set of the polynomial mapping function, F(z), is the boundary between the basin of attraction of the strange attractor at z = 0 and a strange attractor at z = . To illustrate what this means, first consider a case of the simple quadratic function:
F(z) = z2                                [11.12]
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: We can illustrate these three cases with the algorithm Complex_Map.


Complex_Map Algorithm

1. Draw unit circle on complex plane.

2. Request initial point, z1, from the user

3. Repeat until done

Pascal Program Complex_Map

program Complex_Map;

{Program to compute and plot}

{Complex mappings.}


scale = 100.0;

size = 400;


complex = record

r: real;

i: real



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}


c.r := a.r * b.r - a.i * b.i;

c.i := a.r * b.i + a.i * b.r;



rewrite(out, 'Ch11:Complex.dat');


writeln('Initial Value of z: (Zr,Zi)?');

ReadLn(z.r, z.i);


Nit := 1;


sd4 := size div 4;

sd2 := size div 2;

SetRect(circle, sd4, sd4, 3 * sd4, 3 * sd4);



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


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


until button;



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

Figure 11.15

Mapping of initial point with modulus < 1.0. Note effect of strange attractor at z = (0,0).

Figure 11.16

Trajectory of mapping for initial point of modulus > 1.0. Note effect of strange attractor at .


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 systems--the 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 three-body 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 effect--the disturbance caused by a butterfly fluttering its wings in China can, in principle, propagate into a hurricane in the Caribbean Sea!

Interesting Julia Sets

In the above discussion, the unit circle is shown to be the Julia set for the complex transformation, F(z)= z2. 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:

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.



Figure 11.17

Initial 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.2e-1 3.8e-1  

2 7.0e-1 7.1e-1  

3 -8.4e-3 1.0e+0  

4 -1.0e+0 -1.7e-2  

5 1.0e+0 3.3e-2  

6 1.0e+0 6.7e-2  

7 9.9e-1 1.3e-1  

8 9.6e-1 2.6e-1 

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.3e-3 3.3e-3  

31 -9.2e-6 -8.4e-6  

32 1.3e-11 1.5e-10  

33 -2.4e-20 3.9e-21  

34 5.5e-40 -1.8e-40  

35 0.0e+0 -0.0e+0 


Backward mapping

The backward mapping algorithm is suggested by an examination of Figures 11.15-11.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 10-30 iterations and spend the rest of its iteration life hopping from point to point on the set.

Since forward iteration is defined as zn+1 = (zn)2+ c, to reverse the direction of motion along the trajectory, we must write:

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:

Backward mapping algorithm 

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

3. Repeat until set is outlined distinctly

{Julia set will "grow in" point by point.}


Forward mapping

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 filled-in Julia set. The algorithm for forward mapping is a relatively simple extension of the Complex_Map algorithm.


Forward Mapping Algorithm

1. Request starting parameters

2. For each pixel in the 2 ¥ 2 complex plane


Boundary Scanning

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:

Boundary Scanning Algorithm 

1. Select an M ¥ N pixel grid to scan a 2 ¥ 2 complex plane.

2. For each (m,n) pixel on this grid


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 filled-in Julia sets produces more visually appealing images and is implemented below.

Pascal Program Julia 

program Julia;

{Program to compute and plot Julia set.}


scale = 0.01;

R = 10;


complex = record

r: real;

i: real



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}


c.r := a.r * b.r - a.i * b.i;

c.i := a.r * b.i + a.i * b.r;


procedure add (a, b: complex; var c:complex);

{Does complex addition: c = a + b}


c.r := a.r + b.r;

c.i := a.i + b.i;


procedure plot (c, r: integer);

{Procedure to pixel (c.r).}


moveto(c, r);

lineto(c, r);




writeLn('How many iterations at each point?');


writeLn('Value of C: (Cr,Ci)?');

ReadLn(c.r, c.i);



n := 0;

for col := 1 to 400 do

for row := 1 to 400 do


z.r := (col - 200) * scale;

z.i := (row - 200) * scale;


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;




(a) Nit = 100, c = (-1,0)
(b) Nit = 100, c = (-0.9,0.12)
Nit = 100, 
c = (0.25,0) 
(d) Nit = 10, 
c = (0,1) 
(e) Nit = 10, c = (-.4,-.6) 
(f) Nit =20, c = (-.4,-.6) 
(g) Nit = 50,c = (-.4,-.6) 

Figure 11.18

Julia sets (filled-in) 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) = z2 + c. The approximate self-similarity appears in each figure--lobes have sub-lobes which have sub-sublobes, branches have sub-branches which have sub-sub-branches, 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 round-off 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.

Figure 11.19

Julia set for c = (-0.4,0.6) using color-coding, 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.

Figure 11.20

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 filled-in 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

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,

do not cause the complex number zn 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:

Mandelbrot Set Algorithm

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


Implementing this algorithm produces the stark, beetle-like Mandelbrot set representing those complex m values for which the attractor at in complex z-space 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 color-code 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 color-coded Mandelbrot set generator.

Pascal Program Mandelbrot

program Mandelbrot;

{Program to compute and plot the Mandelbrot set.


Nit = 100;

scale = 0.005;

R = 10;


complex = record

r: real;

i: real



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}


c.r := a.r * b.r - a.i * b.i;

c.i := a.r * b.i + a.i * b.r;


procedure sub (a, b: complex; var c: complex);

{Does complex subtraction: c = a - b}


c.r := a.r - b.r;

c.i := a.i - b.i;


procedure plot (c, r, n: integer);

{Procedure to pixel (c.r) in color code n.}


case n of




















moveto(c, r);

lineto(c, r);

r := 400 - r;

moveto(c, r);

lineto(c, r);



for col := 1 to 400 do

for row := 1 to 200 do


z.r := 0.0;

z.i := 0.0;

c.r := (col - 100) * scale;

c.i := (200 - row) * scale;


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;




Figure 11.21

The Mandelbrot Set (black) with color-coded 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.)

Figure 11.22

Segment of the Mandelbrot set with corners at

c1 = (-0.64167, 0.6930) and c2= (-0.453,0.511). 


Figure 11.23

Segment of the Mandelbrot set with corners at

c1 = (-0..5963,0.5821) and c2 = (-0.5550,0.5390). 


Figure 11.24

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

Figure 11.25

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. 


Figure 11.26

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 twenty-five 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) = z2 + 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.

Figure 11.27

The Mandelbrot set as catalog of Julia sets. The Julia set in the lower, right-hand corner of the figure, for example, corresponds to (cr,ci) = (-0.13,0.64). 


Stochastic Fractals

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 self-similar 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 computer-generated (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 fire-storms 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.

Brownian Motion

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 random-walk problem. The random-walk problem, in its simplest 1D form, involves the motion resulting from the flip of a coin. If it's heads, the coin-flipper takes one step ahead; if it's tails, s/he steps back one step. The resulting motion is called a random-walk. 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 random-walk 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.

Random-walk Algorithm

1. Set d(0) = 0 for time t = 0.

2. For t = 1 to 256

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 y-axis and time along the x-axis. A selected example of the output of Random_Walk is shown in Figure 11.28.

Pascal Program Random_Walk

program Random_Walk;

{Program to generate Brownian trace}

{by Random-Walk algorithm.}


xmax = 256;

ymax = 200;

step = 2;


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


ran := random / 32768


begin {Main program}

pensize(2, 2);

d := 0;

{Get "Random" randomized.}


randSeed := time.Minute * 60 + time.Second;

for i := 1 to xmax do


flip := ran;

if flip > 0 then

x := step


x := -step;

d := d + x;

xp := i;

yp := ymax div 2 - d;

drawline(xp, yp, xp, yp);



Figure 11.28

Brownian trace generated by Random_Walk.

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.

Random Midpoint Displacement Algorithm

1. Initialize parameters

2. Repeat for n = 1 to Nit {i.e., subdivide line until we hit pixel level}

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.

Pascal Program MidPoint

program MidPoint;

{Program to generate Brownian trace by}

{method of Random Midpoint Displacement.}


xmax = 256;

ymax = 200;


image = array[0..xmax] of real;


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}


k: integer;

dy, ave: real;


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]));


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


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



Figure 11.29

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 random-walk 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 .



Figure 11.30

Random midpoint displacement algorithm for generating fractal landscapes. Midpoints are displaced randomly in sign and magnitude in the vertical direction--here 2 and 3 upwards and 1 downwards. The process is repeated on each of the four resulting triangles.

Fractal Mountain Landscapes

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.

Fractal Mountain Landscape RMD Algorithm

1. Initialize starting parameters

2. Repeat until desired resolution is obtained

3. Project resulting 3D model onto a 2D image

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 four-triangle 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 left-most 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 side-effect violates the physical constraint that the surface must be continuous.

(a) Nit = 1

(b) Nit = 3

(c) Nit = 5

(d) Nit= 7

(e) Nit = 7 (shaded)

(f) Nit = 7 (shaded, with sea level = 0)

Figure 11.31

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 z-value 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 height-encoded 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.

Fractal Applications

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.

Figure 11.32

Control window for fractal line option of PixelPaint Professional. 

Figure 11.33

Lightning scene generated with fractal line option of PixelPaint Professional. Randomness was set at 82%, and the line width was manually reduced as the "tree" was traversed. 

Fractal Design Tools

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.

Brownian Trace Generator

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







Figure 11.34

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.







Figure 11.35

Fractal trees generated using PixelPaint program. 

Fractal Tree Generators

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.

Figure 11.36

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.

Fractal Image Compression Techniques

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 twenty-eight 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, e-mail. 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 e-mail 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 twenty-eight 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 e-mail to your colleague? The answer to this question involves solving what Barnsley defines as the inverse problem.

The Inverse Problem and the Collage Theorem

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:

Given an arbitrary image, what set of IFS transformations will encode it?

The inverse problem can be conceptualized as the bottom row of Figure 11.37.


(a) The fractal research agenda


(b) The inverse problem

Figure 11.37

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.


Table 11.3


Figure 11.38

Collage Theorem used to encode a square as a 4-transformation 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.

Collage Algorithm

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

The collage, T´, is formally defined as the union of the Wn transformations:



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 2-color level image (i.e., all black interior),

2. Apply an image outline tracer,

3. Convert the bit-mapped 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.2-2.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!





Figure 11.39

The Collage Algorithm. Generating IFS code for a leaf. 

Automating Fractal Image Encoding

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 dual-68020 system. This would appear to rule our fractal image compression as a feasible technique, at least for personal workstations.

Iterated System's P.OEM™

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 fractal-transform 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.

Figure 11.40

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 N-dimensional 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.